Numerical Solution of

Ordinary Differential Equations

(Part - 2)

P. Sam Johnson

May 3, 2020

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 1/55
Overview : Runge-Kutta Methods

The Taylor’s series method of solving differential equations numerically is

handicapped by the problem of finding the higher order derivatives.

Euler’s method is less efficient in practical problems since it requires h to

be small for obtaining reasonable accuracy.

The Range-Kutta methods do not require the calculations of higher

order derivatives and they are designed to give greater accuracy with the
advantage of requiring only the function values at some selected points on
the sub-interval.

These methods agree with Taylor’s series solution upto the terms of hr
where r is the order of the Runge-Kutta method.

In the lecture, we discuss Runge-Kutta methods of orders upto 4.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 2/55
Overview : Predictor-Corrector Methods

A predictor-corrector method is an algorithm that proceeds in two steps.

First, the prediction step calculates a rough approximation of the desired
quantity. Second, the corrector step refines the initial approximation using
another means.

This is not unlike an implementation of Newton-Raphson. In that method

we require an initial guess (we “predict”) and then the Newton-Raphson
approach tells us how to iterate (or “correct”) our latest approximation.

The main difference here is that we have a systematic way of obtaining

the initial prediction.

We discuss some predictor-corrector methods in the lecture.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 3/55
Introduction

Runge-Kutta method is a numerical technique used to solve ordinary

differential equation of the form
dy
= f (x, y ) given y (x0 ) = y0 .
dx
Runge-Kutta methods (also called as RK methods) are the generalization
of the concept used in Modified Euler‘s method. They are iterative
methods to calculate the solution of a differential equation. Starting from
an initial condition, they calculate the solution forward step by step.

The RK methods were first studied by the German mathematicians Carle

Runge and Martin Kutta around 1900. Modern developments are mostly
due to John Butcher in the 1960s.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 4/55
Introduction

The basic reasoning behind Runge-Kutta methods is outlined in the

following.

In Modified Euler’s method, the slope of the solution curve has been
approximated with the slopes of the curve at the end points of the each
sub interval in computing the solution.

The natural generalization of this concept is computing the slope by

taking a weighted average of the slopes taken at more number of
points in each sub interval. However, the implementation of the scheme
differs from Modified Euler’s method so that the developed algorithm is
explicit in nature. The final form of the scheme is of the form

yi+1 = yi + (weighted average of the slopes) for i = 0, 1, 2, . . .

where h is the step length and yi and yi+1 are the values of y at xi and
xi+1 respectively.
P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 5/55
Introduction

In general, the slope is computed at various points x in each sub interval

[xi , xi+1 ] and multiplied them with the step length h and then weighted
average of it is then added to yi to compute yi+1 .

Thus the RK method with r slopes called as r -stage RK method can be

written as
K1 = hf (xi , yi )
K2 = hf (xi + c2 h, yi + a21 K1 )
K3 = hf (xi + c3 h, yi + a31 K1 + a32 K2 )
.. . .. ..
. = .. . .
Kr = hf (xi + cr h, yi + ar 1 K1 + ar 2 K2 + · · · + arr −1 Kr −1 )
and
yi+1 = yi + (W1 K1 + W2 K2 + · · · + Wr Kr ) for i = 0, 1, 2, . . . .
P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 6/55
Runge-Kutta Method of Order p

To determine the parameters

ci ’s for 2 ≤ i ≤ r (coefficients for x)

aij ’s for 1 ≤ j < i ≤ r (coefficients for y )
Wi ’s for 1 ≤ i ≤ r (weights)

in the above equation, yi+1 defined in the scheme is expanded interms of

step length h.

The resultant equation is then compared with Taylor series expansion of

the solution of the differential equation upto a certain number of terms,
say p.

Then the r -stage RK method will be of order p or is an pth order

Runge-Kutta method (called Runge-Kutta method of order p).

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 7/55
The Butcher Tableau

To specify a particular method, we need to provide the parameters :

