Lecture 3.
Second order differential equations
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:
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
Examples and problem from:
Chapter 3 . Reference 1. The problem numbers refer to the textbook listing.
We deal here with some second order differential equations with constant coefficients as they appear in Table III-1 of reference 1.
The analytical approach to the solution is given in Ref. 1. Our purpose here is to present solutions in terms of the finite difference method, already employed in Lecture 2.
DE 6, of Table III-1
d2 C/dt2 = k C (1)
initial conditions (IC) are C(0) , dC(0)/dt .
The analytical approach consists in postulating a solution C ~exp(λt).
Substitution in (1) gives λ2 C = k C , thus λ = + k/12 or λ = - k/12 .
The general solution is
C(t) = A 1 exp(+ k/12 t) + A 2 exp(- k/12 t) . (2)
Some algebraic manipulations are necessary to fit the solutions to the data.
The constants A 1 , A 2 are related to the initial conditions. Notice that if
k<0 the solutions become oscillatory because
exp(+ k/12 t) = exp( i /k/1/2 t ) = cos ( /k/ 1/2 t) + i sin(/k/ 1/2 t) ,
and
exp(-k/12 t) = exp( -i /k/1/2 t ) = cos ( /k/ 1/2 t) - i sin(/k/ 1/2 t).
Applying the IC , we have for example for k>0
C(0) = A1 + A2 (3)
( 1/ k1/2 )dC(0)/dt = A1 - A2
Adding and subtracting yields
A1 = (1/2) { C(0) + 1/ k1/2 )dC(0)/dt } (4)
A2 =(1/2) { C(0) - 1/ k1/2 )dC(0)/dt }
Now we will solve (1) numerically with k=9 and IC C(0)=1 , dC(0)/dt=0.
The second derivative hast the finite difference approximation
d2 C/dt2 ≈ { (C2 – C1 ) /∆t - ( C1 – C0 )/ ∆t } / ∆t
= (C2 – 2C1 + C0 )/ ∆t 2 (5)
C0 is given and the second value is C1 = C0 + ∆t * ( dC(0)/dt ) .
Or in a more general form
d2 C/dt2 ≈ (Cn+1 – 2Cn + Cn-1 )/ ∆t 2 . (6)
So eq (1) d2 C/dt2 = k C , becomes
(Cn+1 – 2Cn + Cn-1 )/ ∆t 2 = k Cn . (7)
The solution is
Cn+1 =2Cn - Cn-1 + ∆t 2 k Cn or,
Cn =2Cn-1 - Cn-2 + ∆t 2 k Cn-1 . (8)
The units of k are ~1/time 2 . It defines a time scale Tscale=1/k1/2 .
When integrating one has to select ∆t << Tscale .
The following plot show the numerical solution along with the
analytical solution , Cexact (t)= (1/2) { exp( k1/2 t) + exp( -k1/2t) }.
Fig 3.1
% d^2c/dt^2= Kc IC c(0)=1 , dc/dt=0;
c(1)=1 ; K= 9 ;dcdt=0 ;
tscale= 1/K^.5 ;tfinal=3*tscale ;nstep=2000;
dt=tfinal/nstep;
c(2)=c(1)+dt*dcdt;
n=1;
for t=[2*dt:dt:tfinal];
n=n+1;
c(n+1)= 2*c(n)-c(n-1) + dt^2*K*c(n);
end
t=[0:dt:tfinal];
cexact=(1/2)*( exp(K^.5*t)+exp(-K^.5*t));
plot(t,c,t,cexact),xlabel('t~sec'),...
ylabel('c,cexact ')
DE 8 of Table III-1 .
d 2C /dt 2 + K1 dC/dt = K2 (9)
Dimensions : K1 ~ 1/time , K2 ~ dimensions of C/ time2 .
The analytic solution given in table III-1 is
C(t) = (K2/K1) t + ( A1/K1) ( 1- exp(-K1t) ) + A2 . (9-a)
The full name for this DE is inhomogeneous differential equation of second order with constant coefficients. The K2 term makes it inhomogeneous.
Two IC are recquired C(0) and dC(0)/dt.
Eq (9) brings up a new feature to our study of DE.
The general solution contains two classes of solution to (9) . First those that satisfy the homogeneous DE
d 2C /dt 2 + K1 dC/dt =0 . Those are called complementary solutions
Ccomplementary = a Y1(t) + b Y2(t) or comparing with (9-a)
= ( A1/K1) ( 1- exp(-K1t) ) + A2 (10)
the constanst a, b are related to the initial conditions. Y1 and Y2 are linearly independent solutions .
There is a third solution,called the particular solution Cp . This one satisfies the inhomogeneous DE
d 2Cp /dt 2 + K1 dCp /dt = K2 .
A little guessing shows that
Cparticular (t) = (K2 /K1 ) t so that
Cgeneral ( t) = Cp + C complementary = Cp + a Y1(t) + b Y2(t) . (11)
There are special “operator methods” for finding the general solution.
Notice that Cp does not introduce arbitray constants. To use (11) in a numerical problem one still has to find a and b from the initial conditions.
That is one establishes two equations
Cgeneral ( t) = Cp (0) + a Y1(0) + b Y2(0) (12)
dCgeneral ( 0)/dt = dCp (0)/dt + a dY1(0)/dt + b dY2(0)/dt ,
and solves for a and b.
Now we will solve (9) in a concrete example without having to invoke any of the above considerations about different solutions and constants of integrations etc.
The finite difference expression for (9) is
(Cn -2 Cn-1 + Cn-2 )/ ∆t2 = - K1 (Cn-1 - Cn-2)/ ∆t + K2
We solve for Cn
Cn = 2Cn-1 - Cn-2 - ∆t K1 (Cn-1 - Cn-2) + ∆t 2 K2 . (13)
Whatever the dimensions of C may be , it is convenient to examine eq (9) by analogy with the dynamics of a particle in one dimension .
In that case , C~x(t) ~ position ~meters. The analogy between the derivatives is
d 2C/dt 2 ~ d2 x/dt2 ~ acceleration ~m/s2
dC/dt ~ dx/dt ~ velocity ~m/s
If K1 > 0 , the term K1 dC/dt is like a friction force/unit mass
If K2 < 0 there is a constant force toward the . i.e. -X axis.
There is a final (terminal speed) dC/dt = K2/K1 . This problem was touched in Prob. 8 , chapter 2 .
http://www.geocities.com/serienumerica4/Lecture2Biolologymedicine.doc
Let C(0)=0 meters , dC(0)/dt=0 m/s , K1 = .05/sec ,
K2= -9.8 meters/s2 . The asymptotic value of the velocity would be
K2/K1 = - 196 m/s
Fig . 3.2 . x(m) vs t(sec) .
In the next figure we provide an approximation to x for
0 <t < Tscale =20 sec= 1/K1 , without solving the DE.
Since vinitial = 0 , and vterminal = K2/K1 we take an average speed (1/2) K2/K1 and x≈ (1/2) (K2/K1) t . Fig 3.3 shows this equation.
Fig 3.3 x= (1/2) vterminal t = (1/2) (K2//K1) t . Compare with Fig 3.2 .
Figure 3.4 v(m/s) vs t(sec)
The acceleration in this problem goes from K2 to approximately zero in a time span of Tscale=20 sec.
So let the average acceleration be (1/2)K2 and approximate the velocity by
v = (1/2)K2 t .
Fig 3.5 Approximation v = (1/2) K2 t . Compare with Fig 3.4 .
MATLAB code
C(1)=0; K1 =.05 ; K2= -9.8 ;
tscale= 1/K1 ;tfinal=tscale ;nstep=3000;
dt=tfinal/nstep; v(1)=0 ;v(2)=0;
C(2) =C(1) +dt*v(1) ;
n=1;
for t=[2*dt:dt:tfinal];
n=n+1;
C(n+1)=2*C(n)-C(n-1) -dt*K1*(C(n)-C(n-1))...
+dt^2*K2 ;
v(n+1)=(C(n+1)-C(n))/dt;
end
t=[0:dt:tfinal];
plot(t,C), xlabel('tsec'),ylabel('x~m');
%plot(t,v), xlabel('t~~sec '),ylabel('v~m/s')
MATALAB code approximations to x and v
% orders of magnitude K1~1/sec K2~ m/sec^2;
K1=.05 ; K2=-9.8 ; tscale=1/K1 ;
vterm=K2/K1;
tfinal=tscale; nstep=100;
dt=tfinal/nstep;
t=[0:dt:tfinal];
x=(1/2)*(K2/K1)*t ; v = (1/2)*K2*t ;
%plot(t,x) ,xlabel('t~ sec'),ylabel('x(m) ');
plot(t,v) ,xlabel('t~ sec'),ylabel('v(m/s)')
DE 9 and 10 of Table III-1
The last two DE In table III -1 recquire from the analytical approach more detailed methods of solution , there are quite a few if’s .
We first post the analytical answers and later in a single numerical stroke deal with the two of them.
d 2C /dt 2 + K1 dC/dt + K2C = 0. (14)
It is a homogenous DE. It has two complementary solutions. Requires two IC , C(0) , dC(0)/dt .
Postulating C ~ exp( λt ) one reaches the quadratic characteristic equation
λ2 + K1 λ + K2 = 0 . It has two roots
λ1 = (1/2) { -K1 + ( K12 – 4 K2)1/2 } ; λ2 = (1/2) { -K1 - ( K12 – 4 K2)1/2}
The general solution is
C(t) = A 1 exp(λ1t ) + A 2 exp(λ2t ) .
If K12 – 4 K2 >0 ,the roots are real and you have growing and decaying exponentials .
If K12 – 4 K2 < 0 then the radicals are imaginary. The solution will be the product of a real exponential exp(-K1/2) and oscillatory sines and cosines.
If K12 – 4 K2 = 0 , λ1 = λ2 = - K1/2 . Another trick has to be invoked to get two independent solutions and the general solution is
C(t) = A1 exp(λ1 t) + A2 t exp(λ1 t) .
Eq 10 of Table III-1.
d 2C /dt 2 + K1 dC/dt + K2C = K3 . (15)
It is a inhomogenous DE. It has two complementary solutions and a particular solution. Requires two IC , C(0) , dC(0)/dt .
The analytical solutions are
Cparticular = K3/K2 , Ccomplementary = = A 1 exp(λ1t ) + A 2 exp(λ2t ) ,
or the complete solution
C(t) = K3/K2 + A 1 exp(λ1t ) + A 2 exp(λ2t ) (16)
when K12 ≠ 4 K2 . If K12 = 4 K2 the solution is
C(t) = K3/K2 + A 1 exp(λ1t ) + A 2 t exp(λ1t ) . (17)
Numerical integration of d 2C /dt 2 + K1 dC/dt + K2C = K3
Let K3 =0 , tscale1=1 sec , ttscale2 =0.1 sec , k1 = 1/tscale1= 1 sec-1 ,
k2 = 1/(tscale2)2 = 100 sec-2
The solution when K3 =0 , was shown to be
C(t) = A 1 exp(λ1t ) + A 2 exp(λ2t ) .
Let A1 = A 2 =1 . Then the initial conditions are ,
C(0) =2 , dC(0)/dt = λ1 + λ2 . (when λ1 ≠ λ2 ) .
The FORTRAN code given next integrates the DE by finite differences.
The solution is (with K3=0) ,
Cn = 2 Cn-1 - C n-2 + (∆t)2{- K1 (Cn-1- Cn-2)/ ∆t -K2 Cn-1 }
Fig 3.6 . Homogenous DE
FORTRAN code
c solutions of 2nd order DE by finite differences
real k1,k2,k3
data tscale1,tscale2 /1., .1 /
complex c0,c1,c2,dcdt, a ,al1,al2
k3(t)=0.
c initial conditions
c0=cmplx(2.,0.)
k1=1./tscale1
k2=1./tscale2**2
c al1 , al2 ~ lambda1 , lambda 2
if(k1**2-4.*k2.ge.0.)al1= .5*(k1+(k1**2-4.*k2)**.5)
if(k1**2-4.*k2.ge.0.)al2= .5*(k1-(k1**2-4.*k2)**.5)
if(k1**2-4.*k2.lt.0.)al1= cmplx(.5*k1,.5*(abs(k1**2-4.*k2))**.5)
if(k1**2-4.*k2.lt.0.)al2= cmplx(.5*k1,-.5*(abs(k1**2-4.*k2))**.5)
dcdt=al1+ al2
dt=amin1(tscale1,tscale2)/100.
tf=2.*amax1(tscale1,tscale2)
c1=c0+dt*dcdt
nstep=int (tf/dt)
kp=int(float(nstep)/80.)
kount=kp
print 100,0.,real(c0),aimag(c0)
do 10 i=2, nstep
t=dt*float(i)
a= k3(t-dt)-k1*dcdt-k2*c1
c2=2.*c1-c0+dt**2*a
if(i.eq.kount)then
print 100 , t ,real(c2) , aimag(c2)
kount=kount+kp
end if
dcdt=(c2-c1)/dt
c0=c1
c1=c2
10 continue
100 format(1x,'t,Real(c2) Aimag(c2) =',3(3x,e10.3))
stop
end
Another example: Oscillations for the inhomogeneous DE with K2=100/s2 , K3= 100m/s2
Fig . 3.7 Inhomogeneous DE. . Lim x t→∞ = K3/K2 =1 .
c solutions Of 2 ordee DE by finite differences
real k1,k2,k3
data tscale1,tscale2 /1., .1 /
complex c0,c1,c2,dcdt ,a
k3(t)=amp
c initial conditions
c0=cmplx(2.,0.)
dcdt=0.
k1=1./tscale1
k2=1./tscale2**2
dt=amin1(tscale1,tscale2)/100.
tf=2.5*amax1(tscale1,tscale2)
c1=c0+dt*dcdt
amp=.5*real(k2*c0)
Print*,'k1,k2, k3=',k1,k2,amp
print*,' '
nstep=int (tf/dt)
kp=int(float(nstep)/80.)
kount=kp
print 100,0.,real(c0),k3(0.)/amp
do 10 i=2, nstep
t=dt*float(i)
a= k3(t-dt)-k1*dcdt-k2*c1
c2=2.*c1-c0+dt**2*a
if(i.eq.kount)then
print 100 , t ,real(c2) , k3(t)/amp
kount=kount+kp
end if
dcdt=(c2-c1)/dt
c0=c1
c1=c2
10 continue
100 format(1x,'t,Real(c2),k3/amp =',3(3x,e10.3))
stop
end
No linear DE
Little if any of the analytical techniques used for the DE in table III is of help in the case of nonlinear equations. However the method of finite difference is universal , no matter what type of DE you confront.
Examples of nonlinear DE
d 2 C/dt2 = k C(t) 2 . ( 18 )
It is nonlinear because the dependent variable , C(t) appears squared.to be specific suppose C ~meters t~ seconds then k~ 1/( meter-sec2)
Requires two IC. C(0) , dC(0)/dt.
The generic solution is
Cn = 2Cn-1 - Cn-2 + (∆t) k ( Cn-1 )2 . (19)
Let C(0)= 1 m and dC(0)/dt=0 m/s , k= 1/(m s2 ) .
Equation (18) represents an acceleration that grows with the square of the position C~x(t).
Fig. 3.8 Solution of nonlinear DE.
MATLAB code
x(1)=1; K =1 ; v0=0 ;
tscale= 1/(K*x(1))^.5 ;tfinal=tscale;nstep=3000;
dt=tfinal/nstep;
x(2) =x(1) + dt*v0 ;
n=1;
for t=[2*dt:dt:tfinal];
n=n+1;
x(n+1)=2*x(n)-x(n-1)+dt^2*K*x(n)^2 ;
end
t=[0:dt:tfinal];
plot(t,x), xlabel('tsec'),ylabel('x~m');