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 { d² x/dt²} = s² 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)│^∞ ₀ + 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 d² x/dt² + k x =0 . (5)

Let m = 0.100 kg , k = 3 kg/s² 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 d² x/dt²} = .100( s² X(s) –s x(0)– dx(0)/dt )

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

The transformed eq(5) is

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

X(s) = (.1*.15) s /(.1*s² +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

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

Section 7. The second compartment…

The Laplace transform for C₂(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 =

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/T² , y(0) =y0=2 , T is a constant. (11)

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

sY(s) –2 – kY(s) = (1/(s*T²))

Y(s)= 2/(s-k) + (1/T²) *( 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/(kT²))*(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

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

L (d²C/dt² – a C(t) ) = s²c(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 d²C/dt² + 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 ω = a^1/2 .

Problem 9.

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

L (d²C/dt² – a C(t) ) = s²c(s) –s C(0)-dC(0)/dt – a c(s) =0

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

c(s) = 1/(s² –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

dY₁/dt = - k Y₂ + R (15-a)

dY₂/dt = k Y₁ . (15-b)

IC Y₁ (0) = Y₂ (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

Y₁(t_n ) = Y₁ (t_n-1) + ( ∆t ){ - k Y₂ (t_n-1 )+ R }

Y₂(t_n ) = Y₂ (t_n-1) + ( ∆t ){ k Y₁ (t_n-1 ) }

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

∆t<< 5 minutes.

According to the textbook , ref .1 ,

Y₁(s) = R/(s^2+k^2)

Using Matlab we obtain y₁(t) as follows ;

MATLAB

>> syms s t ;

k=0.2 ; R=0.1;

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

y₁(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

d²Y₁ /dt² = - k dY₂/dt + 0 (16-a)

d² Y₂/dt² = k dY₁ /dt (16-b)

Inserting 15-b in 16-1 gives

d²Y₁ /dt² = - k dY₂/dt + which is the equation for a harmonic oscillator.

The equation for Y₂ is

d²Y₂ /dt² = - k² Y₂ + kR .

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

Problem 12.

Solve by Laplace transform

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

L{ d ² x/dt² }= s² X(s)-s x(0) –dx(0)/dt = s² 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)/( s² +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