NUMERICAL ANALYSIS I

LECTURE NOTE

PREPARED BY:
KEDIR ALIYI ([Link])

MAY, 2020
AMBO, ETHIOPIA
Contents
1 Basic Concepts in Error Estimation 4
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Theory of Errors or Error Analysis . . . . . . . . . . . . . . . . . . 4
1.3 Common types of errors: . . . . . . . . . . . . . . . . . . . . . . . 4
1.4 Basic sources of Errors . . . . . . . . . . . . . . . . . . . . . . . . 6
1.5 Approximations of Errors . . . . . . . . . . . . . . . . . . . . . . . 6
1.6 Propagation of Errors and Instability . . . . . . . . . . . . . . . . . 9

2 Solutions of non-linear Equations 11

2.1 Locating roots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 Bisection (Bolzano) method . . . . . . . . . . . . . . . . . . . . . 16
2.3 Method of False Position . . . . . . . . . . . . . . . . . . . . . . . 18
2.4 Secant Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.5 Newton-Raphson Method . . . . . . . . . . . . . . . . . . . . . . . 25
2.6 General Iteration Method . . . . . . . . . . . . . . . . . . . . . . . 29

3 System of equation 35
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.2 Direct Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.2.1 Gauss Elimination Method . . . . . . . . . . . . . . . . . 38
3.2.2 Gauss-Jordan Method . . . . . . . . . . . . . . . . . . . . 45
3.2.3 Matrix Decomposition . . . . . . . . . . . . . . . . . . . . 48
3.3 Iterative Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.3.1 Gauss-Jacobi Iteration Method . . . . . . . . . . . . . . . . 56
3.3.2 Gauss-Seidel Iteration Method . . . . . . . . . . . . . . . . 61
3.4 Solving System of non-linear Equations using Newton’s Method . . 66

4 FINITE DIFFERENCES 70
4.1 SHIFT OPERATOR . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.2 FORWARD DIFFERENCE OPERATOR . . . . . . . . . . . . . . 71
4.3 BACKWARD DIFFERENCES . . . . . . . . . . . . . . . . . . . . 75

2
5 Interpolation 77
5.1 Linear Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.2 Lagrange’s Interpolation formula . . . . . . . . . . . . . . . . . . . 78
5.3 Newton Interpolation Formula(forward and backward difference for-
mulas) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.3.1 NEWTON’S FORWARD INTERPOLATION FORMULA . 83
5.3.2 NEWTON– BACKWARD INTERPOLATION FORMULA . . . 87
5.4 Application of interpolation . . . . . . . . . . . . . . . . . . . . . . 90
5.4.1 Derivatives Using Newton’s Forward Difference Formula . . 90
5.4.2 Derivatives using Newton’s Backward Interpolation Formula 92
5.4.3 Numerical Integration . . . . . . . . . . . . . . . . . . . . 94

6 First Order Initial Value Problems 100

6.1 Definition of ODE and examples . . . . . . . . . . . . . . . . . . . 100
6.2 Order of a differential equation, linear and nonlinear ODE . . . . . 100
6.3 Nature of Solutions of ODE: particular and general solutions . . . . 102

3
Chapter 1

Basic Concepts in Error Estimation

1.1 Introduction
What is Numerical Analysis?
The aim of the field numerical analysis is to provide convenient methods for obtain-
ing useful solutions to mathematical problems. Such mathematical problems can be
formulated, for example, in terms of an algebraic or transcendental equation, an or-
dinary or partial differential equation, or an integral equation, or in terms of a set of
such equations.

Example 1.1.1 sinx + x2 ex + lnx = 0

The solution of this in general cannot be obtained so easily by known exact

methods. The numerical analysis will provide solutions which are not exact but
approximate. However, the approximate solutions should be close enough to exact
solution.

1.2 Theory of Errors or Error Analysis

The study of errors is a central concern of numerical analysis. Most numerical
methods give solutions that are only approximate to the desired solutions and it is
important to understand and able to estimate or found the resulting errors.

Error = TrueValue − ApproximateValue

1.3 Common types of errors:

Inherent errors:

• Errors which are already present in the statement of a problem before its use.

4
 • Arise either due to the given data being approximate or due to the limitations
of mathematical tables, calculators, or the digital computers.

• Can be minimized by taking better data or by using high precision computing

aids.

Rounding errors:

• Arise from the process of rounding off the numbers during the computation.

• Are unavoidable in most of the calculations due to the limitations of the com-
puting aids. However, rounding errors can be reduced.

Truncation errors:

• Are caused by using approximate results or on replacing an infinite process

by a finite one.

Example 1.3.1 If we use a decimal computer having a fixed word length of four
digitals, then

• Rounding off of 13.658 gives 13.66; and

• Truncation of 13.658 gives 13.65.

Hence,
Rounding error=13.66-13.658=0.002;and
Truncation error=13.658-13.65=0.008.

Absolute, Relative and Percentage Errors

If X is the true value of a quantity and X is its approximate value, then

• Absolute error, Ea = |X − X|;

Ea |X−X|
• Relative error, Er = |X| = |X| ; and

• Percentage error, E p = Er × 100.

5
If X be such a number that |X − X| ≤ X, then Xis an upper limit on the magnitude
of absolute error and measures the absolute accuracy.
Remarks:

1. The relative and percentage errors are independent of the units used.

1
2. If a number is correct to n decimal places, then the error(absolute) = 2 ×
10−n .

1.4 Basic sources of Errors

• Mathematical statements rarely give exact picture of actual phenomena. Most
of the cases are idealized models (Mathematical modeling).
This is the source of error (errors of the problem).

3 5
• Consider the function sinx = x − x3! + x5! · · ·
In this function, one has to consider a few terms to make computational prac-
tice. Hence, we are forced to stop at some term of the sequence and consider
it to be an approximation. Such a termination of the process gives rise to an
error. This error is called residual error or truncation error.

• Consider π = 22/7 = 3.141516 · · ·

We may have a non-terminating and a non-repeating decimal. It is obvious
that we can use only a finite number of digits in our computation. This is the
source of error called rounding error.

• When we are performing computational with approximate numbers we nat-

urally carry (to some extent) the errors of the original data into final result.
This is the source of error called inherent error.

1.5 Approximations of Errors

Errors in the Approximation of a Function Let y = f (x1 , x2 ) be a function of two
variables x1 and x2 . If δ x1 and δ x2 be the errors in x1 and x2 , respectively, then the

6
error δ y in y is given by

y + δ y = f (x1 + δ x1 , x2 + δ x2 ).

Expanding the right hand side by Taylor’s series, we get:

∂f ∂f
y + δ y = f (x1 , x2 ) + δ x1 + δ x2 +
∂ x1 ∂ x2
terms involving higher powers of

δ x1 andδ x2 · · · (i).

In general, the errorδ y in the function y = f (x1 , x2 , . . . , xn ) corresponding to the er-

ror δ xi in xi (i = 1, 2, . . . , n) is given by:

∂y ∂y ∂y
δy ≈ δ x1 + δ x2 + · · · + δ xn
∂ x1 ∂ x2 ∂ xn
and the relative error in y is given by:

δy ∂ y δ x1 ∂ y δ x2 ∂ y δ xn
Er = ≈ + +···+
y ∂ x1 y ∂ x2 y ∂ xn y

2 3
Example 1.5.1 : Let u = 4xz4y and errors in x, y, z be 0.001. Then compute the
maximum relative error in u when x = y = z = 1.

8xy3 ∂ u 12x2 y2 ∂ u 2 3
Solution: Since ∂u
∂x = ,
z4 ∂ y
= z4
, ∂z = − 16xz5 y , we have

∂u ∂u ∂u 8xy3 12x2 y2 16x2 y3

δu ≈ δx+ δy+ δz = 4 δx+ δ y + − 5 δ z.
∂x ∂y ∂z z z4 z

Since the errors δ x, δ y and δ z may be positive or negative, we take the absolute
values of the terms on the right hand side, giving

8xy3 12x2 y2 16x2 y3

δ umax ≈ | δ x| + | )δ y| + | − δ z|
z4 z4 z5

= 8(0.001) + 12(0.001) + 16(0.001) = 0.036

7
Hence, the maximum relative error in u is

Ermax = δ u/u ≈ 0.036/4 = 0.09.

Example 1.5.2 Find the maximum relative error in the function y = ax1m1 x2m2 · · · xnmn ,
where a, m1 , m2 , · · · , mn are constants.

Solution: lny = lna + m1 lnx1 + m2 lnx2 + · · · + mn lnxn

1 ∂y m1 1 ∂ y
Therefore, y ∂ x1 = x1 , y ∂ x2 = m2
x2 · · · , 1y ∂∂xyn = mn
xn .

∂ y δ x1
Hence, Er ≈ ∂ x1 y + ∂∂xy δ x2
y + · · · + ∂∂xyn δyxn .
2

δ x1 δ x2 δ xn
= m1 + m2 + · · · + mn .
x1 x2 xn
Since the errors δ x1 , δ x2 , · · · , δ xn may be positive or negative, we take the absolute
values of the terms on the right side. This gives:

δ x1 δ x2 δ xn
(Er )max = |m1 | + |m2 | + · · · + |mn |.
x1 x2 xn
Remark: Taking a = 1, m1 = m2 == mn = 1, we have y = x1 x2 · · · xn . Then

δ x1 δ x2 δ xn
Er = + +···+ .
x1 x2 xn

Thus the relative error of a product of n numbers is approximately equal to the

algebraic sum of their relative errors. Errors in a Series Approximation The Taylor’s
series for f(x) at x=a with a remainder after n terms is

(x − a)2 00 (x − a)n−1 (n−1)

f (x) = f (a+x−a) = f (a)+(x−a) f 0 (a)+ f (a)+· · ·+ f (a)+Rn (x)
2! (n − 1)!

n
where Rn (x) = (x−a) n
n! f (θ ), a < θ < x.
If the series is convergent, then Rn (x) → 0 as n → ∞ and hence if f (x) is approxi-
mated by the first n terms of this series, then the maximum error will be given by
the remainder term Rn (x). Moreover, if the accuracy required in a series approxi-

8
mation is pre-assigned, then we can find n, the number of terms which would yield
the desired accuracy.

Example 1.5.3 : Find the number of terms of the exponential series such that their
sum gives the value of ex correct to six decimal places at x = 1.

Solution: The Taylor’s series for ex at x = 0 is

x2 x3 xn−1
ex = 1 + x + + +···+ + Rn (x) · · · (i)
2! 3! (n − 1)!
n
where Rn (x) = xn eθ , 0 < θ < x. Therefore, the maximum absolute error (at θ =
n n
x) = xn! ex ; and the maximum relative error = xn!
Hence, (Er )maxatx=1 = n!1 .
For a six decimal place accuracy at x=1, we have

1 1
< × 10−6
n! 2

⇒ n! > 2 × 106
which gives n=10. Thus we need 10 terms of the series (i) in order that its sum is
correct to 6 decimal places.

1.6 Propagation of Errors and Instability

Algorithm:
An algorithm is a procedure that describes a finite sequence of steps to be performed
in specified order to solve a given problem.

Propagated error:

Propagated error is an error in the succeeding step of an algorithm due to an

error in the initial step.

• If the error E0 is introduced at initial stage of the algorithm and after ’n’
subsequent operation of the algorithm, En is the error resulted, then this is the

9
 propagated error.

Stability:

• If errors are magnified continuously as the algorithm continues, eventually

they will over shadow the true value and hence destroying the validity of the
algorithm. If this is so, the algorithm is instable.

• If errors made at initial step die out (decrease) as the algorithm continuous,
the algorithm is said to be stable.

10
Chapter 2

Solutions of non-linear Equations

2.1 Locating roots
Introduction: A problem of great importance in science and engineering is that of
determining the roots/zeros of an equation of the form

f (x) = 0

A polynomial equation of the form

f (x) = pn (x) = an xn + an−1 xn−1 + +a1 x + a0 = 0

is called an algebraic equation. An equation which contains exponential functions

logarithmic functions, trigonometric functions etc are called a transcendental equa-
tion. For example
3x3 − 2x2 − x − 5 = 0 ,x4 − 3x2 + 1 = 0, x2 − 3x + 1 = 0
Are algebraic (polynomial) equations, and xex − 1 = 0 , cosx − xe2x = 0 , tanx = x
are transcendental equations. Definition 1. A number α ,for which f (α) = 0 is
called a root of the equation f (x) = 0 is called a root of the equation f (x) = 0, or a
zero of f (x) = 0. Geometrically , a root of an equation f (x) = 0 is the value of x at
which the graph of the equation y = f (x) intersects the x-axis.

Figure 2.1: geometrical representation

2 A number α is a simple root of f (x) = 0 if f (α) = 0 and f 0 (α) 6= 0 , Then

we can write f (x) as f (x) = (x − α)g(x), g(α) 6= 0 For example , since x − 1 is

11
a factor of f (x) = x3 + x − 2 = 0 we can write f (x) = (x − 1)(x2 + x + 2) = (x −
1)g(x), g(1) = 4 6= 0 Alternatively, we find f (1) = 0,and f 0 (1) = 4 6= 0 .Hence
,x = 1 is a simple root of f (x) = x3 + x − 2 = 0 .
Multiple Root : A number α is a multiple root, of multiplicity m , of f (x) = 0
if f (α) = 0, f 0 (α) = 0 , f 00 (α) = 0 · · · f m−1 (α) = 0 and f m (α) 6= 0. Then we can
write f (x) as f (x) = (x − α)m g(x), g(α) 6= 0.
For example: consider the equation f (x) = x3 − 3x2 + 4 = 0 .
we find f (2) = 8 − 12 + 4 = 0 f 0 (x) = 3x2 − 6x, f 0 (2) = 12 − 12 = 0

f 00 (x) = 6x − 6, f 00 (2) = 60.

Hence x = 2 is a multiple root of multiplicity 2(double root) of

f (x) = x3 − 3x2 + 4 = 0
We can write f (x) = (x − 2)2 (x + 1) = (x − 2)2 g(x), g(2) = 30
In this chapter , we shall be considering the case of simple roots only.
Remark: 1. A polynomial equation of degree n has exactly n roots , real or complex
,simple or multiple ,where as a transcendental equation may have one root, infinite
number of roots or no root.
We shall derive methods for finding only the real roots. The methods for finding
the roots are classified as (i) direct methods ,and (ii) iterative methods
Direct Method s These methods give the exact value of all the roots in a finite num-
ber of steps( disregarding the round off errors).Therefore , for any direct method
, we can give the total number of operations ( addition, subtraction divisions and
multiplications ) This number is the operational count of the method.
For example ,the roots of the quadratic equation ax2 + bx + c = 0, a = 0 , can be
obtained using the method

1 p
x= [−b ± b2 − 4ac]
2a

For this method we can give the count of total number of operations. There are
direct methods for finding the roots of a cubic and fourth degree polynomial. How-
ever these methods are difficult to use. Direct methods for finding the roots of
polynomial equations of degree greater than 4 or transcendental equations are not
available in literature Iterative Methods: These methods are based on the idea of

12
successive approximatios. We start with one or two initial approximations to the
root and obtain a sequence of approximations x0 , x1 , x2 , · · · , xk , · · · which in the limit
as k → ∞ , converge to the exact root α. An iterative method for finding a root of
the equation f (x) = 0 can be obtained as

xk+1 = ϕ(xk ), k = 0, 1, 2, . . .

This method uses one initial approximation to the root x0 . The sequence of approx-
imations is given by

x1 = ϕ(x0 ), x2 = ϕ(x1 ), x3 = ϕ(x2 ), · · · .

The function ϕ is called the iteration function and x0 is called the initial an approxi-
mation. If the method uses two initial approximations x0 , x1 to the root , then we can
write the method as xk+1 = (xk−1 , xk ) k = 1, 2, · · · Convergence of iterative methods
: The sequence of iterates ,xk , is said to converge to the exact root α , if
Remark 2. Given one or two initial approximations to the root, we require a
suitable iteration function α for a given function f (x) , such that the sequence of
iterates,{xk } , converge to the exact root α. Further , we also require a suitable
criterion to terminate the iteration. Criterion to terminate iteration procedure; Since
we cannot perform infinite number of iterations, we need a criterion to stop the
iterations. We use one or both of the following criterion: The equation f (x) = 0 is
satisfied to a given accuracy or f (xk ) is bounded by an error tolerance ε. | f (xk )| ≤ ε.
The magnitude of the difference between two successive iterates is smaller than a
given accuracy or an error bound ε.

|xk+1 − xk | ≤ ε.

Generally, we use the second criterion. In some very special problems, we re-
quire to use both the criteria.
For example , if we require two decimal place accuracy , then we iterate until
|xk+1 − xk | ≤ 0.005.
If we require three decimal place accuracy, then we iterate until

13
|xk+1 − xk | ≤ 0.0005.
As we have seen earlier , we require a suitable iteration function and suitable
initial approximation(s) to start the iteration procedure. Initial approximation for an
iterative procedure
Descarte’s Rule of iterative procedure : This gives the bound for the number of
positive and negative real root for polynomial equations f (x) = 0;

i) We count the number of changes of signs is the coefficients of pn (x) for the
equation f (x) = pn (x) = 0 . The number of positive roots can not exceed the
number of changes of signs.
For example , If there are 4 changes of signs, then the equation may have 4
positive roots or two(2) positive roots or no positive roots.
If there are 3 changes of signs , then the equation may have 3 positive roots
or definitely one positive root.

ii) We write the equation f (−x) = pn (−x) = 0 and count the number of changes
of signs in the coefficients of pn (−x). The number of negative roots can not
exceed the number of changes of signs. For example , if there are 4 changes
of signs, then the equation may have 4 negative roots or 2. If there are 3
changes of signs , then the equation may have 3 negative roots or definitely
one negative root.

We use the following theorem of calculus to determine an initial approximation. It

is also called the intermediate value theorem.

Theorem 2.1.1 If f (x) is continuous on some interval [a, b] and f (a) f (b) < 0, then
the equation f (x) = 0 has at least one real root or an odd number of real roots in
the interval (a, b).

Example 2.1.1 Determine the maximum number of positive and negative roots and
intervals of length one unit in which the real roots lie for the following equations.

i. 8x3 − 12x2 − 2x + 3 = 0

ii. 3x3 − 2x2 − x − 5 = 0

14
Solutions: Let f (x) = 8x3 − 12x2 − 2x + 3 = 0 the number of changes in signs of
the coefficients (8, −12, −2, 3) is 2 . Therefore , the equation has 2 positive roots
or no positive roots.

i. Now f (−x) = −8x3 − 12x2 + 2x + 3 = 0 The number of changes in signs of

the coefficient (−8, −12, 2, 3) is one . Therefore the equation has one nega-
tive root. We have the following table of values for f (x)

x -2 -1 0 1 2 3
f (x) -105 -15 3 -3 15 105

