Lecture 4.

 

The Laplace Transform

by Reinaldo Baretti Machín

 

www.geocities.com/serienumerica4

www.geocities.com/serienumerica3

www.geocities.com/serienumerica2

www.geocities.com/serienumerica

 

For questions and suggestions write to:

reibaretti2004@yahoo.com

 

 

 

 

 

 

Reference :

 

1. Mathematical Techniques for Biology and Medicine (Dover Publications- Paperback)

by William Simon ; http://www.amazon.com/Mathematical-Techniques-Biology-Medicine-William/dp/0486652475/ref=sr_1_1?ie=UTF8&s=books&qid=1223731765&sr=8-1

 

2. http://geocities.com/serienumerica4/Lecture1Biolologymedicine.doc

3. Mathematical Biology: I. An Introduction (Interdisciplinary Applied Mathematics) by James D. Murray

 

4. http://www.geocities.com/serienumerica4/Lecture3Biolologymedicine.doc

 

5. Schaum's Outline of Advanced Mathematics for Engineers and Scientists by Murray R Spiegel

6. Schaum's Outline of Laplace Transforms by Murray R Spiegel

 

7. http://www.mathworks.com/access/helpdesk/help/toolbox/symbolic/index.html?/access/helpdesk/help/toolbox/symbolic/laplace.html&http://search.yahoo.com/search?p=MATLAB+laplace+transfom&ei=utf-8&fr=b1ie7

 

 

Examples and problems are taken from:

Chapter 4 of  Reference 1. The problem numbers refer to the textbook listing.

 

 

One method of solving linear DE with constants coefficients employs the Laplace transform. There is software to implement the Laplace transform method (see ref . 7   ).

 

Let the independent variable be time –t and the dependent variable be x(t).

 

The transfom of x(t) is called X(s) and is given by

 

 

X(s) =   L {x(t) } ≡ ∫ x(t) exp(-st) dt           ,     0  ≤  t   ≤ ∞  .   (1)

 

X(s) ~ meter-second

 

The dimension of    s ~ 1/time ~ frequency  . It can be shown that the transfoms of the first and second derivatives  are

 

L { dx/dt} = s X(s) – x(0)     ~ meters                                          (2)

 

                                             

L { d2 x/dt2} = s2 X(s) –s x(0)– dx(0)/dt ~ meter/sec     .              (3)  

 

Proof of (2). Using integration by parts

 

  ∫ (dx/dt) exp(-st) dt  =  x(t) exp(-st)│ 0   + s ∫ x(t) exp(-st) dt 

 

                                   = s X(s)  - x(0)                             .   (4)                                        

 

 

Example find the Laplace transform using MATLAB software.

 

MATLAB

 

syms t a s;

f = exp(-a*t);

laplace(f,s)

ans = 1/(s+a)

 

The inverse Laplace command is given by

 

ilaplace(1/(s+a))   , which returns ,

ans =  exp(-a*t)

 

 

Example find the inverse Laplace transform using MATLAB software.

 

 

 

The differential equation along with the initial conditions is transformed into an algebraic expression which can be solved for X(s).

 

We give next a simple example of the harmonic oscillator.

 

 

A harmonic oscillator is governed by the  DE

 

m d2 x/dt2  + k x =0  .                             (5)

 

Let m = 0.100  kg ,  k = 3 kg/s2  with initial conditions x(0)=0.15m ,

dx(0)/dt = 0   ,  find x(t).

 

L{k x(t)}= k X(s)  = 3 X(s)                              

 

L { m d2 x/dt2} = .100( s2 X(s) –s x(0)– dx(0)/dt )

                         =  .100 [ s2 X(s) –.15 s ]        (6)

 

The transformed eq(5) is

0.100 (s2 X(s) –.15 s)   + 3 X(s) =0 , therefore the Laplace transform of the unknown x(t) is  ,

 

X(s) = (.1*.15) s /(.1*s2 +3)                                                           .   (7)

 

To get the inverse transform use the command ilaplace from MATLAB.We have

 

syms s  ;

x=ilaplace( (.1*.15*s)/(.1*s^2+3) )

x = (3/20)*cos(30^(1/2)*t)

 

We proceed to plot it

 

 

Fig 4.1  x(t)  vs t for the simple  harmonic oscillator.

 

Section 6.

The one compartment DE was shown in (2.9 of the textbook)

 to be

 dC/dt = R/V  -(F/V) C(t)       with IC C(0) .                   (8)

 

The Laplace transform is  g = L{ C(t)} 

 

