MAT2437 Matlab Worksheet Laplace Transforms PDF

Preview

Summary

Worksheet on LAPLACE TRANSFORMS...

Description

School of Engineering and Mathematics

MAT2236: Differential Equations Using MATLAB to calculate Laplace transforms and inverse Laplace transforms of functions In this worksheet we consider applications of the Laplace transform and its inverse. The MATLAB command for computing the Laplace transform of a given function is laplace and its inverse is ilaplace. Example 1: Compute the Laplace transform of each of the functions  (i) f (t )  e 4 t  (ii) f (t )  t 3 e 4t (iii) f (t )  u( t)  u( t  3) 0  t  2  sin t if (iv) f (t )    0 else

The function u denotes the unit step function. This function is also known as the Heaviside function and in MATLAB it is called by this name, so f (t )  u( t)  u( t  3)  heaviside(t )  heaviside(t  3) Note that the function in (iv) is piecewise defined and we can rewrite it making use of u. 0  t  2 sin t if f ( t)    ( u(t )  u(t  2 )) sin t  0 else We first define t as a symbolic variable, then define the function whose Laplace transform we wish to find as a function of t and finally apply the MATLAB command laplace. (i)

>> syms t >> f=exp(-4*t); >> laplace(f) ans =1/(s+4)

(ii)

>> f=t^3*exp(-4*t) >> laplace(f) ans = 6/(s+4)^4

(iii)

>> f =heaviside(t)-heaviside(t-3); >> laplace(f) ans =1/s-exp(-3*s)/s (iv)

>> f=sin(t)*(heaviside(t)-heaviside(t-2*pi)); >> laplace(f)

1

ans =1/(s^2+1)-exp(-2*s*pi)/(s^2+1) To find the inverse Laplace transform of a given function we use the command ilaplace. Example 2: Compute the inverse Laplace transform of each of the functions 1 (i) F (s)  s 1 (ii) F (s)  (s  5)2 s 3 (iii) F (s)  2 s  6s  5 s 7  (iv) F (s)  2 e 4s s  14 s  50

We first define s as a symbolic variable, then define the function whose inverse Laplace transform we wish to find as a function of s and finally apply the MATLAB command ilaplace. (i)

>> syms s >> f=1/s; >> ilaplace(f) ans = 1

(ii)

>> f=1/(s-5)^2; >> ilaplace(f) ans = t*exp(5*t)

(iii)

>> f=(s+3)/(s^2+6*s+5); >> ilaplace(f) ans = exp(-3*t)*cosh(2*t)

(iv)

>> f=(s-7)/(s^2-14*s+50)*exp(-4*s); >> ilaplace(f) ans = exp(7*t-28)*cos(t-4)*heaviside(t-4)

The results in (i) and (ii) are straight forward. For (iii) and (iv) this is not quite the case. For s 3 s3 and so the result in (iii) follows using a (iii) observe that F ( s)  2  s  6s  5 ( s  3) 2  4 shifting theorem and the Laplace transform of the function cosh. On the other hand we could have used a partial fraction decomposition first. With MATLAB this is done as follows: We define two vectors, one containing the coefficients of the polynomial in the numerator of the rational function starting with the coefficient of the highest power in s where missing powers get a coefficient of 0 and the other containing the coefficients for the denominator polynomial.

2

s3 , so the numerator is s  3 with corresponding vector b=[1 3] and s  6s  5 the denominator is s 2 6 s 5 . The vector corresponding to it is a=[1 6 5]. The MATLAB command for the decomposition is residue (use the MATLAB help file to get the complete information on residue). Its output are 3 vectors: the vector containing the residues, the vector containing the roots of the denominator polynomial, and a vector containing the coefficients of the direct polynomial term. The latter is 0 unless the degree of the numerator polynomial is greater than that of the denominator. We have F ( s) 

2

>> b=[1 3]; >> a=[1 6 5]; >> [r p k]=residue(b,a) r= 0.5000 0.5000 p= -5 -1 k= []

Now compute the terms in the partial fraction decomposition. >> f1=r(1)/(s-p(1)) f1 = 1/2/(s+5) >> f2=r(2)/(s-p(2)) f2 = 1/2/(s+1)

We can now write down the partial fraction decomposition of F: F ( s)  F1 ( s)  F2 ( s) where 1/ 2 1/ 2 F1 ( s)  and F2 ( s)  . Thus we could have also found that the inverse Laplace s 5 s 1 transform is given by the equivalent expression: >> ilaplace(f1)+ilaplace(f2) ans = 1/2*exp(-5*t)+1/2*exp(-t)

Note that cosh(2 t) 

1 2t e 3 t t ( e  e 2 t )  ( e  e 5t ) . 2 2

We now apply our knowledge of Laplace transforms to the solution of initial value problems.

3

Example 3. Use Laplace transform techniques to determine the solution of each of the following IVPs: (i) y  12 y   40 y  sin(2 t), y(0)  2, y(0)  0 y   2ty   4 y  4, y(0)  0, y(0)  0 (ii) y   16 y  cos 4 t   ( t 1), y(0)  0, y(0)  0 (iii)

legwork. Recall that L( y )  sY  y (0) d d L( y )  s 2Y  sy(0)  y (0) and that L (tf )   L ( f ), so that L( ty )   ( sY  y (0)) . ds ds

Here

we

need

to

first

do

some

and