r (number of stages)
ci ’s for 2 ≤ i ≤ r (coefficients for x)
aij ’s for 1 ≤ j < i ≤ r (coefficients for y )
Wi ’s for 1 ≤ i ≤ r (weights).

These data are usually arranged in a co-called Butcher tableau.

0
c2 a21
c3 a31 a32
.. .. .. ..
. . . .
cr ar 1 ar 2 ... arr −1
W1 W2 ... Wr −1 Wr
P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 8/55
Runge-Kutta Method of Order 1

Consider
dy
= f (x, y ) given y (x0 ) = y0 .
dx
The Euler’s formula for first approximation to the solution of the above
differential equation is given by

yi+1 = yi + hf (xi , yi ) = yi + hy 0 (xi ) [ since y 0 = f (x, y ) ].

Also
h2 00 h2
yi+1 = yi + hy 0 (xi ) +
y (xi ) + y 000 (xi ) + · · · .
2! 2!
Clearly the Euler’s method agrees with the Taylor’s series solution upto the
term in h.

Hence Euler’s method is the Runge-Kutta method of first order.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 9/55
Runge-Kutta Method of Order 2

Now, consider the case r = 2 to derive the 2-stage (second order) RK

method. For this

K1 = hf (xi , yi )
K2 = hf (xi + c2 h, yi + a21 K1 )
yi+1 = yi + (W1 K1 + W2 K2 ) for i = 0, 1, 2, . . . .

Now by Taylor series expansion

h2 00 h2 000
y (xi+1 ) = y (xi ) + hy 0 (xi ) + y (xi ) + y (xi ) + O(h4 )
2! 2!
h2 h3
yi+1 = yi + hfi + (fx + fy f ) + (fxx + 2fxy f + fyy f 2 + fy (fx + fy f )) + O(h4 )
2! 3!

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 10/55
Runge-Kutta Method of Order 2
Also
K1 = hfi
K2 = hf (xi + c2 h, yi + a21 K1 ) = hf (xi + c2 h, yi + a21 hfi )
(c2 h)2 (a21 hfi )2
 
= h fi + c2 hfx + a21 hfi fy + fxx + fyy + (c2 h)(a21 hfi )fxy + O(h4 )
2! 2!
(by Taylor’s series expansion of two variables)
h3
yi+1 = yi + (W1 + W2 )hfi + h2 (W2 c2 fx + W2 a21 ffy ) + W2 (c22 a21 ffxy + a21
2 2
f fyy ) + O(h4 )
2

Now by comparing the equal powers of h in yi+1 of two equations, we get

1 1
W1 + W2 = 1, c 2 W2 = and a21 W2 = .
2 2
The solution of this system is
1 1
a21 = c2 , W2 = and W1 = 1 − .
2c2 2c2
where c2 is any arbitrary constant not equal to zero.
P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 11/55
Runge-Kutta Method of Order 2

Since 2-stage RK method compares with Taylor series upto h2 for any
value of c2 , the 2-stage RK method is of order two and hence this scheme
is denoted in many text books as a second order RK method.

For these values of a21 , W2 , W1 , there exists an infinite number of

2-stage Runge-Kutta methods of order 2.

Now, to give some numerical values to a21 , W2 , W1 ; first the value of c2

has to be fixed.

Generally the value of c2 is fixed such that the values of a21 , W2 , W1 are
integers or some real numbers which easy to remember.
1
Two of such cases, the values of c2 are 2 and 1 respectively.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 12/55
Runge-Kutta Method of Order 2

Let c2 = 1, a21 = 2, W2 = W1 = 12 .

The corresponding 2-stage (second order) RK method is

K1 = hf (xi , yi )
K2 = hf (xi + h, yi + K1 )
K1 + K2
yi+1 = yi + for i = 0, 1, 2, . . . .
2

or equivalently
 
h
yi+1 = yi + f (xi , yi ) + f (xi + h, yi + hf (xi , yi )) for i = 0, 1, 2, . . .
2

