1. Module Title: MATHEMATICS FOR ENGINEERS IV .

2. Level: 2 Semester: 2 Credits: 10.

3. Pre-requisite or co-requisite modules, excluded combinations. Pre-requisite: Mathematics

for Engineers I, II and III

4. Brief descriotion of aims and content

The Module aims to introduce students to the various techniques of numerical analysis applied
to transcendental equations, mathematical functions, matrix algebra and differential equations.
The module will deal with these topics at a basic level, leaving more advanced techniques to
the specialist courses in Engineering.

LEARNING OUTCOMES
[Link] and Understanding

Having successfully completed the module, students should be able to demonstrate knowledge
and understanding of:
1.1 Numerical analysis as applied to simple mathematical topics.

1.2 The implications of basic numerical analysis.

2. Cognitive/Intellectual skills/Apllications of Knowledge

2.1 Present simple arguments and conclusions using numerical analysis.

2.2 Analyse and evaluate problems in mathematics and engineering.

3. Communication /ICT/Numeracy/Analytic Techniques/Practical Skills.

Having successfully completed the module, students should be able to:
3.1. Apply the mathematical knowledge to solve problems in range of Engineering situations.
3.2. Carry out mathematical and numerical manipulation with confidence and accuracy.

4. General transferable skills

Having successfully completed the module, students should be able to:
4.1. Assimilate Abstract Ideas;
4.2. Communicate information having a mathematical content.

1
 INDICATIVE CONTENT

Chapter I: Transcendental, Polynomial Equations and Curve Fitting .

I.1. Bisection and Newton-Raphson method.

I.2. Method of least squares - Fitting a straight line - Fitting a Parabola.
I.3. Solution of a system of linear ordinary differential equations.
I.4. Iterative methods for linear equations.

Chapter II. Approximation of Functions.

[Link] and higher differences: - forward differences and backward differences and central
differences - differences of a polynomial.
II.2. Interpolation: - Newton Gregory Forward and Backward Interpolation formulae - Inverse
Interpolation.
II.3. Numerical Differentiation and integration: Newton forward and backward differences
formula to compute first and higher order derivatives - The Trapezoidal Rule - Simpson’s one
third rule and three eighth rule.

Chapter III. Matrix Algebra

III.1. Solution of simultaneous linear equations.

III.2. Iterative methods for linear equations.

Chapter IV. Numerical Solutions of Differential Equations.

IV.1. Solution of ode’s by Taylor series, Euler’s method, Runge - Kutta methods.
IV.2. Introduction to solution of pde’s: Classification of Partial differential equations of the
second order.
IV.3. Laplace’s equation and its solution.


Quizzes: 10

Assignments: 10
ASSESSMENT STRATEGY:
CATs : 30



Final Exam: 50

REFERENCES:

1. Advanced Engineering Mathematics by E. Kreysig.

2
1 CHAPTER I: Transcendental, Polynomial Equations
and Curve Fitting.
1.1 Introduction: Numerical solutions of non-linear equations
A frequently occuring problem is to find one or more roots of the equations

f (x) = 0. (1)

If f (X) is a quadratic, cubic, or a biquadratic expression, the roots may be evaluated by using
the available algebraic formula. On the other hand, when f (x) is a polynomial of higher degree
or an expression involving transcendental functions,
for example: 1 + cos(x) − 5x, x tan(x) − cosh(x), etc.
Algebraic methods are not available, and rescourse must e taken to find the roots by approx-
imate method. In this section, we describe some numerical methods for the solution of the
equation 1, where f (x) may be algebraic, transcendental or combination of both.

1.2 The Bisection Method

This method is based on an important property of continuous functions, that is:

If a function f (x) is continuous in the closed interval (a, b) and if f (a) and f (b) are of opposite
signs, then there exists at least one root between a and b.
For definition, let f (a) be negative and f (b) be positive. Then the root lies between a and b
and let its approximate value be given by
a+b
x1 = .
2
Now if f (x1 ) = 0, we conclude that x1 is a root of the equation f (x) = 0. Otherwise, the root
lies either between x1 and b, or between a and x1 depending on whether f (x1 ) is negative or
positive. Them we bisect the interval and repeat the process until the root is known to the
desired occuracy.
If f (x1 ) is positive the root lies in between a and x1 . The second approximation to the root is
a + x1
x2 = .
2
If f (x2 ) is negative, the root lies in between x2 and x1 , and the third approximation to the root
is
x1 + x2
x3 = ,
2
and so on. This method is simple but slowly convergent.

Example 1. Find a real root of the equation x3 − x − 1 = 0.

Solution: Let f (x) = x3 − x − 1. Since f (1) = −1 < 0 and f (2) = 5 > 0, a root lies in
between 1 and 2. Therefore, we taken
1+2 3
x1 = = .
2 2

3
 3 3
− 32 − 1 = 7 3


Then f (x1 ) = 2 2
> 0. Hence, the root lies between 1 and 2
= 1.5 and we obtain
1 + 1.5 5
x2 = = 1.25 = .
2 4
19
Now, f (x2 ) = − 64 < 0. Hence the root lies between x1 and x2 , and we obtain
3 5
2
+ 4
x3 = = 1.375.
2
f (x3 ) = 0.224609375 > 0
, hence the root lies between x2 and x3 , and we obtain
x2 + x3 1.25 + 1.375
x3 = = = 1.3125.
2 2
Repeating the process yields the successive approximations until you get the same number
repeating at least four numbers after a comman.
EXERCISES 1:
I.1 Find the root of the equation x3 − x − 11 = 0 correct to the four decimal places using
bisection method.

I.2 Obtain a root of the following equations correct to the three decimal places using bisection
method
(1) x3 − 18 = 0 (2) x3 − x − 4 = 0 (3) x3 − 4x − 9 = 0 (4) x3 + x2 + x + 7 = 0

1.3 The Regular-Falsi Method or Method of False Position

The method of false position is also known as regular falsi Method, is the oldest method of
finding the real root of an equation f (x) = 0 and is somewhat similar to the bisection method.
In this method, we choose two points a and b ( a < b) such that f (a) and f (b) are of opposite
signs. Since the graph of y = f (x) crosses the x-axis between these two points, a root must lies
between these points chosen.
The equation of the chord joining the points (a, f (a)) and(b, f (b)) is given as
y − f (a) f (b) − f (a)
=
x−a b−a
Therefore,
f (b) − f (a)
y − f (a) = (x − a). (2)
b−a
The method consists in replacing the part of the curve bwtween (a, f (a)) and (b, f (b)) by
means of the chord joining these points and taking the intersection of the chord with the x-axis
as an approximation to the root.
In the present case, the point of the intersection is obtained by setting y = 0 in equation 2.
This
f (a)
x=a− (b − a).
f (b) − f (a)
Hence, the first approximation to the root is given by
af (b) − bf (a)
x1 = . (3)
f (b) − f (a)

4
Now if f (x1 ) and f (a) are of opposite signs, then the root lies between a and x1 . So we replace
b by x1 in the equation (3) and get the next approximation x2 given by

f (a) af (x1 ) − x1 f (a)

x2 = a − (x1 − a) = . (∗)
f (x1 ) − f (a) f (x1 ) − f (a)

But if f (x1 ) and f (a) are of the same signs, then f (x1 ) and f (b) will be of the opposite signs
and therefore, the root lies in between x1 and b. Hence, we replace a by x1 in the equation (3)
and get the next approximation x2 .

x1 f (b) − bf (x1 )
x2 = . (∗∗)
f (b) − f (x1 )

The process is to be repeated till the root is found to the desired occuracy.

Example 2. Find a real root of x3 − 2x − 5 = 0 by the method of false position correct to three
decimal places.

Solution: It is easily to see that f (2) = −1 < 0 and f (3) = 16 > 0. Therefore, a root lies
betweeen 2 and 3. It means that a = 2 and b = 3.
The first approximation of the root is x1 and is given by

af (b) − bf (a) 2(16) − 3(−1)

x1 = = = 2.0588.
f (b) − f (a) 16 − (−1)

Now, f (x1 ) = −0.391 < 0. Therefore, the root lies in between x1 = 2.0588 and b = 3. The
second approximation to the root is given by

x1 f (b) − bf (x1 )
x2 = = 2.08125.
f (b) − f (x1 )

Then f (x2 ) = −0.147 < 0, which means that the root of f (x) = 0 lies in between x2 = 2.08125
and b = 3. The third approximation to the root is given by

x2 f (b) − bf (x2 )
x3 = = 2.0896.
f (b) − f (x2 )

Then f (x3 ) = −0.0551 < 0. Hence, the root of f (x) = 0 lies in between x3 = 2.0896 and b = 3.
The fouth approximation to the root is given by

x3 f (b) − bf (x3 )
x4 = = 2.0927.
f (b) − f (x3 )

Now, f (x4 ) = f (2.0927) = −0.0206 < 0. Which implies that the root lies in between x4 =
2.0927 and b = 3. The fifth approximation to the root is given by

x4 f (b) − bf (x4 )
x5 = = 2.0939.
f (b) − f (x4 )

Now f (x5 ) = f (2.0939) = −0.00726 < 0. Therefore, the root lies in between x5 = 2.0939 and
b = 3. The sixth approximation to the root is given by

x5 f (b) − bf (x5 )
x6 = = 2.0943.
f (b) − f (x5 )

5
Now f (x6 ) = f (2.0943) = −0.00286 < 0. Therefore, the root lies in between x6 = 2.0943 and
b = 3. The seventh approximation to the root is given by
x6 f (b) − bf (x6 )
x7 = = 2.0944.
f (b) − f (x6 )
Therefore, the root is 2.094, correct to the three decimal places.

EXERCISES 2:

1. Find the real root of xex = 3 by regular-falsi method correct to three decimal places.

2. Find the root correct to the five places of decimal of the equation 2x − log(x) = 7 lies
between 3.5 and 4.0, using Regula-falsi method.

3. Show that x3 − 2x − 17 = 0 has only one real root and find its value correct to two
decimal places.

4. Find a root of the following equations correct to four decimal places, using the method
of false position.

a) x3 − 4x − 9 = 0 b) x3 + 2x2 + 10x − 20 = 0 c) ex sin(x) = 1 d) x = cos(x)

1.4 Newton-Raphson Method

This method is also known as Newton’s Iteration method.
The solution of f (x) = 0 may often be found by a simple procedure known as the Newton-
Raphson method. This method consists in replacing the part of the curve between a point,
say, (x0 , f (x0 )) and the x-axis by means of the tangent to the curve at the point.
The x-intercept of the tangent is then used as the first approximation. We have
f (x0 )
f 0 (x0 ) = tanθ = (∗)
x0 − x1
and this gives
f (x0 )
x1 = x 0 − . (∗∗)
f 0 (x0 )
Now substituting x1 for x0 and x2 for x1 the next better approximations are given by
f (x1 ) f (x2 )
x 2 = x1 − and x3 = x2 − .
f 0 (x1 ) f 0 f (x2 )
Proceeding in the same way n times, we get the general formulae given as
f (xn ) xn f 0 (xn ) − f (xn )
xn+1 = xn − = . (4)
f 0 (xn ) f 0 (xn )
For n = 0, 1, 2, 3, . . . The formulae (4) can be repeatedly used to find improved approximation
to the real root. But it is obvious from formulae (4) that this method fails if the slope of
the curve becomes zero.
The recurrence relation (4) may also be derived from the Taylor’s series expension of f (x) about
x0 .

6
Let x0 ne am approximate root of f (x) = 0 and let x1 = x0 + h be the correct root so that
f (x1 ) = 0. i.e, f (x0 + h) = 0.
Expanding f (x0 + h) by Taylor’s series, we obtain

h2 00
f (x0 + h) = f (x0 ) + hf 0 (x0 ) + f (x0 ) + . . . = 0.
2!
Neglating the second and higher order derivatives, we have

f (x0 ) + hf 0 (x0 ) = 0.

It means that
f (x0 )
h=− .
f 0 (x0 )
A better approximation than x0 is therefore given by x1 where
f (x0 )
x1 = x 0 − .
f 0 (x0 )
Successive approximations are given by x2 , x3 , x4 , . . . , xn+1 where
f (xn )
xn+1 = xn − . (5)
f 0 (xn )

Which is the same as the formulae (4).

Note: If the initial approximation x0 to the root is not given, choose two values of x, say a
and b, such that f (a) and f (b) are of opposite signs. Therefore, if |f (a)| < |f (b)|, then take
x0 = a.

Geometric Interpretation Let the curve f (x) = 0 meet the x-axs at x = α. It means that α
is the original root of f (x) = 0. Let x0 be the point near the root α of the equation f (x) = 0.
Then, the equation of the tangent at P0 (x0 , f (x0 )) is given by y − f (x0 ) = f 0 (x0 )(x − x0 ) which
cuts the x-axis at
f (x0 )
x1 = x0 − 0 .
f (x0 )
This is the first approximation to the root α.
If P1 (x1 , f (x1 )) is the point corresponding to x1 on the curve, then the tangent at P1 is given
by y − f (x1 ) = f 0 (x1 )(x − x1 ) which cuts the x-axis at

f (x1 )
x2 = x1 − .
f 0 (x1 )
This is the second approximation to the root α.
Repeating this process, we approach to the root α with better approximations quite rapidly.

Example 3. Solve x3 + 2x2 + 10x − 20 = 0 by Newton-Raphson method.

Solution: Let f (x) = x3 + 2x2 + 10x − 20 , then f 0 (x) = 3x2 + 4x + 10. Thus

f (xn ) x3n + 2x2n + 10xn − 20 2 [x3n + x2n + 10]

xn+1 = xn − = X n − = .
f 0 (xn ) 3x2n + 4xn + 10 3x2n + 4xn + 10

7
Now, we can see that f (1) = −7 < 0 and f (2) = 16 > 0. Since |f (1)| < |f (2)|, then x0 = 1.

x0 = 1
2[x30 + x20 + 10]
x1 = = 1.4118
3x20 + 4x0 + 10
2[x31 + x21 + 10]
x2 = = 1.3693
3x21 + 4x1 + 10
2[x32 + x22 + 10]
x3 = = 1.3688
3x22 + 4x2 + 10
2[x33 + x23 + 10]
x4 = 1.3688
3x23 + 4x3 + 10

Hence, the root is 1.3688.

EXERCSISE 3:

1. Using Newton-Raphson method, compute approximately the square root of 24.

√ √
2. Find an iterative formulae for N , where N is a positive number and hence find 12
correct to four decimal places.

3. Solve sin(x) = 1 + x2 using Newton-Raphson method.

1.5 Method of Successive Approximation

This method is also known as Iteration method. Let f (x) = 0 be the given equation whose
roots to be determined. This equation can be written in the form

x = φ(x). (6)

Let x0 be an initial approximation to the actual root, say, α of (6). Then the first root approxi-
mation is x1 = φ(x0 ) and the successive approximation are x2 = φ(x1 ), x3 = φ(x2 ), . . . , xn =
φ(xn−1 ).
If the sequence of approximate roots, x1 , x2 , x3 , . . . , xn converges to α, it is taken as the root
of the equation f (x) = 0.
For convergence purpose, the initial approximation x0 is to be done carefully. The choice of x0
is determined according to the following theorem:

Theorem 1.1. If α is a root of f (x) = 0 which is equivalent to x = φ(x), I is any interval

