Lecture Notes

on

Numerical Methods

Er. Narayan Sapkota, [Link]

Everest Engineering College | Kathmandu University

January 28, 2026

 Disclaimer
This lecture notes have been prepared for the sole purpose of assisting students in their coursework. They
are intended to supplement, not replace, the readings and lectures provided by the instructor. While every
effort has been made to ensure the accuracy and completeness of the information contained in these notes, it
is not guaranteed to be complete, accurate, or up-to-date. Therefore, they should not be used as a substitute
for attending lectures, participating in class discussions, or completing assigned readings.
The author(s) of these notes accept no responsibility for any errors or omissions, or for any actions taken
or not taken based on the information contained herein. Furthermore, the author(s) make no express or
implied warranties of any kind, including but not limited to, warranties of merchant-ability or fitness for a
particular purpose. The reader assumes all responsibility for the use of the information contained in these
notes.
In no event shall the author(s) be liable for any damages whatsoever (including, without limitation, damages
for loss of profits, business interruption, or loss of information) arising out of the use or inability to use these
notes, even if the author(s) have been advised of the possibility of such damages.
These notes are provided on an "as is" basis, and the author(s) make no representations or warranties of
any kind, express or implied, with respect to the contents of these notes. The author(s) reserve the right to
make changes to these notes without notice.
By using these notes, you acknowledge that you have read and understood this disclaimer, and agree to be
bound by its terms and conditions.

i
Prepared By: Er. Narayan Sapkota, [Link]. ii Everest Engineering College
 Contents

0 Preliminaries 1
0.1 Mathematical Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
0.1.1 Derivative and Anti-derivative . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
0.1.2 Limits and Continuity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
0.1.3 Intermediate Value Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
0.1.4 Mean Value Theorem for Derivatives (Lagrange Theorem) . . . . . . . . . . . . . . . . 2
0.1.5 Rolle’s Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
0.1.6 Integrals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
0.1.7 Binomial Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
0.1.8 Taylor’s Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
0.2 Exact and Approximate Numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
0.3 Errors in Numerical Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
0.3.1 Sources of Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
0.3.2 Types of Errors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
0.3.3 Accuracy and Precision Trade-off . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
0.3.4 Significant Digits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
0.4 Methods of Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
0.5 Classification of Iterative Root finding Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 10
0.5.1 Order of Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
0.5.2 Algorithm to Find the Order of Convergence . . . . . . . . . . . . . . . . . . . . . . . 11

1 Solutions of Non-Linear Equations 13

1.1 Bisection Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.2 The Secant Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.3 The Newton-Raphson Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.4 Iteration Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2 Interpolation and Approximation 25

2.1 The Lagrange Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.1.2 Linear Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.1.3 Quadratic Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.1.4 General Lagrange Interpolation Polynomial . . . . . . . . . . . . . . . . . . . . . . . . 27
2.1.5 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.1.6 Inverse Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.1.7 Exercise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.2 Finite Differences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.2.1 Forward Differences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.2.2 Backward Differences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.2.3 Central Differences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.2.4 Detection of Errors by use of Difference Tables . . . . . . . . . . . . . . . . . . . . . . 33
2.3 Error Detection in Difference Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.4 Differences of a Polynomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.5 Newton’s Interpolation using Difference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.5.1 Newton’s Forward Difference Interpolation . . . . . . . . . . . . . . . . . . . . . . . . 35
2.5.2 Newton’s Backward Difference Interpolation . . . . . . . . . . . . . . . . . . . . . . . . 35
2.5.3 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.6 Newton’s Divided Difference Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.6.2 Newton’s Forward Interpolation for Divided Differences . . . . . . . . . . . . . . . . . 38

iii
CONTENTS

2.6.3 Newton’s Backward Interpolation for Divided Differences . . . . . . . . . . . . . . . . 38

2.6.4 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.7 Curve Fitting: Least Squares Lines for Linear and Nonlinear Data . . . . . . . . . . . . . . . 40
2.8 Cubic Spline Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

3 Numerical Differentiation and Integration 41

3.1 Numerical Differentiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2 Newton’s Differentiation Formulas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2.1 Derivatives Using Newton’s Forward Difference Interpolation Formula . . . . . . . . . 41
3.2.2 Derivatives Using Newton’s Backward Difference Interpolation Formula . . . . . . . . 43
3.3 Numerical Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.4 Newton-Cotes Quadrature Formula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.4.1 Trapezoidal Rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.4.2 Simpson’s One-Third Rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.4.3 Simpson’s Three-Eighth Rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.5 Romberg’s Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.6 Gaussian Integration algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

4 Solution of Linear Algebraic Equations 59

4.1 Matrices and Their Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.1.1 Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.1.2 Types of Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.1.3 Matrix Operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.1.4 Positive Definite Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.2 Solution of Linear Simultaneous Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.2.1 Gauss Elimination Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.2.2 Partial Pivoting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.2.3 Complete/Full Pivoting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.2.4 Gauss-Jordan Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.3 Method of Factorization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.3.1 LU Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.3.2 Doolittle Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.3.3 Crout’s Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.3.4 Cholesky’s Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.4 The Inverse of a Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.4.1 Gauss-Jordan Method to Find the Inverse of a Matrix . . . . . . . . . . . . . . . . . . 81
4.5 Ill-Conditioned Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.6 Iterative Methods of Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.6.1 Jacobi’s Iteration Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.6.2 Gauss–Seidel Iteration Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.7 Eigenvalues and Eigenvectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.7.1 Cayley-Hamilton Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
4.7.2 Power Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
4.8 Assignment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.8.1 Gauss Elimination Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.8.2 Gauss-Jordan Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
4.8.3 Factorization Method (LU / Doolittle / Crout) . . . . . . . . . . . . . . . . . . . . . . 93
4.8.4 Gauss-Seidel Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
4.8.5 Jacobi Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.8.6 Gauss-Jordan Method – Inverse of a Matrix . . . . . . . . . . . . . . . . . . . . . . . . 95
4.8.7 Power Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

5 Solution of Ordinary Differential Equations 97

5.1 Review of Ordinary Differential Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

Prepared By: Er. Narayan Sapkota, [Link]. iv Everest Engineering College

CONTENTS

5.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
5.1.2 First-Order ODEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
5.1.3 Second-Order ODEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
5.2 Euler Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
5.3 Runge-Kutta Method for First-Order Differential Equation . . . . . . . . . . . . . . . . . . . 101
5.4 Runge-Kutta Method for Second-Order Differential Equation . . . . . . . . . . . . . . . . . . 106
5.5 Boundary Value Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
5.6 Shooting Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

6 Solutions of PDEs 113

6.1 Classification of PDE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
6.2 Finite Difference Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
6.3 Elliptic Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
6.3.1 Laplace Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
6.3.2 Poisson Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
6.4 Parabolic Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
6.4.1 1D Heat Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
6.5 Exercise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131

Prepared By: Er. Narayan Sapkota, [Link]. v Everest Engineering College

 Chapter 0
Preliminaries

0.1 Mathematical Preliminaries

0.1.1 Derivative and Anti-derivative
Consider the function f (x) = cos(x). Its derivative f ′ (x) = − sin(x), and its anti-derivative F (x) =
sin(x) + C. The former is used to determine the slope f ′ (x0 ) of the curve at a point (x0 , f (x0 )), and
the latter is used to compute the area under the curve for a < x < b.
Z π/2
area = cos(x) dx = F (π/2) − F (0) = sin(π/2) − 0 = 1
0

Figure 1: The slope at the point ( π2 , 0) is m =

Figure 2: The area under the curve for 0 ≤ x ≤ π/2
f ′ ( π2 ) = −1 and can be used to find the tangent
is computed using an integral
line at this point

0.1.2 Limits and Continuity

Limits: The limit of a function f (x) as x approaches a certain value c, denoted by limx→c f (x), represents
the value that f (x) approaches as x gets arbitrarily close to c. Formally, we say that limx→c f (x) = L if for
every positive number ε, there exists a positive number δ such that if 0 < |x − c| < δ, then |f (x) − L| < ε.
This definition captures the intuitive idea that as x gets closer and closer to c, the values of f (x) get closer
and closer to L.

Continuity: A function f (x) is said to be continuous at a point c if the following three conditions are met:

1. f (c) is defined.

2. limx→c f (x) exists.

3. limx→c f (x) = f (c).

In other words, a function is continuous at a point c if its value at c is equal to the limit of the function as
x approaches c. A function is continuous on an interval if it is continuous at every point in that interval.

1
CHAPTER 0. PRELIMINARIES

0.1.3 Intermediate Value Theorem

Let f (x) is a continuous function on the closed interval [a, b] and let u be any number between f (a)
and f (b), we assume f (a) < f (b). Then, there exists at least one number c in the open interval (a, b)
such that f (c) = u.

This has two important corollaries:

1. If a continuous function has values of opposite sign inside an interval, then it has a root in that interval
(Bolzano–Weierstrass Theorem).
2. The image of a continuous function over an interval is itself an interval.

Figure 3: Intermediate Value Theorem

Proof: For Proof click this link

Geometrically, the IVT states that if we have a continuous curve connecting two points (a, f (a)) and (b, f (b)),
and a horizontal line y = u intersects the curve between the two points, then there must be at least one
point c on the curve where the y-coordinate is equal to u.

It ensures the existence of roots for certain equations and guarantees the existence of solutions to various
problems in physics, engineering, and other scientific fields where continuous functions model phenomena.

If a continuous function f on [a, b] satisfies f (a) · f (b) < 0 (i.e., f (a) and f (b) have opposite signs),
then there exists at least one ξ ∈ (a, b) such that f (ξ) = 0.

0.1.4 Mean Value Theorem for Derivatives (Lagrange Theorem)

If f is a continuous function on the closed interval [a, b] and differentiable on the open interval (a, b),
then there exists a point c in (a, b) such that the tangent at c is parallel to the secant line through
the endpoints (a, f (a)) and (b, f (b)), that is,

f (b) − f (a)
f ′ (c) = , a<c<b
b−a

Proof: For proof click this

Geometrically, the MVT says that there is at least one number c in the open interval (a, b) such that the
slope of the tangent line to the graph of y = f (x) at the point (c, f (c)) equals the slope of the secant line
through the points (a, f (a)) and (b, f (b)).

Prepared By: Er. Narayan Sapkota, [Link]. 2 Everest Engineering College

CHAPTER 0. PRELIMINARIES

Figure 4: For any function that is continuous on [a, b] and differentiable on (a, b), there exists some c in the
interval (a, b) such that the secant joining the endpoints of the interval [a, b] is parallel to the tangent at c.

0.1.5 Rolle’s Theorem

If f (x) is a real-valued function that is continuous on the closed interval [a, b], differentiable on the
open interval (a, b), and f (a) = f (b), then there exists at least one point c in the open interval (a, b)
such that f ′ (c) = 0.

Rolle’s Theorem is a special case of the MVT, where the average rate of change of a function over an interval
is zero.

0.1.6 Integrals
First Fundamental Theorem If f is continuous over [a, b] and F is any anti-derivative of f on [a, b], then
Z b
f (x) dx = F (b) − F (a) where F ′ (x) = f (x)
a

Second Fundamental Theorem If f is continuous over [a, b] x ∈ (a, b) then,

Z x
d
f (t) dt = f (x)
dx a

Mean Value Theorem for Integrals Assume that f ∈ C[a, b]. Then there exists a number c, with
c ∈ (a, b), such that
1
Z b
f (x) dx = f (c)
b−a a
The value f (c) is the average value of f over the interval [a, b].

Prepared By: Er. Narayan Sapkota, [Link]. 3 Everest Engineering College

CHAPTER 0. PRELIMINARIES

Figure 5: If a real-valued function f is continuous on a closed interval [a, b], differentiable on the open interval
(a, b), and f (a) = f (b), then there exists a c in the open interval (a, b) such that f ′ (c) = 0.

0.1.7 Binomial Theorem

n  
X n
(x + y) = n
xn−k y k
k
k=0

where n is a non-negative integer, n

represents the binomial coefficient, and x and y are any real numbers.


k

The binomial coefficient n

represents the number of ways to choose k elements from a set of n elements


k
and is defined as:  
n n!
=
k k!(n − k)!

0.1.8 Taylor’s Theorem

Let f be a function with n + 1 continuous derivatives on an interval containing a. Then, for any x in this
interval,
n
X f (k) (a)
f (x) = (x − a)k + Rn (x)
k!
k=0

f ′′ (a) f ′′′ (a) f (n) (a)

f (x) = f (a) + f ′ (a)(x − a) + (x − a)2 + (x − a)3 + · · · + (x − a)n + Rn (x)
2! 3! n!
where the remainder Rn (x) (in Lagrange form) is:

f (n+1) (c)
Rn (x) = (x − a)n+1 , for some c between a and x.
(n + 1)!

The theorem requires:

• f (n+1) must exist (but need not be continuous) on I
• The remainder term shows the error when approximating f by its n-th degree Taylor polynomial
h2 ′′ h3 hn (n)
f (x + h) = f (x) + hf ′ (x) + f (x) + f ′′′ (x) + · · · + f (x) + Rn (x)
2! 3! n!
and
f (n+1) (c) n+1
Rn (x) = h
(n + 1)!

Prepared By: Er. Narayan Sapkota, [Link]. 4 Everest Engineering College

CHAPTER 0. PRELIMINARIES

Taylor Series for Some standard relations

∞
x2 x3 x4 X xn
e =1+x+
x
+ + + ... =
2! 3! 4! n=0
n!

∞
x3 x5 x7 x9 X x2n+1
sin(x) = x − + − + − ... = (−1)n
3! 5! 7! 9! n=0
(2n + 1)!

∞
x2 x4 x6 x8 X x2n
cos(x) = 1 − + − + − ... = (−1)n
2! 4! 6! 8! n=0
(2n)!

0.2 Exact and Approximate Numbers

Exact numbers are values that are known with absolute certainty, having no rounding errors. These
are typically integers, counting numbers, or values defined by rules or constants.

• The number of items in a set (e.g., 5 apples).

• Defined constants (e.g., 100 cm = 1 m).
• Mathematical constants like π (in theoretical contexts).

Approximate numbers are those that arise from measurements or estimations, where the value is
known with some degree of uncertainty or rounding.

• The approximation of π ≈ 3.14 for practical purposes.

• Measurement results (e.g., 5.67 cm, temperature 32.5◦ C).

0.3 Errors in Numerical Calculations

0.3.1 Sources of Errors
When a computational procedure is involved in solving a scientific mathematical problem, errors often will
be involved in the process. A rough classification of the kinds of original errors that might occur is as follows:

Modeling Errors: Mathematical modeling is a process when mathematical equations are used to represent
a physical system. This modeling introduces errors called modeling errors. For example, consider a simple
pendulum system modeled using the differential equation for simple harmonic motion. If the air resistance
or frictional forces are not accurately accounted for in the model, the predictions of the pendulum’s behavior
may deviate from experimental observations.

Blunders and Mistakes: Blunders occur at any stage of the mathematical modeling process and consist
of all other components of error. Blunders can be avoided by sound knowledge of fundamental principles
and by taking proper care in approaching and designing a solution. For instance, in an engineering design
project, incorrectly applying a fundamental equation such as Newton’s second law of motion can lead to
significant errors in predicting the behavior of a mechanical system. Mistakes, on the other hand, are due
to programming errors. For example, a typographical error in a computer program calculating numerical
integrals can result in incorrect results.

Prepared By: Er. Narayan Sapkota, [Link]. 5 Everest Engineering College

CHAPTER 0. PRELIMINARIES

Machine Representation and Arithmetic Errors: These errors are inevitable when using floating-
point arithmetic while using computers or calculators. Examples include rounding and chopping errors. For
instance, when representing the irrational number π in a computer’s finite memory, only a finite number of
digits can be stored. Thus, calculations involving π may introduce rounding errors that accumulate over
successive computations, leading to discrepancies between the computed and exact values.

Mathematical Approximation Errors: This error is also known as truncation error or discretization
error. These errors arise when an approximate formulation is made to a problem that otherwise cannot be
solved exactly. For example, when solving a differential equation numerically using finite difference methods,
the continuous function is discretized into a finite number of points. The error introduced by this discretiza-
tion process is the truncation error, and it affects the accuracy of the numerical solution.

Propagated Error: Propagated error is more nuanced compared to other types of errors. It refers to
errors that arise in subsequent steps of a process due to an earlier mistake, compounding the local errors
present at each stage. This is somewhat similar to errors in initial conditions. Some root-finding techniques
discover additional zeros by modifying the function to eliminate the first root, a process known as reducing
or deflating the equation. In this context, the reduced equations carry forward the errors from the earlier
stages. To address this issue, it is important to verify the later results using the original equation.

In numerical methods discussed in later chapters, propagated error is particularly significant. If errors con-
tinue to increase throughout the process, they can eventually obscure the true value, making the method
invalid; such methods are considered unstable. A stable method—one that is desirable—ensures that errors
introduced at early stages diminish as the process continues.

Each type of error can occur independently of the others, although they may sometimes interact. For example,
round-off error can occur even without truncation error, as in analytical methods. Similarly, truncation errors
can lead to inaccuracies even if calculations are performed with perfect precision. Standard error analysis of
a numerical method often treats truncation error as though such perfect precision exists.

0.3.2 Types of Errors

Absolute and Relative Errors: If XE is the exact or true value of a quantity and XA is its approximate
value, then |XE − XA | is called the absolute error Ea . Therefore, absolute error:

Ea = |XE − XA |

The relative error is defined by

|XE − XA |
Er =
|XE |
provided XE ̸= 0 or XE is not too close to zero.

Percentage relative error is

Ep = Er × 100%.

Suppose a measurement of a length is found to be 100 cm, and the true length is 98 cm. Then the absolute
error is |100 − 98| = 2 cm, and the relative error is 100−98
98 × 100% = 2.04%.

Inherent Errors:: Inherent errors are the errors that preexist in the problem statement itself before its
solution is obtained. These errors exist because of the approximate nature of the data or due to the limita-
tions of the calculations using digital computers. Inherent errors cannot be completely eliminated but can

Prepared By: Er. Narayan Sapkota, [Link]. 6 Everest Engineering College

CHAPTER 0. PRELIMINARIES

be minimized if we select better data or employ high-precision computer computations.

Consider the problem of calculating the area of a circle using the formula A = πr2 . In this formula, π is an
irrational number that cannot be represented exactly in decimal form. Therefore, any calculation involving
π introduces inherent errors due to the finite precision of computer arithmetic.
Suppose we want to calculate the area of a circle with a radius of 2 meters. Using a standard value of π
truncated to four decimal places (3.1416), we get:
A = 3.1416 × 22 = 12.5664 square meters

However, the exact area of the circle, if π were known to infinite precision, would be:
Aexact = π × 22 = π × 4 = 12.5663706 . . . square meters

The inherent error in this calculation arises from the approximation of π. While it may seem small in this
example, for more complex calculations or when high precision is required, the impact of inherent errors can
be significant.

Round-off Errors: Round-off error happens when inaccuracies arise due to a finite number of digits of pre-
cision used to represent numbers. All computers represent numbers, except for integers and some fractions,
with imprecision. Digital computers use floating-point numbers of fixed word length. This type of repre-
sentation will not express the exact or true values correctly. Error introduced by the omission of significant
figures due to computer imperfection is called the round-off error.

Round-off errors are avoidable in most computations. When n digits are used to represent a real number,
one method to avoid the error is to keep the first n digits and chop off all remaining digits. Another method
is to round to the n-th digit by examining the values of the remaining digits.

Methods to Mitigate Round-off Errors

• Truncation: Keep only the first n digits and discard all remaining digits.
• Rounding: Round to the n-th digit by examining the values of the remaining digits. If the (n + 1)-th
digit is 5 or greater, increase the n-th digit by 1; otherwise, leave the n-th digit unchanged.
Using these methods, it’s possible to minimize the impact of round-off errors in computational calculations.
However, it’s important to note that round-off errors can never be completely eliminated, and the choice of
method may depend on the specific requirements of the computation and the level of precision needed.

Truncation Errors: Truncation errors are defined as those errors that result from using an approximation
in place of an exact mathematical procedure. Truncation error results from terminating after a finite number
of terms, known as formula truncation error or simply truncation error.
Example: Consider the function f (x) = ex and its Taylor series expansion around x = 0:
∞
X xn
ex =
n=0
n!

Suppose we want to approximate ex by truncating the series after a finite number of terms. Let’s say we
truncate the series after the n-th term. Then, the truncation error Tn (x) is given by the difference between
the exact function and its truncated approximation:
n
X xk
Tn (x) = e −
x
k!
k=0

Prepared By: Er. Narayan Sapkota, [Link]. 7 Everest Engineering College

CHAPTER 0. PRELIMINARIES

Figure 6: Concepts of accuracy and precisions

As an example, let’s truncate the series expansion of ex after the first two terms (n = 1):

T1 (x) = ex − (1 + x)

This truncation error represents the discrepancy between the exact function ex and its approximation 1 + x
when only the first two terms of the Taylor series expansion are used.

0.3.3 Accuracy and Precision Trade-off

In many scientific and engineering applications, there exists a trade-off between accuracy and precision.
These two terms are often used interchangeably, but they have distinct meanings:
• Accuracy refers to how close a measured or calculated value is to the true or accepted value. It indicates
the absence of systematic errors or biases in the measurement or calculation process. High accuracy
implies that the results are reliable and trustworthy.
• Precision refers to the degree of repeatability or reproducibility of measurements or calculations. It
indicates the level of detail or granularity in the data. High precision implies that the results are
consistent and reproducible, even if they may not be accurate.
In practice, achieving both high accuracy and high precision simultaneously can be challenging, and there
is often a tradeoff between the two:
• High Accuracy, Low Precision: In some cases, it may be possible to obtain highly accurate results by
carefully calibrating instruments, minimizing systematic errors, and using accurate measurement tech-
niques. However, this may come at the cost of reduced precision, as small variations in measurements
may still occur.
• High Precision, Low Accuracy: Conversely, it may be possible to achieve high precision by increasing
the resolution of instruments or using advanced computational techniques to perform calculations with
many significant figures. However, this may come at the cost of reduced accuracy, as systematic errors
or biases in the measurement or calculation process may still exist.

0.3.4 Significant Digits

• Significant digits, also known as significant figures, are digits in a numerical value that carry meaningful
information about the precision of the measurement or calculation.

Prepared By: Er. Narayan Sapkota, [Link]. 8 Everest Engineering College

CHAPTER 0. PRELIMINARIES

• Thus, the numbers 3.1416, 0.66667, and 4.0687 contain five significant digits each. The number 0.00023
has, however, only two significant digits, namely 2 and 3, since the zeros serve only to fix the position
of the decimal point.
• In cases of ambiguity, the scientific notation should be used. For example, in the number 25, 600, the
number of significant figures is uncertain, whereas the numbers 2.56×104 , 2.560×104 , and 2.5600×104
have three, four, and five significant digits, respectively.
The key principles behind significant digits are:
1. Non-zero digits: All non-zero digits are considered significant. For example, in the number 1234, all
four digits (1, 2, 3, and 4) are significant.
2. Leading zeros: Leading zeros (zeros to the left of the first non-zero digit) are not significant. For
example, in the number 0.0023, only the digits 2 and 3 are significant.
3. Captive zeros: Captive zeros (zeros between non-zero digits) are significant. For example, in the
number 506, all three digits (5, 0, and 6) are significant.
4. Trailing zeros: Trailing zeros (zeros to the right of the last non-zero digit) are significant if they are
after the decimal point or if they appear in a number with a decimal point. For example, in the number
4.50, both digits 4 and 5 are significant.
5. Trailing zeros without a decimal point: Trailing zeros at the end of a whole number without a decimal
point are not considered significant. For example, in the number 200, only the digits 2 and 0 are
significant.

0.4 Methods of Solution

We consider the methods for determining the roots of the equation

f (x) = 0 (1)

which may be given explicitly as a polynomial of degree n in x or f (x) may be defined implicitly as a tran-
scendental function. A transcendental (1) may have no root, a finite or an infinite number of real and/or
complex roots while a polynomial equation (1) has exactly n (real and/or complex) roots.

If the function f (x) changes sign in any one of the intervals [x∗ − ϵ, x∗ ], [x∗ , x∗ + ϵ], then x∗ defines an
approximation to the root of f (x) with accuracy ϵ. This is known as the intermediate value theorem. Hence,
if the interval [a, b] containing x∗ and α, where α is the exact root of (1), is sufficiently small, then

|x∗ − α| ≤ b − a

can be used as a measure of the error.

There are two types of methods that can be used to find the roots of the equation (1):
• Direct Methods
– Give the exact value of the roots (in the absence of round-off errors) in a finite number of steps
– Determine all the roots at the same time.
• Iterative Methods
– Based on the idea of successive approximations.
– Starting with one or more initial approximations to the root, obtain a sequence of iterates {xn }
which, in the limit, converges to the root.
– Determine one or two roots at a time.

Prepared By: Er. Narayan Sapkota, [Link]. 9 Everest Engineering College

CHAPTER 0. PRELIMINARIES

0.5 Classification of Iterative Root finding Methods

• Bracketing Methods: Require the limits between which the root lies, Assumption that the function
changes sign in the vicinity of a root.
– Bisection method
– False position method
• Open methods: Require the initial estimation of the solution, Use a single starting value or two values
that do not necessarily bracket the root.
– Newton-Raphson method
– Secant Method
– Fixed-Point Iteration Method

A sequence of iterates {xn } is said to converge to the root ξ if

lim |xn − ξ| = 0.
k→∞

If xn , xn−1 , . . . , xn−m+1 are m approximates to a root, then we write an iteration method in the form

xn+1 = ϕ(xn , xn−1 , . . . , xn−m+1 ) (2)

where we have written the equation (1) in the equivalent form

x = ϕ(x)

The function ϕ is called the iteration function. For m = 1, we get the one-point iteration method

xn+1 = ϕ(xn ), n = 0, 1, . . . (3)

If ϕ(x) is continuous in the interval [a, b] that contains the root and |ϕ′ (x)| ≤ c < 1 in this interval,
then for any choice of x0 ∈ [a, b], the sequence of iterates {xk } obtained from (3) converges to the
root of x = ϕ(x) or f (x) = 0.

Thus, for any iterative method of the form (2) or (3), we need the iteration function ϕ(x) and one or more
initial approximations to the root.

In practical applications, it is not always possible to find ξ exactly. We therefore attempt to obtain an
approximate root xn+1 such that
|f (xn+1 )| < ϵ
and/or
|xn+1 − xn | < ϵ
where xn and xn+1 are two consecutive iterates and ϵ is the prescribed error tolerance.

0.5.1 Order of Convergence

• Measure of how quickly an (iterative) method converges to the solution of a problem.
• Quantifies the rate at which the error decreases as the number of iterations increases.
• How rapidly successive iterates approach the true root of the equation

Prepared By: Er. Narayan Sapkota, [Link]. 10 Everest Engineering College

CHAPTER 0. PRELIMINARIES

Definition An iterative method is said to be of order p or has the rate of convergence p, if p is the
largest positive real number for which
|ϵn+1 | ≤ c|ϵn |p
where ϵn = xn − ξ is the error in the n-th iterate.

Definition If a sequence x1 , x2 , x3 , · · · , xn converges to a value α and if there exist real numbers C > 0
and p > 0 such that

|xn+1 − α|
lim =C (4)
n→∞ |xn − α|p
then we say that p is the rate of convergence of the sequence.

If α is a real root of the equation f (x) = 0 and x1 , x2 , x3 , · · · , xn be the sequence of approximations, then
the error in the nth approximation is defined by |ϵn | = |xn − α|. Now, equation (4) becomes

|ϵn+1 |
lim =C (5)
n→∞ |ϵn |p

When n is sufficiently large, (5) becomes

|ϵn+1 | = C|ϵn |p (6)

The constant C is called the asymptotic error constant. It depends on various order derivatives of f (x)
evaluated at α and is independent of n. The relation

ϵn+1 = Cϵpn + O(ϵp+1

n )

is called the error equation.

By substituting xi = ξ + ϵi for all i in any iteration method and simplifying, we obtain the error equation
for that method. The value of p thus obtained is called the order of this method.

