GE: Numerical Methods

(Class Notes)

Dr. Avnish Kumar Sharma

Assistant Professor
Department of Mathematics
SRCC, DU

August, 2024
Contents

1 Preliminaries 1
1.1 Significant Digits . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 Significant Digits in Rounding Off . . . . . . . . . . . . 2
1.2 Errors in Numerical Methods . . . . . . . . . . . . . . . . . . 3
1.2.1 Types of Errors . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Exercise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2 Roots of Transcendental and Polynomial Equations 7

2.1 Fixed-Point Iteration Method . . . . . . . . . . . . . . . . . . 8
2.2 Bisection Method . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Newton-Raphson Method . . . . . . . . . . . . . . . . . . . . 11
2.4 Regula-Falsi Method . . . . . . . . . . . . . . . . . . . . . . . 13
2.5 Secant Method . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.6 The Order and the Rate of Convergence . . . . . . . . . . . . 14

3 Simultaneous linear equations 17

3.1 Gauss-Jacobi Method . . . . . . . . . . . . . . . . . . . . . . . 18
3.2 Gauss-Seidel Method . . . . . . . . . . . . . . . . . . . . . . . 19
3.3 Gauss Elimination Method with Partial Pivoting . . . . . . . . 20

4 Interpolation 23
4.1 Lagrange Interpolation . . . . . . . . . . . . . . . . . . . . . . 23
4.1.1 Lagrange Interpolating Polynomial . . . . . . . . . . . 23
4.2 Newton’s Form of Interpolation . . . . . . . . . . . . . . . . . 25
4.2.1 Newton’s Divided Difference Formula . . . . . . . . . . 25
4.2.2 Divided Differences . . . . . . . . . . . . . . . . . . . . 26
4.2.3 Newton’s Interpolating Polynomial . . . . . . . . . . . 26
4.3 Finite Difference Operators . . . . . . . . . . . . . . . . . . . 27
4.3.1 Shift Operator E . . . . . . . . . . . . . . . . . . . . . 28
4.3.2 Forward Difference Operator (∆) . . . . . . . . . . . . 28
4.3.3 Backward Difference Operator (∇) . . . . . . . . . . . 28

i
 4.3.4 Central Difference Operator (δ) . . . . . . . . . . . . . 29
4.3.5 Average Difference Operator (µ) . . . . . . . . . . . . . 29

5 Differentiation, Integration and ODE 31

5.1 Differentiation . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
5.1.1 Forward Difference Derivative . . . . . . . . . . . . . . 31
5.1.2 Backward Difference Derivative . . . . . . . . . . . . . 31
5.1.3 Central Difference Derivative . . . . . . . . . . . . . . 32
5.1.4 Three-Point Forward Difference Derivative . . . . . . . 32
5.1.5 Three-Point Backward Difference Derivative . . . . . . 32
5.1.6 General three-point formula for derivative . . . . . . . 33
5.1.7 Formula for Second-Order Derivative . . . . . . . . . . 33
5.2 Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
5.2.1 Trapezoid Rule (Composite form) . . . . . . . . . . . . 34
5.2.2 Simpson’s 1/3rd Rule . . . . . . . . . . . . . . . . . . . 36
5.3 Ordinary Differential Equation . . . . . . . . . . . . . . . . . . 38
5.3.1 Euler’s Method . . . . . . . . . . . . . . . . . . . . . . 39
5.3.2 Second-Order Runge-Kutta Methods (RK2) . . . . . . 39
5.3.3 Special Cases of 2nd-Order Runge-Kutta Methods . . . 40
5.4 Classical fourth-Order Runge-Kutta Method (RK4) . . . . . . 41
Chapter 1

Preliminaries

1.1 Significant Digits

In numerical computation and measurement, precision plays a vital role.
One way of expressing precision is through significant digits (also called
significant figures). These are the digits in a number that carry meaningful
information about its accuracy. The more significant digits a number has,
the higher its precision, whether it comes from a measured quantity or a
computed value.
For example, the number 2.45 is more precise than 2.4, since it conveys
information up to the hundredths place. Similarly, reporting the value of π
as 3.1416 communicates greater accuracy than simply writing 3.14.

Rules for Identifying Significant Digits

To determine how many significant digits a number contains, certain rules
are followed. These rules ensure consistency when reporting results and com-
paring numerical values.

1. Non-zero digits are always significant. Every non-zero digit in a

number contributes to its precision.
Examples: 245 (3 significant digits), 2.45 (3 significant digits).

2. Zeros between non-zero digits are significant. Any zero that lies
between two non-zero digits adds to the count of significant digits.
Example: 2501005 (7 significant digits).

3. Leading zeros are not significant. Zeros placed to the left of the
first non-zero digit (with a decimal point) merely act as placeholders.

1
2 CHAPTER 1. PRELIMINARIES

They do not contribute to the precision.

Example: 0.00456 (3 significant digits).

4. Trailing zeros are significant if a decimal point is present.

When a number contains a decimal point, the zeros at the end are
meaningful and increase the precision.
Example: 12.300 (5 significant digits).

5. Trailing zeros without a decimal point are ambiguous. In whole

numbers where no decimal point is shown, trailing zeros may or may
not be significant. The interpretation often depends on the context,
and scientific notation is preferred for clarity.
Example: 1200 may be taken as: 1.2 × 103 (2 significant digits), 1.20 ×
103 (3 significant digits), 1.200 × 103 (4 significant digits).

Understanding these rules helps ensure that numerical results are commu-
nicated with the appropriate degree of accuracy, avoiding both overstatement
and understatement of precision.

1.1.1 Significant Digits in Rounding Off

In practical computation, it is often necessary to round numbers to a spec-
ified number of significant digits or decimal places. Rounding ensures that
results are reported with an appropriate level of precision without carrying
unnecessary digits. The process follows a few simple rules:

1. If the (n + 1)-th digit is less than 5: Keep the n-th digit unchanged
and simply discard all digits to the right.
Example: Consider 3.1432 rounded to 3 decimal places. The 4th digit
is 2 (less than 5). Hence, the rounded value is 3.143.

2. If the (n + 1)-th digit is greater than 5: Increase the n-th digit by

1 and discard all digits to the right.
Example: For 3.1467 rounded to 3 decimal places, the 4th digit is 7
(greater than 5). Therefore, the rounded value is 3.147.

3. If the (n + 1)-th digit is exactly 5: This is a special case, often han-

dled by the “round half to even” rule (also called banker’s rounding):

• If the n-th digit is odd, increase it by 1.

• If the n-th digit is even, leave it unchanged.
1.2. ERRORS IN NUMERICAL METHODS 3

In both cases, all digits to the right are discarded.

Example 1: For 3.145 rounded to 2 decimal places, the 3rd digit is 5
and the 2nd digit is 4 (even). Hence, the result is 3.14.
Example 2: For 3.135 rounded to 2 decimal places, the 3rd digit is 5
and the 2nd digit is 3 (odd). Therefore, the result becomes 3.14.
These rules help maintain fairness and accuracy in rounding, especially
in large computations, by preventing bias that might occur if numbers were
always rounded up when the (n + 1)-th digit is 5.

1.2 Errors in Numerical Methods

In mathematics, we often come across equations that need to be solved for
practical and theoretical purposes. Some equations can be solved exactly
using algebraic formulas, while others do not yield to such straightforward
methods.
For√instance, consider the quadratic equation x2 − 2 = 0. Its solutions
are ± 2, which can either be identified directly or obtained through the
quadratic formula you learned in high school. However, equations such as
xex − cos x = 0 cannot be solved using simple algebraic techniques. In such
situations, we resort to numerical methods, which provide approximate
solutions through iterative processes.
Numerical methods generate a sequence of approximations that ideally
converge to the exact solution. Naturally, these approximations are not per-
fect, and the discrepancy between the true solution and the computed ap-
proximation is referred to as the error. Formally, we define error as
Error = Exact solution − Approximate solution.
A natural question arises: Why do errors occur in numerical compu-
tation? The sources of these errors lie mainly in three areas:
1. Approximations in mathematical models – Simplifications are of-
ten introduced when formulating real-world problems into mathemati-
cal equations.
2. Inherent limitations of numerical algorithms – Many algorithms
rely on approximations, iterations, or truncations, all of which intro-
duce errors.
3. Finite precision of computers – Digital computers represent real
numbers with a limited number of bits, which restricts accuracy and
leads to rounding effects.
4 CHAPTER 1. PRELIMINARIES