containing the point x = α and |φ0 (x)| < 1 for all x ∈ I, then the sequence of approximations
x0 , x1 , x2 , . . . , xn will converge to the root α provided the initial approximation x0 is
chosen in interval I.

Note: Smaller the value of φ0 (x), more rapid will be the convergence.

Example 4. Find a real root of the equation x3 + x2 − 1 = 0 by iteration method correct to

four decimal places.

Solution: Let f (x) = x3 + x2 − 1. Since f (0) = −1 < 0 and f (1) = 1 > 0, then a real root
lies in between 0 and 1.

8
Now, x3 + x2 − 1 = 0 can be written as
1
x2 (x + 1) = 1 ⇒ x = √
x+1
so
1 1 1
φ(x) = √ and φ0 (x) = − p .
x+1 2 (x + 1)3
Clearly, |φ0 (0)| = 21 < 1 and |φ0 (1)| = 1 √1
2 8
< 1. That is |φ0 (x)| < 1 for all x ∈ [0, 1]. Hence, the
iterative method can be applied.

Let x0 = 0.5 be an initial approximation, then

1
x1 = φ(x0 ) = √ = 0.8165
0.5 + 1
1
x2 = φ(x1 ) = √ = 0.7420
0.8165 + 1
1
x3 = φ(x2 ) = √ = 0.7577
0.7420 + 1
1
x4 = φ(x3 ) = √ = 0.7543
0.7577 + 1
1
x5 = φ(x4 ) = √ = 0.7550
0.7543 + 1
1
x6 = φ(x5 ) = √ = 0.7549
0.7550 + 1
1
x7 = φ(x6 ) = √ = 0.7549
0.7549 + 1

Hence, the root is 0.7549.

Try x0 = 0.65 and find the root correct to four decimal places.

EXERCISE 4:

1. Find a real root of the equation cos(x) = 3x − 1 correct to four decimal places by us-
ing method of successive approximation.

2. Find a root of the following equations correct to three decimal places, using iteration method
1
a) x3 + x2 − 100 = 0, b) x = + sin(x), c) 3x − 6 = log(x), d) xex − cos(x) = 0.
2

1.6 Method of Least Squares: Fitting a straight line and Fitting a

parabola.
In this section, we will fit the following curves:

• A straight line: y = ax + b;

• A parabola: y = ax2 + bx + c.

9
1.6.1 Fitting a straight line
Let (xi , yi ), i = 1, 2, 3, . . . , n be n sets of observations of related data and

y = ax + b (7)

be the straight line to be fitted.

The residual at x = xi is given by

ei = yi − f (xi ) = yi − (axi + b), i = 1, 2, 3, 4, . . . , n.

Therefore,
n
X n
X
E= e2i = [yi − (axi + b)]2 .
i=1 i=1

By the principal of least square, E is minimum. Hence,

∂E ∂E
=0 and = 0.
∂a ∂b

n
∂E X
= −2(yi − axi − b)(xi ) = 0 or
∂a i=1
n
X
∴ −2 (yi − axi − b)(xi ) = 0 or
i=1
n
X
(yi − axi − b) = 0
i=1

Here we have n n n
X X X
yi xi = a x2i +b xi , (8)
i=1 i=1 i=1

and n n
X X
yi = a xi + nb. (9)
i=1 i=1

Sincexi and yi are known in the equations (8) and (9)results in two equations is a and b.
Solving these the best values of a and b can be known and hence

y = ax + b.

Note: Depending on the sufiix i, we can write equations (8) and (9) as
X X X
yx = a x2 + b x and
X X
y = a x + nb.

These are known as Normal equations to the curve y = ax + b.

Example 5. By the method of least squares, find the straight line that best fits the following
data:
x 1 2 3 4 5
y 14 27 40 55 68

10
Solution: Let the straight line of the best fits be y = ax + b (i). The normal equations are
X X X
yx = a x2 + b x (∗) and
X X
y = a x + nb. (∗∗)

x2 ,
P P P P
Here n = 5 (sample size). The values of x, y, and xy are calculated as
shown bellow.
x y x2 xy
1 14 1 14
2 27 4 54
3 40 9 120
4 55 16 220
P 5 P 68 25
P 2 P 240
x = 15 y = 204 x = 55 xy = 748

The equations (*) and (**) becomes:

748 = 55a + 15b, and

204 = 15a + 5b.

Solving these, we have a = 13.6 and b = 0. Putting these values in equation (i), we get the
straight line of the best fits given as
y = 13.6x.
EXERCISE 5:

By the method of least squares, find the straight line that best fits the following data.

a)

x 2 4 6 8 10 12 14 16
y 18 23 20 24 30 28 35 40

b)

x 100 125 140 160 180 200

y 0.45 0.55 0.60 0.70 0.80 0.85

1.6.2 Fitting a Parabola

Let y = ax2 + bx + c be a parabola to be fitted for the data (xi , yi ), i = 1, 2, 3, . . . , n. The
residual at x = xi is
ei = yi − f (xi ) = yi − (ax2i + bxi + c).
Hence,
n n
X X
2
E= e2i = yi − (ax2i + bxi + c) .
i=1 i=1

This is called the sum of square errors. Now E should be minimum for the best values of a, b
and c [Principal of least squares].

∂E ∂E ∂E
= 0, = 0, and = 0.
∂a ∂b ∂c
11
It means that
n
X
2 yi − (ax2i + bxi + c) (−x2i = 0,


i=1
n
X
2 yi − (ax2i + bxi + c) (−xi ) = 0,


i=1
Xn
2 yi − (ax2i + bxi + c) (−1) = 0,


i=1

which results on simplifications after dropping the suffix as

X X X X
x2 y = a x4 + b x3 + c x2 ,
X X X X
xy = a x3 + b x2 + c x,
X X X
y = a x2 + b x + nc.

These are called the Normal equations of the parabola y = ax2 + bx + c.

Solving the normal equations, we get the best values of the constants a, b and c and hence, best
fit for equation
y = ax2 + bx + c.

Example 6. Fit a parabola y = ax2 + bx + c in least square sense to the data

x 10 12 15 23 20
y 14 17 23 25 21
Solution: The normal equations to the curve are
X X X
y = a x2 + b x + 5c
X X X X
xy = a x3 + b x2 + c x
X X X X
x2 y = a x4 + b x3 + c x2

x2 , x3 , x4 , x2 y,
P P P P P P P
Since n = 5 and the values of x, xy, and y
are calculated in the following table

x y x2 x3 x4 xy x2 y
10 14 100 1000 10000 140 1400
12 17 144 1728 20736 204 2448
15 23 225 3375 50625 345 5175
23 25 529 12167 279841 575 13225
P 20 P 21 P 2200 P 38000 P 160000
4
P 420 P 2 8400
x = 80 y = 100 x = 1398 x = 26270 x = 521202 xy = 1684 x y = 30648
Substituting the obtained values from the table in normal equations, we have

100 = 1398a + 80b + 5c

1684 = 26270a + 1398b + 80c

30648 = 521202a + 26270b + 1398c


12
On solving, a = −0.07, b = 3.03, c = −8.89. Therefore,

y = −0.07x2 + 3.03x − 8.89.

EXERCISE 6:

Fit a second degree parabola to the following data using method of least squares.

x 0 1 2 3 4
a)
y 1 1.8 1.3 2.5 6.3

x 1 2 3 4 5 6 7 8 9
b)
y 2 6 7 8 10 11 11 10 9

13
2 Chapter 2: Approximation of Functions
Assume that we have a set of values (xi , yi ), i = 1, 2, 3, . . . , n of functions y = f (x), the value
of x being equally spaced i.e

xi = x0 + ih, i = 1, 2, 3, . . . , n.

x1 = x0 + h, x2 = x0 + 2h, x3 = x0 + 3h, . . . , xn = x0 + nh.

Suppose that we are required to recover the values of f (x) for some intermediate values of x,
or to obtain the derivative of f (x) for some x in the range x0 ≤ x ≤ xn .
The methods for the solution of these problems are based on the concept of the differences of
a function which we now proceed to define.
When the values of x are equally spaced as above, interplation formulae are based on three
types of differences:

1. Forward differences;

2. Backward differences;

3. Central differences.

2.1 Forward Differences

If we subtract from each value of y (except y0 ) the proceding value of y, we get y1 − y0 , y2 −
y1 , y3 −y2 , . . . yn −yn−1 respectively known as the first differences of y. Or if y0 , y1 , y2 , . . . , yn
denotes a set of values of y, then y1 −y0 , y2 −y1 , y3 −y2 , . . . yn −yn−1 are called the differences
of y. Denoting these differences by

4y0 , 4y1 , 4y2 , 4y3 , . . . 4yn .

We have:

4y0 = y1 − y0 , 4y1 = y2 − y1 , 4y2 = y3 − y2 , 4y3 = y4 − y3 , . . . 4yn = yn − yn−1 .

Where 4 is called the forward difference operator, and 4y0 , 4y1 , 4y2 , 4y3 , . . . 4yn
are called the first forward differences. In general, the first forward differences are given by

4yi = yi+1 − yi , i = 0, 1, 2, 3, . . . , n.

The differences of the first forward differences are called the second forward differences and are
denotes by
42 y0 , 42 y1 , 42 y2 , 42 y3 , . . . 42 yn .
Now, the second forward differences are defines as the differences of the first differences, that
is:
42 y0 = 4(y1 − y0 ) = 4y1 − 4y0 = (y2 − y1 ) − (y1 − y0 ) = y2 − 2y1 + y0 .
42 y1 = 4(y2 − y1 ) = 4y2 − 4y1 = (y3 − y2 ) − (y2 − y2 ) = y3 − 2y2 + y1 .
In general
42 yn = 4yn+1 − 4yn = yn+2 − 2yn+1 + yn .

14
Here, 42 is called the second forward difference operator. Similarly, one can define third
forward differences are:

43 y0 = 42 (4y0 ) = 42 y1 − 42 y0 = (y3 − 2y2 + y1 ) − (y2 − 2y1 + y0 ) = y3 − 3y2 + 3y1 − y0 .

43 y1 = 42 (4y1 ) = 42 y2 − 42 y1 = y4 − 3y3 + 3y2 − y1 .

Generally, we have

43 y0 = 42 (4yn ) = 42 yn+1 − 42 yn = yn+3 − 3yn+2 + 3yn+1 − yn .

The fourth forward differences can be defined as

44 y 0 = 43 y 1 − 4 3 y 0
= (y4 − 3y3 + 3y2 − y1 ) − (y3 − 3y2 + 3y1 − y0 )
= y4 − 4y3 + 6y2 − 4y1 + y0

Hence,

44 yn = 43 yn+1 − 43 yn = yn+4 − 4yn+3 + 6yn+2 − 4yn+1 + yn .

It is therefore, clear that any higher order forward differences can easly be expressed in terms
of the ordinated, since the coefficients occuring on the right side are the binomial coefficients.
In general, the nth differences are difined as:

4n yk = 4n−1 yk+1 − 4n−1 yk .

In function notation, the forward differences are as written below:

4f (x) = f (x + h) − f (x),
42 f (x) = f (x + 2h) − 2f (x + h) + f (x),
43 f (x) = f (x + 3h) − 3f (x + 2h) + 3f (x + h) − f (x),

and so on, where h is the interval of differencing.

The following table shows how the forward differences of all orders can be formed:
x y = f (x) 4 42 43 44 45
x0 y0
x1 y1 4y0
x2 y2 4y1 42 y 0
x3 y3 4y2 42 y 1 43 y0
x4 y4 4y3 42 y 2 43 y1 44 y 0
x5 y5 4y4 42 y 3 43 y2 44 y 1 45 y 0

The first term in the table y0 is called the leading term and the differences 4y0 , 42 y0 , 43 y0 , . . .
are called the leading differences. The above differences table is known as forward difference
table or Diagonal difference table.

Properties: The operator 4 satisfies the following properties

• 4(f (x) ± g(x)) = 4f (x) ± 4g(x) i.e 4 is linear.

• 4(αf (x)) = α4f (x), α being a constant.

15
 • 4m 4n f (x) = 4m+n f (x), where m and n are positive integers.
• 4 [f (x).g(x)] 6= f (x).4g(x).
Observation 1: We can express any higher order forward difference of y0 in terms of the entire
y0 , y1 , y2 , . . . , yn . From
4y0 = y1 − yo ,
42 y0 = y2 − 2y1 + y0 ,
43 y 0 = y3 − 3y2 + 3y1 − y0 ,
44 y 0 = y4 − 4y3 + 6y2 − 4y1 + y0 , and so on.
We can see that the coefficients of the entries on the RHS are binomial coefficients. Therefore,
in general      
n n n n
4 yn = yn + yn−1 + yn−2 + yn−3 + . . . + (−1)n y0 (10)
1 2 3
Observation 2: We can express any value of y in terms of leading entry y0 . We know that
y1 − y0 = 4y0 ⇒ y1 = y0 + 4y0
⇒ y1 = (1 + 4)y0 .
Now,
y2 − y1 = 4y1 ⇒ y2 = y1 + 4y1
⇒ y2 = (1 + 4)y1
y2 = (1 + 4)2 y0 .
Similarly, y3 = (1 + 4)3 y0 and so on. In general,
     
n n n 2 n
yn = (1 + 4) y0 = y0 + 4y0 + 4 y0 + 43 y 0 + . . . + 4n y 0 . (11)
1 2 3

2.2 Backward Differences

The differences y1 − y0 , y2 − y1 , y3 − y2 , . . . , yn − yn−1 are called first backward differences if
they are denoted by ∇y1 , ∇y2 , ∇Y3 , . . . , ∇yn , respectively. So that
∇y1 = y1 − y0 ,
∇y2 = y2 − y1 ,
∇y3 = y3 − y2 , . . .
∇yn = yn − yn−1 ,

where ∇ is the backward difference operator.

Now, the second backward differences are defined as the differences of the first differences, i.e
∇2 y2 = ∇(∇y2 ) = ∇(y2 − y1 ) = ∇y2 − ∇y1 = (y2 − y1 ) − (y1 − y0 ) = y2 − 2y1 + y0 .
∇2 y3 = ∇(∇y3 ) = ∇(y3 − y2 ) = ∇y3 − ∇y2 = (y3 − y2 ) − (y2 − y1 ) = y3 − 2y2 + y1 .
In general,
∇n = ∇n−1 yk − ∇n−1 yk−1 . (12)
The backward difference table is given as

16
 x y = f (x) ∇ ∇2 ∇3 ∇4 ∇5
x0 y0
x1 y1 ∇y1
x2 y2 ∇y2 ∇2 y2
x3 y3 ∇y3 ∇2 y3 ∇3 y3
x4 y4 ∇y4 ∇2 y4 ∇3 y4 ∇4 y4
x5 y5 ∇y5 ∇2 y5 ∇3 y5 ∇4 y5 ∇5 y5
In function notation, these are written as

∇f (x) = f (x) − f (x − h),

∇f (x + h) = f (x + h) − f (x),
2
∇ f (x + 2h) = f (x + 2h) − f (x + h) + f (x),
∇3 f (x + 3h) = f (x + 3h) − 3f (x + 2h) + 3f (x + h) − f (x).

And so on, where h is the interval of differencing.