which is nothing but the Modified Euler’s method.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 13/55
Butcher Tableau for Runge-Kutta Method of Order 2

Here are three methods based on convenient choices of c2 . Note that the
first two methods are due to Runge.

The Butcher tableau for Runge-Kutta method of order 2 are as follows:

0 0 0
1 1 2 2
2 2 1 1 3 3

1 1 1 3
0 1 2 2 4 4

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 14/55
Runge-Kutta Method of Order 3

The thrid order Runge-Kutta method is given by

 
1
yi+1 = yi + K1 + 4K2 + K3 for i = 0, 1, 2, . . . ,
6

where

K1 = hf (xi , yi )
h K1
K2 = hf (xi + , yi + )
2 2
K3 = hf (xi + h, yi − K1 + 2K2 )

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 15/55
Butcher Tableau for Runge-Kutta Method of Order 3

The Butcher tableau for Runge-Kutta method of order 3 are as follows:

0 0 0
1 1 1 1
2 2 1 1 3 3
1 1 1
1 -1 2 2 4 4 1 −1 2

1 2 1 1 1 2 3 1
6 3 6 6 6 3 0 4 4

0 0 0
2 2 2 2 2 2
3 3 3 3 3 3
2 2 2 1 1
3 0 3 3 3 3 0 −1 1

1 3 3 1 3 3 1
4 8 8 4 0 4 0 4 4

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 16/55
Runge-Kutta Method of Order 4

The fourth order Runge-Kutta method is given by

 
1
yi+1 = yi + K1 + 2K2 + 2K3 + K4 for i = 0, 1, 2, . . . ,
6

where

K1 = hf (xi , yi )
h K1
K2 = hf (xi + , yi + )
2 2
h K2
K3 = hf (xi + , yi + )
2 2
K4 = hf (xi + h, yi + K3 )

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 17/55
Runge-Kutta Method of Order 4

In each step of fourth-order Runge-Kutta method, the derivative is

evaluated four times:

once at the initial point,

twice at trial midpoints, and
once at a trial endpoint.

From these derivatives the final function value yn+1 (shown as a filled dot)
is calculated.
P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 18/55
Butcher Tableau for Runge-Kutta Method of Order 4

The Butcher tableau for Runge-Kutta method of order 4 are as follows:

0 0 0
1 1 1 1 1 1
2 2 4 4 2 2
1 1 1 1 0 −1 1
2
0 2 2
0 2
3 1
1 0 0 1 1 1 −2 2 1 −1 2 2

1 1 1 1 1 2 1 1 2 1 1
6 3 3 6 6
0 3 6 12 3 12 6

0 0 0
1 1 2 2
1 1 3 3 3 3
1 3 1 2 −1 1 1 1
2 8 8 3 3
1 3 12 4
1 −2 −1 4 1 1 −1 1 1 − 54 1
4
2
1 1 2 1 1 3 3 1 1 3 3 1
6 12 3 12 8 8 8 8 8 8 8 8

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 19/55
Working Rule

For a given point (x0 , y0 ) in the (unknown) solution curve, we are going to
find the ordinate of given x.

Divide the interval connecting x0 and x into n subintervals each of width

h, the value of h is quite small.

To apply one of the Runge-Kutta methods, we should calculate the

weights Ki at (x0 , y0 ) to find y1 = y (x0 + h).

Similarly, the value of y in the second interval is obtained by replacing x0

by x1 and y0 by y1 in the Runge-Kutta method and we obtain y2 .

In general, to find yn substitute xn−1 , yn−1 in the expression for K1 , K2 etc.

To evaluate yn+1 , we need information only at the point yn . Information

at the points yn−1 , yn−2 etc. are not directly required . Hence
Runge-Kutta methods are step methods.
P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 20/55
Exercises

1 Compute y (0.1) and y (0.2) by Runge-Kutta method of fourth order

for the differential equation
dy
= xy + y 2 , y (0) = 1.
dx
2 Use Runge-Kutta method of fourth order to find y (0.1), given that