1.2.1 Types of Errors

Errors encountered in numerical methods can be categorized in several ways.
The following are the most important types:

Absolute Error
The absolute error measures the magnitude of the difference between the
exact value (x) and its approximation (x∗ ). It is
√ defined as EA = |x − x∗ |.
Example: The exact solution of x2 −3 = 0 is 3. Suppose an approximate
solution is√given as 1.731. Then, the absolute error (to three decimal places)
is EA = | 3 − 1.731| = 0.001.

Relative Error
The relative error provides a scale-free measure of accuracy by comparing
|x−x∗ |
the error to the true value: ER = |x| .
Example: Using the previous case, the relative error is ER = 0.001
√ . A
3
smaller value of ER indicates a better approximation.

Percentage Error
The percentage error is simply the relative error expressed as a percentage:
EP = ER ×100. This measure is often used in practical applications to express
error in an intuitive way.

Round-off Error
A round-off error arises because computers cannot represent most real
numbers exactly. Instead, they store numbers using floating-point represen-
tation with finite precision.
Example: The mathematical constant π has the value 3.14159265358979 . . . .
In a computer system, it may be stored as 3.1416. This introduces a round-off
error due to the truncation of digits beyond machine precision.

Truncation Error
A truncation error occurs when an infinite process, such as a series expan-
sion or iterative algorithm, is terminated after a finite number of steps.
Example: The exponential function ex can be expanded as a Taylor series:
2 3
e = 1 + x + x2! + x3! + · · ·
x

If we approximate ex by retaining only the first three terms, then the

neglected higher-order terms introduce a truncation error.
1.3. EXERCISE 5

1.3 Exercise
1. How many significant digits are there in the number 0.007230?

2. Identify the number of significant digits in the following measurement:

4.3050 × 103 .

3. Round the number 12.6759 to 4 significant digits.

4. How many significant digits are present in 0.0020500?

5. Round 27.8571 to 3 significant digits.

6. Calculate the number of significant digits in the value 506000 when it

is written in scientific notation as 5.06 × 105 .

7. Round 0.003056 to 2 significant digits.

8. How many significant digits are in the number 700000? What changes
if it is written as 7.0 × 105 ?

9. Round 145.25 to 2 significant digits and explain the rule applied for
rounding off.

10. The length of an object is measured as 8.3500 cm. How many significant
digits does this measurement have, and what would be the result if it
were rounded to 3 significant digits?

11. If the true value of e is 2.7182818 and its approximate value is given
by 2.7200000, find the absolute and relative errors.
√
12. If the true value of 2 is 1.4142136 and its approximate value is given
by 1.4140000, find the absolute and relative errors.

13. If the true value of the gravitational constant G is 6.67430 × 10−11 and
its approximate value is given by 6.67000 × 10−11 , find the absolute and
relative errors.

14. If the true value of the speed of light c is 299, 792, 458 meters per second
and its approximate value is given by 299, 800, 000 meters per second,
find the absolute and relative errors.

15. If the true value of the Avogadro constant NA is 6.02214076 × 1023 and
its approximate value is given by 6.02000000 × 1023 , find the absolute
and relative errors.
6 CHAPTER 1. PRELIMINARIES
Chapter 2

Roots of Transcendental and

Polynomial Equations

Finding the roots of an equation is one of the most common problems in

mathematics. Broadly, equations can be classified into two categories: alge-
braic equations and transcendental equations.
An algebraic equation is one that can be expressed in terms of polynomi-
als. For example,
3 2 x2 − 2
x − 2x + 5x − 3 = 0, = 5x + 1
x3 + 3x2 − 7
are algebraic equations, since they can ultimately be written in polynomial
form.
On the other hand, a transcendental equation involves transcendental
functions such as trigonometric, logarithmic, or exponential functions. Ex-
amples include
cos x = x, ex = x2 , xex = log(sin x).
For lower-degree polynomials, exact methods are available to determine
the roots. For instance, quadratic equations can be solved using the quadratic
formula. However, when dealing with higher-degree polynomials or transcen-
dental equations, no general formula exists to determine the roots analyti-
cally.
To overcome this limitation, we rely on numerical methods, which pro-
vide approximate solutions with desired accuracy. In this chapter, we will
explore various numerical techniques to find the roots of both algebraic and
transcendental equations.
Numerical methods can be applied only when the given equation is known
to have a solution in a specified interval. This is ensured by the following
result, based on the Intermediate Value Theorem.

7
8CHAPTER 2. ROOTS OF TRANSCENDENTAL AND POLYNOMIAL EQUATIONS

Theorem 2.1 Let f be a continuous real-valued function on [a, b]. If f (a) ·

f (b) < 0, then there exists at least one root of f (x) in (a, b).

From the graph, we observe that at x = a, the function value is negative,

while at x = b, it is positive. Hence, the product f (a) · f (b) is negative. Since
f (x) is continuous on [a, b], the curve must cross the x-axis at some point
c ∈ (a, b), which represents the required root of the equation f (x) = 0.
Let us now focus on the numerical methods for finding the root of the
equation f (x) = 0. These methods yield approximate solutions by generating
a sequence of values through suitable recursive relations, and are especially
useful in cases where analytical techniques are either unavailable or imprac-
tical to apply.

2.1 Fixed-Point Iteration Method

A real number p is called a fixed point of a function g if it satisfies g(p) = p.
In the fixed-point iteration method, the given equation f (x) = 0 is rewritten
in the form g(x) = x, where g(x) is continuous in the specified interval [a, b].
The method then seeks a fixed point of g(x). Since any fixed point of g also
serves as a solution of the original equation f (x) = 0, this approach yields
the required root.
Beginning with an initial approximation x0 , we construct a sequence {xn }
of successive approximations to the root by using the iterative relation

xn+1 = g(xn ), n ≥ 0.
2.1. FIXED-POINT ITERATION METHOD 9

Convergence of Fix-Point Iteration Method

As discussed earlier, the numerical methods generate a sequence of approx-

imate solutions which, intuitively, should converge to the exact root of the
function f (x). By convergence, we mean that as the number of iterations
increases, the approximations move closer to the true root of the equation
f (x) = 0. The natural question, however, is: how can we ensure that the
sequence we generate will indeed converge to the exact root? For this pur-
pose, we require certain convergence criteria for the sequence {xn }.
In the Fixed-Point Iteration Method, the convergence criterion is

|g ′ (x)| < 1, x ∈ (a, b).

To see this, suppose |g ′ (x)| ≤ k < 1 and let p be the exact fix-point of g.
Then, using the definition of the derivative of g(x) at the fixed point p, we
have
g(xn ) − g(p)
|g ′ (p)| = .
xn − p
Hence,

|xn+1 − p| ≤ k|xn − p| ≤ k 2 |xn−1 − p| ≤ · · · ≤ k n |x0 − p|.

Here, |x0 − p| is some fixed positive number, and k is a real constant strictly
less than 1. Therefore, as n → ∞, we obtain

|xn+1 − p| → 0,

which shows that the sequence {xn } converges to the fixed point p. Thus, if
the derivative of g(x) is bounded by 1 in the interval of interest, the sequence
of approximations will converge to the exact solution.

Example 2.1 Find an approximate solution of the equation cos(x) − x = 0

in the interval [0, 1] using Fix-Point Iteration Method.

Solution. We first rewrite the given equation in the form g(x) = x. Here,
g(x) = cos(x). Since the derivative of cos(x) has absolute value less than 1
in the interval (0, 1), the fixed-point iteration is applicable. Let the initial
guess be x0 = 0.5. Using the iteration formula xn+1 = g(xn ), we obtain
x1 = g(x0 ) = cos(0.5) = 0.877582, x2 = g(x1 ) = cos(0.877582) = 0.639012,
and so on. The first ten iterations are summarized in Table 2.1.
10CHAPTER 2. ROOTS OF TRANSCENDENTAL AND POLYNOMIAL EQUATIONS

Iteration n Approximation xn
0 0.500000
1 0.877583
2 0.639012
3 0.802685
4 0.694778
5 0.768196
6 0.719165
7 0.752356
8 0.730081
9 0.745120
10 0.735006

Table 2.1: Iterations for Fixed-Point Iteration Method for f (x) = cos(x) − x

2.2 Bisection Method

