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

y(x₀ + ∆x ) = y(x₀) + ∆x (dy/dx)₀ + (1/2) (∆x)² (d² y/dx²)

+ ..... (1/n!) (∆x)ⁿ (dⁿ y/dxⁿ) . (6.1)

In practice we can do with a few terms.

Example 1. (equation 6.4 from ref 1) Consider

dC₂/dt = (F/V₂)C₁(0) exp(-(F/V₁)t) -(F/V₂) C₂(t) (6.2)

with the parameter values V₁= 10 L , V₂ =20 L , F= 1 L/sec and IC

C₁(0)=.3 mole/L , C₂(0)=0.

Substitution in (6.2) gives

dC₂/dt = (1/20)(.3) exp(-.1 t) –(1/20) C₂(t)

dC₂/dt = .015 exp(-.1t) -.05 C₂(t) (6.3)

The second derivative is

d² C₂/dt ² = -1.5E-3 *exp(-.1 t) -.05 dC₂/dt

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

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

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

C₂(t₀ +∆t) = C₂(t₀) + ∆t{ .015 exp(-.1t₀) -.05 C₂(t₀) }

+(1/2)( ∆t)²{ -1.5E-3 *exp(-.1 t₀) -.05 [ .015 exp(-.1t₀) -.05 C₂(t₀) ] }

(6.5)

In this equation there are two time scales defined by

tscale1= V₁/F = 10 sec ; tscale2 = V₂/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 2^nd order DE

d²y/dt² + b dy/dt + c y =0 , (6.6)

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

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

The solution to second order is,

y(t₀ + ∆t ) = y(t₀) + ∆t (dy/dt)₀ + (1/2) (∆t)² (d² y/dt²) . (6.7)

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

dy( t₀ +∆t )/dt = dy( t₀)/dt + ∆t [d² y(t₀) /dt²] (6.8)

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

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

y(t₀ + ∆t ) = y(t₀) + ∆t (dy/dt)₀ + (1/2) (∆t)² (- (dy/dt)₀ - y₀ ) (6.9)

and

dy( t₀ +∆t )/dt = (dy/dt)₀ + ∆t ( - (dy/dt)₀ - y₀ ) . (6.10)

Fig. 6.2 Solution of d²y/dt² + 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 = - K₁ A + K₃ C (6.11)

dB/dt = - K₂ B + K₁ A

dC/dt = - K₃ A + K₂ 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 K_i.

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

say K₁ (t) = K₁₀+ m t .

The finite difference solutions to (6.11) are

A(t_n ) = A( t _n-1 ) + ∆t { - K₁(t_n-1) A(t_n-1) + K₃ C (t_n-1) } (6.12)

B(t_n ) = B( t _n-1 ) + ∆t { - K₂ B(t_n-1) + K₁(_tn-1) A(t_n-1) }

C(t_n ) = C( t _n-1 ) + ∆t { - K₃ C(t_n-1) + K₂(t _n-1) B(t_n-1) }

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

Let the rate constants be ;

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

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

while K₂= 2/s , K₃ = 1.0/s .

The slope m of K₁ is (9.9/2) /s² = 4.95/s²

The results are plotted in fig 6.4

Fig 6.4 Concentrations for the circular reaction with variable K₁.

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(t_n ) = A( t _n-1 ) + ∆t { R1- aB(t_n-1) – R1’ } (6.15)

and

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

B(t_n ) = B( t _n-1 ) + ∆t { R2 – R2’ + bA(t_n-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 /s² , C(0)=2.

tscale = 1/ a^1/2 = 3.16 s

The finite difference solution is

C(t_n) = C( t_n-1 ) + ∆t { a t_n-1 C(t_n-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(t_n) = y(t_n-1) +∆t { 100 + a y(t_n-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)² .

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

d² y/dt² –λ (1- y²)(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

d² y/dt² +b (dy/dt) + c y = 0 .

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

b ~ –λ (1- y²).

In the damped oscillator the solution is governed by terms

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

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

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

exponentials and therefore no oscillatory behavior.

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

Another possibility is that at an instant (1-y²) < 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-y²) > 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