Lecture 6.  

 

Numerical Methods

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 6 of  Reference 1. The problem numbers refer to the textbook listing.

 

 

We have already employed a number of tools  from numerical methods in solving DE.

 

Finite difference is perhaps the most simple method to get by.

 

Another method of choice to solve DE is Taylor series.

 

If a function and its derivatives are known at a point x0  one can estimate the value of the function at x0 + ∆x , using an approximation to the Taylor series. The Taylor series is an infinite sum

 

 

y(x0 + ∆x ) = y(x0)  +  ∆x  (dy/dx)0 + (1/2) (∆x)2 (d2 y/dx2)

 

                                  + ..... (1/n!) (∆x)n (dn y/dxn)              .     (6.1)

 

In practice we can do with a few terms.

 

Example 1. (equation 6.4 from ref 1) Consider

 

dC2/dt = (F/V2)C1(0) exp(-(F/V1)t)   -(F/V2) C2(t)                    (6.2)

 

with the parameter values V1= 10 L  , V2 =20 L  , F= 1 L/sec and IC

C1(0)=.3 mole/L  , C2(0)=0.

 

Substitution in (6.2) gives

  dC2/dt   =  (1/20)(.3) exp(-.1 t) –(1/20) C2(t)

    

    dC2/dt         = .015 exp(-.1t) -.05 C2(t)                     (6.3)

 

The second derivative is

 

 d2 C2/dt 2  = -1.5E-3 *exp(-.1 t) -.05 dC2/dt  

 

                 =   -1.5E-3 *exp(-.1 t) 

                       -.05 { .015 exp(-.1t) -.05 C2(t)  }         (6.4)

 

The truncated Taylor series to order (∆t)2 is

 

 

C2(t0 +∆t) = C2(t0) + ∆t{ .015 exp(-.1t0) -.05 C2(t0)  }

 

 

   +(1/2)( ∆t)2 { -1.5E-3 *exp(-.1 t0)  -.05 [ .015 exp(-.1t0) -.05 C2(t0) ] }

                                                                                      

                                                                                                   (6.5)

 

 

  In this equation there are two time scales defined by

tscale1= V1/F =  10 sec           ;  tscale2 = V2/F= 20 sec.

 

So one must choose ∆t << 10 sec to apply the Taylor series.

 

 

 

 

Fig 6.1  Taylor series solution of  equation (6.3).

 

FORTRAN code

c solution by Taylor series to 2nd order of

c dC2/dt    = .015 exp(-.1t) -.05 C2(t)

      data c2zero/0./

      tfinal=40.

      dt=10./100.

      nstep=int(tfinal/dt)

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

      kount=kp

      print 100 ,0., c2zero

      do 10 i=1,nstep

      t=dt*float(i)

      c2one=c2zero+dt*(.015*exp(-.1*(t-dt)) -.05*c2zero) +

     $.5*dt**2*(-1.5E-3*exp(-.1*(t-dt))-.05*(.015*exp(-.1*(t-dt)) -

     $  .05*c2zero ) )

      if(i.eq.kount)then

      print 100 ,t, c2one

      kount=kount+kp

      endif

      c2zero=c2one

10    continue

100   format(1x,'t,C2=',2(3x,e10.3))

      stop

      end

     

     

 

 

 

Example 2. Solve by Taylor series the 2nd order DE

 

d2y/dt2 + b dy/dt + c y =0 ,                                           (6.6)

 

where b= 1 /s  , c= 1/s2  and IC

 

y(0)=1 , dy(0)/dt=0. This is the equation of a damped oscillator.

 

The solution to  second order is,

 

 y(t0 + ∆t ) = y(t0)  +  ∆t  (dy/dt)0 + (1/2) (∆t)2 (d2 y/dt2) .  (6.7)

 

One also needs a “Taylor series” for the first derivative which is

 

 

dy( t0 +∆t )/dt =  dy( t0)/dt  +  ∆t  [d2 y(t0) /dt2]                    (6.8)

 

This value (6.7) is needed because the equation is of the second order and the previous  y(t0)  and  (dy/dt)0 are needed.

 

 

Inserting  d2y/dt2 = - dy/dt -  y   in (6.7) and (6.8) the two truncated Taylor series become

 

y(t0 + ∆t ) = y(t0)  +  ∆t  (dy/dt)0 + (1/2) (∆t)2 (- (dy/dt)0 -  y0  )   (6.9)