Thus the Laplace transform of the first DE is 2 2 L( y   12 y   40 y)  s Y  sy(0)  y (0)  12( sY  y(0))  40Y  ( s  12 s  40)Y  2 s  24 and so we need to solve 2 ( s 2 12 s  40)Y  2 s  24  2 s 4 >> syms s >> Y=(2*s-24)/(s^2-12*s+40)+2/(s^2+4)/(s^2-12*s+40); >> y=ilaplace(Y) y = 155/78*exp(6*t)*cos(2*t)-311/52*exp(6*t)*sin(2*t)+1/78*cos(2*t)+1/52*sin(2*t)

For the second IVP we have L ( y   2ty   4y )  s 2Y  sy (0)  y (0)  2

d ( sY  y (0))  4Y  (s 2  6)Y  2sY  ds

and so we need to solve ( s 2  5)Y  sY 

4 s

or equivalently s 3 4 Y   (  )Y   2 s 2 s So this time we need to solve a differential equation: >> dsolve('DY-(s/2-3/s)*Y-4/s^2','s') ans = -8/s^3+exp(1/4*s^2)/s^3*C1

and then set the constant of integration equal to 0 and compute the inverse Laplace transform >> ilaplace(-8/s^3) ans = -4*t^2

In the third case the delta functional appears on the right hand side and the correspondin MATLAB function is dirac.. We have on the left hand side L( y  16 y)  s 2 Y  sy (0)  y(0)  16Y  ( s 2  16)Y We compute the Laplace transform of the right hand side >> syms t 4

>> laplace(dirac(t-pi)+dirac(t-2*pi)) ans = exp(-s*pi)+exp(-2*s*pi)

So we need to find the inverse Laplace transform of Y where ( s 2  16) Y  e  s  e 2 s : >> ilaplace((exp(-s*pi)+exp(-2*s*pi))/(s^2+16)) ans = (1/4*heaviside(t-pi)+1/4*heaviside(t-2*pi))*sin(4*t)

We lastly consider the convolution of two functions. Given f and g we can either compute the convolution of f and g straight from the definition or we can make use of the following property of Laplace transforms: L ( f * g )  L ( f ) L( g ) so that f * g  L1 ( L( f ) L( g )) . We consider two examples:  Example 4 Let f (t )  e t and g (t )  sin t . Determine f * g and plot its graph.

Using f * g  L1 ( L( f ) L( g )) we have >> ilaplace(laplace(exp(-1*t))*laplace(sin(t))) ans = 1/2*exp(-t)-1/2*cos(t)+1/2*sin(t)

Using the definition we have >> syms t >> int('sin(v)*exp(-(t-v))','v',0,t) ans = 1/2*exp(-t)-1/2*cos(t)+1/2*sin(t)

and a graph of the convolution can be obtained using the command >> ezplot('1/2*exp(-t)-1/2*cos(t)+1/2*sin(t)',[-5,10])

Example 4 Let f (t )  u (t )  u (t  2) and g (t )  (u( t) u(t 2 )) sin t . Determine the interval on which f * g is not identically 0, a formula for f * g and plot its graph.

5

We use the definition of the convolution to determine the set of values where f * g is nonzero. t

f  g (t )   f (t  v ) g (v )dv 0

The function value f (t  v ) is obtained from f (v ) though a reflection in the y-axis and a shift by t units in the positive t-direction. The functions f and g are non-zero on [0, 2] , so the function f * g will be non-zero on [0, 4 ] . We will verify this through a series of plots of f (t  v ) and g (v ) for t  , 0, , 2, 3 , 4 , 5  .

6

>> syms t >> laplace((heaviside(t)-heaviside(t-2*pi))*sin(t))

ans =1/(s^2+1)-exp(-2*s*pi)/(s^2+1) >> laplace(heaviside(t)-heaviside(t-2*pi))

ans =1/s-exp(-2*s*pi)/s >>h= laplace(heaviside(t)-heaviside(t-2*pi))*laplace((heaviside(t)-heaviside(t-2*pi))*sin(t))

ans = (1/s-exp(-2*s*pi)/s)*(1/(s^2+1)-exp(-2*s*pi)/(s^2+1)) >> ilaplace(h) ans = (2*cos(t)-2)*heaviside(t-2*pi)+heaviside(t-4*pi)*(-cos(t)+1)-cos(t)+1

The graph of this function is shown in the figure below.

Note that this result appears to have a cosine term without restriction on the domain. This arises because the program cannot distinguish between the Laplace transform of cos t and that of u(t ) cos t . So the result should be ans = (2*cos(t)-2)*heaviside(t-2*pi)+heaviside(t-4*pi)*(-cos(t)+1)+ heaviside(t)*(-cos(t)+1)

The graph is given by

7

We lastly consider the solution of an integro-differential equation. t

Example 5: Determine the solution of

dy  y   y( u) du 1, y(0)  0 . dt 0

t

Y ( s) as we can regard the given integral as the s 0 convolution 1* f . So once the Laplace transform has been carried out we have Y 1 sY  y(0)  Y   s s which is equivalent to s 2Y  sY  Y  1 , that is (s 2  s  1)Y  1 and so 1 Y 2 . s  s 1 So we have The Laplace transform of

 y (u )du

is given by

>> syms s >> ilaplace(1/(s^2+s+1)) ans = 2/3*3^(1/2)*exp(-1/2*t)*sin(1/2*3^(1/2)*t)

8...