Another technique to find the root of an equation f (x) = 0 is the Bisection
Method. It is based on the Intermediate Value Theorem (IVT), which states
that if a continuous function f (x) satisfies the condition f (a) · f (b) < 0, then
there exists at least one root of f (x) in the interval [a, b].
In this method, the interval [a, b] is repeatedly halved, and the subinterval
satisfying the Intermediate Value Theorem is chosen, thereby narrowing down
the root’s location. Let the initial interval be [a0 , b0 ] with a0 = a and b0 = b.
The next approximation of the root is the midpoint
an−1 + bn−1
mn = .
2
The new interval [an , bn ] is then either [an−1 , mn ] or [mn , bn−1 ], depending on
where the IVT holds for f (x). Repeating this process generates a sequence
{mn } of midpoints.

Convergence criteria
At each iteration, the interval is halved. Let the length of the interval at the
nth iteration be ln . Then,
l0
ln = |bn − an | ≤ 12 |bn−1 − an−1 | ≤ · · · ≤ 1
|b
2n 0
− a0 | = 2n
.
This shows that as n increases, ln decreases. Hence, with more iterations, the
midpoint mn approaches the interval’s endpoints. Since the IVT guarantees
that a root of f lies within the interval, the sequence of midpoints {mn }
converges to the solution of f (x) = 0.
2.3. NEWTON-RAPHSON METHOD 11

Example 2.2 Find an approximate solution of the equation cos(x) − x = 0

in the interval [0, 1] using Bisection Method.

Solution. Consider f (x) = cos(x)−x with the initial interval [a0 , b0 ] = [0, 1].
Since f (0) · f (1) < 0, the IVT ensures that a root of f lies within [0, 1]. The
first midpoint is
a0 + b 0
m1 = = 0.5.
2

Here, f (0) · f (0.5) > 0 and f (0.5) · f (1) < 0, so the next interval becomes
[0.5, 1]. Proceeding in this manner, the first 5 approximations are shown in
Table 2.2.

an−1 +bn−1
Iteration n mn = 2
Sign of f (mn ) an bn
0 0 1
1 0.5 + 0.5 1
2 0.75 − 0.5 0.75
3 0.625 + 0.625 0.75
4 0.6857 + 0.685 0.75
5 0.71875 + 0.71875 0.75

Table 2.2: Iterations in Bisection Method for f (x) = cos(x) − x

2.3 Newton-Raphson Method

This method relies on linear approximation and makes use of the derivative of
the function to compute successive estimates. The next approximation xn+1
is obtained as the x-intercept of the tangent drawn at the point (xn , f (xn )) on
the curve f (as illustrated in the figure below). The update formula follows
from the equation of the tangent line:

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

where f ′ (xn ) denotes the derivative of f (x) at x = xn .

12CHAPTER 2. ROOTS OF TRANSCENDENTAL AND POLYNOMIAL EQUATIONS

Conditions for Convergence

The Newton-Raphson Method converges if the function f (x) is differentiable
near the root, the initial guess x0 is sufficiently close to the actual root, and
the derivative f ′ (x) ̸= 0 near the root.

Example 2.3 Find an approximate solution of the equation cos(x) − x = 0

in the interval [0, 1] using Newton-Rapshon Method.

Solution. We now illustrate the Newton-Raphson Method by applying it to

the function
f (x) = cos(x) − x.
In order to proceed, we first compute the derivative of f (x), which is

f ′ (x) = − sin(x) − 1.

Hence, the Newton-Raphson iteration formula takes the form

f (xn ) cos(xn ) − xn
xn+1 = xn − ′
= xn + .
f (xn ) sin(xn ) + 1
Starting with the initial approximation x0 = 0.5, we obtain
cos(0.5) − 0.5
x1 = 0.5 + ≈ 0.725536.
sin(0.5) + 1
Successive approximations obtained by this method are displayed in Table
2.3.
2.4. REGULA-FALSI METHOD 13

Iteration n xn
0 0.500000
1 0.725536
2 0.738369
3 0.739085
4 0.739085

Table 2.3: Iterations using Newton-Raphson Method for f (x) = cos(x) − x

2.4 Regula-Falsi Method

The procedure begins with two initial points x0 and x1 such that f (x0 ) ·
f (x1 ) < 0. By the Intermediate Value Theorem (as in the Bisection Method),
this guarantees that at least one root lies within the interval [x0 , x1 ].
The method makes use of linear interpolation. The next approximation x2
is obtained as the point of intersection of the straight line joining (x0 , f (x0 ))
and (x1 , f (x1 )) with the x-axis. The formula is

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

Once x2 is computed, we determine whether f (x0 ) · f (x2 ) < 0 or f (x2 ) ·

f (x1 ) < 0 to identify the new interval containing the root. For clarity, let
us denote the initial interval as [a0 , b0 ], where a0 = x0 and b0 = x1 . In this
notation, the formula becomes

f (b0 )(b0 − a0 )
x 2 = b0 − ,
f (b0 ) − f (a0 )

and the next interval is chosen as either [a1 , b1 ] = [x0 , x2 ] or [a1 , b1 ] = [x2 , x1 ],
depending on the sign condition.
In general, the iterative step of the method is expressed as

f (bn )(bn − an )
xn+1 = bn − .
f (bn ) − f (an )

The following Table 2.4 presents a few iterations for the function f (x) =
cos(x) − x with initial values x0 = 0 and x1 = 1, computed using the Regula-
Falsi Method.
14CHAPTER 2. ROOTS OF TRANSCENDENTAL AND POLYNOMIAL EQUATIONS

Iteration n xn+1 Sign of f (xn+1 ) an bn

0 0 1
1 0.68507 + 0.68507 1
2 0.73630 + 0.73630 1
3 0.73895 + 0.73895 1
4 0.73908 + 0.73908 1
5 0.73908 + 0.73908 1

Table 2.4: Iterations in Regula-Falsi Method for f (x) = cos(x) − x

2.5 Secant Method

The Secant Method is closely related to the Regula-Falsi Method and employs
the same iterative formula,

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

Unlike Regula-Falsi, however, the Secant Method does not preserve a bracket-
ing interval. Instead, it uses the two most recent approximations to generate
the next one. This often yields faster convergence, though at the cost of
losing the guaranteed bracketing of the root, which may sometimes result in
divergence.

2.6 The Order and the Rate of Convergence

In these iterative methods, we generate a sequence {xn } of approximate so-
lutions to the equation f (x) = 0. Under suitable conditions, this sequence
converges to the exact solution of the equation. A natural question, how-
ever, is how rapidly this convergence takes place and how such speed can
be measured. The rate at which the iterates approach the exact solution
is characterized by two notions: the rate of convergence and the order of
convergence.
Before defining these terms, let us first recall the meaning of convergence.

Definition 2.1 (Limit of a Sequence) A sequence {xn } is said to con-

verge to a real number L if
lim xn = L,
n→∞

or equivalently,
lim |xn − L| = 0.
n→∞
2.6. THE ORDER AND THE RATE OF CONVERGENCE 15

The number L is called the limit of the sequence {xn }. In simple words, as
n becomes larger, the values of xn get closer and closer to L.

Now, let x∗ denote the exact solution to the equation f (x) = 0, and let
{xn } be the sequence of approximations. The error at the nth iteration is
defined as
en := |xn − x∗ |.
With this notation, we are now ready to define the rate and the order of
convergence.

Definition 2.2 Let {xn } be a sequence converging to x∗ . If there exists a

constant C > 0 and a real number p ≥ 1 such that p is the largest real number
for which
|xn+1 − x∗ | en+1
lim = lim p = C,
n→∞ |xn − x∗ |p n→∞ en

then the number p is called the order of convergence, and the constant C
is called the asymptotic error constant. The parameter p characterizes the
efficiency of the iterative method, while C provides a measure of the con-
vergence rate. Larger values of p (and smaller values of C) indicate faster
convergence.

n+2 n+2
Example 2.4 Let xn = . Then limn→∞ xn = limn→∞ = 1, so
n+5 n+5
x∗ = 1. Therefore,

n+2 3
|xn − x∗ | = en = −1 =
n+5 n+5

and
3
en+1 = .
n+6
We compute

en+1 (n + 5)p np−1 (1 + 5/n)p

lim p = lim = lim .
n→∞ en n→∞ 3p−1 (n + 6) n→∞ 3p−1 (1 + 6/n)

