## 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\left(s\right)=\underset{0}{\overset{\infty }{\int }}f\left(t\right){e}^{-ts}dt.$`

### 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

• 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\left(t\right)-\frac{1}{C}Q\left(t\right)\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, I0 and Q0, 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 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/4 F, L = 1.6 H, I1(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]```

## Support

#### Mathematical Modeling with Symbolic Math Toolbox

Get examples and videos