Unit 14 Numerical Solution of Ordinary Differential Equations
Unit 14 Numerical Solution of Ordinary Differential Equations
ORDINARY DIFFERENTIAL
        EQUATIONS
Structure
14.1 Introductio~l
       Objectives
14.2 Basic Co~lcepts
14.3 Taylor Series Method
14.4 Euler's Method
14.5 Richardson's Extrapolation
14.6 Sultnnary
14.7 SolutionslPu~swers
14.1 INTRODUCTION
In the previous two units, you have seen how a colnplicated or tabulated function call be
replaced by an approxunati~~g      poly~~omial so that the fundatnental operations of calculus
v i ~ . ,differentiation and integratio~~
                                        can be performed more easily. 111 this w i t we shall solve
a differential equation, that is, we shall find tl~cUILIUIOWII function which satisties a
colnbinatio~lof the indeyc~ldentvariable, dependenr variable and its derivatives. In physics,
eagiaeering, chetnistry and nlany other disciplines it has become necessary to build
nlathelnatical ~nodelsto represent complicated processes. Differential equations are one of
the rnost hnportant mathematic-al tools uscd hl nlodelling problenls in the engineering and
physical sciences. As it is 1101always possible to obtain the analytical solutio~lof differential
equations recourse must necessarily be niade to numerical rtlethods for solving differential
equations. In this unit, we shall h~troduretwo such methods namely, Euler's ?lethod and
Taylor series method to obain numerical solutio~~      of ordinary differential equations (ODEs).
We shall also introducx Richardsoa's extrapolatio~lmethod to obtain higher order solutio~ls
to ODEs using lower order nlethods. To begin with, we shall recall few basic concepts fro111
the theory of differential equations which we shall be refeml~gquite often.
Objectives
After studying this unit you should be able to :
   identify the initial value problem for the first order ordinary differential equations;
   obtain the solution or the initial value proble~nsby using Taylor series method and
   Euler's method;
   use Richardson's extrapolatio~~  technique. for ilnprovi~lgthe accuracy df the result
   obtained by Euler's method. '
For example,
Numerical Differentiation Integration
and Solution of Differential Equations
                                         Differential equations of the form (I), involving derivatives w.r.t. a single independent
                                         variable are called ordinary differential equations (ODEs) whereas, those involving
                                         derivatives w.r.t. two or more i n d e p e ~ ~ d variables
                                                                                          e~~t      are partial differential equations
                                         (PDEs). Eqn. (2) is an example of PDE.
                                         Definition : The order of a differential equation is the order of the highest order                         1
                                                                                                                                                     !
                                         derivative appearing in the equation a ~ t dits degree is the highest exponent of the
                                                                                                                                                     *
                                                                                                                                                     1
                                         highest order derivative after the equation has been rationalised i.e., after it has been                   I
                                         expressed in the forn~free from radicals and any fractional power of the derivatives or
                                         negative Dower. For cxamale eauation
                                         Definition : When the depew'ent variable and its derivatives occur in the first degree
                                         only and not as higher powers or products, the equation is said to be linear, otherwise
                                         it is nonlinear.                                                                                            1
                                         Equ                      y = x2 is a linear ODE, whereas, ( x + ~ ) '           nonlinear ODE.
