Lecture 3.

 

Second order differential equations – Part 2

 

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

 

 

 

 

This is the continuation of Ref 4.

 

 

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

 

Examples and problem are taken from:

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

 

Section 10. Diffusion between Compartments

 

Fig 3.9   . Diffusion between compartments. The volumes V1 and V2 are not constant but if the initial concentrations C1(0) and C2(0) are very equilibrium can be reached without an appreciable volume change.

 

 

In an interval ∆t a quantity ∆Q (moles) diffuses from left to right proportional to the product of a diffusion constant  (K~volume/time) and the difference in concentrations C2 – C1,

 

∆Q = ∆t K ( C2 – C1  ) ~time (volume/time) ( moles/volume) .

Note that if C2 – C1  < 0 then compartment #1 ,on the left, empties into #2.

 

Dividing by V1 and ∆t gives for C1 = Q(t)/V1 ,

 

d C1/dt =  (K/V1) {C2(t) – C1 (t) }               (20)

 

To obtain the equation for C2 let V= V2 and write instead {C1(t) – C2 (t)}

because if one concentration is increasing the other is decreasing.

 

d C2/dt =  (K/V2) {C1(t) – C2 (t) }     .       (21) 

 

The initial conditions are C1(0) and  C2(0)  ~moles/volume .     

 

Solving (20 ) - (21) by finite differences gives

 

C1( tn ) =  C 1(tn-1 ) +  ∆t (K/V1) {C2(tn-1) – C1 (tn-1) }     (22)          

C2( tn ) =  C 2(tn-1 ) +  ∆t (K/V2) {C1(tn-1) – C2 (tn-1) }   .  (23)

 

There are two time scales

 

tscale 1 = V1/K    and   tscale 2 = V2/K   .

To assume that the volumes are constant would require that

tequilibrium << tscale1  so that  V1 is approximately constant  , or

tequilibrium << tscale2 so that  V2 is approximately constant. 

 

 

Numerical example : Let  V1 = V2 = 10 liters  , K=(1/3) liter/minute

C1(0)=0.3 moles /liter    , C2(0) =0. 

 

Fig. 3.10 Difussion between compartments .Asuming  that V1 and V2 are constants which can not be entirely right if there is a flow from V1 to V2 .

 

Fig 3.11   C1(0)=.3 moles/liter , C2(0)=.2moles/liter  , V1=20 liters , V2=10 liters , K = 2 liter/min

 

Fortran code

c Diffusion between two compartments

      real K

      data v1,v2/10.,10./

      k=1./3.

      c10=0.3

      c20=0.

      ts1=v1/k

      ts2=v2/k

      tlarge=amax1(ts1,ts2)

      tsmall=amin1(ts1,ts2)

      dt=tsmall/200.

      tfinal=4.*tlarge

      nstep=int(tfinal/dt)

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

      kount=kp

      print 100 , 0.,c10,c20

      do 10 i=1,nstep

      t=dt*float(i)

      c11=c10+dt*(k/v1)*(c20-c10)

      c21=c20+dt*(k/v2)*(c10-c20)

      if(i.eq.kount)then

      print 100 , t,c11,c21

      kount=kount+kp

      endif

      c10=c11

      c20=c21

10    continue

100   format(1x,'t,c1,c2~moles/liter=',3(3x,e10.3))

      stop

      end

     

 

 

Section 11. A simple kinetic reaction

 

A reaction has the chemical formula

 

a + b → p                                                                    .  (24)

 

Let the concentration be represented by capital letters A ~ moles/liter ,

B~ moles/liter , P ~ moles/liter.

The concentration P grows at the rate

 

dP/dt = K A B  ~ moles/(liter-time)                           (25)

The contant K which is specific to the particular reaction and thermodynamics conditions, say temperature and pressure, has dimensions

 

K~ volume/(mole-time) . Eq (25) can be written for only one dependent variable , P(t)  , the product concentration.

 

Letting  A(t) = A0 – P(t)     and  B(t) = B0 – P(t)   where A0 ,B0 are initial

 

concentrations. We transform (25) into

 

dP/dt = K ( A0 – P(t) ) ( B0 – P(t))                            ,  (26)

with initial condition P(0)=0.

       

The finite difference solution  of (26) is ,

 

 

P( tn ) =  P (t n-1) + ∆t * K ( A0 – P(tn-1) ) ( B0 – P(tn-1))  

A Fortran code (see below) was employed to reduce the following plot.

A(0)=0.5moles/liter , B(0)=0.2moles/liter ,P(0)=0 , K=1liter/(mole-minute).                        

Fig 3.12

 

FORTRAN code

c Diffusion between two compartments

c a ~moles/liter  time ~ minutes  ,K ~ liters/(mole-minute)

      real K

      data aini , bini , k/.5 ,.2 ,1.   /

      p0=0.

      ts1=1./(aini*k)

      ts2=1./(bini*k)

      tlarge=amax1(ts1,ts2)

      tsmall=amin1(ts1,ts2)

      dt=tsmall/200.

      tfinal=2.*tlarge

      nstep=int(tfinal/dt)

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

      kount=kp

      print 100 , 0.,aini,bini ,p0

      do 10 i=1,nstep

      t=dt*float(i)

      p1=p0+dt*k*(aini-p0)*(bini-p0)

      if(i.eq.kount)then

      print 100 , t,aini-p1,bini-p1, p1

      kount=kount+kp

      endif

      p0=p1

