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:
This is the continuation of Ref 4.
Reference :
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.