Solving ODEs with the Laplace Transform in Matlab

This approach works only for

• linear differential equations with constant coefficients
• right-hand side functions which are sums and products of
  • polynomials
  • exponential functions
  • sine and cosine functions
  • Heaviside (step) functions
  • Dirac (impulse) ``functions''
• initial conditions given at t = 0

The main advantage is that we can handle right-hand side functions which are piecewise defined, and which contain Dirac impulse ``functions''.

You must first save the file Heaviside.m in your directory. (This fixes a bug in Matlab.)

Simple example

Consider the initial value problem

y'' + 3 y' + 2 y = e^-t , y(0) = 4 , y'(0) = 5

Define the necessary symbolic variables:

syms s t Y

Define the right-hand side function and find its Laplace transform:

f = 'exp(-t)'
F = laplace(f,t,s)

Find the Laplace transform of y'(t) : Y₁ = s Y - y(0)

Y1 = s*Y - 4

Find the Laplace transform of y''(t) : Y₂ = s Y₁ - y'(0)

Y2 = s*Y1 - 5

Set the Laplace transform of the left hand side minus the right hand side to zero and solve for Y:

Sol = solve(Y2 + 3*Y1 + 2*Y - F, Y)

Find the inverse Laplace transform of the solution:

sol = ilaplace(Sol,s,t)

Example with piecewise defined right-hand side function

Consider the initial value problem

y'' + 3 y' + 2 y = f(t) , y(0) = 2 , y'(0) = 3

with the right-hand side function

f(t) = 1 for t<3
f(t) = t -2 for 3<t<6
f(t) = 2 for t>6

Define the necessary symbolic variables:

syms s t Y

As the right-hand side function is piecewise defined, rewrite it in terms of the Heaviside function: We have f(t) = f₁(t) + (f₂(t)-f₁(t))H(t-t₁) + (f₃(t)-f₂(t))H(t-t₂)

f = '1 + ((t-2)-1)*Heaviside(t-3) + (2-(t-2))*Heaviside(t-6)'
ezplot(f,[0,10])

Find the Laplace transform of the right hand side function f(t):

F = laplace(f,t,s)

Find the Laplace transform of y'(t) : Y₁ = s Y - y(0)

Y1 = s*Y - 2

Find the Laplace transform of y''(t) : Y₂ = s Y₁ - y'(0)

Y2 = s*Y1 - 3

Set the Laplace transform of the left hand side minus the right hand side to zero and solve for Y:

Sol = solve(Y2 + 3*Y1 + 2*Y - F, Y)

Find the inverse Laplace transform of the solution:

sol = ilaplace(Sol,s,t)

Plot the solution:

ezplot(sol,[0,10])

Example with Dirac ``function''

Consider the initial value problem

y'' + 2 y' + 10 y = 1 + 5 (t-5) , y(0) = 1 , y'(0) = 2

Define the necessary symbolic variables:

syms s t Y

Define the right hand side function:

f = '1 + 5*Dirac(t-5)'

Find the Laplace transform of the right hand side function:

F = laplace(f,t,s)

Find the Laplace transform of y'(t) : Y₁ = s Y - y(0)

Y1 = s*Y - 1

Find the Laplace transform of y''(t) : Y₂ = s Y₁ - y'(0)

Y2 = s*Y1 - 2

Set the Laplace transform of the left hand side minus the right hand side to zero and solve for Y:

Sol = solve(Y2 + 2*Y1 + 10*Y - F, Y)

Find the inverse Laplace transform of the solution:

sol = ilaplace(Sol,s,t)

Plot the solution:

ezplot(sol,[0,10])