Main Content

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)=0f(t)e-tsdt.

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:

  1. Declare equations.

  2. Solve equations.

  3. Substitute values.

  4. Plot results.

  5. 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: R1,R2,R3

  • Currents in ampere: I1,I2,I3

  • Inductance in henry: L

  • Capacitance in farad: C

  • AC voltage source in volts: E(t)

  • Capacitor charge in coulomb: Q(t)

Apply Kirchhoff's voltage and current laws to get the following equations.

I1=I2+I3LdI1dt+I1R1+I2R2=0E(t)+I2R2-QC-I3R3=0

Substitute the relation I3=dQ/dt (which is the rate of the capacitor being charged) to the above equations to get the differential equations for the RLC circuit.

dI1dt-R2LdQdt=-R1+R2LI1dQdt=1R2+R3(E(t)-QC)+R2R2+R3I1

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);        % AC voltage = 1 V

Declare the differential equations.

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

t I1(t)-R2t Q(t)L=-I1(t)R1+R2L

eqn2 = dQ == (1/(R(2)+R(3))*(E-Q/C)) + R(2)/(R(2)+R(3))*I1
eqn2(t) = 

t Q(t)=sin(t)-Q(t)CR2+R3+R2I1(t)R2+R3

Solve Equations

Compute the Laplace transform of eqn1 and eqn2.

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

slaplace(I1(t),t,s)-I1(0)+R2Q(0)-slaplace(Q(t),t,s)L=-R1+R2laplace(I1(t),t,s)L

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

slaplace(Q(t),t,s)-Q(0)=R2laplace(I1(t),t,s)R2+R3+Cs2+1-laplace(Q(t),t,s)CR2+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,LTs-I1(0)+R2Q(0)-QLTsL=-I1,LTR1+R2L

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

QLTs-Q(0)=I1,LTR2R2+R3-QLT-Cs2+1CR2+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 = 

LI1(0)-R2Q(0)+CR2s+Ls2I1(0)-R2s2Q(0)+CLR2s3I1(0)+CLR3s3I1(0)+CLR2sI1(0)+CLR3sI1(0)s2+1R1+R2+Ls+CLR2s2+CLR3s2+CR1R2s+CR1R3s+CR2R3s

Q_LT = 

CR1+R2+Ls+LR2I1(0)+R1R2Q(0)+R1R3Q(0)+R2R3Q(0)+LR2s2I1(0)+LR2s3Q(0)+LR3s3Q(0)+R1R2s2Q(0)+R1R3s2Q(0)+R2R3s2Q(0)+LR2sQ(0)+LR3sQ(0)s2+1R1+R2+Ls+CLR2s2+CLR3s2+CR1R2s+CR1R3s+CR2R3s

Compute I1 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/4F, L=1.6H. Assume that the initial current is I1(0)=2A and the initial charge is Q(0)=2C.

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

200cos(t)8161+405sin(t)8161+16122e-81t40cosh(1761t40)-7425291761sinh(1761t40)141954218161

Qsol = subs(Qsol,vars,values)
Qsol = 

924sin(t)8161-1055cos(t)8161+17377e-81t40cosh(1761t40)+11094251761sinh(1761t40)306008978161

Plot Results

Plot the current I1sol and charge Qsol. Show both the transient and steady state behavior by using two different time intervals: 0t15 and 2t25.

subplot(2,2,1)
fplot(I1sol,[0 15])                      
title('Current')
ylabel('I1(t)')
xlabel('t')

subplot(2,2,2)
fplot(Qsol,[0 15])                        
title('Charge')
ylabel('Q(t)')
xlabel('t')

subplot(2,2,3)
fplot(I1sol,[2 25])                  
title('Current')
ylabel('I1(t)')
xlabel('t')
text(3,-0.1,'Transient')
text(15,-0.07,'Steady State')

subplot(2,2,4)
fplot(Qsol,[2 25])                        
title('Charge')
ylabel('Q(t)')
xlabel('t')
text(3,0.35,'Transient')
text(15,0.22,'Steady State')

Figure contains 4 axes objects. Axes object 1 with title Current contains an object of type functionline. Axes object 2 with title Charge contains an object of type functionline. Axes object 3 with title Current contains 3 objects of type functionline, text. Axes object 4 with title Charge contains 3 objects of type functionline, text.

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);
I1terms = [I1terms{:}];
Qterms = children(Qsol);
Qterms = [Qterms{:}];

figure;
subplot(1,2,1)
fplot(I1terms,[0 15])
ylim([-0.5 2.5])
title('Current terms')

subplot(1,2,2)
fplot(Qterms,[0 15])
ylim([-0.5 2.5])
title('Charge terms')

Figure contains 2 axes objects. Axes object 1 with title Current terms contains 3 objects of type functionline. Axes object 2 with title Charge terms contains 3 objects of type functionline.

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

16122e-81t40cosh(1761t40)-7425291761sinh(1761t40)141954218161

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

(200cos(t)8161405sin(t)8161)

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

17377e-81t40cosh(1761t40)+11094251761sinh(1761t40)306008978161

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

(-1055cos(t)8161924sin(t)8161)