Note: This page has been translated by MathWorks. Click here to see

To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

Solve differential equations by using Laplace transforms in Symbolic Math
Toolbox™ with this workflow. For simple examples on the Laplace transform, see
`laplace`

and `ilaplace`

.

The Laplace transform of a function *f*(*t*) is

$$F\left(s\right)={\displaystyle \underset{0}{\overset{\infty}{\int}}f(t){e}^{-ts}dt}.$$

Symbolic workflows keep calculations in the natural symbolic form instead of numeric form. This approach helps you understand the properties of your solution and use exact symbolic values. You substitute numbers in place of symbolic variables only when you require a numeric result or you cannot continue symbolically. For details, see Choose Symbolic or Numeric Arithmetic. Typically, the steps are:

Declare equations.

Solve equations.

Substitute values.

Plot results.

Analyze results.

You can use the Laplace transform to solve differential equations with initial conditions. For example, you can solve resistance-inductor-capacitor (RLC) circuits, such as this circuit.

Resistances in ohm:

*R*,_{1}*R*,_{2}*R*_{3}Currents in ampere:

*I*,_{1}*I*,_{2}*I*_{3}Inductance in henry:

*L*Capacitance in farad:

*C*Electromotive force in volts:

*E*(*t*)Charge in coulomb:

*Q*(*t*)

Apply Kirchhoff's voltage and current laws to get the differential equations for the RLC circuit.

$$\frac{d{I}_{1}}{dt}+\frac{{R}_{2}}{L}\frac{dQ}{dt}=\frac{{R}_{2}-{R}_{1}}{L}{I}_{1}.$$

$$\frac{dQ}{dt}=\frac{1}{{R}_{3}+{R}_{2}}\left(E(t)-\frac{1}{C}Q(t)\right)+\frac{{R}_{2}}{{R}_{3}+{R}_{2}}{I}_{1}.$$

Declare the variables. Because the physical quantities have positive values,
set the corresponding assumptions on the variables. Let *E*(*t*) be an alternating voltage of 1 V.

syms L C I1(t) Q(t) s R = sym('R%d',[1 3]); assume([t L C R] > 0) E(t) = 1*sin(t); % Voltage = 1 V

Declare the differential equations.

dI1 = diff(I1,t); dQ = diff(Q,t); eqn1 = dI1 + (R(2)/L)*dQ == (R(2)-R(1))/L*I1 eqn2 = dQ == (1/(R(2)+R(3))*(E-Q/C)) + R(2)/(R(2)+R(3))*I1

eqn1(t) = diff(I1(t), t) + (R2*diff(Q(t), t))/L == -(I1(t)*(R1 - R2))/L eqn2(t) = diff(Q(t), t) == (sin(t) - Q(t)/C)/(R2 + R3) + (R2*I1(t))/(R2 + R3)

Assume that the initial current and charge,
*I _{0}* and

`0`

. Declare these initial conditions.cond1 = I1(0) == 0 cond2 = Q(0) == 0

cond1 = I1(0) == 0 cond2 = Q(0) == 0

Compute the Laplace transform of `eqn1`

and
`eqn2`

.

eqn1LT = laplace(eqn1,t,s) eqn2LT = laplace(eqn2,t,s)

eqn1LT = s*laplace(I1(t), t, s) - I1(0) - (R2*(Q(0) - s*laplace(Q(t), t, s)))/L == ... -((R1 - R2)*laplace(I1(t), t, s))/L eqn2LT = s*laplace(Q(t), t, s) - Q(0) == (R2*laplace(I1(t), t, s))/(R2 + R3) + ... (C/(s^2 + 1) - laplace(Q(t), t, s))/(C*(R2 + R3))

The function `solve`

solves only for symbolic variables.
Therefore, to use `solve`

, first substitute
`laplace(I1(t),t,s)`

and
`laplace(Q(t),t,s)`

with the variables
`I1_LT`

and `Q_LT`

.

syms I1_LT Q_LT eqn1LT = subs(eqn1LT,[laplace(I1,t,s) laplace(Q,t,s)],[I1_LT Q_LT])

eqn1LT = I1_LT*s - I1(0) - (R2*(Q(0) - Q_LT*s))/L == -(I1_LT*(R1 - R2))/L

eqn2LT = subs(eqn2LT,[laplace(I1,t,s) laplace(Q,t,s)],[I1_LT Q_LT])

eqn2LT = Q_LT*s - Q(0) == (I1_LT*R2)/(R2 + R3) - (Q_LT - C/(s^2 + 1))/(C*(R2 + R3))

Solve the equations for `I1_LT`

and
`Q_LT`

.

eqns = [eqn1LT eqn2LT]; vars = [I1_LT Q_LT]; [I1_LT, Q_LT] = solve(eqns,vars)