g(s)= C(0) /( s+ (F/V) + (R/V) /( s ( s+(F/V) ) )  . Using MTALAB  we obtain the solution C(t) invoking the command ilaplace.

MATALAB

syms g s C0 F V R;

g = C0/(s+F/V) +(R/V)/(s*(s+(F/V))) ;

ct=ilaplace(g)

 

ct =

 

R/F+exp(-F*t/V)*(C0*F-R)/F

Or

C(t) = (R/F)+ (C(0) –R/F)*exp(-F*t/V)  .

 

Section 7.  The second compartment…

 

The Laplace transform for C2(t) is ,

 

g(s) = (C10*F/V2)*(1/(s+(F/V1)))* (1/(s+(F/V2)))  .

 

Using MATLAB

 syms g s C10 F V1 V2 ;

g = (C10*F/V2)*(1/(s+F/V1))*(1/(s+F/V2));

ct=ilaplace(g)

 

ct =

 

C10*V1/(-V1+V2)*(-exp(-F*t/V1)+exp(-F*t/V2))

 

 

Problems of Chapter 4.

 

 

Problem 4.    Solve the simultaneous system of linear equations ;

 

2x + 4y +6z = 28 

3x - 3y  +3z = 6

3x  + 3y +6 z = 27                                                                   (9)

 

This is of the form A*X=B . Using MATLAB we have

MATLAB

 

A=[2,4,6;3,-3,3;...

 3,3,6 ];

B=[28;6;27];

x=A\B

 

x =

 

     1

     2

     3

i.e  x=1   , y=2  , z=3  .

 

 

 

Problem 6.

 Solve by Laplace transform

 

                    dy/dt –ky = 0     , y(0) =y0=3 .                              (10)

 

L (dy/dt –ky ) = sY(s) –y0 – kY(s)    , where L{y(t)} = Y(s)  .

 

Y(s) = y0/(s-k)      .

MATLAB

syms s  k;

 y0=3;

 yt=ilaplace(y0/(s-k))

 

 yt =  3*exp(k*t)

 

****

 

Problem 7.

 Solve by Laplace transform

 

                    dy/dt –ky = 1/T2     , y(0) =y0=2 , T is a constant.      (11)

 

L (dy/dt –ky ) =L {1/T2}     , where L{y(t)} = Y(s)  .

sY(s) –2 – kY(s) = (1/(s*T2))

Y(s)= 2/(s-k) + (1/T2) *( 1/(s*(s-k))

 

y(t) = L-1 (Y(s) ) 

 

MATLAB

 syms s T  k;

 y0=2;

 yt=ilaplace(y0/(s-k))+...

     (1/T^2)*ilaplace(1/(s*(s-k)))

 

yt =

 

2*exp(k*t)-1/T^2/k*(1-exp(k*t))      ,

i.e.

y(t) = 2 exp(kt) – (1/(kT2))*(1-exp(kt) )

 

******

Problem 8. Solve by Laplace transform

 

dC/dt +(L/V)C = Aexp(-kt) +(R/V)  , C(0)=C0=0  .                      (12)

 

Let L{C(t)} = c (s)

L { dC/dt +(L/V)C }  = L {   Aexp(-kt) +(R/V)  }

 

sc(s)-C0 +(L/V)c(s) = A/(s+k)    +    (R/V)*(1/s)

 

 c(s) = A/[(s+L/V)*(s+k)]  +(R/V)*(1/(s(s+L/V)) )    

 

 MATALAB

syms s A L V R  k;

 y1=A/((s+L/V)*(s+k));

 y2=(R/V)*(1/(s*(s+L/V)));

 yt=ilaplace(y1) + ilaplace(y2)

 

yt =A/(k-L/V)*(exp(-L/V*t)-exp(-k*t)) + R/L*(1-exp(-L/V*t)) or

 

C(t) =(A/(k-L/V))* {exp(-Lt/V) –exp(-kt) }

 

        + (R/L)*{ 1-exp(-Lt/V) }

 

Problem 8. Solve by Laplace transform

 

d2C/dt2 – a C(t) = 0    , IC  C(0)=1   , dC(0)/dt =0.               (13)

 

L (d2C/dt2 – a C(t) ) = s2c(s) –s C(0)-dC(0)/dt – a c(s) =0

 

c(s) = s /(s^2 –a)      .

 

C(t) = L-1 { c(s) }                              .

 

MATLAB

>> syms s a t;

>> cs=s/(s^2-a);

>> ilaplace(cs)

 

ans = cosh(a^(1/2)*t)

 

Changing  the DE to   d2C/dt2 + a C(t) = 0   with the same IC leads to

 

c(s) = s/(s^2+a).

Then

cs=s/(s^2+a);

>> ilaplace(cs)

 

ans =  cos(a^(1/2)*t)  . It is the harmonic oscillator solution with angular frequency  ω = a1/2 .

 

 

Problem 9.

Repeat DE d2C/dt2 – a C(t) = 0    with  IC  C(0)=0   , dC(0)/dt =1.  (14)

 

L (d2C/dt2 – a C(t) ) = s2c(s) –s C(0)-dC(0)/dt – a c(s) =0

 

                              = s2c(s) –0 - 1 – a c(s) = 0.

c(s) = 1/(s2 –a)

MATLAB

cs=1/(s^2-a);

>> ilaplace(cs)

 

ans =  1/a^(1/2)*sinh(a^(1/2)*t)

Problem 11.    Solve the simultaneous system of DE

 

                                  dY1/dt = - k Y2 + R                     (15-a)

 

                                  dY2/dt  =  k Y1  .                         (15-b)

 

IC   Y1 (0) = Y2 (0)= 0 . Let  t~minutes  , R = 0.1 moles/(liter-minute)

 

k = 1/(5minutes)  = 0.2 minute-1 .

 

We will apply the finite difference method.

 

The DE give rise to

       

Y1(tn )  =   Y1 (tn-1) + ( ∆t ){ - k Y2 (tn-1 )+ R  }  

Y2(tn )  =   Y2 (tn-1) + ( ∆t ){  k Y1 (tn-1 )  }  

 

The time sacle of the problem is Tscale= 1/k = 5 miniutes , hence one selects

 

∆t<< 5 minutes.

According to the textbook , ref .1 ,

 

Y1(s) = R/(s^2+k^2) 

 

Using Matlab we obtain  y1(t) as follows ;

 

MATLAB

>>  syms s t ;

k=0.2 ; R=0.1;

ilaplace(R/(s^2+k^2))

 

 y1(t) =  ans =  1/2*sin(1/5*t) .

 

 

Fig 4.2 Concentrations Y1 , Y2.

 

 

 

 FORTRAN code for Problem 11

c prob chapter 4 .11  , Simon math tech biol & medicine

      real k

      data k , r /0.2,0.1/

      y(t)=(1./2.)*sin((1./5.)*t)

      tscale=1./k

      Y10=0.

      Y20=0.

      dt=tscale/200.

      tf=5.*tscale

      nstep=int(tf/dt)

      kp=int(float(nstep)/60.)

      kount=kp

      print 100 , 0. , y10,y(0.), y20

      do 10 i=1,nstep

      t=dt*float(i)

      y11=y10+dt*( - k*y20 + R  )

      y21=y20+dt*k*y10

      if(i.eq.kount)then

      print 100 , t ,y11,y(t),y21

      kount=kount+kp

      endif

      y10=y11

      y20=y21

10    continue

100   format(1x,'t,y1,Y1(exact),y2=',4(3x,e10.3))

      stop

      end

     

Can we see by manipulating the DE that both solutions are oscillators ?

From (15) we get , taking another derivative

d2Y1 /dt2 = - k  dY2/dt + 0                           (16-a)

 

 d2 Y2/dt2  =  k dY1 /dt                                    (16-b)

 

Inserting 15-b in 16-1 gives

 

d2Y1 /dt2 = - k  dY2/dt + which is the equation for a harmonic oscillator.

The equation for Y2 is

 

d2Y2 /dt2 = - k2  Y2 + kR       .

 

The procedure suggests that a system of coupled  first order DE may be substituted by the same number of uncoupled 2nd  orde DE. 

 

 

Problem 12.

Solve   by Laplace transform

 

 d 2 x/dt2  + dx/dt + 0.09 x= 0 , IC  x(0)=1 , dx(0)/dt=0.  (17)

 

 L{ d 2 x/dt2 }= s2 X(s)-s x(0) –dx(0)/dt =   s2 X(s)-s   

 

L{dx/dt} = sX(s)-x(0)                            = s X(s) -1

 

L{.09 x(t)   }   = 0.09 X(s)                      

We have

X(s) = (s+1)/( s2 +s+.09)

 

MATLAB

>> syms s t;

>> ilaplace((s+1)/(s^2+s+.09))

 

ans =  9/8*exp(-1/10*t)-1/8*exp(-9/10*t) , or,

 

x(t) = (9/8)exp(t/10) –(1/8) exp(9t/10)    .

 

END OF LECTURE