Linear Convergence (Order 1 (p=1)): If the error en decreases linearly with each iteration n, the
method exhibits linear convergence. Mathematically, en ≈ C en−1 , where C is a constant.

Super-linear Convergence (Order p > 1): If the error en decreases faster than linearly with each
iteration n, the method exhibits super-linear convergence. Mathematically, en ≈ C (en−1 )p , where p > 1
and C is a constant.

0.5.2 Algorithm to Find the Order of Convergence

1. Choose an initial guess x0 for the root.
2. Iterate the method to compute x1 , x2 , . . . , xn .
3. Compute the absolute error en = |xn − α|, where α is the true root (if known).
en+1
4. Compute the ratio of consecutive errors: rn = en .

5. Repeat steps 2-4 for several iterations to obtain a sequence of ratios r1 , r2 , . . . , rm .

6. Plot the logarithm of the ratios log(rn ) against the iteration number n.
7. Estimate the slope of the resulting curve to determine the order of convergence.

Prepared By: Er. Narayan Sapkota, [Link]. 11 Everest Engineering College

CHAPTER 0. PRELIMINARIES

Prepared By: Er. Narayan Sapkota, [Link]. 12 Everest Engineering College

 Chapter 1
Solutions of Non-Linear Equations

1.1 Bisection Method

After a root of f(x) = 0 has been bracketed in the interval (a, b), bisection method can be used to close in on
it. The bisection method accomplishes this by successfully halving the interval until it becomes sufficiently
small. The bisection method is also known as the interval halving method. The bisection method is not
the fastest method available for finding roots of a function, but it is the most reliable method.

Suppose that f (x) is continuous and real-valued function on an interval a ≤ x ≤ b and f (a) and f (b) are of
opposite signs, that is f (a) × f (b) < 0, then there is at least one real root in the interval between a and b.
Let us define another point xm to be the midpoint between a and b. That is, xm = (a + b)/2. Now, there
exist a following three conditions:
1. If f (xm ) = 0, we have a root at xm
2. If f (xm ) × f (a) < 0, there is a root between xm and a.
3. If f (xm ) × f (b) < 0, there is a root between xm and b.

Figure 1.1: Illustration of Bisection Method

a+b
xm = (1.1)
2

To find out how many iterations will be necessary following equation could be use.

 
|b−a|
log ϵ
n≤ (1.2)
log 2

There are several advantages to the bisection method. The method is guaranteed to converge, and always
converges to an answer, provided a root was bracketed in the interval (a, b) to start with. In addition, the

13
CHAPTER 1. SOLUTIONS OF NON-LINEAR EQUATIONS

error bound, given in (1.1), is guaranteed to decrease by one half with each iteration.

The method may fail when the function is tangent to the axis and does not cross the x-axis at f (x) = 0. The
disadvantage of the bisection method is that it generally converges more slowly than most other methods.
For functions f (x) that have a continuous derivative, other methods are usually faster and these methods
may not always converge. When these methods do converge, they are almost always much faster than the
bisection method.

Example 1

Use the bisection method to find a root of the equation x3 − 4x − 9 = 0 accurate to three decimal
places.

Solution: Here, f (x) = x3 − 4x − 9 = 0

f (2) = 23 − 4 × 2 − 9 = −9 < 0 ∴a=2

f (3) = 33 − 4 × 3 − 9 = 6 > 0 ∴b=3

Hence, a root lies between 2 and 3 and ϵ = 0.0001.

   
|b−a| |2−3|
log ϵ log 0.0001
n≤ = = 9.965 ≈ 10
log 2 log 2
Therefore, total iteration required is 10.
We know that,
a+b
xm =
2
The calculations for the root using bisection method are shown in below table.

f (a) < 0(−ve) f (b) > 0(+ve)

Iteration, n a b xm f (xm )
1 2.000000 3.000000 2.500000 -3.375000
2 2.500000 3.000000 2.750000 0.796875
3 2.500000 2.750000 2.625000 -1.412109
4 2.625000 2.750000 2.687500 -0.339111
5 2.687500 2.750000 2.718750 0.220917
6 2.687500 2.718750 2.703125 -0.061077
7 2.703125 2.718750 2.710938 0.079423
8 2.703125 2.710938 2.707031 0.009049
9 2.703125 2.707031 2.705078 -0.026045
10 2.705078 2.707031 2.706055 -0.008506
11 2.706055 2.707031 2.706543 END

The approximated root of x3 − 4x − 9 = 0 using bisection method is 2.706.

Homework 1

Use the bisection method to find a root of the equation 3x + sin(x) − ex = 0 accurate to four decimal
places.

Prepared By: Er. Narayan Sapkota, [Link]. 14 Everest Engineering College

CHAPTER 1. SOLUTIONS OF NON-LINEAR EQUATIONS

Homework 2

Use the bisection method to find a root of the equation xsin(x) + cos(x) = 0 accurate to four decimal
places.

Homework 3

Use the bisection method to find a root of the equation cos(x) = xex accurate to four decimal places.

Prepared By: Er. Narayan Sapkota, [Link]. 15 Everest Engineering College

CHAPTER 1. SOLUTIONS OF NON-LINEAR EQUATIONS

1.2 The Secant Method

Note
A secant line intersects a curve at two or more points. It represents the average rate of change
between those two points.

A tangent line touches the curve at exactly one point. It represents the instantaneous rate of change
or the slope of the curve at that specific point.

This method is an improvement over the method of false position as it does not require the condition
f (xa ) × f (xb ) < 0 unlike the bisection and false position method. 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 approximation. Additionally, it is not necessary for the interval to contain the root.

Figure 1.2: Illustration of Secant Method

The slope of the secant line passing through (x1 , f (x1 )) and (x3 , f (x3 )) = Slope of the secant line passing
through (x2 , f (x2 )) and (x3 , f (x3 )).

Let x3 be the approximated root of f (x) then f (x3 ) ≈ 0. Therefore,

f (x1 ) f (x2 )
=
x1 − x 3 x2 − x 3
Aftermaths,

x2 − x 1
x3 = x2 − f (x2 ) (1.3)
f (x2 ) − f (x1 )

Similarly,

x1 − x 1
x3 = x1 − f (x1 ) (1.4)
f (x1 ) − f (x2 )

The (1.3) and (1.4) are known as the secant formula. The appropriate value of the root can be refined by
repeating this procedure by replacing x1 by x2 and x2 by x3 , respectively, in the (1.3).

In general form, the above (1.3) becomes,

Prepared By: Er. Narayan Sapkota, [Link]. 16 Everest Engineering College

CHAPTER 1. SOLUTIONS OF NON-LINEAR EQUATIONS

xn+1 − xn
xn+2 = xn+1 − f (xn+1 ) (1.5)
f (xn+1 ) − f (xn )

The rate of convergence of the secant method is 1.61.

Example 1

Use the secant method to find a root of the equation x3 − 4x − 9 = 0 accurate to four decimal places.

Solution: Suppose the initial guess for xn = 1 and xn+1 = 4.

We know that, general form of secant method is,
xn+1 − xn
xn+2 = xn+1 − f (xn+1 )
f (xn+1 ) − f (xn )

Iteration xn f (xn ) xn+1 f (xn+1 ) xn+2 f (xn+2 )

1 1 -12 4 39 1.705882 -10.859353
2 4 39 1.705882 -10.859353 2.205541 -7.093511
3 1.705882 -10.859353 2.205541 -7.093511 3.146719 9.571443
4 2.205541 -7.093511 3.146719 9.571443 2.606157 -1.723461
5 3.146719 9.571443 2.606157 -1.723461 2.688640 -0.318952
6 2.606157 -1.723461 2.688640 -0.318952 2.707372 0.015171
7 2.688640 -0.318952 2.707372 0.015171 2.706521 -0.000123
8 2.707372 0.015171 2.706521 -0.000123 2.706528 -0.000000

The approximated root of x3 − 4x − 9 = 0 using secant method is 2.706.

Homework 1

Use the Secant method to find a root of the equation 3x + sin(x) − ex = 0 accurate to four decimal
places.

Homework 2

Use the Secant method to find a root of the equation xsin(x) + cos(x) = 0 accurate to four decimal
places.

Homework 3

Use the Secant method to find a root of the equation xsin(x) = 1

x accurate to four decimal places.

Prepared By: Er. Narayan Sapkota, [Link]. 17 Everest Engineering College

CHAPTER 1. SOLUTIONS OF NON-LINEAR EQUATIONS

1.3 The Newton-Raphson Method

The Newton-Raphson method is the best-known method of finding roots of a function f (x). The method
is simple and fast. One drawback of this method is that it uses the derivative f ′ (x) of the function as well
as the function f (x). Hence, the Newton-Raphson method is usable only in problems where f ′ (x) can be
computed. Here, again we assume that f (x) is continuous and differentiable and the equation is known to
have a solution near a given point.

Let x0 be an initial approximation to the root of f (x) = 0. Then, the point P (x0 , f (x0 )), where f (x0 ) = f0 ,
lies on the curve. We 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. This 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 , f (x0 )) 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 ) ′
x = x0 − , f (x) ̸= 0
f ′ (x0 )

Figure 1.3: Illustration of Newton-Raphson Method

The next approximation to the root is given by:

f (x0 )
x1 = x0 − , f ′ (x0 ) ̸= 0.
f ′ (x0 )

Again, the next approximation to the root is given by:

f (x1 ) ′
x2 = x1 − , f (x1 ) ̸= 0 (1.6)
f ′ (x1 )

The Equation (1.6) is the Newton-Raphson formula. It is also called the tangent method.

Prepared By: Er. Narayan Sapkota, [Link]. 18 Everest Engineering College

CHAPTER 1. SOLUTIONS OF NON-LINEAR EQUATIONS

We repeat the procedure. The iteration method is defined as:

f (xn )
xn+1 = xn − , f ′ (xn ) ̸= 0 (1.7)
f ′ (xn )

The Newton-Raphson method has second order convergence.

Example 1

Use the Newton-Raphson method to find a root of the equation xx + x − 4 = 0 accurate to three
decimal places.

Solution: Here,
f (x) = xx + x − 4
f ′ (x) = xx (ln(x) + 1) + 1
Suppose the initial guess for x0 = 3.
We know that, general form of Newton-Raphson method is,

f (xn )
xn+1 = xn =
f ′ (xn )

. Therefore,
f (3)
x1 = 3 − = 2.549101
f ′ (3)

Iteration,n xn f (x) f ′ (x) xn+1

1 3.000000 9.411399 22.026616 2.549101
2 2.549101 3.056089 9.646205 2.121827
3 2.121827 0.708684 5.618489 1.805009
4 1.805009 0.065466 4.623144 1.678875
5 1.678875 0.000687 4.526513 1.664715
6 1.664715 0.000000 4.525492 1.664563

The approximated root of xx + x − 4 = 0 is 1.664.

Homework 1

Use the Newton-Raphson method to find a root of the equation xlog10 (x) = 1.2 accurate to four
decimal places.

Homework 2
Use the Newton-Raphson method to find a root of the equation 3x = cosx+1 accurate to four decimal
places.

Homework 3
Use Newton’s method to find the smallest root of the equation ex sin x = 1 to four decimal places.

Prepared By: Er. Narayan Sapkota, [Link]. 19 Everest Engineering College

CHAPTER 1. SOLUTIONS OF NON-LINEAR EQUATIONS

1.4 Iteration Method

The Fixed Point Iteration method is based on the principle of finding a sequence {xn }, each element of which
successively approximates a real root α of an equation. We re-write the equation as
ϕ(x) = 0 so that x = ϕ(x).
Thus, a root of the given equation satisfies
α = ϕ(x)
Therefore, the point x = α remains fixed under the mapping ϕ(x), and so a root of the equation is a fixed
point of ϕ.

The function ϕ(x) is called the iteration function. Here we also assume that ϕ(x) is continuously differentiable
in [a, b]. Let x0 be the initial approximation of r. Thus, x0 satisfies the equation
x0 = g(x0 ).
Putting x0 in the iteration function ϕ, we get the first approximation x1 of α as
x1 = ϕ(x0 ).
Then the successive approximations are calculated as
x2 = ϕ(x1 ), x3 = ϕ(x2 ), ....

The above iteration is generated by the formula

xn+1 = ϕ(xn ) (1.8)

and is called the iteration formula, where xn is the n-th approximation of the root α.
The sequence {xn } of iterations or the successive better approximations may or may not converge to a limit.
If {xn } converges, then it converges to α, and the number of iterations required depends upon the desired
degree of accuracy of the root.
The iteration converges:
lim xn = α
n→∞

if and only if |ϕ′ (x)| < 1 in some neighborhood of α.

We can test this condition using x0 , the initial approximation, before the computations are done.

Remark

Sometimes, it may not be possible to find a suitable iteration function ϕ(x) by manipulating the given
function f (x). In such cases, we may use the following procedure. Write f (x) = 0 as

x = x + αf (x) = ϕ(x),

where α is a constant to be determined. Let x0 be an initial approximation contained in the interval

in which the root lies. For convergence, we require

|ϕ′ (x0 )| = |1 + αf ′ (x0 )| < 1

Simplifying, we find the interval in which α lies. We choose a value for α from this interval and
compute the approximations. A judicious choice of a value in this interval may give faster convergence.

The order of convergence depends on the derivatives of ϕ at α:

Prepared By: Er. Narayan Sapkota, [Link]. 20 Everest Engineering College

CHAPTER 1. SOLUTIONS OF NON-LINEAR EQUATIONS

• Linear convergence if ϕ′ (α) ̸= 0 (and |ϕ′ (α)| < 1)

• Quadratic convergence if ϕ′ (α) = 0 and ϕ′′ (α) ̸= 0
• p-th order convergence if all derivatives up to (p − 1)-th vanish at α

Example 1

Find the root of the equation

f (x) = x3 − 5x + 1 = 0.
by the iteration method, correct to four significant figures.

Solution:
Note: As we do not need the interval where the root lies, but it makes it easier to check the convergence.
The positive root lies in the interval (0, 1).
(i)
x3k + 1
xk+1 = , k = 0, 1, 2, . . .
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 x4 ≈ x5 = 0.20164 is taken as the required approximation to the root.
(ii)
1/3
xk+1 = (5xk − 1) , k = 0, 1, 2, . . .
With x0 = 1, we get the sequence of approximations as

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

which does not converge to the root in (0, 1).

(iii)
5xk − 1
r
xk+1 = , k = 0, 1, 2, . . .
xk
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).

Example 2

Find the root of the equation 3x − cos x − 1 = 0, by the iteration method, correct to four significant
figures.

Solution: Let f (x) = 3x − cos x − 1. As f (0) = −2 < 0, f (1) = 1.45 > 0, there is one root of f (x) = 0
between 0 and 1.

We re-write the equation as

cos x + 1 sin x
x= = ϕ(x), ⇒ ϕ′ (x) = −
3 3

Therefore, |ϕ′ (x)| < 1, as | sin x| < 1 and the convergence criteria is satisfied.

Prepared By: Er. Narayan Sapkota, [Link]. 21 Everest Engineering College

CHAPTER 1. SOLUTIONS OF NON-LINEAR EQUATIONS

n xn ϕ(xn )
0 0 0.6
1 0.6 0.61
2 0.61 0.606
3 0.606 0.6073
4 0.6073 0.60706
5 0.60706 0.60711
6 0.60711 0.60710

As a initial value, we take x0 = 0. The successive approximations of the root are computed in tabular form
as follows:

Remark
Convergence of an iteration method

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

depends on the choice of the iteration function ϕ(x) and a suitable initial approximation x0 to the
root.

Prepared By: Er. Narayan Sapkota, [Link]. 22 Everest Engineering College

CHAPTER 1. SOLUTIONS OF NON-LINEAR EQUATIONS

Assignment 1
Caution 1: log10 x ̸= lne (x). log10 x is commonly used to denote the base 10 logarithm, and lne (x) is used
for the natural logarithm (base e). To convert between the two, one may use:

lne (x)
log10 (x) =
lne (10)

Caution 2: Some functions may have multiple roots; the answer provided is one of them.
1. Use Bisection method to find the root of
(a) x3 − x + 1 = 0
(b) sinx = 1
x

(c) cos x = xex

(d) x log10 x = 1.2
(e) Find the fourth root of 32.
(f) Find the negative root of x3 − x − 11 = 0. [Hint: First find f (−x) then root of f (−x) is the
negative root of f (x)][Answer: -2.3736]
2. Use Secant method to find the root of
(a) x3 − x + 1 = 0
(b) sinx = 1
x

(c) cos x = xex

(d) x log10 x = 1.2
(e) Find the smallest negative root in magnitude of the equation 3x3 − x + 1 = 0 correct to four
decimal places.
(f) Find the smallest positive root of the equation x5 − 64x + 30 = 0 correct to four decimal places.
3. Use Newton Raphson method to find the root of
(a) x3 − x + 1 = 0
(b) sinx = 1
x

(c) cos x = xex

(d) x log10 x = 1.2
(e) xcosx + sinx = 0
(f) Find the positive root, between 0 and 1, of the equation x = e−x to a tolerance of 0.05%.[Answer:
0.567]
4. Use Iteration method to find the root of
(a) x3 − x + 1 = 0
(b) sinx = 1
x

(c) cos x = xex

(d) x log10 x = 1.2
(e) xcosx + sinx = 0

Prepared By: Er. Narayan Sapkota, [Link]. 23 Everest Engineering College

CHAPTER 1. SOLUTIONS OF NON-LINEAR EQUATIONS

Prepared By: Er. Narayan Sapkota, [Link]. 24 Everest Engineering College

 Chapter 2
Interpolation and Approximation

Interpolation refers to the method of predicting the value of a function at any intermediate point along the
independent variable. Extrapolation, on the other hand, involves determining the value of a function for a
point beyond the known range of data.
If (xi , yi ), where i = 0, 1, 2, . . . , represent the (n + 1) given data points of the function y = f (x), then
the process of determining the value of y for any x = xi between x0 and xn is known as interpolation. In
interpolation, we try to construct the polynomial pn (xi ) = yi = f (xi ) ∀i = 0, 1, 2, ..., n. This polynomial
pn (xi ) is called interpolating polynomial.

Interpolation Properties
Interpolation Property: The most fundamental property is that the interpolation polynomial P (x) passes
exactly through all the given data points (xi , yi ) for i = 0, 1, . . . , n. Mathematically, this property is expressed
as P (xi ) = yi for all i.
Degree: The degree of the interpolation polynomial is determined by the number of data points used for
interpolation. If n + 1 points are used, the interpolation polynomial is of degree n or less.
Uniqueness: Under certain conditions (such as when the data points are distinct), the interpolation poly-
nomial is unique for a given set of data points. This means that using the same set of points will always
yield the same interpolation polynomial.

2.1 The Lagrange Interpolation

2.1.1 Introduction
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 (2.1)

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

2.1.2 Linear Interpolation

Linear interpolation is interpolation by the straight line through (x0 , f0 ), (x1 , f1 ); see Fig. 2.1. Thus the
linear Lagrange polynomial p1 is a sum p1 = l0 f0 + l1 f1 with l0 the linear polynomial that is 1 at x0 and 0
at x1 ; similarly, l1 is 0 at x0 and 1 at x1 . Obviously,

x − x1 x − x0
l0 (x) = l1 (x) =
x0 − x 1 x1 − x 0

This gives the linear Lagrange polynomial

25
CHAPTER 2. INTERPOLATION AND APPROXIMATION

x − x1 x − x0
p1 (x) = l0 (x)f0 + l1 (x)f1 = f0 + f1 (2.2)
x0 − x 1 x1 − x 0

Figure 2.1: Linear Interpolation

Find the polynomial f (9.2) for the following data points:

xi 9 9.5
f (xi ) 2.1972 2.2513

Solution:
x − x1 x − 9.5
l0 (x) = = = −2(x − 9.5)
x0 − x 1 9 − 9.5
x − x0 x−9
l1 (x) = = = 2(x − 9)
x1 − x 0 9.5 − 9

l0 (9.2) = 0.6 l1 (9.2) = 0.4

Using formula from (2.2):

x − x1 x − x0
p1 (x) = l0 (x)f0 + l1 (x)f1 = f0 + f1 = 2.2188
x0 − x 1 x1 − x 0

2.1.3 Quadratic Interpolation

Quadratic interpolation is interpolation of given (x0 , f0 ), (x1 , f1 ), (x2 , f2 ) by a second-degree polynomial
p2 (x), which by Lagrange’s idea is

p2 (x) = l0 (x)f0 + l1 (x)f1 + l2 (x)f2

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

Prepared By: Er. Narayan Sapkota, [Link]. 26 Everest Engineering College

CHAPTER 2. INTERPOLATION AND APPROXIMATION

(x − x0 )(x − x2 )
l1 (x) =
(x1 − x0 )(x1 − x2 )
(x − x0 )(x − x1 )
l2 (x) =
(x2 − x0 )(x2 − x1 )

Finally,

(x − x1 )(x − x2 ) (x − x0 )(x − x2 ) (x − x0 )(x − x1 )

p2 (x) = f0 + f1 + f2 (2.3)
(x0 − x1 )(x0 − x2 ) (x1 − x0 )(x1 − x2 ) (x2 − x0 )(x2 − x1 )

Find the polynomial f (9.2) for the following data points:

xi 9 9.5 11
f (xi ) 2.1972 2.2513 2.3972

Solution: The Lagrange basis polynomials and their evaluations at x = 9.2 are:

(x − 9.5)(x − 11.0)
l0 (x) = = x2 − 20.5x + 104.5, L0 (9.2) = 0.54
(9.0 − 9.5)(9.0 − 11.0)

(x − 9.0)(x − 11.0) 1
l1 (x) = =− (x2 − 20x + 99), L1 (9.2) = 0.48
(9.5 − 9.0)(9.5 − 11.0) 0.75

(x − 9.0)(x − 9.5) 1
l2 (x) = = (x2 − 18.5x + 85.5), L2 (9.2) = −0.02
(11.0 − 9.0)(11.0 − 9.5) 3

Using formula from (2.3)

p2 (9.2) = 0.54 2.1972 + 0.48 2.2513 − 0.02 2.3979 = 2.2192

2.1.4 General Lagrange Interpolation Polynomial

Given real values x1 , x2 , x3 , . . . , xn and y1 , y2 , y3 , . . . , yn , there exists a polynomial P with real coefficients
satisfying the conditions P (xi ) = yi ∀i = 0, 1, 2, . . . n, and the degree of polynomial P must be less than
the count of real values, i.e., degree(P ) < n. The Lagrange interpolation polynomial Pn (x) is given by:
n
X
P (x) = yi · li (x)
i=1

Prepared By: Er. Narayan Sapkota, [Link]. 27 Everest Engineering College

CHAPTER 2. INTERPOLATION AND APPROXIMATION

where x1 , x2 , x3 , . . . , xn are distinct real values, and

n
Y x − xj
li (x) =
x
j=1 i
− xj
j̸=i

Proof:
Let’s consider an nth-degree polynomial of the given form,

f (x) = A0 (x − x1 )(x − x2 ) . . . (x − xn ) (2.4)

+ A1 (x − x0 )(x − x2 ) . . . (x − xn )
+ ...
+ An (x − x0 )(x − x1 ) . . . (x − xn−1 )

Substitute observations xi to get Ai :

Put x = x0 to get A0 :
f (x0 ) = y0 = A0 (x0 − x1 )(x0 − x2 ) . . . (x0 − xn )
y0
A0 =
(x0 − x1 )(x0 − x2 ) . . . (x0 − xn )

By substituting x = x1 we get A1 :

f (x1 ) = y1 = A1 (x1 − x0 )(x1 − x2 ) . . . (x1 − xn )

y1
A1 =
(x1 − x0 )(x1 − x2 ) . . . (x1 − xn )

Similarly, by substituting x = xn we get An :

f (xn ) = yn = An (xn − x0 )(xn − x1 ) . . . (xn − xn−1 )

yn
An =
(xn − x0 )(xn − x1 ) . . . (xn − xn−1 )

If we substitute all values of Ai in function f (x) where i = 1, 2, 3, . . . , n then we get the Lagrange interpolation
formula as,

(x − x1 )(x − x2 ) . . . (x − xn )
P (x) = · y0 (2.5)
(x0 − x1 )(x0 − x2 ) . . . (x0 − xn )
(x − x0 )(x − x2 ) . . . (x − xn )
+ · y1
(x1 − x0 )(x1 − x2 ) . . . (x1 − xn )
+ ...
(x − x0 )(x − x1 ) . . . (x − xn−1 )
+ · yn
(xn − x0 )(xn − x1 ) . . . (xn − xn−1 )
= l0 (x) · y0 + l1 (x) · y1 + . . . + ln (x) · yn

n n
X Y x − xj
P (x) = yi · li (x) where li (x) = (2.6)
i=1
x
j=1 i
− xj
j̸=i

Prepared By: Er. Narayan Sapkota, [Link]. 28 Everest Engineering College

CHAPTER 2. INTERPOLATION AND APPROXIMATION

2.1.5 Example
Equal Interval

Find the polynomial f (x) by using Lagrange’s formula and hence find f (3) for following data points:

x 0 1 2 5
f (x) 2 3 12 147

Solution: First, calculate the Lagrange basis polynomials:

(x − 1)(x − 2)(x − 5) (x − 1)(x − 2)(x − 5) x3 − 8x2 + 13x − 10

l0 (x) = = =
(0 − 1)(0 − 2)(0 − 5) (−1)(−2)(−5) −10
(x − 0)(x − 2)(x − 5) (x − 0)(x − 2)(x − 5) −3x3 + 21x2 − 40x + 20
l1 (x) = = =
(1 − 0)(1 − 2)(1 − 5) (1)(−1)(−4) −6
(x − 0)(x − 1)(x − 5) (x − 0)(x − 1)(x − 5) x3 − 6x2 + 5x
l2 (x) = = =
(2 − 0)(2 − 1)(2 − 5) (2)(1)(−3) 6
(x − 0)(x − 1)(x − 2) (x − 0)(x − 1)(x − 2) x3 − 3x2 + 2x
l3 (x) = = =
(5 − 0)(5 − 1)(5 − 2) (5)(4)(3) 18

Now, construct the Lagrange interpolation polynomial P (x):

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

= 2 · l0 (x) + 3 · l1 (x) + 12 · l2 (x) + 147 · l3 (x)
= x3 + x2 − x + 2

Evaluate P (3) to find f (3):

f (3) = P (3) = 33 + 32 − 3 + 2 = 35

Unequal Interval

Find value of log10 (301) by using Lagrange’s interpolation for following data points:

x 300 304 305 307

log10 (x) 2.4771 2.4829 2.4843 2.4871

Solution: Here, x = 301. First, calculate the Lagrange basis polynomials:

(301 − 304)(301 − 305)(301 − 307) 18

l0 (x) = =
(300 − 304)(300 − 305)(300 − 307) 35
(301 − 300)(301 − 305)(301 − 307)
l1 (x) = =2
(304 − 300)(304 − 305)(304 − 307)
(301 − 300)(301 − 304)(301 − 307) −9
l2 (x) = =
(305 − 300)(305 − 304)(305 − 307) 5
(301 − 300)(301 − 304)(301 − 305) 2
l3 (x) = =
(307 − 300)(307 − 304)(307 − 305) 7

Now, construct the Lagrange interpolation polynomial P (x):

P (301) = 2.4771 · l0 (x) + 2.4829 · l1 (x) + 2.4843 · l2 (x) + 2.4871 · l3 (x)

= 2.4786

Prepared By: Er. Narayan Sapkota, [Link]. 29 Everest Engineering College

CHAPTER 2. INTERPOLATION AND APPROXIMATION

2.1.6 Inverse Interpolation

Given a set of values of x and y, the process of finding the value of x for a certain value of y is called inverse
interpolation.

n n
X Y y − yj
P (y) = xi · li (y) where li (y) = (2.7)
i=1
y − yj
j=1 i
j̸=i

Find missing value for following data points:

xi 1 ? 3 4
yi 4 7 14 19

