Solve Differential Equations Using Laplace Transform

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.

Definition: Laplace Transform

The Laplace transform of a function f(t) is

$F (s) = \int_{0}^{\infty} f (t) e^{- t s} d t .$

Concept: Using Symbolic Workflows

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 Numeric or Symbolic Arithmetic. Typically, the steps are:

Declare equations.
Solve equations.
Substitute values.
Plot results.
Analyze results.

Workflow: Solve RLC Circuit Using Laplace Transform

Declare Equations

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₁, R₂, R₃
Currents in ampere: I₁, I₂, I₃
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}}{d t} + \frac{R_{2}}{L} \frac{d Q}{d t} = \frac{R_{2} - R_{1}}{L} I_{1} .$

$\frac{d Q}{d t} = \frac{1}{R_{3} + R_{2}} (E (t) - \frac{1}{C} Q (t)) + \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₀ and Q₀, are both 0. Declare these initial conditions.

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

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

Solve Equations

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₁ and Q by computing the inverse Laplace transform of 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);

Substitute Values

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

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 Results

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')

Analyze Results

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]

Documentation