Since f (−1) f (0) < 0 , there is a root in the interval (−1, 0)

f (0) f (1) < 0 there is a root in the interval (0, 1)
f (1) f (2) < 0 there is a root in the interval (1, 2)
Therefore there are 3 real roots (2 positive and 1 negative ) and the root lies
in the
intervals (−1, 0) , (0, 1) and (1, 2).

ii. Let f (x) = 3x3 − 2x2 − x − 5 = 0 ,The number of changes in signs of the
coefficients (3, −2, −1, −5) is 1 . Therefore the equation has 2 negative roots
or no negative root the table of values for f(x).

x -3 -2 0 1 2 3
f (x) -101 -35 -5 -5 9 55

From the table there is 1 positive real root in the interval (1, 2). The equation has
no negative real root .

Example 2.1.2 Locate the smallest positive real root of the equations

i. xex = cosx

ii. tanx = 2x

Solution

15
 i. Let f (x) = xex − cosx = 0 since f (0) = −1 and

f (1) = e − cos(1) = 2.178,

there is a root in the interval (0, 1).

ii. Let f (x) = tanx − 2x = 0 Then we have the following function values f (0) =
0 , f (0.1) = −0.0997 , f (0.5) = −0.4537 f (1) = −0.4426 f (1.1) = −0.2352
, f (1.2) = 0.1722, since f (1.1) f (1.2) < 0 ,the root lies in the interval (1.1, 1.2)

2.2 Bisection (Bolzano) method

This method is based on the repeated application of the intermediate value theorem.
Let the function f (x) be continuous on [a, b]. For definiteness let f (a) < 0 and
f (b) > 0 then the first approximation to the root.

a+b
x1 =
2
If f (x1 ) = 0 , then x1 is a root of f (x) = 0 .Otherwise, the root lies in (a, x1 ) or
(x1 , b) according as f (x1 ) is positive or negative. Then we bisect the interval as
before and continue the process until the root is found to desired accuracy. In the
figure below, f (x1 ) > 0, so that the root lies in the interval (a, x1 ) . Then the second
approximation to the root is
a + x1
x2 =
2
If f (x1 ) < 0) , then the root lies in (x2 , x1 ) .Then the third approximation to the root
is
x1 + x2
x3 =
2
and so on

Remark Since the new interval containing the root is exactly half the length
of the previous one , the interval width is reduced by a factor of 12 at each step.
Therefore at the end of the nth step , the interval will be of length

b−a
2n
16
 Figure 2.2: Bisection method

If on repeating this process n times , the latest interval is as small as given ε ,then
b−a ln( b−a
ε )
2 n ≤ ε or n ≥ ln2
This gives the number of iterations required for achieving an accuracy ε . In par-
ticular , the minimum number of iterations required for converging to a root in the
interval (0, 1) for a given ε are as under;

ε 10−2 10−3 10−4

n 7 10 14

Example 2.2.1 Find a root of the equation x3 − 4x − 9 = 0 using the bisection

method correct to 3 decimal places.

Solution: Let f (x) = x3 − 4x − 9 = 0. Since f (2) f (3) < 0, the root lies in the in-
terval (2, 3). Therefore, the first approximation to the root is x1 = 2+3 2 = 2.5 and
f (x1 ) = f (2.5) = −3.375.
Since f (2.5) f (3) < 0, the root lies in the interval (2.5, 3). Therefore, the 2nd ap-
proximation to the root is x2 = 2.5+3 2 = 2.75 and f (x2 ) = f (2.75) = 0.7969. Since
f (2.5) f (2.75) < 0, the root lies in the interval (2.5, 2.75). Therefore, the 3rd ap-
proximation to the root is x3 = 2.5+2.75
2 = 2.625 and f (x3 ) = f (2.625) = −1.4121.
Since f (2.625) f (2.75) < 0, the root lies in the interval (2.625, 2.75). Therefore,
the 4th approximation to the root is x4 = 2.625+2.75
2 = 2.6875.
Repeating this process, the successive approximations are
x5 = 2.71875 , x6 = 2.70313
x7 = 2.71094 , x8 = 2.70703
x9 = 2.70508 , x10 = 2.70605
x11 = 2.70654
Since |x11 − x10 | = 0.00012 < 0.0005, the required root is x ≈ x11 = 2.70654.

17
2.3 Method of False Position
The method is also called linear interpolation method or regula-falsi method. At the
start of all iterations of the method, we require the interval in which the root lies.
Let the root of the equation f (x) = 0 lie in the interval (xk−1 , xk ), that is, fk−1 fk < 0,
where f (xk−1 ) = fk−1 , and f (xk ) = fk . Then P(xk−1 , fk−1 ), Q(xk , fk ) are points on
the curve f (x) = 0. Draw a straight line joining the points P and Q From figure
(2.3). The line PQ is taken

Figure 2.3: Metod of false postion

as an approximation of the curve in the interval [xk−1 , xk ]. The equation of the line
fk − fk−1
PQ is given by y− fk
x−xk = xk −xk−1 . The point of intersection of this line PQ with the
x-axis is taken as the next approximation to the root. Setting y = 0, and solving for
x, we get
xk − xk−1
x = xk − ( ) fk .
fk − fk−1
The next approximation to the root is

xk − xk−1
xk+1 = xk − ( ) fk .
fk − fk−1

Simplifying, we can also write the approximation as

xk ( fk − fk−1 ) − (xk − xk−1 ) fk (xk−1 fk − xk fk−1 )

xk+1 = = , k = 1, 2, . . .
fk − fk−1 fk − fk−1

Therefore, starting with the initial interval (x0 , x1 ) in which the root lies, we com-
pute
x0 f1 − x1 f0
x2 = .
f1 − f0
Now, if f (x0 ) f (x2 ) < 0 , then the root lies in the interval (x0 , x2 ). Otherwise, the
root lies in the interval (x2 , x1 ). The iteration is continued using the interval in which

18
the root lies, until the required accuracy criterion is satisfied. Remark: At the start
of each iteration, the required root lies in an interval, whose length is decreasing.
Hence, the method always converges.
Remark: The method of false position has a disadvantage. If the root lies initially
in the interval (x0 , x1 ), then one of the end points is fixed for all iterations. For
example, in Fig.1.2a, the left end point x0 is fixed and the right end point moves
towards the required root. Therefore, in actual computations, the method behaves
like
x0 fk − xk f0
xk+1 = , k = 1, 2, . . .
fk − f0
In Fig.1.2b, the right end point x1 is fixed and the left end point moves towards the
required root. Therefore, in this case, in actual computations, the method behaves
like
x1 fk − xk f1
xk+1 = , k = 2, 3, . . .
fk − f1
Example 2.3.1 Locate the intervals which contain the positive real roots of the
equation x3 − 3x + 1 = 0.
Obtain these roots correct to three decimal places, using the method of false posi-
tion.

Solution: We form the following table of values for the function f(x).

x 0 1 2 3
f (x) 1 -1 3 19

There is one positive real root in the interval (0, 1) and another in the interval (1, 2).
There is no real root for x > 2 as f (x) > 0 for all x > 2. First, we find the root in
(0, 1). We have

x0 = 0, x1 = 1, f0 = f (x0 ) = f (0) = 1, f1 = f (x1 ) = −1

x0 f1 − x1 f0 0(−1) − 1
x2 = = = 0.5, f2 = f (x2 ) = f (0.5) = −0.375 < 0.
f1 − f0 −1 − 1
Since f (0) f (0.5) < 0 ,the root lies in the interval (0, 0.5)

x0 f2 − x2 f0 0 − 0.5
x3 = = = 0.36364.
f2 − f0 −0.375 − 1

19
and
f3 = f (x3 ) = f (0.36364) = −0.04283 < 0

Since f (0) f (0.36364) < 0, the root lies in (0, 0.36364)

x0 f3 − x3 f0 0 − 0.36364
x4 = = = 0.34870
f3 − f0 −0.04283 − 1
and f4 = f (x4 ) = f (0.34870) = −0.00370 < 0
Since f (0) f (0.34870) < 0 the root lies in (0, 0.34870).

x0 f4 − x4 f0 0 − 0.34870(1)
x5 = = = 0.34741
f4 − f0 −0.00370 − 0

and f5 = f (x5 ) = f (0.34741) = −0.00030 < 0

since f (0) f (0.34741) < 0 , the root lies in (0, 0.34741)

x0 f5 − x5 f0 0 − 0.34741
x6 = = = 0.347306
f5 − f0 −0.00030 − 1

Now |x6 − x5 | = |0.347306 − 0.34741| = 0.000104 < 0.0005

Since the root has been computed correct to 3 decimal places the required roots can
be taken as x ≈ x6 = 0.347306 , We may give the result as 0.347 , even though x6 is
more accurate .
Note: The left end point x = 0 is fixed for all iterations. Now we compute the root
in (1, 2).
we have x0 = 1, x1 = 2 , f0 = f (x0 ) = f (1) = −1 and f1 = f (x1 ) = f (2) = 3
Therefore
x0 f1 − x1 f0 3 − 2(−1)
x2 = = = 1.25
f1 − f0 3 − (−1)
f2 = f (x2 ) = f (1.25) = −0.796875
Since f (1.25) f (2) < 0, the root lies in (1.25, 2)

x2 f1 − x1 f2 1.25(3) − 2(−0.796875)
x3 = = = 1.407407
f1 − f2 3 − (−0.796875)

20
 f3 = f (x3 ) = f (1.407407) = −0.434437 Since f (1.407407) f (2) < 0, the root lies
in (1.407407, 2)

x3 f1 − x1 f3 1.407407(3) − 2(−0.434437)
x4 = = = 1.482367
f1 − f3 3 − (−0.434437)

f4 = f (x4 ) = f (1.482367) = −0.189730

Since f (1.482367) f (2) < 0 , the root lies in the interval (1.482367, 2)

x4 f1 − x1 f4 1.482367(3) − 2(−0.189730)
x5 = = = 1.513156
f1 − f4 3 − (−0.189730)

and
f5 = f (x5 ) = f (1.513156) = −0.074884

Since f (1.513156) f (2) < 0 , the root lies in the interval (1.513156, 2)

x5 f1 − x1 f5 1.513156(3) − 2(−0.074884)
x6 = = = 1.525012
f1 − f5 3 − (−0.074884)
and
f6 = f (x6 ) = f (1.525012) = −0.028374

Since f (1.525012) f (2) < 0 ,the root lies in the interval (1.525012, 5)

x6 f1 − x1 f6 1.525012(3) − 2(−0.028374)
x7 = = = 1.529462
f1 − f6 3 − (−0.028374)

f (x7 ) = f (1.529462) = −0.010586

Since, f (1.529462) f (2) < 0, the root lies in the interval (1.529462, 2).

x7 f1 − x1 f7 1.529462(3) − 2(−0.010586)
x8 = = = 1.531116
f1 − f7 3 − (−0.010586)

f (x8 ) = f (1.531116) = −0.003928

Since, f (1.531116) f (2) < 0, the root lies in the interval (1.531116, 2).

x8 f1 − x1 f8 1.531116(3) − 2(−0.003928)
x9 = = = 1.531729
f1 − f8 3 − (−0.003928)

21
 f (x9 ) = f (1.531729) = −0.001454.

Since, f (1.531729) f (2) < 0, the root lies in the interval (1.531729, 2).

x9 f1 − x1 f9 1.531729(3) − 2(−0.001454)
x10 = = = 1.531956
f1 − f9 3 − (−0.001454)

Now, |x10 − x9 | = |1.531956 − 1.53179| ≈ 0.000227 < 0.0005.

The root has been computed correct to three decimal places. The required root can
be taken as xx10 = 1.531956. Note that the right end point x = 2 is fixed for all
iterations.

Example 2.3.2 Find the root correct to two decimal places of the equation xex =
cosx, using the method of false position.

Solution Define f (x) = cosx − xex = 0. There is no negative root for the equation.
We have f (0) = 1, f (1) = cos1 − e = −2.17798.
A root of the equation lies in the interval (0, 1). Let x0 = 0, x1 = 1. Using the method
of false position, we obtain the following results.

x0 f1 − x1 f0 0 − 1(1)
x2 = = = 0.31467
f1 − f0 −2.17798 − 1

and
f (x2 ) = f (0.31467) = 0.51986

Since, f (0.31467) f (1) < 0, the root lies in the interval (0.31467, 1).

x2 f1 − x1 f2 0.31467(−2.17798) − 1(0.51986)
x3 = = = 0.44673
f1 − f2 −2.17798 − 0.51986

f (x3 ) = f (0.44673) = 0.20354.

Since, f (0.44673) f (1) < 0, the root lies in the interval (0.44673, 1).

x3 f1 − x1 f3 0.44673(2.17798) − 1(0.20354)
x4 = = = 0.49402
f1 − f3 −2.17798 − 0.20354

f (x4 ) = f (0.49402) = 0.07079.

22
Since, f (0.49402) f (1) < 0, the root lies in the interval (0.49402, 1).

x4 f1 − x1 f4 0.49402(−2.17798) − 1(0.07079)
x5 = = = 0.50995
f1 − f4 −2.17798 − 0.07079

f (x5 ) = f (0.50995) = 0.02360.

Since, f (0.50995) f (1) < 0, the root lies in the interval (0.50995, 1).

x5 f1 − x1 f5 0.50995(−2.17798) − 1(0.0236)
x6 = = = 0.51520
f1 − f5 −2.17798 − 0.0236

f (x6 ) = f (0.51520) = 0.00776.

Since, f (0.51520) f (1) < 0, the root lies in the interval (0.51520, 1).

x6 f1 − x1 f6 0.5152(−2.17798) − 1(0.00776)
x7 = = = 0.51692.
f1 − f6 −2.17798 − 0.00776

Now, |x7 − x6 | = |0.51692 − 0.51520| ≈ 0.00172 < 0.005.

The root has been computed correct to two decimal places. The required root can
be taken as x ≈ x7 = 0.51692.
Note that the right end point x = 2 is fixed for all iterations

Exercises 2.3.1 Find a root of the equation x3 −4x−9 = 0 using regula-falsi method
correct to 3 decimal places.

2.4 Secant Method

This method is an improvement over the method of false position as it does not
require the condition f (x0 ) f (x1 ) < 0 of the [Link] fig(2.4) Here also the graph
of the function y = f (x) is approximated by a secant line but at each iteration ,two
most recent approximations to the root are used to find the next [Link]
it is not necessary that the interval must contain the root. Taking x0 and x1 as the
intial limits of the interval the equation of the secant line through (x0 , f (x0 )) and
(x1 , f (x1 )) is
f (x1 ) − f (x0 )
y − f (x) = (x1 − x0 )
x1 − x0

23
 Figure 2.4: Secant method

Then the abscissa of the point where it crosses the x-axis (y=0) is given by

x1 − x0
x2 = x1 − ( ) f (x1 )
f (x1 ) − f (x0 )

which is an approximation to the [Link] ,the general formula for successive

approximations is given by

xk − xk−1
xk+1 = xk − ( ) f (xk ) , k = 1, 2, 3, . . .
f (xk ) − f (xk−1 )

Remark: If f (x) = f (xn−1 ) ,then this method fails and shows that it does not con-
verge necessarily. This is a drowback of secant method over the method of false
positon which always converges. But if the secant method once converges,its rate
of convergence is faster than that of the method of false position.

Example 2.4.1 Find the root of the equation cosx = xex using the secant method
correct to 4 decimal places .

Solutions: Let f (x) = cosx−xex = 0 Taking the initial approximations x0 = 0 , x1 =

1 so that f (x0 ) = 1, f (x1 ) = cos1 − e = −2.17798 Then by secant method we have

x1 − x0
x2 = x1 − ( ) f (x1 ) = 0.31467
f (x1 ) − f (x0 )

f (x2 ) = f (0.31467) = 0.51987

x2 − x1
x3 = x2 − ( ) f (x2 ) = 0.44673 , f (x3 ) = f (0.44673) = 0.20354
f (x2 ) − f (x1 )
x3 − x2
x4 = x3 − ( ) f (x3 ) = 0.53171
f (x3 ) − f (x2 )

24
Repeating the process the successive approximations are

x5 = 0.51690 , x6 = 0.51775 , x7 = 0.51776, etc

Hence the root is x ≈ 0.51776 correct to 4 decimal places

2.5 Newton-Raphson Method

This method is also called Newton’s method. This method is also a chord method
in which we approximate the curve near a root, by a straight line. Let x0 be an initial
approximation to the root of f (x) = 0. Then, P(x0 , f0 ), where f0 = f (x0 ), is a point
on the curve. Draw the tangent to the curve at P, (Fig. 1.3). We approximate the
curve in the neighborhood of the root by the tangent to the curve at the point P. The
point of intersection of the tangent with the x-axis is taken as the next approximation
to the root. The process is repeated until the required accuracy is obtained. The
equation of the tangent to the curve y = f (x) at the point P(x0 , f0 ) is given by

y − f (x0 ) = (x − x0 ) f 0 (x0 )

where f 0 (x0 ) is the slope of the tangent to the curve at p.

Figure 2.5: Newton Raphson Method

y = f (x) at the point P(x0 , f0 ) is given by

y − f (x0 ) = (x − x0 ) f (x0 )

where f (x0 ) is the slope of the tangent to the curve at P. Setting y = 0 and solving
for x, we get
f (x0 ) 0
x = x0 − 0 , f (x0 ) 6= 0
f (x0 )

25
 The next approximation to the root is given by

f (x0 )
x1 = x0 − , f 0 (x0 ) 6= 0
f 0 (x0 )

We repeat the procedure. The iteration method is defined as

f (xk )
xk+1 = xk − , f 0 (xk ) 6= 0
f 0 (xk )

This method is called the Newton-Raphson method or simply the Newton’s method.
The method is also called the tangent method. Alternate derivation of the method
Let xk be an approximation to the root of the equation f (x) = 0. Let ∆x be an
increment in x such that xk + ∆x is the exact root, that is f (xk + ∆x) = 0. Expanding
in Taylor’s series about the point xk , we get

(∆x)2 00
f (xk ) + ∆x f 0 (xk ) + f (xk ) + . . . = 0
2!

Neglecting the second and higher powers of ∆x, we obtain

f (xk ) + ∆x f 0 (xk ) ≈ 0,

or
f (xk )
∆x = −
f 0 (xk )
Hence, we obtain the iteration method

f (xk )
xk+1 = xk + ∆x = xk − , f 0 (xk ) 6= 0 , k = 0, 1, 2, . . .
f 0 (xk )
which is same as the method derived earlier.
Remark Convergence of the Newton’s method depends on the initial approxima-
tion to the root. If the approximation is far away from the exact root, the method
diverges (see Example 2.5.1). However, if a root lies in a small interval (a, b) and
x0 ∈ (a, b), then the method converges.
Remark : From Eq.(1.14), we observe that the method may fail when f 0 (x) is close
to zero in the neighborhood of the root. Later, in this section, we shall give the
condition for convergence of the method.