Solution:
y0 = 4, y1 = 14, y2 = 19, y = 7 x0 = 1, x1 = 3, x2 = 4
Using (2.7) for n = 3:
(y − y1 )(y − y2 ) (y − y0 )(y − y2 ) (y − y0 )(y − y1 )
p2 (y) = x0 + x1 + x2
(y0 − y1 )(y0 − y2 ) (y1 − y0 )(y1 − y2 ) (y2 − y0 )(y2 − y1 )
(7 − 14)(7 − 19) (7 − 4)(7 − 19) (7 − 4)(7 − 14)
p2 (7) = 1+ 3+ 4
(4 − 14)(4 − 19) (14 − 4)(14 − 19) (19 − 4)(19 − 14)
= 1.6

2.1.7 Exercise
1. Given the following data points for the function f (x) = sin x + cos x:

x 10◦ 20◦ 30◦

f (x) 1.1585 1.2817 1.3660
Construct the quadratic Lagrange interpolating polynomial that fits the data. Hence, find f ( 12
π
).
Compare with the exact value. [Answer: 1.2247]

2.2 Finite Differences

2.2.1 Forward Differences
If y0 , y1 , y2 , . . . , yn denote a set of values of y, then y1 − y0 , y2 − y1 , . . . , yn − yn−1 are called the differences
of y. Denoting these differences by ∆y0 , ∆y1 , . . . , ∆yn−1 respectively, we have
∆y0 = y1 − y0 , ∆y1 = y2 − y1 , . . . , ∆yn−1 = yn − yn−1 ,
where ∆ is called the forward difference operator and ∆y0 , ∆y1 , . . . , are called first forward differences.
The differences of the first forward differences are called second forward differences and are denoted by
∆2 y0 , ∆2 y1 , . . .. Similarly, one can define third forward differences, fourth forward differences, etc.
Thus,
∆y1 = y1 − y0
∆2 y0 = ∆y1 − ∆y0 = y2 − y1 − (y1 − y0 ) = y2 − 2y1 + y0
∆3 y0 = ∆2 y1 − ∆2 y0 = y3 − 2y2 + y1 − (y2 − 2y1 + y0 ) = y3 − 3y2 + 3y1 − y0
∆4 y0 = ∆3 y1 − ∆3 y0 = y4 − 3y3 + 3y2 − y1 − (y3 − 3y2 + 3y1 − y0 ) = y4 − 4y3 + 6y2 − 4y1 + y0

Prepared By: Er. Narayan Sapkota, [Link]. 30 Everest Engineering College

CHAPTER 2. INTERPOLATION AND APPROXIMATION

Figure 2.2: Forward difference Table

In practical computations, the forward difference table can be formed in the following way. For the data
points (xi , yi ), i = 0, 1, 2, . . . , n and xi = x0 + ih, we have

∆yj = yj+1 − yj , j = 0, 1, . . . , n − 1

In general, nth differences are defined as:

∆n yk = ∆n−1 yk+1 − ∆n−1 yk

2.2.2 Backward Differences

The differences y1 − y0 , y2 − y1 , . . . , yn − yn−1 are called first backward differences if they are denoted by
∇y1 , ∇y2 , . . . , ∇yn respectively, so that
∇y1 = y1 − y0 ,
..
.
∇yn = yn − yn−1 ,
where ∇ is called the backward difference operator. In a similar way, one can define backward differences of
higher orders.
Thus, we obtain
∇y1 = y1 − y0
∇2 y2 = ∇y2 − ∇y1 = y2 − y1 − (y1 − y0 ) = y2 − 2y1 + y0
∇3 y3 = ∇2 y3 − ∇2 y2 = y3 − 3y2 + 3y1 − y0

In general, nth differences are defined as:

∇n yk = ∇n−1 yk+1 − ∇n−1 yk

Prepared By: Er. Narayan Sapkota, [Link]. 31 Everest Engineering College

CHAPTER 2. INTERPOLATION AND APPROXIMATION

Figure 2.3: Backward Difference Table

2.2.3 Central Differences

The central difference operator δ is defined by the relations

y1 − y0 = δy1/2 , y2 − y1 = δy3/2 , . . . , yn − yn−1 = δyn−1/2

Similarly, higher-order central differences can be defined. With the values of x and y as in the preceding two
tables, a central difference can be formed:

Figure 2.4: Central Difference Table

It is clear from all the four tables that in a definite numerical case, the same numbers occur in the same
positions whether we use forward, backward or central differences. Thus, we obtain

∆y0 = ∇y1 = δy1/2 , ∆3 y2 = ∇3 y5 = δ 3 y7/2 , . . .

Prepared By: Er. Narayan Sapkota, [Link]. 32 Everest Engineering College

CHAPTER 2. INTERPOLATION AND APPROXIMATION

2.2.4 Detection of Errors by use of Difference Tables

A difference table is a tabular representation used to compute the differences between successive function
values. It is especially useful when the data is assumed to follow a polynomial pattern. The table is con-
structed by calculating differences between consecutive values, then differences of differences (second, third,
etc.).

The primary use of difference tables is in detecting errors in numerical computations, especially when the
data is expected to follow a known function or pattern.

Given a set of data points (x0 , f (x0 )), (x1 , f (x1 )), . . . , (xn , f (xn )), the difference table is constructed in the
following steps:
1. Step 1: First Differences
The first differences are calculated by subtracting each successive f (x) value from the previous one:

∆f (xi ) = f (xi+1 ) − f (xi )

2. Step 2: Second Differences

The second differences are the differences of the first differences:

∆2 f (xi ) = ∆f (xi+1 ) − ∆f (xi )

3. Step 3: Higher-order Differences

Higher-order differences can be calculated similarly. For third differences, for example:

∆3 f (xi ) = ∆2 f (xi+1 ) − ∆2 f (xi )

2.3 Error Detection in Difference Tables

If the data is consistent and follows a polynomial pattern, the differences should show a predictable, regular
pattern:
• The first differences should increase or decrease linearly for a quadratic function.
• The second differences should be constant for a quadratic function.
• Higher-order differences should be zero for polynomial functions.
If any irregularities or inconsistencies appear in the differences, it is likely that an error exists in the data.
These errors can be due to incorrect data entry, arithmetic mistakes, or rounding errors.

Example

Error Detection in a Quadratic Function

Consider the function f (x) = x2 + 2 and calculate the difference table for the data points from x = 1
to x = 4.
x f (x) = x2 + 2
1 3
2 6
3 11
4 18

Prepared By: Er. Narayan Sapkota, [Link]. 33 Everest Engineering College

CHAPTER 2. INTERPOLATION AND APPROXIMATION

Solution:
The first differences are:
∆f (x0 ) = 6 − 3 = 3

∆f (x1 ) = 11 − 6 = 5

∆f (x2 ) = 18 − 11 = 7

∆f (x) = 3, 5, 7
The second differences are:
∆2 f (x0 ) = 5 − 3 = 2

∆2 f (x1 ) = 7 − 5 = 2

∆2 f (x) = 2, 2
The third differences are:
∆3 f (x0 ) = 2 − 2 = 0

∆3 f (x) = 0
Since the second differences are constant and the third differences are zero, the data follows a quadratic
pattern, and there are no errors.

Error Detection in a Quadratic Function

Let’s consider the function f (x) = x2 + 2 and introduce an error by changing f (3) = 10 instead of
the correct value f (3) = 11.
x f (x)
1 3
2 6
3 10 (Error)
4 18

Solution: The first differences are:

∆f (x0 ) = 6 − 3 = 3

∆f (x1 ) = 10 − 6 = 4 (Error, should be 5)

∆f (x2 ) = 18 − 10 = 8

∆f (x) = 3, 4, 8
The second differences are:
∆2 f (x0 ) = 4 − 3 = 1 (Error, should be 2)

∆2 f (x1 ) = 8 − 4 = 4

∆2 f (x) = 1, 4
The third differences are:
∆3 f (x0 ) = 4 − 1 = 3

∆3 f (x) = 3
The irregularities in the first and second differences indicate that there is an error in the data. Specifically,
the value f (3) = 10 is incorrect. The correct value should be f (3) = 11, which would yield consistent first
and second differences. The difference table has successfully detected the error.

Prepared By: Er. Narayan Sapkota, [Link]. 34 Everest Engineering College

CHAPTER 2. INTERPOLATION AND APPROXIMATION

2.4 Differences of a Polynomial

Let y(x) be a polynomial of the nth degree so that

y(x) = a0 xn + a1 xn−1 + a2 xn−2 + · · · + an

Then we obtain
y(x + h) − y(x) = a0 [(x + h)n − xn ] + a1 [(x + h)n−1 − xn−1 ] + · · ·
= a0 (nh)xn−1 + a′1 xn−2 + · · · + a′n
where a′1 , a′2 , . . . , a′n are the new coefficients.

The above equation can be written as

∆y(x) = a0 (nh)xn−1 + a′1 xn−2 + · · · + a′n

which shows that the first difference of a polynomial of the nth degree is a polynomial of degree (n − 1).
Similarly, the second difference will be a polynomial of degree (n − 2), and the coefficient of xn−2 will be
a0 n(n − 1)h2 .

Thus the nth difference is a0 n!hn , which is a constant. Hence, the (n + 1)th, and higher differences of a
polynomial of nth degree will be zero. Conversely, if the nth differences of a tabulated function are constant
and the (n + 1)th, (n + 2)th, . . ., differences all vanish, then the tabulated function represents a polynomial
of degree n. It should be noted that these results hold good only if the values of x are equally spaced. The
converse is important in numerical analysis since it enables us to approximate a function by a polynomial if
its differences of some order become nearly constant.

2.5 Newton’s Interpolation using Difference

2.5.1 Newton’s Forward Difference Interpolation
Newton’s forward interpolation is a method used to approximate values of a function between known data
points. It’s particularly useful when dealing with equally spaced intervals. It is used to find the value
which is near the begining of the difference table.

The interpolation polynomial in Newton’s forward interpolation formula is given by:

∆2 y 0 ∆3 y0
P (x) = y0 + ∆y0 u + u(u − 1) + u(u − 1)(u − 2) + . . . (2.8)
2! 3!

where
x − x0
u=
h
is the normalized distance.

2.5.2 Newton’s Backward Difference Interpolation

The interpolation polynomial in Newton’s backward interpolation formula is given by:

∇2 y n ∇3 yn
P (x) = yn + ∇yn u + u(u + 1) + u(u + 1)(u + 2) + . . . (2.9)
2! 3!

Prepared By: Er. Narayan Sapkota, [Link]. 35 Everest Engineering College

CHAPTER 2. INTERPOLATION AND APPROXIMATION

where,
x − xn
u=
h
represents the normalized distance. This formula allows us to approximate the value of a function at a given
point x using backward interpolation, which is useful when the point lies near the end of the data range.

2.5.3 Example
Example 1

The table given below gives the value of tan(x) for 0.10 ≤ x ≤ 0.30. Find the value of tan(x) at (i)
x = 0.12 (ii) x = 0.26

xi 0.10 0.15 0.2 0.25 0.30

tan(xi ) 0.1003 0.1511 0.2027 0.2553 0.3093

Solution:

Newton’s forward interpolation (x = 0.12)

Since the value x = 0.12 is near the beginning of the difference table. Therefore, to get the corresponding
value, we use Newton forward interpolation formula.

xi yi ∆y/∇y ∆2 y/∇2 y ∆3 y/∇3 y ∆4 y/∇4 y

0.1 0.1003
0.1511 - 0.1003 = 0.0509
0.15 0.1511 0.0515- 0.0509 = 0.006
0.2027 - 0.1511 = 0.0515 0.0005
0.20 0.2027 0.0526- 0.0515 = 0.0011 -0.0002
0.2553 - 0.2027 = 0.0526 0.0003
0.25 0.2553 0.0540-0.0526= 0.0014
0.3093 - 0.2553 = 0.0540
0.30 0.3093

Since x = 0.12 lies between 0.1 and 0.15, therefore x0 = 0.1 and y0 = 0.1003.

x − x0 0.12 − 0.1
h = 0.15 − 0.1 = 0.05, u= = = 0.4
h 0.05
Now, using newton’s forward interpolation formula,

∆2 y0 ∆3 y 0 ∆4 y0
P (x) = y0 + ∆y0 u + u(u − 1) + u(u − 1)(u − 2) + u(u − 1)(u − 2)(u − 3)
2! 3! 4!
= 0.1205798

Newton’s backward interpolation (x = 0.26)

Since the value x = 0.26 is near the beginning of the difference table, we use Newton backward interpolation
formula to get the corresponding value.
Since x = 0.26 lies between 0.3 and 0.25, therefore xn = 0.3 and yn = 0.3093.

x − xn 0.26 − 0.30
h = 0.30 − 0.25 = 0.05, u= = = −0.8
h 0.05

Prepared By: Er. Narayan Sapkota, [Link]. 36 Everest Engineering College

CHAPTER 2. INTERPOLATION AND APPROXIMATION

Using Newton’s backward interpolation formula,

u(u + 1) 2 u(u + 1)(u + 2) 3 u(u + 1)(u + 2)(u + 3) 4

P (x) = yn + u∆yn + ∆ yn + ∆ yn + ∆ yn
2! 3! 4!
= 0.26602

Example 2

Using Newton’s forward difference formula, find the sum

Sn = 13 + 23 + 33 + · · · + n3

Solution: We have
Sn+1 = 13 + 23 + 33 + · · · + n3 + (n + 1)3
Hence
Sn+1 − Sn = (n + 1)3

or
∆Sn = (n + 1)3

It follows that

∆2 Sn = ∆Sn+1 − ∆Sn = (n + 2)3 − (n + 1)3 = 3n2 + 9n + 7,

∆3 Sn = 3(n + 1)2 + 9(n + 1) + 7 − (3n2 + 9n + 7) = 6n + 12,
∆4 Sn = 6(n + 1) + 12 − (6n + 12) = 6.

Since ∆5 Sn = ∆6 Sn = · · · = 0, Sn is a fourth-degree polynomial in n. Further,

S1 = 1, ∆S1 = 8, ∆2 S1 = 19, ∆3 S1 = 18, ∆4 S1 = 6.

(3.1) gives

(n − 1)(n − 2) (n − 1)(n − 2)(n − 3)

Sn = 1 + (n − 1)(8) + (19) + (18)
2 6
(n − 1)(n − 2)(n − 3)(n − 4)
+ (6)
24
1 1 1
= n4 + n3 + n2
4 2 4
2
n(n + 1)

=
2

2.6 Newton’s Divided Difference Interpolation

2.6.1 Introduction
The application of the Lagrange formula becomes quite lengthy and difficult if the number of values in the
data is more than four. It can be seen that the inclusion of a point in the data list leads to the addition of
one more term in the formula, and the numerator and denominator of each term get changed.

Further, during calculations, there is always a chance of committing some error due to a number of positive
and negative signs in the coefficients of each term. Except for this, there is no way to guess the degree of

Prepared By: Er. Narayan Sapkota, [Link]. 37 Everest Engineering College

CHAPTER 2. INTERPOLATION AND APPROXIMATION

the resulted polynomial before the final result is obtained after lengthy calculations.

For example, suppose the data corresponds to a polynomial of degree 2, but there are 5 data points. We will
have to simplify five terms, with each coefficient being a polynomial of degree 4. Only after simplifying these
five terms can we obtain the final polynomial of degree 2. Hence, it is not advisable to use the Lagrange
formula for more than four points in the data set.

To overcome this problem, we prefer to use the Newton’s Divided Difference formula.

2.6.2 Newton’s Forward Interpolation for Divided Differences

Using Newton’s Divided Difference Method, the interpolating polynomial is given by:

P (x) =f [x0 ] + f [x0 , x1 ] (x − x0 )

+ f [x0 , x1 , x2 ] (x − x0 )(x − x1 )
+ f [x0 , x1 , x2 , x3 ] (x − x0 )(x − x1 )(x − x2 )
+ f [x0 , x1 , x2 , x3 , x4 ] (x − x0 )(x − x1 )(x − x2 )(x − x3 ) + . . . (2.10)

where the divided differences are computed as follows:

f [xi ] = yi
f [xj ] − f [xi ]
f [xi , xj ] =
xj − x i
f [xj , xk ] − f [xi , xj ]
f [xi , xj , xk ] = and so on.
xk − x i

2.6.3 Newton’s Backward Interpolation for Divided Differences

Using Newton’s Divided Difference Method, the interpolating polynomial is given by:

P (x) = f [xn ] + f [xn , xn−1 ] (x − xn )

+ f [xn , xn−1 , xn−2 ] (x − xn )(x − xn−1 )
+ f [xn , xn−1 , xn−2 , xn−3 ] (x − xn )(x − xn−1 )(x − xn−2 ) + . . . (2.11)

2.6.4 Example

Evaluate f (9), using Newton’s divided difference formula from the following data-points.

x 5 7 11 13 17
f (x) 150 392 1452 2366 5202

Solution:

Newton’s forward Divided Difference

The difference table for the given data-points is tabulated below.

Prepared By: Er. Narayan Sapkota, [Link]. 38 Everest Engineering College

CHAPTER 2. INTERPOLATION AND APPROXIMATION

xi yi f [x0 , x1 ] f [x0 , x1 , x2 ] f [x0 , x1 , x2 , x3 ] f [x0 , x1 , x2 , x3 , x4 ]

x0 = 5 150
392−150
7−5 = 121
x1 = 7 392 265−121
11−5 = 24
1452−392
11−7 = 265 32−24
15−5 =1
x2 = 11 1452 457−265
13−7 = 32 0
2366−1452
13−11 = 457 17−7 =
42−32
1
x3 = 13 2366 709−457
17−11 = 42
5202−2366
17−13 = 709
x4 = 17 5202

Using formula for Newton’s Forward Interpolation for Divided Differences:

P (x) =f [x0 ] + f [x0 , x1 ] (x − x0 )

+ f [x0 , x1 , x2 ] (x − x0 )(x − x1 )
+ f [x0 , x1 , x2 , x3 ] (x − x0 )(x − x1 )(x − x2 )
+ f [x0 , x1 , x2 , x3 , x4 ] (x − x0 )(x − x1 )(x − x2 )(x − x3 ) + . . .
= 150 + 121 × (x − 5) + 24 × (x − 5)(x − 7) + 1 × (x − 5)(x − 7)(x − 11)
+ 0 × (x − 5)(x − 7)(x − 11)(x − 13)
= 150 + 121 × (x − 5) + 24 × (x − 5)(x − 7) + 1 × (x − 5)(x − 7)(x − 11)

Therefore,

P (9) = f (9) = 150 + (9 − 5) × 121 + (9 − 5)(9 − 7) × 24

+ (9 − 5)(9 − 7)(9 − 11) × 1
= 150 + 484 + 192 − 16
= 810

Newton’s Backward Divided Difference

xi yi f [x4 , x3 ] f [x4 , x3 , x2 ] f [x4 , x3 , x2 , x1 ] f [x4 , x3 , x2 , x1 , x0 ]

x0 = 5 150
392−150
7−5 = 121
x1 = 7 392 265−121
11−5 = 24
1452−392
11−7 = 265 32−24
15−5 =1
x2 = 11 1452 457−265
13−7 = 32 0
2366−1452
13−11 = 457 17−7 =
42−32
1
x3 = 13 2366 709−457
17−11 = 42
5202−2366
17−13 = 709
x4 = 17 5202

Using formula for Newton’s Backward Interpolation for Divided Differences:

P (x) = f [xn ] + f [xn , xn−1 ] (x − xn ) + f [xn , xn−1 , xn−2 ] (x − xn )(x − xn−1 )

+ f [xn , xn−1 , xn−2 , xn−3 ] (x − xn )(x − xn−1 )(x − xn−2 ) + . . .
= yn + ∆y (x − xn ) + ∆2 y (x − xn )(x − xn−1 )
+ ∆3 y (x − xn )(x − xn−1 )(x − xn−2 ) + . . .

Prepared By: Er. Narayan Sapkota, [Link]. 39 Everest Engineering College

CHAPTER 2. INTERPOLATION AND APPROXIMATION

Therefore,

P (9) = f (9) = 5202 + (9 − 17) × 121 + (9 − 17)(9 − 137) × 24

+ (9 − 17)(9 − 13)(9 − 11) × 1 + (9 − 17)(9 − 13)(9 − 11)(9 − 7) × 0
= 5202 − 5672 + 1344 − 64
= 810

2.7 Curve Fitting: Least Squares Lines for Linear and Nonlinear
Data
2.8 Cubic Spline Interpolation

Prepared By: Er. Narayan Sapkota, [Link]. 40 Everest Engineering College

 Chapter 3
Numerical Differentiation and Integration

3.1 Numerical Differentiation

Numerical differentiation is the technique of estimating the derivative of a function at a given value of x, based
on a set of known data points (xi , yi ). Rather than using the exact function y = f (x), we approximate it with
an interpolating polynomial y = ϕ(x), and then differentiate this polynomial as needed. The appropriate
interpolation formula depends on the location of the point where the derivative is required:

If the values of x are equally spaced and the derivative is needed near the start of the data,
Newton’s forward formula is used. If the derivative is needed near the end, Newton’s
backward formula is applied. When the required point lies near the middle of the table,
Stirling’s or Bessel’s formula is more suitable. If the x-values are not equally spaced, we
use Lagrange’s interpolation formula or Newton’s divided difference formula. Each of
these interpolation methods can be used to derive a corresponding differentiation formula.

Important Note: It’s essential to recognize that the data table only defines the function at specific points
and may not represent the complete function, which might not be differentiable everywhere. Therefore,
numerical differentiation should be applied only when the data indicates that differences of a certain order
remain approximately constant. Otherwise, errors can grow significantly, especially when computing higher-
order derivatives. This occurs because while the interpolating polynomial ϕ(x) might closely match the actual
function f (x) at the given points, the derivative ϕ′ (x) may still differ considerably from f ′ (x).

3.2 Newton’s Differentiation Formulas

3.2.1 Derivatives Using Newton’s Forward Difference Interpolation Formula
Consider the function y = f (x), tabulated at equidistant points xi = x0 + ih, i = 0, 1, 2, . . . , n. Let

x − x0
u=
h

The Newton forward difference interpolation formula is:

∆2 y 0 ∆3 y 0
y(x) = y0 + ∆y0 u + u(u − 1) + u(u − 1)(u − 2)
2! 3!
∆4 y 0
+ u(u − 1)(u − 2)(u − 3)
4! (3.1)
∆5 y 0
+ u(u − 1)(u − 2)(u − 3)(u − 4)
5!
∆ y0
6
+ u(u − 1)(u − 2)(u − 3)(u − 4)(u − 5) + . . .
6!

First Derivative
Differentiating both sides of Equation (3.1) with respect to u:

41
CHAPTER 3. NUMERICAL DIFFERENTIATION AND INTEGRATION

dy ∆2 y 0
= ∆y0 + [(u − 1) + u]
du 2!
∆3 y0
+ [(u − 1)(u − 2) + u(u − 2) + u(u − 1)]
3!
∆4 y0 d
+ [u(u − 1)(u − 2)(u − 3)] (3.2)
4! du
∆ y0 d
5
+ [u(u − 1)(u − 2)(u − 3)(u − 4)]
5! du
∆6 y0 d
+ [u(u − 1)(u − 2)(u − 3)(u − 4)(u − 5)]
6! du

Since u = h ,
x−x0
we have du
dx = h1 . Therefore, by the chain rule:

dy dy du 1 dy
= · = ·
dx du dx h du

Hence, the first derivative of y with respect to x becomes:

"
dy 1 ∆2 y 0 ∆3 y0  2
= ∆y0 + [2u − 1] + 3u − 6u + 2

dx h 2! 3!
∆4 y 0
+ 4u − 18u2 + 22u − 6
 3 
4! (3.3)
∆5 y 0
+ 5u − 40u3 + 105u2 − 100u + 24
 4 
5! #
∆6 y 0
+ 6u − 75u4 + 300u3 − 540u2 + 390u − 120 + · · ·
 5 
6!

At x = x0 , u = 0, so:

1 ∆2 y0 ∆3 y 0 ∆4 y0 ∆5 y 0 ∆6 y0
 
dy
= ∆y0 − + − + − + ··· (3.4)
dx x=x0 h 2 3 4 5 6

Second Derivative
Differentiating Equation (3.2) w.r.t. u:

d2 y 6u − 18u + 11
 2 
= ∆2
y 0 + ∆ 3
y0 (u − 1) + ∆4
y 0
du2 12
2u − 12u + 22u − 10
 3 2

+ ∆ 5 y0 (3.5)
12
30u − 300u3 + 1020u2 − 1350u + 548
4
 
+ ∆ 6 y0 + ···
720

Recalling that
d2 y 1 d2 y
= ,
dx2 h2 du2

Prepared By: Er. Narayan Sapkota, [Link]. 42 Everest Engineering College

CHAPTER 3. NUMERICAL DIFFERENTIATION AND INTEGRATION

we obtain the second derivative with respect to x:

"
d2 y 1 6u − 18u + 11
 2 
= ∆ 2
y 0 + ∆ 3
y0 (u − 1) + ∆4
y0
dx2 h2 12
2u − 12u + 22u − 10
 3 2

+ ∆5 y0 (3.6)
12
#
30u4
300u 3
+ 1020u 2
1350u + 548
 
− −
+ ∆6 y0 + ...
720

At x = x0 (u = 0), this simplifies to:

d2 y 1 11 5 137 6
 
= 2 ∆ y0 − ∆ y0 + ∆4 y0 − ∆5 y0 +
2 3
∆ y0 + · · · (3.7)
dx2 x=x0 h 12 6 180

3.2.2 Derivatives Using Newton’s Backward Difference Interpolation Formula

Consider the function y = f (x), tabulated at equidistant points xi = x0 + ih, i = 0, 1, 2, . . . , n. Let

x − xn
u=
h

The Newton backward difference interpolation formula is:

∇2 y n ∇3 yn
y(x) = yn + ∇yn u + u(u + 1) + u(u + 1)(u + 2)
2! 3!
∇4 yn
+ u(u + 1)(u + 2)(u + 3)
4! (3.8)
5
∇ yn
+ u(u + 1)(u + 2)(u + 3)(u + 4)
5!
∇6 yn
+ u(u + 1)(u + 2)(u + 3)(u + 4)(u + 5) + · · ·
6!

First Derivative

Differentiating Equation (3.8) with respect to u:

dy ∇2 y n ∇3 y n
= ∇yn + (2u + 1) + (3u2 + 6u + 2)
du 2! 3!
∇4 yn
+ (4u3 + 18u2 + 22u + 6)
4!
∇5 yn
+ (5u4 + 40u3 + 105u2 + 100u + 24)
5!
∇6 yn
+ (6u5 + 75u4 + 340u3 + 675u2 + 548u + 120) + · · ·
6!

Prepared By: Er. Narayan Sapkota, [Link]. 43 Everest Engineering College

CHAPTER 3. NUMERICAL DIFFERENTIATION AND INTEGRATION

Using du
dx = h1 , we get:
"
dy 1 2u + 1 2 3u2 + 6u + 2 3
= ∇yn + ∇ yn + ∇ yn
dx h 2! 3!
4u3 + 18u2 + 22u + 6 4
+ ∇ yn
4! (3.9)
5u4 + 40u3 + 105u2 + 100u + 24 5
+ ∇ yn
5! #
6u5 + 75u4 + 340u3 + 675u2 + 548u + 120 6
+ ∇ yn + · · ·
6!

At x = xn , we have u = 0, so:

1 1 2 1 3 1 4 1 5 1 6
 
dy
= ∇yn + ∇ yn + ∇ yn + ∇ yn + ∇ yn + ∇ yn + · · · (3.10)
dx x=xn h 2 3 4 5 6

Second Derivative
Differentiating Equation (3.9) with respect to u:

d2 y ∇2 y n ∇3 yn ∇4 yn
= ·2+ (6u + 6) + (12u2 + 36u + 22)
du 2 2! 3! 4!
∇5 y n
+ (20u3 + 120u2 + 210u + 100)
5!
∇6 y n
+ (30u4 + 300u3 + 1020u2 + 1350u + 548) + · · ·
6!