10    continue

100   format('t~min,a,b,p~moles/liter=',4(3x,e10.3))

      stop

      end

 

 

 

Problems:

 

Find by using finite difference the solution to

 

d2 C/dt2  + a t dC/dt +exp(-bt) C(t)2 = 0 .         (27)

This is a non linear DE with variable coefficients. Table III-1 is not useful in this problem.

Let

a =3./sec2   b= 1/sec  , C(0)=1 , dC(0)/dt=0.

 

The generic solution is

 

C(tn ) = 2C(tn-1) –C(tn-2)  + (∆t)2 { -a tn-1 (C(tn-1) –C(tn-2) ) /∆t

      

         -exp(-btn-1) C 2(tn-1 ) }                               .  (28)

 

Thee are two time scales

 

t1 = (1/a)1/2 = ( 1/3)1/2  and    t2 =1/b  =1    .

 

Select    ∆t << t1 .

 

Fig 3.13 Solution of  eq . (27) .

 

 

MATLAB code

% Chap 3 problems-non linear DE

C(1) =1 ; C(2) =1;

a= 3 ; b=1; ts1=(1/a)^.5 ; ts2=1/b  ;

 

 tf=5.*ts1 ; nstep=5000 ; dt=tf/nstep;

 

k=1;

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

    k=k+1;

     acel=( -a*(t-dt)*( C(k)-C(k-1) )/dt...

         -exp(-b*(t-dt))*C(k)^2 )  ;                           

 

     C(k+1) = 2*C(k)-C(k-1) +(dt^2)*acel;

 end

t=[0:dt:tf];

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

Chapter 3 . Problem 17.

 

 A beaker of water initially at 20 celsius is placed over a hot plate at 80 celsius. Ten minutes later the beaker water reaches 50  celsius. Develop an eqaution for the temperature as  a function of time. Assume the heat loss rate  from the beaker is proportional to the difference in temperature between the beaker and the room.

 

Data: Initial temperature 20 C , T room =20 C  ,final temperature 80 C ,

heat loss rate ~  constant*( Tbeaker (t) – Troom)  ~  constant*( Tbeaker (t) – 20) .

and T(t=10 minutes) = 50 celsius.

 

Guessing the solution.  Suppose the solution includes  an inverted exponential and a constant term.

 

T(t) = A + B ( 1- exp(-λt)        .

We have three data points that allow finding A , B and λ.

T(0)= 20 = A + 0      ;   A=20 celsius

T(t=∞) = 80 = 20 + B   ; B =60 celsius

 

T(t=10min) = 50 = 20 +60 (1-exp(-10λ) )

 

exp(-10λ) = 1/2

λ = ln(2) /10 = 6.93E-2 min-1 .

 

Hence the equation is

T(t) = 20 + 60 {1-exp(-6.93E-2 t) } , t~ minutes

 

From the plot (see below) we see it has reached 65 celsius in 20 minutes.

 

Fig 3.14 Temperature of the beaker.

 

MATLAB code

lambda=6.93E-2;

t=[0:1:25];

Temp= 20 + 60*(1-exp(-lambda*t) );

plot(t,Temp),xlabel('t ~ minutes'),ylabel('Temp')

 

 

 

 Continuation of problem 17. Obtain the DE.

 

Suppose the heat gained by the beaker is proportional to

the difference  (T hot plate – T room )

 

∆ Q gained ~ ( T hot plate – T room ) .

The beaker loses heat  ∆ Q lost ~  -( Tbeaker(t) – Troom ).

The net heat flow is

 

∆ Q ~  T hot plate – T room  - (Tbeaker(t) – Troom)

 

       ~  (T hot plate - Tbeaker(t) )      .

 

The temperature increase in an interval ∆t is  proportional to ∆ Q ,

 

∆ T = ∆t λ λ ~ 1/time . The DE for the beaker temperature  is

 

dT/dt = λ ( 80-T ).

 Integrating gives

- ln( 80 –T) =  λ t + constant . Solving for T(t)

 

80- T(t) = A exp( - λ t)    or

 

T(t)  = 80 + A exp(- λ t) .

 

Applying the IC  T(t=0) =20 determines  the value A= -60 celsius.Then

T(t) = 80 -60 exp(- λ t) or   T(t)= 20 + 60(1-exp(- λ t) .

 

Problem 18.

A patient is given 5 microcuries of   iodine , 131 I. Two hour later 0.5 microcuries had been absorbed by the thyroid. Assume it is a linear process an estimate the absorption at two hours if the initial dose was 15 microcuries. Since the initial dose has been multiplied by three the absorption at two hours increases by the same factor of three. Thus we can expect the thyroid to absorb 3 (0.5) =1.5 microcuries of 131 I .

 

END OF CHAPTER 3.