26
Remark : The computational cost of the method is one evaluation of the function
f (x) and one evaluation of the derivative f 0 (x), for each iteration.

Example 2.5.1 Derive the Newton’s method for finding N1 , where N > 0. Hence,
1
find 17 ,using the initial approximation as (i) 0.05, (ii) 0.15. Do the iterations con-
verge ?

Solution: Let x = N1 or 1x = N Define f (x) = 1x − N.

Then , f 0 (x) = − x12 .
Newton’s method gives

1
f (xk ) xk − N
xk+1 = xk – = xk − = xk + (xk − Nxk2 ) = 2xk − Nxk2
f 0 (xk ) 1
− x2
k

(i) With N = 17, and x0 = 0.05, we obtain the sequence of approximations

x1 = 2x0 − Nx02 = 2(0.05) − 17(0.05)2 = 0.0575

x2 = 2x1 − Nx12 = 2(0.0575) − 17(0.0575)2 = 0.058794.

x3 = 2x2 –Nx22 = 2(0.058794) − 17(0.058794)2 = 0.058823.

x4 = 2x3 − Nx32 = 2(0.058823) − 17(0.058823)2 = 0.058823.

Since, |x4 − x3| = 0, the iterations converge to the root.

The required root is 0.058823.
(ii) With N = 17, and x0 = 0.15, we obtain the sequence of approximations

x1 = 2x0 − Nx02 = 2(0.15) − 17(0.15)2 = −0.0825.

x2 = 2x1 − Nx12 = 2(−0.0825) − 17(−0.8025)2 = −0.280706.

x3 = 2x2 − Nx22 = 2(−0.280706) − 17(–0.280706)2 = −1.900942.

x4 = 2x3 − Nx32 = 2(−1.900942) − 17(−1.900942)2 = −65.23275.

We find that xk → −∞ as k increases.

Therefore, the iterations diverge very fast. This shows the importance of choosing
a proper initial approximation.

27
Example 2.5.2 Derive the Newton’s method for finding the qth root of a positive
number N, N 1/q , where N > 0, q > 0.
Hence, compute 171/3 correct to four decimal places, assuming the initial approxi-
mation as x0 = 2.

Solution : Let x = N 1/q , or xq = N. Define f (x) = xq − N. Then, f 0 (x) = qxq−1 .

Newton’s method gives the iteration
q q q q
xk − N qxk − (xk − N) (q − 1)xk + N
xk+1 = xk − q−1
= q−1
= q−1
qxk qxk qxk

1
For computing 17 3 , we have q = 3 and N = 17. Hence, the method becomes

(3 − 1)xk3 + 17 2xk3 + 17
xk+1 = = k = 0, 1, 2, 3, ...
3xk3−1 3xk2

With x0 = 2, we obtain the following results.

2x03 + 17 2(23 ) + 17
x1 = = = 2.75
3x02 3(2)2

2x13 + 17 2(2.75)3 + 17
x2 = = = 2.582645
3x12 3(2.75)2
2x23 + 17 2(2.582645)3 + 17
x3 = = = 2.571332
3x22 3(2.582645)2
2x33 + 17 2(2.571332)3 + 17
x4 = = = 2.571282.
3x32 3(2.571332)2
Now,
|x4 − x3 | = |2.571282 − 2.571332| = 0.00005.

We may take x ≈ 2.571282 as the required root correct to four decimal places

Example 2.5.3 Perform four iterations of the Newton’s method to find the smallest
positive root of the equation f (x) = x3 − 5x + 1 = 0.

Solution We have f (0) = 1, f (1) = –3. Since, f (0) f (1) < 0 , the smallest positive

28
root lies in the interval (0, 1). Applying the Newton’s method, we obtain

xk3 − 5xk + 1 2xk3 − 1

xk+1 = xk − = 2 , k = 0 , 1 , 2 , ...
3xk2 − 5 3xk − 5

Let x0 = 0.5. We have the following results.

2x03 − 1 2(0.5)3 − 1
x1 = 2 = = 0.176471,
3x0 − 5 3(0.5)2 − 5

2x13 − 1 2(0.176471)3 − 1
x2 = = = 0.201568,
3x12 − 5 3(0.176471)2 − 5
2x23 − 1 2(0.201568)3 − 1
x3 = = = 0.201640,
3x22 − 5 3(0.201568)2 − 5
2x33 − 1 2(0.201640)3 − 1
x4 = 2 = = 0.201640.
3x3 − 5 3(0.201640)2 − 5
Therefore, the root correct to six decimal places is x ≈ 0.201640.

Exercises 2.5.1 Using Newton-Raphson method solve xlog10 x = 12.34 with x0 =

10 correct to four decimal places.

ans:x ≈ 11.594854

2.6 General Iteration Method

The method is also called iteration method or method of successive approximations
or fixed point iteration method. The first step in this method is to rewrite the given
equation f (x) = 0 in an equivalent form as

x = ϕ(x). (1)

There are many ways of rewriting f (x) = 0 in this form. For example, f (x) =
x3 − 5x + 1 = 0, can be rewritten in the following forms.
r
x3 + 1 1 5x − 1
x= , x = (5x − 1) 3 , x = , etc. (2)
5 x

29
Now, finding a root of f (x) = 0 is same as finding a number α such that α = ϕ(α),
that is, a fixed point of ϕ(x). A fixed point of a function ϕ is a point α such that
α = ϕ(α). This result is also called the fixed point theorem. Using Eq.(1), the
iteration method is written as

xk+1 = ϕ(xk ), k = 0, 1, 2, ... (3)

The function ϕ(x) is called the iteration function. Starting with the initial approxi-
mation x0 , we compute the next approximations as

x1 = ϕ(x0 ), x2 = ϕ(x1 ), x3 = ϕ(x2 ), ...

The stopping criterion is same as used earlier. Since, there are many ways of writ-
ing f (x) = 0 as x = ϕ(x), it is important to know whether all or at least one of these
iteration methods converges.

Remark Convergence of an iteration method xk+1 = ϕ(xk ), k = 0, 1, 2, ..., de-

pends on the choice of the iteration function (x), and a suitable initial approximation
x0 , to the root. Consider again, the iteration methods given in Eq.(2), for finding a
root of the equation f (x) = x3 − 5x + 1 = 0. The positive root lies in the interval
(0, 1).
x3 + 1
(i) xk+1 = k k = 0, 1, 2, ... (4)
5

With x0 = 1, we get the sequence of approximations as

x1 = 0.4, x2 = 0.2128, x3 = 0.20193, x4 = 0.20165, x5 = 0.20164.

The method converges and x ≈ x5 = 0.20164 is taken as the required approximation

to the root.
1
(ii) xk+1 = (5xk − 1) 3 , k = 0, 1, 2, ... (5)

With x0 = 1, we get the sequence of approximations as

x1 = 1.5874, x2 = 1.9072, x3 = 2.0437, x4 = 2.0968, ...

30
which does not converge to the root in (0, 1).
s
5xk − 1
xk+1 = , k = 0, 1, 2, ... (6)
kk

With x0 = 1, we get the sequence of approximations as

x1 = 2.0, x2 = 2.1213, x3 = 2.1280, x4 = 2.1284, ...

which does not converge to the root in (0, 1). Now, we derive the condition that the
iteration function ϕ(x) should satisfy in order that the method converges.

Condition of convergence

Theorem 2.6.1 If

i. α be exact root of f (x) = 0

ii. I be any interval containing the points x = α

iii. |φ 0 (x)| < 1 for all x ∈ I

then the sequence of approximations x0 , x1 , x2 , . . . , xk , . . . will converge to the exact

root α provided that the initial approximation x0 is chosen in I.

proof
The iteration method for finding a root of f (x) = 0, is written as

xk+1 = ϕ(xk ) , k = 0, 1, 2, ... (2.1)

Let α be the exact root. That is,

α = ϕ(α) (2.2)

We define the error of approximation at the kth iterate as

εk = xk –,

31
k = 0, 1, 2,... Subtracting (2.2) from (2.1), we obtain

xk+1 − α = ϕ(xk ) − ϕ(α) = (xk − α)ϕ 0 (tk ) (2.3)

(using the mean value theorem) (2.3) or

εk+1 = ϕ(tk )εk , xk < tk < α.

Setting k = k − 1, we get

εk = ϕ 0 (tk−1 )εk−1 , xk−1 < tk−1 < α.

Hence,
εk+1 = ϕ 0 (tk )ϕ 0 (tk−1 )εk−1 .

Using (2.3) recursively, we get

εk+1 = ϕ 0 (tk )ϕ 0 (tk−1 ) . . . ϕ 0 (t0 )ε0 .

The initial error ε0 is known and is a constant. We have

|εk+1 | = |ϕ 0 (tk )||ϕ(tk–1 )|...|ϕ 0 (t0 )||ε0 |.

Let |ϕ 0 (tk )| ≤ c, k = 0, 1, 2, . . . Then,

|εk+1 | ≤ ck+1 |ε0 |. (2.4)

For convergence, we require that |εk+1 | → 0 as k → ∞. This result is possible, if

and only if c < 1. Therefore, the iteration method (2.1) converges, if and only if

|ϕ(xk )| ≤ c < 1, k = 0, 1, 2, . . .

or |ϕ 0 (x)| ≤ c < 1, for all x in the interval (a, b). (2.4) We can test this condition
using x0 , the initial approximation, before the computations are done.
Let us now check whether the methods (4), (5), (6) converge to a root in (0, 1) of
the equation f (x) = x3 − 5x + 1 = 0.
x3 +1 3x2 3x2
i. We have ϕ(x) = 5 , ϕ 0 (x) = 5 and |ϕ 0 (x)| = 5 ≤ 1, for all x in 0 < x <

32
 1. Hence, the method converges to a root in (0, 1).
1
ii. We have ϕ(x) = 5x − 1) 3 , ϕ 0 (x) = 5
3(5x−1)2/3

Now |ϕ 0 (x)| < 1, when x is close to 1 and |ϕ(x)| > 1 in the other part of the
interval. Convergence is not guaranteed.
q
iii. We have ϕ(x) = 5x−1 0 1
x , ϕ (x) = 3 1 .
2x 2 (5x−1) 2

Again, |ϕ 0 (x)| < 1, when x is close to 1 and |ϕ 0 (x)| > 1 in the other part of
the interval. Convergence is not guaranteed.

Example 2.6.1 Find the smallest positive root of the equation x3 −x−10 = 0, using
the general iteration method.

Solution We have

f (x) = x3 − x − 10, f (0) = −10, f (1) = −10,

f (2) = 8 − 2 − 10 = −4, f (3) = 27 − 3 − 10 = 14.

Since, f (2) f (3) < 0, the smallest positive root lies in the interval (2, 3). Write
1
x3 = x + 10, and x = (x + 10) 3 = ϕ(x). We define the iteration method as

1
xk+1 = (xk + 10) 3 .

We obtain
1
ϕ 0 (x) = 2
3(x + 10) 3
We find |ϕ 0 (x)| < 1 for all x in the interval (2, 3). Hence, the iteration converges.
Let x0 = 2.5. We obtain the following results.

1 1
x1 = (12.5) 3 = 2.3208, x2 = (12.3208) 3 = 2.3097,

1 1
x3 = (12.3097) 3 = 2.3090, x4 = (12.3090) 3 = 2.3089.

Since, |x4 − x3 | = |2.3089 − 2.3090| = 0.0001, we take the required root as x ≈

2.3089.

33
Example 2.6.2 Find the smallest negative root in magnitude of the equation 3x4 +
x3 + 12x + 4 = 0, using the method of successive approximations.

Solution We have

f (x) = 3x4 + x3 + 12x + 4 = 0, f (0) = 4, f (−1) = 3 − 1 − 12 + 4 = −6.

Since, f (−1) f (0) < 0, the smallest negative root in magnitude lies in the interval
(−1, 0). Write the given equation as

x(3x3 + x2 + 12) + 4 = 0,

and
4
x=− = ϕ(x)
3x3 + x2 + 12
The iteration method is written as

4
xk+1 = −
3xk3 + xk2 + 12

We obtain
0 4(9x2 + 2x)
ϕ (x) =
(3x3 + x2 + 12)2
We find |ϕ 0 (x)| < 1 for all x in the interval (−1, 0). Hence, the iteration converges.
Let x0 = −0.25. We obtain the following results.

4 4
x1 = − =− = −0.33290,
3x03 + x02 + 12 3(−0.25)3 + (−0.25)2 + 12

4 4
x2 = − =− = −0.33333,
3x13 + x12 + 12 3(−0.33290)3 + (−0.33290)2 + 12
4 4
x3 = − =− = −0.33333.
3x23 + x22 + 12 3(−0.33333)3 + (−0.33333)2 + 12
The required approximation to the root is x ≈ −0.33333.
We conclude that both fixed point iteration and regula-falsi methods converge slowly
as they have only linear rate of convergence. FUrther, Newton’s method converges
at least twice as fast as the fixed point iteration.

34
Chapter 3

System of equation
3.1 Introduction
Consider a system of n linear algebraic equations in n unknowns

a11 x1 + a12 x2 + . . . + a1n xn = b1

a21 x1 + a22 x2 + . . . + a2n xn = b2

..
.

an1 x1 + an2 x2 + . . . + ann xn = bn

where ai j , i = 1, 2, ..., n, j = 1, 2, . . . , n, are the known coefficients, bi , i = 1, 2, . . . , n,

are the known right hand side values and xi , i = 1, 2, . . . , n are the unknowns to be
determined.
In matrix notation we write the system as

Ax = b (3.1)
     
a11 a12 a13 . . . a1n x1 b1
a21 a22 a23 . . . a2n 
     
x2  b2 
Where A = 
 .. .. ..  , x =  ..  and b =  .. 
..     
 . . . .  . .
an1 an2 an3 . . . ann xn bn

The matrix [A|b], obtained by appending the column b to the matrix A is called
the augmented matrix. That is
 
a11 a12 . . . a1n b1
a21 a22 . . . a2n b2 
 

[A|b] =  .. .. .. .. .. 
. . . . . 
 

an1 an2 . . . ann bn
We define the following.

35
 i. The system of equations (3.1) is consistent (has at least one solution), if
rank(A) = rank[A|b] = r. If r = n, then the system has unique solution. If
r < n, then the system has (n − r) parameter family of infinite number of
solutions.

ii. The system of equations (3.1) is inconsistent (has no solution) if

rank(A) 6= rank[A|b].

We assume that the given system is consistent. The methods of solution of the
linear algebraic system of equations (3.1) may be classified as direct and iterative
methods.

(a) Direct methods produce the exact solution after a finite number of steps (dis-
regarding the round-off errors). In these methods, we can determine the total
number of operations (additions, subtractions, divisions and multiplications).
This number is called the operational count of the method.

(b) Iterative methods are based on the idea of successive approximations. We

start with an initial approximation to the solution vector x = x0 , and obtain a
sequence of approximate vectors x0 , x1 , ..., xk , ..., which in the limit as k → ∞
converge to the exact solution vector x.

Now, we derive some direct methods.

3.2 Direct Methods

If the system of equations has some special forms, then the solution is obtained
directly We consider two such special forms.

(a) Let A be a diagonal matrix, A = D. That is, we consider the system of equa-
tions Dx = b as

a11 x1 = b1
a22 x2 = b2
.. ..
. .
ann xn = bn

36
 This system is called a diagonal system of equations. Solving directly, we
obtain
bi
xi = , aii 6= 0, i = 1, 2, ..., n.
aii

(b) Let A be an upper triangular matrix, A = U. That is, we consider the system
of equations Ux = b as

a11 x1 + a12 x2 + · · · + a1n xn = b1

a22 x2 + · · · + a2n xn = b2
... ..
.
an−1,n−1 xn−1 + an−1,n xn = bn−1
ann xn = bn
This system is called an upper triangular system of equations. Solving for the un-
knowns in the order xn , xn−1 , ..., x1 , we get

bn
xn = ann
bn−1 −an−1,n xn
xn−1 = an−1 ,n−1
··· ···
(b1 −∑nj=2 a1 j x j )
x1 = a11

The unknowns are obtained by back substitution and this procedure is called the
back substitution method. Therefore, when the given system of equations is one of
the above two forms, the solution is obtained directly.

Elementary row transformations (operations)

The following operations on the rows of a matrix A are called the elementary row
transformations (operations).

(i) Interchange of any two rows. If we interchange the ith row with the jth row,
then we usually denote the operation as Ri ↔ R j .

(ii) Division/multiplication of any row by a non-zero number p. If the ith row is

multiplied by p, then we usually denote this operation as pRi → Ri .

(iii) Adding/subtracting a scalar multiple of any row to any other row. If all the
elements of the jth row are multiplied by a scalar p and added to the corre-

37
 sponding elements of the ith row, then, we usually denote this operation as
Ri + pR j → Ri .

These row operations change the form of A, but do not change the row-rank of A.
The matrix B obtained after the elementary row operations is said to be row equiv-
alent with A. In the context of the solution of the system of algebraic equations, the
solution of the new system is identical with the solution of the original system. The
above elementary operations performed on the columns of A (column C in place of
row R) are called elementary column transformations (operations). However, we
shall be using only the elementary row operations.

3.2.1 Gauss Elimination Method

The method is based on the idea of reducing the given system of equations Ax = b,
to an upper triangular system of equations Ux = z, using elementary row operations
This reduced system Ux = z, is then solved by the back substitution method to
obtain the solution vector x.
We illustrate the method using the 3 × 3 system

a11 x1 + a12 x2 + a13 x3 = b1

a21 x1 + a22 x2 + a23 x3 = b2 (3.2)
a31 x1 + a32 x2 + a33 x3 = b3

We write the augmented matrix [A — b] and reduce it to the following form

[A|b] Gauss elimination [U|z]

−−−−−−−−−−−−−−→

The augmented matrix of the system (3.2) is

 
a11 a12 a13 b1
 a21 a22 a23 b2  (3.3)
 

a31 a32 a33 b3

First stage of elimination We assume a11 6= 0. This element a11 in the 1 × 1 position
is called the first pivot. We use this pivot to reduce all the elements below this
pivot in the first column as zeros. Multiply the first row in (3.3) by aa2111
and aa11
31

respectively and subtract from the second and third rows. That is, we are performing

38
the elementary row operations R2 − aa21
11
R1 and R3 − aa31
11
R1 respectively. We obtain
the new augmented matrix as
 
a11 a12 a13 b1
(1) (1) (1)
 0 a22 a23 b2 (3.4)
 

(1) (1) (1)
0 a32 a33 b3

where

(1) a21 (1) a21 (1) a21

a22 = a22 − a12 , a23 = a23 − a13 , b2 = b2 − b1
a11 a11 a11

(1) a21 (1) a21 (1) a21

a32 = a32 − a12 , a33 = a33 − a13 , b3 = b3 − b1
a11 a11 a11
Second stage of elimination