If p > 1, the limit diverges since np−1 → ∞. Thus the largest admissible
en+1 n+5
p is 1, and limn→∞ = limn→∞ = 1. Therefore, the order of
en n+6
convergence is 1, with the rate of convergence 1.
16CHAPTER 2. ROOTS OF TRANSCENDENTAL AND POLYNOMIAL EQUATIONS

Example 2.5 Let f (x) = x2 − 2, and we apply Newton-Rapshon

√ Method to
approximate its root. In theory the root of f (x) is 2. The Newton-Rapshon
method gives the following sequence:

x2n 1
xn+1 = + .
2 xn
Here √
√ xn 1 √ (xn − 2)2
|xn+1 − 2| = + − 2= ,
2 xn 2xn
which gives
en+1 1
2
=
en 2xn
So, the limit
en+1 1
lim 2
= √ ,
n→∞ en 2 2
is finite. Thus, the order of convergence is 2 with Asymptotic Error Constant
(rate of convergence) 2√1 2 .

Remark 2.1 The rate of convergence describes the speed at which the se-
quence of errors en approaches zero, whereas the order of convergence reflects
the efficiency of an iterative method. Put differently, the rate of convergence
is a quantitative measure, while the order of convergence is a qualitative mea-
sure.
It is worth noting that not all texts distinguish clearly between these two
notions—some treat them as the same, while others emphasize the difference.
Formally, the rate of convergence is often understood as the limit of the
ratio of the error at the next iteration to that at the current iteration. For a
linearly convergent method, this limit exists as a nonzero constant. However,
for methods of higher order (p > 1), this limit is zero, and hence the definition
loses its meaning in that form. To address this, many books instead define
the asymptotic error constant as the rate of convergence.
Chapter 3

Solution of simultaneous linear

equations

In the previous chapter, we learned how to find the solution of an equa-

tion of the form f (x) = 0, where f (x) may be an algebraic function or a
transcendental function.
In this chapter, we turn our attention to solving a system of linear equa-
tions using numerical methods. In general, an n-dimensional system of linear
equations can be written as
a11 x1 + a12 x2 + · · · + a1n xn = b1 ,
a21 x1 + a22 x2 + · · · + a2n xn = b2 ,
.. ..
. .
an1 x1 + an2 x2 + · · · + ann xn = bn .
Here, x1 , x2 , . . . , xn are the variables, and aij , 1 ≤ i, j ≤ n, are constants.
This system can also be represented in matrix form as
    
a11 a12 · · · a1n x1 b1
 a21 a22 · · · a2n   x2   b2 
..   ..  =  ..  (3.1)
    
 .. .. . .
 . . . .  .   . 
an1 an2 · · · ann xn bn
or, in compact notation,
Ax = b,
where A is the n × n coefficient matrix, and x and b are n × 1 column
matrices.
When n is small, such systems can be solved directly without much dif-
ficulty. However, as n becomes large, solving them by conventional methods
becomes complicated and computationally expensive.

17
18 CHAPTER 3. SIMULTANEOUS LINEAR EQUATIONS

3.1 Gauss-Jacobi Method

The Gauss-Jacobi method is an iterative technique for solving a system of
linear equations such as (5.1). Instead of attempting to solve the equations
directly, this method improves upon an initial guess step by step until the
solution is reached with the desired accuracy.
The general iterative formula is given by:
 
n
(k+1) 1  X (k) 

xi =  b i − a ij x j , k = 0, 1, 2, . . . , 1 ≤ i, j ≤ n,
aii  j=1
j̸=i

where:
(k+1)
• xi represents the new (updated) value of the i-th variable at the
(k + 1)-th iteration.
(k)
• xj is the value of the j-th variable obtained in the k-th iteration.

• aii denotes the diagonal element of the coefficient matrix A correspond-

ing to row i.

• aij is the element of A in row i and column j.

• bi is the i-th entry of the constant vector b.

Convergence Criteria
For the Gauss-Jacobi method to be effective, the coefficient matrix A should
preferably be diagonally dominant. A matrix is diagonally dominant if
n
X
|aii | > |aij |, for all i.
j=1
j̸=i

When this condition holds, convergence is usually guaranteed. If not, the

method may converge very slowly or fail altogether.

Example 3.1 Example. Consider the system:


10x − 2y + z = 7,

−x + 8y − 2z = −4,

2x − y + 9z = 12.

3.2. GAUSS-SEIDEL METHOD 19

To apply the Gauss-Jacobi method, we first rewrite each equation in a

form that isolates the variable on the left-hand side:
1 1 1
x= (7 + 2y − z) , y = (−4 + x + 2z) , z= (12 − 2x + y) .
10 8 9
Starting with the initial guess
x(0) = 0, y (0) = 0, z (0) = 0,
we update the values iteratively. The first few approximations are shown
below:

Iteration x y z
0 0.00000 0.00000 0.00000
1 0.70000 −0.50000 1.3333
2 0.46667 −0.07917 1.1222
3 0.57194 −0.16111 1.2208
4 0.54569 −0.12330 1.1883
5 0.55651 −0.13470 1.1984
6 0.55322 −0.13084 1.1947
7 0.55436 −0.13217 1.1959
8 0.55398 −0.13174 1.1955
9 0.55411 −0.13189 1.1956
10 0.55406 −0.13184 1.1955
After about ten iterations, the solution stabilizes at:
x ≈ 0.55406, y ≈ −0.13184, z ≈ 1.1955.
This demonstrates how the Gauss-Jacobi method gradually refines the values
of the unknowns until convergence is achieved.

3.2 Gauss-Seidel Method

The Gauss-Seidel method is an iterative technique for solving a system of
linear equations(5.1). It can be regarded as an improvement over the Gauss-
Jacobi method. The difference lies in the fact that Gauss-Seidel uses the
newly computed values of the variables immediately, which usually makes
the method converge faster. The convergence criteria is similar to that of
Gauss-Jacobi method.
The iterative formula is given by:
i−1 n
!
(k+1) 1 X (k+1)
X (k)
xi = bi − aij xj − aij xj , k = 0, 1, . . . , 1 ≤ i, j ≤ n.
aii j=1 j=i+1
20 CHAPTER 3. SIMULTANEOUS LINEAR EQUATIONS

Example 3.2 Example. Consider the system:


10x − 2y + z = 7,

−x + 8y − 2z = −4,

2x − y + 9z = 12.


Rewriting each equation to isolate the corresponding variable:

1 1 1
x= (7 + 2y − z) , y= (−4 + x + 2z) , z= (12 − 2x + y) .
10 8 9
Starting with the initial guess
x(0) = 0, y (0) = 0, z (0) = 0,
we iteratively compute new values. The approximations are tabulated below:

Iteration x y z
0 0.000000 0.000000 0.000000
1 0.700000 −0.412500 1.425278
2 0.505972 −0.135366 1.189860
3 0.551184 −0.131721 1.195349
4 0.554043 −0.131849 1.195593
5 0.554065 −0.131850 1.195598
6 0.554067 −0.131850 1.195598
7 0.554067 −0.131850 1.195598
After a few iterations, the solution stabilizes at:
x ≈ 0.554067, y ≈ −0.131850, z ≈ 1.195598.

3.3 Gauss Elimination Method with Partial

Pivoting
The Gauss elimination method is a systematic procedure for solving systems
of linear equations. The central idea is to transform the given system into
an upper triangular form, from which the unknowns can be determined
easily by back substitution.
However, in practical computation, very small pivot elements can cause
large numerical errors. To reduce this risk, we often use a technique known
as partial pivoting, where rows are interchanged so that the largest available
element in a column is chosen as the pivot.
Steps of the method:
3.3. GAUSS ELIMINATION METHOD WITH PARTIAL PIVOTING 21

• Start with the first column and choose the pivot element. If needed,
interchange rows so that the largest element (in magnitude) is placed
in the pivot position.
• Eliminate all entries below the pivot using row operations.
• Repeat the process column by column until the coefficient matrix is in
upper triangular form.
• Solve the resulting triangular system by back substitution.

Example 1. Consider the system

x1 + 2x2 + 3x3 = 9,
2x1 + 3x2 + 4x3 = 13,
3x1 + 4x2 + x3 = 10.
In matrix form, this can be written as
   