Properties:
• ∇(f (x) ± g(x)) = ∇f (x) ± ∇g(x), i.e ∇ is a linear operator.
• ∇(αf (x)) = α∇f (x), α being a constant.
• ∇m ∇n f (x) = ∇m+n f (x), m and n being positive integer.
• ∇(f (x).g(x)) 6= [∇f (x).g(x)].

2.3 Central Differences

The central difference operator δ is defined by the relations

y1 − y0 = δy 1 ,
2
y2 − y1 = δy 3 ,
2
y3 − y2 = δy 5 , . . .
2
yn − yn−1 = δyn− 1 .
2

For the higher order central differences, we have

δy 3 − δy 1 = δ 2 y1 ,
2 2

δy 5 − δy 3 = δ 2 y2 , . . .
2 2

δ 2 y2 − δ 2 y1 = δ 3 y 3 , and so on.
2

The central differences are tabulated below.

x y = f (x) δ δ2 δ3 δ4 δ5
x0 y0
x1 y1 δy 1
2
x2 y2 δy 3 δ 2 y1
2
x3 y3 δy 5 δ 2 y2 δ3y 3
2 2
x4 y4 δy 7 δ 2 y3 δ3y 5 δ 4 y2
2 2
x5 y5 δy 9 δ 2 y4 δ3y 7 δ 4 y3 δ5y 5
2 2 2

17
We can see from the table that central differences on the same horizontal line have the same
suffix. Also, all odd differences have a fractional suffix, and the even differences have an integer
suffix.

Note 1: From the three tables, we can see that only the notation changes, not the differences.
For example,
y1 − y0 = 4y0 = ∇y1 = δy 1 .
2

Note 2: If we write y = f (x) as y = fx or y = yx , then the entries corresponding to x, x +

h, x + 2h, . . . , are yx , yx+h , yx+2h , . . . , respectively, and

4yx = yx+h − yx , 42 yx = 4yx+h − 4yx , and so on.

Similarly,

∇yx = yx − yx−h ,
δyx = yx+ h − yx− h , and so on.
2 2

Example 7. Construct a forward difference table from the following data

x 0 1 2 3 4
y 1 1.5 2.2 3.1 4.6

Evaluate 43 y1 , yx , and y5 .
Solution: The forward difference table is given below.
x y 4 42 43 44
x0 = 0 y0 = 1
x1 = 1 y1 = 1.5 4y0 = 0.5
x2 = 2 y2 = 2.2 4y1 = 0.7 42 y0 = 0.2
x3 = 3 y3 = 3.1 4y2 = 0.9 42 y1 = 0.2 43 y0 = 0
x4 = 4 y4 = 4.6 4y3 = 1.5 42 y2 = 0.6 43 y1 = 0.4 44 y0 = 0.4

Now, 43 y1 = y4 − 3y3 + 3y2 − y1 = 42 y2 − 42 y1 = 4.6 − 3(3.2) + 3(2.2) − 1.5 = 0.4

We known that
       
x x 2 x 3 x
yx = y0 + 4y0 + 4 y0 + 4 y0 + 44 y 0 ,
1 2 3 4
x! x! x! x!
yx = y0 + 4y0 + 42 y 0 + 43 y 0 + 44 y 0 ,
(x − 1)!1! 2!(x − 1)! 3!(x − 3)! 4!(x − 4)!
1 1 1
yx = 1 + x(0.5) + x(x − 1)(0.2) + x(x − 1)(x − 2)(0) + x(x − 1)(x − 2)(x − 3)(0.4)
2 3! 4!
x 1 2 1 4
yx = 1 + + (x − x) + (x − 6x3 + 11x2 − 6x)
2 10 60
1  4
x − 6x3 + 17x2 + 18x + 60

∴ yx =
60
Therefore,
1  4
(5) − 6(5)3 + 17(5)2 + 18(5) + 60 = 7.5.

y5 =
60

18
Example 8. : (Forward differences)
Find the polynomial of degree three which has the values equal to 1, 15, 85 and 259 corresponding
to the values 0, 2, 4 and 6 of the argument.

Solution: Let the variables be x and y. Then we have the following forward difference
table.
x y 4 42 43
0 1
2 15 14
4 85 70 56
6 259 174 104 48
x−x0
Here, h = 2. To use the formulae, we take x0 = 0. Then x = x0 + ph i.e (p = h
) gives
p = x2 , 4y0 = 14, 42 y0 = 56, 43 y0 = 48.

x x x 1 x x x 1
y(x) = 1 + (14) + ( − 1) (56) + ( − 1)( − 2) (48) = x3 + x2 + x + 1
2 2 2 2 2 2 2 6

EXERCISE II.1:
1. Prove the following results:

• 4∇ = ∇4 = 4 − ∇δ 2
4 ∇
• 4+∇= ∇
− 4

2. Tabulate the forward differences for the given data

x 1 2 3 4 5 6 7 8 9
y 1 8 27 64 125 216 343 512 729

3. Form a backward differences table of the function f (x) = x3 − 3x2 − 5x − 7 for x =

−1, 0, 1, 2, 3, 4, 5.

4. Show that

• y3 = y2 + 4y1 + 42 y0 + 43 y0 .

• 42 y8 = y8 − 2y7 + y6 .

5. If y0 = 3, y1 = 12, y2 = 81, y3 = 2000, y4 = 100 show that 44 y0 = −7459.

19
2.4 INTERPOLATION WITH EQUAL INTERVALS
Interpolation is a technique of obtaining the value of a function for any intermediate values of
the independent variable i.e argument within an interval, when the values of the arguments are
given. Here, x(argument): x0 , x1 = x0 + h, x2 = x0 + 2h, . . . , xn and yx (entry): y0 , y1 , y2 , . . . , yn .

The process of finding the value of y corresponding to any value of x = xi between x0 and xn
is called interpolation. The process of finding the value of a function outside the given range
of arguments is called extrapolation. However, the term interpolation is applied to both
process.

2.4.1 Gregory-Newton Forward Interpolation Formula

It also known as Newton’s Forward Interpolation Formula. Let y = f (x) be a function
which takes the values y0 , y1 , y2 , . . . , yn for (n+1) values of x0 , x1 , x2 , . . . , xn , of the independent
variable x (argument).
Let these values of x be equidistant, i.e xi = x0 + ih, i = 0, 1, 2, 3, . . . , n and let y(x) be the
polynomial in x of nth degree, such that yi = f (xi ), i = 0, 1, 2, 3, . . . , n. Suppose we need to
evaluate y(x) near the beginning of table of values. Therefore,

y(x) = A0 +A1 (x−x0 )+A2 (x−x0 )(x−x1 )+A3 (x−x0 )(x−X1 )(x−x2 )+. . .+An (x−x0 )(x−x1 ) . . . (x−xn ).
(13)
Putting x = x0 , x1 , x2 , . . . , xn successively in Equation (13), we get y0 = A0 , y1 = A0 +
A1 (x1 − x0 ), y2 = A0 + A1 (x2 − x0 ) + A2 (x2 − x1 )(x2 − x1 ) and so on. From these