(1) (1)
We assume a22 6= 0 . This element a22 in the 2 × 2 position is called the second
pivot. We use this pivot to reduce the element below this pivot in the second column
(1)
a32
as zero. Multiply the second row in (3.4) by (1) and subtract from the third row.
a22
(1)
a
That is, we are performing the elementary row operation R3 − ( a32
22
)R2 . We obtain
the new augmented matrix as
 
a11 a12 a13 b1
(1) (1) (1)
 0 a22 a23 b2 (3.5)
 

(2) (2)
0 0 a33 b3

Where
(1) (1)
(2) (1) a32 (1) (2) (1) a32 (1)
a33 = a33 − a , b3 = b3 −
(1) 23
b
(1) 2
a22 a22
(2)
The element a33 6= 0 is called the third pivot. This system is in the required upper
triangular form [U|z]. The solution vector x is now obtained by back substitution.
From the third row, we get x3 = b(2)
3
a33
(1) (2)
(b2 −a13 x3 )
From the second row, we get x2 = (2) .
a(12)
(b1 −a12 x2 −a13 x3 )
From the first row, we get x1 = a11 .

39
 In general, using a pivot, all the elements below that pivot in that column are
made zeros. Alternately, at each stage of elimination, we may also make the pivot
as 1, by dividing that particular row by the pivot.
Remark : When does the Gauss elimination method as described above fail? It
fails when any one of the pivots is zero or it is a very small number, as the elimi-
nation progresses. If a pivot is zero, then division by it gives over flow error, since
division by zero is not defined. If a pivot is a very small number, then division by it
introduces large round-off errors and the solution may contain large errors.

Pivoting Procedures How do we avoid computational errors in Gauss elimina-

tion? To avoid computational errors, we follow the procedure of partial pivoting. In
the first stage of elimination, the first column of the augmented matrix is searched
for the largest element in magnitude and brought as the first pivot by interchang-
ing the first row of the augmented matrix (first equation) with the row (equation)
having the largest element in magnitude. In the second stage of elimination, the
second column is searched for the largest element in magnitude among the n − 1
elements leaving the first element, and this element is brought as the second pivot
by interchanging the second row of the augmented matrix with the later row hav-
ing the largest element in magnitude. This procedure is continued until the upper
triangular system is obtained. Therefore, partial pivoting is done after every stage
of elimination. There is another procedure called complete pivoting. In this proce-
dure, we search the entire matrix A in the augmented matrix for the largest element
in magnitude and bring it as the first pivot.
This requires not only an interchange of the rows, but also an interchange of the
positions of the variables. It is possible that the position of a variable is changed a
number of times during this pivoting. We need to keep track of the positions of all
the variables. Hence, the procedure is computationally expensive and is not used in
any software.

Remark Gauss elimination method is a direct method. Therefore, it is possible

to count the total number of operations, that is, additions, subtractions, divisions
and multiplications. Without going into details, we mention that the total number
of divisions and multiplications (division and multiplication take the same amount

40
 2
of computer time) is n( n +3n−1
3 ). The total number of additions and subtractions
(addition and subtraction take the same amount of computer time) is n( (n−1)(2n+5)
6 ).

Remark When the system of algebraic equations is large, how do we conclude

that it is consistent or not, using the Gauss elimination method? A way of deter-
mining the consistency is from the form of the reduced system (3.5). We know that
if the system is inconsistent then rank(A) 6= rank[A|b]. By checking the elements
of the last rows, conclusion can be drawn about the consistency or inconsistency.
(2) (2)
Suppose that in (3.5), a33 6= 0 and b3 6= 0. Then, rank(A) = rank[A|b] = 3. The
system is consistent and has a unique solution. Suppose that we obtain the reduced
system as
 
a11 a12 a13 b1
(1) (1) (1) 
 0 a22 a23 b2 

(2)
0 0 0 b3
Then, rank(A) = 2, rank[A|b] = 3 and rank(A)rank[A|b]. Therefore, the system
is inconsistent and has no solution. Suppose that we obtain the reduced system as
 
a11 a12 a13 b1
(1) (1) (1) 
 0 a22 a23 b2 


0 0 0 0

Then, rank(A) = rank[A|b] = 2 < 3. Therefore, the system has 3 − 2 = 1 parameter

family of infinite number of solutions.

Example 3.2.1 Solve the system of equations

x1 + 10x2 –x3 = 3

2x1 + 3x2 + 20x3 = 7

10x1 − x2 + 2x3 = 4

using the Gauss elimination with partial pivoting.

41
Solution: We have the augmented matrix as
 
1 10 −1 3
 2 3 20 7 
 

10 −1 2 4

We perform the following elementary row transformations and do the eliminations.

 
10 1 2 4
R1 ↔ R3 :  2 3 20 7 
 

1 10 1 3

R2 − (R1 /5), R3 − (R1 /10) :

 
10 1 2 4
 0 3.2 19.6 6.2 
 

0 10.1 −1.2 2.6

R2 ↔ R3 :  
10 1 2 4
 0 10.1 −1.2 2.6 
 

0 3.2 19.6 6.2

R3 − (3.2/10.1)R2 :
 
10 1 2 4
 0 101 12 26 
 

0 0 19.98020 537624

Back substitution gives the solution.

5.37624
Third equation gives x3 = 19.98020 = 0.26908

Second equation gives

1 1
x2 = (2.6 + 1.2x3) = (2.6 + 1.2(0.26908)) = 0.28940
10.1 10.1

42
First equation gives

1 1
x1 = (4 + x2 − 2x3 ) = (4 + 0.2894 − 2(0.26908)) = 0.37512.
10 10

Example 3.2.2 Solve the system of equations

2x1 + x2 + x3 − 2x4 = −10

4x1 + 2x3 + x4 = 8

3x1 + 2x2 + 2x3 = 7

x1 + 3x2 + 2x3 − x4 = −5

using the Gauss elimination with partial pivoting.

Solution: The augmented matrix is given by

 
2 1 1 −2 −10
 

 4 0 2 1 8 

 3 2 2 0 7 
1 3 2 −1 −5

We perform the following elementary row transformations and do the elimina-

tions.  
4 0 2 1 8
 
 2 1 1 −2 −10 
R1 ↔ R2 : 
 3 2 2 0

 7 

1 3 2 −1 −5

 
4 0 2 1 8
1 0 − 52 −14 
 
 0
R2 − (1/2)R1 , R3 − (3/4)R1 , R4 − (1/4)R1 :  

 0 2 12 − 43 1 
3 5
0 3 2 − 4 −7

43
  
4
0 2 1 8
 3 5 
 3 0 2 − 4 −7 
R2 ↔ R4 :  

 2 21 − 34
0 1  
5
1 0 − 2 −14
0
 
4 0 2 1 8
 3 5 
 0 3 2 − 4 −7 
R3 − (2/3)R2 , R4 − (1/3)R2 : 
 0 0 −1 1 17


 2 12 3 
0 0 − 12 − 25 2 − 35
3
 
4 0 2 1 8
 3 5 
 0 3 2 − 4 −7 
R4 − R3 , :  0 0 −1 1 17 

 2 12 3 
0 0 0 − 6 − 52 13
3

Using back substitution, we obtain

52 6 17 1 17 1
x4 = (− )(− ) = 8 , x3 = −2( − )x4 = −2( − ( )8) = −10
3 13 3 12 3 12

1 3 5 1 3 5
x2 = [−7 − ( )x3 + ( )x4 ] = [−7 − ( )(−10) + ( )8] = 6
3 2 4 3 2 4
1 1
x1 = [8 − 2x3 − x4 ] = [8 − 2(−10) − 8] = 5.
4 4

Exercises 3.2.1 Solve the system of equations

3x1 + 3x2 + 4x3 = 20

2x1 + x2 + 3x3 = 13

x1 + x2 + 3x3 = 6

using the Gauss elimination method.

Example 3.2.3 Test the consistency of the following system of equations

x1 + 10x2 –x3 = 3

2x1 + 3x2 + 20x3 = 7

44
 9x1 + 22x2 + 79x3 = 45

using the Gauss elimination method.

Solution: We have the augmented matrix as

 
1 10 1 3
 2 3 20 7 
 

9 22 79 45

We perform the following elementary row transformations and do the eliminations.

 
1 10 1 3
R2 − 2R1 , R3 − 9R1 :  0 17 22 1 
 

0 68 88 18
.  
1 10 1 3
R3 − 4R2 :  0 17 22 1 
 

0 0 0 14
Now, rank [A] = 2, and rank[A|b] = 3. Therefore, the system is inconsistent and has
no solution.

3.2.2 Gauss-Jordan Method

The method is based on the idea of reducing the given system of equations Ax = b,
to a diagonal system of equations Ix = d, where I is the identity matrix, using ele-
mentary row operations. We know that the solutions of both the systems are iden-
tical. This reduced system gives the solution vector x. This reduction is equivalent
to finding the solution as x = A−1 b.

[A|b] Gauss-Jordan method [I|X]

−−−−−−−−−−−−−−−→

In this case, after the eliminations are completed, we obtain the augmented ma-
trix for a 3 × 3 system as  
1 0 0 d1
 0 1 0 d2  (3.6)
 

0 0 1 d3

45
and the solution is xi = di , i = 1, 2, 3.

Elimination procedure
The first step is same as in Gauss elimination method, that is, we make the elements
below the first pivot as zeros, using the elementary row transformations. From the
second step onwards, we make the elements below and above the pivots as zeros
using the elementary row transformations. Lastly, we divide each row by its pivot
so that the final augmented matrix is of the form (3.6). Partial pivoting can also
be used in the solution. We may also make the pivots as 1 before performing the
elimination.

Example 3.2.4 Solve the following system of equations

x1 + x2 + x3 = 1

4x1 + 3x2 − x3 = 6

3x1 + 5x2 + 3x3 = 4

using the Gauss-Jordan method (i) without partial pivoting, (ii) with partial pivot-
ing.

Solution: We have the augmented matrix as

 
1 1 1 1
 4 3 −1 6 
 

3 5 3 4

(i) We perform the following elementary row transformations and do the elimi-
nations.  
1 1 1 1
R2 − 4R1 , R3 − 3R1 :  0 1 5 2 
 

0 2 0 1
 
1 0 4 3
R1 + R2 , R3 + 2R2 :  0 1 5 2 
 

0 0 10 5

46
  
1 0 0 1
4 5
R1 − R3 , R2 − R3 :  0 −1 0 − 12 
 
10 10
0 0 −10 5
R3
Now, making the pivots as 1, ((−R2 ), −10 ) we get
 
1 0 0 1
 0 1 0 12 
 

0 0 1 − 21

Therefore, the solution of the system is x1 = 1, x2 = 12 , x3 = − 12 . (ii) We perform the

following elementary row transformations and do the elimination.
 
4 3 −1 6
R1 ↔ R2 :  1 1 1 1 
 

3 5 3 4

1 43 − 41 32
 
R1 
: 1 1 1 1 

4
3 5 3 4

1 34 − 14 32
 

R2 − R1 , R3 − 3R1 :  0 41 5
− 12 
 
4
0 114
15
4 − 12

1 43 − 14 3
 
2
4R2  15 2 
:  0 1 11 − 11
11

1 5 1
0 4 4 −2
14 18
 
1 0 − 11 11
3 1
R1 − R2 , R3 − R2 :  0 1 15 − 112 

4 4 11 
10 5
0 0 11 − 11

1 0 − 14 18
 
11 11
R3 /(10/11) :  0 1 15 − 2 

11 11 
0 0 1 − 12

47
  
1 0 0 1
R1 + (14/11)R3 , R2 − (15/11)R3 :  0 1 0 12 
 

0 0 1 − 21

Therefore, the solution of the system is x1 = 1, x2 = 21 , x3 = − 12 .

Inverse of a Matrix by Gauss-Jordan Method

We start with the augmented matrix of A with the identity matrix I of the same
order. When the Gauss-Jordan procedure is completed, we obtain

[A|I] Gauss-Jordan method [I|A−1 ]

−−−−−−−−−−−−−−−→

since, AA−1 = I.
Remark Partial pivoting can also be done using the augmented matrix [A—I].
However, we cannot first interchange the rows of A and then find the inverse. Then,
we would be finding the inverse of a different matrix.

Exercises 3.2.2 Find the inverse of the matrix

 
1 1 1
 4 3 −1 
 

3 5 3

using the Gauss-Jordan method (i) without partial pivoting, and (ii) with partial
pivoting.

3.2.3 Matrix Decomposition

LU decomposition
Suppose we have the system of equations

AX = B.

The motivation for an LU decomposition is based on the observation that systems

of equations involving triangular coefficient matrices are easier to deal with.

48
Definition: An LU decomposition of a matrix A is the product of a lower triangular
matrix and an upper triangular matrix that is equal to A.

It turns out that we need only consider lower triangular matrices L that have 1’s
down the diagonal. Here is an example, let
 
1 2 4
A =  3 8 14  = LU
 

2 6 13
   
1 0 0 U11 U12 U13
where L =  L21 1 0  and U =  0 U22 U23  Multiplying out LU and
   

L31 L32 1 0 0 U33

setting the answer equal to A gives
   
U11 U12 U13 1 2 4
LU =  L21U11 L21U12 +U22 L21U13 +U23  =  3 8 14 
   

L31U11 L31U12 + L32U22 L31U13 + L32U23 +U33 2 6 13

Now we have to use this to find the entries in L and U. Fortunately this is not nearly
as hard as it might at first seem. We begin by running along the top row to see that

U11 = 1 , U12 = 2 , U13 = 4.

Now consider the second row

L21U11 = 3 ⇒ L21 × 1 = 3 ∴ L21 = 3,

L21U12 +U22 = 8 ⇒ 3 × 2 +U22 = 8 ∴ U22 = 2,

L21U13 +U23 = 14 ⇒ 3 × 4 +U23 = 14 ∴ U23 = 2

Notice how, at each step, the equation in hand has only one unknown in it, and other
quantities that we have already found. This pattern continues on the last row

L31U11 = 2 ⇒ L31 × 1 = 2 ⇒ L31 = 2,

49
 L31U12 + L32U22 = 6 ⇒ 2 × 2 + L32 × 2 = 6 ⇒ L32 = 1,

L31U13 + L32U23 +U33 = 13 ⇒ (2 × 4) + (1 × 2) +U33 = 13 ⇒ U33 = 3.

We have shown that

    
1 2 4 1 0 0 1 2 4
A =  3 8 14  =  3 1 0   0 2 2 
    

2 6 13 2 1 1 0 0 3

and this is an LU decomposition of A.

Using an LU decomposition to solve systems of equations

Once a matrix A has been decomposed into lower and upper triangular parts it is
possible to obtain the solution to AX = B in a direct way. The procedure can be
summarised as follows

Figure 3.1: Lu decomposition

• Given A, find L and U so that A = LU. Hence LUX = B.

• Let D = UX so that LD = B. Solve this triangular system for Y .

• Finally solve the triangular system UX = D for X.

The benefit of this approach is that we only ever need to solve triangular systems.
The cost is that we have to solve two of them.

50
       
x1 1 2 4 x1 3
Example 3.2.5 Find the solution of X =  x2  of  3 8 14   x2  =  13 
      

x3 2 6 13 x3 4

Solution
• The first step is to calculate the LU decomposition of the coefficient matrix on the
lefthand side. In this case that job has already been done since this is the matrix we
considered earlier. We found that
   
1 0 0 1 2 4
L= 3 1 0  , U = 0 2 2 
   

2 1 1 0 0 3
 
d1
The next step is to solve LD = B for the vector D =  d2  That is we consider
 

d3
    
1 0 0 d1 3
LD =  3 1 0   d2  =  13  = B
    

2 1 1 d3 4

which can be solved by forward substitution. From the top equation we see that
d1 = 3. The middle equation states that 3d1 + d2 = 13 and hence d2 = 4. Finally
the bottom line says that 2d1 + d2 + d3 = 4 from which we see that d3 = 6. Now
that we have found D we finish the procedure by solving UX = D for X. That is we
solve     
1 2 4 x1 3
UX =  0 2 2   x2  =  4  = D
    

0 0 3 x3 6
by using back substitution. Starting with the bottom equation we see that 3x3 = −6
so clearly x3 = −2. The middle equation implies that 2x2 + 2x3 = 4 and it follows
that x2 = 4. The top equation states that x1 + 2x2 + 4x3 = 3 and consequently x1 = 3.

51
Therefore we have found that the solution to the system of simultaneous equations
      
1 2 4 x1 3 3
 3 8 14   x2  =  13  is X = 4 
      

2 6 13 x3 4 −2

QR Decomposition
If A is a m × n matrix with linearly independent columns, then A can be decom-
posed as A = QR, where Q is a m × n matrix whose columns form an orthonormal
basis for the column space of A and R is an nonsingular upper triangular matrix.

Theorem 3.2.1 : If A is a m × n matrix with linearly independent columns, then A

can be decomposed as A = QR , where Q is a m × n matrix whose columns form
an orthonormal basis for the column space of A and R is an nonsingular upper
triangular matrix.

Proof: Suppose A = [u1 | u2 | . . . | un ] and rank(A) = n. Apply the Gram-Schmidt

process to {u1 , u2 , ..., un } and the orthogonal vectors v1 , v2 , . . . , vn are

hui , v1 i hui , v2 i hui , vi−1 i

vi = ui − 2
v1 − 2
v2 − · · · − vi−1
kv1 k kv2 k kvi−1 k2
Let qi = kvvii k for i = 1, 2, . . . , n. Thus q1 , q2 , . . . , qn form a orthonormal basis for
the column space of A.
Now,

hui , v1 i hui , v2 i hui , vi−1 i

ui = vi + 2
v1 + 2
v2 + · · · + vi−1
kv1 k kv2 k kvi−1 k2
⇒ ui = kvi kqi + hui , q1 iq1 + hui , q2 iq2 + · · · + hui , qi−1 iqi−1

i.e. ui ∈ span{v1 , v2 , v3 , . . . , vi } = span{q1 , q2 , q3 , . . . , qi }

Thus ui is orthogonal to q j for j > i

u1 = kv1 kq1

u2 = kv2 kq2 + hu2 , q1 iq1

u3 = kv3 kq3 + hu3 , q1 iq1 + hu3 , q2 iq2

52
 ..
.

un = kvn kqn + hun , q1 iq1 + hun , q2 iq2 + · · · + hun , qn−1 iqn−1

Let Q = [q1 q2 ... qn ] , so Q is a m × n matrix whose columns form an orthonormal

basis for the column space of A . Now
 
kv1 k hu2 , q1 i hu3 , q1 i
. . . hun , q1 i
 

 0 kv2 k hu3 , q2 i
. . . hun , q2 i 

A = [u1 u2 . . . un ] = [q1 q2 q3 · · · qn ] 
 0 0 kv3 k