and

 

dy( t0 +∆t )/dt =  (dy/dt)0  +  ∆t  ( - (dy/dt)0 -  y0  )   .          (6.10)       

 

 

Fig. 6.2 Solution of  d2y/dt2 + dy/dt +  y =0 , IC y(0)=1, dy(0)/dt=0.

 

FORTRAN code

c taylor series fro 2nd order DE

      data y0,dydt0/1.,0./

      dydt2(dydt,y)= -dydt-y

      dt=1./100.

      tfinal=8.

      nstep=int(tfinal/dt)

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

      kount=kp

      print 100 ,0.,y0

      do 10 i=1,nstep

      t=dt*float(i)

c Taylor series up to dt**2 for the position

      y1=y0 +dt*dydt0 +.5*dt**2*dydt2(dydt0,y0)

c Taylor series up to dt for the velocity

      dydt1=dydt0+dt*dydt2(dydt0,y0)

      if(i.eq.kount)then

      print 100 ,t, y1

      kount=kount+kp

      endif

      y0=y1

      dydt0=dydt1

10    continue

100   format(1x,'t,y=',2(3x,e10.3))

      stop

      end

 

 

3. Numerical solution of a circular chain reaction with changing rate constant

Fi 6.3 Simplified representation of a circular reaction.

 

The DE are

 

dA/dt = - K1 A + K3 C                                            (6.11)

dB/dt = - K2 B + K1 A

dC/dt = - K3 A + K2 B

 

with specified initial conditions A(0) , B(0)  and  C(0) ~ moles

 

In chapter 2 the circular reaction DE equations were solved with constant Ki.

 

Now we take one of the rate constants to be variable ,

 

say  K1 (t) = K10 + m t     .

 

The finite difference solutions to (6.11) are

 

 

 

A(tn ) = A( t n-1 ) + ∆t { - K1(tn-1 ) A(tn-1) + K3 C (tn-1)  }          (6.12)                                                                

 

 

B(tn ) = B( t n-1 ) + ∆t { - K2 B(tn-1) + K1(tn-1 ) A(tn-1)  }                          

 

 

C(tn ) = C( t n-1 ) + ∆t { - K3 C(tn-1) + K2(t n-1 ) B(tn-1)  }      

 

Let    A(0) = 300 micro mol, B(0)= 60 micro moles  and  C(0)=30 micro moles. 

Let the  rate constants be ;

K1  = 0.1 + [(10-.1)/2]  t       for t ≤ 2 s ,

K1 = .1+9.9=10.0/s         for  t  ≥    2 s.

 

 

  while   K2 = 2/s       , K3 = 1.0/s  .

 

The slope m of K1 is  (9.9/2) /s2 = 4.95/s2        

 

The results are plotted in fig 6.4

 

Fig 6.4 Concentrations for the circular reaction with variable K1.

 

FORTRAN code

c Problem from Math Tech for BIOL & MED page 119 Simon

      real k1,k2,k3

      data k2,k3/2.,1./

      data a0,b0,c0 /300.,60.,30./

      tsmall=amin1(1./10.,1./k2,1./k3)

      tf=4.

      dt=tsmall/100.

      nstep=int(tf/dt)

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

      kount=kp

      print 100 ,0., a0,b0,c0

      do 10 n=1,nstep

      t=dt*float(n)

      a1=a0+dt*(-k1(t-dt)*a0+k3*c0)

      b1=b0+dt*(-k2*b0 +k1(t-dt)*a0)

      c1=c0+dt*(-k3*c0+k2*b0)

      if(n.eq.kount)then

      print 100 , t,a1,b1,c1

      kount=kount+kp

      endif

      a0=a1

      b0=b1

      c0=c1

10    continue

100   format('t(s),A(micromol),B,C=',4(3x,e10.3))

      stop

      end

 

 

 

      function k1(t)

      real k1

      k1=0.1+ (9.9/2.)*t

      return

      end

 

 

 

Section  4.

 

Consider a first order pair of coupled DE for two metabolic substances A and B given by

 

dA/dt = (R1 – a B) –  R1’                                             (  6.13  )

 

dB/dt = R2  - ( R2’ – b A )                         .                    ( 6.14   )

 

