Lecture 2.
First order differential equations in dilution problems and other situations
by Reinaldo Baretti Machνn
http://www1.uprh.edu/rbaretti/methodsoftheoreticalphysics.htm
http://www1.uprh.edu/rbaretti/methodsoftheoreticalphysicsPart2.htm
http://www1.uprh.edu/rbaretti/methodsoftheoreticalphysicsPart3.htm
Para preguntas y sugerencias escriba a :
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
Examples and problem from:
Chapter 2 . Reference 1. The problem numbers refer to the textbook listing.
Section 1. One compartment dilution problem.
Fig 1.
A bucket contains say V=10 liters of water. The volume is kept constant becaues the inflow F =2 liter/minute is adjusted is equal to the outflow F at the bottom. Thres moles of a tracer Q(0)= 3 moles are injetct at t=0 and we assume it is instantly mixed throughout the volume V.
At any instant ,the concentration C(t) (moles/liter) = Q(t) /V .
The initial concentration C(0) = 3 moles/10 liters =0.3 moles/liter
Establish a first order differential equation and solved it.
In an interval ∆t (minutes) the number of moles diminishes by ∆Q .
A finite difference equation is
∆Q = - ∆t F C(t) (1)
moles ~ minutes * (liters/minute) * (moles/liter)
We can write an equation for either Q or C.
Substitute ∆Q = V ∆C above an get
∆C = - ∆t ( F/V) C(t) . The differential equation follows from invoking the
limit (∆t→0) (∆C / ∆t ) = - (F/V) C
i.e dC/dt = -(F/V) C (2)
This is of the form dy/dt = -k y (see chapter 1) and the solution is
y(t) = y(0)exp(-kt) . Hence the solution is,
C(t) = C(0) exp( -(F/V) t) . The numerical values are
C(0)= 0.3 moles/liter , F/V =k or λ = 2(liters/min)/10(liters) = .2 /minute
Tscale = 1/λ = 5 minutes The plot of C(t) is shown next.
Plots of C(t) analytical and the numerical result from an integration.
The differential equation
dC/dt = -(F/V) C leads to the finite difference equation
(Cn Cn-1 )/∆t = -(F/V) Cn-1 wrere Cn = C( n ∆t)
Solving for Cn
Cn =Cn-1 - ∆t (F/V) Cn-1 (3)
C0 = C(t=0) is given , then the solution is constructed from the previous
values of C which are known.
C1 =C0 - ∆t (F/V) C0 ,
C2 =C1 - ∆t (F/V) C1
Now ∆t has to be small relative to the time scale 1/λ .
The following MATLAB code integrates numerically the DE.
The value of ∆t is chosen to be small in comparison with Tscale.
C analytical i.e C(0) exp(-λt) and C numerical are plotted together. They are practically identical within the resolution of the graph.
tscale = 5 minutes , dt= 2.0 E-2 minutes
MATLAB code
%numerical solution of dC/dt= -(F/V)*C(t)
F=2 ; V=10 ; Q0= 3; c(1)=Q0/V ;
lambda=F/V; tscale=1/lambda ,
cexact=c(1)*exp(-lambda*t);
ti=0 ; tf = 4*tscale ; nstep=1000;
dt=(tf-ti)/nstep ,
% k is the index kounter it starts with k=1 at t=0
% c(1) is given as initial condition , we wil find c*2) ,
% c(3)...etc
k=0 ;
for x=[ti+dt:dt:tf]
k=k+1;
c(k+1)=c(k) +dt*(-(F/V)*c(k));
end
t=[ti:dt:tf];
plot(t,c,t,cexact), xlabel('t(minutes)'),ylabel('C,cexact')
Other examples of exponential decay.
Section 2.
Consider a metabolite of Q0 moles to which a small quantity of radioactively
labeled q(0) is added at t=0. ( i.e q(0) << Q0 ) . The total quantity is
Qtotal = Q0 + q(t) ≈Q0 .
q decays at a fixed rate of R ~ moles/time
Before writing an equation for ∆q lets rewrite eq (1)
∆Q = - ∆t F C(t) as
change in q =(moles)= - ∆t (time) R (moles/time) (radioactive moles/total moles)
∆q = - ∆t R q(t)/Q0 from which we can pass to the DE
dq/dt = -k q k= R/Q0 ~ 1/time (4)
Radioactive decay
It is empirically found that from a total of N atoms, there is a decay , ∆N<0 , in an interval ∆t such that
∆N =- λ N ∆t , λ~ 1/time . The corresponding DE is
dN/dt = - λ N(t) . Let the initial condition be N(t=0) = N0 , then
N = N0 exp( - λt)
Absorption of X- rays
The intensity I ~ energy/(time-area) , of X-rays after penetrating a material slab of thickness ∆x is diminished acording to
∆I = - D I(x) ∆x , D ~ 1/length
this leads to the DE dI/dx = - D I . Let the initital condition
be I(x=0) = I0 , then I (x) = Io exp(-Dx)
Sec 5. One compartment with constant injection (see figure)
Assume the intial concentration of substance Q is zero , C(0)=0.
A volume V~ liter is kept fixed by a constant inflow flow and outflow of F ~
liters /time . The quantity of Q(t) moles of a tracer suffers in an interval
a loss
∆Qloss = - ∆t F Q/V ~ - ∆t ( volume/time) * (moles/volume).
At the same time a quantity ∆t R ~ time *(moles/time) is injected instantly.
R is the rate of injection in moles/time .
∆Qgain = ∆t R . Adding the loss and gain parts of Q ,
∆Q = ∆t ( R - F Q/V) = ∆t ( R - F C(t) ).
Dividing by ∆t V and taking the limit ∆t→0 , results in the DE
d C /dt = ( R/V (F/V) C(t) ) , (5)
with initial condition C(t=0) =0.
Remember R , V , F are parameters. They assume constant values in a given example. The independent variable is time , the dependent variable is the concentration C(t) ~ moles/volume.
The DE equation (5) is of a type called separable. One can separate the variables C and t , in two different integrals.
∫ dC /(b -a C(t) ) = ∫ dt , where , b=R/V a= - F/V .
If you go to a table of integrals you will find the indefinite integral in a form like this ;
∫ dx /(b+ax) = (1/a) ln (b+ax) . Hence
∫ dC / { (R/V) (F/V) C(t) } =
-(1/(F/V) ) ln [ (R/V) (F/V) C(t)] .
Equating this to the right hand side integral gives
the indefinite integral
-(1/(F/V) ) ln [ (R/V) (F/V) C(t)] = t + A
A is the integration constant that can be related to the initial condition as we shall show.
It is convenient to rewrite the integral as
ln [ (R/V) (F/V) C(t)] = - (F/V) t + A* where A* = (F/V) A is still a constant. To simplify matters call A all constants that originate with the original A.
Then
(R/V) (F/V) C(t) = A exp(-(F/V)t ) , (another A)
and solving for C(t)
C(t) = (V/F) { (R/V) - A exp(-(F/V)t ) } (6).
Apply now the initial condition to determine A
C(0)= 0 = (V/F) { (R/V) A } , thus A =R/V and finally
C(t) = (R/F) { 1- exp( -(F/V) t) } (7)
which dimensionally checks ,
moles/volume = ( moles/time)( time/volume) .
At t= infinty C(∞) = R/F = constant .
The time scale of the problem is T scale = V/F
Example : Plot of eq (7) . Such plot is called an inverted exponential.
set % R=.1 mol sec , F =2 liter/(60 sec) , V=10 liters
Matlab code
% One compartment with constant injection
% R=.1 mol sec , F =2 liter/(60 sec) , V=10 liters
R=0.1 , F =2/60 , V=10 ,
tscale=V/F; tfinal=3*tscale ;dt = tfinal/100;
t=[0:dt:tfinal];
c=(R/F)*(1-exp(-t/tscale));
plot(t,c) , xlabel('t(sec)') ,...
ylabel('C(t)~moles/liter ' )
The two compartment series dilution
The one compartment dilution problem have already been solved
and the solution is C1(t)=C1(0) exp( -(F/V1) t) where the subscript 1 refers to compartment #1. The initial condition was C1(0)= Q0 /V1 i.e a total amount of tracer Q0 is injected and mixed instantly at t=0.
Volumes V1 and V2 are kept constant by insuring that F = F1 = F2 .
The problem is to set up a differential equation for C2(t) and to find its solution. Remember C2(t) = Q2(t) /V2 and Q2 ~ moles in the second compartment. The initial condition is C2(t) =0.
The quantity ∆Q2 ~ moles ,in compartment 2 within an interval ∆t is increased by the amount
∆Q2 gain = inflow rate * C1 * ∆t~ (liter/sec)(moles/liter) (sec)
= F C1 ∆t
The loss is ( negative quantity)
∆Q2 loss = - outflow rate * C2 * ∆t~ (liter/sec)(moles/liter) (sec).
= - F C2 (t) ∆t
So the net change is
∆Q2 = F C1 ∆t - F C2 (t) ∆t and the corresponding DE is
d Q2 /dt = F C1 - F C2 (t) divide by V2 to obtain the DE for C2 ,
dC2/dt = (F/V2 ) ( C1(t) C2 (t) ) . (8)
Rearranging (8) , and introducing the solution for C1 get
dC2/dt +( F/V2 ) C2 (t) = (F/V2 ) C1(0) exp( -(F/V1) t ) . (9)
This is a first order linear DE of the form
dy/dx + P(x) y = Q(x) (10)
which has the general solution
y = exp(-I) ∫ Q(x) exp(I) dx + B exp(-I) (11)
B is an integration constant and I = ∫ P(x) dx .
To avoid the burden of too many symbols lets put some data in (9).
Set , C1 (0)=0.3moles/liter , V1= 10 liters , V2= 5 liters ,
F= 2liters/min= 2liters/(60sec) = 3.33E-2 liters/sec.
There are two time scales ; (F/V1) = 3.33E-3 sec-1 ; Tscale1 = V1/F= 300 sec and (F/V2) = 6.66E-3 sec-1 ; Tscale2 = V2/F = 150 sec.
With the proposed data and the value of C1(0)=.3 moles/liter , eq (9) becomes
dC2/dt + 6.67E-3 C2 (t) = 6.67E-3 *(.3) * exp( -3.33E-3 t ) . (12)
Identify P(t) = 6.67E-3 ~ 1/sec and
Q(t) = 2.00E-3 *exp(-3.33E-3 t) ~ moles/(liter-sec)
Proceed to implement the method of eq (11)
I = ∫ P(t) dt = 6.67E-3 * t (13)
∫Q(x) exp(I) dx = ∫ 2.00E-3 *exp(-3.33E-3 t)*exp(6.67E-3 t) dt
= ∫ 2.00E-3 *exp(+3.33E-3 t) dt
= ( 2.00E-3/3.33E-3) exp( +3.33E-3 t)
= 6.01E-1 exp( 3.33E-3 t) . (14)
Inserting (13), (14) in (11) results in
C2(t) = exp(-6.67E-3 t) * 6.01E-1 exp( 3.33E-3 t) + Bexp(-6.67E-3 t)
= 6.01E-1 exp(-3.33E-3t) + B exp(-6.67E-3 t) .
Now invoke the initail condition
C2(0)= 0 = 6.01E-1 + B ; thus B= -6.01E-1 and
C2(t) = 6.01E-1 { exp(-3.33E-3t) exp (-6.67E-3 t) } . (15)
Applying (11) but keeping the data (the parameters) unspecified gives,
C2(t) = [ V1 C1(0)/(V1-V2 )] {exp(-(F/V1)t ) exp(-(F/V2)t ) } . (16)
Plot of C1 , C2
Matlab code
% two compartment dilution
tscale=1/6.67E-3 ; tf=5*tscale;
dt=tf/200;
t=[0:dt:tf];
c1=.3*exp(-3.33E-3*t);
c2=6.01E-1*(exp(-3.33E-3*t)-exp(-6.67E-3*t) );
plot(t,c2,t,c1), xlabel('t~sec') ,...
ylabel('C1 ,C2~moles/liter')
Roughly speaking without solving the DE we can anticipate the following
results. The 3 moles are emptied from compartment one in roughly 300 seconds ( tscale 1). If there is no loss in compartment #2 the concentration would be 3 moles/5 liters = 0.6 moles/liter but in 300 sec there are two half lifes of tscale2 (tscale2=150 sec) . Thus the concentratio reduces to (1/2)2=1/4 . So at t=300 sec the concentration at compartment # 2 ,is at a maximum and roughly equals (1/4)(.6moles/liter) = 0.15 moles/liter .
A check of the plot above bears out this estimate.
Numerical integration of eq.(8).
We expand next on the procedure already used ( in these notes) to integrate numerically a DE by finite differences.
The equation
dC2/dt +( F/V2 ) C2 (t) = (F/V2 ) C1(0) exp( -(F/V1) t ) (17)
has initial condition C2(t=0) = 0 . Suppose you want the solution in intervals
of size ∆t where , ∆t = ( tfinal tinitial) /nsteps , and ∆t << tscale1 or
∆t << tscale2 , whichever is the samller one.
The finite difference generic solution is ,
C2 (tn) = C2 (tn-1) +
∆t {- ( F/V2 ) C2 (tn-1) + (F/V2 ) C1(0) exp( -(F/V1) tn-1 ) } . (18)
note that every term in the right hand is known.
The Fortran given below uses this method. The plot is given next , which , of course is identical to the previous one.
Same initial data C1(0)= 0.3 moles/liter , c2(0)=0, but V1=V2=10 liters.
Fortran code (two compartment series dilution)
c the two compartment series dilution
c F~liters/min ,V ~liters
data c10 /.3 /
data F, V1,V2/2. ,10., 5./
c1(t) = c10*exp(-F*t/v1)
c f is changed to liters /second by the next declaration
f=f*(1./60.)
tscale1=V1/f
tscale2=V2/f
tsmall=amin1(tscale1,tscale2)
tlarge=amax1(tscale1,tscale2)
dt= tsmall/50.
nstep=int(3.*tlarge/dt)
c20=0.
kp=int(float(nstep)/70.)
kount=kp
print 100 , 0. , c1(0.) ,c20
do 10 i=1,nstep
t=dt*float(i)
c21= c20 + dt*(-(F/v2)*c20 +(f/v2)*c1(t-dt) )
if(i.eq.kount)then
print 100 ,t ,c1(t),c21
kount=kount+kp
endif
c20=c21
10 continue
100 format(1x,'t(sec),c1, c2 ~moles/liter=',3(3x,e10.3))
stop
end
Section 8. Mother Daughter Radioactive Decay
This problem is analogous to the two compartment dilution series.
A mother (M) nuclei species decays like
dM/dt = - λ1 M ( 19 ).
The daughter nuclei ( D) are born at the rate + λ1 M and decay at a rate
- λ2 D , therefore the DE for the daughter nucei is
d D/dt = λ1 M(t) - λ2 D (t) (20)
The equation is analogous ( though not the same) to
dC2/dt = (F/V2 ) ( C1(t) C2 (t) ) if we make C1=M , C2= D
and λ1 = λ2 = (F/V2) . But in general λ1 ≠ λ2 .
We need two initial conditions to solve (19) and (20) , set M(0)= 1mol
D(0)=0 and tscale1 =1/ λ1 = 5 minutes , tscale2 = 1/λ2 =10 minutes.
We can modify slightly the above FORTRAN code to integrate , both equations simultaneously. The finite difference solutions to (19)-(20) are
M(tn ) = M(tn-1) - ∆t λ1 M(t n-1) (21)
D( tn ) = D(tn-1) + ∆t ( λ1 M(tn-1 ) - λ2 D (tn-1) ) . (22)
Solution of (21)-(22) with the above parameters
FORTRAN code
c Fortran code for mother daughter radioactive decay
real m0 ,m1 ,lambda1, lambda2
data m0 ,D0 /1.,0. /
c time scales are in minutes
data tscale1 , tscale2/5.,10. /
lambda1=1./tscale1
lambda2=1./tscale2
tsmall=amin1(tscale1,tscale2)
tlarge=amax1(tscale1,tscale2)
dt= tsmall/50.
nstep=int(3.*tlarge/dt)
kp=int(float(nstep)/70.)
kount=kp
print 100 , 0. , m0 ,d0
do 10 i=1,nstep
t=dt*float(i)
M1 = M0 - lambda1*dt*M0
D1=d0+ dt*(lambda1*m0-lambda2*d0)
if(i.eq.kount)then
print 100 ,t ,m1,d1
kount=kount+kp
endif
m0=m1
d0=d1
10 continue
100 format(1x,'t(min),mother,daughter~moles=',3(3x,e10.3))
stop
end
Section 9. Circular reactions
The concentration A has a gain in an interval ∆t of
∆t k3 C C (third line) and a loss of - ∆t k1 A A (first line). The balance is
∆A = ∆t k3 C C - ∆t k1 A A , which leads to the DE
dA/dt = k3 C C - k1 A A (23)
the ki ~ ( 1/time)*(1/concentration) .
Similarly
dB/dt = k1 A A - k2 B B (24)
dC/dt = k2 B B - k3 C C . (25)
To simplify matters suppose A , B and C are practically constants and can be absorbed into new constants , capital K1 , K2 , K3 .
The capital Ks are ; K1 ≡ k1 A , K2 ≡ k2 B , K3 ≡ k3 C ~1/time . The
equations take the new form
dA/dt = K3 C - K1 A , 26-a
dB/dt= K1A -K2 B , 26-b
dC/dt = K2 B K3 C . 26-c
We will solve numerically this system of three first order coupled DE.
Let tscale1= 10 minutes =1/K1 , tscale2=5minutes =1/K2 ,
Tscale3= 3minutes=1/K3 . A(0)=1.0 moles , B(0)=0 moles , C(0)=0 moles
The finite difference solutions to eqs 26 a to 26 c are
A( tn ) =A(tn-1) + ∆t { K3 C (tn-1) - K1 A(tn-1 ) }
B( tn ) =B(tn-1) + ∆t { K1 A (tn-1) - K2 B(tn-1 ) }
C( tn ) =C(tn-1) + ∆t { K2 B (tn-1) - K3 C(tn-1 ) }
and ∆t has to be small relative to the smallest of the Tscale. In the present
case ∆t << Tscale 3.
The next FORTRAN code expands the mother daughter code to deal with the circular reaction. The results are plotted next.
We see that in about 10 minutes a steady state is reached with constant concentrations
A=.56 , B= .28 , C= .16 .
Problem 3. A substance Q~ moles , diffuses from the compartment of volume V at a rate proportional to the product of the diffusion constant K ~ volume/time and the difference in concentration with the environment (C- Cenvironment) ~ moles/volume. The initial condition is C(0)= Q0/V.
Also assume C > Cenvironment hence there is a loss of Q in the compartment,
V is kept constant , the mechanism is not shown in the drawing.
Then in an interval ∆t , a loss occurs
∆Q = - ∆t K (C- Cenvironment) ~ time( volume/time)*(moles/volume)
Dividing by the volume V and by ∆t we obtain the differential equation
dC/dt = -(K/V) (C(t) - Cenvironment) .
Problem 4. Assume that in the compartement of Prob. 3 a quantity ∆Qin
is pumped given by flow rate R(moles/time) * ∆t .
The net ∆Q is now
∆Q = - ∆t K (C- Cenvironment) + R ∆t and the DE is
dC/dt = - (K/V)(C- CE ) + R/V
Example : Solve the DE dC/dt = - (K/V)(C- CE ) + R/V .
Let K= F=2liter/minute , V=10 liter ,C(0)=0.3 moles/liter ,
Cenvironmnet=.05 moles/liter , R =.1mole/minute.
Tscale = 1/(F/V) =10/2 = 5 minutes .
The finite difference solution is
C ( tn) = C(tn-1) + ∆t { - (K/V)(C(tn-1) - CE ) + R/V }.
FORTRAN code
real k
data c0 , cE/.3, .10/
data k,v, R/2.,10.,0.1/
c time scales are in minutes
tscale=V/K
dt= tscale/100.
nstep=int(6.5*tscale/dt)
kp=int(float(nstep)/70.)
kount=kp
print 100 , 0. , c0
do 10 i=1,nstep
t=dt*float(i)
c1 = c0 +dt*( -(K/v)*(c0-ce) +R/V)
if(i.eq.kount)then
print 100 ,t ,c1
kount=kount+kp
endif
c0=c1
10 continue
100 format(1x,'t(min),C~ moles/liter=',2(3x,e10.3))
stop
end
Problem 5 . In a one compartment dilution the concentration falls with time according to the data :
|
t(sec) |
C (moles/liter) |
|
0 |
.024 |
|
1 |
.011 |
|
2 |
4.8E-3 |
|
3 |
2.4E-3 |
|
4 |
1.0E-3 |
If Q(0)=0.1 mole , find F ~liters/sec and V ~liters
C(0) =.024 ≡ Q(0)/V , hence V= .1/.024 = 4.17 liters
If we assume a perfect fit , which may not be the case ,
then C(t) = C(0) exp(- (F/4.17) t) . Any value from the table will serve to solve for F. At t=1 sec
.011 = .024 exp( -.240F) . Taking the natural log
ln(.011/.024) = -.780 = -.240F
F = 3.25 liters/second
Problem 6.
In a radioactive process half the atoms decayed in one day.What fraction existed at 10 hours. We know that the half life (t1/2 ), is related to the rate λ by
λ = .693 / t1/2 = .693 /24 hours = 2.89E-2 hour-1 .
After 10 hours the fraction of atoms left is
N/N0 = exp( - λ t) = exp( -2.89E-2 (10) ) = .749 .
Problem 8. Substance A is created at a rate R ~moles/sec and destroyed at the rate - K(1/sec) * A(moles) ~ moles/sec. Let the initial condition be A(0) =0.
In an interval ∆t the net change in A is
∆A = R ∆t K A ∆t
The corresponding DE is
dA/dt = R K A .
A steady state is reached when dA/dt =0 . This gives A(t=∞) = R/K.
Again the solution is an inverted exponential as mentioned in a previous example of the one compartment dilution with constant injection.
Data: R= .5 moles/sec , K = 0.1/sec , A(0)=0.
The asymptotic value of A is A (t=∞) = R/K=5 moles
MATLAB code
% R~moles/sec K ~(moles/mole-sec)
A(1)=0; R =.5 ; K= .1 ;
tscale= 1/K ;tfinal=3*tscale ;nstep=2000;
dt=tfinal/nstep;
n=0;
t=[dt:dt:tfinal];
n=n+1;
A(n+1)=A(n) + dt*(R-K*A(n));
t=[0:dt:tfinal];
plot(t,A),xlabel('t~sec'),ylabel('A~moles')
END OF CHAPTER II