Using du
dx = h1 , we get:
"
d2 y 1 12u2 + 36u + 22 4
= 2 ∇2 yn + (u + 1)∇3 yn + ∇ yn
dx 2 h 24
20u3 + 120u2 + 210u + 100 5
+ ∇ yn (3.11)
60 #
30u4 + 300u3 + 1020u2 + 1350u + 548 6
+ ∇ yn + · · ·
720

At x = xn (u = 0):

d2 y 1 11 4 5 5 137 6
 
= ∇ 2
y n + ∇ 3
yn + ∇ yn + ∇ yn + ∇ yn + · · · (3.12)
dx2 x=xn h2 12 3 180

Example 1

Given
x 1.0
1.1 1.2 1.3 1.4 1.5 1.6
y 7.989
8.403 8.781 9.129 9.451 9.750 10.031
   2 
dy d y
Using the above table, compute the and at x = 1.1 and x = 1.6.
dx dx2

Prepared By: Er. Narayan Sapkota, [Link]. 44 Everest Engineering College

CHAPTER 3. NUMERICAL DIFFERENTIATION AND INTEGRATION

Solution: Given that the values of x are equidistant with step size:h = 1.1 − 1.0 = 0.1
(a) At x = 1.1 we, use Newton’s Forward Difference Formula

x y ∆y ∆2 y ∆3 y ∆4 y ∆5 y ∆6 y
1.0 7.989
0.414
1.1 8.403 −0.036
0.378 0.006
1.2 8.781 −0.030 −0.002
0.348 0.004 0.001
1.3 9.129 −0.026 −0.001 0.002
0.322 0.003 0.003
1.4 9.451 −0.023 0.002
0.299 0.005
1.5 9.750 −0.018
0.281
1.6 10.031

Table 3.1: Difference Table for given data

Using Newton’s Forward Difference formula for the first derivative:

"
dy 1 ∆2 y0 ∆3 y0  2
= ∆y0 + [2u − 1] + 3u − 6u + 2

dx h 2! 3!
∆4 y 0
+ 4u − 18u2 + 22u − 6
 3 
4! (3.13)
∆5 y 0
+ 5u − 40u3 + 105u2 − 100u + 24
 4 
5! #
∆6 y 0
+ 6u − 75u4 + 300u3 − 540u2 + 390u − 120
 5 
6!

For x = 1.1, u = x−x0

h = 1.1−1.0
0.1 =1

dy
∴ = 3.947
dx
Now, the second derivative is:
"
d2 y 1 6u − 18u + 11
 2 
= ∆ y0 + ∆ y0 (u − 1) + ∆ y0
2 3 4
dx2 h2 12
2u − 12u2 + 22u − 10
 3 
+ ∆5 y0
12
#
30u4 − 300u3 + 1020u2 − 1350u + 548

+ ∆ y0
6
720

d2 y
∴ = −3.583
dx2
(b) At x = 1.6 we use Newton’s Backward Difference Formula
We use the above difference table and the backward difference operator ∇ instead of ∆.
For x = 1.6, u = x−xn
h = 1.6−1.6
0.1 =0

Prepared By: Er. Narayan Sapkota, [Link]. 45 Everest Engineering College

CHAPTER 3. NUMERICAL DIFFERENTIATION AND INTEGRATION

Using Newton’s Backward Difference formula for the first derivative:

"
dy 1 2u + 1 2 3u2 + 6u + 2 3
= ∇yn + ∇ yn + ∇ yn
dx h 2! 3!
4u3 + 18u2 + 22u + 6 4
+ ∇ yn
4!
5u4 + 40u3 + 105u2 + 100u + 24 5
+ ∇ yn
5!
6u5 + 75u4 + 340u3 + 675u2 + 548u + 120 6
+ ∇ yn +]
6!

dy
∴ = 2.75
dx
Now, the second derivative:

d2 y 1 11 4 5 5 137 6
 
= 2 ∇ yn + ∇ yn + ∇ yn + ∇ yn +
2 3
∇ yn + · · ·
dx2 x=xn h 12 3 180

d2 y
∴ = −0.715
dx2

Homework 1
Given time and velocity in the table below. Find the initial acceleration.
t (sec) 0 5 10 15 20
v (m/s) 0 3 14 69 228

Homework 2

Find the value of cos(1.74) from the following table:

x 1.7 1.74 1.78 1.82 1.86
sin x 0.9916 0.9857 0.9781 0.9691 0.9584

Prepared By: Er. Narayan Sapkota, [Link]. 46 Everest Engineering College

CHAPTER 3. NUMERICAL DIFFERENTIATION AND INTEGRATION

3.3 Numerical Integration

The process of evaluating a definite integral from a set of tabulated values of the integrand f (x) is called
numerical integration. This process, when applied to a function of a single variable, is known as quadrature.
The problem of numerical integration, like that of numerical differentiation, is solved by representing f (x)
by an interpolation formula and then integrating it between the given limits. In this way, we can derive
quadrature formulae for approximate integration of a function defined by a set of numerical values only.

3.4 Newton-Cotes Quadrature Formula

Let Z b
I= f (x) dx
a

where f (x) takes the values y0 , y1 , y2 , . . . , yn for x = x0 , x1 , x2 , . . . , xn .

Divide the interval (a, b) into n sub-intervals of width h, so that x0 = a, x1 = x0 + h, x2 = x0 + 2h, . . . , xn =
x0 + nh = b.

We want to evaluate the definite integral:

Z x0 +nh
I= f (x) dx
x0

Let
x = x0 + rh ⇒ dx = h dr
As x goes from x0 to x0 + nh, r goes from 0 to n. So,
Z x0 +nh Z n
I= f (x) dx = h f (x0 + rh) dr
x0 0

We approximate f (x0 + rh) using Newton’s forward interpolation:

r(r − 1) 2 r(r − 1)(r − 2) 3

f (x0 + rh) = y0 + r∆y0 + ∆ y0 + ∆ y0 + · · ·
2! 3!

Substituting into the integral:

r(r − 1) 2 r(r − 1)(r − 2) 3

Z n 
I=h y0 + r∆y0 + ∆ y0 + ∆ y0 + · · · dr
0 2! 3!

Prepared By: Er. Narayan Sapkota, [Link]. 47 Everest Engineering College

CHAPTER 3. NUMERICAL DIFFERENTIATION AND INTEGRATION

Each term is a polynomial in r, so we integrate term by term.

Z n
y0 dr = y0 · n
0

n
n2
Z
r∆y0 dr = ∆y0 ·
0 2
r(r − 1) 2
n
∆2 y 0 n 2
Z Z
∆ y0 dr = (r − r) dr
0 2 2 0

∆2 y0 n3 n2
 
= −
2 3 2
n
r(r − 1)(r − 2) 3 ∆3 y0 n 3
Z Z
∆ y0 dr = (r − 3r2 + 2r) dr
0 6 6 0

∆3 y 0 n 4
 
= −n +n
3 2
6 4
Putting all terms together:

n2 1 n3 n2 1 n4
     
I = h ny0 + ∆y0 + − ∆2 y 0 + − n 3 + n 2 ∆3 y 0 + · · ·
2 2 3 2 6 4
Factor out nh:
n(2n − 3) 2 n(n − 2)(2n − 1) 3
 
n
I = nh y0 + ∆y0 + ∆ y0 + ∆ y0 + · · ·
2 12 72

This is the general Newton–Cotes quadrature formula, from which various rules such as the trapezoidal
rule, Simpson’s 1/3 rule, Simpson’s 3/8 rule, and others are derived by choosing specific values of n and
neglecting higher-order differences.

3.4.1 Trapezoidal Rule

For n = 1, using first-order polynomial:

Z x6
h
f (x) dx ≈ [(y0 + y6 ) + 2(y1 + y2 + y3 + y4 + y5 )] (3.14)
x0 2

General form: Z x0 +nh

h
f (x) dx ≈ [y0 + 2(y1 + y2 + · · · + yn−1 ) + yn ]
x0 2

3.4.2 Simpson’s One-Third Rule

For n = 2, using a parabola (second-order polynomial):

Z x6
h
f (x) dx ≈ [(y0 + y6 ) + 4(y1 + y3 + y5 ) + 2(y2 + y4 )] (3.15)
x0 3

General form (n even):

Z x0 +nh
h
f (x) dx ≈ [y0 + yn + 4(y1 + y3 + · · · ) + 2(y2 + y4 + · · · )]
x0 3

Prepared By: Er. Narayan Sapkota, [Link]. 48 Everest Engineering College

CHAPTER 3. NUMERICAL DIFFERENTIATION AND INTEGRATION

3.4.3 Simpson’s Three-Eighth Rule

For n = 3, using a cubic polynomial:

x6
3h
Z
f (x) dx ≈ [(y0 + y6 ) + 3(y1 + y2 + y4 + y5 ) + 2y3 ] (3.16)
x0 8

General form:
 
n−1 n−3
xn
3h 
Z X X
f (x) dx ≈ y0 + yn + 3 yi + 2 yi 

x0 8 i=1 i=3
i mod 3̸=0 i mod 3=0

Example 1
R6 1
Evaluate 0 x+1 dx using
1. Trapezoidal rule
1
2. Simpson’s rule
3
3
3. Simpson’s rule
8
4. Compare the results with its actual value.

Solution: Let h = 1, and define f (x) = x+1 .

1
Then the interval [0, 6] is divided into n = 6 sub-intervals.

The nodes are:

x0 = 0, x1 = 1, x2 = 2, x3 = 3, x4 = 4, x5 = 5, x6 = 6

Function values:
1
y0 = f (0) = =1
1
1
y1 = f (1) =
2
1
y2 = f (2) =
3
1
y3 = f (3) =
4
1
y4 = f (4) =
5
1
y5 = f (5) =
6
1
y6 = f (6) =
7

(1) Trapezoidal Rule

Z 6
h
f (x) dx ≈ [y0 + 2(y1 + y2 + y3 + y4 + y5 ) + y6 ]
0 2

1 1 1 1 1 1 1
   
≈ 1+2 + + + + + ≈ 2.0214
2 2 3 4 5 6 7

Prepared By: Er. Narayan Sapkota, [Link]. 49 Everest Engineering College

CHAPTER 3. NUMERICAL DIFFERENTIATION AND INTEGRATION

1
(2) Simpson’s Rule
3
Z 6
h
f (x) dx ≈ [y0 + 4(y1 + y3 + y5 ) + 2(y2 + y4 ) + y6 ]
0 3

1 1 1 1 1 1 1
     
= 1+4 + + +2 + + ≈ 1.9587
3 2 4 6 3 5 7

3
(3) Simpson’s Rule
8
6
3h
Z
f (x) dx ≈ [y0 + 3(y1 + y2 + y4 + y5 ) + 2y3 + y6 ]
0 8

3 1 1 1 1 1 1
     
= 1+3 + + + +2 + ≈ 1.9661
8 2 3 5 6 4 7

(4) Actual Value

6
1
Z 6
dx = ln(x + 1) = ln(7) − ln(1) = ln(7) ≈ 1.9459
0 x+1 0

Method Absolute Error

Trapezoidal Rule
Simpson’s 31 Rule
Simpson’s 83 Rule

Prepared By: Er. Narayan Sapkota, [Link]. 50 Everest Engineering College

CHAPTER 3. NUMERICAL DIFFERENTIATION AND INTEGRATION

3.5 Romberg’s Method

We have derived approximate quadrature formulae with the help of finite differences method. Romberg’s
method provides a simple modification to these quadrature formulae for finding their better approximations.
As an illustration, let us improve upon the value of the integral
Z b
I= f (x) dx,
a

by the Trapezoidal rule. If I1 , I2 are the values of I with sub-intervals of width h1 , h2 and E1 , E2 their
corresponding errors, respectively, then

(b − a)h21 ′′ (b − a)2 h22 ′′ ¯

E1 = − f (X̄), E2 = − f (X̄ )
12 12

¯ ) are very
Since f ′′ (X̄) is also the largest value of f ′′ (x), we can reasonably assume that f ′′ (X̄) and f ′′ (X̄
nearly equal.

E1 h2 E1 h2
∴ = 12 or = 2 1 2 (1)
E2 h2 E2 − E1 h2 − h1

Now since I = I1 + E1 = I2 + E2 ,
⇒ E2 − E1 = I1 − I2 (2)

From (1) and (2), we have:

h21
E1 = (I1 − I2 )
h22 − h21

Hence,
h21 I1 h22 − I2 h21
I = I1 + E1 = I1 + (I1 − I2 ) i.e., I= (3)
h22 − h21 h22 − h21
which is a better approximation of I.
To evaluate I systematically, we take h1 = h and h2 = h
2 so that (3) gives

I1 (h/2)2 − I2 h2 4I2 − I1
I= =
(h/2) − h
2 2 3

1
i.e., I(h, h/2) = [4I(h/2) − I(h)] (4)
3

Now we use the trapezoidal rule several times, successively halving h and apply (4) to each pair of values as
per the following scheme:

I(h)
I(h, h/2)
I(h/2) I(h, h/2, h/4)
I(h/2, h/4) I(h, h/2, h/4, h/8)
I(h/4) I(h/2, h/4, h/8)
I(h/4, h/8)
I(h/8)

The computation is continued until successive values are close to each other. This method is called Richard-
son’s deferred approach to the limit and its systematic refinement is called Romberg’s method.

Prepared By: Er. Narayan Sapkota, [Link]. 51 Everest Engineering College

CHAPTER 3. NUMERICAL DIFFERENTIATION AND INTEGRATION

Example
R6
Evaluate 1
0 x+1
dx using Romberg Integration.

Solution:
When n=8:
b−a 6−0
a = 0, b = 6, n = 8, h= = = 0.75
n 8

i 0 1 2 3 4 5 6 7 8
xi 0 0.75 1.5 2.25 3.0 3.75 4.5 5.25 6.0
yi = xi1+1 1 0.571429 0.4 0.307692 0.25 0.210526 0.181818 0.16 0.142857

Z 6
h
I(h/8) = f (x) dx ≈ [y0 + 2(y1 + y2 + y3 + y4 + y5 + y6 + y7 ) + y8 ] ≈ 1.9897
0 2
When n=4:
b−a 6−0
a = 0, b = 6, n = 4, h= = = 1.5
n 4

i 0 1 2 3 4
xi 0 1.5 3.0 4.5 6.0
yi = xi1+1 1 0.4 0.25 0.181818 0.142857

Z 6
h
I(h/4) = f (x) dx ≈ [y0 + 2(y1 + y2 + y3 ) + y4 ] ≈ 2.1049
0 2

When n=2:
b−a 6−0
a = 0, b = 6, n = 2, h= = =3
n 2

i 0 1 2
xi 0 3.0 6.0
yi = xi1+1 1 0.25 0.142857

Z 6
h
I(h/2) = f (x) dx ≈ [y0 + 2(y1 ) + y2 ] ≈ 2.4643
0 2
When n=1:
b−a 6−0
a = 0, b = 6, n = 1, h= = =6
n 1

i 0 1
xi 0 6.0
yi = xi1+1 1 0.142857

Z 6
h
I(h/1) = f (x) dx ≈ [y0 + y1 ] =≈ 3.4286
0 2

Prepared By: Er. Narayan Sapkota, [Link]. 52 Everest Engineering College

CHAPTER 3. NUMERICAL DIFFERENTIATION AND INTEGRATION

As we know,
1
I(h, h/2) = [4I(h/2) − I(h)]
3
1
I(h/2, h/4) = [4I(h/4) − I(h/2)]
3
1
I(h/4, h/8) = [4I(h/8) − I(h/4)]
3
1
I(h, h/2, h/4) = [4I(h/2, h/4) − I(h, h/2)]
3
1
I(h/2, h/4, h/8) = [4I(h/4, h/8) − I(h/2, h/4)]
3
1
I(h, h/2, h/4, h/8) = [4I(h/2, h/4, h/8) − I(h, h/2, h/4)]
3
Therefore,

I(h) = 3.4286
I(h, h/2) = 2.1429
I(h/2) = 2.4643 I(h, h/2, h/4) = 1.9325
I(h/2, h/4) = 1.9851 I(h, h/2, h/4, h/8) = 1.9425
I(h/4) = 2.1049 I(h/2, h/4, h/8) = 1.9400
I(h/4, h/8) = 1.9513
I(h/8) = 1.9897

R6
∴ 1
0 x+1
dx = 1.9425

Prepared By: Er. Narayan Sapkota, [Link]. 53 Everest Engineering College

CHAPTER 3. NUMERICAL DIFFERENTIATION AND INTEGRATION

3.6 Gaussian Integration algorithm

Gaussian integration, also known as Gaussian quadrature, is a numerical method for approximating the
definite integral of a function. This method is highly efficient and accurate for polynomials and functions
that can be well-approximated by polynomials over the integration interval.
The integral
Z b
I= f (x) dx
a
is approximated by
n
X
I≈ wi × f (xi )
i=1

where xi are the nodes (or abscissas) and wi are the corresponding weights. These nodes and weights are
chosen such that the approximation is exact for polynomials of degree 2n − 1 or less. The general idea of
Gaussian quadrature is to approximate the integral of a function f (x) over an interval [a, b] by a weighted
sum of the function’s values at specified points within the interval.

While only defined for the interval [−1, 1], this is actually a universal function, because we can convert the
limits of integration for any interval [a, b] to the Legendre-Gauss interval [−1, 1]:

Z b Z 1 n
X
f (x) dx = f (ξi ) dξ ≈ wi × f (ξi ) (3.17)
a −1 i=1

where, x = b−a
2 ξ + b+a
2

Number of points, n Points, ξi Weights, wi

1 0 2
2 ± √13 = ±0.57735 1
3 q 0 8
9 = 0.889
± = ±0.774597
3
5 9 =
5
0.556
r q √
4 ± 37 − 27 65 = ±0.339981 18+ 30
36 = 0.652
r q √
± 37 + 27 65 = ±0.86113 18− 30
36 = 0.347
5 r 0 225 = 0.568889
128
q √
± 13 5 − 2 10
7 = ±0.538469
322+13 70
900 = 0.4786
r q √
± 13 5+2 10
7 = ±0.90618 322−13 70
900 = 0.23692

Example

Given the integral Z 1.3

5xe−2x dx
0.1

a) Use the one-point Gauss quadrature rule to estimate the value of the integral.
b) Use the two-point Gauss quadrature rule to estimate the value of the integral.
c) Use the three-point Gauss quadrature rule to estimate the value of the integral.
d) Find the true error for part (a), (b) and (c).
e) Find the absolute relative true error for part (a), (b) and (c).

Prepared By: Er. Narayan Sapkota, [Link]. 54 Everest Engineering College

CHAPTER 3. NUMERICAL DIFFERENTIATION AND INTEGRATION

Solution: The Gaussian Quadrature nodes and weights are defined on the interval [−1, 1]. We change
variables from x ∈ [0.1, 1.3] to ξ ∈ [−1, 1] via

b−a a+b
x= ξ+ = 0.6ξ + 0.7,
2 2
where a = 0.1, b = 1.3.
Then, the integral becomes
1.3 1 1
b−a
Z Z Z
I= 5xe−2x dx = f (ξ) dξ = 0.6 f (ξ) dξ,
0.1 2 −1 −1

where
f (ξ) = 5 (0.6ξ + 0.7) e−2(0.6ξ+0.7) .

One-point Gauss Quadrature (n = 1)

The one-point Gauss quadrature has
ξ1 = 0, w1 = 2.

The estimate is
I1 = 0.6 × w1 × f (ξ1 ) = 1.2 × f (0).

Calculate f (0):

f (0) = 5 × 0.7 × e−2×0.7 = 5 × 0.7 × e−1.4 ≈ 5 × 0.7 × 0.2466 = 0.8631.

Therefore,
I1 ≈ 1.2 × 0.8631 = 1.0357.

Two-point Gauss Quadrature (n = 2)

The two-point Gauss quadrature has nodes and weights

1
ξ1,2 = ± √ ≈ ±0.57735, w1 = w2 = 1.
3

The estimate is
I2 = 0.6 × [f (−0.57735) + f (0.57735)] .

Calculate

f (−0.57735) = 5 × (0.6 × −0.57735 + 0.7) × e−2(0.6×−0.57735+0.7) = 5 × 0.3464 × e−0.6928 ≈ 0.8676,

f (0.57735) = 5 × (0.6 × 0.57735 + 0.7) × e−2(0.6×0.57735+0.7) = 5 × 1.0536 × e−2.1072 ≈ 0.6401.

f (−0.57735) + f (0.57735) = 0.8676 + 0.6401 = 1.5077.

Therefore,
I2 = 0.6 × 1.5077 = 0.9046.

Prepared By: Er. Narayan Sapkota, [Link]. 55 Everest Engineering College

CHAPTER 3. NUMERICAL DIFFERENTIATION AND INTEGRATION

Three-point Gauss Quadrature (n = 3)

The three-point Gauss quadrature has
ξ = {0, ±0.7746},
with weights
8 5
w0 = ≈ 0.8889, w± = ≈ 0.5556.
9 9
The estimate is
I3 = 0.6 × [w0 f (0) + w+ f (0.7746) + w− f (−0.7746)] .

Calculate
f (0) = 0.8631,
f (0.7746) = 5 × (0.6 × 0.7746 + 0.7) × e−2(0.6×0.7746+0.7) = 5 × 1.1648 × e−2.3296 ≈ 0.5687,
f (−0.7746) = 5 × (0.6 × −0.7746 + 0.7) × e−2(0.6×−0.7746+0.7) = 5 × 0.2332 × e−0.4664 ≈ 0.731.

Sum weighted terms:

0.8889 × 0.8631 + 0.5556 × 0.5687 + 0.5556 × 0.731 = 0.767 + 0.316 + 0.406 = 1.489.

Therefore,
I3 = 0.6 × 1.489 = 0.893.

True value and errors

Using a numerical tool, the true value is
Z 1.3
Itrue = 5xe−2x dx ≈ 0.8922.
0.1

Calculate the true errors:

Error = |Itrue − In |.

n In Error = |Itrue − In |
1 1.0357 0.1435
2 0.9046 0.0124
3 0.8930 0.0008

Calculate the absolute relative true error in percentage:

|Itrue − In |
ϵrel = × 100%.
|Itrue |

n ϵrel (%)
1 16.07%
2 1.39%
3 0.09%

Prepared By: Er. Narayan Sapkota, [Link]. 56 Everest Engineering College

CHAPTER 3. NUMERICAL DIFFERENTIATION AND INTEGRATION

Assignment 3
1. Given the following table of values of x and y:

x 1.00 1.05 1.10 1.15 1.20 1.25 1.30

y 1.000 1.025 1.049 1.072 1.095 1.118 1.140

dy d2 y
Use finite differences to approximate the first derivative and the second derivative at:
dx dx2
(a) x = 1.05
(b) x = 1.25
(c) x = 1.15
2. From the following table of values of x and y, find the values of

dy d2 y
and
dx dx2
at x = 2.03.
x 1.96 1.98 2.00 2.02 2.04
y 0.7825 0.7739 0.7651 0.7563 0.7473

3. The table below gives the velocity v of a body at various times t. Using finite differences, find the
acceleration at t = 1.1.

t 1.0 1.1 1.2 1.3 1.4

v 43.1 47.7 52.1 56.4 60.8
R 1.4
4. Evaluate 0.2
(sin x + log x − ex ) dx using
(a) Trapezoidal rule
1
(b) Simpson’s rule
3
3
(c) Simpson’s rule
8
(d) Compare the results with its actual value.
R 0.6
5. Evaluate 0 x2 e−x dx using
(a) Trapezoidal rule
1
(b) Simpson’s rule
3
3
(c) Simpson’s rule
8
(d) Compare the results with its actual value.
6. Solve Question no. 4 and 5 using Romberg Integration method.
7. Solve Question no. 4 and 5 using Gaussian Integration method.(use n= 2 and n=3).

Prepared By: Er. Narayan Sapkota, [Link]. 57 Everest Engineering College

CHAPTER 3. NUMERICAL DIFFERENTIATION AND INTEGRATION

Prepared By: Er. Narayan Sapkota, [Link]. 58 Everest Engineering College

 Chapter 4
Solution of Linear Algebraic Equations

4.1 Matrices and Their Properties

4.1.1 Matrix
A matrix is a rectangular array of elements (numbers or variables) arranged in rows and columns. If a matrix
has m rows and n columns, its order is m × n.

1 2
 
A=
3 4
Here, A is a matrix of order 2 × 2.

4.1.2 Types of Matrices

Row Matrix (Row Vector)
A matrix consisting of a single row is called a row matrix.

A= 2 4 6
 

Column Matrix (Column Vector)

A matrix consisting of a single column is called a column matrix.

1
 

B = 3
5

Submatrix
A matrix obtained by deleting one or more rows or columns from a given matrix is called a submatrix.

1 2
 
A= ⇒ 1
 
3 4

Square Matrix
A matrix having the same number of rows and columns is called a square matrix.

2 1
 
A=
0 3

Upper Triangular Matrix

A square matrix in which all elements below the principal diagonal are zero is called an upper triangular
matrix.
1 2 3
 

A = 0 4 5
0 0 6

59
CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

Lower Triangular Matrix

A square matrix in which all elements above the principal diagonal are zero is called a lower triangular
matrix.
1 0 0
 

A = 2 3 0
4 5 6

Diagonal Matrix
A square matrix in which all non-diagonal elements are zero is called a diagonal matrix.

2 0 0
 

A = 0 5 0
0 0 7

Identity Matrix
An identity matrix is a diagonal matrix in which all diagonal elements are equal to one.

1 0
 
I=
0 1

For any square matrix A, AI = IA = A.

Symmetric and Skew-Symmetric Matrices

A square matrix [aij ] is said to be symmetric if

aij = aji for all i and j.

A square matrix [aij ] is said to be skew-symmetric if

aij = −aji for all i and j,

so that all the leading diagonal elements are zero.

A symmetric matrix:  
a h g
h b f
g f c

A skew-symmetric matrix:
0
 
h g
−h 0 f
−g −f 0

4.1.3 Matrix Operations

Addition of Matrices
Two matrices can be added only if they have the same order. The addition is performed element-wise.

1 2 5 6
   
A= , B=
3 4 7 8

6 8
 
A+B =
10 12

Prepared By: Er. Narayan Sapkota, [Link]. 60 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

Subtraction of Matrices
Subtraction of matrices is also defined only for matrices of the same order.

 
−4 −4
A−B =
−4 −4

Multiplication of Matrices
Let A be a matrix of order m × n and B be a matrix of order n × p. Then the product AB is a matrix of
order m × p.

1 2 2 0
   
A= , B=
3 4 1 2

(1)(2) + (2)(1) (1)(0) + (2)(2) 4 4

   
AB = =
(3)(2) + (4)(1) (3)(0) + (4)(2) 10 8

Transpose of a Matrix
The matrix obtained from a given matrix A by interchanging its rows and columns is called the transpose
of A and is denoted by AT .
Let
1 2 3
 
A= .
4 5 6

Then the transpose of A is

1 4
 

AT = 2 5 .
3 6

Inverse of a Matrix
If A is a non-singular square matrix of order n, then a square matrix B of the same order such that

AB = BA = I,

is called the inverse of A, where I is the unit matrix.

The inverse of A is denoted by A−1 .
Let
1 2
 
A= .
3 4

The determinant of A is
|A| = (1)(4) − (2)(3) = −2 ̸= 0,

so A is non-singular.
The inverse of A is
1 4 1
   
−2 −2
A−1 = = 3 .
|A| −3 1 2 − 12

Prepared By: Er. Narayan Sapkota, [Link]. 61 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

Rank of a Matrix
If we select any r rows and r columns from any matrix A, deleting all other rows and columns, then the
determinant formed by these r × r elements is called a minor of A of order r. Clearly, there will be a number
of different minors of the same order, obtained by deleting different rows and columns from the same matrix.
Definition: A matrix is said to be of rank r if
1. it has at least one non-zero minor of order r, and
2. every minor of order higher than r is zero.

4.1.4 Positive Definite Matrix

A real symmetric matrix A ∈ Rn×n is called positive definite if

xT Ax > 0 for all x ∈ Rn , x ̸= 0.