1 2 3 9
A = 2 3 4 , b = 13 .
3 4 1 10
In the first column, the largest element is 3 (in row 3). By swapping row
1 and row 3, we place this value in the pivot position. We then eliminate
entries below the pivot and continue the process column by column, applying
pivoting whenever necessary. Finally, the system is transformed into upper
triangular form and solved by back substitution.
Why do we need partial pivoting? The importance of partial pivoting
becomes clear when small pivot elements are present. For instance, consider
the system
0.001x + y = 3,
x + 2y = 5.
If we proceed without pivoting, the small coefficient 0.001 in the first
equation leads to significant numerical error in the elimination process. The
solution obtained, rounding off to 2 decimal digits, is x = 0, y = 3, which is
not very accurate.
On the other hand, if we apply partial pivoting by interchanging the two
equations first, the pivot element becomes 1 instead of 0.001. Eliminating as
before, the system gives the more accurate solution x = −1, y = 3.
This shows that partial pivoting plays a crucial role in ensuring numerical
stability and reliability of the Gauss elimination method, especially when the
coefficient matrix contains very small or very large values.
22 CHAPTER 3. SIMULTANEOUS LINEAR EQUATIONS
Chapter 4

Interpolation

4.1 Lagrange Interpolation

Lagrange interpolation is a polynomial interpolation method used to estimate
the value of a function given a set of data points. The idea is to construct
a polynomial that passes through all the given points. This method is espe-
cially useful when you have a small set of known data points and want to
estimate values in between them.

4.1.1 Lagrange Interpolating Polynomial

Given a set of n + 1 data points (x0 , y0 ), (x1 , y1 ), . . . , (xn , yn ), the Lagrange
interpolating polynomial P (x) is defined as follows:

n
X
P (x) = yi Li (x),
i=0

where Li (x) are the Lagrange basis polynomials given by:

Y x − xj
Li (x) = .
0≤j≤n
x i − x j
j̸=i

Explanation

The Lagrange basis polynomial Li (x) has the property that it equals 1 at
x = xi and 0 at all other xj (where j ̸= i). This property ensures that the
interpolating polynomial P (x) passes through each of the given data points.

23
24 CHAPTER 4. INTERPOLATION

Steps to Construct the Lagrange Interpolating Polyno-

mial
1. Identify the Data Points: Determine the set of known data points
(x0 , y0 ), (x1 , y1 ), . . . , (xn , yn ).

2. Construct Lagrange Basis Polynomials: For each data point i,

construct the Lagrange basis polynomial Li (x) using the formula:
Y x − xj
Li (x) = .
0≤j≤n
x i − x j
j̸=i

3. Form the Interpolating Polynomial: Combine the basis polynomi-

als to form the interpolating polynomial:
n
X
P (x) = yi Li (x).
i=0

4. Simplify (if possible): Sometimes, the resulting polynomial can be

simplified for easier computation or analysis.

Example 1
Suppose you have three data points: (1, 2), (2, 3), and (4, 5). The Lagrange
interpolating polynomial P (x) would be:

(x − 2)(x − 4) (x − 2)(x − 4)
L0 (x) = = ,
(1 − 2)(1 − 4) 3
(x − 1)(x − 4) (x − 1)(x − 4)
L1 (x) = =− ,
(2 − 1)(2 − 4) 2
(x − 1)(x − 2) (x − 1)(x − 2)
L2 (x) = = .
(4 − 1)(4 − 2) 6

Then, the polynomial P (x) is:

P (x) = 2L0 (x) + 3L1 (x) + 5L2 (x).

4.2. NEWTON’S FORM OF INTERPOLATION 25

Example 2
Given the data points (0, 1), (1, 2), and (2, 0), the Lagrange interpolating
polynomial P (x) is:
2
X
P (x) = yi Li (x),
i=0

where the Lagrange basis polynomials are:

(x − 1)(x − 2) (x − 1)(x − 2)
L0 (x) = = ,
(0 − 1)(0 − 2) 2
(x − 0)(x − 2)
L1 (x) = = −(x)(x − 2),
(1 − 0)(1 − 2)
(x − 0)(x − 1) (x)(x − 1)
L2 (x) = = .
(2 − 0)(2 − 1) 2
Substituting the values into P (x):

P (x) = 1 · L0 (x) + 2 · L1 (x) + 0 · L2 (x),

(x − 1)(x − 2)
P (x) = − 2x(x − 2).
2
Simplifying this expression gives the quadratic polynomial:
3 5
P (x) = − x2 + x + 1.
2 2

4.2 Newton’s Form of Interpolation

Newton’s interpolation provides an efficient way to compute the interpolating
polynomial for a given set of data points. It uses divided differences to
construct the polynomial incrementally, making it useful when adding new
data points.

4.2.1 Newton’s Divided Difference Formula

Given a set of n + 1 distinct data points (x0 , y0 ), (x1 , y1 ), . . . , (xn , yn ), the
interpolating polynomial in Newton’s form is expressed as:

P (x) = a0 +a1 (x−x0 )+a2 (x−x0 )(x−x1 )+· · ·+an (x−x0 )(x−x1 ) . . . (x−xn−1 ),
26 CHAPTER 4. INTERPOLATION

where the coefficients a0 , a1 , . . . , an are determined using divided differences.

4.2.2 Divided Differences

Divided differences are computed recursively based on the values of the data
points. The first divided difference is simply the value yi :

f [xi ] = yi .
The first-order divided difference for two points xi and xi+1 is:

f [xi+1 ] − f [xi ]
f [xi , xi+1 ] = .
xi+1 − xi
The second-order divided difference for three points xi , xi+1 , xi+2 is:

f [xi+1 , xi+2 ] − f [xi , xi+1 ]

f [xi , xi+1 , xi+2 ] = .
xi+2 − xi
In general, the k-th order divided difference for k+1 points xi , xi+1 , . . . , xi+k
is:

f [xi+1 , . . . , xi+k ] − f [xi , . . . , xi+k−1 ]

f [xi , xi+1 , . . . , xi+k ] = .
xi+k − xi

4.2.3 Newton’s Interpolating Polynomial

Using the divided differences, the polynomial P (x) can be written as:

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

· · · + f [x0 , x1 , . . . , xn ](x − x0 )(x − x1 ) . . . (x − xn−1 ).

Each term in the polynomial uses an additional factor of (x − xi ) as the

degree increases.
4.3. FINITE DIFFERENCE OPERATORS 27

Example
Consider the following data points:

(x0 , y0 ) = (1, 1), (x1 , y1 ) = (2, 4), (x2 , y2 ) = (3, 9), (x3 , y3 ) = (4, 16).

We will now compute the divided differences in tabular form.

xi yi = f [xi ] f [xi , xi+1 ] f [xi , xi+1 , xi+2 ] f [xi , xi+1 , xi+2 , xi+3 ]
1 1 4−1
2−1
=3 5−3
2 4 3−1
=1 0−0
4−1
=0
3 9 9−4
3−2
=5
4 16

Explanation
• The first column contains the values of xi .
• The second column contains the values of yi = f (xi ), which are the
given data points.
• The third column contains the first-order divided differences:
4−1 9−4 16 − 9
f [x0 , x1 ] = = 3, f [x1 , x2 ] = = 5, f [x2 , x3 ] = = 7.
2−1 3−2 4−3

• The fourth column contains the second-order divided differences:

5−3 7−5
f [x0 , x1 , x2 ] = = 1, f [x1 , x2 , x3 ] = = 1.
3−1 4−2

• The fifth column contains the third-order divided difference:

1−1
f [x0 , x1 , x2 , x3 ] = = 0.
4−1

Thus, the Newton interpolating polynomial is:

P (x) = 1 + 3(x − 1) + 1(x − 1)(x − 2) = x2 .

4.3 Finite Difference Operators

In this section, we will discuss various finite difference operators that are
widely used in numerical analysis and discrete calculus. Let the tabular
points x0 , x1 , . . . , xn be equally spaced, that is, xi = x0 + ih, i = 0, 1, . . . , n.
We now define the following operators.
28 CHAPTER 4. INTERPOLATION

4.3.1 Shift Operator E

The shift operator, denoted as E, shifts a function f (xi ) by h units along
the x-axis:

Ef (xi ) = f (xi + h)
For any integer n, the shift operator can shift by nh units:

E n f (xi ) = f (xi + nh)

4.3.2 Forward Difference Operator (∆)

The forward difference operator with step size h, denoted as ∆, calculates
the difference between the function value at xi + h and xi :

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

For small h, the forward difference operator approximates the derivative:

∆f (xi ) ≈ hf ′ (xi )
The n-th order forward difference can be given as:
n  
k n
X
n
∆ f (xi ) = (−1) f (xi + (n − k)h).
k=0
k