A and B represent the quantity of moles above or below the equilibrium quantities A(equilibrium) and  B(equilibrium ). Thus both A and B may oscillate between positive and negative values.

 

The  figure serves to illustrate the kinetics involved.

 

Fig. 6.5 .    If B > 0  the term –aB inhibits the formation of A. On the contrary if B < 0, the term –aB augments the formation of A.

Like wise a positive A inhibits the loss of B and a negative value of A increases the loss of B.

 

 

 

The dimensions could be  the following.

A, B ~ moles , R ~moles/time , a~1/time , b~ 1/time.

The rate constans R1 , R1’ ,R2, R2’are independent of the concentrations A or B .

 

Modified example from page Ref 1, 121.

 

Take the set of values

A(0)= 300 micro moles           a= 1/s

B(0)= 200 μ mol                   b= 10/s

 

R1= 200 μ mol/s                   R1’= 100 μ mol/s                  

R2= 20 μ mol/s                     R2’ = 10  μ mol/s                  

 

We convert  (6.13 ) – (6.14) into the finite difference expressions

 

 

 

A(tn ) = A( t n-1 ) + ∆t {  R1- aB(tn-1) – R1’ }      (6.15)

and

 

dB/dt = R2  - ( R2’ – b A )                        

 

B(tn ) = B( t n-1 ) + ∆t {  R2 – R2’ + bA(tn-1)  }   

 

There are two time scales in this DE ,

tscale1 =1/a= 1.0 s   , tscale2=1/b =0.1 s .

One has to take    ∆t <<  tscale2.

 

 

Fig 6.6-  Plot of A and B in a reaction with mutual feedback. Oscillations with a 2.0 second period.

 

FORTRAN code

c Simon section 4 page 121

      data r1, r1prime ,r2,r2prime/200.,100.,20.,10. /

      data smalla,smallb/1.,10. /

      data a0,b0/300.,200. /

      tsmall=amin1(1./smalla,1./smallb)

c      tfinal=10.*amax1(1./smalla,1./smallb)

      tfinal=4.

      dt=tsmall/100.

      nstep=int(tfinal/dt)

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

      kount=kp

      print 100 ,0.,a0,b0

      do 10 i=1,nstep

      t=dt*float(i)

      p= r1-smalla*b0

      q= r2prime-smallb*a0

c      if(p.lt.0.)p=0.

c      if(q.lt.0.)q=0.

      a1=a0+dt*(p-r1prime)

      b1=b0+dt*(r2-q)

c      if(a1.lt.0.)a1=0.

c      if(b1.lt.0.)b1=0.

      if(i.eq.kount)then

      print 100 ,t, a1,b1

      kount=kount+kp

      endif

      a0=a1

      b0=b1

10    continue

100   format('t(sec),A,B=',3(3x,e10.3))

      stop

      end

     

 

Problem 1 , page 124.

 

 Solve numerically

 

 dC/dt = a t C(t)     , a =0.1 /s2  , C(0)=2.

 

tscale = 1/ a1/2 = 3.16 s

 

The finite difference solution is

 

C(tn) = C( tn-1 )  + ∆t { a tn-1 C(tn-1 )   }  .

 

MATLAB code

% Simon page 124

a=.1; c(1)=2. ; tscale=1/a^(1/2) ;

dt =tscale/100 ; tfinal=3*tscale ;

  

k=1;

 for  t=[dt:dt:tfinal];

k=k+1;

 c(k)=c(k-1)+dt*a*(t-dt)*c(k-1) ;

end

 

t=[0:dt:tfinal];

plot(t,c),xlabel ('time '), ylabel('C(t)')

 

 

Fig 6.7  Solution of dC/dt = a t C(t)    .

 

 

Problem 3.

Solve the following DE numerically.

 

 dy/dt = 100 + y      , IC   y(0)=0.

 The analytical answer is readily found

 

∫ dy/(100+y) = ∫ dt + C

 

ln ( 100+y ) = t+c

100+ y = exp ( c ) exp( t)

 

y= 100- exp(c) exp(t)  .

Applying the IC  y(0)=0 = 100 – exp(c)  , gives exp(c) =100. finally

y(t) = 100(1-exp(t) ).

 

To solve it numerically we write the finite difference solution

 

 y(tn) = y(tn-1) +∆t { 100 + a y(tn-1 ) }.

 

We have introduced a “factor a =1/s ” in order to define the time scale,