I1_LT = (R2*Q(0) + L*I1(0) - C*R2*s + L*s^2*I1(0) + R2*s^2*Q(0) + C*L*R2*s^3*I1(0) + ... C*L*R3*s^3*I1(0) + C*L*R2*s*I1(0) + C*L*R3*s*I1(0))/((s^2 + 1)*(R1 - R2 + L*s + ... C*L*R2*s^2 + C*L*R3*s^2 + C*R1*R2*s + C*R1*R3*s - C*R2*R3*s)) Q_LT = (C*(R1 - R2 + L*s + L*R2*I1(0) + R1*R2*Q(0) + R1*R3*Q(0) - R2*R3*Q(0) + ... L*R2*s^2*I1(0) + L*R2*s^3*Q(0) + L*R3*s^3*Q(0) + R1*R2*s^2*Q(0) + R1*R3*s^2*Q(0) - ... R2*R3*s^2*Q(0) + L*R2*s*Q(0) + ... L*R3*s*Q(0)))/((s^2 + 1)*(R1 - R2 + L*s + C*L*R2*s^2 + C*L*R3*s^2 + ... C*R1*R2*s + C*R1*R3*s - C*R2*R3*s))

Calculate *I _{1}* and

`I1_LT`

and `Q_LT`

. Simplify the
result. Suppress the output because it is long.I1sol = ilaplace(I1_LT,s,t); Qsol = ilaplace(Q_LT,s,t); I1sol = simplify(I1sol); Qsol = simplify(Qsol);

Before plotting the result, substitute symbolic variables by the numeric
values of the circuit elements. Let *R**1*
= 4 Ω , *R2* = 2 Ω,
*R3* = 3 Ω, *C* =
1/4 F, *L* = 1.6 H, *I _{1}*(0) = 15 A , and

vars = [R L C I1(0) Q(0)]; values = [4 2 3 1.6 1/4 15 2]; I1sol = subs(I1sol,vars,values) Qsol = subs(Qsol,vars,values)

I1sol = 15*exp(-(51*t)/40)*(cosh((1001^(1/2)*t)/40) - (293*1001^(1/2)*sinh((1001^(1/2)*t)/40))/21879) - (5*sin(t))/51 Qsol = (4*sin(t))/51 - (5*cos(t))/51 + (107*exp(-(51*t)/40)*(cosh((1001^(1/2)*t)/40) + (2039*1001^(1/2)*sinh((1001^(1/2)*t)/40))/15301))/51

Plot the current `I1sol`

and charge
`Qsol`

. Show both the transient and steady state behavior
by using two different time intervals: 0 ≤ `t`

≤ 10 and 5 ≤ `t`

≤ 25.

subplot(2,2,1) fplot(I1sol,[0 10]) title('Current') ylabel('I1(t)') xlabel('t') subplot(2,2,2) fplot(Qsol,[0 10]) title('Charge') ylabel('Q(t)') xlabel('t') subplot(2,2,3) fplot(I1sol,[5 25]) title('Current') ylabel('I1(t)') xlabel('t') text(7,0.25,'Transient') text(16,0.125,'Steady State') subplot(2,2,4) fplot(Qsol,[5 25]) title('Charge') ylabel('Q(t)') xlabel('t') text(7,0.25,'Transient') text(15,0.16,'Steady State')

Initially, the current and charge decrease exponentially. However, over the long term, they are oscillatory. These behaviors are called "transient" and "steady state", respectively. With the symbolic result, you can analyze the result's properties, which is not possible with numeric results.

Visually inspect `I1sol`

and `Qsol`

. They
are a sum of terms. Find the terms by using `children`

.
Then, find the contributions of the terms by plotting them over ```
[0
15]
```

. The plots show the transient and steady state
terms.

I1terms = children(I1sol); Qterms = children(Qsol); subplot(1,2,1) fplot(I1terms,[0 15]) ylim([-2 2]) title('Current terms') subplot(1,2,2) fplot(Qterms,[0 15]) ylim([-2 2]) title('Charge terms')

The plots show that `I1sol`

has a transient and steady
state term, while `Qsol`

has a transient and two steady
state terms. From visual inspection, notice `I1sol`

and
`Qsol`

have a term containing the
`exp`

function. Assume that this term causes the
transient exponential decay. Separate the transient and steady state terms
in `I1sol`

and `Qsol`

by checking terms
for `exp`

using `has`

.

I1transient = I1terms(has(I1terms,'exp')) I1steadystate = I1terms(~has(I1terms,'exp'))

I1transient = 15*exp(-(51*t)/40)*(cosh((1001^(1/2)*t)/40) - (293*1001^(1/2)*sinh((1001^(1/2)*t)/40))/21879) I1steadystate = -(5*sin(t))/51

Similarly, separate `Qsol`

into transient and steady state
terms. This result demonstrates how symbolic calculations help you analyze
your problem.

Qtransient = Qterms(has(Qterms,'exp')) Qsteadystate = Qterms(~has(Qterms,'exp'))

Qtransient = (107*exp(-(51*t)/40)*(cosh((1001^(1/2)*t)/40) + (2039*1001^(1/2)*sinh((1001^(1/2)*t)/40))/15301))/51 Qsteadystate = [ -(5*cos(t))/51, (4*sin(t))/51]