4.3.3 Backward Difference Operator (∇)

The backward difference operator with step size h, denoted as ∇, calcu-
lates the difference between the function value at xi and xi − h:

∇f (xi ) = f (xi ) − f (xi − h)

For small h, the backward difference operator approximates the derivative:

∇f (xi ) ≈ hf ′ (xi )
The n-th order backward difference can be given as:
n  
k n
X
n
∇ f (xi ) = (−1) f (xi − kh).
k=0
k
4.3. FINITE DIFFERENCE OPERATORS 29

4.3.4 Central Difference Operator (δ)

The central difference operator with step size h, denoted as δ, is defined
as the difference between the function values symmetrically around xi :

δf (xi ) = f (xi + h/2) − f (xi − h/2)

This operator gives a more accurate approximation of the derivative com-
pared to the forward and backward difference operators, reducing the trun-
cation error:

δf (xi ) ≈ f ′ (xi ) + O(h2 )

The n-th order central difference can be given as:
n  
k n
X
n
δ f (xi ) = (−1) f (xi + (n/2 − k)h).
k=0
k

4.3.5 Average Difference Operator (µ)

The average difference operator, denoted as µ, is defined as follows:

f (xi + h/2) + f (xi − h/2)

µf (xi ) =
2
This operator provides a balanced approximation by smoothing out varia-
tions.

Example 4.1 Show that

(i) ∆ = E − 1 and ∇ = 1 − E −1

(ii) ∆ + 1 = (1 − ∇)−1

(iii) δ = E 1/2 − E −1/2 and µ = 21 (E 1/2 + E −1/2 )

(iv) ∆ = ∇(1 − ∇)−1

q
2
(v) µ = 1 + δ4

(vi) ∆n = (E − 1)n and ∇n = (1 − E −1 )n

30 CHAPTER 4. INTERPOLATION
Chapter 5

Numerical Differentiation,
Integration, and ODE

5.1 Differentiation
Numerical differentiation is a technique used to approximate the derivative
of a function when only discrete data points are available. Three common
methods for first-order numerical differentiation are the forward difference,
backward difference, and central difference methods.

5.1.1 Forward Difference Derivative

The forward difference method approximates the derivative at a point xi by
using the function values at xi and xi+1 . Mathematically, it is expressed as:

f (xi+1 ) − f (xi )
f ′ (xi ) ≈
h
where h is the step size between xi and xi+1 .
The error in the forward difference method is of order O(h), meaning the
error decreases linearly with decreasing step size. Thus, it is a first-order
accurate approximation.

5.1.2 Backward Difference Derivative

The backward difference method approximates the derivative at a point xi
by using the function values at xi and xi−1 . The formula is given by:

f (xi ) − f (xi−1 )
f ′ (xi ) ≈
h
31
32 CHAPTER 5. DIFFERENTIATION, INTEGRATION AND ODE

Here, h is the step size between xi and xi−1 .

Similar to the forward difference, the backward difference has an order
O(h) error, making it a first-order accurate method.

5.1.3 Central Difference Derivative

The central difference method uses both the forward and backward points
around xi to approximate the derivative. The formula gives it:

f (xi+1 ) − f (xi−1 )
f ′ (xi ) ≈
2h
This method computes the slope of the secant line between the points (xi−1 , f (xi−1 ))
and (xi+1 , f (xi+1 )), centered at xi .
The central difference method has an error of order O(h2 ), meaning it
is a second-order accurate approximation. The error decreases quadratically
as the step size decreases, making it more accurate than the forward and
backward difference methods.
For higher accuracy in numerical differentiation, we can use three points
instead of just two. This results in higher-order approximations for both
the forward and backward difference methods. The three-point forward and
backward difference formulas give better estimates of the derivative at a point
xi .

5.1.4 Three-Point Forward Difference Derivative

The three-point forward difference derivative uses the values of the function
at xi , xi+1 , and xi+2 . The formula is:

−f (xi+2 ) + 4f (xi+1 ) − 3f (xi )

f ′ (xi ) ≈
2h
where h is the step size between consecutive points. The error in this ap-
proximation is of the order O(h2 ), which means that the error decreases
quadratically as the step size h decreases.

5.1.5 Three-Point Backward Difference Derivative

Similarly, the three-point backward difference derivative uses the values of
the function at xi , xi−1 , and xi−2 . The formula is:

3f (xi ) − 4f (xi−1 ) + f (xi−2 )

f ′ (xi ) ≈
2h
5.1. DIFFERENTIATION 33

Here, h is the step size between consecutive points. This method also provides
a higher-order approximation with an order O(h2 ) error, making it more
accurate than the two-point backward difference method.

5.1.6 General three-point formula for derivative

Lagrange interpolation provides a method to approximate a function by a
polynomial passing through given points. Using three points, we can con-
struct a second-degree polynomial and differentiate it to approximate the
derivative at one of the points. In this section, we derive the general three-
point formula for numerical differentiation using Lagrange interpolation.
Given three points (x1 , y1 ), (x2 , y2 ), and (x3 , y3 ), the Lagrange interpola-
tion polynomial P (x) is:

(x − x2 )(x − x3 )y1 (x − x1 )(x − x3 )y2 (x − x1 )(x − x2 )y3

P (x) = + +
(x1 − x2 )(x1 − x3 ) (x2 − x1 )(x2 − x3 ) (x3 − x1 )(x3 − x2 )

This polynomial passes through the three points x1 , x2 , and x3 . To

approximate the derivative at any of these points, we differentiate P (x).

Derivative at x1 :
(2x1 − x2 − x3 )y1 (x1 − x3 )y2 (x1 − x2 )y3
P ′ (x1 ) ≈ + +
(x1 − x2 )(x1 − x3 ) (x2 − x1 )(x2 − x3 ) (x3 − x1 )(x3 − x2 )

Derivative at x2 :
(x2 − x3 )y1 (2x2 − x1 − x3 )y2 (x2 − x1 )y3
P ′ (x2 ) ≈ + +
(x1 − x2 )(x1 − x3 ) (x2 − x1 )(x2 − x3 ) (x3 − x1 )(x3 − x2 )

Derivative at x3 :
(x3 − x2 )y1 (x3 − x1 )y2 (2x3 − x1 − x2 )y3
P ′ (x3 ) ≈ + +
(x1 − x2 )(x1 − x3 ) (x2 − x1 )(x2 − x3 ) (x3 − x1 )(x3 − x2 )

5.1.7 Formula for Second-Order Derivative

The second-order derivative at xi can be approximated using the Taylor
expansion around xi . We have the following formula:

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

f ′′ (xi ) ≈
h2
34 CHAPTER 5. DIFFERENTIATION, INTEGRATION AND ODE

5.2 Integration
In this section, we will discuss the Trapezoid Rule and Simpson’s Rule for
numerical integration.

5.2.1 Trapezoid Rule (Composite form)

The Trapezoid Rule is a numerical method used to approximate the definite
integral of a function. It is beneficial when the exact integral cannot be found
easily. The general formula for the Trapezoid Rule is given as:
Z b " n−1
#
h X
f (x) dx ≈ f (x0 ) + 2 f (xi ) + f (xn ) ,
a 2 i=1

where h = b−a n
is the width of each sub-interval, and n is the number of
sub-intervals. The points that divide the interval [a, b] are:

x0 = a, x1 = a + h, x2 = a + 2h, ..., xn = b.
In particular, if we take n = 1, then x0 = a and xn = x1 = b, and hence the
Basic Trapezoid Rule will be
Z b
h
f (x) ≈ [f (a) + f (b)].
a 2
Note: If n is not given in the question, then one should apply the basic
Trapezoid rule.

Geometrical Interpretation
The idea behind the Trapezoid Rule is to approximate the region under the
curve as a series of trapezoids. Each trapezoid’s area is determined by the
function values at the endpoints of the sub-interval. Summing the areas of
these trapezoids gives an approximation of the integral.

Error in the Trapezoid Rule

The error in the Trapezoid Rule is given by:

(b − a)3 ′′
ET = − f (ξ),
12n2
where ξ is some point in the interval [a, b]. This shows that the error decreases
as n increases.
5.2. INTEGRATION 35

Example 1: f (x) = x2
Let’s apply the Trapezoid Rule to approximate the integral of f (x) = x2 over
the interval [0, 2]. We will divide the interval into n = 4 sub-intervals.
b−a
a = 0, b = 2, h= = 0.5
4
The points dividing the interval are:

x0 = 0, x1 = 0.5, x2 = 1, x3 = 1.5, x4 = 2.
Evaluating the function at these points:

f (0) = 02 = 0, f (0.5) = (0.5)2 = 0.25, f (1) = 12 = 1,

f (1.5) = (1.5)2 = 2.25, f (2) = 22 = 4.

Now applying the Trapezoid Rule formula:

Z 2
0.5
x2 dx ≈ [0 + 2(0.25 + 1 + 2.25) + 4]
0 2
=0.25 × [0 + 2(3.5) + 4]
=0.25 × 11
=2.75.

The exact value of the integral is:

Z 2
x3 2 8
x2 dx = = ≈ 2.6667.
0 3 0 3
Thus, the Trapezoid Rule gives an approximation of 2.75, which is close to
the exact value.
2
Example 2: f (x) = e−x
2
Now, let’s approximate the integral of f (x) = e−x over the interval [0, 1]
using the Trapezoid Rule with n = 4 sub-intervals.
b−a
a = 0, b = 1, h= = 0.25
4
The points dividing the interval are:

x0 = 0, x1 = 0.25, x2 = 0.5, x3 = 0.75, x4 = 1.

36 CHAPTER 5. DIFFERENTIATION, INTEGRATION AND ODE

Evaluating the function at these points:

2 2 2
f (0) = e−(0) = 1, f (0.25) = e−(0.25) ≈ 0.9394, f (0.5) = e−(0.5) ≈ 0.7788,
2 2
f (0.75) = e−(0.75) ≈ 0.5698, f (1) = e−(1) ≈ 0.3679.
Now applying the Trapezoid Rule formula:

Z 1
2 0.25
e−x dx ≈ [1 + 2(0.9394 + 0.7788 + 0.5698) + 0.3679]
0 2
=0.125 × [1 + 2(2.2879) + 0.3679]
=0.125 × [1 + 4.5758 + 0.3679]
=0.125 × 5.9437
≈0.74296.
The exact value of the integral is approximately 0.7468, so the Trapezoid
Rule gives a good approximation of 0.74296.

5.2.2 Simpson’s 1/3rd Rule

Simpson’s Rule or Simpson’s 1/3rd Rule is another numerical method for
approximating definite integrals. It generally provides more accurate results
than the Trapezoid Rule, especially for smooth functions. It approximates
the region under the curve using parabolic segments instead of straight lines
or trapezoids.
The formula for Simpson’s Rule is:

" n−1 n−2

#
Z b
h X X
f (x) dx ≈ f (x0 ) + 4 f (xi ) + 2 f (xi ) + f (xn ) ,
a 3 i=1,3,5,... i=2,4,6,...

where h = b−an
is the width of each sub-interval, and n is the number of sub-
intervals. It is important to note that for Simpson’s Rule n must be
even. The points that divide the interval [a, b] are:
x0 = a, x1 = a + h, x2 = a + 2h, ..., xn = b.
The Basic Simpson’s Rule can be obtained by taking n = 2 in the above
formula. For n = 2, x0 = a, x1 = a + h and xn = x2 = b, and hence the
formula will be
Z b
h
f (x) ≈ [f (x0 ) + 4f (x1 ) + f (x2 )].
a 3
5.2. INTEGRATION 37

Note: If n is not given in the question, then one should apply the Basic
Simpson’s Rule.

Geometrical Interpretation
Simpson’s Rule fits a quadratic polynomial (parabola) to every pair of adja-
cent intervals and approximates the area under the curve using this parabola.
This makes Simpson’s Rule more accurate for smooth and continuous func-
tions.

Error in Simpson’s Rule

The error in Simpson’s Rule is related to the fourth derivative of the function
f (x). For a smooth function, the error ES is approximately:

(b − a)5 (4)
ES = − f (ξ),
180n4
where ξ is some point in the interval [a, b]. This shows that Simpson’s Rule
is generally more accurate than Trapezoid.

Example 1: f (x) = x2
Let’s apply Simpson’s Rule to approximate the integral of f (x) = x2 over
the interval [0, 2], using n = 4 sub-intervals.
b−a
a = 0, b = 2, h= = 0.5
4
The points dividing the interval are:

x0 = 0, x1 = 0.5, x2 = 1, x3 = 1.5, x4 = 2.
Evaluating the function at these points:

f (0) = 02 = 0, f (0.5) = (0.5)2 = 0.25, f (1) = 12 = 1,

f (1.5) = (1.5)2 = 2.25, f (2) = 22 = 4.
Now applying Simpson’s Rule formula:

Z 2
0.5 0.5
x2 dx ≈ [0 + 4(0.25 + 2.25) + 2(1) + 4] = × [0 + 4(2.5) + 2 + 4]
0 3 3
0.5 0.5 8
= × [10 + 2 + 4] = × 16 = ≈ 2.6667.
3 3 3
38 CHAPTER 5. DIFFERENTIATION, INTEGRATION AND ODE

The exact value of the integral is:

Z 2
2 x3 2 8
x dx = = ≈ 2.6667.
0 3 0 3
Thus, Simpson’s Rule gives the exact value in this case.
2
Example 2: f (x) = e−x
2
Now, let’s approximate the integral of f (x) = e−x over the interval [0, 1]
using Simpson’s Rule with n = 4 sub-intervals.
b−a
a = 0, b = 1, h= = 0.25
4
The points dividing the interval are:

x0 = 0, x1 = 0.25, x2 = 0.5, x3 = 0.75, x4 = 1.

Evaluating the function at these points:

2 2 2
f (0) = e−(0) = 1, f (0.25) = e−(0.25) ≈ 0.9394, f (0.5) = e−(0.5) ≈ 0.7788,
2 2
f (0.75) = e−(0.75) ≈ 0.5698, f (1) = e−(1) ≈ 0.3679.
Now applying Simpson’s Rule formula:
Z 1
2 0.25
e−x dx ≈ [1 + 4(0.9394 + 0.5698) + 2(0.7788) + 0.3679]
0 3
0.25
= × [1 + 4(1.5092) + 2(0.7788) + 0.3679]
3
0.25 0.25
= × [1 + 6.0368 + 1.5576 + 0.3679] = × 8.9623 ≈ 0.7470.
3 3
The exact value of the integral is approximately 0.7468, so Simpson’s Rule
gives an excellent approximation of 0.7470.

5.3 Ordinary Differential Equation

We discuss three methods for solving ordinary differential equations (ODEs):
Euler’s method, the second-order Runge-Kutta method (RK2), and the clas-
sical fourth-order Runge-Kutta method (RK4). These methods approximate
the solution of the initial value problem:
dy
= f (t, y), y(t0 ) = y0 . (5.1)
dt
5.3. ORDINARY DIFFERENTIAL EQUATION 39

5.3.1 Euler’s Method

Euler’s method is the simplest numerical method, using a linear approxima-
tion for yn+1 . The formula for the next approximation is:

yn+1 = yn + h · f (tn , yn )

Example

Consider the ODE:

dy
= t + y, y(0) = 1
dt
With step size h = 0.1, applying Euler’s method for one step gives:

y1 = y0 + 0.1 · f (0, 1) = 1 + 0.1 · (0 + 1) = 1.1

Similarly, you can calculate y2 , y3 , ....

5.3.2 Second-Order Runge-Kutta Methods (RK2)

The general second-order Runge-Kutta method is a family of numerical tech-
niques for solving first-order ordinary differential equations 5.1.
In a general second-order Runge-Kutta method, the next approximation
yn+1 at tn+1 is given by:

yn+1 = yn + h [a1 k1 + a2 k2 ] ,

where:

• h is the step size.

• k1 and k2 are the slope approximations, given by:

k1 = f (tn , yn ), k2 = f (tn + p1 h, yn + u11 hk1 ).

• a2 , p1 , and u11 are coefficients that define a specific second-order Runge-

Kutta method satisfying the conditions a1 + a2 = 1, a2 p1 = 1/2, and
a2 u11 = 1/2.
40 CHAPTER 5. DIFFERENTIATION, INTEGRATION AND ODE

5.3.3 Special Cases of 2nd-Order Runge-Kutta Meth-

ods
1. Modified Euler’s Method (Mid-point Method)

The Modified Euler’s Method (Mid-Point method), also called the second-
order Runge-Kutta method, is another variant, defined by:

a1 = 0, a2 = 1, p1 = 1/2, u11 = 1/2.