tscale= 1/a = 1.0 s

 

MATLAB code

% Simon Prob 3 -page 125

a=1; y(1)=0. ; tscale=1/a ;

dt =tscale/100 ; tfinal=3*tscale ;

  

k=1;

 for  t=[dt:dt:tfinal];

k=k+1;

 y(k)=y(k-1)+dt*(100+y(k-1)) ;

end

 

t=[0:dt:tfinal];

plot(t,y),xlabel ('time '), ylabel('Y(t)')

 

 

Fig 6.8 Solution of dy/dt = 100 + y      , IC   y(0)=0 .

 

 

 

 

Problem 4. The logistic equation  (LE).

 

                                      \[ 
\frac{dN}{dt} = r N \frac{(K - N)}{K} 
\]

 

 

Let N(t) be a population at time t  , r a rate of growth ~ 1/time , K some

 

population “limit”.  

If     K>>N     , dN/dt = r N , and the population grows as N=N(0)exp(rt).

 

This answer is unrealistic in the long run. The logistic equation sets a limit to growth.

 

The logistic equation is non linear

                 

  dN(t)/dt = r N(t) – (r/K) N(t)2 .

 

We may view the DE as containing two time scales

first  tscale1 = 1/r   but a larger time scale is represented by

 

tscale2=K/(rN) = (K/N)tscale1

 

Let N(0)=1.E 9  K=6.5 E9 and  r=.02/year

 

 

We solve the LE using the following code with two different values of the limiting population.

 

MATLAB code

 

% Simon Prob 4 -page 125

r=.02; N(1)=1.e9 ; K=6.5E9 ; tscale=1/r ;

dt =tscale/100 ; tfinal=6*tscale ;

  

k=1;

 for  t=[dt:dt:tfinal];

k=k+1;

 N(k)=N(k-1)+dt*(r*N(k-1)-(r/K)*N(k-1)^2) ;

end

 

t=[0:dt:tfinal];

plot(t,N),xlabel ('t~years'), ylabel('Population')

 

 

Figure  6.9  Population growth with K= 3 E9.

 

 

Figure 6.10  . Population growth with K= 6 E9.

 

 

 

 

Problem 5 . Van der Pol equation (VdP).

 

Solve the non linear equation

 

 d2 y/dt2λ (1- y2)(dy/dt) +  y = 0          .

 

with y adimensional  , λ = 0.4 /s  nad IC y(0)=0.1 , dy(0)/dt = 0 .

 

To know what to expect , recall first  the  equation of the damped oscillator

 

d2 y/dt2 +b (dy/dt) +  c y = 0      .

 

When comparing with VdP  we see that c=1 and b is variable like

   

    b ~   –λ (1- y2).

 

In the damped oscillator the solution is governed by terms

 

   y ~  exp(-bt/2)exp( (+/-)ω t)    , where  ω = (1/2) ( b2 -4 c) 1/2  .

 

It may happen that  (1-y2) > 0  , therefore  b is negative. It would tend  increase the amplitude  y through the first factor

 

exp(-bt/2) , if at the same time (b2 -4 c ) > 0 there are no imaginary

 

exponentials and therefore no oscillatory behavior.

 

If b < 0 and  (b2 -4 c ) < 0 there can be oscillatory motion.

 

Another possibility is that at an instant (1-y2) < 0  , therefore b is positive

and the amplitude y would decrease if there is an oscillatory motion.

 

But if y decreases sufficiently it send the solution to the first instance where

 

(1-y2) > 0  .

 

Fig 6.11. An oscillation with increasing amplitude is produced.

 

 

 

 

Fig 6.12. The same lambda but different initial conditions y(0)=1.5 , dy/dt=0

 

 

MATLAB code for problem 5

% Simon Prob 5 -page 126

lambda=.4; y(1)=.1 ; y(2)=.1  ; tscale=1/lambda ;

dt =tscale/400 ; tfinal=8*tscale ;

 

  

k=1;

 for  t=[2*dt:dt:tfinal];

k=k+1;

dydt=(y(k)-y(k-1))/dt;

y(k+1)=2*y(k)-y(k-1) + dt^2*(+lambda*(1-y(k)^2)*dydt...

    -y(k));

end

t=[0:dt:tfinal];

plot(t,y),xlabel('t~ sec'), ylabel('y ')

END OF LECTURE 6