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:

reibaretti2004@yahoo.com

 

 

 

 

 

 

 

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

 

 

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');