dy 1
= , y (0) = 1.
dx x +y
3 Given y 0 = x 2 − y , y (0) = 1, find y (0.1) using Runge-Kutta method
of fourth order.
4 Using 4th order Runge-Kutta method, find y (0.1), y (0.2) and y (0.3),
given that
dy
= 1 + xy , y (0) = 2.
dx
P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 21/55
Exercises

5 Using 4th order Runge-Kutta method, evaluate the value of y when

x = 1.1, given that
dy y 1
+ = 2, y (1) = 1.
dx x x
6 Apply third order Runge-Kutta method to find an approximate value
of y when x = 0.2, given that
dy
= x + y, y (0) = 1.
dx
7 Using Runge-Kutta method of fourth order, solve

dy y2 − x2
= 2
dx y + x2
with y (0) = 1 at x = 0.2, 0.4.
P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 22/55
Exercises

8 Apply Runge-Kutta method to find an approximate value of y when

x = 0.2 in steps of 0.1 if
dy
= x + y2
dx
given that y = 1 where x = 0.
9 Using Runge-Kutta method of fourth order, solve for y at x = 1.2, 1.4
from
dy 2xy + e x
= 2
dx x + xe x
given x0 = 1, y0 = 0.
10 Given
dy
=y −x
dx
where y (0) = 2, find y (0.1) and y (0.2) correct to 4 decimal places by
I Runge-Kutta second order formula
I Runge-Kutta fourth order formula.
P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 23/55
Exercises

11 Given
dy
= 1 + y2
dx
where y (0) = 0, find y (0.2), y (0.4) and y (0.6) by Runge-Kutta fourth
order formula.
12 Taking h = 0.5, solve the initial value problem for x = 0.05 of the
differential equation
dy y
= 3x + , y (0) = 1.
dx 2

I Euler’s method
I Modified Euler’s method
I Runge-Kutta method of order 4.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 24/55
Solving Simultaneous First Order Differential Equations :
Runge-Kutta Methods

Consider the simultaneous differential equations of the type