. . . hun , q3 i 

 .. .. .
.. .. .
..


 . . . 

0 0 . . . kvn k
0
(3.7)
Thus A can be decomposed as A = QR, where R is an upper triangular and nonsin-
gular matrix.

Example 3.2.6 Find the QR decomposition of

 
1 −1 −1
 
 1 0 0 
A= 

 1 −1 0 

0 0 −1

Solution: Applying Gram-Schmidt process of computing QR decomposition

Step 1
√
r11 = ka1 k = 3
 1 
√
3
√1
 
1 
3

q1 = a1 =  
ka1 k  √1

 3 
0
step 2
2
r12 = qT1 a2 = − √
3

53
step 3

√1
     
−1 3 − 13
√1
     2 
0 2
qb2 = a2 − q1 qT1 a2
     3 
= a2 − q1 r12 =   − (− √ )  3
= 
 −1  3  √1
  −1 
     3
3 
0 0 0

r
2
r22 = kqb2 k =
3
 
− √16
 q 
2 
1 
3 
q2 = qb2 = 

kqb2 k

1
 − √6 
 
0
step 4
1
r13 = qT1 a3 = − √
3
step 5
1
r23 = qT2 a3 = √
6
step 6

− 12
 

qb3 = a3 − q1 qT1 a3 − q2 qT2 a3 = a3 − r13 q1 − r23 q2 =  0 21 

 

−1
√
6
r33 = kqb3 k =
2
− √16
 
 
1  0 
q3 = qb3 = 
 √1

kqb3 
 6 
− √26

54
Therefore, A = QR
 
√1 − √16 − √16
 
1 −1 −1  √ 
3 √2 √1
   3 − 3 − 3

√1 √2 0

 1 0 0   
= 3 6
 0 √2 √1
  
 1 −1 0  
  √1 − √16 √1  6 √6 
 3 6 6
  0 0
0 0 −1 0 0 − √26 2

Uses: QR decomposition is widely used in computer codes to find the eigenvalues

of a matrix, to solve linear systems, and to find least squares approximations.

3.3 Iterative Methods

general linear iterative method for the solution of the system of equations Ax = b,
can be written in matrix form as

x(k+1) = Hx(k) + c, k = 0, 1, 2, . . .

respectively. H is called the iteration matrix, which depends on A and c is a column

vector, which depends on A and b. When to stop the iteration We stop the iteration
procedure when the magnitudes of the differences between the two successive iter-
ates of all the variables are smaller than a given accuracy or error tolerance or an
error bound ε, that is,
|xik+1 − xik | ≤ ε, for all i (3.8)

Now, we derive two iterative methods for the solution of the system of algebraic
equations
a11 x1 + a12 x2 + a13 x3 = b1
a21 x1 + a22 x2 + a23 x3 = b2 (3.9)
a31 x1 + a32 x2 + a33 x3 = b3

55
3.3.1 Gauss-Jacobi Iteration Method

Sometimes, the method is called Jacobi method. We assume that the pivots aii 6= 0,
for all i. Write the equations as

a11 x1 = b1 − (a12 x2 + a13 x3 )

a22 x2 = b2 − (a21 x1 + a23 x3 )
a33 x3 = b3 − (a31 x1 + a32 x2 )

The Jacobi iteration method is defined as

(k+1) 1 k k
x1 = a11 [b1 − (a12 x2 + a13 x3 )]
(k+1) 1 k k
a22 [b2 − (a21 x1 + a23 x3 )]
x2 = k = 0, 1, 2, ... (3.10)
(k+1) 1 k k
x3 = a33 [b3 − (a31 x1 + a32 x2 )]

Since, we replace the complete vector x(k) in the right hand side of (3.10) at the end
of each iteration, this method is also called the method of simultaneous displace-
ment.
Remark A sufficient condition for convergence of the Jacobi method is that the
system of equations is diagonally dominant, that is, the coefficient matrix A is di-
agonally dominant. We can verify that
n
|aii | ≥ ∑ |ai j |.
j=1,i6= j

If the system is not diagonally dominant, we may exchange the equations, if pos-
sible, such that the new system is diagonally dominant and convergence is guar-
anteed. This implies that convergence may be obtained even if the system is not
diagonally dominant.
Remark If the system is diagonally dominant, then the iteration converges for any
initial solution vector. If no suitable approximation is available, we can choose
x = 0, that is xi = 0 for all i. Then, the initial approximation becomes xi = abiii , for
all i.

56
Example 3.3.1 Solve the system of equations

4x1 + x2 + x3 = 2
x1 + 5x2 + 2x3 = −6
x1 + 2x2 + 3x3 = −4

using the Jacobi iteration method. Use the initial approximations as

(i) xi = 0, i = 1, 2, 3 , (ii) x1 = 0.5, x2 = −0.5, x3 = −0.5.
Perform five iterations in each case.

Solution: Note that the given system is diagonally dominant. Jacobi method gives
the iterations as
(k+1) (k) (k)
x1 = 0.25[2 − (x2 + x3 )]
(k+1) (k) (k)
x2 = 0.2[−6 − (x1 + 2x3 )]
(k+1) (k) (k)
x3 = 0.33333[−4 − (x1 + 2x2 )], k = 0, 1, ...
(0) (0) (0)
We have the following results. (i) x1 = 0, x2 = 0, x3 = 0 First iteration

(1) (0) (0)

x1 = 0.25[2 − (x2 + x3 )] = 0.5

(1) (0) (0)

x2 = 0.2[−6 − (x1 + 2x3 )] = –1.2
(1) (0) (0)
x3 = 0.33333[−4 − (x1 + 2x2 )] = –1.33333.

Second iteration

(2) (1) (1)

x1 = 0.25[2 − (x2 + x3 )] = 0.25[2 − (−1.2 − 1.33333)] = 1.13333

(2) (1) (1)

x2 = 0.2[−6 − (x1 + 2x3 )] = 0.2[−6 − (0.5 + 2(−1.33333))] = −0.76668,
(2) (1) (1)
x3 = 0.33333[−4 − (x1 + 2x2 )] = 0.33333[−4 − (0.5 + 2(−1.2))] = −0.7.

Third iteration

(3) (2) (2)

x1 = 0.25[2 − (x2 + x3 )] = 0.25[2 − (−0.76668 − 0.7)] = 0.86667,

(3) (2) (2)

x2 = 0.2[−6 − (x1 + 2x3 )] = 0.2[−6 − (1.13333 + 2(−0.7))] = −1.14667,
(3) (2) (2)
x3 = 0.33333[−4−(x1 +2x2 )] = 0.33333[−4−(1.13333+2(−0.76668))] = −1.19998.

57
Fourth iteration

(4) (3) (3)

x1 = 0.25[2 − (x2 + x3 )] = 0.25[2–(–1.14667–1.19999)] = 1.08666,

(4) (3) (3)

x2 = 0.2[−6 − (x1 + 2x3 )] = 0.2[−6 − (0.86667 + 2(−1.19998))] = −0.89334,
(4) (3) (3)
x3 = 0.33333[−4−(x1 +2x2 )] = 0.33333[−4−(0.86667+2(−1.14667))] = −0.85777.

Fifth iteration

(5) (4) (4)

x1 = 0.25[2 − (x2 + x3 )] = 0.25[2 − (−0.89334 − 0.85777)] = 0.93778,

(5) (4) (4)

x2 = 0.2[−6 − (x1 + 2x3 )] = 0.2[−6 − (1.08666 + 2(−0.85777))] = −1.07422,
(5) (4) (4)
x3 = 0.33333[−4−(x1 +2x2 )] = 0.33333[−4−(1.08666+2(−0.89334))] = −1.09998.

It is interesting to note that the iterations oscillate and converge to the exact solution

x1 = 1.0, x2 = –1, x3 = −1.0.

(0) (0) (0)

(ii) x1 = 0.5, x2 = −0.5, x3 = −0.5.

First iteration

(1) (0) (0)

x1 = 0.25[2 − (x2 + x3 )] = 0.25[2 − (−0.5 − 0.5)] = 0.75,

(1) (0) (0)

x2 = 0.2[−6 − (x1 + 2x3 )] = 0.2[−6 − (0.5 + 2(−0.5))] = −1.1,
(1) (0) (0)
x3 = 0.33333[−4 − (x1 + 2x2 )] = 0.33333[−4 − (0.5 + 2(−0.5))] = −1.16667.

Second iteration

(2) (1) (1)

x1 = 0.25[2 − (x2 + x3 )] = 0.25[2 − (−1.1 − 1.16667)] = 1.06667,

(2) (1) (1)

x2 = 0.2[−6 − (x1 + 2x3 )] = 0.2[−6 − (0.75 + 2(−1.16667))] = −0.88333,
(2) (1) (1)
x3 = 0.33333[−4−(x1 +2x2 )] = 0.33333[−4−(0.75+2(−1.1))] = −0.84999.

Third iteration

(3) (2) (2)

x1 = 0.25[2 − (x2 + x3 )] = 0.25[2 − (−0.88333 − 0.84999)] = 0.93333,

58
 (3) (2) (2)
x2 = 0.2[−6 − (x1 + 2x3 )] = 0.2[−6 − (1.06667 + 2(−0.84999))] = −1.07334,
(3) (2) (2)
x3 = 0.33333[−4−(x1 +2x2 )] = 0.33333[−4−(1.06667+2(−0.88333))] = −1.09999.

Fourth iteration

(3) (3)
x1 (4) = 0.25[2 − (x2 + x3 )] = 0.25[2 − (−1.07334 − 1.09999)] = 1.04333,

(4) (3) (3)

x2 = 0.2[−6 − (x1 + 2x3 )] = 0.2[−6 − (0.93333 + 2(−1.09999))] = −0.94667,
(4) (3) (3)
x3 = 0.33333[−4−(x1 +2x2 )] = 0.33333[−4−(0.93333+2(−1.07334))] = −0.92887.

Fifth iteration

(5) (4) (4)

x1 = 0.25[2 − (x2 + x3 )] = 0.25[2 − (−0.94667 − 0.92887)] = 0.96889,

(5) (4) (4)

x2 = 0.2[−6 − (x1 + 2x3 )] = 0.2[−6 − (1.04333 + 2(−0.92887))] = −1.03712,
(5) (4) (4)
x3 = 0.33333[−4−(x1 +2x2 )] = 0.33333[−4−(1.04333+2(−0.94667))] = −1.04999.

Example 3.3.2 Solve the system of equations

26x1 + 2x2 + 2x3 = 12.6

3x1 + 27x2 + x3 = −14.3
2x1 + 3x2 + 17x3 = 6.0

using the Jacobi iteration method. Obtain the result correct to three decimal places.

Solution: The given system of equations is strongly diagonally dominant. Hence,

we can expect faster convergence. Jacobi method gives the iterations as

(k+1) (k) (k)

x1 = [12.6 − (2x2 + 2x3 )]/26

(k+1) (k) (k)

x2 = [−14.3 − (3x1 + x3 )]/27
(k+1) (k) (k)
x3 = [6.0 − (2x1 + 3x2 )]/17 k = 0, 1, ...

Choose the initial approximation as

(0) (0) (0)

x1 = 0, x2 = 0, x3 = 0.

59
We obtain the following First iteration

(1) 1 (0) (0) 1

x1 = [12.6 − (2x2 + 2x3 )] = [12.6] = 0.48462
26 26

(1) 1 (0) (0) 1

x2 = [−14.3 − (3x1 + x3 )] = [−14.3] = −0.52963,
27 27
(1) 1 (0) (0) 1
x3 = [6.0 − (2x1 + 3x2 )] = [6.0] = 0.35294.
17 17
Second iteration

(2) 1 (1) (1) 1

x1 = [12.6 − (2x2 + 2x3 )] = [12.6 − 2(−0.52963 + 0.35294)] = 0.49821,
26 26

(2) 1 (1) (1) 1

x2 = [−14.3−(3x1 +x3 )] = [−14.3−(3(0.48462)+0.35294)] = −0.59655,
27 27
(2) 1 (1) (1) 1
x3 = [−6.0−(2x1 +3x2 )] = [6.0−(2(0.48462)+3(−0.52963))] = 0.38939
17 17
Third iteration

(3) 1 (2) (2) 1

x1 = [12.6 − (2x2 + 2x3 )] = [12.6 − 2(−0.59655 + 0.38939)] = 0.50006,
26 26

(3) 1 (2) (2) 1

x2 = [−14.3−(3x1 +x3 )] = [−14.3−(3(0.49821)+0.38939)] = −0.59941,
27 27
(3) 1 (2) (2) 1
x3 = [−6.0−(2x1 +3x2 )] = [6.0−(2(0.49821)+3(−0.59655))] = 0.39960
17 17

Fourth iteration

(4) 1 (3) (3) 1

x1 = [12.6 − (2x2 + 2x3 )] = [12.6 − 2(−0.59941 + 0.39960)] = 0.50000
26 26

(4) 1 (3) (3) 1

x2 = [−14.3−(3x1 +x3 )] = [−14.3−(3(0.50006)+0.39960)] = −0.59999
27 27
(4) 1 (3) (3) 1
x3 = [−6.0−(2x1 +3x2 )] = [6.0−(2(0.50006)+3(−0.59941))] = 0.39989
17 17
We find
(4) (3)
|x1 − x1 | = |0.5 − 0.50006| = 0.00006,
(4) (3)
|x2 − x2 | = | − 0.59999 + 0.59941| = 0.00058,

60
 (4) (3)
|x3 − x3 | = |0.39989 − 0.39960| = 0.00029.

Three decimal places of accuracy have not been obtained at this iteration. Fifth
iteration

(5) 1 (4) (4) 1

x1 = [12.6 − (2x2 + 2x3 )] = [12.6 − 2(−0.59999 + 0.39989)] = 0.50001,
26 26

(5) 1 (4) (4) 1

x2 = [−14.3−(3x1 +x3 )] = [−14.3−(3(0.50000)+0.39989)] = −0.60000
27 27
(5) 1 (4) (4) 1
x3 = [−6.0−(2x1 +3x2 )] = [6.0−(2(0.50000)+3(−0.59999))] = 0.40000
17 17
. We find
(5) (4)
|x1 − x1 | = |0.50001 − 0.5| = 0.00001
(5) (4)
|x2 − x2 | = | − 0.6 + 0.59999| = 0.00001,
(5) (4)
|x3 − x3 | = |0.4 − 0.39989| = 0.00011.

Since, all the errors in magnitude are less than 0.0005, the required solution is

x1 = 0.5, x2 = –0.6, x3 = 0.4.

3.3.2 Gauss-Seidel Iteration Method

we use the updated values of x1 , x2 , ..., xi−1 in computing the value of the variable
xi . We assume that the pivots aii 6= 0, for all i.
We write the equations as

a11 x1 = b1 − (a12 x2 + a13 x3 )

a22 x2 = b2 − (a21 x1 + a23 x3 )
a33 x3 = b3 − (a31 x1 + a32 x2 )

The Gauss-Seidel iteration method is defined as

(k+1) 1 (k) (k)
x1 =a11 [b1 − (a12 x2 + a13 x3 )]
(k+1) (k+1) (k)
x2 = a122 [b2 − (a21 x1 + a23 x3 )] k = 0, 1, 2, ... (3.11)
(k+1) (k+1) (k+1)
x3 = a133 [b3 − (a31 x1 + a32 x2 )]

This method is also called the method of successive displacement.

61
 If both the Gauss-Jacobi and Gauss-Seidel methods converge, then Gauss-Seidel
method converges at least two times faster than the Gauss-Jacobi method.

Example 3.3.3 Find the solution of the system of equations

45x1 + 2x2 + 3x3 = 58

−3x1 + 22x2 + 2x3 = 47

5x1 + x2 + 20x3 = 67

correct to three decimal places, using the Gauss-Seidel iteration method.

Solution: The given system of equations is strongly diagonally dominant. Hence,

we can expect fast convergence. Gauss-Seidel method gives the iteration

(k+1) 1 (k) (k)

x1 = (58 − 2x2 − 3x3 ),
45

(k+1) 1 (k+1) (k)

x2 =(47 + 3x1 − 2x3 ),
22
(k+1) 1 (k+1) (k+1)
x3 = (67 − 5x1 − x2 ).
20
(0) (0) (0)
Starting with x1 = 0, x2 = 0, x3 = 0, we get the following results.

First iteration

(1) 1 (0) (0) 1

x1 = (58 − 2x2 − 3x3 ) = (58) = 1.28889,
45 45

(1) 1 (1) (0) 1

x2 = (47 + 3x1 − 2x3 ) = (47 + 3(1.28889) − 2(0)) = 2.31212,
22 22
(1) 1 (1) (1) 1
x3 = (67–5x1 − x2 ) = (67 − 5(1.28889) − (2.31212)) = 2.91217
20 20
Second iteration

(2) 1 (1) (1) 1

x1 = (58 − 2x2 − 3x3 ) = (58 − 2(2.31212) − 3(2.91217)) = 0.99198,
45 45

(2) 1 (2) (1) 1

x2 = (47 + 3x1 − 2x3 ) = (47 + 3(0.99198) − 2(2.91217)) = 2.00689,
22 22
62
 (2) 1 (2) (2) 1
x3 = (67 − 5x1 − x2 ) = (67 − 5(0.99198) − (2.00689)) = 3.00166.
20 20
Third iteration

(3) 1 (2) (2) 1

x1 = (58 − 2x2 − 3x3 ) = (58 − 2(2.00689) − 3(3.00166) = 0.99958,
45 45

(3) 1 (3) (2) 1

x2 = (47 + 3x1 − 2x3 ) = (47 + 3(0.99958) − 2(3.00166)) = 1.99979,
22 22
(3) 1 (3) (3) 1
x3 = (67 − 5x1 − x2 ) = (67 − 5(0.99958) − (1.99979)) = 3.00012.
20 20
Fourth iteration

(4) 1 (3) (3) 1

x1 = (58 − 2x2 − 3x3 ) = (58 − 2(1.99979) − 3(3.00012)) = 1.00000,
45 45

(4) 1 (4) (3) 1

x2 = (47 + 3x1 − 2x3 ) = (47 + 3(1.00000) − 2(3.00012)) = 1.99999,
20 22
(4) 1 (4) (4) 1
x3 = (67 − 5x1 − x2 ) = (67 − 5(1.00000) − (1.99999)) = 3.00000.
20 20
We find
(4) (3)
|x1 − x1 | = |1.00000 − 0.99958| = 0.00042,
(4) (3)
|x2 − x2 | = |1.99999 − 1.99979| = 0.00020,
(4) (3)
|x3 − x3 | = |3.00000 − 3.00012| = 0.00012.

Since, all the errors in magnitude are less than 0.0005, the required solution is

x1 = 1.0, x2 = 1.99999, x3 = 3.0.

Rounding to three decimal places, we get x1 = 1.0, x2 = 2.0, x3 = 3.0.

Example 3.3.4 Computationally show that Gauss-Seidel method applied to the

system of equations
3x1 − 6x2 + 2x3 = 23

−4x1 + x2 − x3 = −8

x1 − 3x2 + 7x3 = 17

63
diverges. Take the initial approximations as x1 = 0.9, x2 = −3.1, x3 = 0.9. Inter-
change the first and second equations and solve the resulting system by the Gauss-
Seidel method. Again take the initial approximations as x1 = 0.9, x2 = –3.1, x3 =
0.9, and obtain the result correct to two decimal places. The exact solution is
x1 = 1.0, x2 = –3.0, x3 = 1.0.

Solution: Note that the system of equations is not diagonally dominant. Gauss-
Seidel method gives the iteration

(k+1) 1 (k) (k)

x1 = [23 + 6x2 − 2x3 )]
3
(k+1) (k+1) (k)
x2 = [−8 + 4x1 + x3 ]
(k+1) 1 (k+1) (k+1)
x3 = [17 − x1 + 3x2 ].
7
Starting with the initial approximations x1 = 0.9, x2 = −3.1, x3 = 0.9, we obtain
the following results.

