Documentation Center

  • Trial Software
  • Product Updates

Compute Laplace and Inverse Laplace Transforms

The Laplace transform of a function f(t) is defined as

while the inverse Laplace transform (ILT) of f(s) is

where c is a real number selected so that all singularities of f(s) are to the left of the line s =  c. The notation L[f] indicates the Laplace transform of f at s. Similarly, L–1[f] is the ILT of f at t.

The Laplace transform has many applications including the solution of ordinary differential equations/initial value problems. Consider the resistance-inductor-capacitor (RLC) circuit below.

Let Rj and Ij, j = 1, 2, 3 be resistances (measured in ohms) and currents (amperes), respectively; L be inductance (henrys), and C be capacitance (farads); E(t) be the electromotive force, and Q(t) be the charge.

By applying Kirchhoff's voltage and current laws, Ohm's Law, and Faraday's Law, you can arrive at the following system of simultaneous ordinary differential equations.

Solve this system of differential equations using laplace. First treat the Rj, L, and C as (unknown) real constants and then supply values later on in the computation.

clear E
syms R1 R2 R3 L C real;
syms I1(t) Q(t) s;
dI1(t) = diff(I1(t), t);
dQ(t) = diff(Q(t),t);
E(t) = sin(t);  % Voltage
eq1(t) = dI1(t) + R2*dQ(t)/L - (R2 - R1)*I1(t)/L;
eq2(t) = dQ(t) - (E(t) - Q/C)/(R2 + R3) - R2*I1(t)/(R2 + R3);

At this point, you have constructed the equations in the MATLAB® workspace. An approach to solving the differential equations is to apply the Laplace transform, which you will apply to eq1(t) and eq2(t). Transforming eq1(t) and eq2(t)

L1(t) = laplace(eq1,t,s)
L2(t) = laplace(eq2,t,s)

returns

L1(t) = 
s*laplace(I1(t), t, s) - I1(0) 
+ ((R1 - R2)*laplace(I1(t), t, s))/L 
- (R2*(Q(0) - s*laplace(Q(t), t, s)))/L

L2(t) =
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))

Now you need to solve the system of equations L1 = 0, L2 = 0 for laplace(I1(t),t,s) and laplace(Q(t),t,s), the Laplace transforms of I1 and Q, respectively. To do this, make a series of substitutions. For the purposes of this example, use the quantities R1 = 4 Ω (ohms), R2 = 2 Ω, R3 =  3 Ω, C = 1/4 farads, L = 1.6 H (henrys), I1(0) = 15 A (amperes), and Q(0) = 2 A*sec. Substituting these values in L1

syms LI1 LQ
NI1 = subs(L1(t),{R1,R2,R3,L,C,I1(0),Q(0)}, ...
      {4,2,3,1.6,1/4,15,2})

returns

NI1 =
s*laplace(I1(t), t, s) + (5*s*laplace(Q(t), t, s))/4
 + (5*laplace(I1(t), t, s))/4 - 35/2

The substitution

NQ = subs(L2,{R1,R2,R3,L,C,I1(0),Q(0)},{4,2,3,1.6,1/4,15,2})

returns

NQ(t) =
s*laplace(Q(t), t, s) - 1/(5*(s^2 + 1)) -...
(2*laplace(I1(t), t, s))/5 + (4*laplace(Q(t), t, s))/5 - 2

To solve for laplace(I1(t),t,s) and laplace(Q(t),t,s), make a final pair of substitutions. First, replace the strings laplace(I1(t),t,s) and laplace(Q(t),t,s) by the sym objects LI1 and LQ, using

NI1 =...
 subs(NI1,{laplace(I1(t),t,s),laplace(Q(t),t,s)},{LI1,LQ})

to obtain

NI1 =
(5*LI1)/4 + LI1*s + (5*LQ*s)/4 - 35/2

Collecting terms

NI1 = collect(NI1,LI1)

gives

NI1 =
(s + 5/4)*LI1 + (5*LQ*s)/4 - 35/2

A similar string substitution

NQ = ...
subs(NQ,{laplace(I1(t),t,s), laplace(Q(t),t,s)}, {LI1,LQ})

yields

NQ(t) =
(4*LQ)/5 - (2*LI1)/5 + LQ*s - 1/(5*(s^2 + 1)) - 2

which, after collecting terms,

NQ = collect(NQ,LQ)

gives

NQ(t) =
(s + 4/5)*LQ - (2*LI1)/5 - 1/(5*(s^2 + 1)) - 2

Now, solving for LI1 and LQ

[LI1, LQ] = solve(NI1, NQ, LI1, LQ)

you obtain

LI1 = 
(5*(60*s^3 + 56*s^2 + 59*s + 56))/((s^2 + 1)*(20*s^2 + 51*s + 20))
 
LQ =
(40*s^3 + 190*s^2 + 44*s + 195)/((s^2 + 1)*(20*s^2 + 51*s + 20))

To recover I1 and Q, compute the inverse Laplace transform of LI1 and LQ. Inverting LI1

I1 = ilaplace(LI1, s, t)

produces

I1 =
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

Inverting LQ

Q = ilaplace(LQ, s, t)

yields

Q = 
(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

Now plot the current I1(t) and charge Q(t) in two different time domains, 0 ≤ t ≤ 10 and 5 ≤ t ≤ 25. The following statements generate the desired plots.

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

Note that the circuit's behavior, which appears to be exponential decay in the short term, turns out to be oscillatory in the long term. The apparent discrepancy arises because the circuit's behavior actually has two components: an exponential part that decays rapidly (the "transient" component) and an oscillatory part that persists (the "steady-state" component).

Was this topic helpful?