A0 = y0
y1 − A0 y1 − y0 4y0
A1 = = =
x1 − x0 x1 − x0 h
y2 − A0 − A1 (x2 − x0
A2 =
(x2 − x0 )(x2 − x1
y2 − y0 − A1 (2h)
=
2h2
y2 − y0 − 4y h
0
(2h)
= 2
2h
y2 − y0 − 24y0
=
2h2
y2 − y0 − 2(y1 − y0
=
2h2
y2 − 2y1 + y0
=
2h2
y2 − 2y1 + y0 1
∴ A2 = = 42 y0
2h2 2!h2
Similarly,
1
A3 =
3
43 y0 , and so on.
3!h
Putting these values in Equation (13), we get

4y0 42 y0 43 y0
y(x) = y0 + (x − x0 ) + (x − x 0 )(x − x 1 ) + (x − x0 )(x − x1 )(x − x2 ) + . . . (14)
h 2!h2 3!h3

20
 x−x0
Putting p = h
, i.e x = x0 + ph. where p is a real number, Equation (14) takes the form

p(p − 1) 2 p(p − 1)(p − 2) 3 p(p − 1)(p − 2) . . . [P − (n − 1)] n

yp = y0 +p4y0 + 4 y0 + 4 y0 +. . .+ 4 y0 .
2! 3! n!
(15)
Where yp = y(x0 + ph). Equation (15) is known as Newton’s Forward Interpolation
formula. The first two terms of (15) will be lineae interpolation and the first three terms will
be parabolic interpolation and so on.

2.4.2 Gregory-Newton Backward Interpolation Formula

It is also known as Newton’s Backward Interpolation Formula.
Let y = f (x) be a function which takes the values y0 , y1 , y2 , . . . , yn for (n + 1) values of
x0 , x1 , x2 , . . . , xn , of the independent variable x.
Let these values of x be equidistant, i.e xi = x0 + ih, i = 0, 1, 2, 3, . . . , n and let y(x) be the
polynomial in x of nth degree, such that yi = f (xi ), i = 0, 1, 2, 3, . . . , n. Suppose that it is
required to evaluate y(x) near the end of the table of values, then we can assume that

y(x) = A0 + A1 (x − xn ) + A2 (x − xn )(x − xn−1 ) + . . . + An (x − xn )(x − xn−1 ) . . . (x − x1 ). (16)

Putting x = xn , xn−1 , xn−2 , . . . , x1 successively in (16), we get

yn = y(xn ) = A0 ,
yn−1 = y(xn−1 ) = A0 + A1 (xn−1 − xn )
yn−2 = A0 + A1 (xn−2 − xn ) + A2 (xn−2 − xn )(cn−2 − xn−1 ),

and so on. These equations give

A 0 = yn ,
yn−1 − A0 yn−1 − yn yn − yn−1 ∇yn
A1 = = = = ,
xn−1 − xn xn−1 − xn xn − xn−1 h
yn−2 − yn − A1 (xn−2 − xn )
A2 =
(xn−2 − xn )(xn−2 − xn−1 )
yn−2 − yn − A1 (−2h)
=
2h2
yn−2 − yn − ∇yh n (−2h)
=
2h2
yn − 2yn−1 + yn−2
=
2h2
1
∴ A2 = ∇2 yn .
2!h2
Similarly,
1
A3 = ∇3 yn , and so on.
3!h3
Putting these values in Equation (16), we have
1 1 1
y(x) = yn + ∇yn (x−xn )+ 2 ∇2 yn (x−xn )(x−xn−1 )+ 3 ∇3 yn (x−xn )(x−xn−1 )(x−xn−2 )+. . .
h 2!h 3!h
(17)

21
 x−xn
Let p = h
, i.e x = xn + ph, where p is a real number. Then Equation (17) takes the form

p(p + 1) 2 p(p + 1)(p + 2) 3 p(p + 1)(p + 2) . . . (p + n − 1) n

yp = yn +p∇yn + ∇ yn + ∇ yn +. . .+ ∇ yn ,
2! 3! n!
(18)
where yp = y(xn + ph), Equation (18) is known as Newton’s Backward Interpolation
Formula.

Note: Since the formula involves the backward differences, it is called backward interpolation
formula sn is used to interpolate the values of y near the end of a set of tabular values.

Example 9. The following data give I, the indicated HP and V, the speed in knots developed
by ship.
V 8 10 12 14 16
I 1000 1900 3250 5400 8950

Find I when V = 9, using Newton’s forward interpolation formula.

Solution: We note that v = 9 is near to the beginning of the table. Hence, to get the corre-
sponding I, we use Newton’s forward interpolation formula.

The forward differences are calculated and tabulated as follows:

V I 4 42 43 44
8 1000
10 1900 900
12 3250 1350 450
14 5400 2150 800 350
16 8950 3550 1400 600 250

Here, V0 = 8, I0 = 1000, 4I0 = 900, 42 I0 = 450, 43 I0 = 350, 44 I0 = 250. Hence, the

interpolation polynomial will be of degree 4. That is

p(p − 1) 2 p(p − 1)(p − 2) 3 p(p − 1)(p − 2)(p − 3) 4

I = I0 + p4I0 + 4 I0 + 4 I0 + 4 I0
2! 3! 4!
Let Ip be the value of I when V = 9. Then

V − V0 9−8 1
p= = = = 0.5.
h 2 2
Therefore,

(0.5)(0.5 − 1) (0.5)(0.5 − 1)(0.5 − 2)

Ip = 1000 + (0.5)(900) + (450) + (350)
2! 3!
(0.5)(0.5 − 1)(0.5 − 2)(0.5 − 3)
+ (250) = 1405.8594
4!
Example 10. The amount A of a substance remaining in a reaction system after an interval
of time t in a certain chemical experiment is tabulated below:
t(min) 2 5 8 11
A(gm) 94.8 87.9 81.3 75.1

22
Obtain the value of A where t = 9 using Newton’s backward interpolation formula.
Solution: Since the value of t = 9 is near the end of the table, to get the corresponding value
of A we use Newton’s backward interpolation formula.

The backward differences are calculated abd tabulated below.

t A ∇ ∇2 ∇3
2 94.8
5 87.9 -6.9
8 81.3 -6.6 0.3
11 75.1 -6.2 0.4 0.1
Here,
tn = 11, An = 75.1, ∇An = −6.2, ∇2 An = 0.4, ∇3 An = 0.1
Hence, the interpolation polynomial will be of degree 3. That is
p(p + 1) 2 p(p + 1)(p + 2) 3
A = An + p∇An + ∇ An + ∇ An ,
2! 3!
let Ap be the value of A when t = 9. Then
t − tn 9 − 11 2
p= = =− .
h 3 3
Therefore,
        
−2 1 −2 −2 1 −2 −2 −2
Ap = 75.1 + (−6.2) + + 1 (0.4) + +1 + 2 (0.1)
3 2! 3 3 3! 3 3 3
Ap = 79.183951.

Find a polynomial which takes the following data

x 1 3 5 7 9 11
and hence compute yx at x = 2, 12.
y 3 14 19 21 23 28
Solution: The forward differences table is given by
x y 4 42 43 44
1 3
3 14 11
5 19 5 -6
7 21 2 -3 3
9 23 2 0 3 0
11 28 5 3 3 0
Taking
x − x0 x−1
x0 = 1, y0 = 3, p = = .
h 2
Using Newton’s forward interpolation formula, we get
1 1 1
yp = y0 + p4y0 + p(p − 1)42 y0 + p(p − 1)(p − 2)43 y0 + p(p − 1)(p − 2)(p − 3)44 y0
2!  3!  4!   
x−1 1 x−1 1 11 1 1
= 3+ (11) + (x − 1) − 1) (−6) + (x − 1) (x − 1) − 1 (x − 1) − 2 (3)
2 2! 2 2 3! 2 2 2
11 3 1
= 3 + (x − 1) − (x2 − 4x + 3) + (x3 − 9x2 + 23x − 15)
2 4 16
1  3
x − 21x2 + 159x − 91 .

∴ yp =
16
23
Again take
x − 11
xn = 11, yn = 28, p = .
2
Using Newton’s backward interpolation formula,
1 1
yp = yn + p∇yn + p(p + 1)∇2 yn + p(p + 1)(p + 2)∇3 yn
2! 3!
5 1 1 1 1
= 28 + (x − 11) + 2
(x − 11)(x − 9)(3) + (x − 11)(x − 9)(x − 7)(3)
2 2! 2 3! 23
5 1
= 28 + (x − 11) + (x − 11)(x − 9)(x − 1)
2 16
1
x3 − 21x2 + 159x − 91 .


∴ yp =
16
So we can use any one of the formula to find the polynomial. Therefore,
1
x3 − 21x2 + 159x − 91 .


yx =
16
Now,
1
(2)3 − 21(2)2 + 159(2) − 91 = 9.4375.


y2 =
16
1
(12)3 − 21(12)2 + 159(12) − 91 = 32.5625.


y12 =
16

2.5 EQUIDISTANT TERMS WITH ONE OR MORE MISSING

VALUES
When one or more of the values of the function y = f (x) corresponding to the equidistant
values of x are missing. We can find these missing values using finite difference operator E and
4. The method is best illustrated by the following example.
Example 11. Find the missing value in the following table.
x 16 18 20 22 24 26
y 43 89 - 155 268 388
Solution: Since five values are given it is possible to express y as a polynomial of fourth degree.
Hence, the fifth differences of y are zeros. Taking the origin for x at 16, from the given table
we have y0 = 43, y1 = 89, y3 = 155, y4 = 268, y5 = 388 and we have to find y2 . We know that
45 y0 = 0 for all values of x,
45 y0 = 0 i.e (E − 1)5 y0 = 0.
         
5 5 4 5 3 5 2 5
i.e E − E + E − E + E − 1 yo = 0
1 2 3 4
E 5 − 5E 4 + 10E 3 − 10E 3 + 5E − 1 y0 = 0


or
i.e E 5 y0 − 5E 4 y0 + 10E 3 y0 − 10E 2 y0 + 5Ey0 − y0 = 0.
Hence
y5 − 5y4 + 10y3 − 10y2 + 5y1 − y0 = 0
Substituting the given values,
388 − 5(268) + 10(155) − 10y2 + 5(89) − 43 = 0 ⇒ y2 = 100.

24
Example 12. Find the missing values in the following table of values of x and y.
x 0 1 2 3 4 5 6
y -4 -2 - - 220 546 1148

Hint: There being given five values and two missing values, we may have 45 y0 = 0 and 46 y0 =
0.

EXERCISE II. 2

1. From the following data find y at x = 43 using Newton’s forward interpolation formula.
x 40 50 60 70 80 90
y 184 204 226 250 176 304

2. The population of a certain town in decennial census was as given below. Estimate the
population for the year 1895.
Year (x) 1891 1901 1911 1921
Population in thousands (y) 46 66 81 101

3. The area A of a circle of diameter d is given for the following values

d 80 85 90 95 100
A 5026 5674 6362 7088 7854

Calculate the area of a circle of diameter d = 105.

4. From the following table, estimate the values of f (22) and f (42).

x 20 25 30 35 40 45
f (x) 354 332 291 260 231 204

5. Find the polynomial which takes the following data

x 4 6 8 10
Hence, calculate y at x = 5.
y 1 3 8 16

6. Obtain the estimate of the missing value in the following table

x 1 2 3 4 5
y 2 5 7 - 32

7. Given y0 = 3, y1 = 12, y2 = 81, y3 = 200, y4 = 100. Find 44 y0 without forming the

difference table.

25
2.6 CENTRAL DIFFERENCE INTERPOLATION FORMULA
The central difference formula are most suited for interpolation near the middle for a
tabulated set. The most important central difference formula are those due to Stirling,
Bassel and Everett.

For convenience, we state the central difference formula by taking the central ordinate as y0
corresponding to x = x0 :
x y 4 42 43 44
x−2 y−2
x−1 y−1 4y−2
x0 y0 4y−1 42 y−2
x1 y1 4y0 42 y−1 43 y−2
x2 y2 4y2 42 y0 43 y−1 44 y−2

2.6.1 Stirling’s Formula

4y−1 + 4y0 1 2 2 p(p1 − 12 ) (43 y−1 + 43 y−2 ) p2 (p2 − 12 ) 4
yp = y0 + p + p 4 y−1 + + 4 y−2
2 2! 3! 2 4!
p(p2 − 12 )(p2 − 22 ) (45 y−2 + 45 y−3 ) p2 (p2 − 12 )(p2 − 22 ) 6
+ + 4 y−3 + . . .
5! 2 6!

2.6.2 Bassel’s Formula

p(p − 1) 42 y−1 + 42 y0 (p − 12 )p(p − 1) 3
 
1 1
yp = (y0 − y1 ) + (p − )4y0 + + 4 y−1
2 2 2! 2 3!
(p + 1)p(p − 1)(p − 2) 44 y−2 + 44 y−1
 
(p + 1)p(p − 1)(p − 2) 3
+ 4 y−1 + + ...
3! 4! 2

2.6.3 Everett’s Formula

1 1
yp = qy0 + q(q 1 − 12 )42 y−1 + g(q 2 − 12 )(q 2 − 22 )44 y−2 + . . .
3! 5!
1 1
+py1 + p(p2 − 12 )42 y0 + p(p2 − 12 )(p2 − 22 )44 y−1 + . . .
3! 5!
where q = 1 − p and x = x0 + ph.

26
2.7 INTERPOLATION WITH UNEQUAL INTERVALS
2.7.1 Newton’s Divided Difference Formula
Let (xi , yi ), i = 0, 1, 2, 3, . . . , n be the set of values of x and y satisfying the function y = f (x),
the values of xi being not necessarily equally spaced.
f (x1 )−f (x0 )
The first divided difference of f (x) for the arguments x0 and x1 is defined as x1 −x0
and
denoted by [x0 , x1 ] so that
f (x1 ) − f (x0 )
f (x0 , x1 ) = [x0 , x1 ] = .
x1 − x0
Similarly, we have
f (x2 ) − f (x1 )
f (x1 , x2 ) = [x1 , x2 ] = .
x2 − x1
In general, we have
f (xn ) − f (xn−1 )
f (xn−1 , xn ) = [xn−1 , xn ] = .
xn − xn−1
Again, the second divided difference of f (x) for the arguments x0 , x1 and x2 is defined as
[x1 , x2 ] − [x0 , x1 ]
f (x0 , x1 , x2 ) = [x0 , x1 , x2 ] = .
x2 − x0
[x2 , x3 ] − [x1 , x2 ]
f (x1 , x2 , x3 ) = [x1 , x2 , x3 ] = .
x3 − x1
The third divided difference is given as
[x1 , x2 , x3 ] − [x0 , x1 , x2 ]
f (x0 , x1 , x2 , x3 ) = [x0 , x1 , x2 , x3 ] =
x3 − x0
In general, we have
[x1 , . . . xn−1 ] − [x0 , x1 , . . . , xn−2 ]
f (x0 , x1 , . . . xn−1 ) = [x0 , x1 , . . . xn−1 ] = . (19)
xn−1 − x0

Which is the nth order divided difference of f (x) for the arguments x0 , x1 , . . . , xn−1 .

The divided difference table

x Entry f (x) 1st divided difference 2nd divided difference 3rd divided difference
x0 f (x0 )
x1 f (x1 ) f (x0 , x1 )
x2 f (x2 ) f (x1 , x2 ) f (x0 , x1 , x2 )
x3 f (x3 ) f (x2 , x3 ) f (x1 , x2 , x3 ) f (x0 , x1 , x2 , x3 )
x4 f (x4 ) f (x3, x4 ) f (x2 , x3 , x4 ) f (x1 , x2 , x3 , x4 )

We note the following properties of divided differences:

•
f (x1 − f (x − 0) f (x0 ) − f (x1 )
[x0 , x1 ] = = = [x1 , x0 ].
x1 − x0 x0 − x 1
•
[x1 , x2 ] − [x0 , x1 ] [x2 , x1 ] − [x1 , x0 ] [x1 , x0 ] − [x2 , x1 ]
[x0 , x1 , x2 ] = = = ,
x 2 − x0 x2 − x0 x0 − x2

27
which shows that the divided difference is symmetric with respect to its arguments.

To derive a formula involving a function y = f (x) and its divided differences, we observe that
y − y0
[x, x0 ] = [x0 , x] =
x − x0
⇒ y = y0 + (x − x0 )[x, x0 ] (∗)

[x0 , x1 ] − [x, x0 ] [x, x0 ] − [x0 , x1 ]

[x, x0 , x1 ] = =
x1 − x x − x1
⇒ [x, x0 ] = [x0 , x1 ] + (x − x0 )[x, x0 , x1 ] (∗∗)

Substituting (**) into (*), we obtain

y = y0 + (x − x0 )[x0 , x1 ] + (x − x0 )(x − x1 )[x, x0 , x1 ] (∗ ∗ ∗)

Further, we have
[x, x0 , x1 ] − [x0 , x1 , x2 ]
[x, x0 , x1 , x2 ] =
x − x2
⇒ [x, x0 , x1 ] = [x0 , x1 , x2 ] + (x − x2 )[x, x0 , x1 , x2 ] (∗ ∗ ∗∗)

Putting (****) in (***), we obtain

y = y0 + (x − x0 )[x0 , x1 ] + (x − x0 )(x − x1 ) ([x0 , x1 , x2 ] + (x − x2 )[x, x0 , x1 , x2 ])

y = y0 + (x − x0 )[x0 , x1 ] + (x − x0 )(x − x1 )[x0 , x1 , x2 ] + (x − x0 )(x − x1 )(x − x2 )[x, x0 , x1 , x2 ].

Proceeding in this way, we obtain the formula

y = y0 + (x − x0 )[x0 , x1 ] + (x − x0 )(x − x1 )[x0 , x1 , x2 ] + (x − x0 )(x − x1 )(x − x2 )[x, x0 , x1 , x2 ]

+ . . . + (x − x0 )(x − x1 ) . . . (x − xn )[x, x0 , x1 , . . . , xn ]. (i)

Where the last term is the reminder term in the formula. The formula (i) above is called
Newton’s Divided interpolation formula.
Example 13. Construct a divided difference table for the following data.
x 4 5 7 10 11 13
f (x) 48 100 294 900 1210 2028
Solution: The divided difference table is givena as follows
x f (x) 1st divided diff 2nd divided diff 3rd divided diff 4th divided diff
4 48
100−48
5 100 5−4
= 52
294−100 97−52
7 294 7−5
= 97 7−4
= 15
900−294 202−97 21−15
10 900 10−7
= 202 10−5
= 21 10−4
=1
1210−900 310−202 27−21
11 1210 11−10
= 310 11−7
= 27 11−5
=1 0
2028−1210 409−310 33−27
13 2028 13−11
= 409 13−10
= 33 13−7
=1 0
Example 14. Using Newton’s divided difference formula, find the cubic function from the
following table of values:
x 0 2 3 5
y 1 15 40 156

28
Solution: We form the table of divided differences as follows:
x y 1st divided difference 2nd divided difference 3rd divided difference
0 1
2 15 7
3 40 25 6
5 156 58 11 1
Hence,

y = y0 + (x − x0 )[x0 , x1 ] + (x − x0 )(x − x1 )[x0 , x1 , x2 ] + (x − x0 )(x − x1 )(x − x2 )[x0 , x1 , x2 , x3 ]

y = 1 + (x − 0)(7) + (x − 0)(x − 2)(6) + (x − 0)(x − 2)(x − 3)(1)
∴ y = x3 + x2 + x + 1.

EXERCISE II. 3

1. Find a polynomial satisfied by the following data.

x -4 -1 0 2 5
f (x) 1245 33 5 9 1335

2. Find f (4) from the following table of values

x 0 1 3 5
f (x) 8 11 35 123

3. Find a third order polynomial satisfying the following date

x 0 2 3 6
f (x) 659 705 729 804

29
2.7.2 Lagrange’s Interpolation Formula
Let y = f (x) be a function which takes the (n + 1) values y0 , y1 , y2 , . . . , yn polynomial of nth
degree in x. Let this polynomial be the form

y = f (x) = a0 (x − x1 )(x − x2 ) . . . (x − xn ) + a1 (x − x0 )(x − x2 ) . . . (x − xn )

+a2 (x − x0 )(x − x1 )(x − x3 ) . . . (x − xn ) + . . . + an (x − x0 )(x − x1 ) . . . (x − xn ). (20)

Putting x = x0 , y = y0 in equation (20), we get

y0 = a0 (x0 − x1 )(x0 − x2 )(x0 − x3 ) . . . (x0 − xn )

y0
∴ a0 =
(x0 − x1 )(x0 − x2 )(x0 − x3 ) . . . (x0 − xn )

Again putting x = x1 and y = y1 in equation (20), we get

y1
a1 = .
(x1 − x0 )(x1 − x2 )(x1 − x3 ) . . . (x1 − xn )

Proceeding in the same way, we obatin

y2
a2 = .
(x2 − x0 )(x2 − x1 )(x2 − x3 ) . . . (x2 − xn )

Until
yn
an = .
(xn − x0 )(xn − x1 )(xn − x2 ) . . . (xn − xn−1 )

Substituting the values of a0 , a1 , a2 , . . . , an in equation (20), we get

(x − x1 )(x − x2 )(x − x3 ) . . . (x − xn )
y = f (x) = y0
(x0 − x1 )(x0 − x2 )(x0 − x3 ) . . . (x0 − xn )
(x − x0 )(x − x2 )(x − x3 ) . . . (x − xn )
+ y1
(x1 − x0 )(x1 − x2 )(x1 − x3 ) . . . (x1 − xn )
+.......................................
(x − x0 )(x − x1 )(x − x2 ) . . . (x − xn−1 )
+ yn . (21)
(xn − x0 )(xn − x1 )(xn − x2 ) . . . (xn − xn−1 )

This is knwon as Lagrange’s Interpolation formula.

Note:

1. This formula can be used irrespectively of whether the values x0 , x1 , x2 , . . . , xn are equally
spaced or not.

2. It is simple and easy to remember but its application is not speedly.

3. The main drawback of it is that if another interpolation value is inserted, then the inter-
polation coefficients are required to be recalculated.

Example 15. Use Lagrange’s interpolation formula to find the value of y when x = 10, if the
values of x and y are given as below:
x 5 6 9 11
y 12 13 14 16

30
Solution: Here,
x0 = 5, x1 = 6, x2 = 9, x3 = 11, y0 = 12, y1 = 13, y2 = 14, y3 = 16
. By Lagrange’s formula in equation (21)
(x − 6)(x − 9)(x − 11) (x − 5)(x − 9)(x − 11)
y= (12) + (13)
(5 − 6)(5 − 9)(5 − 11) (6 − 5)(6 − 9)(6 − 11)
(x − 5)(x − 6)(x − 11) (x − 5)(x − 6)(x − 9)
+ (14) + (16)
(9 − 5)(9 − 6)(9 − 11) (11 − 5)(11 − 6)(11 − 9)
Putting x = 10 in the above equation, we get
(10 − 6)(10 − 9)(10 − 11) (10 − 5)(x − 9)(10 − 11)
y= (12) + (13)
(5 − 6)(5 − 9)(5 − 11) (6 − 5)(6 − 9)(6 − 11)
(10 − 5)(10 − 6)(10 − 11) (10 − 5)(10 − 6)(10 − 9)
+ (14) + (16) = 14.6666667.
(9 − 5)(9 − 6)(9 − 11) (11 − 5)(11 − 6)(11 − 9)
EXERCISE II.3

1. Use Lagrange’s interpolation formula to find the form of f (x), given

x 0 2 3 6
y 648 704 729 792
2. By means of Lagrange’s formula, prove that
y1 = y3 − 0.3(y5 − y−3 ) + 0.2(y−3 − y−5 )

2.8 INVERSE INTERPOLATION

So far, given a set of values of x and y, we were required to find the value of y corresponding
to a value of x. Sometimes, we may require to find the value of x corresponding to a certain
value of y. The process of finding such a value of x is called inverse interpolation. In this
section, we shall study two methods of inverse interpolation:
1. Lagrange’s method.
2. Iterative method.
We apply Lagrange’s method when the arguments (xi ) are unequally spaced and the iterative
method when the arguments are equally spaced.

2.8.1 Lagrange’s Method

Choosing y as an independent variable and x and a dependent variable in Lagrange’s interpo-
lation.
(y − y1 )(y − y2 )(y − y3 ) . . . (y − yn )
x = f (y) = x0
(y0 − y1 )(y0 − y2 )(y0 − y3 ) . . . (y0 − yn )
(y − y0 )(y − y2 )(y − y3 ) . . . (y − yn )
+ x1
(y1 − y0 )(y1 − y2 )(y1 − y3 ) . . . (y1 − yn )
+.......................................
(y − y0 )(y − y1 )(y − y2 ) . . . (y − yn−1 )
+ xn . (22)
(yn − y0 )(yn − y1 )(yn − y2 ) . . . (yn − yn−1 )

31
Which is used for the inverse interpolation.

Example 16. The following table gives the values of x and y:

x 30 35 40 45 50
y 15.9 14.9 14.1 13.3 12.5
Find the value of x corresponding to y = 13.6.

Solution: Here,

y0 = 15.9, y1 = 14.9, y2 = 14.1, y3 = 13.3, y4 = 12.5, x0 = 30, x1 = 35, x2 = 40, x3 = 45, x4 = 50.

Taking y = 13.6 in Equation (22) above, we get

(13.6 − 14.9)(13.6 − 14.1)(13.6 − 13.3)(13.6 − 12.5)
x= (30)
(15.9 − 14.9)(15.9 − 14.1)(15.9 − 13.3)(15.9 − 12.5)
(13.6 − 15.9)(13.6 − 14.1)(13.6 − 13.3)(13.6 − 12.5)
+ (35)
(14.9 − 15.9)(14.9 − 14.1)(14.9 − 13.3)(14.9 − 12.5)
(13.6 − 15.9)(13.6 − 14.9)(13.6 − 13.3)(13.6 − 12.5)
+ (40)
(14.1 − 15.9)(14.1 − 14.9)(14.1 − 13.3)(14.1 − 12.5)
(13.6 − 15.9)(13.6 − 14.9)(13.6 − 14.1)(13.6 − 12.5)
+ (45)
(13.3 − 15.9)(13.3 − 14.9)(13.3 − 14.1)(13.3 − 12.5)
(13.6 − 15.9)(13.6 − 14.9)(13.6 − 14.1)(13.6 − 13.3)
+ (50)
(12.5 − 15.9)(12.5 − 14.9)(12.5 − 14.1)(12.5 − 13.3)
∴ x = 43.14185 =≈ 43.1.

2.8.2 Iterative Method

Let us consider Netwoton’s forward interpolation formula:
x(x − 1) 2 x(x − 1)(x − 2) 3
yx = y0 + x4y0 + 4 y0 + 4 y0 + . . . (23)
2! 3!
From this we get
 
1 x(x − 1) 2 x(x − 1)(x − 2) 3
x = yx − y0 − 4 y0 − 4 y0 − . . . . (24)
4y0 2! 3!
Negleting the second and higher differences, we get a first approximation for x and this we
write as:
1
x(1) = [yx − y0 ] . (25)
4y0
Substituting this in RHS of Equation (24), the second approximation is given by

x(1) (x(1) − 1) 2 x(1) (x(1) − 1)(x(1) − 2) 3

 
(2) 1
x = yx − y0 − 4 y0 − 4 y0 − . . . . (26)
4y0 2! 3!

The Third approximation is obtained by substituting x = x(2) in the RHS of Equation (24).

x(2) (x(2) − 1) 2 x(2) (x(2) − 1)(x(2) − 2) 3

 
(3) 1
x = yx − y0 − 4 y0 − 4 y0 − . . . . (27)
4y0 2! 3!

32
Similarly, further approximations can be obtained ad the process us to be continued until
two successive approximations of x agree with other to the required accuracy. If the final
approximation is x∗ , then the required is

x = x0 + x∗ h.

Note:

1. This method can be applied equally well by starting with any other interpolation formula.

2. It is powerful iterative procedure to find the root of an equation to a good degree of

accuracy.

Example 17. From the following table

x 1.8 2.0 2.2 2.4 2.6
Find x when y = 5, using iterative method.
y 2.9 3.6 4.4 5.5 6.7

Solution: The forward differences table is given as follows

x y 4 42 43 44
1.8 2.9
2.0 3.6 0.7
2.2 4.4 0.8 0.1
2.4 5.5 1.1 0.3 0.2
2.6 6.7 1.2 0.1 -0.2 -0.4

Here, yx = 5, y0 = 2.9, y1 = 3.6, y2 = 4.4, y3 = 5.5, y4 = 6.7, 4y0 = 0.7, 42 y0 = 0.1, 43 y0 =

0.2, 44 y0 = −0.4.

Let us consider Equation (25), we get

1 1
x(1) = [yx − y0 ] = (5 − 2.9) = 3.
4y0 0.7

Using Equation (26), then the second approximation is given as

x(1) (x(1) − 1) 2 x(1) (x(1) − 1)(x(1) − 2) 3

 
(2) 1
x = yx − y0 − 4 y0 − 4 y0 − . . .
4y0 2! 3!
 
1 3(3 − 1) 3(3 − 1)(3 − 2) 3(3 − 1)(3 − 2)(3 − 3)
= 5 − 2.9 − (0.1) − (0.2) − (−0.4)
0.7 2 3! 4!
1
= (5 − 2.9 − 0.3 − 0.2)
0.7
∴ x(2) = 2.2857143.

33
Proceeding on as explained above, we get

x(3) = 2.4724147
x(4) = 2.64364724
x(5) = 2.5373869
x(6) = 2.5985062
x(7) = 2.5611579
x(8) = 2.5841148
x(9) = 2.5700665
x(10) = 2.5786872
x(11) = 2.573406
x(12) = 2.5766447
x(13) = 2.5746598
x(14) = 2.5758768
x(15) = 2.575625.

Therefore ∴ x∗ = 2.575, correct to three decimal places. Hence, the value of x corresponding
to to yx = 5 is
x = x0 + x∗ h = 1.8 + (2.575)(0.2) = 2.315.
EXERCISE

1. Given values f (14) = 68.7, f (17) = 64, f (31) = 44, and f (35) = 39. Find f (27) using
Lagrange’s formula.

Given u1 = 22, u2 = 30, u3 = 82, u4 = 106 and u5 = 206. Find u6 using Lagrange’s interpola-
tion formula.

3. Apply Lagrange’s formula inversely to obtain the root of f (x) = 0, given that f (300 =
−30, f (34) = −13, f (38) = 3, and f (42) = 18.

4. Find x when f (x) = 0.163, given that

x 80 82 84 86 88
f (x) 0.131 0.154 0.176 0.200 0.221

34
2.9 NUMERICAL DIFFERENTIATION
2.9.1 Introduction
In this section, we consider the general polynomial of interpolation which states that given a
set of values (xi , yi ), i = 0, 1, 2, 3, . . . of two variables x and y, we obtain a polynomial φ(x),
called interpolation polynomial, such that φ(x) agrees with y(x) at all xi using Lagrange’s
formula and Newton’s formula.

In this section, we shall be concerned with the problem of numerical differentiation which states
that given a set of values of x and y(x) as previously, we wish to derive formula for computing
dy d2 y
the approximate values of dx , dx2 , . . . , for any x ∈ [x0 , xn ].

2.9.2 Numerical Differentiation Formula

The general method for deriving the numerical differentiation formula is to differentiate the
interpolating polynomial. Hence, corresponding to each of the formula derived in the section
on interpolation; We may derive a formula for the derivative.

We illustrate the derivation with Newton’s forward interpolation formula only, the method of
derivation being the same in each case.

(a) Newton’s forward interpolation formula is

1 1
y = y0 + p4y0 + p(p − 1)42 y0 + p(p − 1)(p − 2)43 y0 + . . . , (28)
2! 3!
x−x0
where x = x0 + ph, i.e p = h
. Then

dy dy dp
= By Chain rule.
dx dp dx

dy 1 1 1 1 1
= 4y0 + (2p − 1)42 y0 + [(p − 1)(p − 2) + p(p − 1) + p(p − 2)] 43 + . . .
dx  h 2! h 3! h
3p2 − 6p + 2 3 4p3 − 18p2 + 22p − 6 4

dy 1 2p − 1 2
∴ = 4y0 + 4 y0 + 4 y0 + 4 y0 + . . . (29)
.
dx h 2 6 24
dy
Expression (29) can be used for computing the value of dx for any non-tabular value of x in
[x0 , xn ]. For tabular values of x, the formula takes a simpler form. By setting x = x0 we obtain
p = 0 from p = x−x h
0
and hence
   
dy 1 1 2 1 3 1 4
= 4y0 − 4 y0 + 4 y0 − 4 y0 + . . . . (30)
dx x=x0 h 2 3 4

Differentiating (29) again, we obtain

d2 y 12p2 − 36p + 22 4
 
1 2 6p − 6 3
= 2 4 y0 + 4 y0 + 4 y0 + . . . . (31)
dx2 h 6 24

From which we obtain

d2 y
   
1 2 3 11 4
= 2 4 y0 − 4 y0 + 4 y0 + . . . . (32)
dx2 x=x0 h 12

35
Formula for computing higer derivatives may be obtained by successive differentiation. In a
similar way, different differention formula can be derived by starting with other interpolation
formula.

(b) Newton’s Backwork Interpolation formula gives

   
dy 1 1 2 1 3
= ∇yn + ∇ yn + ∇ yn + . . . . (33)
dx x=xn h 2 3

And
d2 y
   
1 2 3 11 4 5 5
= 2 ∇ yn + ∇ yn + ∇ yn + ∇ yn + . . . . (34)
dx2 x=xn h 12 6

Example 18. The population of a certain town (as obtain from census data) is shown in the
following table
year (x) 1951 1961 1971 1981 1991
Population in thousands (y) 19.96 36.65 58.81 77.21 94.61

Find the rate of grouth of the population in the year 1981.

Solution: We have to find the derivative at 1981 which is near the end of the table. Hence,
we use the derivative of Newton’s backward interpolation formula. The table of differences is
given as follows
x y ∇y ∇2 y ∇3 y ∇4 y
1951 19.96
1961 36.65 16.69
1971 58.81 22.16 5.47
1981 77.21 18.40 -3.76 -9.23
1991 94.61 17.40 -1 2.76 11.99

Here, h = 10, xn = 1991, ∇yn = 17.4, ∇2 yn = −1, ∇3 yn = 2.76, ∇4 yn = 11.99. So,

3p2 + 6p + 2 3 2p3 + 9p2 + 11p + 3 4

 
dy 1 2p + 1 2
= ∇yn + ∇ yn + ∇ yn + ∇ yn .
dx h 2 6 12

Now, we have to find out the rate of growth of the population in the year 1981. Then,
1981 − xn 1981 − 1991
xn + ph = 1981 ⇒ p = = = −1.
h 10
 
2(−1)+1 3(−1)2 +6(−1)+2 2(−1)3 +9(−1)2 +11(−1)+3
 dy  1
dx x=1981
= 10
17.4 + 2
(−1) + 6
(2.76) + 24
(11.99) = 1.6440833

dy d2 y
Example 19. Find the values of dx
and dx2
at x = 1.2 from the following table.

x 1.0 1.2 1.4 1.6 1.8 2.0 2.2

y 2.7183 3.3201 4.0552 4.9530 6.0496 7.3891 9.025

Soluation: The difference table is

36
 x y 4y 42 y 43 y 44 y 45 y 46 y
1.0 2.7183
1.2 3.3201 0.6018
1.4 4.0552 0.7351 0.1333
1.6 4.9530 0.8978 0.1627 0.0294
1.8 8.0496 1.0966 0.1988 0.0361 0.0067
2.0 7.3891 1.3395 0.2429 0.0441 0.0080 0.0013
2.2 9.025 0.2964 0.2964 0.0535 0.0094 0.0014 0.0001
Here,
x − x0 1.2 − 1
x = 1.2, h = 0.2, p = = = 1.
h 0.2

3p2 − 6p + 2 3 2p3 − 9p2 + 11p − 3 4

 
dy 1 2p − 1 2
= 4y0 + 4 y0 + 4 y0 + 4 y0 + . . .
dx x=1.2 h 2 6 12
 
1 1 1 1 1 1
= 0.618 + (0.1333) − (0.0294) + (0.0067) + (0.0013) + (0.0001)
0.2 2 6 12 20 30
= 3.2305.

On the other hand, by setting x0 = 1.2 and y0 = 3.3201, we get

 
dy 1 1 2 1 3 1 4
= 4y0 − 4 y0 + 4 y0 − 4 y0 + . . .
dx x=x0 h 2 3 4
 
dy 1 1 1 1
= 3.3201 − (0.1627) + (0.0361) −
dx x=1.2 0.2 2 3 4(0.0080) + 15 (0.0014)
= 3.3205

d2 y
 
1 2 3 11 4 5 5
= 2 4 y0 − 4 y0 + 4 y0 − 4 y0 + . . .
dx2 x=x0 h 12 6
2
 
dy 1 11 5
= 0.1627 − 0.361 + (0.0080) − (0.0014)
dx2 x=1.2 0.22 22 6
= 3.318.

EXERCISES

1. The table given below reveals velocity v of a body during time t specified. Find its acceler-
ation at t = 1.1.
t 1.0 1.1 1.2 1.3 1.4
v 43.1 47.7 52.1 56.4 60.8

2. For the following pairs of values of x and y, find numerically the first derivative at x = 4.
x 1 2 4 8 10
y 0 1 5 21 27

37
2.9.3 Maximum and Minimum values of a tabulated function
The maximum and minimum values of a function can be found by equating the first derivative
to zero and solving the variable. The same procedure can be applied to determine the maximum
and minimum of a tabulated function.

We consider Newton’s forward difference formula:

p(p − 1) 2 p(p − 1)(p − 2) 3
y = y0 + p4y0 + 4 y0 + 4 y0 + . . .
2 6
Differentiating with respect to p, we obtain
dy 2p − 1 2 3p2 − 6p + 2 3
= 4y0 + 4 y0 + 4 y0 + . . .
dp 2 6
dy
For maxima or minima, dp
= 0.

Determinating the right-hand side after third order difference and equating it to zero, we obtain
 
1 2 1 3
1
4y0 − 4 y0 + 4 y0 + 42 y0 − 43 y0 + (42 y0 )p2 = 0
2 3 2
which is a quadratic equation in p and can be solved. The values of x can be found from the
relation
x = x0 + ph.
Example 20. Given the following data, find the maximum value of y.
x -1 1 2 3
y -21 15 12 3
Solution: Since the arguments are not equispaced, we will form the divided difference table as
follows
x y 1st divided difference 2nd divided difference 3rd divided difference
-1 -21
1 15 18
2 12 -3 -7
3 3 -9 -3 1
Using Newton’s divided difference formula, we get
y = y0 + (x − x0 )[x0 , x1 ] + (x − x0 )(x − x1 )[x0 , x1 , x2 ] + (x − x0 )(x − x1 )(x − x2 )[x0 , x1 , x2 , x3 ]
y = −21 + (x + 1)(18) + (x + 1)(x − 1)(−7) + (x + 1)(x − 1)(x − 2)(1)
∴ y = x3 − 9x2 + 17x + 6
For maximum,
dy d2 y
= 0 and < 0.
dx dx2
That is
( (
3x2 − 18x + 17 = 0 x = 4.8257 or x = 1.743
⇒ ⇒ x = 1.743
6x − 18 < 0 x<3
Thus, the maximum value is
y(1.743) = (1.743)3 − 9(1.743)2 + 17(1.743) + 6 = 15.171612.

38
Example 21. Determine the maxima and minima of the function y = f (x) tabulated below.

x 0 1 2 3
y -5 -7 -3 13
Solution: The difference table is
x y 4 42 43
0 -5
1 -7 -2
2 -3 4 6
3 13 16 12 6
x−x0
Here, x0 = 0, h = 1, then p = h
= x. Now,

p(p − 1) 2 p(p − 1)(p − 2) 3

y = y0 + p4y0 + 4 y0 + 4 y0 + . . . ,
2 6
Hence
x(x − 1) x(x − 1)(x − 2)
y = −5 + x(−2) + (6) + (6)
2 6
y = x3 − 3x − 5.

In order to find the maxima and minima of the function, we have to solve the following equations
( (
dy dy
dx
= 0 dx
= 3x2 − 3 = 0
2
d y ⇒ 2
d y ⇒ x = ±1
dx2 dx2
= 6x

So there is a maxima at x = −1 and a minima at x = 1.

y(−1) = (−1)3 − 3(−1) − 5 = −3 Maxima

3
y(1) = (1) − 3(1) − 5 = −7 Minima

EXERCISES

1. Find the maximum and minimum values of the function from the following table
x 0 1 2 3 4 5
y 0 0.25 0 2.25 16.00 56.25
dy d2 y
2. Evaluate dx
and dx2
for the following table

x 0 1 3 6
y 18 10 -18 90
3. For the following data gives corresponding values of pressure and specific volume of super-
heated steam.
v 2 4 6 8 10
p 105 42.07 25.3 16.7 13

(i) Find the rate of change of pressure with respect to volume when v = 2.

(ii) Find the rate of change of volume with respect to pressure when p = 105.

39
2.10 NUMERICAL INTEGRATION
We are required to compute the value of the definite integral
Z b
I = y dx, (35)
a

where the function y = f (x) is defined by the set of data points (xi , yi ), i = 0, 1, 2, 3, . . . , n.
To derive a general formula for numerical integration, we substitute Newton’s forward
difference formula for y(x) in Equation (35), where the interval [a, b] is divided into n equal
sub-intervals such that
a = x0 < x1 < x2 < . . . < xn = b.
We obtain
Z xn  
1 2 1 3
I = y0 + p4y0 + p(p − 1)4 y0 + p(p − 1)(p − 2)4 y0 + . . . dx, (36)
x0 2 6
where putting x = x0 + ph. So dx = h dp, we get
Z n  
1 2 1 3
I = h y0 + p4y0 + p(p − 1)4 y0 + p(p − 1)(p − 2)4 y0 + . . . dp,
0 2 6
2 2
n2 (n − 2)2 3
 
n n (2n − 3) 2
= h ny0 + 4y0 + 4 y0 + 4 y0 + . . . .
2 12 24
Therefore,
xn
n(n − 2)2 3
 
n(2n − 3) 2
Z
n
y dx = nh y0 + 4y0 + 4 y0 + 4 y0 + . . . . (37)
x0 2 12 24

The equation (37) is the general quadratic formula and is called Newton’s-cote’s formula.

Different numerical integration can be derived by setting n = 1, 2, 3, . . . , in Equation (37)

2.10.1 For n = 1: Trapezoidal Rule

We define that differences of second and higher order vanish, and we obtain
Z x1  
1
y dx = h y0 + 4y0
x0 2
h
= (2y0 + 4y0 )
2
h
= (2y0 + (y1 − y0 ))
2
= (y0 + y1 )
Z x1
h
∴ y dx = (y0 + y1 ).
x0 2
For subsequent interval, we similary obtain
Z x2
h
y dx = (y1 + y2 )
x1 2
Z x3
h
y dx = (y2 + y3 )
x2 2

40
Until
Z xn
h
y dx = (yn−1 + yn )
xn−1 2

Combining the above formulae, we obtain the rule

Z xn
h
y dx = [y0 + 2(y1 + y2 + y3 + . . . + yn−1 ) + yn ] , (38)
x0 2

which is known as the Trapezoidal rule.

1
2.10.2 For n = 2: Simpson 3
Rule
Setting n = 2 in Equation (37), we observe that differences of third and higher order vanish,
and we obtain
Z x2  
2 2(2.2 − 3 2
y dx = 2h y0 + 4y0 + 4 y0
x0 2 12
 
1 2
= 2h y0 + 4y0 + 4 y0
6
h
6y0 + 64y0 + 42 y0
 
=
3
h
= [6y0 + 6(y1 − y0 ) + (y2 − 2y1 + y0 )]
Z x2 3
h
∴ y dx = [y0 + 4y1 + y2 ] .
x0 3

Similarly,
Z x4
h
y dx = [y2 + 4y3 + y4 ]
x2 3
.
.
Z xn .
h
y dx = [yn−2 + 4yn−1 + yn ] .
xn−2 3

Summing up the above results, we obtain

Z xn
h
y dx = [y0 + 4(y1 + y3 + . . . + yn−1 ) + 2(y2 + y4 + . . . + yn−2 ) + yn ] , (39)
x0 3
1
which is known as Simpson’s 3
rule, and requires an even number of subintervals of width h.

3
2.10.3 For n = 3: Simpson’s 8
Rule
Setting n = 3 in Eqaution (37) and negleting all differences above the third order, we get
Z x3  
3 3 2 1 3
y dx = 3h y0 + 4y0 + 4 y0 + 4 y0
x0 2 4 8
3h
= [y0 + 3y1 + 3y2 + y3 ] .
8
41
Similarly,
Z x6
3h
y dx = [y3 + 3y4 + 3y5 + y6 ]
x3 8
.
.
.
Z xn
3h
y dx = [yn−3 + 3yn−2 + 3yn−1 + yn ]
xn−3 8

Adding up all these integrals, we obtain

Z xn
3h
y dx = [y0 + 3(y1 + y2 + y4 + y5 + . . . + yn−2 + yn−1 ) + 2(y3 + y6 + . . . + yn−3 ) + yn ] ,
x0 8
(40)

which is called the Simpson’s 83 rule.

Example 22. Evaluate Z 10

dx
by using
0 1 + x2
1. Trapezoidal rule.
1
2. Simpson’s 3
rule.
3
3. Simpson’s 8
rule.

Solution: Diving the whole range of integration [0, 10] into 10 equal parts; h = 1 and the value
of integral for each point of subdivision are given below.
x y
x0 = 0 y0 = 1
x1 = 1 y1 = 12 = 0.5
x2 = 2 y2 = 15 = 0.2
1
x3 = 3 y3 = 10 = 0.1
1
x4 = 4 y4 = 17 = 0.0588235
1
x5 = 5 Y5 = 26 = 0.0384615
1
x6 = 6 y6 = 37 = 0.027027
1
x7 = 7 y7 = 50 = 0.02
1
x8 = 8 y8 = 65 = 0.0153846
1
x9 = 9 y9 = 82 = 0.0121951
1
x10 = 10 y10 = 101

1. By Trapezoidal Rule
Z 10
h
y dx = [y0 + 2(y1 + y2 + y3 + . . . + y9 ) + y10 ] ,
x0 2
 
1 1
= 1 + 2(0.5 + 0.2 + 0.1 + . . . + 0.0121951) +
2 101
= 1.4768422.

42
 1
2. By Simpson’s 3
rule
Z 10
h
y dx = [y0 + 4(y1 + y3 + . . . + y9 ) + 2(y2 + y4 + . . . + y8 ) + y10 ] ,
0 3
= 1.4316659.

3
3. Simpson’s 8
rule
Z 10
3h
y dx = [y0 + 3(y1 + y2 + y4 + y5 + y7 + y8 ) + 2(y3 + y6 + . . . + y9 ) + y10 ] ,
0 8
= 1.4198828.

Example 23. The velocity v of a particle at distance s from a point on its path is given by the
following table:
s[f t] 0 10 20 30 40 50 60
v[t/s] 47 58 64 65 61 52 38
1
Estimate the time taken to travel 60f t using Simpson’s 3
rule. Compare the result with Simp-
son’s 83 rule.

Solution:
ds
v= = rate of displacement is velocity
dt
So Z 60 Z 60
1
t= ds = y dx = time taken to travel 60 ft.
0 v 0
1
It means that y = v
and s = x.

s[x] y = 1/v
x0 = 0 y0 = 14 = 0.0212765
1
x1 = 10 y1 = 58 = 0.0172413
1
x2 = 20 y2 = 64 = 0.015625
1
x3 = 30 y3 = 65 = 0.0153846
1
x4 = 40 y4 = 61 = 0.0163934
1
x5 = 50 y5 = 52 = 0.192307
1
x6 = 60 y6 = 38 = 0.0263157
1
Here h = 10; By simpson’s 3
rule, we get
Z 60
h
y dx = (y0 + 2(y2 + y4 ) + 4(y1 + y3 + y5 ) + y6 )
0 3
1.063502.
3
By Simpson’s 8
rule, we get
Z 60
3h
y dx = (y0 + 3(y1 + y2 + y4 + y5 ) + 2(y3 ) + y6 )
0 8
= 1.0643723.

43
 EXERCISES:
R2
1. Evaluate 0
y dx from the following table using Trapezoidal rule

x 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8

y 1.21 1.37 1.46 1.59 1.67 2.31 2.91 3.83 4.01 4.79
1
2. Find an approximate value of ln 5 by calculating to four decimal places by Simpson’s 3
rule
R 5 dx
the integral 0 4x+5 dividing the range into 10 equal parts.

[Link] the approximate distance traveled by a train between 11 : 50 a.m to 12 : 30 p.m

from the following data using Simpson’s 1/3 rule.

time 11:50 12:00 12:10 12:20 12:30

speed 24.2 35.0 41.3 42.8 39.2

4. Whe an a train is moving at 30 miles an hour, steam is shut off and brakes are applied. The
speed of the train in miles per hour after t seconds is given by
t 0 5 10 15 20 25 30 35 40
v 30 24 19.5 16 13.6 11.7 10.0 8.5 7.0

Determine how far the train has moved in the 40 seconds.

44
3 Numerical Solution to Ordinary Differential Equations
3.1 Power Series Solution
Consider the differential equation
dy
= y 0 = f (x, y), (41)
dx
subject to the condition y(x0 ) = y0 .
This is an initial value problem. Now we can expand the solution to equation (41) i.e y(x) in
the neighborhood of x = x0 in power series as

y(x) = A0 + A1 (x − x0 ) + A2 (x − x0 )2 + . . . , (42)

where A0 , A1 , A2 , . . . are constant to be determineed such that Equation (42) satisfies Equation
(41) subject to the initial value.

To expand y(x) in power series, generally, Taylor’s Theorem or Maclaurin’s Theorem is used.
If the boundary condition x0 (6= 0), then we use Taylor’s series to expand y(x) about x = x0
and it is
(x − x0 ) 0 (x − x0 )2 00
y(x) = y(x0 ) + y (x0 ) + y (x0 ) + . . . (43)
1! 2!
If the boundary condition is at x0 = 0, then we use Maclaurin’s series to expand y(x) about
x = 0 and it is
x 0 x2 x3 000
y(x) = y(0) + y (0) + y 00 (0) + !y (0) + . . . (44)
1! 2! 3
The method is best illustrated by the following example.
Example 24. Evaluate the solution to the differential equation y 0 = (1 + x)xy 2 subject to
y(0) = 1 by taking five terms in Maclaurin’s series for x = 0 (0.1) 0.4. Compare with the exact
solution.

Solution: Maclaurin’s series is

x 0 x2 x3 000 x4
y(x) = y(0) + y (0) + y 00 (0) + !y (0) + y (4) + . . . (45)
1! 2! 3 4!
The given equation is
y 0 (x) = (1 + x)xy 2 , y(0) = 1.

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

y 00 (x) = 2xyy 0 + y 2 + 2xy 2 + 2x2 yy 0 , y 00 (0) = 1
000
y (x) = 4yy 0 + 8xyy 0 + 2xy 02 + 2y 2 + 2xyy 00 + 2x2 yy 00 + 2x2 y 02 , y 000 (0) = 2
y (4) = 2yy 02 + 6y 02 + 6xy 0 y 00 + 6yy 00 + 12yy 0 + 14xyy 00 + 12yy 02 + 6x2 y 0 y 00 + 2x2 yy 00 , y (4) (0) = 6

Putting these values in Equation (45), we get

x x2 x3 x4
y(x) = 1 + (0) + (1) + (2) + (6)
1 2 3! 4!
x2 x3 x4
∴ y(x) = 1 + + + + .... (46)
2 3 4
45
Putting x = 0.1, 0.2, 0.3, 0.4 successively in Equation (46), we get

(0.1)2 (0.1)3 0.1)4

y1 = y(0.1) = 1 + + + = 1.0053583.
2 3 4
(0.2)2 (0.2)3 (o.2)4
y2 = y(0.2) = 1 + + + = 1.0230667.
2 3 4
(0.3)2 (0.3)3 (0.3)4
y3 = y(0.3) = 1 + + + = 1.056025.
2 3 4
(0.4)2 (0.4)3 (0.4)4
y4 = y(0.4) = 1 + + + = 1.1077333.
2 3 4
Exact Solution: Given that
dy
= (1 + x)xy = (x + x2 )y 2 .
dx
It means that
dy
= (x + x2 )dx by separation of variables
y2
 2
x3

1 x
⇒− = + + c by integration
y 2 3

Using the initial condition y(0) = 1 i.e y = 1 at x = 0, in the above equation, we get

−1 = c ⇒ c = −1

Therefore,

x2 x3
 
1
− = + −1
y 2 3
3x2 + 2x3
= −1
6
3x2 + 2x3 − 6
=
6
6
∴y = . (47)
6 − 3x2 − 2x3
Now putting x = 0.1, 0.2, 0.3, 0.4 successively in Equation (47), we get

y1 = y(0.1) = 1.0053619.
y2 = y(0.2) = 1.0231924.
y3 = y(0.3) = 1.0570825.
y4 = y(0.4) = 1.1127596.

Exercise: Using four terms of the Maclaurin’s series, find y at x = 0 (0.1) 0.6 given that
2y 0 = (1 + x)y 2 , y(0) = 1. Compare the results with the exact solution.

46
3.2 Euler’s Method
Consider the differential equation
dy
= f (x, y), (48)
dx
where y(x0 ) = y0 .
Suppose that we wish to find successively y1 , y2 , y3 , . . . , ym where ym os the value of y corre-
sponding to x = xm , where xm = x0 + mh, m = 1, 2, 3, . . . , h being small. Here we use the
property that in a small interval, a curve is nearly a straight line.

Thus, in the interval x0 to x1 of x, we approximate the curve by the tangent at the point
(x0 , y0 ). Therefore, the equation of the tangent at (x0 , y0 ) is
 
dy
y − y0 = (x − x0 )
dx (x0 , y0 )
= f (x0 , y0 )(x − x0 ). [From the equation (48)]
∴ y = y0 + (x − x0 )f (x0 , y0 ).
Hence, the value of y corresponding to x = x1 is
y1 = y0 + (x1 − x0 )f (x0 , y0 )
∴ y1 = y0 + h f (x0 , y0 ). (49)
Similarly, approximating the curve in the next interval [x1 , y1 ] by a line through (x1 , y1 ) with
the slope f (x1 , y1 ), we get
y2 = y1 + h f (x1 , y1 ). (50)
Proceeding on; In general it can be shown that
ym+1 = ym + h f (xm , ym ). (51)
Remark: The process is very slow and to obtain it with reasonable accuracy using
Euler Method. We have to use h very small.

3.3 Improved Euler’s Method

Here we consider a line passing through A(x0 , y0 ) whose slope is the average of the slopes at
(1)
A(x0 , y0 ) and P (x1 , y1 ) such that
(1)
y1 = y0 + h f (x0 , y0 ).
Therefore, the equation of the tangent line is given as
1 (1)

y − y0 = (x − x0 ) f (x0 , y0 ) + f (x1 , y1 ) . (52)
2
As we assuming that the line passes through the point (x0 , y0 ) and (x1 , y1 ). Therefore, the
equation is
1 
(1)

y1 = y0 + (x1 − x0 ) f (x0 , y0 ) + f (x1 , y1 )
2
h (1)

or y1 = y0 + f (x0 , y0 ) + f (x1 , y1 ) .
2
47
Therefore, we have
h
y1 = y0 + (f (x0 , y0 ) + f [x0 + h, y0 + h f (x0 , y0 )]) . (53)
2
In general, we have the formula
h
ym+1 = ym + (f (xm , ym ) + f [xm + h, ym + h f (xm , ym )]) . (54)
2
Where xm − xm−1 = h.

3.4 Modified Euler’s Method

In this method the curve in the interval [x0 , y0 ], where x1 = x0 + h is approximated by the line
through (x0 , y0 ) with slope
 
h h
f x0 + , y0 + f (x0 , y0 . (55)
2 2

That is, the slope at the middle point whose abscisse is the average of x0 and x0 + h2 . The
equation of line L is
  
h h
y − y0 = (x − x0 ) f x0 + , y0 + f (x0 , y0 . (56)
2 2

Putting x = x1 , we get
  
h h
y1 = y0 + (x1 − x0 ) f x0 + , y0 + f (x0 , y0
2 2
  
h h
∴ y1 = y0 + h f x0 + , y0 + f (x0 , y0 . (57)
2 2

Proceeding in the same way, the general formula can be show that
  
h h
ym+1 = ym + h f xm + , ym + f (xm , ym . (58)
2 2
dy
Example 25. Solve dx = 1 − y, y(0) = 0 in the range 0 ≤ x ≤ 0.3 using the following
methods by choosing h = 0.1

1. Euler’s Method.

2. Improved Euler’s Method.

3. Modified Euler’s Method by.

Compare the answers with the exact solution.

Solution: Given that

dy
= 1 − y, y(0) = 0, h = 0.1.
dx
Now, we have to find out the solutions at x = 0.1, x = 0.2 and x = 0.3.

48
 dy
1. Euler’s Method: The algorithm is if dx
= f (x, y), y(x0 ) = y0 , then

ym+1 = ym + hf (xm , ym ).

Here f (x, y) = 1 − y, h = 0.1, x0 = 0 and y(0) = y0 = 0. Hence

ym+1 = ym + (0.1)(1 − ym )
⇒ ym+1 = 0.1 + 0.9ym .

Putting m = 0, 1, 2 successively, we get

y1 = 0.1 + 0.9y0 = 0.1 + 0.9(0) = 0.1

y2 = 0.1 + 0.9y1 = 0.1 + 0.9(0.1) = 0.19
y3 = 0.1 + 0.9y2 = 0.1 + 0.9(0.19) = 0.271.

2. Improved Euler’s Method: Here, the formula is

h (1)

ym+1 = ym + f (xm , ym ) + f [xm+1 , ym+1 )] ,
2
where
(1)
ym+1 = f [xm + h, ym + h f (xm , ym )].
Here, f (x, y) = 1 − y, therefore f (xm , ym ) = 1 − ym . Hence,
   
(1)
f xm+1 , ym+1 = 1 − ym + hf (xm , ym )

= 1 − ym − h(1 − ym )
 
(1)
∴ f xm+1 , ym+1 = (1 − h)(1 − ym )

So
h
ym+1 = ym + [1 − ym + (1 − h)(1 − ym )]
2
h
= ym + [(1 − ym )(1 + 1 − h)]
2
h
⇒ ym+1 = ym + [(2 − h)(1 − ym )] .
2
Here, h = 0.1

ym+1 = ym + 0.095(1 − ym )
∴ ym+1 = 0.095 + 0.905ym .

Putting m = 0, 1, 2 successively in the above equation, we get

y1 = 0.095 + 0.905y0 = 0.095 + 0.905(0) = 0.095

y2 = 0.095 + 0.905y1 = 0.095 + 0.905(0.095) = 0.180975
y3 = 0.095 + 0.905y2 = 0.095 + 0.905(0.180975) = 0.2587823.

49
 3. Modified Euler’s Method: The formula is
  
h h
ym+1 = ym + h f xm + , ym + f (xm , ym
2 2
  
h
= ym + h 1 − ym + f (xm , ym )
2
 
h
= ym + h 1 − ym − (1 − ym )
2
 
h h
= ym + h 1 − ym − + ym
2 2
 
h h
= ym + h (1 − ) − (1 − )ym
2 2
 
h
= ym + h 1 − (1 − ym )
2
 
0.1
= ym + 0.1 1 − (1 − ym )
2
∴ ym+1 = 0.095 + 0.905ym .

Which is identical to the equation of the Improved Euler’s Method. Putting successively
m = 0, 1, 2, we get

y1 = 0.095 + 0.905y0 = 0.095 + 0.905(0) = 0.095

y2 = 0.095 + 0.905y1 = 0.095 + 0.905(0.095) = 0.180975
y3 = 0.095 + 0.905y2 = 0.095 + 0.905(0.180975) = 0.2587823.

The exact solution: we have

dy
= 1 − y Given differential equation
dx
dy
⇒ = dx by separation of variables
1−y
⇒ − log(1 − y) + log c = x by integration
c
⇒ ex = by separating the variables
1−y
⇒ ex (1 − y) = c by cross product.

By using the given initial condition, we have y = 0 at x = 0. so c = 1, therefore

y = 1 − ex .
y1 = y(0.1) = 0.0951625.
y2 = y(0.2) = 0.1812691.
y3 = y(0, 3) = 0.2591817.

EXERCSISE: Solve
dy 2x
=y− , y(0) = 1
dx y
in the range 0 ≤ x ≤ 0.2 using (i) Euler’s Method, (ii) Improved Euler’s Method, and
(iii) Modified Euler’s Method. Take h = 0.1.

50
3.5 Runge - Kutta Method
Runge - Kutta methods do not required the calculations of higher order derivatives and
give great accuracy.

3.5.1 First order Runge-Kutta method

We have see that Euler’s method gives

y1 = y + 0 + hf (x0 , y0 ) = y0 + hy00 [y 0 = f (x, y)]. (59)

Now y1 = f (x0 + h. Expanding it by Taylor’s series, we get

h 0 h2 00
y1 = y0 + y + y + ...
1! 0 2! 0
This implies that Euler’s method agrees with Taylor’s series solution up to term in h. Hence,
Runge-Kutta method of first order is Euler’s method

3.5.2 Runge-Kutta method of second order

The improved Euler’s method gives
h (1)

y1 = y0 + f (x0 , y0 ) + f [x0 + h, y1 )] ,
2
where
(1)
y1 = f [x0 + h, y0 + h f (x0 , y0 )].
That is
h
y1 = y0 + (f (x0 , y0 ) + f [x0 + h, y0 + hf (x0 , y0 )]) , (60)
2
Now, by Taylor’s series

h2 00 h3 000
y1 = y(x0 + h) = y0 + hy00 + y + y0 + . . . (61)
2! 0 3!
Expanding f((x0 + h, y0 + hf (x0 , y0 )) by Taylor’s series for a function of two variables, we get
 
h ∂f ∂


y1 = y0 + 2 f (x0 , y0 ) + f (x0 , y0 ) + h ∂x x0 ,y0 + 1 + hf (x0 , y0 ) ∂y
x0 ,y0
+terms in 2nd and higher power of h
"     !#
1 ∂f ∂
y1 = y0 + 2hf (x0 , y0 ) + h2 + f (x0 , y0 ) + O(h3 )
2 ∂x x0 ,y0 ∂y x0 ,y0
2
h 0
= y0 + hf (x0 , y0 ) + f (x0 , y0 ) + O(h3 )
2
h2 00
= y0 + hf (x0 , y0 ) + y0 + O(h3 )
2
Then,

h2 00
y1 = y0 + hf (x0 , y0 ) + y + O(h3 ). (62)
2 0
51
Therefore, the Improved Euler’s method is the Runge-Kutta method of second order.

The algorithm to Runge - Kutta method of second order.

k1 = hf (x0 , y0 )
k2 = hf (x0 + h, y0 + k1 )
1
k = (K1 + k2 )
2
∴ y1 = y0 + k.

3.5.3 Higher order Runge - Kutta Methods

Runge - Kutta method of third order: It is defined by the following equations
k1 = hf (x0 , y0 )
k2 = hf (x0 + h, y0 + k1 )
k3 = hf (x0 + h, y0 + k2 )
h k1
k4 = hf (x0 + , y0 + )
2 2
1
k = (k1 + 4k4 + k3 )
6
∴ y1 = y0 + k.

3.5.4 Runge - Kutta Method of fourth order

It is most commonly known as Runge - Kutta method and the working procedure is as
follows: Consider the following equation
dy
= f (x, y), y(x0 ) = y0 .
dx
To compute y1 , we calculate successively
k1 = hf (x0 , y0 )
 
h k1
k2 = hf x0 + , y0 +
2 2
 
h k2
k3 = hf x0 + , y0 +
2 2
k4 = hf (x0 + h, y0 + k3 )
1
k = (k1 + 2k2 + 3k3 + k4 )
6
∴ y1 = y0 + k and x1 = x0 + h.
The increment in y in second interval is computed in a similar manner by means of the formula
k1 = hf (x1 , y1 )
 
h k1
k2 = hf x1 + , y1 +
2 2
 
h k2
k3 = hf x1 + , y1 +
2 2
k4 = hf (x1 + h, y1 + k3 )
1
k = (k1 + 2k2 + 3k3 + k4 )
6
∴ y2 = y1 + k and x2 = x1 + h.

52
And so on for succeeding intervals.

Example 26. Given y 0 = x2 − y, y(0) = 1, find y(0.1), y(0.2) using

1. Second order Runge - Kutta method

2. Third orde Runge - Kutta method

3. Fourth order Runge - Kutta method

Solution: Given
y 0 = x2 − y, x0 = 0, y0 = 1, ∴ f (x0 , y0 ) = −1.
Let h = 0.1

1. Runge - kutta method of second order: Here,

k1 = hf (x0 , y0 ) = 0.1(−1) = −0.1

k2 = hf (x0 + h, y0 + k1 ) = 0.1f (0.1, 0.9) = 0.1[(0.1)2 − 0.9] = −0.089
1 1
k = (K1 + k2 ) = (−0.1 − 0.089) = −0.0945.
2 2
∴ y1 = y0 + k = 1 − 0.0945 = 0.9055.

Again, taking x1 = 0.1, y1 = 0.9055 in place of x0 , y0 and repeating the process, we get

k1 = hf (x1 , y1 ) = h(x21 − y) = 0.1[(0.1)2 − 0.9055) = −0.08955

k2 = hf (x1 + h, y1 + k1 ) = hf (0.2, 0.81595) = (0.1)[(0.2)2 − 0.81595) = −0.077595
1 1
k = (K1 + k2 ) = (−0.08955 − 0.077595) = −0.0835725
2 2
∴ y2 = y1 + k = 0.9055 − 0.0835725 = 0.8219275.

2. Runge - Kutta method of Third order: Here, we have

k1 = hf (x0 , y0 ) = −0.1
k2 = hf (x0 + h, y0 + k1 ) = −0.089
k3 = hf (x0 + h, y0 + k2 ) = −0.0901
h k1
k4 = hf (x0 + , y0 + ) = −0.09475
2 2
1
k = (k1 + 4k4 + k3 ) = −0.9485
6
∴ y1 = y0 + k = 0.90515.

Taking x1 = 0.1, y1 = 0.90515, h = 0.1 in place of x0 , y0 and repeating the process,

we get

k1 = hf (x1 , y1 ) = −0.089515
k2 = hf (x1 + h, y1 + k1 ) = −0.0775635
k3 = hf (x1 + h, y1 + k2 ) = −0.0787586
h k1
k4 = hf (x1 + , y1 + ) = −0.0837892
2 2
1
k = (k1 + 4k4 + k3 ) = −0.0839051
6
∴ y2 = y1 + k = 0.82124490

53
 3. Runge - Kutta method of fourth order: Here,

k1 = hf (x0 , y0 ) = −0.1
 
h k1
k2 = hf x0 + , y0 + = −0.09475
2 2
 
h k2
k3 = hf x0 + , y0 + = −0.0550125
2 2
k4 = hf (x0 + h, y0 + k3 ) = −0.0894987
1
k = (k1 + 2k2 + 3k3 + k4 ) = −0.0948372
6
∴ y1 = y0 + k = 0.9051627.

Taking x1 = 0.1, y1 = 0.9051627 in place of x0 , y0 and repeating the process, we get

k1 = hf (x1 , y1 ) = −0.0895162
 
h k1
k2 = hf x1 + , y1 + = −0.0837904
2 2
 
h k2
k3 = hf x1 + , y1 + = −0.0840767
2 2
k4 = hf (x1 + h, y1 + k3 ) = −0.0781085
1
k = (k1 + 2k2 + 3k3 + k4 ) = −0.083893
6
∴ y2 = y1 + k = 0.8212696.

Exercise:

1. Using Runge -Kutta method of fourth order, solve for y(0.1), y(0.2) and y(0.3) given
that
y 0 = xy + y 2 , y(0) = 1.

54
4 Simultaneous Linear Algebraic Equations
4.1 Introduction
Simultaneous linear algebraic equations are very common in various fields of engineering and
science. We use matrix inversion method or cramer’s rule to solve these equations in general.
To solve the equations contain a large number of unknowns; Particularly suited for
computer operations. These are of two types: Direct and Iterative methods.

4.2 Direct Method

4.2.1 Gauss Elimination Method
In this method, the unknowns are eliminated successively by transforming the given matrix
into an equivalent system with upper triangular coefficient matrix (i.e a matrix in which all the
elements below the principal diagonal are zero) by means of elementary row operations, from
which the unknowns are found by back substitution. Here, we shall explain it by considering
of three equations in three unknowns.

Consider a system

a1 x + b1 y + c1 z = d1

a2 x + b 2 y + c 2 z = d 2 (63)

a3 x + b 3 y + c 3 z = d 3


where x, y, z are unknowns. The system in matrix form is

AX = B,

where
     
a1 b 1 c 1 x d1
A = a2 b2 c2  , X = y
 and B = d2  .

a3 b3 c13 z d3

Consider the augmanted matrix [A|b],

 
a1 b1 c1 d1
[A|B] = a2 b2 c2 d2  (64)
a3 b3 c13 d3

Now, Equation (64) is to be reduced to an upper triangular matrix, let a1 6= 0. Then

 
(
a2 a1 b 1 c 1 d 1
R2 ⇒ R2 − a1 R1
a3
⇔  0 b02 c02 d02  (65)
R3 ⇒ R3 − a1 R1 0 b0 c 0 d 0 3 13 3

Here, a1 is called the first pivot and b02 , c02 , b03 , c03 , d02 , d02 are transformed elements. Now take b02
as the pivot (b0 2 6= o), then
 
0 a1 b 1 c 1 d 1
b
R3 ⇒ R3 − 30 R2 ⇔  0 b02 c02 d02  (66)
b2 00 00
0 0 c13 d3

55
Now, if c003 6= 0, from Equation (66). The given system of linear equation is equivalent to

a1 x + b1 y + c1 z = d1

0 0 0
b2 y + c 2 z = d 2 (67)
 00
 00
c3 z = d3 .

Using back substitution,

00

d3


 z= 00
c3

1 0 00 0 00

y= 0 00
b2 c3
d2 c 3 − c 2 d3 (68)
 0 00 00 00 0 00 0 00
1


x =

0 00
a1 b2 c3
d 1 b 2 c 3 − b 1 d 2 c3 + b1 c 2 d 3 − b2 c 1 d 3

0 00
Note: This method fails if any one of the pivots a1 , b2 or c3 becomes zero. In such
case, by interchanging the rows we can get the non-zero pivots.
Example 27. Solve the system of equations

3x + y − z = 3

2x − 8y + z = −5 ,

x − 2y + 9z = 8


using Gauss elimination method.

Solution: The given system is equivalent to
    
3 1 −1 x 3
2 −8 1  y  = −5
1 −2 9 z 8

The augmented matrix is

 
3 1 −1 3
[A|B] = 2 −8 1 −5
1 −2 9 8

Now we will make A as upper triangular choosing 3 as a pivot,

 
(
2 3 1 −1 3
R2 ⇒ R2 − 3 R1
1
⇔ 0 − 26 3
5
3
−7
R3 ⇒ R3 − 3 R1 0 −3 7 28
7
3

Now choosing − 26
3
as the pivot from the second column, we get
 
3 1 −1 3
7
R3 ⇒ R3 − ⇔ 0 − 263
5
3
−7
26 693 231
0 0 78 26

From this, we get

  
3x + y − z = 3
 z = 1
 z = 1

− 26
3
y + 53 z = −7 − 26
3
y = −7 + 35 = − 26
3
⇒y=1 ∴ y=1

 693  
z = 231 1 1
x = 3 (3 − y + z) = 3 (3 − 1 + 1) = 1 x=1
 
78 26

56
 EXERCISES

1. Solve the system of equations


28x + 4y − z = 32

x + y + 10z = 24 By Gauss elimination method.

2x + 17y + 4z = 35


2. Using Gauss elimiantion method. Solve the system of equations


3.15x − 1.96y + 3.85z = 12.95

2.13x + 5.12y − 2.89z = −8.61

5.92x + 3.05y + 2.15z = 6.88


3. Solve the system of equations



 x1 + x 2 + x3 + x4 = 2

x + x + 3x − 2x = −6
1 2 3 4


 2x1 + 3x2 − x3 + 2x4 = 7
x1 + 2x2 + x3 − x4 = −2


4.2.2 Gauss - Jordan Method

This method is a modified form of Gauss elimination method. In this method, the coefficient
matrix A of AX = B is reduced to a diagonal matrix or unit matrix by making all the elemets
above and below to the principal diagonal of A zero. The labour of back substutution is saved
here even though it involves additional computation.

Example 28. Solve the system of equations


10x + y + z = 12

2x + 10y + z = 13 by Gauss - Jordan method.

x + y + 5x = 7


Solution: The given system in matrix form is

      
10 1 1 x 12 10 1 1 12
 2 10 1 y  = 13 ∴ [A|B] =  2 10 1 13 .
1 1 5 z 7 1 1 5 7

 
1 −8 −44 −51
R1 → R1 − 9R3 ⇔ 2 10 1 13 
1 1 5 7

 
( 1 −8 −44 −51
R2 → R2 − 2R1
⇔ 0 26 89 115 
R3 → R3 − R1 0 9 49 58

57
  
1 −8 −44 −51
R2 → −(R2 − 3R3 ) ⇔ 0 1 58 59 
0 9 49 58

 
( 1 0 420 421
R1 → R1 + 8R2
⇔ 0 1 58 59 
R3 → R3 − 9R2 0 0 −473 −473

 
1 0 420 421
1
R3 → − R3 ⇔ 0 1 58 59 
473
0 0 1 1

 
(
420 1 0 0 1
R1 → R1 + 473 R3
⇔ 0 1 0 1
R2 → R2 − 58R3 0 0 1 1

Therefore, the system AX = B reduced to the form


x = 1
    
1 0 0 x 1 
0 1 0 y  = 1 ∴ y=1
0 0 1 z 1

z=1


EXERCISES

1. Solve the system of equations


10x1 + x2 + x3 = 12,

x1 + 10x2 − x3 = 10, byGauss − Jordanmethod.

x1 − 2x2 + 10x3 = 9,


2. Solve the following system of equations by Gauss - Jordan method



 x + 2y + z − w = −2,

2x + 3y − z + 2w = 7,


 x + y + 3z − 2w = −6,

x + y + z + w = 2.

58
4.3 Iterative Method
The method of iteration is not applicable to all systems of equations. For this, each equa-
tion of the system must contain one large coefficient and large coefficient must be
attached to a different unknown in that equation. “ The solution to a system of linear
equations will exist by iterative procedure if the absolute value of largest coefficient is greater
than the sum of the absolute values of all remaining coefficients in each equation (condition
for convergence)”. Now let us see some of such method in details:

4.3.1 Jacobi Method of iteration

This method is also known as Gauss-Jacobi method. Here, consider the following system of
linear equations

a1 x + b 1 y + c 1 z = d 1
a2 x + b2 y + c2 z = d2 (69)
a3 x + b 3 y + c 3 z = d 3

Let
|a1 | > |b1 | + |c1 |; |b2 | > |a2 | + |c2 |; |c3 | > |a3 | + |b3 |.
That is, in each equation the coefficients of diagonal terms are large. Hence the system (70) is
ready for iteration. Solving for x, y and z respectively, we get
1
x = (d1 − b1 y − c1 z)
a1
1
y = (d2 − a2 x − c2 z)
b2
1
z = (d3 − a3 x − b3 y) (70)
c3
Let x0 , y0 and z0 be the initial approximations of the unknowns x, y and z. Substituting these
values on RHS of Equation (70), the first approximations are given by
1
x1 = (d1 − b1 y0 − c1 z0 )
a1
1
y1 = (d2 − a2 x0 − c2 z0 )
b2
1
z1 = (d3 − a3 x0 − b3 y0 ) (71)
c3
Substituting the values of x1 , y1 and z1 in the RHS of Equation (70), the second approximations
are given by
1
x2 = (d1 − b1 y1 − c1 z1 )
a1
1
y2 = (d2 − a2 x1 − c2 z1 ) (72)
b2
1
z2 = (d3 − a3 x1 − b3 y1 )
c3

59
Proceeding in the same way, if xr , yr , zr are the rth iterates, then
1
xr+1 = (d1 − b1 yr − c1 zr )
a1
1
yr+1 = (d2 − a2 xr − c2 zr )
b2
1
zr+1 = (d3 − a3 xr − b3 yr ) (73)
c3
The process continued till convergency is secured.

Note: In the absence of any better estimates, the initial approximations are taken
as x0 = y0 = z0 = 0.
Example 29. Solve by Jacobi iteration method the system

8x − 3y + 2z = 20;
6x + 3y + 12z = 35;
4x + 11y − z = 33.

Solution: Consider the given system as

8x − 3y + 2z = 20;
4x + 11y − z = 33;
6x + 3y + 12z = 35.

So that the diagonal elements are dominant in the coefficient matrix. Now we write the
equations in the form
1
x = (20 + 3y − 2z)
8
1
y = (33 − 4x + z)
11
1
z = (35 − 6x − 3y) (74)
2
We start from an approximations x0 = y0 = z0 = 0. Substituting these values on RHS of
equation (74).

First approximation is given as:

1
x1 = (20 + 3y0 − 2z0 ) = 2.5
8
1
y1 = (33 − 4x0 + z0 ) = 3
11
1
z1 = (35 − 6x0 − 3y0 ) = 2.9166667
2
Second approximation: substituting x1 , y1 , z1 on RHS of Equations (74), we get
1
x2 = (20 + 3y1 − 2z1 ) = 2.895833
8
1
y2 = (33 − 4x1 + z1 ) = 2.3560606
11
1
z2 = (35 − 6x1 − 3y1 ) = 0.9166666
2
60
Third approximation: Substituting x2 , y2 , z2 on RHS of Equation (74), we get
1
x3 = (20 + 3y2 − 2z2 ) = 3.1543561
8
1
y3 = (33 − 4x2 + z2 ) = 2.030303
11
1
z3 = (35 − 6x2 − 3y2 ) = 0.8797348
2
Proceeding in the same way, we get the following results
x4 = 3.0419299 x5 = 3.0168733 x6 = 3.0104409 x7 = 3.0157694
y4 = 1.9329373 y5 = 1.9696539 y6 = 1.9859295 y7 = 1.9885503
z4 = 0.8319128 z5 = 0.9127173 z6 = 0.9158165 z7 = 0.9149638
Until
x11 = 3.0167368 x12 = 3.0167424
y11 = 1.9858681 y12 = 1.9858987
z11 = 0.9118326 x12 = 0.9118312
Therefore, from 11th and 12th approximations, the values of x, y, z are same correct to four
decimal places. Stopping at this stage, we get
x = 3.0167;
y = 1.9858;
z = 0.9118.
Exercise: Solve the following system of equations by Jacobi iteration method:
3x + 4y + 15z = 54.8;
x + 12y + 3z = 39.66;
10x + y − 2z = 7.74.

4.3.2 GAUSS - SEIDEL Iteration Method

This is a modification of Gauss - Jacobi method. As before, the system of linear equations
a1 x + b 1 y + c 1 z = d 1
a2 x + b 2 y + c 2 z = d 2
a3 x + b 3 y + c 3 z = d 3 (75)
is written as
1
x = (d1 − b1 y − c1 z) (a)
a1
1
y = (d2 − a2 x − c2 z) (b)
b2
1
z = (d3 − a3 x − b3 y) (c) (76)
c3
and we start with the initial approximations x0 , y0 , z0 . Substituting y0 and z0 in equation (a),
we get
1
x1 = (d1 − b1 y0 − c1 z0 ).
a1

61
Now substituting x = x1 , z = z0 in Equation (b), we get
1
y1 = (d2 − a2 x1 − c2 z0 ).
b2
Now substituting x = x1 , y = y1 in Equation (c), we get
1
z1 = (d3 − a3 x1 − b3 y1 ).
c3
This process is continued till the values of x, y, z are obtained to the desired degree of accuracy.
The general algorithm is as follows:
1
xk+1 = (d1 − b1 yk − c1 zk )
a1
1
yk+1 = (d2 − a2 xk+1 − c2 zk )
b2
1
zk+1 = (d3 − a3 xk+1 − b3 yk+1 ) (77)
c3
Since the current values of the unknows at each stage of iteration are used in preceeding to the
next stage of iteration, this method is more rapid in convergence than Gauss- jacobi method.
Example 30. Solve by Gauss-Seidal iteration method the system
8x − 3y + 2z = 20;
4x + 11y − z = 33;
6x + 3y + 12z = 35.
So that the diagonal elements are dominant in the coefficient matrix. Now we write the
equations in the form
1
x = (20 + 3y − 2z) (a)
8
1
y = (33 − 4x + z) (b)
11
1
z = (35 − 6x − 3y) (c)
2
Putting y = 0, z = 0 in RHS of (a), we get
20
x1 = = 2.5.
8
Putting x = 2.5, z = 0 in RHS of (b), we get
1
y1 = (33 − 4(2.5)) = 2.0909091.
11
Putting x = 2.5, y = 2.0909091 in RHS of (c), we get
1
z1 = (35 − 6(2.5) − 3(2.0909091)) = 1.1439394.
12
For second approximation:
1
x2 = (20 + 3y1 − 2z1 ) = 2.9981061
8
1
y2 = (33 − 4x2 + z1 ) = 2.0137741
11
1
z2 = (35 − 6x2 − 3y2 ) = 0.9141701
2
62
For third approximation:
1
x3 = (20 + 3y2 − 2z2 ) = 3.0266228
8
1
y3 = (33 − 4x3 + z2 ) = 1.9825163
11
1
z3 = (35 − 6x3 − 3y3 ) = 0.9077262
2
For fourth approximation:
1
x4 = (20 + 3y3 − 2z3 ) = 3.0165121
8
1
y4 = (33 − 4x4 + z3 ) = 1.9856071
11
1
z4 = (35 − 6x4 − 3y4 ) = 0.9120088
2
For fifth approximation:
1
x5 = (20 + 3y4 − 2z4 ) = 3.0166005
8
1
y5 = (33 − 4x5 + z4 ) = 1.9859643
11
1
z5 = (35 − 6x5 − 3y5 ) = 0.9118753
2
For sixth approximation:
1
x6 = (20 + 3y5 − 2z5 ) = 3.0167678
8
1
y6 = (33 − 4x6 + z5 ) = 1.9858913
11
1
z6 = (35 − 6x6 − 3y6 ) = 0.9118099
2
For seventh approximation:
1
x7 = (20 + 3y6 − 2z6 ) = 3.0167568
8
1
y7 = (33 − 4x7 + z6 ) = 1.9858894
11
1
z7 = (35 − 6x7 − 3y7 ) = 0.9118159
2
Since at the sixth and seventh approximations, the values of x, y, z are the same, correct to four
decimal places; We can stop the iteration process. Therefore

x = 3.0167
y = 1.9858
z = 0.9118

We find that 12 iterations are necessary in Gauss - Jacobi method to get the same accuracy as
achieved by 7 iterations in Gauss - Seidel method.

63
 EXERCISES

1. Solve by Gauss- Seidel method the following system of equations

28x + 4y − z = 32
x + 3y + 10z = 24
2x + 17y + 4z = 35

2. Using Gauss-Seidel iteration method, find the solution of the following system of equations
correct to four decimal places.

10x − 2y − z − w = 3
−2x + 10y − z − w = 15
−x − y + 10z − 2w = 27
−x − y − 2z + 10w = −9

3. Solve the following system of linear equations by

1. Gauss-Jacobi method

2. Gauuss-seidel iteration method

(i)

2x + y + z = 4
x + 2y + z = 4
x + y + 2z = 4

(ii)

8x + y + z = 8
2x + 4y + z = 4
x + 3y + 5z = 5

64