First iteration

(1) 1 (0) (0) 1

x1 = [23 + 6x2 − 2x3 ] = [23 + 6(−3.1) − 2(0.9)] = 0.8667,
3 3
(1) (1) (0)
x2 = [−8 + 4x1 + x3 ] = [−8 + 4(0.8667) + 0.9] = −3.6332
(1) 1 (1) (1) 1
x3 = [17 − x1 + 3x2 ] = [17 − (0.8667) + 3(−3.6332)] = 0.7477
7 7
Second iteration

(2) 1 (1) (1) 1

x1 = [23 + 6x2 − 2x3 ] = [23 + 6(−3.6332) − 2(0.7477)] = −0.0982,
3 3
(2) (2) (1)
x2 = [−8 + 4x1 + x3 ] = [−8 + 4(−0.0982) + 0.7477] = −7.6451,
(2) 1 (2) (2) 1
x3 = [17 − x1 + 3x2 ] = [17 + 0.0982 + 3(−7.6451)] = −0.8339.
7 7
Third iteration

(3) 1 (2) (2) 1

x1 = [23 + 6x2 − 2x3 ] = [23 + 6(−7.6451) − 2(−0.8339)] = −7.0676,
3 3

64
 (3) (3) (2)
x2 = [−8 + 4x1 + x3 ] = [−8 + 4(−7.0676) − 0.8339] = −37.1043,
(3) 1 (3) (3) 1
x3 = [17 − x1 + 3x2 ] = [17 + 7.0676 + 3(−37.1043)] = −12.4636.
7 7
It can be observed that the iterations are diverging very fast.
Now, we exchange the first and second equations to obtain the system

−4x1 + x2 − x3 = −8

3x1 − 6x2 + 2x3 = 23

x1 − 3x2 + 7x3 = 17.

The system of equations is now diagonally dominant. Gauss-Seidel method gives

iteration
(k+1) 1 (k) (k)
x1 = [8 + x2 − x3 ]
4
(k+1) (k+1) (k)
x2 = −[23 − 3x1 − 2x3 ]
(k+1) 1 (k+1) (k+1)
x3 = [17 − x1 + 3x2 ]
7
Starting with the initial approximations x1 = 0.9, x2 = –3.1, x3 = 0.9, we obtain
the following results. First iteration

(1) 1 (0) (0) 1

x1 = [8 + x2 − x3 ] = [8 − 3.1 − 0.9] = 1.0,
4 4

(1) −1 (1) (0) −1

x2 = [23 − 3x1 –2x3 ] = [23 − 3(1.0) − 2(0.9)] = −3.0333,
6 6
(1) 1 (1) (1) 1
x3 = [17 − x1 + 3x2 ] = [17 − 1.0 + 3(−3.0333)] = 0.9857
7 7
Second iteration

(2) 1 (1) (1) 1

x1 = [8 + x2 − x3 ] = [8 − 3.0333 − 0.9857] = 0.9953,
4 4

(2) 1 (2) (1) 1

x2 = − [23 − 3x1 − 2x3 ] = − [23 − 3(0.9953) − 2(0.9857)] = −3.0071,
6 6
(2) 1 (2) (2) 1
x3 = [17 − x1 + 3x2 ] = [17 − 0.9953 + 3(−3.0071)] = 0.9976.
7 7

65
Third iteration

(3) 1 (2) (2) 1

x1 = [8 + x2 − x3 ] = [8 − 3.0071 − 0.9976] = 0.9988,
4 4

(3) 1 (3) (2) 1

x2 = − [23 − 3x1 − 2x3 ] = − [23 − 3(0.9988) − 2(0.9976)] = −3.0014,
6 6
(3) 1 (3) (3) 1
x3 = [17 − x1 + 3x2 ] = [17 − 0.9988 + 3(−3.0014)] = 0.9996.
7 7
Fourth iteration

(4) 1 (3) (3) 1

x1 = [8 + x2 − x3 ] = [8 − 3.0014 − 0.9996] = 0.9998,
4 4

(4) 1 (4) (3) 1

x2 = − [23 − 3x1 − 2x3 ] = − [23 − 3(0.9998) − 2(0.9996)] = −3.0002,
6 6
(4) 1 (4) (4) 1
x3 = [17 − x1 + 3x2 ] = [17 − 0.9998 + 3(−3.0002)] = 0.9999.
7 7
We find
(4) (3)
|x1 − x1 | = |0.9998 − 0.9988| = 0.0010,
(4) (3)
|x2 − x2 | = | − 3.0002 + 3.0014| = 0.0012,
(4) (3)
|x3 − x3 | = |0.9999 − 0.9996| = 0.0003.

Since, all the errors in magnitude are less than 0.005, the required solution is x1 =
0.9998, x2 = −3.0002, x3 = 0.9999. Rounding to two decimal places, we get
x1 = 1.0, x2 = −3.0, x3 = 1.0.

3.4 Solving System of non-linear Equations using New-

ton’s Method

Now we consider the solution of simultaneous non-linear equations by Newton’s

method.
Consider the system
f (x, y) = 0
(3.12)
f (x, y) = 0

66
involving two non-linear equations. Let (x0 , y0 ) be an initial approximation to the
root of the system, and (x0 + h, y0 + k) be the root of the system given by (3.12).
Then we must have
f (x0 + h, y0 + k) = 0
(3.13)
g(x0 + h, y0 + k) = 0
Let us assume that f and g are differentiable expanding (3.13) by Taylor’s series, we
obtain
f (x0 + h, y0 + k) = f0 + h ∂∂xf + k ∂∂yf + . . . = 0
∂g
0
∂g
0 (3.14)
g(x0 + h, y0 + k) = g0 + h ∂ x + k ∂ y + . . . = 0
0 0

neglecting, the second and higher order terms and retaining only the linear terms of
(3.14), we obtain
h ∂∂xf + k ∂∂yf = − f0
∂g
0
∂g
0 (3.15)
h ∂ x + k ∂ y = −g0
0 0

where f0 = f (x0 , y0 ), ∂∂xf = ( ∂∂ xf )x=x0 ; ∂∂yf = ( ∂∂ yf )y=y0 solving (3.14) for h and k,
0 0
the next approximation of the root is given by

x1 = x0 + h

and
y1 = y0 + k

The above process is repeated to desired degree of accuracy.

Example 3.4.1 Solve

x2 − y2 = 4

and
x2 + y2 = 16

By Newton–Raphson Method

Solution: To obtain the initial approximation we replace the first equation by its
asymptote y = x, which gives
√
2x2 = 16 ⇒ x = 2 2
√ √
let x0 = 2 2, y0 = 2 2, and (x0 , y0 ) be the initial approximation to the root of the

67
system.

We have f = x2 − y2 − 4 ⇒ f0 = −4 and g = x2 + y2 − 16 ⇒ g0 = 0 differenti-

ating partially, we obtain
∂f ∂f
= 2x , = −2y
∂x ∂y
∂g ∂g
= 2x , = 2y
∂x ∂y
∂f √ √
so that ∂ x0 = 2x0 = 4 2 , ∂∂yf = −4 2
0

∂g √ ∂g √
= 2x0 = 4 2 , = 2y0 = 4 2
∂ x0 ∂ y0

The system of linear equations can be written as

∂f ∂f √ √
h +k = − f0 ⇒ h(4 2) − k(4 2) = −(−4)
∂ x0 ∂ y0

⇒ h − k = 0.7072

and
∂g ∂g √ √
h +k = −g0 ⇒ h(4 2) + k(4 2) = 0
∂ x0 ∂ y0
⇒ h+k = 0

so that
h–k = 0.7072 (i)

h+k = 0 (ii)

solving we get h = 0.3536, k = –0.3536

The second approximation to the root is given by
√
x1 = x0 + h = 2 2 + 0.3536 = 3.1820
√
y1 = y0 + k = 2 2 − 0.3536 = 2.4748

The process can be repeated.

68
Example 3.4.2 Solve

f (x, y) = x2 + y − 20x + 40 = 0

g(x, y) = x + y2 − 20y + 20 = 0

Solution: Let x0 = 0, y0 = 0 be the initial approximation to the root

∂f ∂f
f = x2 + y − 20x + 40 ⇒ = 2x − 20 , =1
∂x ∂y

∂g ∂g
g = x + y2 − 20y + 20 ⇒ =1 , = 2y − 20
∂x ∂y
and f0 = 40, g0 = 20

So that
∂f ∂f
= −20, =1
∂ x0 ∂ y0
∂g ∂g
=1 , = −20
∂ x0 ∂ y0
the linear equations are

∂f ∂f
h +k = − f0 ⇒ 20h + k = −40
∂ x0 ∂ y0

∂g ∂g
h +k = −g0 ⇒ h − 20k = −20
∂ x0 ∂ y0
Solving, we get

h = 2.055, k = 1.103
The next approximation is given by

x1 = x0 + h = 2.055

y1 = y0 + k = 1.103.

69
Chapter 4

FINITE DIFFERENCES
Introduction

Numerical Analysis is a branch of mathematics which leads to approximate so-

lution by repeated application of four basic operations of Algebra. The knowledge
of finite differences is essential for the study of Numerical Analysis. In this section
we introduce few basic operators.

4.1 SHIFT OPERATOR

Let y = f (x) be function of x and x, x + h, x + 2h, x + 3h, . . . , etc., be the consecutive
values of x, then the operator E is defined as

E f (x) = f (x + h),

E is called shift operator. It is also called displacement operator.

Note: E is only a symbol but not an algebraic sum. E 2 f (x) means the operator E is
applied twice on f (x), i.e.,
E 2 f (x) = E[E f (x)]

= E f (x + h)

= f (x + 2h)

Similarly
E n f (x) = f (x + nh)

and
E −n f (x) = f (x − nh).

The operator E has the following properties:

1. E( f1 (x) + f2 (x) + . . . + fn (x)) = E f1 (x) + E f2 (x) + . . . + E fn (x)

70
 2. E(c f (x)) = cE f (x) (where c is constant)

3. E m (E n f (x)) = E n (E m f (x)) = E m+n f (x) where m, n are positive integers

4. If n is positive integer E n [E −n f (x)] = f (x)

Alternative notation If y0 , y1 , y2 , . . . , yn , . . . , etc are consecutive values of the func-

tion y = f (x) corresponding to equally spaced values x0 , x1 , x2 , . . . , xn , etc., of x then
in alternative notation
Ey0 = y1

Ey1 = y2

...

E 2 y0 = y2

...

and in general E n y0 = yn .

4.2 FORWARD DIFFERENCE OPERATOR

Let y = f (x) be any function given by the values y0 , y1 , y2 , . . . , yn , which it takes
for the equidistant values x0 , x1 , x2 , . . . , xn , of the independent variable x, then y1 −
y0 , y2 − y1 , . . . , yn − yn–1 are called the first differences of the function y. They are
denoted by y0 , y1 , . . . , etc. We have

∆y0 = y1 − y0

∆y1 = y2 − y1

...

∆yn−1 = yn − yn−1

The symbol ∆ is called the difference operator. The differences of the first differ-
ences denoted by ∆y and ∆2 y0 , ∆2 y1 , . . . , ∆2 yn are called second differences, where

∆2 y0 = ∆[∆y0 ]

71
 = ∆[y1 − y0 ]

= ∆y1 − ∆y0

= y2 − 2y1 + y0

Similarly
∆r yk = ∆r−1 yk+1 − ∆r−1 yk

Alternative Notation
Let the functions y = f (x) be given at equal spaces of the independent variable x,
say at x = a, a + h, a + 2h, . . . , etc., and the corresponding values of f (a), f (a +
h), f (a + 2h), . . . , etc. The independent variable x is often called the argument and
the corresponding value of the dependent variable is of the function at x = a, and is
denoted by ∆ f (a). Thus we have

∆ f (a) = f (a + h) − f (a),

writing the above definition we can write

∆ f (a + h) = f (a + h + h) − f (a + h) = f (a + 2h) − f (a + h)

Similarly

∆2 f (a) = ∆(∆ f (a)) = ∆( f (a + h) − f (a)) = f (a + 2h) − 2 f (a + h) + f (a)

Note: The operator ∆ is called forward difference operator and in general it is de-
fined as
∆ f (x) = f (x + h) − f (x),

where h is called the interval of differencing. Difference Table

It is a convenient method for displaying the successive differences of a function.

72
The following table is an example to show how the differences are formed.

x y ∆y ∆2 y ∆3 y ∆4 y ∆5 y
x0 y0
∆y0
x1 y1 ∆2 y0
∆y1 ∆3 y0
x2 y2 ∆2 y1 ∆4 y
∆y2 ∆3 y1 ∆5 y0
x3 y3 ∆2 y2 ∆4 y1
∆y3 ∆3 y1
x4 y4 ∆2 y3
∆y4
x5 y5

The above table is called a diagonal difference table. The first term in the table is
y0 . It is called the leading term. The differences ∆y0 , ∆2 y0 , ∆3 y0 , . . . ∆n y, are called
the leading differences. The differences ∆n yn with a fixed subscript are called for-
ward differences. In forming such a difference table care must be taken to maintain
correct sign.
Properties of the Operator

1. If c is a constant then ∆c = 0.

2. ∆ is distributive, i.e., ∆[ f (x) ± g(x)] = ∆ f (x) ± ∆g(x).

3. If c is a constant then ∆c f (x) = c∆ f (x).

4. If m and n are positive integers then ∆m ∆n f (x) = ∆m+n f (x).

5. ∆[ f1 (x) + f2 (x) + . . . + fn (x)] = ∆ f1 (x) + ∆ f2 (x) + . . . + ∆ fn (x).

6. ∆ f (x)g(x) = f (x)∆g(x) + g(x)∆ f (x).

f (x) g(x)∆ f (x)− f (x)∆g(x)
7. ∆[ g(x) ]= g(x)g(x+h) .

Example 4.2.1 Construct a difference table for y = f (x) = x3 + 2x + 1 for x =

1, 2, 3, 4, 5.

73
Solution:
x y ∆y ∆2 y ∆3 y
1 4
9
2 13 12
21 6
3 34 18
39 6
4 73 24
63
5 136

Theorem 4.2.1 : The nth forward differences of a polynomial of the nth degree
are constant when the values of independent variable are at equal intervals. i.e
∆ f (x) = n!hn where f (x) is nth degree polynomial.

Proof Exercise

Example 4.2.2 By constructing a forward difference table and taking the second
order forward differences as constant find the sixth term of the series 8, 12, 19, 29, 42, . . . .

Solution: Let K be the sixth term of the series. The difference table is

x y ∆ ∆2
1 8
4
1 12 3
7
3 19 3
10
4 29 3
13
5 42 k − 55
k−4
6 k

74
The second differences are constant.

⇒ K − 55 = 3

⇒ K = 58.

The sixth term of the series is 58.

4.3 BACKWARD DIFFERENCES

Let y = f (x) be a function given by the values y0 , y1 , . . . yn which it takes for the
equally spaced values x0 , x1 , . . . , xn of the independent variable x. Then y1 − y0 , y2 −
y1 , . . . , yn − yn−1 are called the first backward differences of y = f (x). They are
denoted by 5y0 , 5y1 , . . . , 5yn respectively. Thus we have

y1 − y0 = 5y1

y2 − y1 = 5y2
..
.

yn − yn−1 = 5yn

where 5 is called the backward difference operator.

x y 5 52 y 53 y 54 y
x0 y0
5y0
x1 y1 52 y0
5y1 53 y0
x2 y2 52 y1 54 y
5y2 53 y1
x3 y3 52 y2
5y3
x4 y4

75
Note: In the above table the differences 5n y with a fixed subscript i, lie along the
diagonal upward sloping.
Alternative notation: Let the function y = f (x) be given at equal spaces of the
independent variable x at x = a, a + h, a + 2h, . . . then we define 5 f (a) = f (a) −
f (a − h) where 5 is called the backward difference operator, h is called the interval
of differencing.
In general we can define

5 f (x) = f (x) − f (x − h).

We observe that
5 f (x + h) = f (x + h) − f (x) = ∆ f (x)

5 f (x + 2h) = f (x + 2h) − f (x + h) = ∆ f (x + h)

...

5 f (x + nh) = f (x + nh) − f (x + (n − 1)h) = ∆ f (x + (n − 1)h)

Similarly we get
52 f (x + 2h) = 5[5 f (x + 2h)]