Sylvester’s Criterion: A real symmetric matrix is positive definite if and only if all its leading principal
minors are positive.
For a 3 × 3 matrix  
a11 a12 a13
A = a12 a22 a23  ,
a13 a23 a33
the conditions are:
1. det([a11 ]) > 0
 
a11 a12
2. det >0
a12 a22
3. det(A) > 0

Example
Consider the matrix
2 0
 
−1
A = −1 2 −1 .
0 −1 2

Step 1: Symmetry

AT = A,
so the matrix is symmetric.
Step 2: First Leading Principal Minor

det([2]) = 2 > 0.

Step 3: Second Leading Principal Minor

2
 
−1
det = (2)(2) − (−1)(−1) = 3 > 0.
−1 2

Prepared By: Er. Narayan Sapkota, [Link]. 62 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

Step 4: Third Leading Principal Minor

2 −1 0
det(A) = −1 2 −1
0 −1 2
Expanding along the first row:
2 −1 −1 −1
=2 − (−1)
−1 2 0 2

= 2(3) + 1(−2) = 4 > 0.

Since all leading principal minors are positive, the matrix A is positive definite.

4.2 Solution of Linear Simultaneous Equations

Simultaneous linear equations occur quite often in engineering and science. The analysis of electronic circuits
consisting of invariant elements, the analysis of a network under sinusoidal steady-state conditions, the
determination of the output of a chemical plant, and the finding of the cost of chemical reactions are some
of the exercises which depend on the solution of simultaneous linear algebraic equations.
The solution of such equations can be obtained by direct iterative methods. We describe below some of these
methods of solution.

4.2.1 Gauss Elimination Method

In this method, the unknowns are eliminated successively and the system is reduced to an upper triangular
system, from which the unknowns are found by back substitution. The method is quite general and is well-
adapted for computer operations. Here we explain it by considering a system of three equations for the sake
of clarity.
Consider the system of equations:
a1 x + b1 y + c1 z = d1 ,
a2 x + b2 y + c2 z = d2 ,
a3 x + b3 y + c3 z = d3 .

Step I: Elimination of x from the second and third equations

a2
Assuming a1 ̸= 0, we eliminate x from the second equation by subtracting times the first equation from
a1
a3
the second equation. Similarly, we eliminate x from the third equation by subtracting times the first
a1
equation from the third equation. Thus, we obtain the new system:
a1 x + b1 y + c1 z = d1 ,
b′2 y + c′2 z = d′2 ,
b′3 y + c′3 z = d′3 .

Here, the first equation is called the pivotal equation and a1 is called the first pivot.
Step II: Elimination of y from the third equation
b′3
Assuming b′2 ̸= 0, we eliminate y from the third equation by subtracting times the second equation from
b′2
the third equation. We then obtain:
a1 x + b1 y + c1 z = d1 ,
b′2 y + c′2 z = d′2 ,
c′′3 z = d′′3 .

Prepared By: Er. Narayan Sapkota, [Link]. 63 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

Here, the second equation is the pivotal equation and b′2 is the new pivot.
Step III: Evaluation of the unknowns
The values of x, y, and z are obtained from the reduced system by back substitution.

Example 1

Apply the Gauss elimination method to solve the equations

x + 4y − z = −5 (i)
x + y − 6z = −12 (ii)
3x − y − z = 4 (iii)

Solution:
Step I: Elimination of x
Subtract equation (i) from equation (ii), and subtract 3 times equation (i) from equation (iii).

(ii) − (i) : −3y − 5z = −7 (iv)

(iii) − 3(i) : −13y + 2z = 19 (v)

Step II: Elimination of y

Eliminate y from equations (iv) and (v) by operating

13
(v) − (iv) :
3

71 148
z= (vi)
3 3

Step III: Back Substitution

From equation (vi),
148
z= = 2.0845
71

From equation (iv),

−3y − 5z = −7
81
y=− = −1.1408
71

From equation (i),

x + 4y − z = −5
117
x= = 1.6479
71

Hence, the solution is

x = 1.6479, y = −1.1408, z = 2.0845

Prepared By: Er. Narayan Sapkota, [Link]. 64 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

Example 1

Using the Gauss elimination method, solve the equations

x + 2y + 3z − u = 10,
2x + 3y − 3z − u = 1,
2x − y + 2z + 3u = 7,
3x + 2y − 4z + 3u = 2.

Solution:
The system can be written in matrix form as
1 2 3 10
    
−1 x
2 3 −3 −1 y  =  1  .
   
2 −1 2 3  z   7 


3 2 −4 3 u 2

Step I: Apply the row operations

R2 → R2 − 2R1 , R3 → R3 − 2R1 , R4 → R4 − 3R1 .

We obtain
1 2 3 −1 | 10
 
0 −1 −9 1 | −19.
0 5 | −13

−5 −4
0 −4 −13 6 | −28

Step II: Apply the row operations

R3 → R3 − 5R2 , R4 → R4 − 4R2 .

This gives
1 2 3 −1 | 10
 
0 −1 −9 1 | −19.
0 0 41 0 | 82 


0 0 23 2 | 48

Step III: Back Substitution

From the third row:
41z = 82 ⇒ z = 2.

From the fourth row:

23z + 2u = 48 ⇒ 46 + 2u = 48 ⇒ u = 1.

From the second row:

−y − 9z + u = −19 ⇒ −y − 18 + 1 = −19 ⇒ y = 2.

From the first row:

x + 2y + 3z − u = 10 ⇒ x + 4 + 6 − 1 = 10 ⇒ x = 1.

Hence, the solution is

x = 1, y = 2, z = 2, u = 1.

Prepared By: Er. Narayan Sapkota, [Link]. 65 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

4.2.2 Partial Pivoting

Pivoting is a technique used in numerical linear algebra to improve the stability of Gaussian elimination. It
involves rearranging the rows and/or columns of a matrix to avoid division by very small numbers, which
can lead to large numerical errors.
Definition: Partial pivoting involves interchanging rows only so that the largest absolute value in the
current column (from the current row downward) is used as the pivot element.
Why: This reduces numerical errors caused by dividing by very small numbers. At step k, the pivot element
akk is chosen such that
|akk | = max |aik |
i=k,...,n

This means: "Choose the element with the largest absolute value in column k, from row k downward, as the
pivot."

Example