dy
= f (x, y , z) (1)
dx
and
dz
= g (x, y , z) (2)
dx
with initial conditions y (x0 ) = y0 and z(x0 ) = z0 can be solved by
Runge-Kutta methods.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 25/55
Starting at (xi , yi , zi ) and taking the step-sizes for x, y , z to be h, k, `
respectively, we have the following for the Runge-Kutta methods of fourth
order.
K1 = h f (xi , yi , zi ) J1 = h g (xi , yi , zi )
h K1 J1 h K1 J1
K2 = h f (xi + , yi + , zi + ) J2 = h g (xi + , yi + , zi + )
2 2 2 2 2 2
h K2 J2 h K2 J2
K3 = h f (xi + , yi + , zi + ) J3 = h g (xi + , yi + , zi + )
2 2 2 2 2 2
K4 = h f (xi + h, yi + K3 , zi + J3 ) J4 = h g (xi + h, yi + K3 , zi + J3 )
 
1
yi+1 = yi + K1 + 2K2 + 2K3 + K4 for i = 0, 1, 2, . . . ,
6
and  
1
zi+1 = zi + J1 + 2J2 + 2J3 + J4 for i = 0, 1, 2, . . . .
6

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 26/55
Exercises

13 Solve the system of differential equations

dy dz
= xz + 1, = −xy
dx dx
for x = 0.3(0.3)0.9 using Runge-Kutta fourth order formula. Initial
values are x = 0, y = 0, z = 1.
14 Using Runge-Kutta method of order 4, find the approximate values of
x and y at t = 0.2 for the following system:
dx dy
= 2x + y , = x − 3y
dt dt
given that when t = 0, x = 0, y = 0.5.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 27/55
Solving Second Order Differential Equations by
Runge-Kutta Methods

Let us describe the fourth order Runge-Kutta method for solving the
second order differential equation

d 2y
= y 00 = f (x, y , y 0 ). (3)
dx 2
If we put y 0 = z, then y 00 = z 0 . So the differential equation becomes

z 0 = f (x, y , z).

We have now the simultaneous differential equations

dy dz
=z (that is, f1 ≡ z) and z = = f (x, y , z) (that is, f2 ≡ f )
dx dx
where f1 means f1 (x, y , z) and f2 denotes f2 (x, y , z).

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 28/55
Exercises

15 Given
d 2y
− y 3 = 0, y (0) = 10, y 0 (0) = 5.
dx 2
Evaluate y (0.1) using Runge-Kutta method.
16 Use the Runge-Kutta method with fourth order accuracy to determine
the approximate value of y at x = 0.1 if y satisfies the differential
equation
d 2y dy
2
− x2 − 2xy = 1, y (0) = 1, y 0 (0) = 0.
dx dx
17 Using Runge-Kutta method, solve
 2
d 2y dy
2
=x − y2
dx dx
for x = 0.2 correct to 4 decimal places. Initial conditions are
y (0) = 1, y 0 (0) = 0.
P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 29/55
Introduction : Predictor-Corrector Methods

A predictor-corrector method is an algorithm that proceeds in two steps.

First, the prediction step calculates a rough approximation of the desired
quantity. Second, the corrector step refines the initial approximation using
another means. This involves

The predictor step. We use an explicit method to obtain an

p
approximation yn+k to yn+k .
The corrector step. We use an implicit method, but with the
p
predicted value yn+k on the right-hand side in the evaluation of fn+k .
p
We use fn+k to denote this approximate (predicted) value of fn+k .
We can then go on to correct again and again. At each step we put
the latest approximation to yn+k in the right-hand side of the scheme
to generate a new approximation from the left-hand side.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 30/55
Introduction : Predictor-Corrector Methods

What has been done in Newton-Raphson method?

This is not unlike an implementation of Newton-Raphson. In that method

we require an initial guess (we “predict”) and then the Newton-Raphson
approach tells us how to iterate (or “correct”) our latest approximation.

The main difference here is that we have a systematic way of obtaining

the initial prediction.

It is sufficient for our purposes to illustrate the idea of a

predictor-corrector method using the simplest possible pair of methods.

We use Euler’s method to predict and the trapezium method to correct.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 31/55
Predictor-Corrector Methods

Consider the differential equation

dy
= f (x, y ) given y (x0 ) = y0 .
dx
We divide the range for x into a number of step sizes of equal width h. If
xi and xi+1 are two consecutive points, then xi+1 = xi + h.

Euler’s formula for the above differential equation

yi+1 = yi + hf (xi , yi ), i = 0, 1, 2, . . . . (4)

The modified Euler’s formula is

 
h
yi+1 = yi + f (xi , yi ) + f (xi+1 , yi+1 ) , i = 0, 1, 2, . . . . (5)
2

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 32/55
Predictor-Corrector Methods

The value of yi+1 is first estimated by (4) and this value is inserted on the
right side of (5) to get a better approximation for yi+1 .

This value of yi+1 is again substituted in (5) to find a still better

approximation of yi+1 .

This process is repeated until two consecutive values of yi+1 are almost
equal.

This technique of refining an initially crude estimate by means of more

accurate formula is known as predictor-corrector methods. Equation (4)
is called the predictor and (5) is called corrector.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 33/55
Predictor-Corrector Methods

In the method described so far to solve a differential equation over an

interval only the value of y at the beginning of the interval was required.

In the predictor and corrector methods, four prior values of y are

needed to evaluate the value of y at xi+1 .

A predictor formula is used to predict the value yi+1 of y at xi+1 and

then a corrector formula is used to improve the value of yi+1 .

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 34/55
Milne’s Method

Consider the first order differential equation

dy
= f (x, y ) given y (x0 ) = y0 .
dx
Newton’s forward difference formula can be written as
n(n − 1) 2 n(n − 1)(n − 2) 3
f (x, y ) = f0 + n∆f0 + ∆ f0 + ∆ f0 + · · · . (6)
2! 3!
Substituting this in the relation
Z x0 +4h
y4 = y0 + f (x, y )dx
x0

we get
x0 +4h  
n(n − 1) 2
Z
y4 = y0 + f0 + n∆f0 + ∆ f0 + · · · dx.
x0 2!

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 35/55
Milne’s Method
Let x = x0 + nh. Therefore
Z 4 n(n − 1) 2

y4 = y0 + h f0 + n∆f0 + ∆ f0 + · · · dn
0 2!
 
20 2 8
= y0 + h 4f0 + 8∆f0 + ∆ f0 + ∆3 f0 + · · ·
3 3
 
0 0 20 2 8
= y0 + h 4y0 + 8(E − 1)y0 + (E − 2E + 1)y00 + (E 3 − 3E 2 + 3E − 1)y00
3 3
(neglecting fourth and higher order differences)
 
20 0 8
= y0 + h 4y00 + 8(y10 − y00 ) + (y2 − 2y10 + y00 ) + (y30 − 3y20 + 3y10 − y00 )
3 3
 
8 0 4 0 8 0
= y0 + h y − y2 + y3 .
3 1 3 3

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 36/55
Milne’s Method

Thus
4h  0
2y1 − y20 + 2y30 .

y4 = y0 +
3
If x0 , x1 , . . . , x4 are any 5 consecutive values of x, then the above equation
can be written as
4h  0 0
+ 2yn0 .

yn+1,p = yn−3 + 2yn−2 − yn−1
3

This is called Milne’s Predictor formula (the subscript p indicates that it

is a predicted value).

This formula can be used to predict the value of y4 when those of

y0 , y1 , y2 , y3 are known.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 37/55
Milne’s Method

To get a corrector formula we substitute Newton’s formula (6) in the

relation Z x0 +2h
y2 = y0 + f (x, y )dx
x0
and we get
Z x0 +2h

n(n − 1) 2

y2 = y0 + f0 + n∆f0 + ∆ f0 + · · · dx
x0 2!
Z 2 
n(n − 1) 2
= y0 + h f0 + n∆f0 + ∆ f0 + · · · dn putting x = x0 + nh
0 2!
 
1 2
= y0 + h 2f0 + 2∆f0 + ∆ f0 + · · ·
3
 
1
= y0 + h 2y00 + 2(E − 1)y00 + (E 2 − 2E + 1)y00
3
(neglecting higher order differences)
 
1
= y0 + h 2y00 + 2(y10 − y00 ) + (y20 − 2y10 + y00 )
3

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 38/55
Milne’s Method

Thus
h 0
y0 + 4y10 + y20 .

y2 = y0 +
3

If x0 , x1 , x2 are any three corrective values of x, the above relation can be

written as
h 0
yn−1 + 4yn0 + yn+1
0

yn+1,c = yn−1 +
3

This is known Milne’s corrector formula where the suffix c stands for
corrector.
0
An improved value of yn+1 is computed and again the corrector formula is
applied until we get yn+1 to the derived accuracy.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 39/55
Adams-Bashforth Method

Consider
dy
= f (x, y ) with y (x0 ) = y0 .
dx
Newton’s backward interpolation formula can be written as

n(n + 1) 2 n(n + 1)(n + 2) 3

f (x, y ) = f0 + n∇f0 + ∇ f0 + ∇ f0 + · · ·
2! 6
x−xn
where n = h and f0 = f (x0 , y0 ).

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 40/55
Adams-Bashforth Method

Substituting this
Z x0 +h
y1 = y0 + f (x, y )dx, (7)
x0
we get
Z x0 +h

n(n + 1) 2

y1 = y0 + f0 + n∇f0 + ∇ f0 + · · · dx
x0 2!
Z 1 
n(n + 1) 2
= y0 + f0 + n∇f0 + ∇ f0 + · · · dn (putting x = x0 + nh)
0 2!
 
1 5 2 3
= y0 + h f0 + ∇f0 + ∇ f0 + ∇3 f0 + · · ·
2 12 8

Neglecting fourth and higher order differences and expressing ∇f0 , ∇2 f0 ,

∇3 f0 in terms of function values we get
h 
55y00 − 59y−1
0 0 0

y1 = y0 + + 37y−2 − 9y−3 .
24
P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 41/55
Adams-Bashforth Formula

This can also be written as

h
55yn0 − 59yn−1
0 0 0

yn+1,p = y0 + + 37yn−2 − 9yn−3 .
4

This is called Adams-Bashforth formula and is used as a predictor

formula. The notation p indicates that it is a predicted value.

A corrector formula can be derived in a similar manner by using Newton’s

backward difference formula at f1 .

That is,
n(n + 1) 2 n(n + 1)(n + 2) 3
f (x, y ) = f1 + n∇f1 + ∇ f1 + ∇ f1 + · · ·
2! 3!

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 42/55
Adams-Moulton Corrector Formula
Substituting this in (7), we get,
Z  x1 n(n + 1) 2

y1 = y0 + f1 + n∇f1 + ∇ f1 + · · · dx
x0 2!
Z 0  
n(n + 1) 2
= y0 + h f1 + n∇f1 + ∇ f1 + · · · dn (putting x = x1 + nh)
−1 2!
 
1 1 2 1 3
= y0 + h f1 − ∇f1 − ∇ f1 + ∇ f1 − · · ·
2 12 24

Neglecting 4th and higher order differences and expressing ∇f1 , ∇2 f1 ,

∇3 f1 in terms of function values, we get
 
h
y1 = y0 + 9f1 + 19f0 − 5f−1 + f−2 .
24

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 43/55
Adams-Moulton Corrector Formula

This can also be written as

 
h p
yn+1,c = yn + 9f + 19fn − 5fn−1 + fn−2 .
24 n+1

This is called Adams-Moulton corrector formula.

The subscript c indicates that the value obtained in the corrected value
and the superscripts p on the right indicates that the predicted value of
yn+1 should be used for computing the value of f (xn−1 , yn−1 ).

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 44/55
Important Note

To apply a predictor-corrector method, we need four starting values of y

which can be calculated by any of the numerical methods.

In practice, the 4th order Runge-Kutta method is found to be most useful.

Once the predictor-corrector method starts, it is very effective, as past
information is effectively used.

To apply Milne’s predictor-corrector method, we require 4 prior values of y .

If these values are not given, we can determine them by using Taylor’s
series, Euler’s or Runge-Kutta methods.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 45/55
Exercises

18 Given
dy 1
= , y (0) = 2.
dx x +y
If y (0.2) = 2.09, y (0.4) = 2.17 and y (0.6) = 2.24, find y (0.8) using
Milne’s method.
19 Using Milne’s predictor-corrector formula, find y (0.4), for the
differential equation
dy
= 1 + xy , y (0) = 2.
dx
20 Given
dy 1
= (1 + x 2 )y 2
dx 2
and y (0) = 1, y (0.1) = 1.6, y (0.2) = 1.12, y (0.3) = 1.21. Evaluate
y (0.4) by Milne’s predictor-corrector method.
P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 46/55
Exercises

21 Find by Milne’s method for the equation

y 0 = y − x 2, y (0) = 1,

by obtaining the starting values by Taylor’s series method.

22 Using Adams-Bashforth method find y (4.4) given

5xy 0 + y 2 = 2, y (4) = 1, y (4.1) = 1.0049, y (4.2) = 1.0097

and y (4.3) = 1.0143.

23 Using Adams-Bashforth method determine y (1.4) given that

y 0 − x 2y = x 2, y (1) = 1.

Obtain the starting values from the Euler’s method.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 47/55
Exercises

24 Using Adams-Bashforth method find y (0.4) given that

y 0 = 1 + xy , y (0) = 2.

25 Apply Milne’s method to find y (0.4) of the initial value problem

y 0 = x − y 2, y (0) = 1.

Starting solutions required are to be obtained using Runge-Kutta

method of order 4 using step value h = 0.1.
26 Using Milne’s method, find y (4.5) given

5xy 0 + y 2 − 2 = 0

with y (4) = 1, y (4.1) = 1.0049, y (4.2) = 1.0097, y (4.3) = 1.0143,

y (4.4) = 1.0187.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 48/55
Exercises

27 Given
y 0 = x(x 2 + y 2 )e −x , y (0) = 1,
find y at x = 0.1, 0.2, and 0.3 by Taylor’s series method and compute
y (0.4) by Milne’s method.
28 Using Runge-Kutta method of order 4, find y for x = 0.1, 0.2, 0.3
given that
dy
= xy + y 2 , y (0) = 1.
dx
Continue the solution at x = 0.4 using Milne’s method.
29 If
dy
= 2e x y , y (0) = 2,
dx
find y (4) using Adams predictor-corrector formula by calculating
y (1), y (2) and y (3) using Euler’s modified formula.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 49/55
Exercises

30 Given
y 00 + xy 0 + y = 0, y (0) = 1, y 0 (0) = 0
obtain y for x = 0.1, 0.3 by any method. Further, continue the
solution by Milne’s method to calculate y (0.4).
31 Given
dy
= 1 + y2
dx
where y = 0 when x = 0. Find y (0.8) by Adams-Bashforth formula.
Find y (0.2), y (0.4), y (0.6) by fourth order Runge-Kutta-method.
32 Given
dy
= 1 + y2
dx
where y = 0 when x = 0. Find y (0.8) and y (1.0) by Milne’s formula.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 50/55
Exercises

33 The differential equation y 0 = x 2 + y 2 − 2 satisfies the following data:

x y
-0.1 1.0900
0 1.0000
0.1 0.8900
0.2 0.7605

Use Milne’s method to obtain the value of y (0.3).

34 Using Adams-Bashforth predictor-corrector formulae, evaluate y (1.4)
if y satisfies
dy y 1
+ = 2
dx x x
and y (1) = 1, y (1.1) = 0.996, y (1.2) = 0.986, y (1.3) = 0.972.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 51/55
Exercises

35 Find y (2) by Milne’s method if y (x) is the solution of

dy 1
= (x + y )
dx 2
assuming y (0) = 2, y (0.5) = 2.636, y (1) = 3.595, y (1.5) = 4.968.
36 Tabulate by Milne’s method the numerical solution of
dy
=x +y
dx
from x = 0.2 to x = 0.3 given that

x y
0 1
0.05 1.0526
0.1 1.1104
0.15 1.1737
P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 52/55
Exercises

37 Solve the differential equation

dy
= x2 + y2 − 2
dx
given y (−0.1) = 1.09, y (0) = 1, y (0.1) = 0.89. Find y (0.2) by series
expansion and find y (0.3) by Milne’s method.
38 Solve the initial value problem
dy
= 1 + xy 2 , y (0) = 1
dx
for x = 0.4, 0.5 by using Milne’s method when it is given that

x 0.1 0.2 0.3

y 1.105 1.223 1.355

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 53/55
Exercises

39 Using the Adams method, solve the simultaneous differential

equations
dy
= xy + z
dx
dz
= y −z
dx
with y (0) = 0, z(0) = 1.
40 Use Milne’s method to solve the simultaneous differential equations
dy
= x +z
dx
dz
= −xy
dx
with y (0) = 1, z(0) = 0.
P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 54/55
References

1 Richard L. Burden and J. Douglas Faires, Numerical Analysis -

Theory and Applications, Cengage Learning, Singapore.
2 Kendall E. Atkinson, An Introduction to Numerical Analysis, Wiley
India.
3 David Kincaid and Ward Cheney, Numerical Analysis -
Mathematics of Scientific Computing, American Mathematical
Society, Providence, Rhode Island.
4 S.S. Sastry, Introductory Methods of Numerical Analysis, Fourth
Edition, Prentice-Hall, India.

P. Sam Johnson (NITK) Numerical Solution of Ordinary Differential Equations (Part - 2) May 3, 2020 55/55