= 5[∆ f (x + h) = ∆[∆ f (x)] = ∆2 f (x)

...

5n f (x + nh) = ∆n f (x).

Relation between E and 5 :

5 f (x) = f (x) − f (x − h) = f (x) − E −1 f (x)

⇒ 5 = 1 − E −1

or
E −1
5=
E

76
Chapter 5

Interpolation
Introduction:
In this chapter, we discuss the problem of approximating a given function by
polynomials. There are two main uses of these approximating polynomials. The
first use is to reconstruct the function f (x) when it is not given explicitly and only
values of f (x) and/ or its certain order derivatives are given at a set of distinct
points called nodes or tabular points. The second use is to perform the required
operations which were intended for f (x), like determination of roots, differentiation
and integration etc. can be carried out using the approximating polynomial P(x).
The approximating polynomial P(x) can be used to predict the value of f (x) at a
nontabular point. The deviation of P(x) from f (x), that is f (x) − P(x), is called
the error of approximation. Let f (x) be a continuous function defined on some
interval [a, b], and be prescribed at n + 1 distinct tabular points x0 , x1 , . . . , xn such
that a = x0 < x1 < x2 < ... < xn = b. The distinct tabular points x0 , x1 , ..., xn may
be non-equispaced or equispaced, that is xk+1 − xk = h, k = 0, 1, 2, . . . , n − 1. The
problem of polynomial approximation is to find a polynomial Pn (x), of degree ≤ n,
which fits the given data exactly, that is,

Pn (xi ) = f (xi ), i = 0, 1, 2, . . . , n. (5.1)

The polynomial Pn (x) is called the interpolating polynomial. The conditions given
in (5.1) are called the interpolating conditions

Remark 1 Through two distinct points, we can construct a unique polynomial

of degree 1 (straight line). Through three distinct points, we can construct a unique
polynomial of degree 2 (parabola) or a unique polynomial of degree1 (straight line).
That is, through three distinct points, we can construct a unique polynomial of de-
gree ≤ 2. In general, through n + 1 distinct points, we can construct a unique poly-
nomial of degree ≤ n. The interpolation polynomial fitting a given data is unique.
We may express it in various forms but are otherwise the same polynomial

77
5.1 Linear Interpolation
Given two points (x0 , y0 ) and (x1 , y1 ), the linear polynomial passing through the
two points is the equation of the line passing through the points. One way to write
its formula is
x1 − x x − x0
P1 (x) = y0 + y1
x1 − x0 x1 − x0

Example 5.1.1 For the data points (2, 3) and (5, 7) find P1 (x).

Solution: let x0 = 2, y0 = 3 , x1 = 5 , y1 = 7

x1 − x x − x0 2−x x−2 7
P1 (x) = y0 + y1 = 3( ) + 7( ) = (5 − x) + (x − 2)
x1 − x0 x1 − x0 5−2 5−2 3

In general, if y0 = f (x0 ) and y1 = f (x1 ) for some function f , then

−x
P1 (x) = f (x0 )( xx11−x0
) + f (x1 )( xx−x 0
1 −x0
) is a linear approximation of f (x) for all x ∈
[x0 , x1 ].

5.2 Lagrange’s Interpolation formula

Let the data

x x0 x1 x2 ... xn
f (x) f (x0 ) f (x1 ) f (x2 ) . . . f (xn )

be given at distinct unevenly spaced points or non-uniform points x0 , x1 , . . . , xn.

This data may also be given at evenly spaced points. For this data, we can fit a
unique polynomial of degree ≤ n. Since the interpolating polynomial must use all
the ordinates f (x0 ), f (x1 ), . . . , f (xn ), it can be written as a linear combination of
these ordinates. That is, we can write the polynomial as

Pn (x) = l0 (x) f (x0) + l1 (x) f (x1 ) + ... + ln (x) f (xn )

(5.2)
= l0 (x) f0 + l1 (x) f1 + ... + ln (x) fn

where f (xi ) = fi and li (x), i = 0, 1, 2, . . . , n are polynomials of degree n. This

polynomial fits the data given in (5.1) exactly.

78
At x = x0 , we get

f (x0 ) = Pn (x0 ) = l0 (x0 ) f (x0 ) + l1 (x0 ) f (x1 ) + . . . + ln (x0 ) f (xn ).

This equation is satisfied only when l0 (x0 ) = 1 and li (x0 ) = 0, i 6= 0. At a general

point x = xi , we get f (xi )Pn (xi ) = l0 (xi ) f (x0 ) + . . . + li (xi ) f (xi ) + . . . + ln (xi ) f (xn ).
This equation is satisfied only when li (xi ) = 1 and l j (xi ) = 0, i j. Therefore, li (x),
which are polynomials of degree n, satisfy the conditions
(
0 if i 6= j
li (x j ) = (5.3)
1 if i = j

Since, li (x) = 0 at x = x0 , x1 , . . . , xi−1 , xi+1 , . . . , xn , we know that

(x − x0 ), (x − x1 ), . . . , (x − xi−1 ), (x − xi+1 ), . . . , (x − xn )

are factors of li (x). The product of these factors is a polynomial of degree n. There-
fore, we can write

li (x) = C(x − x0 )(x − x1 )...(x − xi−1 )(x − xi+1 )...(x − xn )

where C is a constant.
Now, since li (xi ) = 1, we get

li (xi ) = 1 = C(xi − x0 )(xi − x1 )...(xi − xi−1 )(xi − xi+1 )...(xi − xn )

Hence
1
C=
(xi − x0 )(xi − x1 )...(xi − xi−1 )(xi − xi+1 )...(xi − xn )
Therefore,

(x − x0 )(x − x1 )...(x − xi−1 )(x − xi+1 )...(x − xn )

li (x) = (5.4)
(xi − x0 )(xi − x1 )...(xi − xi−1 )(xi − xi+1 )...(xi − xn )

Note that the denominator on the right hand side of li (x) is obtained by setting x = xi
in the numerator. The polynomial given in (5.2) where li (x) are defined by (5.4)
is called the Lagrange interpolating polynomial and li (x) are called the Lagrange

79
fundamental polynomials.
We can write the Lagrange fundamental polynomials li (x) in a simple notation.
Denote
w(x) = (x − x0 )(x − x1 ) . . . (x − xn )

which is the product of all factors. Differentiating w(x) with respect to x and sub-
stituting x = xi we get

w0 (xi ) = (xi − x0 )(xi − x1 ) . . . (xi − xi−1 )(xi − xi+1 ) . . . (xi − xn )

since all other terms vanish. Therefore, we can also write li (x) as

w(x)
li (x) = .
(x − xi )(w0 (x)

Linear interpolation

For n = 1, we have the data

x x0 x1
f (x) f (x0 ) f (x1 )

The Lagrange fundamental polynomials are given by

x − x1
l0 (x) = ,
x0 − x1

x − x0
l1 (x) =
x1 − x0
The Lagrange linear interpolation polynomial is given by

P1 (x) = l0 (x) f (x0 ) + l1 (x) f (x1 ).

Quadratic interpolation

For n = 2, we have the data

x x0 x1 x2
f (x) f (x0 ) f (x1 ) f (x2 )

80
The Lagrange fundamental polynomials are given by

(x − x1 )(x − x2 )
l0 (x) = ,
(x0 − x1 )(x0 − x2 )

(x − x0 )(x − x2 )
l1 (x) =
(x1 − x0 )(x1 − x2
(x − x0 )(x − x1 )
l2 (x) =
x2 − x0 )(x2 − x1 )
The Lagrange quadratic interpolation polynomial is given by

P1 (x) = l0 (x) f (x0 ) + l1 (x) f (x1 ) + l2 (x) f (x2 ).

Example 5.2.1 Use Lagrange’s formula, to find the quadratic polynomial that takes
the values

x 0 1 3
y 0 1 0

Solution Since f0 = 0 and f2 = 0, we need to compute l1 (x) only. We have

(x − x0 )(x − x2 ) x(x − 3) 1
l1 (x) = = = (3x − x2 )
(x1 − x0 )(x1 − x2 ) (1)(−2) 2

The Lagrange quadratic polynomial is given by

1 1
f (x) = l1 (x) f (x1 ) = (3x − x2 )(1) = (3x − x2 )
2 2

Example 5.2.2 Construct the Lagrange interpolation polynomial for the data

x −1 1 4 7
f (x) −2 0 63 342

Hence, interpolate at x = 5.

81
Solution The Lagrange fundamental polynomials are given by

(x − x1 )(x − x2 )(x − x3 ) (x − 1)(x − 4)(x − 7)

l0 (x) = =
(x0 − x1 )(x0 − x2 )(x0 − x3 ) (−1 − 1)(−1 − 4)(−1 − 7)

1 3
=− (x − 12x2 + 39x − 28).
80
(x − x0 )(x − x2 )(x − x3 ) (x + 1)(x − 4)(x − 7)
l1 (x) = =
(x1 − x0 )(x1 − x2 )(x1 − x3 ) (1 + 1)(1 − 4)(1 − 7)
1 3
= (x − 10x2 + 17x + 28).
36
(x − x0 )(x − x1 )(x − x3 ) (x + 1)(x − 1)(x − 7)
l2 (x) = =
(x2 − x0 )(x2 − x1 )(x2 − x3 ) (4 + 1)(4 − 1)(4 − 7)
1 3
=− (x − 7x2 − x + 7).
45
(x − x0 )(x − x1 )(x − x2 ) (x + 1)(x − 1)(x − 4)
l3 (x) = =
(x3 − x0 )(x3 − x1 )(x3 − x2 ) (7 + 1)(7 − 1)(7 − 4)
1 3
=(x − 4x2 − x + 4).
144
Note that we need not compute l1 (x) since f (x1 ) = 0. The Lagrange interpolation
polynomial is given by

P3 (x) = l0 (x) f (x0 ) + l1 (x) f (x1 ) + l2 (x) f (x2 ) + l3 (x) f (x3 )

1 3 1 1 3
=− (x −12x2 +39x−28)(−2)− (x3 −7x2 −x+7)(63)+ (x −4x2 −x+4)(342)
80 45 144
= x3 − 1.

Hence, f (5) = P3 (5) = 53 − 1 = 124

Exercises 5.2.1 Using Lagrange’s interpolation formula find a polynomial which

passes the points (0, –12), (1, 0), (3, 6), (4, 12). Hence find the value of y at x = 6.

82
5.3 Newton Interpolation Formula(forward and back-
ward difference formulas)

5.3.1 NEWTON’S FORWARD INTERPOLATION FORMULA

Let y = f (x) be a function which takes the values y0 , y1 , y2 , . . . , yn corresponding to

the (n + 1) values x0 , x1 , x2 , . . . , xn of the independent variable x. Let the values x be
equally spaced, i.e.,
xr = x0 + rh, r = 0, 1, 2, . . . ,

where h is the interval of differencing. Let φ (x) be a polynomial of the nth degree in
x taking the same values as y corresponding to x = x0 , x1 , . . . , xn , then, φ (x) repre-
sents the continuous function y = f (x) such that f (xr ) = φ (xr ) for r = 0, 1, 2, . . . , n
and at all other points f (x) = φ (x) + R(x) where R(x) is called the error term (Re-
mainder term) of the interpolation formula. Ignoring the error term let us assume

f (x) ≈ φ (x) ≈ a0 +a1 (x−x0 )+a1 (x−x0 )(x−x1 ) . . . an (x−x0 )(x−x1 ) . . . (x−xn−1 )
(5.5)
the constants a0 , a1 , a2 , . . . , an can be determine as follows. Putting x = x0 in
(5.5) we get
f (x0 ) ≈ φ (x0 ) = a0

⇒ y0 = a0

putting x = x1 in (5.5) we get

f (x1 ) ≈ φ (x1 ) = a0 + a1 (x1 − x0 ) = y0 + a1 h

∴ y1 = y0 + a1 h
y1 − y0 ∆y0
⇒ a1 = =
h h
Putting x = x2 in (5.5) we get

f (x2 ) ≈ φ (x2 ) = a0 + a1 (x2 − x0 ) + a2 (x2 − x0 )(x2 − x1 )

∆y0
∴ y2 = a0 + a1 (2h) + a2 (2h)(h) = y0 + (2h) + a2 (2h)(h)
h

83
 ⇒ y2 = y0 + 2(y1 − y0 ) + a2 (2h2 )
y2 − 2y1 + y0 ∆2 y0
⇒ a2 = =
2h2 2!h2
Similarly by putting x = x3 , x = x4 , . . . , x = xn in (5.5) we get

∆3 y0 ∆4 y0 ∆n y0
a3 = , a 4 = , . . . , an =
3!h3 4!h4 n!hn

putting the values of a0 , a1 , . . . , an in (5.5) we get

∆y0 ∆2 y0 ∆n y0
f (x) ≈ φ (x) = y0 + (x−x0 )+ (x−x0 )(x−x1 )+· · ·+ (x−x0 )(x−x1 ) . . . (x−xn−1 )
h 2!h2 n!hn
(5.6)

x−x0
Writing p = h , we get x − x0 = ph

x − x1 = x − x0 − (x1 − x0 ) = ph − h = (p − 1)h

x − x2 = x − x1 + x1 − x2 = (x − x1 ) − (x2 − x1 ) = (p − 1)h − h = (p − 2)h

Similarly
x − xn−1 = (p − n + 1)h

Equation (5.6) can be written as

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

φ (x) = y0 + ∆y0 + ∆ y0 + · · · + ∆ y0
1! 2! n!

The above formula is called Newton’s forward interpolation formula.

Note:

1. Newton forward interpolation formula is used to interpolate the values of y

near the beginning of a set of tabular values.

2. y0 may be taken as any point of the table, but the formula contains only those
values of y which come after the value chosen as y0 .

84
Example 5.3.1 Evaluate y = e2x for x = 0.05 using the following table

x 0.00 0.10 0.20 0.30 0.40

y = e2x 1.000 1.2214 1.4918 1.8221 2.255

Solution The difference table is

x y = e2x ∆y ∆2 y ∆3 y ∆4 y
0.000 1.0000
0.2214
0.10 1.2214 0.0490
0.2704 0.0109
0.20 1.4918 0.0599 0.0023
0.3303 0.0132
0.30 1.8221 0.0731
0.4034
0.4 2.2255

We have x0 = 0.00, x = 0.05, h = 0.1.

x − x0 0.5 − 0.00
∴ p= = = 0.5
h 0.1

Using Newton’s forward formula

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

f (x) = y0 + ∆y0 + ∆ y0 + ∆ y0 + ∆ y0
1! 2! 3! 4!

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

f (0.05) = 1.0000+0.5×0.2214+ (0.0490)+ (0.0109)
2 6
0.5(0.5 − 1)(0.5 − 2)(0.5 − 3)
+ (0.0023)
24
= 1.000 + 0.1107 − 0.006125 + 0.000681 − 0.000090 = 1.105166

∴ f (0.05) ≈ 1.052.

Exercises 5.3.1 In an examination the number of candidates who obtained marks

85
between certain limits were as follows:

Marks 30 − 40 40 − 50 50 − 60 60 − 70 70 − 80
No.o f Students 31 42 51 35 31

Find the number of candidates whose scores lie between 45 and 50.
Solution: Fist of all we construct a cumulative frequency table for the given data.

U pper limits o f the class intervals 40 50 60 70 80

Cumulative f requency 31 73 124 159 190
let x = U pper limits o f the class intervals(marks) and y = Cumulative f requency
The difference table is

x y ∆y ∆2 y ∆3 y ∆4 y
40 31
42
50 73 9
51 −25
60 124 −16 37
35 12
70 159 −4
31
80 190

we have x0 = 40, x = 45, h = 10

x − x0 45 − 40
p= = = 0.5
h 10
and
y0 = 73, ∆y0 = 42, ∆2 y0 = 9, ∆3 y0 = −25, ∆4 y0 = 37

From Newton’s forward interpolation formula

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

f (x) = y0 + ∆y0 + ∆ y0 + ∆ y0 + ∆ y0
1! 2! 3! 4!

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

= 73+0.5(42)+ (9)+ (−25)+ (37)
2! 3! 4!

86
 = 31 + 21 − 1.125 − 1.5625 − 1.4452 = 47.8673

= 48 (approximately)

The number of students who obtained marks less than 45 = 48, and the number of
students who scored marks between 45 and 50 = 73 - 48 = 25.

Exercises 5.3.2 A second degree polynomial passes through the points (1, -1), (2,
-1), (3, 1), (4, 5). Find the polynomial.

5.3.2 NEWTON– BACKWARD INTERPOLATION FORMULA

Newton’s forward interpolation formula cannot be used for interpolating a value

of y near the end of a table of values. For this purpose, we use another formula
known as Newton–backward interpolation formula. It can be derived as follows.
Let y = f (x) be a function which takes the values y0 , y1 , y2 , . . . , yn corresponding
to the values x0 , x1 , x2 , . . . , xn of the independent variable x. Let the values of x be
equally spaced with h as the interval of differencing, i.e.,
Let xr = x0 + rh, r = 0, 1, 2, . . . , n
Let φ (x) be a polynomial of the nth degree in x taking the same values as y
corresponding to x = x0 , x1 , . . . , xn , i.e., φ (x) represents y = f (x) such that f (xr ) =
φ (xr ), r = 0, 1, 2, . . . , we may write φ (x) as

f (x) ≈ φ (x) = a0 +a1 (x−xn )+a2 (x−xn )(x−xn−1 )+· · ·+an (x−xn )(x−xn−1 ) . . . (x−x1 )
(5.7)
Putting x = xn is (5.7) we get

f (xn ) ≈ φ (xn ) = a0 .

⇒ yn = a0 .

Putting x = xn−1 in (5.7) we get

f (xn−1 ≈ φ (x) = a0 + a1 (xn−1 − xn )

⇒ yn−1 = yn + a1 (−h)

⇒ a1 h = yn − yn−1

87
 yn − yn−1 5yn
⇒ a1 = =
h 1!h
Putting x = xn−2 , we get

f (xn−2) ≈ φ (xn−2 ) = a0 + a1 (xn−2 − xn ) + a2 (xn−2 − xn )(xn−2 − xn−1 )

yn − yn−1
⇒ yn−2 = yn + (−2h) + a2 (−2h)(−h)
h
⇒ yn−2 = yn − 2yn + 2yn−1 + 2(h)2 a2
yn − 2yn−1 + yn−2 52 yn
⇒ a2 = =
2h2 2!h2
similarly putting x = xn−3 , x = xn−4 , x = xn−5 , . . . we get

53 yn 54 yn 5n yn
a3 = , a 4 = , . . . , a n =
3!h3 4!h4 n!hn

substituting these values in (5.7)

5yn 52 yn 5n yn
f (x) ≈ φ (x) = yn + (x−xn )+ (x−xn )(x−xn−1 )+. . .+ (x−xn )(x−xn−1 ) . . . (x−x1 )
1!h 2!h2 n!hn
(5.8)
x−xn
writing p = h we get x − xn = ph

∴ x − xn−1 = x − xn + xn − xn−1 = ph + h = (p + 1)h

⇒ x − xn−2 = (p + 2)h, . . . , (x − x1 ) = (p + n − 1)h

∴ The equation (5.8) may be written as

5yn P(p + 1) 2 p(p + 1)(p + 2) 3

f (x) ≈ φ (x) = yn + p + 5 yn + 5 yn
1! 2! 3!

p(p + 1)(p + 2) . . . (p + n − 1) n
+···+ 5 yn
n!
The above formula is known as Newton’s backward interpolation formula.

Example 5.3.2 The following data gives the melting point of an alloy of lead and
zinc, where t is the temperature in degrees c and P is the percentage of lead in the

88
alloy. P is the percentage of lead in the alloy.

P 40 50 60 70 80 90
t 180 204 226 250 276 304

Find the melting point of the alloy containing 84 percent lead.

Solution The value of 84 is near the end of the table, therefore we use the Newton’s
backward interpolation formula. The difference table is

x y 5y 52 y 53 y 54 y 55 y
40 184
20
50 204 2
22 0
60 226 2 0
24 0 0
70 250 2 0
26 0
80 276 2
28
90 304

We have xn = 90, x = 84, h = 10 ,tn = yn = 304, 5 tn = 5yn = 28, 52 yn = 2

53 yn = 54 yn = 55 yn = 0

x − xn 84 − 90
p== = −0.6
h 10
From Newton’s backward formula

5tn P(p + 1) 2
f (84) = tn + p + 5 tn + . . .
1! 2!
(−0.6)(−0.6 + 1)
f (84) = 304 + (−0.6)(28) + (2)
2!
= 304 − 16.8 − 0.24 = 286.96.

89
Exercises 5.3.3 Calculate the value of f (7.5) using table below

x 1 2 3 4 5 6 7 8
f (x) 1 8 27 64 125 216 343 512

5.4 Application of interpolation

5.4.1 Derivatives Using Newton’s Forward Difference Formula

Consider the data (xi , f (xi )) given at equispaced points xi = x0 + ih, i = 0, 1, 2, . . . , n

where h is the step length. The Newton’s forward difference formula is given by

∆y0 ∆2 y0 ∆n y0
f (x) = y0 + (x−x0 )+ (x−x0 )(x−x1 )+· · ·+ (x−x0 )(x−x1 ) . . . (x−xn−1 )
h 2!h2 n!hn

Set x = x0 + ph. we get

f (x) = f (x0 + ph)

p p(p−1) 2 p(p−1)(p−2)...(p−n+1) n (5.9)
= y0 + 1! ∆y0 + 2! ∆ y0 + · · · + n! ∆ y0

Differentiating (5.9) with respect to x, we get

df d f dp 1 d f
= =
dx d p dx h d p

dp 1
Since dx = h

df 1 1 1 1
= [∆ f0 + (2p−1)∆2 f0 + (3p2 −6p+2)∆3 f0 + (4p3 −18p2 +22p−6)∆4 f0 +. . .]
dx h 2 6 24
(5.10)
0
At x = x0 , that is, at s = 0, we obtain the approximation to the derivative f (x) as

1 1 1 1
f 0 (x) = [∆ f0 − ∆2 f0 + ∆3 f0 − ∆4 f0 + . . .]
h 2 3 4

Differentiating (5.10) with respect to x, we get

d2 f
dx2
= h1 ddp ( dd pf ) ddxp
1 1
(5.11)
= h2
[∆2 f0 + (p − 1)∆3 f0 + 12 (6p2 − 18p + 11)∆4 f0 + · · · ]

90
putting p = 0 in (5.11) we have

d2y 1 11
( 2
)x=x0 = 2 [∆2 f0 − ∆3 f0 − ∆4 f0 + · · · ] (5.12)
dx h 12

and so on similarly we get

d3y 1 3
( 3
)x=x0 = 3 [∆3 f0 + ∆4 f0 + · · · ] (5.13)
dx h 2

Example 5.4.1 Find dy/dx at x = 1 from the following table of values

x 1 2 3 4
f (x) 1 8 27 64

Solution We have the following forward difference table.

Forward difference table

x y ∆y ∆2 y ∆3 y
1 1
7
2 8 12
9 6
3 27 18
37
4 64
We have h = 1, x0 = 1, and x = x0 + sh = 1 + s. For x = 1, we get s = 0. Therefore,

dy 1 1 1
(1) = (∆ f0 − ∆2 f0 + ∆3 f0 )
dx h 2 3

1 1
= 7 − (12) + (6) = 3
2 3

Exercises 5.4.1 The following data gives the velocity of a particle for 8 seconds at
an interval of 2 seconds. Find the initial acceleration using the entire data.

time(sec) 0 2 4 6 8
velocity(m/sec) 0 172 1304 4356 10288

91
5.4.2 Derivatives using Newton’s Backward Interpolation For-
mula

Consider the Newton’s backward interpolation formula

y = yn + p 5y n
1! +
P(p+1)
2! 52 yn + p(p+1)(p+2)
3! 53 yn
(5.14)
+ · · · + p(p+1)(p+2)...(p+n−1)
n! 5n yn

where p = x−x h
n
⇒ ddxp = 1h h being the interval of differencing
Differentiating (5.14) w.r.t. x we get

df d f dp 1 d f
= =
dx d p dx h d p

df
dx = 1h [5 fn + 12 (2p + 1) 52 fn + 16 (3p2 + 6p + 2) 53 fn + 24
1
(4p3 + 18p2 + 22p + 6) 54 fn

1
+ 120 (5p4 + 40p3 + 105p2 + 100p + 24) 54 fn + · · · ]
(5.15)
At x = xn , we get s = 0. Hence, we obtain the approximation to the first derivative
f (xn ) as

1 1 1 1 1
f 0 (xn ) = [5 fn + 52 fn + 53 fn + 54 fn + 55 fn + · · · ]
h 2 3 4 5

Differentiating (5.15) with respect to x again, we get

d2 f
   
1 d df dp 1 d df
2
= = 2
dx h dp dp dx h d p dp

d2 f 1 1
dx2
= h2
[52 fn + (p + 1) 53 fn + 12 (6p2 + 18p + 11) 54 fn
(5.16)
1
+ 24 (4p3 + 24p2 + 42p + 20) 54 fn + · · · ]

92
At x = xn , that is, at s = 0, we obtain the approximation to the second derivative
f 00 (x) as
f 00 (xn ) = h12 [52 fn + 53 fn + 11 4
12 5 f n
(5.17)
+ 56 54 fn + 137 5
180 5 f n + · · · ]

Remark We use the backward difference formulas for derivatives, when we need
the values of the derivatives near the end of table of values.

Example 5.4.2 Find f 0 (3) using the Newton’s backward difference formula, for the
data

x 1.0 1.5 2.0 2.5 3.0

f (x) −1.5 −2.875 −3.5 −2.625 0.5

Solution: The step length is h = 0.5 and x = xn + sh = 3.0 + s(0.5). For x = 3.0,
we get s = 0. We have the following backward difference table.
Backward difference table

x y 5y 52 y 53 y 54 y
1.0 −1.5
−1.375
1.5 −2.875 0.75
−0.625 0.75
2.0 −3.5 1.5 0
0.875 0.75
2.5 −2.625 2.25
3.125
3.0 0.5

From the formula

1 1 1 1 1
f 0 (xn ) = [5 fn + 52 fn + 53 fn + 54 fn + 55 fn + · · · ]
h 2 3 4 5

we obtain
1 1 1
f 0 (3.0) = [3.125 + (2.25) + (0.75)] = 9
0.5 2 3

93
5.4.3 Numerical Integration

Introduction
The problem of numerical integration is to find an approximate value of the integral
Z b
I= w(x) f (x)dx
a

where w(x) > 0 in (a, b) is called the weight function. The function f (x) may be
given explicitly or as a tabulated data. We assume that w(x) and f (x) are integrable
on [a, b]. The limits of integration may be finite, semi-infinite or infinite. The
integral is approximated by a linear combination of the values of f (x) at the tabular
points as
Z b n
I= w(x) f (x)dx = ∑ λk f (xk )
a k=0

Integration Rules Based on Uniform Mesh Spacing

When w(x) = 1 and the nodes xk0 s are prescribed and are equispaced with x0 = a,
xn = b, where h = (b−a)/n, the methods are called Newton-Cotes integration rules.
The weights λk0 s are called Cotes numbers. We shall now derive some Newton-
Cotes formulas. That is, we derive formulas of the form
Z b n
I= f (x)dx = ∑ λk f (xk ) (5.18)
a k=0

We note that, ab f (x)dx defines the area under the curve y = f (x), above the x-axis,
R

between the lines x = a, x = b.

Trapezium Rule
This rule is also called the trapezoidal rule. Let the curve y = f (x), a ≤ x ≤ b, be
approximated by the line joining the points P(a, f (a)), Q(b, f (b)) on the curve (see
5.1).
Using the Newton’s forward difference formula, the linear polynomial approxima-
tion to f (x), interpolating at the points P(a, f (a)), Q(b, f (b)), is given by

1
f (x) = f (x0 ) + (x − x0 )∆ f (x0 )
h

94
 Figure 5.1: Trapizium rule

where x0 = a, x1 = b and h = b − a. Substituting in

we obtain
Z b Z x1 Z x1 Z x1
1
I= f (x)dx = f (x)dx = f (x0 ) dx + [ (x − x0 )dx]∆ f0
a x0 x0 h x0

1 1
= (x1 − x0 ) f (x0 ) + [ (x − x0 )2 ]xx10 ∆ f (x0 )
h 2
1
= (x1 − x0 ) f (x0 ) + [ f (x1 ) − f (x0 )](x1 –x0 )2
2h
h
= h f (x0 ) + [ f (x1 ) − f (x0 )]
2
h (b − a)
= [ f (x1 ) + f (x0 )] = [ f (b) + f (a)].
2 2
The trapezium rule is given by
Z b
h b−a
I= f (x)dx = [ f (x1 ) + f (x0 )] = [ f (b) + f (a)] (5.19)
a 2 2

Hence, the trapezium rule integrates exactly polynomial of degree ≤ 1, and the
method is of order 1.
Composite trapezium rule
Let the interval [a, b] be subdivided into N equal parts of length h. That is,
h = b−a
N . The nodal points are given by

a = x0 , x1 = x0 + h, x2 = x0 + 2h, ..., xN = x0 + Nh = b.

We write
Z b Z xN Z x1 Z x2 Z xN
f (x)dx = f (x)dx = f (x)dx + f (x)dx + · · · + f (x)dx
a x0 x0 x1 xN−1

95
There are N integrals. Using the trapezoidal rule to evaluate each integral, we get
the composite trapezoidal rule as
Rb
af (x)dx = 2h [{ f (x0 ) + f (x1 )} + { f (x1 ) + f (x2 )} + · · · + { f (xN−1 ) + f (xN )}]
= h2 [ f (x0 ) + 2{ f (x1 ) + f (x2 ) + ... + f (xN−1 )} + f (xN )].
(5.20)
The composite trapezium rule is also of order 1.
The formula (5.20) is called Trapezoidal rule for numerical integration. The error
committed in this formula is given by

h3 00 (b − a)3 00
E ≈ − f (ξ ) = − f (ξ ), (5.21)
12 12n2

where
a = x0 < ξ < xn = b.

Example 5.4.3 Find the approximate value of I = 01 1+x 1 R

dx , using the trapezium
rule with 2, 4 and 8 equal subintervals. Using the exact solution, find the absolute
errors.

Solution With N = 2, 4 and 8, we have the following step lengths and nodal points.
1
N = 2, h = b−a
N = 2 the nodes are 0, 0.5, 1.0
1
N = 4, h = b−a
N = 4 the nodes are 0, 0.25, 0.5, 0.75, 1.0
1
N = 8, h = b−a
N = 8 the nodes are 0, 0.125, 0.25, 0.375, 0.5, 0.675, 0.75, 0.875, 1.0
We have the following tables of values.

N=2:x 0 0.5 1.0

f (x) 1.0 0.666667 0.5

N = 4: We require the above values. The additional values required are the
following:
x 0.25 0.75
f (x) 0.8 0.571429
N = 8: We require the above values. The additional values required are the follow-
ing:
x 0.125 0.375 0.625 0.875
f (x) 0.888889 0.727273 0.615385 0.533333

96
Now, we compute the value of the integral

h
N =2: I1 = [ f (0)+2 f (0.5)+ f (1.0)] = 0.25[1.0+2(0.666667)+0.5] = 0.708334.
2

h
N=4: I2 = [ f (0) + 2 f (0.25) + f (0.5) + f (0.75) + f (1.0)]
2
= 0.125[1.0+20.8 + 0.666667 + 0.571429+0.5] = 0.697024.
h
N=8: I3 = [ f (0) + 2( f (0.125) + f (0.25) + f (0.375) + f (0.5)
2
+ f (0.625) + f (0.75) + f (0.875)) + f (1.0)]

= 0.0625[1.0+2(0.888889+0.8+0.727273+0.666667+0.615385

+ 0.571429 + 0.533333) + 0.5] = 0.694122.

The exact value of the integral is I = ln2 = 0.693147. The errors in the solutions
are the following:

|Exact − I1 | = |0.693147 − 0.708334| = 0.015187

|Exact − I2 | = |0.693147 − 0.697024| = 0.003877

|Exact − I3 | = |0.693147 − 0.694122| = 0.000975.

Example 5.4.4 Evaluate I = 12 5+3x

1
R
dx with 4 and 8 subintervals using the trapez-
ium rule. Compare with the exact solution and find the absolute errors in the solu-
tions. Comment on the magnitudes of the errors obtained. Find the bound on the
errors.

Solution: With N = 4 and 8, we have the following step lengths and nodal points.
1
N = 4 : h = b−aN = 4 . The nodes are 1, 1.25, 1.5, 1.75, 2.0.

1
N = 8 : h = b−aN = 8 . The nodes are 1, 1.125, 1.25, 1.375, 1.5, 1.675, 1.75, 1.875, 2.0.
We have the following tables of values.
N = 4:
x 1.0 1.25 1.5 1.75 2.0
f (x) 0.125 0.11429 0.10526 0.09756 0.09091

97
 N = 8: We require the above values. The additional values required are the
following
x 1.125 1.375 1.625 1.875
f (x) 0.11940 0.10959 0.10127 0.09412
Now, we compute the value of the integral.

h
N=4: I1 = [ f (1) + 2( f (1.25) + f (1.5) + f (1.75)) + f (2.0)]
2

= 0.125[0.125+2(0.11429+0.10526+0.09756)+0.09091] = 0.10627.
h
N=8: I2 = [ f (1) + 2( f (1.125) + f (1.25) + f (1.375) + f (1.5)
2
+ f (1.625) + f (1.75) + f (1.875)) + f (2.0)]

= 0.0625[0.125 + 2(0.11940 + 0.11429 + 0.10959 + 0.10526 + 0.10127

+ 0.09756 + 0.09412) + 0.09091] = 0.10618.

The exact value of the integral is

1 1
I = [5 + 3x]21 = [ln(11) − ln(8)] = 010615.
3 3

The errors in the solutions are the following:

|Exact − I1 | = |0.10615 − 0.10627| = 0.00012.

|Exact − I2 | = |0.10615 − 0.10618| = 0.00003.

Exercises 5.4.2 Using the trapezium rule, evaluate the integral I = 01 x2 +6x+10
dxR
,
with 2 and 4 subintervals. Compare with the exact solution. Comment on the
magnitudes of the errors obtained.

Exercises 5.4.3 The velocity of a particle which starts from rest is given by the
following table.

t(sec) 0 2 4 6 8 10 12 14 16 18 20
v( f t/sec) 0 16 29 40 46 51 32 18 8 3 0

Evaluate using trapezium rule, the total distance travelled in 20 seconds.

98
Note: Trapezoidal rule can be applied to any number of subintervals odd or even.

Simpson’s 1/3 Rule

In the previous section, we have shown that the trapezium rule of integration in-
tegrates exactly polynomials of degree 1, that is, the order of the formula is 1.
In many science and engineering applications, we require methods which produce
more accurate results. One such method is the Simpson’s 1/3 rule. Let the interval
[a, b] be subdivided into two equal parts with step length h = b−a 2 . We have three
a+b
abscissas x0 = a, x1 = 2 , and x2 = b. Then, P(x0 , f (x0 )), Q(x1 f (x1 )), R(x2 , f (x2 ))
are three points on the curve y = f (x). We approximate the curve y = f (x), axb,
by the parabola joining the points P, Q, R, that is, we approximate the given curve
by a polynomial of degree 2. Using the Newton’s forward difference formula, the
quadratic polynomial approximation to f(x), interpolating at the points P(x0 , f (x0 )), Q(x1 f (x1 )), R(x2 , f (x
is given by

1 1
f (x) = f (x0 ) + (x − x0 )∆ f (x0 ) + 2 (x − x0 )(x − x1 )∆ f (x0 ).
h 2h

Substituting in (5.18), we obtain

Z b Z x2 Z x2
1 1
f (x)dx = f (x)dx = [ f (x0 )+ (x−x0 )∆ f (x0 )+ 2 (x−x0 )(x−x1 )∆ f (x0 )]dx
a x0 x0 h 2h
(5.22)

99
Chapter 6

First Order Initial Value Problems

6.1 Definition of ODE and examples
Definition: An ordinary differential equation(ODE) is a relation between a function
y(t), its derivatives y0 , y00 , · · · , y(m) where m denotes the order of the highest deriva-
tive, and the variable t upon which the function y and all its derivatives depend. The
most general form of an ODE is given by

Φ(t, y, y0 , y00 , · · · , y(m) ) = 0. (6.1)

6.2 Order of a differential equation, linear and non-

linear ODE
Definition: The order of the differential equation (DE) is the order of its highest
derivative.

Example 6.2.1 A first order ordinary differential equation in the unknown y is

y0 (t) = f (t, y(t)) (6.2)

where y : R → R is the unknown function and f : R2 → R is a given function. The

equation
in (6.1) is called linear iff the function with values f (t, y) is linear on its second
argument; that is, there exist functions a, b : R → R such that

y0 (t) = a(t)y(t) + b(t), f (t, y) = a(t)y + b(t).

If P = P(x) and Q = Q(x) are functions of x only, then

dx/dy + Py = Q

100
is called a linear differential equation order 1.
We can solve these linear DEs using an integrating factor. For linear DEs of
order 1, the integrating factor is:

R
P dx
e

Example 6.2.2 Solve

dy 3
− y=7
dx x

Solution:
Here,
3
P(x) = −
x
and
Q(x) = 7

Now for the integrating factor:

− 3x dx x−3
R R
IF = e P dx
=e = e−3 ln x = eln = x−3

We need to apply the following formula:

R Z  R 
P dx P dx
ye = Qe dx

Example 6.2.3 A second order differential equation is an equation of the form

F(x, y, y0 , y00 ) = 0.

Definition: The degree of the DE is the degree of the derivative of the highest order
after the equation has been rationalized in derivatives. An ODE is of two types
namely linear and non-linear.
Definition:An ODE is said to be linear if no product of the dependent variable y(t)
with itself or any one of its derivatives occur. Otherwise, it is non-linear. A linear
DE of order m can be expressed in the form
m
∑ Φ p(t)y p(t) = r(t)
p=0

101
in which Φ p (t) are known functions. If the general non-linear DE (6.1) of order m
can be written as
y(m) = F(t, y, y0 , . . . , y(m−1) ) (6.3)

then the equation in (6.3) is called a canonical (or natural) representation of the DE
(6.1). In such form, the highest order derivative is expressed in terms of the function
y, the lower order derivatives and the independent variable t.

Example 6.2.4 A linear differential equations are

a. a(x)y0 + b(x)y = c(x) is first order ,

b. a(x)y00 + b(x)y0 + c(x)y = d(x) is second order,

where a(x), b(x), c(x), d(x) are known functions of x

Example 6.2.5 Some non-linear differential equations

a. (y00 )2 = xy

b. y0 y = −x

c. y0 = siny

6.3 Nature of Solutions of ODE: particular and gen-

eral solutions
{x, yz

102