Solve the system: (

0.0003x + 1.567y = 1.568
1.234x + 2.345y = 3.579

Solution:
Step 1: Write the augmented matrix

0.0003 1.567 1.568

 

1.234 2.345 3.579

Step 2: Find the pivot in column 1

Look at absolute values in column 1: |0.0003|, |1.234| → largest is 1.234 (row 2).
Swap row 1 and row 2:
1.234 2.345 3.579
 
0.0003 1.567 1.568

Step 3: Proceed with Gaussian elimination

Eliminate the element below the pivot:
0.0003
R2 → R2 − R1
1.234

1.234 2.345 3.579

 
0 1.56643 1.56713

Then back-substitute to find y and x. This is partial pivoting because we only swapped rows.

Example

Solve the system with partial pivoting:

x+y+z =9
2x − 3y + 4z = 13
3x + 4y + 5z = 40

Prepared By: Er. Narayan Sapkota, [Link]. 66 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

Solution:
Step 1: Form the augmented matrix

1 1 1 9
 
 2 −3 4 13 
3 4 5 40

Step 2: Pivoting (rows only)

- Column 1 pivot: choose largest absolute value in column 1.
- Max value is 3 in row 3 =⇒ swap row 1 and row 3:

3 4 5 40
 
 2 −3 4 13 
1 1 1 9

Step 3: Eliminate below pivot

2 1
R2 → R2 − R1 , R3 → R3 − R1
3 3

3 4 5 40
 
 0 − 17 2
− 86 
3 3 3
0 − 13 − 32 7
3

Step 4: Pivot column 2

- Max absolute value in column 2 (rows 2-3) = 17
3 in row 2 (no swap needed)
Step 5: Eliminate below pivot

−1/3 1
R3 → R3 − R2 = R3 − R2
−17/3 17

3 4 5 40
 
 0 − 17 2
− 86 
3 3 3
0 0 − 36
17
145
17

Step 6: Back-substitution

145/17 145
z= =−
−36/17 36

−86/3 − (2/3)z −86/3 − (2/3)(−145/36)

y= = =5
−17/3 −17/3

40 − 4y − 5z
x= =3
3

145
x = 3, y = 5, z = −
36

Prepared By: Er. Narayan Sapkota, [Link]. 67 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

4.2.3 Complete/Full Pivoting

Complete (or full) pivoting involves interchanging both rows and columns so that the largest absolute value
in the remaining submatrix is chosen as the pivot.
At step k, the pivot element apq is chosen such that

|apq | = max |aij |

i=k,...,n
j=k,...,n

This means: "Choose the element with the largest absolute value in the remaining submatrix as the pivot."
Why: It gives better numerical stability than partial pivoting, especially for ill-conditioned matrices.

Example

Consider the matrix:

0 2 3 1
   

A = 4 5 6 , b = 2
7 8 9 3

Solution:
Step 1: Find the largest absolute value in the matrix
The largest element is 9 (row 3, column 3).
Step 2: Swap row 1 ↔ row 3 and column 1 ↔ column 3
New augmented matrix:
9 8 7 3
 
 6 5 4 2 
3 2 0 1

Now proceed with Gaussian elimination as usual.

Step 1: Form augmented matrix

1 1 1 9
 
 2 −3 4 13 
3 4 5 40

Step 2: Choose largest element in entire matrix

- Max absolute value = 5 at row 3, column 3 - Swap row 1 ↔ row 3, column 1 ↔ column 3

5 4 3 40
 
 4 −3 2 13 
1 1 1 9

Step 3: Proceed with Gaussian elimination (eliminate below pivot)

4 1
R2 → R2 − R1 , R3 → R3 − R1
5 5

5 4 3 40
 
 0 − 31 − 25 −19 
5
0 1
5
2
5 1

Prepared By: Er. Narayan Sapkota, [Link]. 68 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

- Continue elimination and back-substitution as usual

Complete pivoting gives same solution (after correcting column order):

145
x = 3, y = 5, z = −
36

4.2.4 Gauss-Jordan Method

The Gauss–Jordan method is an extension of the Gauss elimination method. In this method, the unknowns
are eliminated successively so as to reduce the given system to a diagonal system. The values of the unknowns
are then obtained directly without the need for back substitution. The method is systematic and well-suited
for computer implementation. For the sake of clarity, we explain the method by considering a system of
three equations.
Consider the system of equations:
a1 x + b1 y + c1 z = d1 ,
a2 x + b2 y + c2 z = d2 ,
a3 x + b3 y + c3 z = d3 .

The augmented matrix corresponding to the above system is

 
a1 b1 c1 d1
 a2 b2 c2 d2  .
a3 b3 c3 d3

Step I: Elimination of x from the second and third equations

a2
Assuming a1 ̸= 0, we eliminate x from the second and third equations by subtracting times the first
a1
a3
equation from the second equation and times the first equation from the third equation, respectively.
a1
Thus, we obtain the system:
a1 x + b1 y + c1 z = d1 ,
b′2 y + c′2 z = d′2 ,
b′3 y + c′3 z = d′3 .

Here, the first equation is called the pivotal equation and a1 is called the first pivot.
Step II: Elimination of y from the first and third equations
b1
Assuming b′2 ̸= 0, we eliminate y from the first and third equations by subtracting times the second
b′2
b′3
equation from the first equation and times the second equation from the third equation. We then obtain:
b′2

a′1 x + c′1 z = d′1 ,

b′2 y + c′2 z = d′2 ,
c′′3 z = d′′3 .

Here, the second equation is the pivotal equation and b′2 is the new pivot.
Step III: Elimination of z from the first and second equations

Prepared By: Er. Narayan Sapkota, [Link]. 69 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

c′1
Assuming c′′3 ̸= 0, we eliminate z from the first and second equations by subtracting times the third
c′′3
c′2
equation from the first equation and times the third equation from the second equation. This yields:
c′′3

a′′1 x = d′′1 ,
b′′2 y = d′′2 ,
c′′3 z = d′′3 .

Here, the third equation is the pivotal equation and c′′3 is the final pivot.
Step IV: Evaluation of the unknowns

a′′1 0 0 d′′1
 
 0 b′′2 0 d′′2  .
0 0 c′′3 d′′3

Dividing each row by its pivot, we obtain:

1 0 0
 
x
 0 1 0 y .
0 0 1 z

The values of x, y, and z are obtained directly from the diagonal system by simple division. The augmented
matrix corresponding to the given system is

1 1 1 9
 
 2 −3 4 13  .
3 4 5 40

Step I: Make the first pivot equal to 1 and eliminate x

R2 → R2 − 2R1 , R3 → R3 − 3R1

1 1 1 9
 
 0 −5 2 −5  .
0 1 2 13

Step II: Make the second pivot equal to 1

1
R2 → − R2
5

1 1 1 9
 
 0 1 −2 1 .
5
0 1 2 13

Step III: Eliminate y from the first and third rows

R1 → R1 − R2 , R3 → R3 − R2

Prepared By: Er. Narayan Sapkota, [Link]. 70 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

1 0 7
8
 
5
 0 1 − 25 1 .
0 0 12
5 12

Step IV: Make the third pivot equal to 1

5
R3 → R3
12

1 0 7
8
 
5
 0 1 −2 1 .
5
0 0 1 5

Step V: Eliminate z from the first and second rows

7 2
R1 → R1 − R3 , R2 → R2 + R3
5 5

1 0 0 1
 
 0 1 0 3 .
0 0 1 5

x = 1, y = 3, z = 5.

4.3 Method of Factorization

The method of factorization is a numerical technique used to solve a system of linear equations of the
form
AX = B

by decomposing the coefficient matrix A into a product of simpler matrices. This approach reduces compu-
tational effort and improves numerical stability compared to direct methods.

4.3.1 LU Decomposition
In LU decomposition, the matrix A is factorized as

A = LU

where L is a lower triangular matrix and U is an upper triangular matrix.

The solution procedure consists of:

1. Solving LY = B using forward substitution

2. Solving U X = Y using backward substitution

Prepared By: Er. Narayan Sapkota, [Link]. 71 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

4.3.2 Doolittle Method

In the Doolittle method, the matrix A is factorized such that the lower triangular matrix L has unit
diagonal elements.

1 0 0
   
u11 u12 u13
L = l21 1 0 , U = 0 u22 u23 
l31 l32 1 0 0 u33

The elements of L and U are computed using the following formulae:

i−1
X
uij = aij − lik ukj , i≤j
k=1

j−1
!
1 X
lij = aij − lik ukj , i>j
ujj
k=1

Algorithm:
1. Set diagonal elements of L equal to 1
2. Compute elements of U
3. Compute sub-diagonal elements of L
4. Solve LY = B
5. Solve U X = Y

Example

Apply the Do Little method to solve the equations:


3x + 2y + 7z = 4

2x + 3y + z = 5
3x + 4y + z = 7



Solution: Given the system of linear equations:

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

This can be written in matrix form as:

AX = B
where
3 2 7 4
     
x
A = 2 3 1 , X = y  , B = 5
3 4 1 z 7

Step 1: LU Factorization (Doolittle Method)

In Doolittle method,
A = LU

Prepared By: Er. Narayan Sapkota, [Link]. 72 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

where L has unit diagonal elements.

Let
1 0 0
   
u11 u12 u13
L = l21 1 0 , U = 0 u22 u23 
l31 l32 1 0 0 u33

Equating A = LU :

3 2 7 1 0 0
      
u11 u12 u13 u11 u12 u13
2 3 1 = l21 1 0  0 u22 u23  = l21 u11 l21 u12 + u22 l21 u13 + u23 
3 4 1 l31 l32 1 0 0 u33 l31 u11 l31 u12 + l32 u22 l31 u13 + l32 u23 + u33

Step 2: Determine Elements of U and L

From the first row:
u11 = 3, u12 = 2, u13 = 7

From the second row:

2
l21 =
3
2 5
u22 = 3 − (2) =
3 3
2 11
u23 = 1 − (7) = −
3 3
From the third row:
3
l31 = =1
3
4 − (1)(2) 6
l32 = =
5
3
5
6 11 3
  
u33 = 1 − (1)(7) + − =
5 3 5

Step 3: Write L and U

1 0 0 3 2 7
   

L =  23 1 0 , U = 0 5
3 − 11
3

1 6
5 1 0 0 3
5

Step 4: Solve LY = B
Let  
y1
Y = y2 
y3
Then,
1 0 0 y1 4
    

LY = B ⇒  32 1 0 y2  = 5
1 6
5 1 y3 7
y1 = 4
2 7
y1 + y 2 = 5 ⇒ y 2 =
3 3
6 1
y1 + y2 + y3 = 7 ⇒ y3 =
5 5

Prepared By: Er. Narayan Sapkota, [Link]. 73 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

Thus,
4
 

Y =  73 
1
5

Step 5: Solve U X = Y

3 2 7 4
    
x
U X = Y ⇒ 0 5
3 − 11
3
 y  =  7 
3
0 0 3
5 z 1
5

3 1 1
z= ⇒z=
5 5 3
5 11 7 6
y− z= ⇒y=
3 3 3 5
4
3x + 2y + 7z = 4 ⇒ x =
15

4 6 1
x= , y= , z=
15 5 3

4.3.3 Crout’s Method

In Crout’s method, the upper triangular matrix U has unit diagonal elements.

0 0 1
   
l11 u12 u13
L = l21 l22 0 , U = 0 1 u23 
l31 l32 l33 0 0 1

The elements are computed using:

j−1
X
lij = aij − lik ukj , i≥j
k=1

i−1
!
1 X
uij = aij − lik ukj , i<j
lii
k=1

Algorithm:

1. Set diagonal elements of U equal to 1

2. Compute columns of L

3. Compute rows of U

4. Solve LY = B

5. Solve U X = Y

Prepared By: Er. Narayan Sapkota, [Link]. 74 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

Example

Apply the Crout’s method to solve the equations:


3x + 2y + 7z = 4

2x + 3y + z = 5
3x + 4y + z = 7



Solution:
The given system of equations is 
3x + 2y + 7z = 4

2x + 3y + z = 5
3x + 4y + z = 7



This can be written in matrix form as

AX = B
where
3 2 7 4
     
x
A = 2 3 1 , X = y  , B = 5
3 4 1 z 7

Step 1: Crout’s Factorization

In Crout’s method,
A = LU
where L is a lower triangular matrix and U is an upper triangular matrix with unit diagonal elements.
Let
0 0 1
   
l11 u12 u13
L = l21 l22 0 , U = 0 1 u23 
l31 l32 l33 0 0 1

Then
0 0 1
    
l11 u12 u13 l11 l11 u12 l11 u13
A = LU = l21 l22 0  0 1 u23  = l21 l21 u12 + l22 l21 u13 + l22 u23 
l31 l32 l33 0 0 1 l31 l31 u12 + l32 l31 u13 + l32 u23 + l33

Step 2: Determine the elements of L and U

From the first column:
l11 = 3, l21 = 2, l31 = 3

From the first row:

2 7
u12 = , u13 =
3 3
From the second column:
2 5
 
l22 = 3 − 2 =
3 3
2
 
l32 =4−3 =2
3

From the second row:

1−2 7
11


u23 = 3
=−
5
3
5

Prepared By: Er. Narayan Sapkota, [Link]. 75 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

From the third column:

7 11 3
    
l33 =1− 3 +2 − =
3 5 5

Step 3: Write the matrices L and U

3 0 0 1 2 7
   
3 3
L = 2 5
3 0 , U = 0 1 − 11
5

3 2 3
5 0 0 1

Step 4: Solve LY = B (Forward Substitution)

Let  
y1
Y = y2 
y3

Then,
3 0 0 4
    
y1
LY = B ⇒ 2 5
3 0  y2  = 5
3 2 3
5 y3 7

4
3y1 = 4 ⇒ y1 =
3
5 7
2y1 + y2 = 5 ⇒ y2 =
3 5
3 1
3y1 + 2y2 + y3 = 7 ⇒ y3 =
5 3

Thus,
4
3
Y =  75 
1
3

Step 5: Solve U X = Y (Backward Substitution)

1
 2 7
  4
3 3 x 3
U X = Y ⇒ 0 1 − 11
5
 y  =  75 
0 0 1 z 1
3

1
z=
3
11 7 6
y− z= ⇒y=
5 5 5
2 7 4 4
x+ y+ z = ⇒x=
3 3 3 15

4 6 1
x= , y= , z=
15 5 3

Prepared By: Er. Narayan Sapkota, [Link]. 76 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

4.3.4 Cholesky’s Method

Cholesky’s method is applicable when the matrix A is symmetric and positive definite. The factorization
is given by

A = LLT

where L is a lower triangular matrix and LT is its transpose.

0 0
 
l11
L = l21 l22 0
l31 l32 l33

The elements of L are computed as:

v
u
u i−1
X
lii = taii − 2
lik
k=1

j−1
!
1 X
lij = aij − lik ljk , i>j
ljj
k=1

Algorithm:
1. Verify that A is symmetric and positive definite
2. Compute the matrix L
3. Solve LY = B
4. Solve LT X = Y

Example

Apply Cholesky’s method to solve the system


4x + 2y + 2z = 6

2x + 5y + 4z = 3
2x + 4y + 6z = 10



Solution:
Step 1: Matrix Form

AX = B
where
4 2 2 6
     
x
A = 2 5 4 , X = y  , B =3
2 4 6 z 10

Step 2: LU Factorization with U = LT

Assume
A = LU, U = LT

Prepared By: Er. Narayan Sapkota, [Link]. 77 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

with
0 0
   
l11 l11 l21 l31
L = l21 l22 0 , U =0 l22 l32 
l31 l32 l33 0 0 l33

Then
2
4 2 2
   
l11 l11 l21 l11 l31
LU = l21 l11 2
l21 + l222
l21 l31 + l22 l32  = 2 5 4
l31 l11 l31 l21 + l32 l22 2
l31 + l32
2
+ l33
2
2 4 6

Step 3: Solve for L elements step by step

From the (1,1) entry:
2
l11 = 4 ⇒ l11 = 2

From the (2,1) and (3,1) entries:

l21 l11 = 2 ⇒ l21 = 1, l31 l11 = 2 ⇒ l31 = 1

From the (2,2) entry:

2
l21 + l22
2
= 5 ⇒ 1 + l22
2
= 5 ⇒ l22 = 2

From the (3,2) entry:

3
l31 l21 + l32 l22 = 4 ⇒ 1 ∗ 1 + l32 ∗ 2 = 4 ⇒ l32 =
2

From the (3,3) entry:

 2 √
3 11 11
2
l31 + 2
l32 + 2
l33 =6⇒1+ + l33 = 6 ⇒ l33 =
2 2
⇒ l33 =
2 4 2

Thus,
2 0 0 2 1 1
   

L = 1 2 0 , U = LT = 0 2
3 
√ √2
1 3
2
11
2 0 0 11
2

Step 4: Solve LY = B (Forward Substitution)

2 0 0 6
    
y1
LY = B ⇒ 1 2 √
0  y2  =  3 
1 3
2 2
11 y3 10

2y1 = 6 ⇒ y1 = 3
y1 + 2y2 = 3 ⇒ y2 = 0
√
3 11 14
y1 + y2 + y3 = 10 ⇒ y3 = √
2 2 11

3
 

Y = 0 
√14
11

Step 5: Solve U X = Y (Backward Substitution)

Prepared By: Er. Narayan Sapkota, [Link]. 78 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

2 1 1 3
    
x
UX = Y ⇒ 0 2 = 0 
3  

√2 y
0 0 11 z √14
2 11

√
11 14 28
z= √ ⇒z=
2 11 11
3 21
2y + z = 0 ⇒ y = −
2 11
13
2x + y + z = 3 ⇒ x =
11

Final Solution:

13 21 28
x= , y=− , z=
11 11 11

Example

Apply the Cholesky’s method to solve the equations:


3x + 2y + 7z = 4

2x + 3y + z = 5
3x + 4y + z = 7



Solution: Remark: Cholesky’s method is applicable only to matrices that are symmetric and positive
definite. The coefficient matrix of the given system is

3 2 7
 

A = 2 3 1
3 4 1

Since A ̸= AT , Cholesky decomposition cannot be applied directly. Hence, we form the normal equations:

AT AX = AT B

Step 1: Form the Normal Equations

The given system is not symmetric, so we form the normal equations:

AT AX = AT B

3 2 3
 

AT = 2 3 4
7 1 1

22 24 26 43
   

AT A = 24 29 21 , AT B = 51

26 21 51 40

So the normal equations become:

Prepared By: Er. Narayan Sapkota, [Link]. 79 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

22 24 26 43
    
x
24 29 21 y  = 51
26 21 51 z 40

Step 2: LU Factorization (Crout / Doolittle style)

We set AT A = LU with U = LT , i.e.,

0 0
   
l11 l11 l21 l31
L = l21 l22 0 , U = LT =  0 l22 l32 
l31 l32 l33 0 0 l33

Multiplying LU elementwise gives:

2
22 24 26
   
l11 l11 l21 l11 l31
LU = l21 l11
 2
l21 + l222
l21 l31 + l22 l32  = 24 29 21
l31 l11 l31 l21 + l32 l22 2
l31 + l32
2
+ l33
2
26 21 51

Step 3: Solve for the elements of L

√
2
l11 = 22 ⇒ l11 = 22

24 26
l21 l11 = 24 ⇒ l21 = √ , l31 l11 = 26 ⇒ l31 = √
22 22

2
24 61
 r
2
l21 + 2
l22 = 29 ⇒ √ + 2
l22 = 29 ⇒ l22 =
22 11

21 − 24·26
l21 l31 + l22 l32 = 21 ⇒ l32 = q 22
61
11

q
2
l31 + l32
2
+ l33
2
= 51 ⇒ l33 = 51 − (l31
2 + l2 )
32

Hence,

√
22 q0 0


L =  √2422 61
0 , U = LT
 
11
√26 l32 l33
22

Step 4: Solve LY = AT B (Forward Substitution)

43
 

LY = 51
40

Prepared By: Er. Narayan Sapkota, [Link]. 80 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

43
l11 y1 = 43 ⇒ y1 = √
22
√
7 11
l21 y1 + l22 y2 = 51 ⇒ y2 =
5
1
l31 y1 + l32 y2 + l33 y3 = 40 ⇒ y3 =
3
 
√43
 √22 
Y =  7 11 
 5 
1
3

Step 5: Solve U X = Y (Backward Substitution)

UX = Y

1 1
l33 z =⇒z=
3√ 3
7 11 6
l22 y + l32 z = ⇒y=
5 5
43 4
l11 x + l21 y + l31 z = √ ⇒ x =
22 15

4 6 1
x= , y= , z=
15 5 3

Conclusion: Although the original matrix is not symmetric, Cholesky’s method can be applied successfully
by first converting the system into its normal equations.

4.4 The Inverse of a Matrix

4.4.1 Gauss-Jordan Method to Find the Inverse of a Matrix
Let  
a11 a12 a13
A = a21 a22 a23  .
a31 a32 a33

Step 1: Write the augmented matrix

1 0 0
 
a11 a12 a13
[A|I] =  a21 a22 a23 0 1 0 
a31 a32 a33 0 0 1

Step 2: Apply row operations to transform A into I

1. Make the pivot in row 1 equal to 1: R1 → R1 /a11 .
2. Eliminate elements below the pivot in column 1: R2 → R2 − a21 R1 , R3 → R3 − a31 R1 .
3. Make the pivot in row 2 equal to 1 and eliminate above/below: R2 → R2 /a22 , R1 → R1 − a12 R2 ,
R3 → R3 − a32 R2 .

Prepared By: Er. Narayan Sapkota, [Link]. 81 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

4. Make the pivot in row 3 equal to 1 and eliminate above: R3 → R3 /a33 , R1 → R1 − a13 R3 , R2 →
R2 − a23 R3 .
Step 3: Result
At the end, the augmented matrix becomes

1 0 0
 
∗ ∗ ∗
[I|A−1 ] =  0 1 0 ∗ ∗ ∗ ,
0 0 1 ∗ ∗ ∗

where the right-hand side gives A−1 .

Example

Find the inverse of the following matrix using Gauss Jordan Method:

2 1 1
 

A = 1 3 2 .
1 0 0

Solution:
Step 1: Write the augmented matrix

2 1 1 1 0 0
 

[A|I] =  1 3 2 0 1 0 
1 0 0 0 0 1

Step 2: Make the pivot in row 1 equal to 1

1 1 1 1
0 0
 
R1 2 2 2
R1 → ⇒ 1 3 2 0 1 0 
2
1 0 0 0 0 1

Step 3: Eliminate the elements below the pivot in column 1

R2 → R2 − R1 , R3 → R3 − R1

1 1 1 1
0 0
 
2 2 2
 0 5 3
− 21 1 0 
2 2
0 − 12 − 21 − 21 0 1

Step 4: Make the pivot in row 2 equal to 1

1 1 1 1
0 0
 
2 2 2 2
R2 → R2 ⇒  0 1 3
− 15 2
0 
5 5 5
0 − 21 − 12 − 12 0 1

Step 5: Eliminate elements above and below the pivot in column 2

1 1
R1 → R1 − R2 , R3 → R3 + R2
2 2

Prepared By: Er. Narayan Sapkota, [Link]. 82 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

1 0 1 11
− 15 0
 
5 20
 0 1 3
− 15 2
0 
5 5
0 0 − 10
7
− 11
20
1
5 1

Step 6: Make the pivot in row 3 equal to 1

1 0 1 11
− 51 0
 
10 5 20
R3 → − R3 ⇒  0 1 3
− 51 2
0 
7 5 5
0 0 1 11
14 − 72 − 10
7

Step 7: Eliminate elements above the pivot in column 3

1 3
R1 → R1 − R3 , R2 → R2 − R3
5 5

1 0 0 0 0 1
 
 0 1 0 −2 1 3 
0 0 1 11
14 − 27 − 10
7

Step 8: Read off the inverse

0 0 1
 

A−1 = −2 1 3
3 −1 −5

4.5 Ill-Conditioned Systems

In numerical analysis, the concept of conditioning of a system of linear equations plays an important role
in determining the accuracy and reliability of computed solutions.
Consider a system of linear equations
AX = B,
where A is a square matrix, X is the vector of unknowns, and B is the constant vector.
Well-Conditioned and Ill-Conditioned Systems
A system is said to be well-conditioned if a small change in the coefficients of A or in the vector B produces
only a small change in the solution X.
On the other hand, a system is called ill-conditioned if a very small change in the coefficients or in the
right-hand side leads to a large change in the solution. In such systems, numerical solutions may be highly
unreliable due to amplification of rounding and truncation errors.
Consider the following two systems of equations:
System I:
x + y = 2,
x + 1.01y = 2.01.

System II:
x + y = 2,
x + 1.001y = 2.001.

Prepared By: Er. Narayan Sapkota, [Link]. 83 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

Although the coefficients differ only slightly, the solutions of these systems differ significantly, illustrating
the sensitivity of ill-conditioned systems.
Condition Number
The degree of conditioning of a matrix A is measured by its condition number, denoted by κ(A).

κ(A) = ∥A∥ · ∥A−1 ∥

where ∥ · ∥ denotes a suitable matrix norm.

• If κ(A) is small, the system is well-conditioned.
• If κ(A) is large, the system is ill-conditioned.
A large condition number indicates that small errors in data or arithmetic operations can result in large
errors in the computed solution.
Effects of Ill-Conditioning
Ill-conditioned systems have the following effects:
• Amplification of round-off errors
• Loss of significant digits in the solution
• Poor convergence of iterative methods
• Numerical instability in direct methods
Detection of Ill-Conditioning
Some common indicators of ill-conditioned systems include:
• Nearly linearly dependent rows or columns of the coefficient matrix
• Very small determinant of the coefficient matrix
• Large condition number
Remedies
To reduce the effects of ill-conditioning, the following measures may be adopted:
• Scaling of equations
• Use of pivoting techniques
• Employing higher precision arithmetic
• Reformulating the problem, if possible
• Using stable numerical algorithms

4.6 Iterative Methods of Solution

The methods discussed earlier for solving systems of simultaneous linear equations are known as direct
methods, since they yield the solution after a fixed number of computations. In contrast, an iterative
method begins with an initial approximation to the solution and produces successively better approxima-
tions through repeated computational cycles until the desired accuracy is achieved.
Thus, in iterative methods, the amount of computation depends on the degree of accuracy required. For
large systems of equations, iterative methods are often faster than direct methods. Moreover, round-off
errors are usually smaller because iteration is a self-correcting process; errors introduced at any stage tend
to be corrected in subsequent steps.

Prepared By: Er. Narayan Sapkota, [Link]. 84 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

Simple iterative methods are particularly effective when the coefficients of the leading diagonal are large
compared to the other coefficients. We now discuss two important iterative methods:
1. Jacobi’s Iteration Method
2. Gauss–Seidel Iteration Method

4.6.1 Jacobi’s Iteration Method

Consider the system of linear equations:
a1 x + b1 y + c1 z = d1 ,
a2 x + b2 y + c2 z = d2 , (4.1)
a3 x + b3 y + c3 z = d3 .

If the coefficients a1 , b2 , c3 are large compared to the other coefficients, we solve each equation for its leading
variable:
1
x= (d1 − b1 y − c1 z) ,
a1
1
y= (d2 − a2 x − c2 z) , (4.2)
b2
1
z= (d3 − a3 x − b3 y) .
c3
Let the initial approximations be
x(0) , y (0) , z (0) .

The first approximations are obtained by substituting these initial values into the right-hand sides:
1  
x(1) = d1 − b1 y (0) − c1 z (0) ,
a1
1  
y (1) = d2 − a2 x(0) − c2 z (0) , (4.3)
b2
1  
z (1) = d3 − a3 x(0) − b3 y (0) .
c3

Similarly, the second approximations are

1  
x(2) = d1 − b1 y (1) − c1 z (1) ,
a1
1  
y (2) = d2 − a2 x(1) − c2 z (1) , (4.4)
b2
1  
z (2) = d3 − a3 x(1) − b3 y (1) .
c3

This process is continued until the difference between two successive approximations is negligible.
Observation: If no better initial estimates are available, one may take
x(0) = y (0) = z (0) = 0.

Solve the system of equations 

10x + y − z = 11.19

x + 10y + z = 28.08
−x + y + 10z = 35.61



using Jacobi iteration method.

Prepared By: Er. Narayan Sapkota, [Link]. 85 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

Step 1: Rewrite each equation in Jacobi form

1
x= (11.19 − y + z)
10
1
y= (28.08 − x − z)
10
1
z= (35.61 + x − y)
10

Step 2: Choose initial approximation

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

First Iteration
1
x(1) = (11.19) = 1.119
10
1
y (1) = (28.08) = 2.808
10
1
z (1) = (35.61) = 3.561
10

Second Iteration
1
x(2) = (11.19 − 2.808 + 3.561) = 1.19
10
1
y (2) = (28.08 − 1.119 − 3.561) = 2.34
10
1
z (2) = (35.61 + 1.119 − 2.808) = 3.39
10

Third Iteration
1
x(3) = (11.19 − 2.34 + 3.39) = 1.22
10
1
y (3) = (28.08 − 1.19 − 3.39) = 2.35
10
1
z (3) = (35.61 + 1.19 − 2.34) = 3.45
10

Fourth Iteration
1
x(4) = (11.19 − 2.35 + 3.45) = 1.23
10
1
y (4) = (28.08 − 1.22 − 3.45) = 2.34
10
1
z (4) = (35.61 + 1.22 − 2.35) = 3.45
10

Fifth Iteration
x(5) = 1.23, y (5) = 2.34, z (5) = 3.45

Final Solution (correct to two decimal places):

x = 1.23, y = 2.34, z = 3.45

Prepared By: Er. Narayan Sapkota, [Link]. 86 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

Example

Solve the system using Jacobi’s method:

20x + y − 2z = 17,
3x + 20y − z = −18,
2x − 3y + 20z = 25.

Solution:
Rewriting the equations:
1
x= (17 − y + 2z),
20
1
y= (−18 − 3x + z),
20
1
z= (25 − 2x + 3y).
20

Starting with x(0) = y (0) = z (0) = 0, successive iterations yield:

(x(1) , y (1) , z (1) ) = (0.85, −0.9, 1.25),

(x(2) , y (2) , z (2) ) = (1.02, −0.965, 1.03),

(x(3) , y (3) , z (3) ) ≈ (1.001, −1.001, 1.003).

Further iterations show no significant change. Hence, the solution is

x = 1, y = −1, z=1.

4.6.2 Gauss–Seidel Iteration Method

The Gauss–Seidel method is a modification of Jacobi’s method. The system

a1 x + b1 y + c1 z = d1 ,
a2 x + b2 y + c2 z = d2 ,
a3 x + b3 y + c3 z = d3

is rewritten as
1
x= (d1 − b1 y − c1 z) ,
a1
1
y= (d2 − a2 x − c2 z) ,
b2
1
z= (d3 − a3 x − b3 y) .
c3

Unlike Jacobi’s method, newly computed values are used immediately.

Starting with initial approximations x(0) , y (0) , z (0) :
• Compute x(1) using y (0) , z (0) ,
• Compute y (1) using x(1) , z (0) ,
• Compute z (1) using x(1) , y (1) .

Prepared By: Er. Narayan Sapkota, [Link]. 87 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

This process is repeated until convergence.

Example

Solve using Gauss–Seidel Iteration Method:

20x + y − 2z = 17,
3x + 20y − z = −18,
2x − 3y + 20z = 25.

Solution:
Rewriting:
1
x= (17 − y + 2z),
20
1
y= (−18 − 3x + z),
20
1
z= (25 − 2x + 3y).
20

Starting with x(0) = y (0) = z (0) = 0:

First iteration:
x(1) = 0.85, y (1) = −1.0275, z (1) = 1.0109.

Second iteration:
x(2) = 1.0025, y (2) = −0.9998, z (2) = 0.9998.

Third iteration:
x(3) = 1.0000, y (3) = −1.0000, z (3) = 1.0000.

Since successive approximations agree closely, we stop.

x = 1, y = −1, z=1.

4.7 Eigenvalues and Eigenvectors

Let A be a square matrix of order n with elements aij . We can find a column vector X and a scalar λ such
that
AX = λX or (A − λI)X = 0,
where I is the identity matrix of order n.
This represents a system of n homogeneous linear equations:

(a11 − λ)x1 + a12 x2 + · · · + a1n xn = 0

a21 x1 + (a22 − λ)x2 + · · · + a2n xn = 0
..
.
an1 x1 + an2 x2 + · · · + (ann − λ)xn = 0

A non-trivial solution exists only if the determinant of the coefficient matrix vanishes:

det(A − λI) = 0

Prepared By: Er. Narayan Sapkota, [Link]. 88 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

This determinant gives an n-th degree equation in λ, called the characteristic equation of A. Its roots
λi , i = 1, 2, . . . , n are called the eigenvalues of A. For each eigenvalue λi , a corresponding non-zero solution
Xi = [x1 , x2 , . . . , xn ]T is called an eigenvector.

4.7.1 Cayley-Hamilton Theorem

Every square matrix satisfies its own characteristic equation. If the characteristic equation of an n×n matrix
A is
det(A − λI) = (−1)n λn + k1 λn−1 + · · · + kn = 0,
then
(−1)n An + k1 An−1 + · · · + kn I = 0.

Example

Find the eigenvalues and eigenvectors of

5 4
 
A= .
1 2

Solution:
The characteristic equation is
5−λ 4
det(A − λI) = =0
1 2−λ

(5 − λ)(2 − λ) − 4 = λ2 − 7λ + 6 = 0
⇒ λ = 6, 1.
For λ = 6:
4
  
−1 x
(A − 6I)X = =0 ⇒ x + 4y = 0 ⇒ x = −4y.
1 −4 y
Hence, an eigenvector is (4, 1).
For λ = 1:
4 4 x
  
(A − I)X = =0 ⇒ x + y = 0.
1 1 y
Hence, an eigenvector is (1, −1).

Example

Find the eigenvalues and eigenvectors of

8 6 3
 

A = 6 7 4 .
2 4 3

Solution:
The characteristic equation is

det(A − λI) = 0 ⇒ (λ − 0)(λ − 3)(λ − 15) = 0

⇒ λ = 0, 3, 15.
Eigenvectors:

Prepared By: Er. Narayan Sapkota, [Link]. 89 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

• For λ = 0:
(A − 0I)X = 0 ⇒ X = (1, 2, 2)

• For λ = 3:
X = (2, 1, −2)

• For λ = 15:
X = (2, −2, 1)

4.7.2 Power Method

In many engineering problems, it is required to compute the numerically largest eigenvalue and the corre-
sponding eigenvector. The Power Method is an iterative procedure well-suited for this purpose and for
machine computations.
Let X1 , X2 , . . . , Xn be the eigenvectors corresponding to eigenvalues λ1 , λ2 , . . . , λn . An arbitrary vector X
can be expressed as
X = k1 X1 + k2 X2 + · · · + kn Xn

Multiplying by A:
AX = k1 λ1 X1 + k2 λ2 X2 + · · · + kn λn Xn

Similarly,
A2 X = k1 λ21 X1 + k2 λ22 X2 + · · · + kn λ2n Xn
Ar X = k1 λr1 X1 + k2 λr2 X2 + · · · + kn λrn Xn

If |λ1 | > |λ2 | > · · · > |λn |, then λ1 is the dominant eigenvalue, and repeated multiplication of X by A makes
it increasingly aligned with X1 , the corresponding eigenvector.
Procedure
1. Start with an initial approximation X (0) to the eigenvector.
2. Compute AX (r) and normalize it to obtain X (r+1) .
3. Approximate the eigenvalue as
X (r+1)T AX (r+1)
λ(r+1) =
X (r+1)T X (r+1)

4. Repeat until [X (r+1) − X (r) ] is negligible.

Note: Using A−1 instead of A in the method allows computation of the smallest eigenvalue.

Example

Find the largest eigenvalue and the corresponding eigenvector of

2 1 0
 

A = 1 2 1
0 1 2

using the initial vector X (0) = [1, 0, 0]T .

Solution:

1
 

X (0) = 0
0

Prepared By: Er. Narayan Sapkota, [Link]. 90 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

2 1 0 1 2
    

AX (0) = 1 2 1 0 = 1

0 1 2 0 0

Normalize by the largest component (2):

1
 

X (1) = 0.5 , λ(1) = 2

0

Next iteration:
2 1 0 1 2.5
    

AX (1) = 1 2 1 0.5 =  2 
0 1 2 0 0.5

Normalize by the largest component (2.5):

1
 

X (2) = 0.8 , λ(2) = 2.5

0.2

Repeating the process:

X (3) = [1, 0.65, 0.54]T , λ(3) = 3.43

X (4) = [1, 0.76, 0.80]T , λ(4) = 3.41
X (5) = [1, 0.67, 0.74]T , λ(5) = 3.41

Since X (5) ≈ X (4) and λ(5) ≈ λ(4) , we stop.

Largest eigenvalue: λmax ≈ 3.41, Corresponding eigenvector: Xmax ≈ [0.74, −1, 0.67]T

Example

Obtain the numerically dominant eigenvalue and eigenvector of the matrix

15 4 3
 

A = 10 12 6
20 4 2

using the initial vector X (0) = [1, 1, 1]T .

Solution:
Iteration 1:
15 4 3 1 22
    

AX (0) = 10 12 6 1 = 28

20 4 2 1 26

Normalize by the largest component (28):

22 0.786
   
1   
X (1) = 28 ≈ 1 , λ(1) ≈ 28
28
26 0.929

Prepared By: Er. Narayan Sapkota, [Link]. 91 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

Iteration 2:
15 4 3 0.786 22.57
    

AX (1) = 10 12 6  1  ≈  27.9 

20 4 2 0.929 25.5

Normalize by the largest component (27.9):

0.809
 

X (2) ≈  1 , λ(2) ≈ 27.9

0.914

Iterations 3–8: Repeating this process:

X (3) ≈ [0.815, 1, 0.907]T , λ(3) ≈ 28.1

X (4) ≈ [0.818, 1, 0.904]T , λ(4) ≈ 28.05
..
.
X (8) ≈ [−1, 0.5, 1]T , λ(8) ≈ 20

Conclusion: Since X (8) ≈ X (7) and λ(8) ≈ λ(7) , the dominant eigenvalue and corresponding eigenvector
are:

λmax ≈ 20, Xmax ≈ [−1, 0.5, 1]T

4.8 Assignment
4.8.1 Gauss Elimination Method
Question 1: Solve the following system of equations by Gauss elimination:

x + 2y + z = 7


2x − y + 3z = 12

3x + y − z = 2



Question 2: Solve:


2x + y + z = 10

x + 3y − z = 5

3x − 2y + 2z = 7



Question 3: Solve:

x + y + 2z = 9


2x − y + z = 8

x + 3y + z = 12



Prepared By: Er. Narayan Sapkota, [Link]. 92 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

4.8.2 Gauss-Jordan Method

Question 1: Solve the following system by Gauss-Jordan elimination:


 x+y+z =6

2x − y + 3z = 14

3x + 4y + z = 20



Question 2: Solve: 
2x + 3y + z = 10


x − y + 2z = 5

3x + y + 4z = 15



Question 3: Solve: 

 x + 2y + 3z = 14

2x + y − z = 3

x−y+z =2



4.8.3 Factorization Method (LU / Doolittle / Crout)

Question 1: Solve the following system using factorization (LU decomposition):


 x+y+z =6

2x + 3y + z = 14

x + 4y + 3z = 20



Question 2: Solve: 

3x + y + 2z = 10

x + 2y + z = 7

2x + 3y + 4z = 20



Question 3: Solve: 

 2x + 3y + z = 9

x + 4y + 2z = 10

3x + y + 3z = 12



4.8.4 Gauss-Seidel Method

Question 1: Solve the following system by Gauss-Seidel method:


 2x + y + 6z = 9

8x + 3y + 2z = 13

x + 5y + z = 7



Question 2: Solve: 
28x + 4y − z = 32


x + 3y + 10z = 24

2x + 17y + 4z = 35



Prepared By: Er. Narayan Sapkota, [Link]. 93 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

Question 3: Solve:

10x + y + z = 12


2x + 10y + z = 13

2x + 2y + 10z = 14



Question 4: Solve:

7x + 52x2 + 13x3 = 104
 1


83x1 + 11x2 − 4x3 = 95

3x1 + 8x2 + 29x3 = 71



Question 5: Solve:

3x − 0.1x2 − 0.2x3 = 7.85
 1


0.1x1 + 7x2 − 0.3x3 = −19.3

0.3x1 − 0.2x2 + 10x3 = 71.4



4.8.5 Jacobi Method

Question 1: Solve the following system by Jacobi iteration:


 2x + y + 6z = 9

8x + 3y + 2z = 13

x + 5y + z = 7



Question 2: Solve:


 28x + 4y − z = 32

x + 3y + 10z = 24

2x + 17y + 4z = 35



Question 3: Solve:


 10x + y + z = 12

2x + 10y + z = 13

2x + 2y + 10z = 14



Question 4: Solve:

7x + 52x2 + 13x3 = 104
 1


83x1 + 11x2 − 4x3 = 95

3x1 + 8x2 + 29x3 = 71



Question 5: Solve:

3x1 − 0.1x2 − 0.2x3 = 7.85


0.1x1 + 7x2 − 0.3x3 = −19.3

0.3x1 − 0.2x2 + 10x3 = 71.4



Prepared By: Er. Narayan Sapkota, [Link]. 94 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

4.8.6 Gauss-Jordan Method – Inverse of a Matrix

Question 1: Use the Gauss-Jordan method to find the inverse of the following matrix:

2 3 4
 

A = 4 3 1
1 2 4

Question 2: Use the Gauss-Jordan method to find the inverse of the following matrix:

0 1 2
 

B = 1 2 3
3 1 1

Question 3: Use the Gauss-Jordan method to find the inverse of the following matrix:

8 4 3
 

C = 2 1 1
1 2 1

4.8.7 Power Method

Question 1: Find the largest eigenvalue and the corresponding eigenvector of

2 1 0
 

A = 1 2 1 .
0 1 2

Question 2: Find the largest eigenvalue and the corresponding eigenvector of

3 1 0
 

B = 1 3 1 .
0 1 3

Question 3: Find the largest eigenvalue and the corresponding eigenvector of

4 2 0
 

C = 2 4 2 .
0 2 4

Question 4: Find the largest eigenvalue and the corresponding eigenvector of

1 2 0
 

D = 2 1 2 .
0 2 1

Prepared By: Er. Narayan Sapkota, [Link]. 95 Everest Engineering College

CHAPTER 4. SOLUTION OF LINEAR ALGEBRAIC EQUATIONS

Prepared By: Er. Narayan Sapkota, [Link]. 96 Everest Engineering College

 Chapter 5
Solution of Ordinary Differential Equations

5.1 Review of Ordinary Differential Equations

5.1.1 Introduction
Many problems in science, engineering, and technology can be formulated as differential equations. These
equations describe the relationship between a function and its derivatives, representing physical phenomena
such as motion, heat transfer, population growth, or electrical circuits. Analytical solutions are only possible
for a limited class of differential equations, and many practical problems lead to equations that do not belong
to these standard types. In such cases, numerical methods provide an essential tool for approximating
solutions. With the widespread availability of computers, these numerical methods have become even more
practical, allowing us to tackle complex problems efficiently.
Numerical methods allow us to obtain approximate solutions when exact solutions are impossible or cum-
bersome. They provide values of the unknown function at discrete points, often with controllable accuracy.
Depending on the method, solutions may be expressed as a series expansion or computed step-by-step
over a range of the independent variable. Understanding the properties and limitations of these methods is
crucial for applying them effectively in scientific and engineering problems.
The solution of an ODE is an expression for y in terms of a finite number of elementary functions of x. This
is called a closed-form solution. When such a solution is not obtainable, we rely on numerical methods
to approximate the solution.

5.1.2 First-Order ODEs

Consider a first-order ODE:

dy
= f (x, y), y(x0 ) = y0
dx
In most numerical methods, the differential equation is replaced by a difference equation, which is then
solved to yield either:
1. Series solutions: Methods such as Picard or Taylor series approximate y by a truncated series in
x. These methods use information at a single point and do not iterate to find the next value. They
are called single-step methods.
2. Step-by-step solutions: Methods such as Euler, Runge-Kutta, Milne, and Adams-Bashforth
compute the solution iteratively, moving from one point to the next in small steps. These methods are
called multi-step or step-by-step methods because each subsequent value of y depends on previous
points.
A boundary value problem (BVP) for a first-order ODE is less common but occurs when the value of
the function is specified at two different points:

y(a) = α, y(b) = β

In such cases, numerical methods like the shooting method are used, where an initial guess for y ′ (a) (if
needed) is iteratively adjusted to satisfy the boundary condition at x = b.
• Euler and Runge-Kutta methods are particularly suited for computing y over a limited range of
x-values.

97
CHAPTER 5. SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS

• Milne and Adams methods are more appropriate for computing over a wider range, but they require
starting values obtained from single-step methods like Runge-Kutta.

5.1.3 Second-Order ODEs

A second-order ODE has the general form:

d2 y
= f (x, y, y ′ )
dx2

Second-order ODEs frequently appear in mechanical vibrations, electrical circuits, and motion under
force. To solve them numerically, they are often reduced to a system of first-order ODEs:

dy dy1 dy2
y1 = y, y2 = ⇒ = y2 , = f (x, y1 , y2 )
dx dx dx

This allows systematic computation of both y and y ′ simultaneously.

Initial Conditions for Second-Order ODEs

For a second-order ODE, the initial value problem (IVP) requires:

y(x0 ) = y0 , y ′ (x0 ) = y0′

This specifies the values of the function and its first derivative at a single point x0 . Numerical methods such
as Euler and Runge-Kutta can then compute the solution over a range of x.

Boundary Conditions for Second-Order ODEs

For a boundary value problem (BVP), the solution must satisfy conditions at two or more points. For
example:

y(a) = α, y(b) = β

Here, the values of the function are prescribed at two different points a and b. BVPs are commonly solved
using shooting methods or finite difference methods, often starting with initial guesses and iterating
until the boundary conditions are satisfied.

5.2 Euler Method

Consider the first-order ordinary differential equation (ODE):

dy
= f (x, y), y(x0 ) = y0 (5.1)
dx
As we have seen, the Taylor-series method requires higher-order derivatives for good accuracy of the solution,
and in many cases, it may be awkward to apply if the derivatives become complicated. Further, we know
that the error in a Taylor series will be small if the step size h is small. In fact, if we make h small enough, we
may only need a few terms of the Taylor-series expansion for good accuracy. The Euler method follows this
idea to the extreme for first-order differential equations—it uses only the first two terms of the Taylor series.
Suppose that we have chosen h small enough so that we may truncate the series after the first derivative
term. Then:

Prepared By: Er. Narayan Sapkota, [Link]. 98 Everest Engineering College

CHAPTER 5. SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS

h2 ′′
y(x0 + h) = y(x0 ) + hy ′ (x0 ) + y (ξ); x0 ≤ ξ ≤ x0 + h
2!

where the last term is the usual error term for the truncated Taylor series. Adopting a subscript notation
for the successive y-values and representing the error by the order relation, we may write the algorithm for
Euler’s method in the form:
yn+1 = yn + hy ′ + O(h2 ) (5.2)
Finally,

yn+1 = yn + hf (xn , yn ), n = 0, 1, 2, 3 (5.3)

This error is just the local error, i.e., the error in a single step. Over many steps, the global error becomes
of O(h).

Example

Solve dy
dx = −(2x + y), y(0) = −1 find y(0.4) taking h = 0.1.

Solution: The formula of the Euler method from (5.3) is

yn+1 = yn + hf (xn , yn ), n = 0, 1, 2, 3, · · ·

and
x0 = 0, x1 = x0 + h = 0.2, x2 = x1 + h = 0.3, x3 = x2 + h = 0.4

At x0 = 0:
y0 = −1, f (x0 , y0 ) = −(2(0.0) + (−1)) = 1, hf (x0 , y0 ) = 0.1 × 1 = 0.1
y1 = y0 + hf (x0 , y0 ) = −1 + 0.1 = −0.9

The remaining calculations are shown in table below:

n xn yn f (xn , yn ) hf (xn , yn ) yn+1 = yn + hf (xn , yn )
0 0.0 -1 1 0.1 -0.9
1 0.1 -0.9 0.7 0.07 -0.83
2 0.2 -0.83 0.43 0.043 -0.787
3 0.3 -0.787 0.187 0.0187 -0.7683
4 0.4 -0.7683 END END END
Since the analytical solution is:

y(x) = −2x + 2 − 3e−x ∴ y(0.4) = −0.81096

The analytical solution is −0.81096. Thus, the error is −0.04266.

Here, we see that there is only one-decimal-place accuracy, even though we have advanced the solution only
for four steps. To gain four-decimal-place accuracy, one has to repeat the process more than 400 times,
taking a very small value of h, say h < 0.00005.

Alternative Method:
We observe that,
yn+1 = yn + hf (xn , yn ), n = 0, 1, 2, 3

Prepared By: Er. Narayan Sapkota, [Link]. 99 Everest Engineering College

CHAPTER 5. SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS

Put f (xn , yn ) = −(2xn + yn ), we get

yn+1 = yn + h(−2 xn − yn ), n = 0, 1, 2, 3

Put values of n we can get:

n = 0, y1 = y0 + h(−2 x0 − y0 ) = −0.9
n = 1, y2 = y1 + h(−2 x1 − y1 ) = −0.83
n = 2, y3 = y2 + h(−2 x2 − y2 ) = −0.787
n = 3, y4 = y3 + h(−2 x3 − y3 ) = −0.7683

Prepared By: Er. Narayan Sapkota, [Link]. 100 Everest Engineering College
CHAPTER 5. SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS

5.3 Runge-Kutta Method for First-Order Differential Equation

Consider the first-order ordinary differential equation (ODE):

dy
= f (x, y), y(x0 ) = y0 (5.4)
dx

Runge-Kutta First Order

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

yn+1 = yn + hf (xn , yn ), n = 0, 1, 2, 3 (5.5)

Runge-Kutta 2nd Order

k1 = f (xn , yn )
k2 = f (xn + h, yn + hk1 )
1
k = (k1 + k2 )
2
yn+1 = yn + h · k (5.6)

Runge-Kutta 4th Order

k1 = f (xn , yn )
 
h k1
k2 = f xn + , yn + h
2 2
 
h k2
k3 = f xn + , yn + h
2 2
k4 = f (xn + h, yn + hk3 )
1
k = (k1 + 2(k2 + k3 ) + k4 )
6
yn+1 = yn + h · k (5.7)

Example

Solve dy
dx = −(2x + y), y(0) = −1 for y(0.4) taking h = 0.1. (b) RK1 (b) RK2 (c) RK4

Solution: We are given the differential equation:

dy
= −(2x + y), y(0) = −1 ∴ f (x, y) = −(2x + y)
dx
and we will compute y from x = 0 to x = 0.4, with step size h = 0.1, so we need 4 steps:

x = 0.1, 0.2, 0.3, 0.4

Since the analytical solution is:

y(x) = −2x + 2 − 3e−x ∴ y(0.4) = −0.81096

(a) Using RK1 Method [Euler Method]

Prepared By: Er. Narayan Sapkota, [Link]. 101 Everest Engineering College
CHAPTER 5. SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS

The RK1 formula is:

yn+1 = yn + hf (xn , yn ), n = 0, 1, 2, 3

Step 1: x0 = 0, y0 = −1

f (x0 , y0 ) = −[2(0) + (−1)] = 1

y1 = −1 + 0.1 · 1 = −0.9

Step 2: x1 = 0.1, y1 = −0.9

f (x1 , y1 ) = −[2(0.1) + (−0.9)] = −[0.2 − 0.9] = 0.7

y2 = −0.9 + 0.1 · 0.7 = −0.83

Step 3: x2 = 0.2, y2 = −0.83

f (x2 , y2 ) = −[2(0.2) + (−0.83)] = −[0.4 − 0.83] = 0.43

y3 = −0.83 + 0.1 · 0.43 = −0.787

Step 4: x3 = 0.3, y3 = −0.787

f (x3 , y3 ) = −[2(0.3) + (−0.787)] = −[0.6 − 0.787] = 0.187

y4 = −0.787 + 0.1 · 0.187 = −0.7683

y(0.4) ≈ −0.7683

(b) Using RK2 Method

The RK2 formulas are:

k1 = f (xn , yn )
k2 = f (xn + h, yn + hk1 )
1
k = (k1 + k2 )
2
yn+1 = yn + h · k

Step 1: x0 = 0, y0 = −1

k1 = f (0, −1) = −[2(0) + (−1)] = 1

k2 = f (0.1, −1 + 0.1 · 1) = f (0.1, −0.9) = −[2(0.1) + (−0.9)] = −[0.2 − 0.9] = 0.7
1
k = (1 + 0.7) = 0.85
2
y1 = −1 + 0.1 · 0.85 = −0.915

Step 2: x1 = 0.1, y1 = −0.915

k1 = f (0.1, −0.915) = −[2(0.1) + (−0.915)] = −[0.2 − 0.915] = 0.715

k2 = f (0.2, −0.915 + 0.1 · 0.715) = f (0.2, −0.8435) = −[0.4 − 0.8435] = 0.4435
1
k = (0.715 + 0.4435) = 0.57925
2
y2 = −0.915 + 0.1 · 0.57925 = −0.857075

Prepared By: Er. Narayan Sapkota, [Link]. 102 Everest Engineering College
CHAPTER 5. SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS

Step 3: x2 = 0.2, y2 = −0.857075

k1 = f (0.2, −0.857075) = −[0.4 − 0.857075] = 0.457075

k2 = f (0.3, −0.857075 + 0.1 · 0.457075) = f (0.3, −0.8113675) = −[0.6 − 0.8113675] = 0.2113675
1
k = (0.457075 + 0.2113675) = 0.33422125
2
y3 = −0.857075 + 0.1 · 0.33422125 = −0.823652875

Step 4: x3 = 0.3, y3 = −0.823652875

k1 = f (0.3, −0.823652875) = 0.223652875

k2 = f (0.4, −0.823652875 + 0.1 · 0.223652875) = f (0.4, −0.8012875875) = 0.0012875875
1
k = (0.223652875 + 0.0012875875) = 0.11247023125
2
y4 = −0.823652875 + 0.1 · 0.11247023125 = −0.812405851875

y(0.4) ≈ −0.8124

(c) Using RK4 Method

The RK4 method formula is given by:

k1 = f (xn , yn )
 
h h
k2 = f xn + , yn + k1
2 2
 
h h
k3 = f xn + , yn + k2
2 2
k4 = f (xn + h, yn + hk3 )
1
k = (k1 + 2k2 + 2k3 + k4 )
6
yn+1 = yn + h · k

We will now apply this method to find y(0.4).

Step 1: x0 = 0, y0 = −1

k1 = f (0, −1) = −(2(0) + (−1)) = 1,

0.1 0.1
 
k2 = f 0 + , −1 + · 1 = f (0.05, −0.95) = −(2(0.05) + (−0.95)) = −(0.1 − 0.95) = 0.85,
2 2
0.1 0.1
 
k3 = f 0 + , −1 + · 0.85 = f (0.05, −0.9575) = 0.8575,
2 2
k4 = f (0.1, −1 + 0.1 · 0.8575) = f (0.1, −0.91425) = 0.71425,
1 1
k = (1 + 2(0.85 + 0.8575) + 0.71425) = (1 + 2(1.7075) + 0.71425) ≈ 0.8549,
6 6
y1 = y0 + 0.1 × 0.8549 = −1 + 0.08549 = −0.91451.

Prepared By: Er. Narayan Sapkota, [Link]. 103 Everest Engineering College
CHAPTER 5. SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS

Step 2: x1 = 0.1, y1 = −0.91451

k1 = f (0.1, −0.91451) = −(2(0.1) + (−0.91451)) = 0.71451,

k2 = f (0.15, −0.91451 + 0.05 · 0.71451) = f (0.15, −0.87829) = 0.57829,
k3 = f (0.15, −0.91451 + 0.05 · 0.57829) = f (0.15, −0.88597) = 0.58597,
k4 = f (0.2, −0.91451 + 0.1 · 0.58597) = f (0.2, −0.85591) = 0.45591,
1
k = (0.71451 + 2(0.57829 + 0.58597) + 0.45591) ≈ 0.58316,
6
y2 = −0.91451 + 0.1 × 0.58316 = −0.91451 + 0.05832 = −0.85619.

Step 3: x2 = 0.2, y2 = −0.85619

k1 = f (0.2, −0.85619) = −(2(0.2) + (−0.85619)) = 0.45619,

k2 = f (0.25, −0.85619 + 0.05 · 0.45619) = f (0.25, −0.83349) = 0.33349,
k3 = f (0.25, −0.85619 + 0.05 · 0.33349) = f (0.25, −0.83945) = 0.33945,
k4 = f (0.3, −0.85619 + 0.1 · 0.33945) = f (0.3, −0.82224) = 0.22224,
1
k = (0.45619 + 2(0.33349 + 0.33945) + 0.22224) ≈ 0.33739,
6
y3 = −0.85619 + 0.1 · 0.33739 = −0.85619 + 0.03374 = −0.82245.

Step 4: x3 = 0.3, y3 = −0.82245

k1 = f (0.3, −0.82245) = −(2(0.3) + (−0.82245)) = 0.22245,

k2 = f (0.35, −0.82245 + 0.05 · 0.22245) = f (0.35, −0.81132) = 0.11132,
k3 = f (0.35, −0.82245 + 0.05 · 0.11132) = f (0.35, −0.81683) = 0.11683,
k4 = f (0.4, −0.82245 + 0.1 · 0.11683) = f (0.4, −0.81077) = 0.01077,
1
k = (0.22245 + 2(0.11132 + 0.11683) + 0.01077) ≈ 0.11492,
6
y4 = −0.82245 + 0.1 · 0.11492 = −0.82245 + 0.01149 = −0.81096.

Thus, we obtain:
y(0.4) ≈ −0.81096

Example

Using the Runge-Kutta method of fourth order, solve

dy y 2 − x2
= 2 , y(0) = 1
dx y + x2
at x = 0.2, 0.4.

Solution:
We are given the differential equation:

dy y 2 − x2 y 2 − x2
= 2 , y(0) = 1 ∴ f (x, y) =
dx y + x2 y 2 + x2

To find y(0.2) and y(0.4), we use the 4th order Runge-Kutta method with step size h = 0.2. The given
differential equation is:

Prepared By: Er. Narayan Sapkota, [Link]. 104 Everest Engineering College
CHAPTER 5. SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS

dy y 2 − x2
= 2 , y(0) = 1
dx y + x2
The Runge-Kutta formulas are:
k1 = f (xn , yn )
 
h h
k2 = f xn + , yn + k1
2 2
 
h h
k3 = f xn + , yn + k2
2 2
k4 = f (xn + h, yn + hk3 )
1
k = (k1 + 2k2 + 2k3 + k4 )
6
yn+1 = yn + h · k
Step 1: Solving for x = 0.2 We start with x0 = 0 and y0 = 1.
y 2 − x2
f (x, y) =
y 2 + x2
At x0 = 0, y0 = 1:
12 − 02 1
k1 = f (0, 1) = 2 = =1
1 + 02 1
0.2 0.2 1.12 − 0.12
 
k2 = f 0 + ,1 + · 1 = f (0.1, 1.1) = ≈ 0.9836
2 2 1.12 + 0.12
0.2 0.2 1.09842 − 0.12
 
k3 = f 0 + ,1 + · 0.9836 = f (0.1, 1.0984) = ≈ 0.9835
2 2 1.09842 + 0.12
1.19672 − 0.22
k4 = f (0.2, 1 + 0.2 · 0.9835) = f (0.2, 1.1967) = ≈ 0.9465
1.19672 + 0.22
1
k = (1 + 2(0.9836 + 0.9835) + 0.9465) ≈ 0.9801
6
y1 = y0 + 0.2 · 0.9801 = 1.1960

Step 2: Solving for x = 0.4 Now, we use x1 = 0.2 and y1 = 1.1960 to calculate the next step:
1.19602 − 0.22
k1 = f (0.2, 1.1960) = ≈ 0.9463
1.19602 + 0.22
0.2 0.2 1.46932 − 0.32
 
k2 = f 0.2 + , 1.1960 + · 0.9463 = f (0.3, 1.4693) = ≈ 0.9211
2 2 1.46932 + 0.32
0.2 0.2 1.46812 − 0.32
 
k3 = f 0.2 + , 1.1960 + · 0.9211 = f (0.3, 1.4681) = ≈ 0.9211
2 2 1.46812 + 0.32
1.48022 − 0.42
k4 = f (0.4, 1.1960 + 0.2 · 0.9211) = f (0.4, 1.4802) = ≈ 0.8632
1.48022 + 0.42
1
k = (0.9463 + 2(0.9211 + 0.9211) + 0.8632) ≈ 0.9157
6
y2 = y1 + 0.2 · 0.9157 = 1.3791

The values of y are:

y(0.2) ≈ 1.1960

y(0.4) ≈ 1.3791

Prepared By: Er. Narayan Sapkota, [Link]. 105 Everest Engineering College
CHAPTER 5. SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS

5.4 Runge-Kutta Method for Second-Order Differential Equation

Consider the second-order ODE:

d2 y dy dy
2
= f (x, y, ), y(x0 ) = y0 , (x0 ) = y0′ (5.1)
dx dx dx

We introduce a substitution to reduce it to first-order equations by letting z = dx .

dy
Thus, we have:

dz d2 y
= 2 = f (x, y, z)
dx dx

Therefore, the second-order equation is reduced to the system of two first-order equations:

dy

= z = ϕ(x, y, z)

dx (5.2)
dz
= f (x, y, z)


dx

This is the system of first-order simultaneous differential equations.

Runge-Kutta First Order

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

yn+1 = yn + hϕ(xn , yn , zn ), n = 0, 1, 2, 3
′
yn+1 = zn+1 = zn + hf (xn , yn , zn ), n = 0, 1, 2, 3

Runge-Kutta 2nd Order

k1 = ϕ(xn , yn , zn ) = zn l1 = f (xn , yn , zn )
k2 = ϕ (xn + h, yn + hk1 , zn + hl1 ) = zn + hl1 l2 = f (xn + h, yn + hk1 , zn + hl1 )
1 1
k = (k1 + k2 ) l = (l1 + l2 )
2 2
yn+1 = yn + hk
′
yn+1 = zn+1 = zn + hl

Runge-Kutta 4th Order

k1 = ϕ(xn , yn , zn ) = zn l1 = f (xn , yn , zn )
 
h k1 l1 
h k1 l1

k2 = ϕ xn + , yn + h , zn + h l2 = f xn + , yn + h , zn + h
2 2 2 2 2 2
 
h k2 l2 
h k2 l2

k3 = ϕ xn + , yn + h , zn + h l3 = f xn + , yn + h , zn + h
2 2 2 2 2 2
k4 = ϕ (xn + h, yn + hk3 , zn + hl3 ) l4 = f (xn + h, yn + hk3 , zn + hl3 )
1 1
k = (k1 + 2k2 + 2k3 + k4 ) l = (l1 + 2l2 + 2l3 + l4 )
6 6
yn+1 = yn + hk ′
yn+1 = zn+1 = zn + hl

Prepared By: Er. Narayan Sapkota, [Link]. 106 Everest Engineering College
CHAPTER 5. SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS

Example

Using the Runge-Kutta method, solve the second-order differential equation

y ′′ = xy ′2 − y 2

for x = 0.2 and h = 0.2

The initial conditions are: x = 0, y = 1, y′ = 0

Solution:

We first convert the second-order ODE into a system of first-order equations by introducing a new variable:

z = y′ , so that z ′ = f (x, y, z) = xz 2 − y 2

The system becomes:

dy
= z = ϕ(x, y, z)
dx
dz
= xz 2 − y 2 = f (x, y, z)
dx

(a) Using RK1 Method [Euler Method]

Using the first-order Runge-Kutta method (Euler’s method), the update formulas are:

yn+1 = yn + h · ϕ(xn , yn , zn ), n = 0, 1, 2, . . .
zn+1 = zn + h · f (xn , yn , zn ), n = 0, 1, 2, . . .

Step 1: At x0 = 0

y0 = 1, z0 = 0, h = 0.2
ϕ(x0 , y0 , z0 ) = z0 = 0
f (x0 , y0 , z0 ) = x0 z02 − y02 = 0 · 02 − 12 = −1
y1 = y0 + h · ϕ(x0 , y0 , z0 ) = 1 + 0.2 · 0 = 1.0
z1 = z0 + h · f (x0 , y0 , z0 ) = 0 + 0.2 · (−1) = −0.2

Step 2: At x1 = 0.2

y1 = 1.0 z1 = −0.2, h = 0.2

ϕ(x1 , y1 , z1 ) = z1 = −0.2
f (x1 , y1 , z1 ) = x1 z12 − y12 = 0.2 · (−0.2)2 − (1)2 = 0.008 − 1 = −0.992

y(0.2) ≈ 1.0

y ′ (0.2) ≈ −0.2

(b) Using RK2 Method

Initial values: x0 = 0, y0 = 1, z0 = 0

Prepared By: Er. Narayan Sapkota, [Link]. 107 Everest Engineering College
CHAPTER 5. SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS

Step 1:

k1 = z0 = 0
l1 = f (0, 1, 0) = 0 · 02 − 12 = −1
k2 = z0 + hl1 = 0 + 0.2 · (−1) = −0.2
l2 = f (0.2, 1 + 0.2 · 0, −0.2) = 0.2 · 0.04 − 1 = 0.008 − 1 = −0.992
1
k = (0 + (−0.2)) = −0.1
2
1
l = (−1 + (−0.992)) = −0.996
2
y1 = y0 + hk = 1 − 0.02 = 0.98
z1 = z0 + hl = 0 + 0.2 · (−0.996) = −0.1992

y(0.2) ≈ 0.98, y ′ (0.2) ≈ −0.1992

(c) Using RK4 Method

Let: x0 = 0, y0 = 1, z0 = 0, h = 0.2

Step 1: Compute k1 , l1

k1 = z0 = 0, l1 = f (x0 , y0 , z0 ) = 0 · 02 − 12 = −1

Step 2: Compute k2 , l2

h
x0 + = 0.1
2
h
y0 + k1 = 1 + 0.1 · 0 = 1
2
h
z0 + l1 = 0 + 0.1 · (−1) = −0.1
2
k2 = ϕ(0.1, 1, −0.1) = −0.1
l2 = f (0.1, 1, −0.1) = 0.1 · 0.01 − 1 = 0.001 − 1 = −0.999

Step 3: Compute k3 , l3

h
y0 + k2 = 1 + 0.1 · (−0.1) = 0.99
2
h
z0 + l2 = 0 + 0.1 · (−0.999) = −0.0999
2
k3 = ϕ(0.1, 0.99, −0.0999) = −0.0999
l3 = f (0.1, 0.99, −0.0999) = 0.1 · (−0.0999)2 − (0.99)2
= 0.00099801 − 0.9801 = −0.97910199

Prepared By: Er. Narayan Sapkota, [Link]. 108 Everest Engineering College
CHAPTER 5. SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS

Step 4: Compute k4 , l4

x0 + h = 0.2
y0 + hk3 = 1 + 0.2 · (−0.0999) = 0.98002
z0 + hl3 = 0 + 0.2 · (−0.97910199) = −0.195820398
k4 = ϕ(0.2, 0.98002, −0.195820398) = −0.195820398
l4 = f (0.2, 0.98002, −0.195820398)
= 0.2 · (−0.1958)2 − (0.98002)2
= 0.2 · 0.03837 − 0.9604392 = 0.007674 − 0.9604392 = −0.9527652

Step 5: Compute RK4 Updates

1
k= (k1 + 2k2 + 2k3 + k4 )
6
1
= (0 + 2(−0.1) + 2(−0.0999) + (−0.19582))
6
−0.59562
= = −0.09927
6
1
l= (l1 + 2l2 + 2l3 + l4 )
6
1
= (−1 + 2(−0.999) + 2(−0.9791) + (−0.9528))
6
−5.909
= = −0.9848
6

y1 = y0 + h · k = 1 + 0.2 · (−0.09927) = 0.98015

z1 = z0 + h · l = 0 + 0.2 · (−0.9848) = −0.19696

Final Answer:
y(0.2) ≈ 0.9802, y ′ (0.2) ≈ −0.1970

Prepared By: Er. Narayan Sapkota, [Link]. 109 Everest Engineering College
CHAPTER 5. SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS

5.5 Boundary Value Problems

Such a problem requires the solution of a differential equation in region R subject to the various conditions
on the boundary of R. Practice applications give rise to many such problems. We shall discuss two-point
linear boundary value problems of the following types:

d2 y dy
(i) 2
+ λ(x) + µ(x)y = γ(x)
dx dx
with the conditions y(x0 ) = a, and y(xn ) = b.

d4 y
(ii) + λ(x)y = µ(x)
dx4

with the conditions y(x0 ) = y ′ (x0 ) = a and

y(xn ) = y ′ (xn ) = b.

There exist two numerical methods for solving such boundary value problems. The first one is known as the
finite difference method which makes use of finite difference equivalents of derivatives. The second one is
called the shooting method which makes use of the techniques for solving initial value problems.

5.6 Shooting Method

Prepared By: Er. Narayan Sapkota, [Link]. 110 Everest Engineering College
CHAPTER 5. SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS

Assignment 5
[A] Solve the following first order differential equation with initial conditions using (a)RK1 method (b) RK2
method (c) RK4 method.

Note: You can solve any 4 from the following questions, but question number 1 is compulsory.
1. Using the Runge-Kutta method, solve for y at x = 1.2 and x = 1.4, given the differential equation:
dy y + ex
=
dx 1 + xex
with initial condition:
x0 = 1, y0 = 0

2. Using the Runge-Kutta method, solve for y at x = 1.2 and x = 1.4, given the differential equation:
dy y + ex
=
dx 1 + xex
with initial condition:
x0 = 1, y0 = 0

3. Use Runge’s method to approximate y when x = 1.1, given that y = 1.2 when x = 1 and dy
dx = 3x + y 2 .
4. Using the Runge-Kutta method, find y(0.2) given that dy
dx = 3x + y 2 , y(0) = 1, taking h = 0.1.
5. Using the Runge-Kutta method, compute y(0.2) and y(0.4) from dy
dx = 2x + y 2 , y(0) = 1, taking
h = 0.1.
6. Use the Runge-Kutta method to find y when x = 1.2 in steps of 0.1, given that dy
dx = x2 + y 2 and
y(1) = 1.5.
7. Given dy
dx = x3 + y, y(0) = 2, compute y(0.2), y(0.4), and y(0.6) by the Runge-Kutta method.
8. Find y(0.1) and y(0.2) using the Runge-Kutta formula, given that y ′ = x2 − y and y(0) = 1.
9. Using the Runge-Kutta method, solve the following equation, taking each step of h = 0.1, given
y(0) = 3. dx
dy
= 4x
y − xy. Calculate y for x = 0.1 and x = 0.2.

[B] Solve the following second order differential equation with initial conditions using (a)RK1 method (b)
RK2 method (c) RK4 method.

Note: You can solve any 4 from the following questions, but question number 1 and 2 is compulsory.
1. Use the Runge-Kutta method to solve the second-order ODE:

d2 y
= −y
dx2
with initial conditions y(0.1) = 0, y ′ (0) = 1, and step size h = 0.1, to compute y(0.2) and y(0.3).
2. Solve the differential equation for y(0.25) using the RRunge-Kutta method:
2
d2 y

dy
10 2 + + 6x = 0
dx dx
with initial conditions:
y(0) = 1, y ′ (0) = 0, h = 0.25
Answer: y(0.25) ≈ 0.9937, y ′ (0.25) ≈ −0.0188

Prepared By: Er. Narayan Sapkota, [Link]. 111 Everest Engineering College
CHAPTER 5. SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS

3. Solve using the RK method:

d2 y dy
+2 +y =0
dx2 dx
Given: y(0) = 1, y ′ (0) = 0, use step size h = 0.2 to compute y(0.2) and y(0.4).
4. Using the Runge-Kutta method, solve:

d2 y dy
2
=x· −y
dx dx
With initial values: y(0) = 2, y ′ (0) = −1, and step size h = 0.1, find y(0.1) and y(0.2).
5. Apply the Runge-Kutta method to solve:

d2 y
= ex − y
dx2
with y(0) = 1, y ′ (0) = 0, and step size h = 0.2. Compute y(0.2) and y(0.4).
6. Use the Runge-Kutta method to solve:

d2 y dy
= sin(x) −
dx2 dx
Given: y(0) = 0, y ′ (0) = 1, with step size h = 0.1, find y(0.1), y(0.2), and y(0.3).

Prepared By: Er. Narayan Sapkota, [Link]. 112 Everest Engineering College
 Chapter 6
Solutions of PDEs

6.1 Classification of PDE

A Partial Differential Equation (PDE) is a mathematical equation that involves an unknown function of
multiple variables and its partial derivatives with respect to those variables. Formally, a PDE can be written
as:

∂ku
 
∂u ∂u
F x1 , x2 , . . . , xn , u(x1 , x2 , . . . , xn ), , ,..., k = 0
∂x1 ∂x2 ∂xn

Here:

• u(x1 , x2 , . . . , xn ) is the unknown function of the independent variables x1 , x2 , . . . , xn , which could

represent physical quantities like temperature, pressure, displacement, etc.

• ∂xi and higher-order

∂u
partial derivatives represent the rates of change of u with respect to the variables
x1 , x2 , . . . , xn .

• F is a general function that relates u, its partial derivatives, and possibly the independent variables.

A general second-order partial differential equation (PDE) in two variables x and y can be written as:

∂2u ∂2u ∂2u

A + B + C + lower-order terms = 0
∂x2 ∂x∂y ∂y 2

Here:

• u(x, y) is the unknown function.

• A, B, C are constants (coefficients) or functions of x only or y only, or both.

This equation is classified based on the discriminant D = b2 − 4ac, where a = A, b = B, and c = C. The
classification is as follows:

• Elliptic: When D < 0, the equation is classified as elliptic. The solution tends to describe steady-
state or equilibrium problems, such as Laplace’s equation or Poisson’s equation. The spatial variable
solutions in this case are typically smooth and continuous.

• Parabolic: When D = 0, the equation is classified as parabolic. The solution typically describes
diffusion or heat flow, such as in the heat equation. Parabolic equations represent processes that
smooth out over time.

• Hyperbolic: When D > 0, the equation is classified as hyperbolic. The solution typically represents
wave propagation or transmission problems, such as the wave equation. Hyperbolic equations describe
systems where information or disturbances propagate with finite speed.

113
CHAPTER 6. SOLUTIONS OF PDES

Example

Classify the following PDEs:

2
∂2u 2
1. ∂∂xu2 − 2 ∂x∂y + ∂∂yu2 = 0
∂2u ∂2u
2. ∂x2 + ∂y 2 + ∂x = 0
∂u

∂2u 2
3. ∂x2 − ∂y 2 + 2 ∂x = 0
∂ u ∂u

∂2u
4. ∂x2 + ∂x + 2 ∂y = 0
∂u ∂u

Solution:
1.
∂2u ∂2u ∂2u
− 2 + =0
∂x2 ∂x∂y ∂y 2
Compare with the general form

Auxx + Buxy + Cuyy + (lower-order terms) = 0

Here:
A = 1, B = −2, C=1
D = B 2 − 4AC = (−2)2 − 4(1)(1) = 4 − 4 = 0
The given PDE is Parabolic

2.
∂ 2 u ∂ 2 u ∂u
+ 2 + =0
∂x2 ∂y ∂x
Here:
A = 1, B = 0, C=1
D = B 2 − 4AC = 02 − 4(1)(1) = −4
The given PDE is Elliptic

3.
∂2u ∂2u ∂u
2
− 2 +2 =0
∂x ∂y ∂x

Here:
A = 1, B = 0, C = −1
D = B 2 − 4AC = 02 − 4(1)(−1) = 4
The given PDE is Hyperbolic

4.
∂ 2 u ∂u ∂u
+ +2 =0
∂x2 ∂x ∂y

Here, there is no second derivative with respect to y:

A = 1, B = 0, C=0

D = B 2 − 4AC = 02 − 4(1)(0) = 0
The given PDE is Parabolic

Prepared By: Er. Narayan Sapkota, [Link]. 114 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

6.2 Finite Difference Method

If y(x) and its derivatives are single-valued continuous functions of x, then by Taylor’s expansion, we have

h2 ′′ h3
y(x + h) = y(x) + hy ′ (x) + y (x) + y ′′′ (x) + . . . (1)
2! 3!
and

h2 ′′ h3
y(x − h) = y(x) − hy ′ (x) + y (x) − y ′′′ (x) + . . . (2)
2! 3!
Equation (1) gives:

y(x + h) − y(x)
y ′ (x) = + O(h) (3)
h
which is the forward difference approximation of y ′ (x) with an error of the order h.
Similarly, equation (2) gives:

y(x) − y(x − h)
y ′ (x) = + O(h) (4)
h
which is the backward difference approximation of y ′ (x) with an error of the order h.
Subtracting equation (2) from equation (1), we obtain:

y(x + h) − y(x − h)
y ′ (x) = + O(h2 ) (5)
2h
which is the central difference approximation of y ′ (x) with an error of the order h2 . Clearly, this central
difference approximation is better than the forward or backward difference approximations and should be
preferred.
Adding equations (1) and (2), we get:

y(x + h) − 2y(x) + y(x − h)

y ′′ (x) = + O(h2 ) (6)
h2
which is the central difference approximation of y ′′ (x). Similarly, we can derive central difference approxi-
mations to higher derivatives.
Hence, the working expressions for the central difference approximations are as follows:
For the first derivative:

yi+1 − yi−1
yi′ =
2h
For the second derivative:

yi+1 − 2yi + yi−1

yi′′ =
h2
Observation: The accuracy of this method depends on the size of the sub-interval h and the order of
approximation. As we reduce h, the accuracy improves, but the number of equations to be solved also
increases.

Prepared By: Er. Narayan Sapkota, [Link]. 115 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

6.3 Elliptic Equations

The Laplace equation and the Poisson’s equation are examples of elliptic partial differential equations.
The Laplace equation arises in steady-state flow and potential problems, and is given by:

∇2 u = uxx + uyy = 0, (6.1)

where u(x, y) is the unknown function of the spatial coordinates x and y.

The Poisson equation arises in fluid mechanics, electricity and magnetism, and torsion problems, and is
given by:
∇2 u = uxx + uyy = f (x, y), (6.2)
where f (x, y) is a known source function.
The solution of these equations is a function u(x, y) which is satisfied at every point of a region R subject
to certain boundary conditions specified on the closed curve C (Figure 6.1).
In general, problems concerning steady viscous flow, equilibrium stresses in elastic structures, etc., lead to
elliptic type of equations.

Figure 6.1

6.3.1 Laplace Equation

The Laplace equation is given by:

∂2u ∂2u
∇2 u(x, y) = 0 or + 2 =0 (6.3)
∂x2 ∂y

Discretization using the 5-Point Formula

In this problem, both are space variables. Now, we approximate this PDE by the central difference approx-
imation. Then the finite difference approximation of the above equation is,

∂2u u(i + 1, j) − 2u(i, j) + u(i − 1, j)

≈
∂x2 h2

∂2u u(i, j + 1) − 2u(i, j) + u(i, j − 1)

≈
∂y 2 k2

Prepared By: Er. Narayan Sapkota, [Link]. 116 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

Figure 6.3: New 5-point configuration

Figure 6.2: Old 5-point configuration

Substituting these into the Laplace equation gives:

u(i + 1, j) − 2u(i, j) + u(i − 1, j) u(i, j + 1) − 2u(i, j) + u(i, j − 1)

+ =0
h2 k2

It is assumed that the length of the subintervals along x and y directions are equal, i.e., h = k. Then, the
above equations becomes

1
ui,j = (ui+1,j + ui−1,j + ui,j+1 + ui,j−1 ) (6.4)
4

This equation represents the discretized form of the Laplace equation using the 5-point standard formula,
which can be used to iteratively solve for u(i, j) at each grid point.

Jacobi Iteration Method

(r+1) 1  (r) (r) (r) (r)


ui,j = ui+1,j + ui−1,j + ui,j+1 + ui,j−1 (6.5)
4

Gauss Seidal Method

(r+1) 1  (r) (r+1) (r) (r+1)


ui,j = ui+1,j + ui−1,j + ui,j+1 + ui,j−1 (6.6)
4

Example 1

Consider a steel plate of size 15 cm × 15 cm, where two sides are held at 100◦ C and the other two
sides are held at 0◦ C. Find the steady-state temperature at the interior points assuming a grid size
of 5 cm × 5 cm.

Solution:

Prepared By: Er. Narayan Sapkota, [Link]. 117 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

100◦ C
15

u1 u2
10

0◦ C 0◦ C
u3 u4
5

0
0 5 100◦ C 10 15

The steady–state temperature distribution satisfies Laplace’s equation. Using the finite difference approxi-
mation,
1
ui,j = (ui+1,j + ui−1,j + ui,j+1 + ui,j−1 ) .
4

With grid spacing h = 5 cm, the plate is divided into 4 × 4 nodes. There are four interior points, which we
denote as

u1 = u(5, 10),
u2 = u(10, 10),
u3 = u(5, 5),
u4 = u(10, 5).

The boundary temperatures are:

u(x, 15) = 100◦ C,
u(x, 0) = 100◦ C,
u(0, y) = 0◦ C,
u(15, y) = 0◦ C.

Applying the finite difference formula at each interior point gives:

1
u1 = (100 + u3 + 0 + u2 ),
4

1
u2 = (100 + u4 + u1 + 0),
4

1
u3 = (100 + 0 + u4 + u1 ),
4

1
u4 = (100 + 0 + u3 + u2 ).
4

By symmetry of the problem,

u1 = u2 = u3 = u4 = u.

Prepared By: Er. Narayan Sapkota, [Link]. 118 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

Substituting into any equation,

1
u= (100 + u + u),
4

4u = 100 + 2u,

2u = 100, ⇒ u = 50◦ C.

Therefore, the steady–state temperature at all interior points is

xy 0 5 10 15
15 100 100 100 100
10 0 50 50 0
5 0 50 50 0
0 100 100 100 100

Example 2

Consider the following Laplace equation:

uxx + uyy = 0

with the boundary conditions:

u(x, 0) = 0, u(0, y) = 0, u(x, 1) = 5x, u(1, y) = 5y.

Find the approximate values at the interior meshes by dividing the square region [0, 1] × [0, 1] into a
4 × 4 grid.

Solution:
Step 1: Grid and mesh size
Dividing [0, 1] × [0, 1] into a 4 × 4 grid gives
1
h= .
4
The grid points are
xi = ih, yj = jh, i, j = 0, 1, 2, 3, 4.
There are 9 interior points:

(1, 1), (2, 1), (3, 1), (1, 2), (2, 2), (3, 2), (1, 3), (2, 3), (3, 3).

Step 2: Boundary values

u(x, 0) = 0, u(0, y) = 0

u(0.25, 1) = 1.25,

u(x, 1) = 5x ⇒ u(0.5, 1) = 2.5,
u(0.75, 1) = 3.75,



u(1, 0.25) = 1.25,

u(1, y) = 5y ⇒ u(1, 0.5) = 2.5,
u(1, 0.75) = 3.75.



Prepared By: Er. Narayan Sapkota, [Link]. 119 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

Step 3: Finite difference equation

For Laplace’s equation,
1
ui,j = (ui+1,j + ui−1,j + ui,j+1 + ui,j−1 ) .
4

Step 4: Label interior points

u13 u23 u33
u12 u22 u32
u11 u21 u31

Step 5: Difference equations

1
u11 = (u21 + u12 )
4
1
u21 = (u11 + u31 + u22 )
4
1
u31 = (u21 + u32 + 1.25)
4
1
u12 = (u11 + u22 + u13 )
4
1
u22 = (u12 + u21 + u23 + u32 )
4
1
u32 = (u22 + u31 + u33 + 2.5)
4
1
u13 = (u12 + u23 + 1.25)
4
1
u23 = (u13 + u22 + u33 + 2.5)
4
1
u33 = (u23 + u32 + 3.75 + 3.75)
4

Prepared By: Er. Narayan Sapkota, [Link]. 120 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

Step 6: Gauss–Jacobi iteration The Jacobi iteration formula is

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

ui,j = ui+1,j + ui−1,j + ui,j+1 + ui,j−1 .
4

Assume the initial guess:

(0)
ui,j = 0 for all interior points.

Iterating until convergence gives the approximate solution.

Step 7: Converged numerical values

u11 ≈ 0.469, u21 ≈ 0.938, u31 ≈ 1.641,

u12 ≈ 0.938, u22 ≈ 1.875, u32 ≈ 2.813,
u13 ≈ 1.641, u23 ≈ 2.813, u33 ≈ 3.516.

Example 3

Consider a steel plate of size 20 cm × 15 cm, where upper and lower two sides are held at 80◦ C and the
other two sides are held at 20◦ C. Find the steady-state temperature at the interior points assuming
a grid size of 4 cm × 5 cm.

Solution:
Using central finite differences with unequal mesh widths h = 4 cm (along x-direction) and k = 5 cm (along
y-direction), we have
ui+1,j − 2ui,j + ui−1,j ui,j+1 − 2ui,j + ui,j−1
+ = 0.
h2 k2
Solving for ui,j ,
k 2 ui+1,j + ui−1,j + h2 ui,j+1 + ui,j−1



ui,j = .
2(k 2 + h2 )

Substituting h = 4 and k = 5,
25 ui+1,j + ui−1,j + 16 ui,j+1 + ui,j−1



ui,j = .
82

Grid arrangement:
The plate dimensions give
20 15
=5 intervals in x, =3 intervals in y.
4 5

Thus, there are 4 × 2 = 8 interior grid points, denoted as

u12 u22 u32 u42

u11 u21 u31 u41

Boundary conditions:
u(x, 0) = 80◦ C,
u(x, 15) = 80◦ C,
u(0, y) = 20◦ C,
u(20, y) = 20◦ C.

Prepared By: Er. Narayan Sapkota, [Link]. 121 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

Finite difference equations:

25(u21 + 20) + 16(u12 + 80)

u11 = ,
82
25(u31 + u11 ) + 16(u22 + 80)
u21 = ,
82
25(u41 + u21 ) + 16(u32 + 80)
u31 = ,
82
25(20 + u31 ) + 16(u42 + 80)
u41 = ,
82

25(u22 + 20) + 16(80 + u11 )

u12 = ,
82
25(u32 + u12 ) + 16(80 + u21 )
u22 = ,
82
25(u42 + u22 ) + 16(80 + u31 )
u32 = ,
82
25(20 + u32 ) + 16(80 + u41 )
u42 = .
82

Solution of the system:

Solving the above system of linear equations (using Gauss–Seidel iteration), the steady–state temperatures
at the interior points are approximately

u11 ≈ 58.6, u21 ≈ 64.8, u31 ≈ 64.8, u41 ≈ 58.6,

u12 ≈ 58.6, u22 ≈ 64.8, u32 ≈ 64.8, u42 ≈ 58.6.

6.3.2 Poisson Equation

The Poisson equation is given by:
∇2 u(r) = g(r)
where:
• ∇2 is the Laplacian operator,
• u(r) is the unknown scalar function (the potential, for example),
• g(r) is a known source term.
In two-dimensional Cartesian coordinates, the equation becomes:

∂ 2 u(x, y) ∂ 2 u(x, y)
+ = g(x, y)
∂x2 ∂y 2

Discretization
To apply the Finite Difference Method, we discretize the spatial domain (grid) and approximate the deriva-
tives by finite differences.
Consider a grid with Nx points in the x-direction and Ny points in the y-direction. Let hx and hy represent
the grid spacing in the x- and y-directions, respectively. Using a central difference scheme for the second
derivative, we approximate the Laplacian as:

Prepared By: Er. Narayan Sapkota, [Link]. 122 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

∂ 2 u(x, y) u(x + hx , y) − 2u(x, y) + u(x − hx , y)

2
≈
∂x h2x

∂ 2 u(x, y) u(x, y + hy ) − 2u(x, y) + u(x, y − hy )

≈
∂y 2 h2y

Thus, the Poisson equation becomes:

u(x + hx , y) − 2u(x, y) + u(x − hx , y) u(x, y + hy ) − 2u(x, y) + u(x, y − hy )

+ = g(x, y)
h2x h2y

Discretized Equation at each grid point (xi , yj ), the equation becomes:

ui+1,j − 2ui,j + ui−1,j ui,j+1 − 2ui,j + ui,j−1

2
+ = gi,j
hx h2y

where ui,j is the value of u(xi , yj ) at the grid point (xi , yj ), and gi,j is the value of the source term g(xi , yj )
at that point.
It is assumed that the length of the subintervals along x and y directions are equal, i.e., hx = hy = h. Then,
the above equations becomes,

1
ui,j = (ui+1,j + ui−1,j + ui,j+1 + ui,j−1 ) − h2 gi,j

4

Boundary Conditions
To fully solve the Poisson equation, boundary conditions are required. These are often given as Dirichlet
(specified u values) or Neumann (specified ∂u/∂n) boundary conditions.
For example:
• Dirichlet boundary conditions: u(x, y) is specified on the boundary.
• Neumann boundary conditions: The gradient of u(x, y) is specified on the boundary.

Example 1

Solve the equation ∇2 u = −10(x2 + y 2 + 10) over the square with sides x = 0, y = 0, x = 3, y = 3,
with u = 0 on the boundary and mesh length h = 1.

Solution:
The Poisson equation in two dimensions is

uxx + uyy = −10(x2 + y 2 + 10).

Using central finite difference approximations,

ui+1,j − 2ui,j + ui−1,j ui,j+1 − 2ui,j + ui,j−1
uxx ≈ , uyy ≈ .
h2 h2

Thus,
ui+1,j + ui−1,j + ui,j+1 + ui,j−1 − 4ui,j
= −10(x2i + yj2 + 10).
h2

Prepared By: Er. Narayan Sapkota, [Link]. 123 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

Since h = 1,
ui+1,j + ui−1,j + ui,j+1 + ui,j−1 − 4ui,j = −10(x2i + yj2 + 10).

The grid points are x, y = 0, 1, 2, 3. and the interior points are (1, 1), (2, 1), (1, 2), (2, 2).
Let u11 = u(1, 1), u21 = u(2, 1), u12 = u(1, 2), u22 = u(2, 2).
Finite difference equations:
At (1, 1):
u21 + u12 − 4u11 = −10(12 + 12 + 10) = −120,
4u11 − u21 − u12 = 120.

At (2, 1):
u11 + u22 − 4u21 = −10(22 + 12 + 10) = −150,
4u21 − u11 − u22 = 150.

At (1, 2):
u11 + u22 − 4u12 = −10(12 + 22 + 10) = −150,
4u12 − u11 − u22 = 150.

At (2, 2):
u12 + u21 − 4u22 = −10(22 + 22 + 10) = −180,
4u22 − u12 − u21 = 180.

System of equations:
4u11 − u21 − u12 = 120,
−u11 + 4u21 − u22 = 150,
−u11 + 4u12 − u22 = 150,
−u12 − u21 + 4u22 = 180.

Solving the above system of linear equations, we obtain

u11 = 70, u21 = 95, u12 = 95, u22 = 110.

Tabular form of the solution:

y\x 0 1 2 3
3 0 0 0 0
2 0 95 110 0
1 0 70 95 0
0 0 0 0 0

Hence, the approximate solution at the interior mesh points is

u(1, 1) = 70, u(2, 1) = 95, u(1, 2) = 95, u(2, 2) = 110.

Example 2

Solve the Poisson equation ∇2 u = 2x2 y 2 over the square domain 0 < x < 3 and 0 < y < 3, with
Dirichlet boundary conditions u(x, y) = 0 on the boundary. The grid spacing is h = k = 1, and the
Gauss-Seidel method is to be used for the solution.

Solution:

Prepared By: Er. Narayan Sapkota, [Link]. 124 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

Step 1: Finite difference approximation

The Poisson equation in 2D is
uxx + uyy = 2x2 y 2 .

Using central differences with uniform mesh h = k = 1:

ui+1,j − 2ui,j + ui−1,j ui,j+1 − 2ui,j + ui,j−1
+ = 2x2i yj2 .
h2 k2

Since h = k = 1:
ui+1,j + ui−1,j + ui,j+1 + ui,j−1 − 4ui,j = 2x2i yj2 .

Solving for ui,j :

1 
ui,j = ui+1,j + ui−1,j + ui,j+1 + ui,j−1 − 2x2i yj2 .
4
Step 2: Grid and interior points
The domain is x, y = 0, 1, 2, 3.
Interior points: (1, 1), (2, 1), (1, 2), (2, 2).
Let
u11 = u(1, 1), u21 = u(2, 1), u12 = u(1, 2), u22 = u(2, 2).

Step 3: Finite difference equations at interior points

1  1
u11 = u21 + u12 + u01 + u10 − 2(12 · 12 ) = (u21 + u12 + 0 + 0 − 2)
4 4
1 1
u21 = (u31 + u11 + u22 + u20 − 2(22 · 12 )) = (0 + u11 + u22 + 0 − 8)
4 4
1 1
u12 = (u22 + u02 + u11 + u13 − 2(12 · 22 )) = (u22 + 0 + u11 + 0 − 8)
4 4
1 1
u22 = (u32 + u12 + u23 + u21 − 2(22 · 22 )) = (0 + u12 + 0 + u21 − 32)
4 4

Boundary points are zero (u = 0).

Step 4: Gauss-Seidel iteration formula
(0)
Starting with an initial guess uij = 0, the Gauss-Seidel updates are:

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


u11 = u21 + u12 − 2 ,
4
(k+1) 1  (k+1) (k)

u21 = u11 + u22 − 8 ,
4
(k+1) 1  (k+1) (k)

u12 = u11 + u22 − 8 ,
4
(k+1) 1  (k+1) (k+1)

u22 = u12 + u21 − 32 .
4
Step 5: Iterative solution
After several Gauss-Seidel iterations, the solution converges to approximately:

u11 ≈ −1.93, u21 ≈ −2.49, u12 ≈ −2.49, u22 ≈ −8.48.

Prepared By: Er. Narayan Sapkota, [Link]. 125 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

Step 6: Tabular form

y\x 0 1 2 3
3 0 0 0 0
2 0 −2.49 −8.48 0
1 0 −1.93 −2.49 0
0 0 0 0 0

Example 3

Solve the Poisson equation

uxx + uyy = −81xy, 0 < x < 1, 0 < y < 1,

with boundary conditions

u(0, y) = 0, u(x, 0) = 0, u(1, y) = 100, u(x, 1) = 100,

using a uniform mesh size h = 13 .

Solution:
Here, h = 13 . The standard five-point finite difference formula for the Poisson equation is:

ui−1,j + ui+1,j + ui,j−1 + ui,j+1 − 4ui,j

= f (xi , yj ),
h2
or equivalently,
−4ui,j + ui−1,j + ui+1,j + ui,j−1 + ui,j+1 = h2 f (xi , yj ).

1 2
Here, f (x, y) = −81xy, and xi = ih, yj = jh. Since h2 = = 19 , we have


3

1
h2 f (xi , yj ) = (−81xi yj ) = −9xi yj .
9

Let the interior points be labeled u1 , u2 , u3 , u4 as:

u3 u4
u1 u2

The finite difference equations at each interior point are:

For u1 (i = 1, j = 1) : 0 + u2 + 0 + u3 − 4u1 = −9(1 · 1) = −9

=⇒ −4u1 + u2 + u3 = −9,

For u2 (i = 2, j = 1) : u1 + 100 + 0 + u4 − 4u2 = −9(2 · 1) = −18

=⇒ u1 − 4u2 + u4 = −118,

For u3 (i = 1, j = 2) : 0 + u4 + u1 + 100 − 4u3 = −9(1 · 2) = −18

=⇒ u1 − 4u3 + u4 = −218,

For u4 (i = 2, j = 2) : u3 + 100 + u2 + 100 − 4u4 = −9(2 · 2) = −36

=⇒ u2 + u3 − 4u4 = −236.

Prepared By: Er. Narayan Sapkota, [Link]. 126 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

From these equations, we solve step by step:

Subtract the first equation from the third: (u1 − 4u3 + u4 ) − (−4u1 + u2 + u3 ) = −218 − (−9)
=⇒ 5u1 − u2 − 5u3 + u4 = −209.

(Step-by-step Gaussian elimination proceeds similarly, ultimately giving:)

u1 = u4 = 51.0833, u2 = 76.5477, u3 = 25.7916.

Thus, the approximate solution at interior points is:

u3 = 25.792 u4 = 51.083
u1 = 51.083 u2 = 76.548

6.4 Parabolic Equations

6.4.1 1D Heat Equation
The 1D heat equation is a partial differential equation (PDE) given by:

∂u(x, t) ∂ 2 u(x, t) k
= α2 and α2 = (6.1)
∂t ∂x2 ρcp

where:
• u(x, t) is the temperature at position x and time t,
• k is the thermal conductivity,
• ρ is the density of the material,
• cp is the specific heat capacity at constant pressure,
• α is the thermal diffusivity constant,
∂ 2 u(x,t)
• ∂x2 is the spatial second derivative of the temperature.

The solution is not defined in a closed form (as in elliptic equations) but propagates in an open-ended region
from initial values satisfying the prescribed boundary conditions (see Figure).

The Schmidt method is an explicit finite difference method used to solve the heat equation numerically. This
involves discretizing the spatial and time domains.

Discretization of the Heat Equation

Let the spatial domain x ∈ [0, L] be divided into N equal intervals with grid spacing ∆x, and the time
domain t ∈ [0, T ] be divided into intervals with time step ∆t.
• The spatial grid points are xi = i∆x for i = 0, 1, 2, . . . , N ,
• The time levels are tn = n∆t for n = 0, 1, 2, . . . , M .

The finite difference approximations for the derivatives are:

Prepared By: Er. Narayan Sapkota, [Link]. 127 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

• Time Derivative:
∂u(xi , tn ) un+1 − uni
≈ i
∂t ∆t
• Second Spatial Derivative:
∂ 2 u(xi , t) un − 2uni + uni−1
≈ i+1
∂x 2 (∆x)2

Finite Difference Scheme for the Schmidt Method

Substituting these approximations into the heat equation:

un+1 − uni un − 2uni + uni−1

i
= α2 i+1
∆t (∆x)2

Rearranging the equation:

α2 ∆t n
un+1 = uni + ui+1 − 2uni + uni−1


i
(∆x)2

α2 ∆t
Put λ = (∆x)2 , above equation becomes

un+1
i = λuni−1 + (1 − 2λ)uni + λuni+1 (6.2)

This is the Schmidt Explict formula for solving the 1D heat equation. This formula is valid only for 0 < λ ≤ 21 .
The temperature at the next time step un+1i is computed based on the values at the current time step.

Bendre-Schmidt Recurrence Relation

when λ = 1/2 Eq. 6.2 reduces to

1 n
un+1 = ui−1 + uni+1 (6.3)

i
2
This shows that the value of u at xi at time tn+1 is the mean of the u-values at xi−1 and xi+1 at time tn .
This relation, known as the Bendre-Schmidt recurrence relation, gives the values of u at the internal mesh

Prepared By: Er. Narayan Sapkota, [Link]. 128 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

points with the help of boundary conditions.

Initial and Boundary Conditions

1. Initial Condition: The initial temperature distribution is specified as:
u(x, 0) = f (x), for x ∈ [0, L]
which gives u0i for each spatial grid point.
2. Boundary Conditions: The boundary conditions specify the temperature at the boundaries:
• Dirichlet Boundary Condition: u(0, t) = u(L, t) = 0, which means the temperature at both ends of the
bar is fixed.
• Neumann Boundary Condition: ∂u
∂x = 0, meaning no heat flux at the boundaries.

Example 1

Solve the boundary value problem ut = uxx under the conditions u(0, t) = 0, u(1, t) = 0, and
u(x, 0) = sin(πx), where 0 ≤ x ≤ 1, using the Bendre-Schmidt method. (Take h = 0.2.

Solution:
Step 1: Identify parameters Compare
ut = uxx with ut = α2 uxx .
We have α = 1 and h = ∆x = 0.2.
The stability parameter is
α2 ∆t
λ= .
(∆x)2
Taking λ = 0.1 (typical small value for stability in explicit methods),
∆t = λ(∆x)2 = 0.02.
Step 2: Discretization using Bender-Schmidt
The Bender-Schmidt (explicit) scheme is
1 n
un+1 = ui−1 + uni+1 , i = 1, 2, . . . , 4


i
2
with boundary points
un0 = 0, un5 = 0.

Step 3: Grid points and initial condition

Divide 0 ≤ x ≤ 1 into 5 intervals (h = 0.2):

x0 = 0, x1 = 0.2, x2 = 0.4, x3 = 0.6, x4 = 0.8, x5 = 1

Initial condition u(x, 0) = sin(πx):

u00 = 0,
u01 = sin(π · 0.2) = 0.5875,
u02 = sin(π · 0.4) = 0.9511,
u03 = sin(π · 0.6) = 0.9511,
u04 = sin(π · 0.8) = 0.5875,
u05 = sin(π · 1) = 0.

Prepared By: Er. Narayan Sapkota, [Link]. 129 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

Step 4: First time step using Bender-Schmidt

For n = 0 → n = 1, compute interior points:

1 0 1
u11 = (u + u02 ) = (0 + 0.9511) = 0.4756,
2 0 2
1 1
u12 = (u01 + u03 ) = (0.5875 + 0.9511) = 0.7693,
2 2
1 1
u13 = (u02 + u04 ) = (0.9511 + 0.5875) = 0.7693,
2 2
1 1
u14 = (u03 + u05 ) = (0.9511 + 0) = 0.4756.
2 2

Boundary points remain:

u10 = 0, u15 = 0.

Step 5: Tabular form of solution (first time step)

x 0 0.2 0.4 0.6 0.8 1

u(x, 0) 0 0.5875 0.9511 0.9511 0.5875 0
u(x, ∆t) 0 0.4756 0.7693 0.7693 0.4756 0

Step 6: Subsequent time steps

Repeat the same Bender-Schmidt iteration for n = 1, 2, . . . until the desired time is reached.

1 n
un+1 = ui−1 + uni+1 , i = 1, . . . , 4.


i
2

This completes the solution using the Bender-Schmidt explicit method.

Prepared By: Er. Narayan Sapkota, [Link]. 130 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

6.5 Exercise
1. Solve the partial differential equation uxx + uyy = 0 over a square domain of side 4 units, subject to
the following boundary conditions:

u(0, y) = 0, 0 ≤ y ≤ 4,
u(4, y) = 12 + y, 0 ≤ y ≤ 4,
u(x, 0) = 3x, 0 ≤ x ≤ 4,
u(x, 4) = x ,2
0 ≤ x ≤ 4.

2. Solve the Poisson equation: ∇2 f = 2x2 + y over the square domain 1 ≤ x ≤ 4 and 1 ≤ y ≤ 4, with the
boundary condition f = 0 on the boundary. The grid spacing is h = k = 1.
3. Consider a steel plate of size 20 cm × 15 cm, where upper and lower two sides are held at 80◦ C and the
other two sides are held at 20◦ C. Find the steady-state temperature at the interior points assuming a
grid size of 4 cm × 5 cm. The Poisson equation is ∇2 u = −10(x2 + y 2 + 10).
4. Solve the heat equation
ut = uxx , 0 ≤ x ≤ 1, t > 0,
subject to the boundary conditions

u(0, t) = 0, u(1, t) = 0,

and the initial condition

u(x, 0) = sin(πx), 0 ≤ x ≤ 1,
using the Schmidt method with
1
h = 0.2, α= .
2
5. Solve the Laplace equation
uxx + uyy = 0
over a rectangular plate of size 2 m × 1 m, with boundary conditions:

u(0, y) = 0, u(2, y) = 100, u(x, 0) = 0, u(x, 1) = 50.

Assume a uniform grid spacing of h = 0.5 in both x and y directions.

6. Solve the Poisson equation
∇2 u = −50(x + y)
on the square domain 0 ≤ x ≤ 2, 0 ≤ y ≤ 2, with boundary condition u = 0 on all sides. Use a grid
spacing of h = k = 1.
7. A metal plate of size 10 cm × 10 cm has its edges maintained at the following temperatures:

u(x, 0) = 0, u(x, 10) = 100, u(0, y) = 50, u(10, y) = 50.

Solve for the steady-state temperature distribution using the finite difference method with a grid size
of 2 cm × 2 cm.
8. Solve the one-dimensional heat equation

ut = αuxx , 0 ≤ x ≤ L, t > 0,

subject to
u(0, t) = u(L, t) = 0, u(x, 0) = f (x),
. Let ∆x = 0.1, ∆t = 0.01, and α = 0.5.

Prepared By: Er. Narayan Sapkota, [Link]. 131 Everest Engineering College
CHAPTER 6. SOLUTIONS OF PDES

9. Consider a square plate of size 1 m × 1 m with the following boundary conditions:

u(x, 0) = 0, u(x, 1) = 100, u(0, y) = 0, u(1, y) = 100.

Solve the Laplace equation

uxx + uyy = 0
using the finite difference method with a grid spacing of h = k = 0.25.
10. Consider a rectangular metal plate of size 3 m × 2 m. The left and right edges are held at 0◦ C, and the
top and bottom edges are held at 100◦ C. Solve for the steady-state temperature distribution u(x, y)
using the finite difference method. Assume a uniform grid spacing of 1 m in both directions.

Prepared By: Er. Narayan Sapkota, [Link]. 132 Everest Engineering College