The formula for the Mid-Point method is:

k1 = f (tn , yn ),
 
h h
k2 = f tn + , yn + k1 ,
2 2
yn+1 = yn + hk2 .

2. Heun’s Method

Heun’s method is a common variant of the second-order Runge-Kutta method,

characterized by:
a1 = 1/2, a2 = 1/2,

and hence
p1 = 1, u11 = 1.

The formula for Heun’s method is:

k1 = f (tn , yn ),
k2 = f (tn + h, yn + hk1 ),
h
yn+1 = yn + (k1 + k2 ).
2

3. Ralston’s Method

Ralston’s method aims to minimize the error coefficient and is characterized

by:
1 3 2 2
a1 = , a2 = , p1 = , u11 = .
4 4 3 3
The formula for Ralston’s method is:
5.4. CLASSICAL FOURTH-ORDER RUNGE-KUTTA METHOD (RK4)41

k1 = f (tn , yn ),
 
2h 2h
k2 = f tn + , yn + k1 ,
3 3
h
yn+1 = yn + (k1 + 3k2 ).
4

5.4 Classical fourth-Order Runge-Kutta Method

(RK4)
The RK4 method uses four evaluations of the function:
 
h h
k1 = f (xn , yn ), k2 = f xn + , yn + k1
2 2
 
h h
k3 = f xn + , yn + k2 , k4 = f (xn + h, yn + h · k3 )
2 2
The solution is updated using:
h
yn+1 = yn + (k1 + 2k2 + 2k3 + k4 )
6

Example
dy
For the same ODE, dx
= x + y, y(0) = 1, and step size h = 0.1:

k1 = f (0, 1) = 0 + 1 = 1

k2 = f (0.05, 1 + 0.05 · 1) = 0.05 + 1.05 = 1.1

k3 = f (0.05, 1 + 0.05 · 1.1) = 0.05 + 1.055 = 1.105
k4 = f (0.1, 1 + 0.1 · 1.105) = 0.1 + 1.1105 = 1.2105
0.1
y1 = 1 + (1 + 2 · 1.1 + 2 · 1.105 + 1.2105) = 1.11035
6
Similarly, one can calculate y2 , y3 , ..... etc.
42 CHAPTER 5. DIFFERENTIATION, INTEGRATION AND ODE

Miscellaneous Exercise
1. A real root of the equation x3 − 5x + 1 = 0 lies in the interval (0, 1).
Perform three iterations of the Regula Falsi Method to obtain the root.
2. Perform three iterations of the Bisection Method to find the root of the
equation Cos(x) − xex in (0, 1).
3. Discuss the order of convergence of the Secant method and give the
geometrical interpretation of the method.
√
4. (a) Verify x = a is a fixed point of the function h(x) = 12 x + xa2 .

Determine the√ order of convergence of the sequence pn = h(pn−1 )

towards x = a.
(b) Use the Secant method to find the root of 3x + sin(x) − ex = 0 in
the interval ]1, 10[. Perform three iterations.
5. Prove that Newton’s Method is of order two using x3 + 2x2 − 3x − 1 = 0
and the initial approximation x0 = 2.
6. Define a lower and an upper triangular matrix. Solve the system of
equations:
3x1 + 2x2 − x3 = −1
6x1 + 8x2 + x3 = 14
x1 + 2x2 + 7x3 = 1
by Gauss-Jacobi method. Take the initial approximation as X (0) =
(0, 0, 0) and do three iterations.
7. Set up the Gauss-Seidel iteration scheme to solve the system of equa-
tions:
4x1 + 2x2 − x3 = 1
2x1 + 4x2 + x3 = −1
−x1 + x2 + 4x3 = 1
Take the initial approximation as X (0) = (0, 0, 0) and do three itera-
tions.
8. (a) Construct the Lagrange form of the interpolating polynomial from
the following data:
x f (x) = ln(x)
1 ln(1)
2 ln(2)
3 ln(3)
5.4. CLASSICAL FOURTH-ORDER RUNGE-KUTTA METHOD (RK4)43

(b) Prove that for n + 1 distinct nodal points x0 , x1 , . . . , xn , there

exists a unique interpolating polynomial of at most degree n.
(c) Find the maximum value of the step size h that can be used in
the tabulation of f (x) = ex on the interval [0, 1] so that the error
in the linear interpolation of f (x) is less than 5 × 10−4 .

9. (a) Define the backward difference operator ∇ and the Newton di-
vided difference. Prove that:
∇n f
f [x0 , x1 , . . . , xn ] =
n!hn
where h = xi+1 − xi .
(b) Construct the divided difference table for the following data set
and then write out the Newton form of the interpolating polyno-
mial:
x y
−7 10
−5 5
−4 2
−1 10
Find the approximation of y for x = −3.
(c) Use the formula

f (x0 + h) − f (x0 )
f ′ (x0 ) =
h

to approximate the derivative of f (x) = 1+x+x3 at x0 = 1 taking

h = 1, 0.1, 0.001. What is the order of approximation?
R1
10. (a) Approximate the value of 0 e−x dx using the Trapezoidal rule and
verify that the theoretical error bound holds for the same.
Rb
(b) State Simpson’s 1/3rd rule for the evaluation of a f (x) dx and
prove that it has a degree of precision 3.
(c) Use Euler’s method to approximate the solution of the initial value
problem:
dx 1 + x2
= , x(1) = 0, 1 ≤ t ≤ 4
dt t
taking 5 steps.
44 CHAPTER 5. DIFFERENTIATION, INTEGRATION AND ODE

11. Define the rate of convergence of an iterative scheme {xn }. Use the
bisection method to determine the root of the equation:

x5 + 2x − 1 = 0 on (0, 1).

Further, compute the theoretical error bound at the end of the fifth
iteration and the next enclosing (bracketing) interval.

12. Differentiate between the method of false position and the secant method.
Apply the method of false position to solve:

cos(x) − x = 0

to determine an approximation to the root lying in the interval (0, 1)

until the absolute error is less than 10−5 (p = 0.739085).

13. (a) Using scaled partial pivoting during the factoring step, find ma-
trices L, U , and P such that LU = P A, where:
 
4 2 −1
A= 1 2 3 
−1 1 4

(b) Let g be a continuous function on the closed interval [a, b]. Sup-
pose g is differentiable on the open interval (a, b) and there exists a
positive constant k < 1 such that |g ′ (x)| ≤ k < 1 for all x ∈ (a, b).
Then:
i. The sequence {pn } generated by pn = g(pn−1 ) converges to
the fixed point p for any p0 ∈ [a, b].
ii. Find the approximated root of f (x) = x3 + 2x2 − 3x − 1 using
the secant method, taking p0 = 2 and p1 = 1, until |pn − pn−1 |
is sufficiently small.
(c) Determine the spectral radius ρ of the matrix:
 
−1 0 −2
A =  3 2 −1
1 3 0


14
Hence, solve the system Ax = b, given b = −3.
7
5.4. CLASSICAL FOURTH-ORDER RUNGE-KUTTA METHOD (RK4)45

14. Use the Jacobi method to solve the following system of linear equations:

2x1 − x2 = 3
−x1 + 4x2 + x3 = 7
x2 + 5x3 = 1
 
0
(0)
Use the initial approximation x = 0 and perform three iterations.

0

15. (a) If x0 , x1 , x2 , . . . , xn are n + 1 distinct points and f is defined at

x0 , x1 , . . . , xn , prove that the interpolating polynomial Pn of de-
gree at most n is unique.
(b) Define the shift operator E and the central difference operator ∆.
Prove that:
f (x0 + h) − f (x0 )
f ′ (x0 ) =
h
approximates the derivative of f (x) for small h, and verify that
the order of approximation is O(h).
(c) Derive the second-order forward difference approximation to the
first-order derivative of a function.

16. Approximate the value of the integral:

Z 1
e−x dx
0

using Simpson’s rule. Further, verify the theoretical error bound.

17. Apply Euler’s method to approximate the solution of the initial value
problem:
x′ = 1 + x2 , x(1) = 0, 1 ≤ t ≤ 4
taking 5 steps.

18. For the function f (x) = ex , construct the Lagrange form of the inter-
polating polynomial passing through the points (−1, e−1 ), (0, 1), and
(1, e). Estimate f (x) using the polynomial and determine the error in
the approximation. Verify that the theoretical error bound is satisfied.
46 CHAPTER 5. DIFFERENTIATION, INTEGRATION AND ODE