The general form ,of a linear ODE of order n can be expressed in the form
                                         ~ [ y =] a, (t)     (t) + a l (t) y(n-')(t) + . . . . . . + a,-, ( 9 Y (t) + an (t) Y (t) = r(t)      (5)
                                         where r(t), a, (t), i = 1, 2, . . . . . . ., n are known functions of t and
                                                                                                         I
                                                                                                             d
                                                          (t) -1
                                                               dtn-
                                                                                     ..
                                                                    + . . . . . . . + a n - l ( t ) z + an(t),
                                         is the linear differential operator. The general nonlinear ODE of order                 11   can be
                                         written as
                                         The general solution of an nth order ODE contains n arbitrary constants. In order to
                                         determine these arbitrary constants, we require n conditions. If these conditions are
                                         given at one point, then these conditions are known as initial conditions a ~ l dthe
                                         differential equation together with the initial conditions is called an initial value
                                         problem (IVP). The nth order IVP can be written as
    Heocc, it is sufficient to study numerical ~rlcthodsfor the solution of the first order IVP.
                 Y' = f(t, Y) ~ ( t " )= Yo
                                          9                                                       (10)
    The vector fonlr of thrse neth hods can their be used to solve Eqn. (9). Before attempting to
    obtain ~lunleriralsoIutions to Eqn, (lo), we lust make sure (hat the proble~nhas a unique
    solution. Thc Lollowhg (heorem cllsures the existence and uniqueiress ot (he solution to IVP (10).
    then for any yo, the IVP (10) has a unique solution. This condition is called the
    Lipschitz condition and L is called the 1,ipschitz constant.
    W e assunre the existence and uniqueness of the solution and also that f(t, y) has
    continuous partial derivatives w.r.t. t and y of a s high order as w e desire.
                                          I
    Let us assume that t b b e an interval over which the solution of the IVP (10) is
                                    [o,
    required. If w e subiivide the interval t b into
                                                                     [ n    I        11   subintervals using a steysize
. where tn = b, we obtain (he mesh points or grid points to, t,, t2 ...... 7 tn
a s shown in Fig. 1.
                                                                            Fig. 1
    W e can then write t, = to + kh, k = 0, 1, . . . . . n. A numerical method for the solution
    of the IVP ( l o ) , will produce approximate values y, at the grid points tk.
Numerical    :renliation Integration   Remember that the approximate values yk may contain the truncation and round-off errors.
and Solutb    Differential Equations
                                       We shall now discuss the cot~structionof numerical methods and related basic concepts
                                       with reference to a simple ODE.
                                       where     = a and b + Nh = b.                                                                  1
                                       Separating the variables and integrating, we find that the exact solution of Eqn. (11) is
                                                          eY' - I',)
                                              ~ ( t =) Y(Q                                                                    (12)
                                       In order to obtain a relation connecting two successive solution values, we set t = t,
                                       and gtl in Eqn. (12). Thus we get
y ( t , + ~= y(b) eYt- 1
Dividing, we get
Hence we have
                                                                y(t,,),n=O, 1,.....,N-1
                                                           Ah
                                               y(t,+J=e
Eqn. (13) gives the required relation between y(t,,) and y(t,,,).
                                                                  .
                                       Setting 11 = 0,1, 2, . . . ., N-1, successively, we can find y(tl), Y(~J,,.. . ., Y ( ~ N )
                                       from the given value y(fo).
                                                                                                                                      I
                                       An approximate methgd or a numerical method can be obtained by approximating ehh in
                                       Eqn. (13). For example, we may use the followiilg polynomial ayroxhations.                     I
and so on.
                                       Let us retain (p+l) terms in the expansion of ehy and denote the approximation to ehy
                                       by E(Ah). Tbe ~~umerical method for obtaining the approximate values y, of y(t,) can
                                       then be written as
                                               TE = ~ ( h t , -Yntr
                                                              )
                                       Sinc'e (p+l) tenns are re.tained in the expa~lsionof ehh,we have
                                                                                                     Numerical Solution of Ordinary
                                                                                                              Differential Equations
The TE is of order p+l. The integer p is then called the order of the method.
We say that a ~~umerical          method is stable if the error at ally stage, i.e. y, - y($) = E,
remains bounded as n + m. Let us examine the stability of the numerical method (17).
Putting Y , + ~ = ~ ( t , + ~+ )E , + ~ and y, = y(t,) + en in Eqn. (17), we have
We note from Eqn. (18) that the error at t,,, consists of two parts. The first part E [ U ]
- eAhis the local truncation error and can be made as small as we like by suitably
deterttlit~illgE [ U ] . The second part ( E ( U )1 E, is the propagation error from the
previous step t, to $+, and will not grow if IE(&h)l < 1. If I E ( U ) ( c 1, then as
11 + m the 'propagation error tends to zero and method is said to be absolutely stable.
Formally we give the followil~gdefinition.
Defintion : A numerical method (17) is called absolutely stable if ( E ( U ) ( s 1.
You may also observe here that the exact value y(t,) given by Eqn. (13) increases if h
> 0 and decreases if h c 0, with the growth factor eAh.The approximate value yn given
by Eqn. (17) grows or decreases with the factor IE(U) 1. Thus, in order to have
mea~li~lgful~~umericalresults, it is necessary that the growth fact-r 6f the nu~tierical
method should not increase faster than the growth factor of exact solution when h > 0
and should decay at least as fast as the growth factor of the exact solution when h c 0.
Accordilrgly, we give here the following definitioa.
Defil~itiol~
           : A numerical method is said to he relatively stal,le if IE ( U ) l s eAh,h > 0.
Firstorder: I l + h h l 4 1
or-1s U s 1
or-2s hhs 2
Second order :       I   1 + hh   +-
                                   h:2(sl
                                        The second condition gives -2 a hh. Hence the right inequality gives -2    5 Ah   a 0 . The
                                        left inequality gives
                                                                       h2h2 h3h3
                                        Thirdorder:        ( + A h + - + -2l a 1 6
-2.5rhhsO.
                                        Numerical methods for finding the solution of IVP given by Eqn. (10) may be broadly
                                        classified as
                                          i) Singlestep methods
                                         ii) Multistep methods
                                        Singlestep methods ehable us to find y,+,, an approximation to y(t,+,), if y,, y,' and h
                                        are known.
                                        If yn+l can be determined from Eqn. (19) by evaluatilig the right hand side, then the
                                        singlestep method is known as an explicit method, otherwise it is known as an implicit
                                        method. The local truncation error of the method (19) is defined by
                                        Let us now take up an example to understand how the singlestep method works.
                                        Example 1 : find the solution of the IVP y' = hy, y(0) = 1 in 0 < t a 0.5, using the
                                        first order method
                                                          -
                                               Y , + ~ = (1 t Ah) yn with h = 0.1 and h = 1.
                                                                                               0.5   0.5
                                        Solution : Here the number of intervals are N      = -=      -=
                                                                                                h    0.1
                                               We have yo = 1
       11, =(1 t Ah) y,=(l t A h ) = ( l tO.1A)                                            Numerical Solution of Ordinary
                                                                                                    Differential Equations
       y 2 = (1 t Ah) y, =(1 tAh12 =(1 t 0 . 1 ~ ) ~
'We now give in Table 1 the values of y, for h = * 1 together with exact values.
In the same way you can obtain the solution using the second order method and
compare the results obtained in the two cases.
               (
       y,+, = 1 t Ah t - y, with h = d.1 and A = 1.
                       *
We are now prepared to consider numerical methods for integrating differential
equations. The first method we discuss is the Taylor series method. It is not strictly a
numerical method, but it is the most fundamental method to which every numerical
method must compare.
                                                                        h
                                          where   + (tk, yk, h) = y; + 21 yrlk+ . . . + -yk
                                                                                        h ~ - l (P)
                                                                                         P!
                                          This is called the Taylor Series method of order p. The truncation error of the lllethod
                                          is given by
                                                         hp+l
                                                      --
                                                      -           y(P+l)   (tk + Oh), 0 < 0 < 1
                                                        (p+l) !
                                                  Y' =f(t,y)
                                                  yl'=f,+f$
                                                            h3
                                          with t h e T E = - y U ' ( a ) , t,, < a < 4,'
                                                             6
                                          The Taylor series method of 0(h3), (p=3) is
                                                             h4
                                           with the TE =          y(P) (a), t,, c a <   C+,.
                                           Example 2 : Using the third order Taylor series method find the solution of the
                                           differential equation
                                           Solution : We have the derivatives and their values at x=2, y=2 as follows :
                                                   yr=l-Y                                         y' (2) = 0
                                                             X
                                                                                     Numerical Solution of Ordinary
                                                                                              Differential Equations
Example 3 : Solve the equation x2y' = 1 - xy - x2y2, y(1) = - 1 from x=l to x=2 by
using Taylor series method of 0(h2) with h =      ln
                                              and 114 and find the actual error at
x=2 if the exact solution is y = - llx.
Since the exact value is y(2) = -0.5, we have the actual errors as
                            1
       e, = 0.0583 with h = -
                            3
                            1
       e2 = 0.0321 with h = -
                            4
- - - - - - -
Write the Taylor series method of order four and solve the IVPs E2) and E3).
E4) Using second order Taylor series method solve the IVP
                                           Notice that though the Taylor series method of order p gives us results of desired accuracy
                                           in a few number of steps, it requires evaluation of the higher order derivati. and becomes
                                           tedious to apply if the various derivatives are complicated. Also, it is d f l .*'t to determine
                                           the error in such cases. We now consider a method, the Eukr's metbod -.,hichcan be
                                           regarded as Taylor series method of order one and avoids tbese difficulties.
                                           Let
                                                 [t o ,bI be the interval over which the solution of the given IVP is to be
                                           determined. Let h be the steplength. Then the nodal points are defined by tk = to + kh,
                                           k = 0 , 1 , 2 ,........, N w i t h t N = t O + N h = b .
Fii. ;I
Yk+l = Yk +h ~ ' k
                                           with       TE =
                                                                 (TI
                                                                 - y"(a),       tk < a < t,
                                           Eqn. (32) is known as the Euler's method and it calculates recursively the solution at
                                           the nodal points &, k = 0, 1,      ., N... . ..
                                           Since the truncation error (31) is of order b2, Euler's method is of first order. It is also
                                           called an 0@)method.
      Let us now see the geometrioal representation of the Euler's method.                                   Numerical Solution of Ordinary
                                                                                                                      Differential Equations
      Geometrical Interpretation
      Let y(t) be the solution of the given IVP, Integrating
      we get
                                                             dt
                                                                         *
                                                                = f(t, y) from tk to tk+l,
      (y$    dt =   j'.'f(t,y) dt
                        tk
                                    =~    ( t ~-+y(tk)
                                                   ~)                                                 (33)
      We know that geometrically f(t, y) represents the slope of the curve y(t). Let us
      approximate the slope of the curve between $ and tk+,by the slope at $ only. If we
      approximate y(tk+J and y(tJ by yk+, and yk respectively, then we have
      Example 4 : Use Euler method to find the solution of y' = t              + ly 1,   given y(0) = 1.
      Find the solution on [O, 0.81 with h = 0.2.
                    .
      Solution : We have
Numerical Differentiation Integration
and Solution of Differential Equations     y(0.81= Y, = Y, + (0.2) f,
                                                      = 1.856 + (0.2) [0.6 + 1.856)
                                                      = 2.3472
                                         Example 5 : Solve the differential equation y' = t+y, y(0) = 1. t E [0,1] by Euler's
                                         method using h = 0.1. If the exact value is y(1) = 3.436564, find the exact muT.
                                         Solution : Euler's inethod is
                                                              = Yn + hy',
                                                           Yn+,
                                         For the given problem, we have
                                                               = (1 + h)yn + ht,
                                                          h = O.l,y(O) = 1,
                                                          y1 = yo = (1 + 0.1) + (0.1) (0) = 1.1
                                                          y2=(1.1)(1.1)+(0.1)(0.1)= 1.22, y,= 1.362
                                                          y, = 1.5282, y, = 1.72102, y, = 1.943122,
                                                          y, = 2.197434, y, = 2.487178, y, = 2.815895
                                                          ylo = 3.187485 = y(1)
                                                          actual error = y(1) - ylo = 3.436564 - 3.187485 = 0.2491.
                                         Remark':    Since  Euler's method is of O(h), it requires h to be very small to attain the
                                         desired accuracy. Hence, very often, the number of steps to be carried out becomes very
                                         large. In such cases, we need higher order methods to obtain the required accuracy in a
                                         limited number of steps.
                                         This equation is called the difference equation associated with Euler's method. A
                                                                                                        . . . . . . ., Y , + ~ soine
                                         difference equation of order N is a relation involving y,, Y,+~,
                                         simple difference equations are
where n is an integer.
                                         where the coefficients aN-,, aN,,   .....   . . ., a. and b may be functions of n but not of
                                         y. All the Eqns. (35) are linear. It is easy to solve the difference Eqn. (36), when the
                                         coefficients are constant or a function of n say linear of a quadratic function of n.
I
    setting y,(p) = A (a c o ~ ~ s t a ~inl tEqn.
                                             )    (36) and detrrmillil~gthe value of A. For detail,
    you call refer to ele~tlentarynumerical analysis by Coate- deBoor. We illustrate this
    method by considering a few examples.
    Example 6 : Find the solution of the initial-value difference equations
                  Ynt2-4~,+, + 3yn = 2", Yo = 0, Y, = 1
I
t
        .-. y, (C) = C, (1)" + C, (3)"
    This gives
                 = C,   + 3"C2
    For obtaining the particular solutio~lwe try y,(p) = A2".
I or, A = -1
          ydc) = C (1 + 3h)'.
    For obtaining the particular solution we try yk(p) = Ah.
    This give
Numerical   DifferentiationIntegration   Therefore, the general solution of the given problem is
and Solution of Differential Equations
                                                                 5
                                               y k = ~+3l1)~--.
                                                          (l
                                                                 3
                                         Using the condition y(0) = 1, we obtain C = 8/3.
                                         Thus
                                                     8             5
                                                y k = y (1 +3i1)~--.
                                                                   3
Yk+l = + 3h) yk + 5h
                                         E6)    y' =   -
                                                       x2 - 4y
                                                              Ly(4) ,
                                                                    = 4. Find y(4.l) taking h = 0.1
E8) y' = 1 + y2, y(0) = 1. Find y(0.6) taking h = 0.2 and h = 0.1.
                                         which is au O(h) method. Let F(h) and F(h/2) be the solutions obtained by using step
                                         lengths h and h/2 respectively.
                                         Recall that the Richardson's extrapolation method of combining two computed values with two
                                         different step sizes, to obtain a higher order method is given by (ref.. Formula (54) Unit 12)
Let y1(4) and y2(t,J be the two values obtained by a numerical method of order p with
step sizes hl and b. If e: and 5 are the corresponding errors, then
                                    P
in the interval (0, I ] taking h = 0.2, 0.1. Using Richardson's extrapolation technique
obtain the improved value at t = 1.
                                        Table 2 : h    = 0.2
                                         -    -    -
Similarly, starting with t, = 0, yo = 1, we obtain the following table of values for h = 0.1.
Table 3 : h = 0.1
                                                               2F(0.1) - F(0.2)
                                                     (0.1) =
                                                                   1
                                                = 2(0.50364) - (0.50706)
                                                = 0.50022
                                         Example 7 : Use Euler's method to solve numerically the initial value: proble
                                         y' = t t y , y(0) = 1 with h = 0.2, 0.1 and 0.05 in the interval [0, 0.61. Apply
                                         Richardson's extrapolation technique to compute y(0.6).
Table 4 : h = 0.2
Table 5 : h = 0.1
Table 6 : h = 0.05
                      = 2.043649
        The exact solution is y = - (l+t) + 2et
i
I
               = 0.000589
        Aad now a few exercises for you
E9) The IW
                 is given. Find (0.6) with h = 0.2 and h = 0.1, using EulerTs method and extrapolate
                 the value y(0.6). Compare with the exact solution.
1       E10) Extrapolate the value y(0.6) obtained in B).
i
We now end this unit by giving a summary of what we have cwered in it.
is given by
              2) Euler's method is the Taylor series method of order one. The steps involved in
                 solving the IVP given by (10) by Euler's method are as follows :
                 Step 1 : Evaluate f(t,,, yo)
                 Step 2 : Find y, = yo + h f(t, yo)
                 Step 3 : If t,, < b, change b to b + h and yo to y, and repeat steps 1 and 2
Numerical Differentiation Integration
and Solution o f Differential Equations         Step 4 : If t, = b, write the value of y
                                            3) Richardson's extrapolation method given by Eqn. (42) can be used to improve the
                                               values sf the function evaluated by the Euler's method.
                                                       y2= (1.105)~
                                                       ------------------
                                                    ys = (1.105)~
                                                Table giving the values of y, together with exact values is
                                                                                   Table 7
                                                 t                 Second order method         Exact solution
                                                 0                 1                            1
                                                 0.1               1.105                        1.10517
                                                 0.2               1.22103                      1.22140
                                                 0.3               1.34923                      1.34986
                                                 0.4               1.49090                     .1.49182
                                                 0.5      s        X 34745                    ' 1.64872
Substituting
Table 8 : h = 0.2
Table 9 : h = 0.1
                                    Table 10 : h = 0.2
Numerical Diff6rentiation Integration                                         Table 11 : h   = 0.1
and Solution of Differential Equations