NUMERICAL METHODS KIT

FOR MATLAB, SCILAB AND OCTAVE USERS

Rohan Verma
MSc (Physics), Indian Institute of Technology Kanpur, India
BSc (Physics), University of Delhi, India
© Rohan Verma 2024
All rights reserved
 Preface
I take great pleasure in presenting to the readers this book entitled
“Numerical Methods Kit”. The book has been designed for Science,
Engineering, Mathematics and Statistics students.
A look at the contents of the book will give the reader a clear idea of the
variety of numerical methods discussed and analysed.
The book has been written in a concise and lucid style with proper
explanation of Mathematics involved in each method.
Each method is explained with solved examples, computer programs and
their results like a screenshot of the graphic window and console window.
The careful organisation of figures, solved examples, codes, graphic window
and console window help the students grasp quickly.
Despite careful editing, there is zero probability that this book is error-free. If
anything looks amiss, please send that part to my email id given below.

Email: rohanv.i341@[Link] Rohan Verma

 CONTENTS
𝟏. Solutions of Algebraic and Transcendental Equations 1
1.1 Bisection Method 1
1.2 Newton-Raphson Method 4
1.3 Regula-Falsi Method 8
1.4 Secant Method 14

2. Interpolation 19
2.1 Finite Difference Operator 19
2.2 Newton’s Gregory Forward Interpolation Method 20
2.3 Newton’s Gregory Backward Interpolation Method 23
2.4 Lagrange’s Interpolation Method 26

3. Curve Fitting 28
3.1 Line Fitting 28
3.2 Parabola Fitting 30

4. Numerical Differentiation 33
4.1 Equal Interval 33
4.2 Unequal Interval 34

5. Numerical Integration 35
5.1 Newton-Cotes Formula 36
5.2 Trapezoidal Rule 37
5.3 Simpson's 1/3 rule 39
5.4 Simpson's 3/8 rule 40
5.5 Monte Carlo method 42

6. Ordinary Differential Equation 45

6.1 Classification of differential equations 46
6.2 Euler Method 47
6.3 Runge–Kutta method 61
7. Linear Algebra 63
7.1 Gaussian Elimination Method 67
7.2 Gauss-Jordan Elimination Method 73
7.3 Jacobi Method 77
7.4 Gauss-Seidel Method 81

8. Finite Difference Method 85

8.1 Classification of PDE 85
8.2 Finite Difference Fundamentals 87
8.3 2D Laplace's Equation 92
8.4 1D Heat Equation 96
8.5 1D Wave Equation 100
CHAPTER
Solutions of Algebraic and Transcendental
Equations
𝑓(𝑥) = 0 is called equation and the values of 𝑥 for which this equation
satisfied is called roots of the equation.
Algebraic Equation - If 𝑓(𝑥) is purely polynomial. E.g. 𝑥 − 𝑥 + 1 = 0
Transcendental Equations - If 𝑓(𝑥) containing logarithm, trigonometric,
exponential function, etc. E.g. 𝑥𝑒 − 𝑐𝑜𝑠 𝑥 = 0

𝟏. 𝟏 Bisection Method
The bisection method applies to a continuous function. Let our equation is
𝑓(𝑥) = 0. Here (Figure 1.1) we are going to consider two values of 𝑥, 𝑎 and 𝑏
such that 𝑓(𝑎) < 0 and 𝑓(𝑏) > 0. Clearly, we can see in Figure 1.1 that the root
must lie somewhere between 𝑎 and 𝑏.

Figure 𝟏. 𝟏

Take 𝑐 =

There could be three possibilities.

(a) 𝑓(𝑐) > 0

replace 𝑏 with 𝑐

(b) 𝑓(𝑐) < 0

1
replace 𝑎 with 𝑐

(c) 𝑓(𝑐) = 0

𝑐 is the root of equation 𝑓(𝑥) = 0. We are replacing 𝑎 or 𝑏 with 𝑐 to make a

smaller interval for finding the root. In this way, we keep bisecting the
interval for finding the root. Eventually, the interval will coincide with the 𝑐
for which 𝑓(𝑐) = 0.

Working rule

(a) Find 𝑎 and 𝑏 such that 𝑓(𝑎) < 0 and 𝑓(𝑏) > 0 or {𝑓(𝑎)𝑓(𝑏)} < 0.
(b) Find the first approximate root using.
𝑎+𝑏
𝑥 =
2
Now calculate 𝑓(𝑥 ) and Examine its sign. If 𝑓(𝑥 ) < 0 it implies that root lies
between 𝑥 and 𝑏. 2nd approximate root is given by 𝑥 = .
If 𝑓(𝑥 ) > 0 it implies that root lies between 𝑎 and 𝑥 .Then 2nd approximate
root is given by 𝑥 = . Calculate 𝑓(𝑥 ). Repeat the (b) step until the
required accuracy of root or 𝑓(𝑥 ) → 0.
Example 𝟏. 𝟏 Find the 5th approximate real root of the equation 𝑥 −1=0
using the Bisection Method.
Let 𝑓(𝑥) = 𝑥 − 1 = 0
take 𝑎 = 0 .8 and 𝑏 = 1.1
𝑓(0.8) = −0.36 < 0
𝑓(1.1) = 0.21 > 0
First approximate root is
0.8 + 1.1
𝑥 = = 0.95
2
𝑓(0.95) = −0.0975
New value of 𝑎 = 0.95
Second approximate root is
0.95 + 1.1
𝑥 = = 1.025
2
𝑓(1.025) = 0.050625
New value of 𝑏 = 1.025
Third approximate root is
0.95 + 1.025
𝑥 = = 0 ⋅ 9875
2
𝑓(0.9875) = −0.0248437

2
New value of 𝑎 = 0.9875
Fourth approximate root is
0.9875 + 1.025
𝑥 = = 1.00625
2
𝑓(1.00625) = 0.0125391
New value of 𝑏 = 1.00625
Fifth approximate root is
0.9875 − 1 ⋅ 00625
𝑥 = = 0.996875
2
𝑓(0.996875) = −0.0062402
In the previous example, we can see that 𝑥 → 𝑟, where 𝑟 is the root of
equation 𝑓(𝑥) = 0.
Scilab Code 𝟏. 𝟏
//bisection method
clc
clear
function z=f(x)
z=x^3
endfunction
n=input("input the interval where you want to find the roots of given
equation")
a=n(1)
b=n(2)
c=(a+b)/2
tol=0.00000000001
if f(a)*f(b)<0 then

disp("convergence")
while abs(f(c))>tol
if f(a)*f(c)<0 then
b=c
else
a=c
end

disp(c)
c=(a+b)/2
end
disp(""+string(c)+" is the root of eqution")
else
disp('root does not exist in this interval')
end

3
Console Window 𝟏. 𝟏

𝟏. 𝟐 Newton-Raphson Method
The Newton - Raphson Method applies to continuous and differentiable
function. Let our equation is 𝑓(𝑥) = 0. In Newton - Raphson Method First we
need to consider some initial guess let take 𝑥 , and then we draw a tangent at
𝑥 , 𝑓(𝑥 ) (Figure 1.2). This tangent intersects at x-axis at 𝑥 . from Figure 1.2
we can see that 𝑥 is closer to root. To get closer values (𝑥 , 𝑥 … 𝑥 ), we repeat
the same process for 𝑥 and so on until 𝑥 → 𝑟, where 𝑟 is the root of equation
𝑓(𝑥) = 0. The rate of convergence of the Newton - Raphson Method is faster
than the Bisection method.

4
 Figure 𝟏. 𝟐

Working rule

(a) Consider equation 𝑓(𝑥) = 0 to find root by Newton - Raphson method

And take initial guess root 𝑥 .
(b) Find 1st approximate root 𝑥 by constructing a tangent equation at
𝑥 , 𝑓(𝑥 ) and putting 𝑦 = 0.
𝑦 − 𝑓(𝑥 ) = 𝑓 (𝑥 )(𝑥 − 𝑥 )
Putting 𝑦 = 0
𝑓(𝑥 )
𝑥 =𝑥 −
𝑓 (𝑥 )
(c) 2 approximate root
nd

𝑓(𝑥 )
𝑥 =𝑥 −
𝑓 (𝑥 )
(d) repeat the (b) step until 𝑥 → 𝑟, where 𝑟 is the root of equation 𝑓(𝑥) = 0.

General Formula
𝑓(𝑥 )
𝑥 =𝑥 − , 𝑛 = 1,2,3, … (1.1)
𝑓 (𝑥 )

Example 𝟏. 𝟐 Find the 6 approximate real root of the equation 𝑥 − 2𝑥 +

1 = 0 by Newton - Raphson Method.

Let 𝑓(𝑥) = 𝑥 − 2𝑥 + 1

5
𝑓’(𝑥) = 2𝑥 − 2
Consider initial guess 3
1st approximate root
𝑓(3)
𝑥 =3−
𝑓 (3)
𝑥 =2
𝑓(2) = 1
2nd approximate root
𝑓(2)
𝑥 =2−
𝑓 (2)
𝑥 = 1.5
𝑓(1.5) = 0.25
3rd approximate root
𝑓(1.5)
𝑥 = 1.5 −
𝑓 (1.5)
𝑥 = 1.25
𝑓(1.25) = 0.0625
4th approximate root
𝑓(1.25)
𝑥 = 1.25 −
𝑓 (1.25)
𝑥 = 1.125
𝑓(1.125) = 0.015625
5th approximate root
𝑓(1.125)
𝑥 = 1.125 −
𝑓 (1.125)
𝑥 = 1.0625
𝑓(1.0625) = 0.0039062
6th approximate root
𝑓(1.0625)
𝑥 = 1.0625 −
𝑓 (1.0625)

6
 𝑥 = 1.03125
𝑓(1.03125) = 0.0009766
6th approximate root is 1.03125 for equation 𝑥 − 2𝑥 + 1 = 0.

Scilab Code 𝟏. 𝟐

//Newton-Raphson Method
clc
clear
x0=input("guess the root of equation:-")
function z=f(x)
z=x^2-2*x+1
endfunction
disp("convergence")
tol=0.0000000001
while abs(f(x0))>tol
x1=x0-f(x0)/numderivative(f,x0)
disp(x1)
x0=x1
end
disp("root of equation is "+string(x0)+"")

7
Console Window 𝟏. 𝟐

𝟏. 𝟑 Regula-Falsi Method
Let our equation is 𝑓(𝑥) = 0. First, we need to consider a root-finding interval
(𝑎, 𝑏). Such that 𝑓(𝑎) < 0 and 𝑓(𝑏) > 0 or {𝑓(𝑎)𝑓(𝑏)} < 0. It is quite obvious
that root must lie somewhere between 𝑎 and 𝑏 (see Figure 1.3). Here we can
see that it is looking very similar to the Bisection Method till now. But in the
Regula Falsi Method instead of the bisecting interval, we draw a line segment
joining (𝑎, 𝑓(𝑎)) and (𝑏, 𝑓(𝑏)) see Figure 1.3, this line segment intersects at x-

8
axis at 𝑥 .This 𝑥 is closer to the root of equation 𝑓(𝑥) = 0, we can say that 𝑥
is our first approximate root of the equation.

Figure 1.3

Now for 2nd approximate root, we check the sign of 𝑓(𝑥 )

If 𝑓(𝑥 ) > 0 replace 𝑏 with 𝑥 .

If 𝑓(𝑥 ) < 0 replace 𝑎 with 𝑥 .

Keep doing this process until 𝑥 → 𝑟, where 𝑟 is the root of equation 𝑓(𝑥) = 0.

Working rule

(a) Let 𝑓(𝑥) = 0 be the equation to solve. Take 𝑎 and 𝑏 such that 𝑓(𝑎) < 0
and 𝑓(𝑏) > 0.

(b) Find 1st approximate root 𝑥 by constructing equation of line segment

Joining (𝑎, 𝑓(𝑎)) and (𝑏, 𝑓(𝑏)) and putting 𝑦 = 0.
𝑓(𝑏) − 𝑓(𝑎)
𝑦 − 𝑓(𝑎) = (𝑥 − 𝑎)
𝑏−𝑎
Putting 𝑦 = 0

9
 𝑎𝑓(𝑏) − 𝑏𝑓(𝑎)
𝑥 = (1.2)
𝑓(𝑏) − 𝑓(𝑎)
(c) Check the sign of 𝑓(𝑥 ).

If 𝑓(𝑥 ) > 0 replace 𝑏 with 𝑥

2nd approximate root

𝑎𝑓(𝑥 ) − 𝑥 𝑓(𝑎)
𝑥 =
𝑓(𝑥 ) − 𝑓(𝑎)
If 𝑓(𝑥 ) < 0 replace 𝑎 with 𝑥

2nd approximate root

𝑥 𝑓(𝑏) − 𝑏𝑓(𝑥 )
𝑥 =
𝑓(𝑏) − 𝑓(𝑥 )
(c) repeat the third step until convergence reached.

Example 1.3 Find the 7th approximate root of 𝑥𝑒 = 𝑐𝑜𝑠(𝑥) by Regula Falsi
method.

Let 𝑓(𝑥) = 𝑥𝑒 − 𝑐𝑜𝑠(𝑥) = 0 and 𝑎 = 0, 𝑏 = 1

𝑓(0) = −1 < 0 and 𝑓(1) = 2.1779795 > 0

1st approximate root

0𝑓(1) − 1𝑓(0)
𝑥 =
𝑓(1) − 𝑓(0)
𝑥 = 0.3146653
𝑓(0.3146653) = −0.5198712 < 0
New value of 𝑎 = 0.3146653

2nd approximate root

0.3146653𝑓(1) − 1𝑓(0.3146653)
𝑥 =
𝑓(1) − 𝑓(0.3146653)
𝑥 = 0.4467281
𝑓(0.4467281) = −0.2035448 < 0
New value of 𝑎 = 0.4467281

10
3rd approximate root
0.4467281 (1) − 1𝑓(0.4467281 )
𝑥 =
𝑓(1) − 𝑓(0.4467281 )
𝑥 = 0.4940153
𝑓(0.4940153) = −0.0708023 < 0
New value of 𝑎 = 0.4940153

4th approximate root

0.4940153 𝑓(1) − 1𝑓(0.4940153 )
𝑥 =
𝑓(1) − 𝑓(0.4940153 )
𝑥 = 0.5099461
𝑓(0.5099461) = −0.0236077 < 0
New value of 𝑎 = 0.5099461

5th approximate root

0.5099461 𝑓(1) − 1𝑓(0.5099461 )
𝑥 =
𝑓(1) − 𝑓(0.5099461 )
𝑥 = 0.515201
𝑓(0.515201) = −0.0077601 < 0
New value of 𝑎 = 0.515201

6th approximate root

0.515201 𝑓(1) − 1𝑓(0.515201 )
𝑥 =
𝑓(1) − 𝑓(0.515201)
𝑥 = 0.5169222
𝑓(0.5169222) = −0.0025389 < 0
New value of 𝑎 = 0.5169222

7th approximate root

0.5169222 𝑓(1) − 1𝑓(0.5169222 )
𝑥 =
𝑓(1) − 𝑓(0.5169222)
𝑥 = 0.5174847
𝑓(0.5174847) = −0.0008294

11
7th approximate root is 0.5174847 for equation 𝑥𝑒 − 𝑐𝑜𝑠(𝑥) = 0.

Scilab Code 𝟏. 𝟑

//Regula-Falsi Method
clc
clear
clf
x=[-1:0.1:1]
function z=f(x)
z=x*exp(x)-cos(x)
endfunction
plot(x,f,'c')
xgrid(3)
xlabel("x","Fontsize",4)
ylabel("y","Fontsize",4)
title("y=x*exp(x)-cos(x)","Fontsize",4)
a=0
b=1
for i=1:10
x1=(a*f(b)-b*f(a))/(f(b)-f(a))
disp('root after '+string(i)+' iteration is '+string(x1)+'')
if f(x1)< 0 then
a=x1
end
if f(x1)> 0 then
b=x1
end
if f(x1)==0 then
break
end
end
disp(''+string(x1)+' is root of equation' )
disp('value of f('+string(x1)+') is = '+string(f(x1))+'')

12
Graphic Window 𝟏. 𝟏

13
Console Window 𝟏. 𝟑

𝟏. 𝟒 Secant Method
Secant method is an improved form of Regula Falsi method. In this method,
we don’t need to check the sign of the approximate root. Like we have done in
the Regula Falsi method. Let our equation is 𝑓(𝑥) = 0 with 𝑓(𝑥 ) < 0 and
𝑓(𝑥 ) > 0. So that mean root must lie between 𝑥 and 𝑥 . First approximate
root 𝑥 is given by
𝑥 𝑓(𝑥 ) − 𝑓(𝑥 )𝑥
𝑥 =
𝑓(𝑥 ) − 𝑓(𝑥 )
Similarly, Second approximate root 𝑥 is given without checking the sing of
𝑓(𝑥 ).
𝑥 𝑓(𝑥 ) − 𝑓(𝑥 )𝑥
𝑥 =
𝑓(𝑥 ) − 𝑓(𝑥 )
General Formula
𝑥 𝑓(𝑥 ) − 𝑥 𝑓(𝑥 )
𝑥 = ,𝑛 ≥ 1 (1.3)
𝑓(𝑥 ) − 𝑓(𝑥 )

14
Keep finding 𝑛th approximate root until 𝑥 → 𝑟 where 𝑟 is the root of equation
𝑓(𝑥) = 0.

Example 1.4 Find the 7th approximate real root of the equation 2 −𝑥 =0
in interval −1 ≤ 𝑥 ≤ 1 by secant method.

Let 𝑓(𝑥) = 2 − 𝑥 = 0 and 𝑥 = −1, 𝑥 = 1

𝑓(−1) = −0.5 < 0 and 𝑓(1) = 1 > 0

1st approximate root

−1𝑓(1) − 1𝑓(−1)
𝑥 =
𝑓(1) − 𝑓(−1)
𝑥 = −0.3333333
𝑓(−0.3333333) = 0.6825894
2nd approximate root
1𝑓(−0.3333333) − (−0.3333333)𝑓(1)
𝑥 =
𝑓(−0.3333333) − 𝑓(1)
𝑥 = −3.2006581
𝑓(−3.2006581) = −10.135443
3rd approximate root
−0.3333333 𝑓(−3.2006581) − (−3.2006581)𝑓(−0.3333333 )
𝑥 =
𝑓(−3.2006581) − 𝑓(−0.3333333 )
𝑥 = −0.514254
𝑓(−0.514254) = 0.4356977
4th approximate root
−3.2006581 𝑓(−0.514254) − (−0.514254)𝑓(3.2006581 )
𝑥 =
𝑓(−0.5142541) − 𝑓(−3.2006581 )
𝑥 = −0.6249762
𝑓(−0.6249762) = 0.2578352
5th approximate root
−0.514254 𝑓(−0.6249762) − 𝑓(−0.6249762)(−0.514254 )
𝑥 =
𝑓(−0.6249762) − 𝑓(−0.514254)

15
 𝑥 = −0.7854827
𝑓(−0.7854827) = −0.0368252

6th approximate root

−0.6249762 𝑓(−0.7854827) − (−0.7854827)𝑓(−0.6249762 )
𝑥 =
𝑓(−0.7854827) − 𝑓(−0.6249762)
𝑥 = −0.7654234
𝑓(−0.7654234) = 0.0024078
7th approximate root
−0.7854827 𝑓(−0.7654234) − (−0.7654234)𝑓(−0.7854827 )
𝑥 =
𝑓(−0.7654234) − 𝑓(−0.7854827)
𝑥 = −0.7666544
𝑓(−0.7666544) = 0.0000199
7th approximate root is −0.7666544 for equation 2 − 𝑥 = 0.

Scilab Code 𝟏. 𝟒

//Secant Method
clc
clear
clf
x=[-1:0.1:1]
function z=f(x)
z=2^x-x^2
endfunction
plot(x,f,'c')
xgrid(3)
xlabel("x","Fontsize",4)
ylabel("y","Fontsize",4)
title("y=2^x-x^2","Fontsize",4)
a=-1
b=1
for i=1:10
x1=(a*f(b)-b*f(a))/(f(b)-f(a))
disp('root after '+string(i)+' iteration is '+string(x1)+'')
disp(f(x1))

16
 a=b
b=x1
if f(x1)==0 then
break
end
end
disp(''+string(x1)+' is root of equation' )
disp('value of f('+string(x1)+') is = '+string(f(x1))+'')

Graphic Window 𝟏. 𝟐

17
Console Window 𝟏. 𝟒

18
CHAPTER 2
Interpolation
Interpolation is the way of estimating function 𝑓(𝑥) for given data (𝑥 , 𝑦 )
within the range. Whereas Extrapolation is about estimating the function
outside the given data.

There are many methods for Interpolation and Extrapolation. A few of them
are listed below.

(a) Newton’s Gregory Forward Interpolation Method

(b) Newton’s Gregory Backward Interpolation Method

(c) Lagrange’s Interpolation Method

(d) Newton’s Divided Difference Interpolation method

Newton’s forward and backward interpolation method are used for constant
step size (|𝑥 − 𝑥 | = 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡).

Lagrange and Newton’s Divided Difference Interpolation method are used for
varying step size.

It’s good to use Newton’s forward method for better estimation 𝑓(𝑥 ) when 𝑥
is closer to the first value of range and vice versa. It’s ok to use Newton’s
Forward or Backward Interpolation Method to estimate 𝑓(𝑥 ). Where 𝑥 may
lie inside or outside the range. Before getting into details of methods, first, we
need to learn Finite Difference operators.

𝟐. 𝟏 Finite Difference Operator

For 𝑦 = 𝑓(𝑥) given data (𝑥 , 𝑥 , 𝑥 … 𝑥 ) and corresponding 𝑦 = 𝑓(𝑥 )
(𝑦 , 𝑦 , 𝑦 … 𝑦 )

such that 𝑥 = 𝑥 + 𝑛ℎ, where 𝑛 ≥ 1 and ℎ is the length of the interval which is
constant.
Shifting Operator (E)
It defines as 𝐸𝑓(𝑥) = 𝑓(𝑥 + ℎ) similarly 𝐸 𝑓(𝑥) = 𝐸 𝐸𝑓(𝑥) = 𝐸 𝑓(𝑥 + ℎ) =
𝑓(𝑥 + 2ℎ), 𝐸 𝑓(𝑥) = 𝑓(𝑥 − 2ℎ),𝐸 𝑓(𝑥) = 𝑓 𝑥 − .
General Formula
𝐸 𝑓(𝑥) = 𝑓(𝑥 + 𝑛ℎ) (2.1)

19
Forward Difference Operator (𝜟)
It defines as 𝛥𝑓(𝑥) = 𝑓(𝑥 + ℎ) − 𝑓(𝑥).
𝛥 𝑓(𝑥) = 𝛥 𝛥𝑓(𝑥)
= 𝛥 𝑓(𝑥 + ℎ) − 𝑓(𝑥)
= 𝛥𝑓(𝑥 + ℎ) − 𝛥𝑓(𝑥)
= 𝑓(𝑥 + 2ℎ) − 𝑓(𝑥 + ℎ) − ( 𝑓(𝑥 + ℎ) − 𝑓(𝑥))
= 𝑓(𝑥 + 2ℎ) − 2𝑓(𝑥 + ℎ) + 𝑓(𝑥)

𝛥𝑦 = 𝛥𝑓(𝑥 )
= 𝑓(𝑥 + ℎ) − 𝑓(𝑥 )
= 𝑓(𝑥 ) − 𝑓(𝑥 )
𝛥𝑦 = 𝑦 − 𝑦
Similarly
𝛥𝑦 = 𝑦 − 𝑦
𝛥𝑦 = 𝑦 − 𝑦

General Formula
𝛥𝑦 = 𝑦 −𝑦 (2.2)

Backward Difference Operator (𝜵)

It defines as 𝛻𝑓(𝑥) = 𝑓(𝑥) − 𝑓(𝑥 − ℎ).

𝛻𝑦 = 𝛻𝑓(𝑥 )
= 𝑓(𝑥 ) − 𝑓(𝑥 − ℎ)
= 𝑓(𝑥 ) − 𝑓(𝑥 )
= 𝛻𝑦 = 𝑦 − 𝑦
General Formula
𝛻𝑦 = 𝑦 − 𝑦 (2.3)
The relation between Finite Difference Operator
𝛥=𝐸−1 (2.4)
𝛻 = 1−𝐸 (2.5)

𝟐. 𝟐 Newton’s Gregory Forward Interpolation Method

Let for 𝑦 = 𝑓(𝑥) a given data point is (𝑦 , 𝑦 , 𝑦 … . 𝑦 ), (𝑥 , 𝑥 , 𝑥 … . 𝑥 ) . In
this method we estimate 𝑓(𝑥 ) or 𝑦 where 𝑥 𝜖 (𝑥 , 𝑥 ). Data points of the
independent variable must be equally spaced (|𝑥 − 𝑥 | = 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡). So 𝑥 =
𝑥 + 𝑛ℎ , 𝑛 = 0,1,2,3, … and ℎ is some constant. Finding 𝑦 = 𝑓(𝑥) for
𝑥 𝜖 (𝑥 , 𝑥 ).

20
Let 𝑥 = 𝑥 + 𝑢ℎ or 𝑢 = (2.6)

𝐸 𝑓(𝑥) = 𝑓(𝑥 + 𝑢ℎ)

𝐸 𝑓(𝑥 ) = 𝑓(𝑥 + 𝑢ℎ)
Using (2.6)
𝑦 = 𝑓(𝑥) = 𝐸 𝑦
Using (2.4)
𝑦 = (1 + 𝛥) 𝑦 (2.7)
Binomial expansion of (1 + 𝛥) is
𝑢𝛥 𝑢(𝑢 − 1)𝛥 𝑢(𝑢 − 1)(𝑢 − 2)𝛥
(1 + 𝛥) = 1 + + + +⋯ (2.8)
1! 2! 3!
From 2.7 and 2.8 equation

𝑢𝛥𝑦 𝑢(𝑢 − 1)𝛥 𝑦 𝑢(𝑢 − 1)(𝑢 − 2)𝛥 𝑦

𝑦=𝑦 + + + +⋯ (2.9)
1! 2! 3!

Our equation (2.9) is Newton’s Forward Interpolation formula, where 𝑢 =

.

Formation of Forward difference table using 𝜟𝒚𝒏 = 𝒚𝒏 𝟏 − 𝒚𝒏 .

Start the calculations from the first data point.

𝒙 𝒚 𝜟𝒚 𝜟𝟐 𝒚 𝜟𝟑 𝒚

𝑥 𝑦 𝛥𝑦 = 𝑦 − 𝑦 𝛥 𝑦 = 𝛥𝑦 − 𝛥𝑦 𝛥 𝑦 = 𝛥 𝑦 − 𝛥 𝛥𝑦

𝑥 𝑦 𝛥𝑦 = 𝑦 − 𝑦 𝛥 𝑦 = 𝛥𝑦 − 𝛥𝑦

𝑥 𝑦 𝛥𝑦 = 𝑦 − 𝑦

𝑥 𝑦

Table 2.1

Example 2.1 Find the cubic polynomial which takes the following values.
𝑥 0 1 2 3
𝑓(𝑥) 1 2 1 10
Find 𝑓(4) and 𝑓 (4).

The Forward difference table is given from Table 2.1

21
 𝒙 𝒚 𝜟𝒚 𝜟𝟐 𝒚 𝜟𝟑 𝒚

0 1 𝛥𝑦 =1 𝛥 𝑦 = −2 𝛥 𝑦 = 12

1 2 𝛥𝑦 = −1 𝛥 𝑦 = 10

2 1 𝛥𝑦 = 9

3 10

Table 2.2

Here ℎ = 1 ,𝑥 = 0

Using (2.6) ,(2.9) and Table 2.2

𝑢 = 𝑥 and 𝑓(𝑥) = 2𝑥 − 7𝑥 + 6𝑥 + 1
𝑓 (𝑥) = 6𝑥 − 14𝑥 + 6
𝑓(4) = 41
𝑓 (4) = 46
Scilab Code 𝟐. 𝟏

//NEWTON'S FORWARD INTERPOLATION

clc
clear
//matrix u store table 2.1
n=input(' input the number of data point ')
u=zeros(n,n+1)
u(:,1)=input(' input x data')
u(:,2)=input(' input y data')
h=abs(u(2,1)-u(1,1))//step size
x=input(' input interpolating value ')

a=(x-u(1,1))/h

for j=3:n+1//Formation of Forward difference table

for i=1:(n-j+2)
u(i,j)=u(i+1,j-1)-u(i,j-1)
end
end
disp(u)
//s variable is for summation of terms in equation (2.9)
22
s=u(1,2)
t=1
for i=0:n-2
t=(a-i)*t
s=s+t*u(1,i+3)/factorial(i+1)
end
disp(' value of y at x='+string(x)+' is '+string(s)+'')
Console Window 𝟐. 𝟏

𝟐. 𝟑 Newton’s Gregory Backward Interpolation Method

Let for 𝑦 = 𝑓(𝑥) a given data point is (𝑦 , 𝑦 , 𝑦 … . 𝑦 ), (𝑥 , 𝑥 , 𝑥 … . 𝑥 ) . In
this method we estimate 𝑓(𝑥 ) or 𝑦 where 𝑥 𝜖 (𝑥 , 𝑥 ). Data point of the
independent variable must be equally spaced (|𝑥 − 𝑥 | = 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡). So
𝑥 = 𝑥 + 𝑛ℎ , 𝑛 = 0,1,2,3, … and ℎ is some constant. Finding 𝑦 = 𝑓(𝑥) for
𝑥 𝜖 (𝑥 , 𝑥 ). Let 𝑥 = 𝑥 + 𝑢ℎ or 𝑢 = (2.10)
𝐸 𝑓(𝑥) = 𝑓(𝑥 + 𝑢ℎ)
𝐸 𝑓(𝑥 ) = 𝑓((𝑥 + 𝑢ℎ)
Using (2.10)
𝑓(𝑥) = 𝐸 𝑦
𝑓(𝑥)=(𝐸 ) 𝑦
Using (2.5)
𝑓(𝑥) = (1 − 𝛻) 𝑦 (2.11)
Binomial expansion of (1 − 𝛻)
𝑢𝛻 𝑢(𝑢 + 1)𝛻 𝑢(𝑢 + 1)(𝑢 + 2)𝛻
(1 − 𝛻) = 1 + + + +⋯ (2.12)
1! 2! 3!

23
Using (2.11) and (2.12)

𝑢𝛻𝑦 𝑢(𝑢 + 1)𝛻 𝑦 𝑢(𝑢 + 1)(𝑢 + 2)𝛻 𝑦

𝑓(𝑥) = 𝑦 + + + +⋯ (2.13)
1! 2! 3!

Our equation (2.13) is Newton’s Backward interpolation formula.

Where 𝑢 = .

Formation of Backward difference table using 𝜵𝒚𝒏 = 𝒚𝒏 − 𝒚𝒏 𝟏

Starts the calculations from the last data point.

𝒙 𝒚 𝜵𝒚 𝜵𝟐 𝒚 𝜵𝟑 𝒚

𝑥 𝑦

𝑥 𝑦 𝛻𝑦 = 𝑦 − 𝑦

𝑥 𝑦 𝛻𝑦 = 𝑦 − 𝑦 𝛻 𝑦 = 𝛻𝑦 − 𝛻𝑦

𝑥 𝑦 𝛻𝑦 = 𝑦 − 𝑦 𝛻 𝑦 = 𝛻𝑦 − 𝛻𝑦 𝛻 𝑦 =𝛻 𝑦 −𝛻 𝑦

Table 2.3

Example 2.2 Using Newton’s Backward Interpolation formula calculate

the population of the city in year 1925 from given data.

Year (𝑥) 1891 1901 1911 1921 1931

Population of
city in 46 66 81 93 101
thousands (𝑦)

Table 2.5
Backward difference table from Table 2.5
𝒙 𝒚 𝜵𝒚 𝜵𝟐 𝒚 𝜵𝟑 𝒚 𝜵𝟒 𝒚
1891 46

1901 66 𝛻𝑦 = 20

24
 1911 81 𝛻𝑦 = 15 𝛻 𝑦 = −5
1921 93 𝛻𝑦 = 12 𝛻 𝑦 = −3 𝛻 𝑦 =2
1931 101 𝛻𝑦 = 8 𝛻 𝑦 = −4 𝛻 𝑦 = −1 𝛻 𝑦 = −3
Table 2.4
Using Table 2.4 and equation (2.13)
8𝑢 4𝑢(𝑢 + 1) 𝑢(𝑢 + 1)(𝑢 + 2) 3𝑢(𝑢 + 1)(𝑢 + 2)(𝑢 + 3)
𝑓(𝑥) = 1 + − − −
1! 2! 3! 4!
Here ℎ = 10 and 𝑥 = 1931
Using (2.10)
𝑢= ,𝑓(1925) = 96.8368
Scilab Code 𝟐. 𝟐
//NEWTON'S Backward INTERPOLATION
clc
clear
n=input('input the number of data point')
//matrix u store table 2.3
u=zeros(n,n+1)
u(:,1)=input(' input x data')
u(:,2)=input(' input y data')
h=abs(u(2,1)-u(1,1))//step size
x=input(' input interpolating value ')

a=(x-u(n,1))/h

for j=3:n+1//Formation of Backward difference table

for i=j-1:n
u(i,j)=u(i,j-1)-u(i-1,j-1)
end
end
disp(u)
//s variable is for summation of terms in equation (2.13)
s=u(n,2)
t=1
for i=0:n-2
t=(a+i)*t
s=s+t*u(n,i+3)/factorial(i+1)
end
disp(' value of y at x='+string(x)+' is '+string(s)+'')

25
Console Window 𝟐. 𝟐

𝟐. 𝟒 Lagrange’s Interpolation Method

Consider a set of 𝑘 + 1 data points of the unequal interval. ( 𝑥 − 𝑥 =
𝑛𝑜𝑡 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡)
(𝑥 , 𝑦 ), … . 𝑥 , 𝑦 … . (𝑥 , 𝑦 )
Where no two 𝑥 are the same. The interpolating function is given by

𝑓(𝑥) = 𝑦 𝑙 (𝑥)

Where
𝑥−𝑥
𝑙 (𝑥) =
𝑥 −𝑥

Let’s take 3 data (𝑥 , 𝑦 ), (𝑥 , 𝑦 ), (𝑥 , 𝑦 )

(𝑥 − 𝑥 )(𝑥 − 𝑥 ) (𝑥 − 𝑥 )(𝑥 − 𝑥 ) (𝑥 − 𝑥 )(𝑥 − 𝑥 )
𝑓(𝑥) = 𝑦 + 𝑦 + 𝑦
(𝑥 − 𝑥 )(𝑥 − 𝑥 ) (𝑥 − 𝑥 )(𝑥 − 𝑥 ) (𝑥 − 𝑥 )(𝑥 − 𝑥 )
Inverse Lagrange’s Interpolation method
This method is used to calculate 𝑥 for given 𝑓(𝑥) where 𝑥 𝜖 (𝑥 , 𝑥 )

𝑥= 𝑥 𝑚 (𝑦)

Where

26
 𝑦−𝑦
𝑚 (𝑦) =
𝑦 −𝑦

Example 𝟐. 𝟑 Compute 𝑓(0.3) for the data

𝑥 0 1 3 4 7
𝑓 1 3 49 129 813
by Lagrange's interpolation formula using Scilab.
Scilab Code 𝟐. 𝟑
//Lagrange's Interpolation Method
clc
clear
k=input(' input number of data point ')
x=input(' input x data points ')
y=input(' input y data points ')
a=input(' input interpolating value ')
s=0
for j=1:k
l(j)=1
for i=1:k
if j~=i then
l(j)=((a-x(i))/(x(j)-x(i)))*l(j)
end
end
s=s+y(j)*l(j)
end
disp(' value of y at x='+string(a)+' is '+string(s)+'')

Console Window 𝟐. 𝟑

27
CHAPTER
Curve Fitting
We do experiments and observe data in the form of (𝑥 , 𝑓(𝑥 )) most of the
time. Function 𝑓(𝑥) is required for further calculation.

Let’s consider a simple experiment which is about finding the resistance of

ohmic material. Where we are observing voltage for different current. We
know that relation between 𝑉 and 𝐼 is a straight line and slope corresponds to
the resistance of ohmic material. Now we need to find 𝑓(𝑥) or a straight for
finding resistance.

We can connect all data point (𝑉 , 𝐼 ) to get a straight line. But in reality,
connecting all data point does not give a straight line as there are many
sources of error which do not lead to a straight line. So, we need to find the
approximate straight line and the method by which we find the approximate
straight line from given (𝑥 , 𝑓(𝑥 )) data is called Curve Fitting or Line fitting
(in this case we are fitting into line curve). The 2nd kind of curve fitting is
parabola fitting for a physical problem which is described by the parabolic
equation. E.g. Length vs Time period in simple pendulum [ 𝐿 = ].

𝟑. 𝟏 Line Fitting
Let we have (𝑥 , 𝑦 ) data points for finding the best fit line 𝑦 = 𝑎𝑥 + 𝑏. To find
𝑎 and 𝑏 we need to solve the system of the linear equation such that
𝑎∑𝑥 + 𝑛𝑏 = ∑𝑦
𝑎∑𝑥 + 𝑏∑𝑥 = ∑𝑥𝑦
Where 𝑛 is no. of data points.
Scilab Code 𝟑. 𝟏
//Line Fitting
clc
clear
clf
n=input(' input number of data points ')
x=input(' input x data points')
y=input(' input y data points')
plot(x,y,'*r')
for i=1:n
z(i)=x(i)*y(i)

28
end
sx=sum(x)
sy=sum(y)
sz=sum(z)
sx2=sum(x^2)
A=[sx n;sx2 sx]
B=[sy;sz]
C=inv(A)*B
a=C(1,1)//a is the slope of the line
b=C(2,1)//b is the y-intercept of the line
disp(''+string(a)+' is the slope of the line')
disp(''+string(b)+' is the y-intercept of the line')
for i=1:n
t(i)=a*x(i)+b
end
plot(x,t,'m')
xlabel('X','Fontsize',4)
ylabel('Y','Fontsize',4)
legend('Data points','Line fit curve',-4)
xgrid(3)

Console Window 𝟑. 𝟏

29
Graphic Window 𝟑. 𝟏

𝟑. 𝟐 Parabola Fitting
Let we need to find best fit parabola 𝑦 = 𝑎𝑥 + 𝑏𝑥 + 𝑐 for given data point
(𝑥 , 𝑦 ) . To find 𝑎, 𝑏 and 𝑐 we need to solve the system of the linear equation
such that
𝑎∑𝑥 + 𝑏∑𝑥 + 𝑛𝑐 = ∑𝑦
𝑎∑𝑥 + 𝑏∑𝑥 + 𝑐∑𝑥 = ∑𝑥𝑦
𝑎∑𝑥 + 𝑏∑𝑥 + 𝑐∑𝑥 = ∑𝑥 𝑦
Where 𝑛 is no. of data points.
Scilab Code 𝟑. 𝟐
//Parabola Fitting
clc
clear
clf
n=input(' input number of data points ')
x=input(' input x data points ')
y=input(' input y data points ')
plot(x,y,'*b')
for i=1:n
z1(i)=x(i)*y(i)

30
 z2(i)=x(i)^2*y(i)
end
sx=sum(x)
sx2=sum(x^2)
sy=sum(y)
sx3=sum(x^3)
sx4=sum(x^4)
sxy=sum(z1)
sx2y=sum(z2)
A=[sx2 sx n;sx3 sx2 sx;sx4 sx3 sx2]
B=[sy;sxy;sx2y]
C=inv(A)*B
a=C(1,1)//y=ax^2+bx+c
b=C(2,1)
c=C(3,1)
disp('coefficient of x^2 '+string(a)+'')
disp('coefficient of x '+string(b)+'')
disp('constant term '+string(c)+'')
for i=1:n
t(i)=a*x(i)^2+b*x(i)+c
end
plot(x,t,'c')
xlabel('X','Fontsize',4)
ylabel('Y','Fontsize',4)
legend('Data points','Parabola fit curve',-4)
xgrid(3)

Console Window 𝟑. 𝟐

31
Graphic Window 𝟑. 𝟐

32
CHAPTER
Numerical Differentiation
Numerical Differentiation is about finding 𝑦 (𝑥) for given (𝑥 , 𝑦 ) data points.

𝟒. 𝟏 Equal Interval
When |𝑥 − 𝑥 | = 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡, Newton Forward and Newton Backward
interpolation formula are used to construct 𝑓(𝑥) for given (𝑥 , 𝑦 ) data points.
So, we can deduce the Numerical Differentiation formula from Newton
Forward and Newton Backward interpolation formula. Instead of
constructing 𝑓(𝑥) and differentiating it to find 𝑓 (𝑥) we will construct a
formula for 𝑓 (𝑥) from interpolating formula.
Forward Numerical Differentiation Formula
From Newton Forward interpolation formula equation (2.9)
𝑢𝛥𝑦 𝑢(𝑢 − 1)𝛥 𝑦 𝑢(𝑢 − 1)(𝑢 − 2)𝛥 𝑦
𝑦(𝑥) = 𝑦 + + + +⋯
1! 2! 3!
Where 𝑢 = ,𝑥 = 𝑥 + 𝑢ℎ

Substituting 𝑥 with 𝑥 + 𝑢ℎ in (2.9)

𝑢𝛥𝑦 (𝑢 − 𝑢)𝛥 𝑦 (𝑢 − 3𝑢 + 2𝑢)𝛥 𝑦
𝑦(𝑥 + 𝑢ℎ) = 𝑦 + + + +⋯ (4.1)
1! 2! 3!
Differentiating (4.1) with respect to 𝑢

𝑑𝑦 𝛥𝑦 (2𝑢 − 1)𝛥 𝑦 (3𝑢 − 6𝑢 + 2)𝛥 𝑦

𝑦 (𝑥) = = + + +⋯ (4.2)
𝑑𝑥 1! ℎ 2! ℎ 3! ℎ

Equation (4.2) is our formula for Numerical Differentiation from the Newton
Forward Interpolation formula.

Backward Numerical Differentiation Formula

From Newton Backward Interpolation formula equation (2.13)

𝑢𝛻𝑦 𝑢(𝑢 + 1)𝛻 𝑦 𝑢(𝑢 + 1)(𝑢 + 2)𝛻 𝑦
𝑦(𝑥) = 𝑦 + + + +⋯
1! 2! 3!

Where 𝑢 = ,𝑥 = 𝑥 + 𝑢ℎ

Substituting 𝑥 with 𝑥 + 𝑢ℎ in (2.13)

33
 𝑢𝛻𝑦 (𝑢 + 𝑢)𝛻 𝑦 (𝑢 + 3𝑢 + 2𝑢)𝛻 𝑦
𝑦(𝑥 + 𝑢ℎ) = 𝑦 + + + +⋯ (4.3)
1! 2! 3!
Differentiating (4.3) with respect to 𝑢

𝑑𝑦 𝛻𝑦 (2𝑢 + 1)𝛻 𝑦 (3𝑢 + 6𝑢 + 2)𝛻 𝑦

𝑦 (𝑥) = = + + +⋯ (4.4)
𝑑𝑥 1! ℎ 2! ℎ 3! ℎ

Equation (4.4) is our formula for Numerical Differentiation from the Newton
Backward Interpolation formula.

𝟒. 𝟐 Unequal Interval
We have seen in chapter 2 Lagrange’s interpolation method is used to
construct 𝑓(𝑥) for an unequal interval (|𝑥 − 𝑥 | = 𝑛𝑜𝑡 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡) for given
data point (𝑥 , 𝑦 ). In unequal interval data points, we don’t derive a separate
formula for numerical differentiation. We just differentiate 𝑓(𝑥) by
constructing it from interpolating formula to get 𝑓 (𝑥).

34
CHAPTER
Numerical Integration
Numerical Integration is about solving definite integral by discretization of
integral limit or finding ∫ 𝑦 𝑑𝑥 for given (𝑥 , 𝑦 ).
Before jumping into Numerical Integration. First, we must have to understand
the basic definition of Integration.
We are given to evaluate 𝐼 = ∫ 𝑦 𝑑𝑥 , 𝑦 = 𝑓(𝑥). That means we need to
calculate area bounded by 𝑦 = 𝑓(𝑥), 𝑦 = 0, 𝑥 = 𝑎 and 𝑥 = 𝑏.
Solving integral by slicing bounded area into Rectangular strip (Figure 5.1).
Then calculating the individual area of rectangular strip and adding them to
get the total area. So, we need to choose the width of the rectangular strip (ℎ =
|𝑥 − 𝑥 |)

Obviously, the length of the rectangle is 𝑓(𝑥 ) at 𝑥 (Figure 5.1).

Considering a smaller value of ℎ leads to more accuracy (to make a perfect

rectangle).

Number of discrete point in x-axis is 𝑛 = + 1.

Figure 5.1
There are many methods for Numerical Integration but we will focus on the
listed below:
(a) Trapezoidal rule

35
(b) Simpson's 1/3 rule
(c) Simpson's 3/8 rule

All three formulas are deduced from Newton Cotes formula and this formula
is derived from the Newton Forward Interpolation formula.

𝟓. 𝟏 Newton-Cotes Formula
Let 𝑦 = 𝑓(𝑥) be a function such that 𝑦 = 𝑓(𝑥 ) , 𝑦 = 𝑓(𝑥 ) … .., 𝑦 = 𝑓(𝑥 ).

Where 𝑥 = 𝑥 + 𝑛ℎ, ℎ is the step size.

From Newton’s Forward Interpolation formula equation (2.9)

𝑢𝛥𝑦 𝑢(𝑢 − 1)𝛥 𝑦 𝑢(𝑢 − 1)(𝑢 − 2)𝛥 𝑦
𝑓(𝑥) = 𝑦 = 𝑦 + + + +⋯
1! 2! 3!
Where 𝑢 = ,𝑥 = 𝑥 + 𝑢ℎ

Now integrate ∫ 𝑦 𝑑𝑥 where 𝑎 = 𝑥 , 𝑏 = 𝑥

∫ 𝑦 𝑑𝑥 = ∫ 𝑦 𝑑𝑥 (5.1)

Substituting 𝑥 with 𝑥 + 𝑢ℎ in (2.9)

𝑑𝑥 = ℎ𝑑𝑢 and 𝑢 = 0 to 𝑢 = 𝑛

From (2.9) and (5.1)

𝑦 𝑑𝑥 = ℎ 𝑦(𝑥 + 𝑢ℎ) 𝑑𝑢

𝑢𝛥𝑦 (𝑢 − 𝑢)𝛥 𝑦 (𝑢 − 3𝑢 + 2𝑢)𝛥 𝑦

=ℎ (𝑦 + + + … ) 𝑑𝑢
1! 2! 3!

𝑛 𝛥𝑦 (2𝑛 − 3𝑛 )𝛥 𝑦 (𝑛 − 4𝑛 + 4𝑛 )𝛥 𝑦
𝑦 𝑑𝑥 = ℎ[(𝑛𝑦 + + + … )] (5.2)
2 12 24

The above equation (5.2) is our Newton-Cotes formula as a parent formula for
Numerical Integration methods.

36
𝟓. 𝟐 Trapezoidal Rule
Take 𝑛 = 1(this mean only two data point(𝑥 , 𝑦 ), (𝑥 , 𝑦 )). Substitute into
newton-cotes formula (5.2)

𝒙 𝒚 𝜟𝒚

𝑥 𝑦 𝛥𝑦 = 𝑦 − 𝑦

𝑥 𝑦

From the above table, we can see that only 𝛥𝑦 is defined. So, terms containing
𝛥 𝑦 , 𝛥 𝑦 …are redundant.

1 1 ℎ
𝑦 𝑑𝑥 = ℎ 1𝑦 + 𝛥𝑦 = ℎ 𝑦 + (𝑦 − 𝑦 ) = [𝑦 + 𝑦 ]
2 2 2

Similarly,

∫ 𝑦 𝑑𝑥 = [𝑦 + 𝑦 ],∫ 𝑦 𝑑𝑥 = [𝑦 + 𝑦 ]….∫ 𝑦 𝑑𝑥 = [𝑦 +𝑦 ]
∫ 𝑦 𝑑𝑥 =∫ 𝑦 𝑑𝑥 + ∫ 𝑦 𝑑𝑥 + ⋯ ∫ 𝑦 𝑑𝑥

= [𝑦 + 𝑦 + 𝑦 + 𝑦 𝑦 … … . 𝑦 +𝑦 +𝑦 ]

ℎ
𝑦 𝑑𝑥 = [𝑦 + 2(𝑦 + 𝑦 … . +𝑦 )+𝑦 ] (5.3)
2

The above equation (5.3) is the Trapezoidal formula for Numerical

Integration. In this formula, any value of 𝑛 is allowed.

37
Scilab Code 𝟓. 𝟏

//Numerical Integration by Trapezoidal rule

clc
clear
function z=f(x)
z=1/(1+x^2)
endfunction
disp('f(x)=1/(1+x^2)')
a=input(' enter intial value ')
b=input(' enter final value ')
n=input(' input number of discrete points ')
h=(b-a)/(n-1)//step size
x(1)=a
y(1)=f(a)
for i=2:n
x(i)=x(i-1)+h
y(i)=f(x(i))
end
I=(sum(y)*h)-(h*0.5*(y(1)+y(n)))
disp('value of ∫1/(1+x^2)dx limit '+string(a)+' to '+string(b)+' by Trapezoidal
rule is '+string(I)+'')
z=intg(a,b,f)
disp('exact value of ∫1/(1+x^2)dx limit '+string(a)+' to '+string(b)+' =
'+string(z)+'' )

Console Window 𝟓. 𝟏

38
𝟓. 𝟑 Simpson's 𝟏/𝟑 rule
Take 𝑛 = 2 and substitute into Newton-Cotes formula (5. 2)

data points (𝑥 , 𝑦 ), (𝑥 , 𝑦 ), (𝑥 , 𝑦 )

𝒙 𝒚 𝜟𝒚 𝜟𝟐 𝒚
𝑥 𝑦 𝛥𝑦 = 𝑦 − 𝑦 𝛥 𝑦 = 𝛥𝑦 − 𝛥𝑦
𝑥 𝑦 𝛥𝑦 = 𝑦 − 𝑦
𝑥 𝑦
From the above table, we can see that term containing 𝛥 𝑦 , 𝛥 𝑦 … .. are
redundant in newton cote’s formula.

4 1
𝑦 𝑑𝑥 = ℎ 2𝑦 + 𝛥𝑦 + (16 − 12)𝛥 𝑦
2 12

1
= ℎ 2𝑦 + 2(𝑦 −𝑦 ) + 𝛥(𝑦 − 𝑦 )
3
1
= ℎ 2𝑦 + 2(𝑦 −𝑦 ) + (𝑦 − 2𝑦 + 𝑦 )
3
ℎ
= (𝑦 +4𝑦 + 𝑦 )
3
Similarly,

∫ 𝑦 𝑑𝑥 = (𝑦 +4𝑦 + 𝑦 ) ,∫ 𝑦 𝑑𝑥 = (𝑦 +4𝑦 + 𝑦 ) …

ℎ
𝑦 𝑑𝑥 = (𝑦 +4𝑦 +𝑦 )
3

∫ 𝑦 𝑑𝑥 =∫ 𝑦 𝑑𝑥 + ∫ 𝑦 𝑑𝑥 + ⋯ ∫ 𝑦 𝑑𝑥

(5.4)
ℎ
𝑦 𝑑𝑥 = (𝑦 + 4(𝑦 + 𝑦 + 𝑦 … 𝑦 ) + 2(𝑦 + 𝑦 + 𝑦 … 𝑦 )+𝑦 )
3

The above equation (5.4) is Simpson's 1/3 rule for Numerical Integration. In
this formula, only multiple of 2 is allowed for 𝑛.

39
𝟓. 𝟒 Simpson's 𝟑/𝟖 rule
By substituting 𝑛 = 3 in Newton-Cotes formula (5. 2). We get the formula for
Simpson’s 3/8 rule

∫ 𝑦 𝑑𝑥 = (𝑦 + 2(𝑦 + 𝑦 + 𝑦 … 𝑦 ) + 3(𝑦 + 𝑦 + 𝑦 … 𝑦 )+𝑦 ) (5.5)

2nd term in (5.5) containing 𝑦 where 𝑖 is multiple of 3.

3rd term in (5.5) contains 𝑦 other than 𝑦 , 𝑦 and 𝑦 from 2nd term.

Only multiple of 3 is allowed for 𝑛 in (5.5)

Example 5.1 Evaluate by using (i) Trapezoidal rule (ii) Simpson's

1/3 rule (iii) Simpson's 3/8 rule and verify analytically.

Take 𝑛 = 6 , ℎ = =1

𝒙𝒊 = 𝒙𝒊 + 𝒊𝒉 𝟏
𝒚𝒊 =
𝟏 + 𝒙𝟐𝒊
𝑥 =0 𝑦 =1
𝑥 =1 𝑦 = 0.5
𝑥 =2 𝑦 = 0.2
𝑥 =3 𝑦 = 0.1
𝑥 =4 𝑦 = 0.0588
𝑥 =5 𝑦 = 0.0385
𝑥 =6 𝑦 = 0.027

(a)By Trapezoidal rule

= [𝑦 + 2(𝑦 + 𝑦 + 𝑦 +𝑦 + 𝑦 ) + 𝑦 ]

= 1.4108
(b) By Simpson's 1/3 rule

= (𝑦 + 4(𝑦 + 𝑦 + 𝑦 ) + 2(𝑦 + 𝑦 ) + 𝑦 )

= 1.3662

40
(c) By Simpson's 3/8 rule

= (𝑦 + 2𝑦 + 3(𝑦 + 𝑦 + 𝑦 + 𝑦 ) + 𝑦 )

= 1.3571
Solving analytically

= [tan (𝑥)] = 1.4056

Example 5.2 A river is 80 𝑚 wide. The depth 𝑑 in meters at a distance 𝑥

from one bank is given by the following table. Find the approximate area of
the cross-section.

𝑥(𝑚) 0 10 20 30 40 50 60 70 80
𝑑(𝑚) 0 4 7 9 12 15 14 08 03

Figure 5.2

In this problem cross section area= ∫ 𝑦𝑑𝑥 , ℎ = 10 and 𝑛 = 8

Let’s solve above integral by Simpson’s 1/3 rule (5.4).

∫ 𝑦𝑑𝑥 = (0 + 4(4 + 9 + 15 + 8) + 2(7 + 12 + 14) + 3)

=710𝑚

41
𝟓. 𝟓 Monte Carlo method
Solving integral by Monte Carlo method is based on probability.

Let’s consider an example ∫ 𝑥 𝑑𝑥 to see how this method works

Figure 5.3

Solving ∫ 𝑥 𝑑𝑥 this mean calculating area under curve see Figure 5.3
((0,0), (2,0), (2,4))

Steps to implement Monte Carlo method

(a) Consider a Rectangle of size(2 × 4) Figure 5.3 ((0,0), (2,0), (2,4), (0,4).

Pour grain over entire rectangular region ((0,0), (2,0), (2,4), (0,4) uniformly.

(b) Count total no. of grains poured and grains inside the region
((0,0), (2,0), (2,4))

(c) Calculate proportion of the total grains under the curve

. (( , ),( , ),( , ))
(0,0), (2,0), (2,4)) = 𝑝 =
. .

(d)Calculating integral ∫ 𝑥 𝑑𝑥 = 𝑝 × 𝑎𝑟𝑒𝑎 𝑜𝑓 𝑟𝑒𝑐𝑡𝑎𝑛𝑔𝑢𝑙𝑎𝑟 𝑟𝑒𝑔𝑖𝑜𝑛

A higher number of pouring grains leads to more accurate result. Solving

integral by the Monte Carlo method manually is not a good idea. So, in

42
computers instead of pouring grains, we generate random numbers and count
them.

Scilab Code 𝟓. 𝟐

//Solving Integral by Monte Carlo Method

clc
clear
clf
inside = 0
n = 1000//total number of random variables
xgrid(4)
xlabel('x','Fontsize',4)
ylabel('y','Fontsize',4)
title('Solving Integral by Monte Carlo Method','Fontsize',3)
for i = 1:n
x = 2*rand()
y = 4*rand()
if y <= x^2 then
plot(x,y,'g.')
inside = inside + 1
else
plot(x,y,'y.')
end
end
a=(inside/n)*8
disp('value of ∫x^2dx limit 0 to 2 is '+string(a)+'')

Console Window 𝟓. 𝟐

43
Graphic Window 𝟓. 𝟐

44
CHAPTER
Ordinary Differential Equation
Many physical problems are expressed in mathematical function and equation
involving many variables. These Expressions are called a mathematical model
for that physical system, any such mathematical model which involves
function and their derivatives are called differential equation.

Many real-world problems are modeled into the differential equation. A few
of them are listed below.

(a) Radioactive decay of Radioactive element such as uranium.

(b) Motion of body

(c) population Models

(d) prediction of weather

(e) prediction of stock prices

(f) fluid flow (flow in the river), Motion of airplane, blood flow, swimming of
animals, air flow in lungs etc.

(g) Atmospheric pollution

(h) Electric field and potentials

(i) Diffusion of Materials

(k) heat flow and distribution

(l) deformation of the body, elasticity

(m) Modeling cancer growth

(n) pattern formation (e.g. strips in zebra, spots in leopard)

(o) drug delivery to part of the body

(p) propagation of sound wave through medium

(q) Enzyme kinetics

(r) Epidemic models

45
Versatile use of differential equations motivates us to learn how to solve the
differential equation numerically since many differential equations cannot be
solved analytically.

𝟔. 𝟏 Classification of differential equations

(a) ordinary differential equation (ODE) is an equation that contains one or
more derivatives of an unknown function of a single independent variable.

E.g. (a) + =𝑥 (b) + 4𝑦 = 0

(b) Partial differential equation (PDE) is an equation involving one or more

partial derivates of an unknown function that depends on two or more
independent variables.
E.g. (a) + = 0 (b) =𝑐

The details of PDE and their numerical methods to solve are discussed in
chapter 8.

This chapter is dedicated for numerical methods to solve ODE.

General Terminology of ODE

Order- Highest order derivative present in ODE

E.g. 𝑥 − 3𝑥 + 6𝑥 − 6𝑦 = 0 is of order 3.

Degree-Power of highest order derivative in ODE.

E.g. +𝑥 = 0 is of order 2 and degree 2.

Linearity and Homogeneity of 𝒏𝒕𝒉 Order ODE

An 𝑛 order ODE in unknown function 𝑦(𝑥) is said to be linear if it can be

written in the form

𝑦( )
+𝑝 (𝑥)𝑦 ( )
+𝑝 (𝑥)𝑦 ( )
+ ⋯ 𝑝 (𝑥)𝑦 = 𝑟(𝑥)
𝑝 (𝑥), 𝑝 (𝑥), … 𝑝 (𝑥) are continuous in the defined domain. If 𝑟(𝑥) = 0 for
the defined domain then the above equation becomes Homogeneous
otherwise it is non-Homogeneous.

46
𝟔. 𝟐 Euler Method
Let’s consider 1st order ODE = 𝑓(𝑥, 𝑦) with an initial condition, 𝑦(𝑥 ) = 𝑦 .

The formula for the Euler method is

𝑦 = 𝑦 + 𝛥𝑥𝑓(𝑥 , 𝑦 ) (6.1)
We can calculate 𝑦 from previously calculated 𝑦 using (6.1). Smaller 𝛥𝑥
leads to higher accuracy for 𝑦 .

Example 6.1 Find 𝑦(2.2) from differential equation = −𝑥𝑦 with the
initial condition,𝑦(2) = 1.

Here, 𝑓(𝑥, 𝑦) = −𝑥𝑦 ,𝑥 = 2, 𝑦 = 1

Taking 𝛥𝑥 = 0.05
𝑦 = 𝑦 + 𝛥𝑥𝑓(𝑥 , 𝑦 )
𝑦 = 𝑦(2.05) = 0.9
𝑦 = 𝑦 + 𝛥𝑥𝑓(𝑥 , 𝑦 )
𝑦 = 𝑦(2.1) = 0.817
𝑦 = 𝑦 + 𝛥𝑥𝑓(𝑥 , 𝑦 )
𝑦 = 𝑦(2.15) = 0.747
𝑦 = 𝑦 + 𝛥𝑥𝑓(𝑥 , 𝑦 )
𝑦 = 𝑦(2.2) = 0.687
Higher order ODE and System of ODE can be solved using the Euler method.
The basic idea behind this is converting ODE into a system of first order ODE
and solving them using the Euler method.

2nd Order ODE

2nd order ODE 𝑥̈ (𝑡) = 𝑓(𝑡, 𝑥, 𝑥̇ ) (6.2)

where 𝑥̈ (𝑡) = , 𝑥̇ (𝑡) = with initial condition 𝑥(𝑡 ) = 𝑥 , 𝑥̇ (𝑡 ) = 𝑢

Taking new notation

𝑥̇ (𝑡) = 𝑢(𝑡) (6.3)
Using 6.2 and 6.3

47
𝑢̇ (𝑡) = 𝑓(𝑡, 𝑥, 𝑢) (6.4)
6.3 and 6.4 is system of 1st order ODE and can be solved using (6.1)

1st data point

𝑡 = 𝑡 + 𝛥𝑡
𝑥 = 𝑥 + 𝛥𝑡𝑢
𝑢 = 𝑢 + 𝛥𝑡𝑓(𝑡 , 𝑥 , 𝑢 )
2nd data point

𝑡 = 𝑡 + 𝛥𝑡
𝑥 = 𝑥 + 𝛥𝑡𝑢
𝑢 = 𝑢 + 𝛥𝑡𝑓(𝑡 , 𝑥 , 𝑢 )
.

.
Nth data point

𝑡 =𝑡 + 𝛥𝑡
𝑥 =𝑥 + 𝛥𝑡𝑢
𝑢 =𝑢 + 𝛥𝑡𝑓(𝑡 ,𝑥 ,𝑢 )

Example 𝟔. 𝟐 Differential equation for the driven damped harmonic

oscillator is given below. Describe the position and velocity of a particle by a
computer program. 𝑥̈ + 0.5𝑥̇ + 3𝑥 = 0.4 cos(3𝑡) with initial condition 𝑥(0) =
10, 𝑢(0) = 0
Constructing the system of 1st order ODE
𝑥̇ (𝑡) = 𝑢(𝑡)
𝑢̇ (𝑡) = 0.4 cos(3𝑡) − 3𝑥 − 0.5𝑢
Taking 𝛥𝑡 = 0.01

48
Scilab Code 𝟔. 𝟏

//Example 6.2
clc
clear
clf
function z=f(t, x, u)
z=0.4*cos(3*t)-0.5*u-3*x
endfunction
t=[0:0.01:20]'//time(s)
x(1)=10//intial postion (m)
u(1)=0//intial velocity(m/s)
for i=1:2000//solving system of 1st order ODE
x(i+1)=x(i)+0.01*u(i)
u(i+1)=u(i)+0.01*f(t(i),x(i),u(i))
end
subplot(2,1,1)
plot(t,x,'b')
xgrid(0)
xlabel('t(s)','Fontsize',3)
ylabel('x(m)','Fontsize',3)
title('Postion vs time ','Fontsize',3)
subplot(2,1,2)
plot(t,u,'r')
xgrid(0)
xlabel('t(s)','Fontsize',3)
ylabel('v(m/s)','Fontsize',3)
title('Velocity vs time ','Fontsize',3)

49
Graphic Window 6.1

System of 𝟏st order ODE

System of 1st order ODE 𝑥̇ (𝑡) = 𝑓 (𝑡, 𝑥, 𝑦)

𝑦̇ (𝑡) = 𝑓 (𝑡, 𝑥, 𝑦)
with the initial condition, 𝑥(𝑡 ) = 𝑥 𝑎𝑛𝑑 𝑦(𝑡 ) = 𝑦

Where 𝑥̇ (𝑡) = and 𝑦̇ (𝑡) =

1st data point

𝑡 = 𝑡 + 𝛥𝑡
𝑥 = 𝑥 + 𝛥𝑡𝑓1 (𝑡 , 𝑥 , 𝑦 )
𝑦 = 𝑦 + 𝛥𝑡𝑓2 (𝑡 , 𝑥 , 𝑦 )
2nd data point

𝑡 = 𝑡 + 𝛥𝑡
𝑥 = 𝑥 + 𝛥𝑡𝑓1 (𝑡 , 𝑥 , 𝑦 )
50
𝑦 = 𝑦 + 𝛥𝑡𝑓2 (𝑡 , 𝑥 , 𝑦 )

.
Nth data point

𝑡 =𝑡 + 𝛥𝑡
𝑥 =𝑥 + 𝛥𝑡𝑓1 (𝑡 ,𝑥 ,𝑦 )
𝑦 =𝑦 + 𝛥𝑡 𝑓2 (𝑡 ,𝑥 ,𝑦 )

Example 𝟔. 𝟑 Solve the given initial value problem using the Euler method
in Scilab.
𝑑𝑥
= 3𝑥 − 5𝑦
𝑑𝑡
𝑑𝑦
= 𝑥−𝑦
𝑑𝑡
𝑥(0) = 3, 𝑦(0) = 1
Taking 𝛥𝑡 = 0.001

𝑥 = 3, 𝑦 = 1, 𝑡 = 0
𝑖 = 0,1,2,3 …
𝑡 = 𝑡 + 0.001
𝑥 = 𝑥 + 0.001 ⋅ (3𝑥 − 5𝑦 )
𝑦 = 𝑦 + 0.001 ⋅ (𝑥 − 𝑦 )

51
Scilab Code 𝟔. 𝟐

//Example 6.3
clc
clear
clf
t=[0:0.001:20]'//t
//intial value
x(1)=3//x(0)=3
y(1)=1//y(0)=1
for i=1:20000//solving system of 1st order ODE
x(i+1)=x(i)+0.001*(3*x(i)-5*y(i))
y(i+1)=y(i)+0.001*(x(i)-y(i))
end
subplot(2,1,1)
plot(t,x,'r')
xgrid(0)
xlabel('t','Fontsize',3)
ylabel('x','Fontsize',3)
subplot(2,1,2)
plot(t,y,'b')
xgrid(0)
xlabel('t','Fontsize',3)
ylabel('y','Fontsize',3)

52
Graphic Window 𝟔. 𝟐

Higher order ODE

Consider the 3rd order ODE

𝑦 (𝑡) = 𝑓(𝑡, 𝑦, 𝑦 , 𝑦 ) (6.5)

with initial condition 𝑦(𝑡0 ) = 𝑦0 , 𝑦 (𝑡0 ) = 𝑢0 , 𝑦 (𝑡0 ) = 𝑎0

Where 𝑦 (𝑡) = , 𝑦 (𝑡) = , 𝑦 (𝑡) =

Taking new notation

𝑦 (𝑡) = 𝑢(𝑡) (6.6)
𝑦 (𝑡) = 𝑢 (𝑡) = 𝑎(𝑡) (6.7)

Using 6.5, 6.6 and 6.7

𝑎 (𝑡) = 𝑓(𝑡, 𝑦, 𝑢, 𝑎) (6.8)
(6.6), (6.7) and (6.8) is system of 1st order ODE and can be solved using (6.1)

53
1st data point

𝑡 = 𝑡 + 𝛥𝑡
𝑦 = 𝑦 + 𝛥𝑡𝑢
𝑢 = 𝑢 + 𝛥𝑡𝑎
𝑎 = 𝑎 + 𝛥𝑡𝑓(𝑡 , 𝑦 , 𝑢 , 𝑎 )
2nd data point

𝑡 = 𝑡 + 𝛥𝑡
𝑦 = 𝑦 + 𝛥𝑡𝑢
𝑢 = 𝑢 + 𝛥𝑡𝑎
𝑎 = 𝑎 + 𝛥𝑡𝑓(𝑡 , 𝑦 , 𝑢 , 𝑎 ).

.
Nth data point

𝑡 =𝑡 + 𝛥𝑡
𝑦 =𝑦 + 𝛥𝑡𝑢
𝑢 =𝑢 + 𝛥𝑡𝑎
𝑎 =𝑎 + 𝛥𝑡𝑓(𝑡 ,𝑦 ,𝑢 ,𝑎 )
We can extend the above notion to solve higher order ODE.

Example 𝟔. 𝟒 Solve the given initial value problem using the Euler method
in Scilab.

𝑦 − 𝑦 + 100𝑦 − 100𝑦 = 0 , 𝑦(0) = 4, 𝑦 (0) = 11, 𝑦 (0) = −299

Constructing the system of 1st order ODE

𝑦 (𝑡) = 𝑢(𝑡)
𝑢 (𝑡) = 𝑎(𝑡)
𝑎 (𝑡) = 𝑎 − 100𝑢 + 100𝑦

54
𝑦(0) = 4, 𝑢(0) = 11, 𝑎(0) = −299
Taking 𝛥𝑡 = 0.001

Scilab code 𝟔. 𝟑

//Example 6.4
clc
clear
clf
t=[0:0.001:5]'//t
//intial value
y(1)=4//y(0)=4
u(1)=11//u(0)=1
a(1)=-299// a(0)=-299
for i=1:5000//solving system of 1st order ODE
y(i+1)=y(i)+0.001*u(i)
u(i+1)=u(i)+0.001*a(i)
a(i+1)=a(i)+0.001*(a(i)-100*u(i)+100*y(i))
end
subplot(3,1,1)
plot(t,y,'b')
xgrid(0)
xlabel('t','Fontsize',3)
ylabel('y','Fontsize',3)
subplot(3,1,2)
plot(t,u,'r')
xgrid(0)
xlabel('t','Fontsize',3)
ylabel('u','Fontsize',3)
subplot(3,1,3)
plot(t,a,'r')
xgrid(0)
xlabel('t','Fontsize',3)
ylabel('a','Fontsize',3)

55
Graphic Window 𝟔. 𝟑

56
System of Higher order ODE

Let’s consider the system of second order ode

𝑥̈ (𝑡) = 𝑓 (𝑡, 𝑥, 𝑥̇ , 𝑦, 𝑦̇ )
𝑦̈ (𝑡) = 𝑓 (𝑡, 𝑥, 𝑥̇ , 𝑦, 𝑦̇ )

Where 𝑥̈ (𝑡) = , 𝑦̈ (𝑡) = , 𝑥̇ (𝑡) = , 𝑦̇ (𝑡) =

With initial condition 𝑥(𝑡 ) = 𝑥 , 𝑥̇ (𝑡) = 𝑢 , 𝑦(𝑡 ) = 𝑦 , 𝑦̇ (𝑡) = 𝑣

Considering new variable 𝑢 and 𝑣 to make system of first order ode

𝑥̇ (𝑡) = 𝑢(𝑡)
𝑦̇ (𝑡) = 𝑣(𝑡)
𝑢̇ (𝑡) = 𝑓 (𝑡, 𝑥, 𝑢, 𝑦, 𝑣)
𝑣̇ (𝑡) = 𝑓 (𝑡, 𝑥, 𝑢, 𝑦, 𝑣)

Applying Euler method in above system of first ODE

1st data point

𝑡 = 𝑡 + 𝛥𝑡
𝑥 = 𝑥 + 𝛥𝑡𝑢
𝑦 = 𝑦 + 𝛥𝑡𝑣
𝑢 = 𝑢 + 𝛥𝑡𝑓 (𝑡 , 𝑥 , 𝑢 , 𝑦 , 𝑣 )
𝑣 = 𝑣 + 𝛥𝑡𝑓 (𝑡 , 𝑥 , 𝑢 , 𝑦 , 𝑣 )
2nd data point
𝑡 = 𝑡 + 𝛥𝑡
𝑥 = 𝑥 + 𝛥𝑡𝑢
𝑦 = 𝑦 + 𝛥𝑡𝑣
𝑢 = 𝑢 + 𝛥𝑡𝑓 (𝑡 , 𝑥 , 𝑢 , 𝑦 , 𝑣 )
𝑣 = 𝑣 + 𝛥𝑡𝑓 (𝑡 , 𝑥 , 𝑢 , 𝑦 , 𝑣 )
.

57
.

Nth data point

𝑡 =𝑡 + 𝛥𝑡
𝑥 =𝑥 + 𝛥𝑡𝑢
𝑦 =𝑦 + 𝛥𝑡𝑣
𝑢 =𝑢 + 𝛥𝑡𝑓 (𝑡 ,𝑥 ,𝑢 ,𝑦 ,𝑣 )
𝑣 =𝑣 + 𝛥𝑡𝑓 (𝑡 ,𝑥 ,𝑢 ,𝑦 ,𝑣 )

We can extend the above notion to solve the system of higher order
ODE.

Example 𝟔. 𝟓 Equation of 2 mass coupled oscillator is given below with

initial condition, describe the position and velocity of both masses using
Scilab.

𝑥̈ (𝑡) = −6𝑥 + 2𝑦 , 𝑦̈ (𝑡) = 2𝑥 − 6𝑦 , 𝑥(0) = 1, 𝑥̇ (0) = 0, 𝑦(0) = 0, 𝑦̇ (0) = 0

Constructing system of 1st order ODE

𝑥̇ (𝑡) = 𝑢(𝑡)
𝑦̇ (𝑡) = 𝑣(𝑡)
𝑢̇ (𝑡) = −6𝑥 + 2𝑦
𝑣̇ (𝑡) = 2𝑥 − 6𝑦
𝑥(0) = 1, 𝑢(0) = 0, 𝑦(0) = 0, 𝑣(0) = 0 , Taking 𝛥𝑡 = 0.001

58
Scilab code 𝟔. 𝟒

//Example 6.5
clc
clear
clf
t=[0:0.001:50]'//time(s)
x(1)=1//intial postion(m) of first mass
y(1)=0//intial postion(m) of second mass
u(1)=0//intial velocity(m/s) of first mass
v(1)=0//intial velocity(m/s) of second mass
for i=1:50000//solving system of 1st order ODE
x(i+1)=x(i)+0.001*u(i)
y(i+1)=y(i)+0.001*v(i)
u(i+1)=u(i)+0.001*(-6*x(i)+2*y(i))
v(i+1)=v(i)+0.001*(2*x(i)-6*y(i))
end
subplot(2,1,1)
plot(t,x,'b')
plot(t,u,'c')
xgrid(3)
xlabel('t','Fontsize',3)
ylabel('1st mass','Fontsize',3)
title('Coupled oscillator','Fontsize',3)
legend('position','velocity',-4)
subplot(2,1,2)
plot(t,y,'r')
plot(t,v,'m')
xgrid(3)
xlabel('t','Fontsize',3)
ylabel('2nd mass','Fontsize',3)
legend('position','velocity',-4)

59
Graphic Window 𝟔. 𝟒

60
𝟔. 𝟑 Runge–Kutta method
Order 2

Let 1st order ODE = 𝑓(𝑥, 𝑦), 𝑦(𝑥 ) = 𝑦 , ℎ is the step size.

𝑦 = 𝑦 + (𝑘 + 𝑘 )

Where 𝑘 = ℎ𝑓(𝑥 , 𝑦 )
𝑘 = ℎ𝑓(𝑥 ,𝑦 + 𝑘 )

Order 4

𝑦 = 𝑦 + (𝑘 + 2𝑘 + 2𝑘 + 𝑘 )

Where 𝑘 = ℎ𝑓(𝑥 , 𝑦 )

𝑘 = ℎ𝑓 𝑥 + , 𝑦 +

𝑘 = ℎ𝑓 𝑥 + , 𝑦 +

𝑘 = ℎ𝑓(𝑥 + ℎ, 𝑦 + 𝑘 )

Scilab code 𝟔. 𝟓

//Runge–Kutta Method order 4

clc
clear
clf
function z=f(x, y)
z=x^2
endfunction
//initial condition
x(1)=0
y(1)=1
h=0.1//step size
for i=1:100
k1(i)=h*f(x(i),y(i))
k2(i)=h*f(x(i)+h/2,y(i)+k1(i)/2)
k3(i)=h*f(x(i)+h/2,y(i)+k2(i)/2)
k4(i)=h*f(x(i)+h,y(i)+k3(i))
y(i+1)=y(i)+(k1(i)+2*k2(i)+2*k3(i)+k4(i))/6

61
 x(i+1)=x(i)+h
end
plot(x,y,'b')
xgrid(3)
xlabel('x','Fontsize',4)
ylabel('y','Fontsize',4)

Graphic Window 𝟔. 𝟓

62
CHAPTER
Linear Algebra
A linear equation in 𝑛 variables 𝑥 , 𝑥 , … , 𝑥 is an equation of the form

𝑎 𝑥 + 𝑎 𝑥 + ⋯ + 𝑎 𝑥 = 𝑏.

System of Linear Equations

A system of linear equation in 𝑛 variables 𝑥 , 𝑥 , … , 𝑥 is a finite collection of

linear equations in these variables. A system of m linear equation in these n
variables can be written as
𝑎 𝑥 + 𝑎 𝑥 + ⋯+ 𝑎 𝑥 = 𝑏
𝑎 𝑥 + 𝑎 𝑥 + ⋯+ 𝑎 𝑥 = 𝑏
. . . (7.1)

. . .
𝑎 𝑥 +𝑎 𝑥 + ⋯+ 𝑎 𝑥 =𝑏
A set of values of unknowns 𝑥 , 𝑥 , … , 𝑥 which simultaneously satisfy the
system of linear equation are called solution. The set of all solutions of a
system of linear equation is called the solution set. If 𝑥 , 𝑥 , … , 𝑥 are all zero,
then the solution is said to be a trivial solution Otherwise, the solution is said
to be non-trivial.

A system of linear equations is said to be consistent if it has at least one

solution.

A system of linear equations is said to be inconsistent if it has no solution.

Given any system of linear equations precisely one of the following three is
true:

(a) The system has exactly one solution.

(b) The system has infinitely many solutions.

(c) The system has no solution.

63
The system of linear equations (7.1) may be written as
𝐴𝑋 = 𝐵
𝑎 𝑎 ⋯ 𝑎 𝑏 𝑥
𝑎 𝑎 ⋯ 𝑎 𝑏 𝑥
Where 𝐴 = ⋮ ⋮ ⋮ ⋮ ,𝐵= ,𝑋= ⋮
⋮
𝑎 𝑎 ⋯ 𝑎 𝑏 𝑥

Augmented matrix of the system of linear equations (7.1)

𝑎 𝑎 ⋯ 𝑎 𝑏
𝑎 𝑎 ⋯ 𝑎 𝑏
[𝐴|𝐵] = ⋮ ⋮ ⋮ ⋮ ⋮
𝑎 𝑎 ⋯ 𝑎 𝑏
ELEMENTARY ROW TRANSFORMATIONS

The following three operations applied on the rows of matrix are called
elementary row operations:

(a) Interchange of any two rows,𝑅 ⇔ 𝑅 .

(b) Multiplying the elements of a row by a non-zero scalar, 𝑅 → 𝑘𝑅 .

(c) Adding to the elements of a row, the constant times the corresponding
elements of another row, 𝑅 → 𝑅 + 𝑘𝑅 .

Row Echelon Form of a Matrix

A matrix is said to be in row echelon form if it has the following three

properties:

1. If a row does not consist entirely of zeros, then the first non-zero entry in
the row should be 1, called the leading 1(pivot entry).

2. If any two successive rows that do not consist entirely of zeros, then the
leading 1 in the lower row occurs farther to the right than the leading 1 in the
upper row, i.e., the leading 1 appears farther to the right as we move down
the rows of the matrix.

3. If there are any rows that consist entirely of zeros, then they are grouped
together at the bottom of the matrix.

A column containing a leading 1 is called a pivot column, while a column

having no leading 1 called a non-pivot column.

64
Transforming a Matrix to Row Echelon Form Using Elementary Row
Operations

1. Interchange the top row with another row, if necessary, to bring a non-zero
entry at the top of the column.

[Link] the entry in the first row and first column 𝑎 . divide the first row by 𝑎
in order to get a leading 1.

3. Add suitable multiple of the top row to the rows below so that all entries
below the leading 1 become zeros.

4. Now cover the top row of the matrix and start again with step 1 Applied to
the submatrix. Continue in this way until the entire matrix is in row echelon
form.
0 1 3 −4
Example 𝟕. 𝟏 Reduce the matrix 𝐴 = 2 4 8 6 to row echelon form
2 −4 0 1
by using elementary row operations.
𝑅 ⇔𝑅
2 4 8 6
0 1 3 −4
2 −4 0 1
1
𝑅 → 𝑅
2
1 2 4 3
0 1 3 −4
2 −4 0 1
𝑅 → 𝑅 − 2𝑅
1 2 4 3
0 1 3 −4
0 −8 −8 −5
𝑅 → 𝑅 + 8𝑅
1 2 4 3
0 1 3 −4
0 0 16 −37
1
𝑅 → 𝑅
16

65
 1 2 4 3
0 1 3 −4
0 0 1 −37/16
The last matrix is row echelon form.

Reduced Row Echelon Form

A matrix is said to be in reduced row echelon form if it satisfies the following

properties:

1. It is in row echelon form.

2. Each column that contains a leading 1 has zeros everywhere else in that
column.

Transforming a Matrix to Reduced Row Echelon Form Using Elementary

Row Operations

To reduce a matrix to reduced row echelon form, we first reduce it to row

echelon form, we then apply the additional step:

Beginning with last non-zero row and working upward we add suitable
multiples of each row to the rows above to get zeros above the leading 1’s.
1 2 −1 6
Example 𝟕. 𝟐 Reduce the matrix 𝐴 = 3 8 9 10 to reduced row
2 −1 2 −2
echelon form by using elementary row operations.
𝑅 → 𝑅 − 3𝑅
𝑅 → 𝑅 − 2𝑅
1 2 −1 6
0 2 12 −8
0 −5 4 −14
1
𝑅 → 𝑅
2
1 2 −1 6
0 1 6 −4
0 −5 4 −14
𝑅 → 𝑅 + 5𝑅
1 2 −1 6
0 1 6 −4
0 0 34 −34

66
 1
𝑅 → 𝑅
34
1 2 −1 6
0 1 6 −4
0 0 1 −1
𝑅 → 𝑅 +𝑅
𝑅 → 𝑅 − 6𝑅
1 2 0 5
0 1 0 2
0 0 1 −1
𝑅 → 𝑅 − 2𝑅
1 0 0 1
0 1 0 2
0 0 1 −1
The last matrix is in reduced row echelon form.

There are many methods to get solution of system of linear equation out of
which four are listed below.

(a) Gaussian Elimination Method

(b) Gauss-Jordan Elimination Method

(c) Jacobi Method

(d) Gauss-Seidel Method

(a) and (b) are the direct method which gives exact solution of linear system
𝐴𝑋 = 𝐵.
(c) and (d) are iterative method which gives approximate solution of a linear
system 𝐴𝑋 = 𝐵.

𝟕. 𝟏 Gaussian Elimination Method

Let’s consider 𝐴𝑋 = 𝐵, The steps for solving the system of linear equations by
Gaussian Elimination Method is given below:

1. Write the augmented matrix [𝐴|𝐵] of the system.

2. Use the elementary row operations to reduce the augmented matrix [𝐴|𝐵] to
a matrix [𝐶|𝐷] in row echelon form.
67
3. Write the linear system corresponding to the row echelon matrix [𝐶|𝐷] and
use back substitution to obtain the solution.

Types of Solution

(a) Inconsistent solution – If any row of the form [0 0 ⋯ 0|𝑐] occur at any
stage of the Gaussian Elimination method.

(b) Unique solution- when every column of 𝐶 is pivot column.

(c) Infinitely many solution-when non-pivot column are present in 𝐶. The

variable corresponding to non-pivot column are called independent that can
take on any real values, while those corresponding to pivot columns are called
dependent variable the values of the dependent variable are then determined
in term of independent variable by making use of back substitution.

Example 𝟕. 𝟑 Solve the following system of linear equations using

Gaussian Elimination Method.
𝑥 + 2𝑦 + 3𝑧 = 1
𝑥 + 3𝑦 + 5𝑧 = 2
2𝑥 + 5𝑦 + 9𝑧 = 3
1 2 3 1
[𝐴|𝐵 ] = 1 3 5 2
2 5 9 3
Transforming [𝐴|𝐵 ] into row echelon form
𝑅 → 𝑅 −𝑅
𝑅 → 𝑅 − 2𝑅
1 2 3 1
0 1 2 1
0 1 3 1
𝑅 → 𝑅 −𝑅
1 2 3 1
0 1 2 1
0 0 1 0
linear system corresponding to the above row echelon matrix
𝑥 + 2𝑦 + 3𝑧 = 1
𝑦 + 2𝑧 = 1

68
 𝑧=0
Getting values by back substitution

𝑧 = 0 , 𝑦 = 1, 𝑥 = −1.

Example 𝟕. 𝟒 Solve the following system of linear equations using

Gaussian Elimination Method.
3𝑥 − 6𝑥 + 3𝑥 = 9
−2𝑥 + 4𝑥 + 2𝑥 − 𝑥 = −11
4𝑥 − 8𝑥 + 6𝑥 + 7𝑥 = −5
3 −6 0 3 9
[𝐴|𝐵 ] = −2 4 2 −1 −11
4 −8 6 7 −5
Transforming [𝐴|𝐵 ] into row echelon form
1
𝑅 → 𝑅
3
1 −2 0 1 3
−2 4 2 −1 −11
4 −8 6 7 −5
𝑅 → 𝑅 + 2𝑅
𝑅 → 𝑅 − 4𝑅
1 −2 0 1 3
0 0 2 1 −5
0 0 6 3 −17
1
𝑅 → 𝑅
2
1 −2 0 1 3
0 0 1 1/2 −5/2
0 0 6 3 −17
𝑅 → 𝑅 − 6𝑅
1 −2 0 1 3
0 0 1 1/2 −5/2
0 0 0 0 −2
The last row of the above row echelon matrix is of the form [0 0 ⋯ 0|𝑐] .
So, this linear system of equations has no solution.

linear system corresponding to the above row echelon matrix

69
 𝑥 − 2𝑥 + 𝑥 = 3

𝑥 + 𝑥 =−

0 = −2

Example 7.5 Solve the following system of linear equations using Gaussian
Elimination Method.
3𝑥 + 𝑥 + 7𝑥 + 2𝑥 = 13
2𝑥 − 4𝑥 + 14𝑥 − 𝑥 = −10
5𝑥 + 11𝑥 − 7𝑥 + 8𝑥 = 59
2𝑥 + 5𝑥 − 4𝑥 − 3𝑥 = 39
3 1 7 2 13
[𝐴|𝐵 ] = 2 −4 14 −1 −10
5 11 −7 8 59
2 5 −4 −3 39
1
𝑅 → 𝑅
3
1 1/3 7/3 2/3 13/3
2 −4 14 −1 −10
5 11 −7 8 59
2 5 −4 −3 39
𝑅 → 𝑅 − 2𝑅
𝑅 → 𝑅 − 5𝑅
𝑅 → 𝑅 − 2𝑅
1 1/3 7/3 2/3 13/3
0 −14/3 28/3 −7/3 −56/3
0 28/3 −56/3 14/3 112/3
0 13/3 −26/3 −13/3 91/3
3
𝑅 →− 𝑅
14
1 1/3 7/3 2/3 13/3
0 1 −2 1/2 4
0 0 0 0 0
0 0 0 −13/2 13
𝑅 ⇔𝑅

70
 1 1/3 7/3 2/3 13/3
0 1 −2 1/2 4
0 0 0 −13/2 13
0 0 0 0 0
2
𝑅 →− 𝑅
13
1 1/3 7/3 2/3 13/3
0 1 −2 1/2 4
0 0 0 1 −2
0 0 0 0 0
linear system corresponding to the above row echelon matrix

𝑥 + 𝑥 + 𝑥 + 𝑥 =

𝑥 − 2𝑥 + 𝑥 = 4

𝑥 = −2
0=0

We can ignore the last equation, since it provides us with no information

regarding the values of the variables. The third equation gives 𝑥 = −2.
Column 3 of the matrix is a no-pivot column, so 𝑥 is the independent
variable. Therefore, we can let 𝑥 take on any value, say 𝑥 = 𝑐. Substituting
this into the second equation yields

𝑥 − 2𝑐 + (−2) = 4 ⇒ 𝑥 = 5 + 2𝑐

Then from the first equation, we obtain

𝑥 + (5 + 2𝑐) + (𝑐) + (−2) = ⇒ 𝑥 = 4 − 3𝑐

Therefore, the complete solution set is {(4 − 3𝑐, 5 + 2𝑐, 𝑐, −2): 𝑐𝜖ℝ}

Particular solutions can be found by choosing values for 𝑐. For example,

choosing 𝑐 = 1, yields

𝑥 = 1, 𝑥 = 7, 𝑥 = 1 and 𝑥 = −2.

Example 𝟕. 𝟔 Solve the following system of linear equations by Gaussian

Elimination Method using Scilab.
2𝑥 + 3𝑧 = 1
𝑦+𝑧=2

71
 𝑥+𝑦+𝑧 = 3
2 0 3 1
[𝐴|𝐵 ] = 0 1 1 2
1 1 1 3
Scilab code 7.1

//Example 7.6
clc
clear
C=[2 0 3 1;0 1 1 2;1 1 1 3]//C is augmented Matrix
disp("Augmented Matrix:-")
disp(C)
for i=1:3//Elementary Row operation
for k=(i+1):3
d=C(i,i)
t=C(k,i)
for j=1:4
C(i,j)=C(i,j)/d
C(k,j)=C(k,j)-t*C(i,j)
end
end

if (i==3) then

C(3,4)=C(3,4)/C(3,3)
C(3,3)=1
end

end
disp("Row Echelon Matrix:-")
disp(C)
for i=3:-1:1//back substitution
x(i,1)=C(i,4)
for j=(i+1):3
x(i,1)=x(i,1)-C(i,j)*x(j,1)
end
end
disp('Ans:-')
disp('x='+string(x(1,1))+'')
disp('y='+string(x(2,1))+'')
disp('z='+string(x(3,1))+'')

72
Console window 7.1

𝟕. 𝟐 Gauss-Jordan Elimination Method

Let’s consider 𝐴𝑋 = 𝐵 linear system, the steps for solving 𝐴𝑋 = 𝐵 by Gauss-
Jordan Elimination Method is given below:

1. Write the augmented matrix [𝐴|𝐵 ] of the system.

2. Use the elementary row operations to reduce the augmented matrix [𝐴|𝐵] to
a matrix [𝐶|𝐷 ] in reduced row echelon form.

3. Write the linear system corresponding to the reduced row echelon matrix
[𝐶|𝐷].

Types of Solution

(a) Inconsistent solution – If any row of the form [0 0 ⋯ 0|𝑐] occur at any
stage of the Gaussian -Jordan Elimination method.

(b) Unique solution- when every column of 𝐶 is pivot column.

(c) Infinitely many solution-when non-pivot column are present in 𝐶. The

variable corresponding to non-pivot column are called independent that can
take on any real values, while those corresponding to pivot columns are called

73
dependent variable the values of the dependent variable are then determined
in term of independent variable.

Example 7.7 Solve the following system of linear equations using the
Gauss-Jordan Method:
𝑥 + 2𝑦 + 3𝑧 = 6
2𝑥 − 3𝑦 + 2𝑧 = 14
3𝑥 + 𝑦 − 𝑧 = −2
1 2 3 6
[𝐴|𝐵] = 2 −3 2 14
3 1 −1 −2
𝑅 → 𝑅 − 2𝑅
𝑅 → 𝑅 − 3𝑅
1 2 3 6
0 −7 −4 2
0 −5 −10 −20
𝑅 ⇔𝑅
1 2 3 6
0 −5 −10 −20
0 −7 −4 2
1
𝑅 →− 𝑅
5
1 2 3 6
0 1 2 4
0 −7 −4 2
𝑅 → 𝑅 − 2𝑅
𝑅 → 𝑅 + 7𝑅
1 0 −1 −2
0 1 2 4
0 0 10 30
1
𝑅 → 𝑅
10
1 0 −1 −2
0 1 2 4
0 0 1 3
𝑅 → 𝑅 +𝑅

74
 𝑅 → 𝑅 − 2𝑅
1 0 0 1
0 1 0 −2
0 0 1 3
The last matrix is in reduced row echelon form. The linear system
corresponding to this final matrix is:

𝑥 = 1, 𝑦 = −2 and 𝑧 = 3.

Example 7.8 Solve the following system of linear equations by Gauss-

Jordan Elimination Method using Scilab
𝑥+𝑦+𝑧 = 3
2𝑥 + 3𝑦 + 7𝑧 = 0
𝑥 + 3𝑦 − 2𝑧 = 17
1 1 1 3
[𝐴|𝐵] = 2 3 7 0
1 3 −2 17
Scilab code 7.2

//Example 7.8
clc
clear
C=[1 1 1 3;2 3 7 0;1 3 -2 17]//C is augmented Matrix
disp("Augmented Matrix:-")
disp(C)
for i=1:3// for row echelon matrix
for k=(i+1):3
d=C(i,i)
t=C(k,i)
for j=1:4
C(i,j)=C(i,j)/d
C(k,j)=C(k,j)-t*C(i,j)
end
end
if (i==3) then
C(3,4)=C(3,4)/C(3,3)
C(3,3)=1
end
end
for i=3:-1:1// for reduced row echelon matrix
for k=(i-1):-1:1
75
 t=C(k,i)
for j=1:4
C(k,j)=C(k,j)-t*C(i,j)
end
end
end
disp(" Reduced Row Echelon Matrix:-")
disp(C)
disp('Ans:-')
disp('x='+string(C(1,4))+'')
disp('y='+string(C(2,4))+'')
disp('z='+string(C(3,4))+'')
Console window 7.2

76
𝟕. 𝟑 Jacobi Method
The Jacobi method is an iterative algorithm for determining the solutions of a
strictly diagonally dominant system of linear equations (𝐴 must be strictly
diagonally dominant matrix for system of linear equation 𝐴𝑋 = 𝐵).

A square matrix is said to be diagonally dominant if, for every row of the
matrix, the magnitude of the diagonal entry in a row is larger than or equal
to the sum of the magnitudes of all the other (non-diagonal) entries in that
row.

𝑎 ≤ |𝑎 |, 𝑖 = 1,2,3 … , 𝑛
,

In case of strictly diagonally dominant matrix.

𝑎 < |𝑎 |, 𝑖 = 1,2,3 … , 𝑛
,

Let 𝐴𝑋 = 𝐵 is square system of 𝑛 linear equation. 𝐴 can be decomposed such

that
𝐴=𝐷−𝐶
𝑎 𝑎 ⋯ 𝑎 𝑥 𝑏
𝑎 𝑎 ⋯ 𝑎 ⋮ ⋮
Where 𝐴 = ⋮ ⋮ ⋮ ⋮ ,𝑋= ⋮ ,𝐵 =
⋮
𝑎 𝑎 ⋯ 𝑎 𝑥 𝑏
𝑎 0 ⋯ 0 0 −𝑎 ⋯ −𝑎
0 𝑎 ⋯ 0 −𝑎 0 ⋯ −𝑎
𝐷= ,𝐶 =
⋮ ⋮ ⋱ ⋮ ⋮ ⋮ ⋱ ⋮
0 0 ⋯ 𝑎 −𝑎 −𝑎 ⋯ 0
Now, 𝐴𝑋 = 𝐵 can be Re-written as
(𝐷 − 𝐶)𝑋 = 𝐵
𝐷𝑋 = 𝐵 + 𝐶𝑋 (7.2)
Diagonal system of linear equation can be easily solved E.g.
𝑎 0 0 𝑥 𝑏
0 𝑎 0 𝑥 = 𝑏
0 0 𝑎 𝑥 𝑏

77
From above system of linear equation:

𝑥 = ,𝑥 = ,𝑥 = .

The right-hand side vector of (7.2) cannot be a known quantity because it

involves the unknown vector 𝑋. Rather, if we choose (arbitrarily) some
specific value for 𝑋, say 𝑋 = 𝑋 ( ) , then the resulting system 𝐷𝑋 = 𝐶𝑋 ( ) + 𝐵
can be readily solved. Let us call the solution of this system as 𝑋 ( ) . That is,
𝐷𝑋 ( ) = 𝐶𝑋 ( ) + 𝐵.Now taking 𝑋 = 𝑋 ( ) on the right-hand side of (7.2) we
can obtain the solution of this system, which we denote as 𝑋 ( ) and repeat this
procedure to get a general iterative procedure as

𝐷𝑋 ( )
= 𝐶𝑋 ( )
+ 𝐵, 𝑘 = 0,1,2,··· (7.3)

If 𝐷 is invertible, then the above iterative procedure can be written as

𝑋 ( ) = 𝑏𝑋 ( ) + 𝑐, 𝑘 = 0,1,2,··· (7.4)
where 𝑏 = 𝐷 𝐶 and 𝑐 = 𝐷 𝐵.

(7.3) can be Re-written as

( )
1 ( )
𝑥 = (𝑏 − 𝑎 𝑥 ) , 𝑖 = 1,2,3, …
𝑎

The iterative procedure (7.4) is called the Jacobi method.

The convergence of (7.3) to the solution of 𝐴𝑋 = 𝐵, only guaranteed when 𝐴 is

strictly diagonally dominant matrix.

Example 𝟕. 𝟗 Solve the following system of linear equations by the Jacobi

method.
2𝑥 + 𝑥 = 11
5𝑥 + 7𝑥 = 13

2 1 11 𝑥
𝐴= ,𝐵 = ,𝑋 = 𝑥
5 7 13
𝐴 is strictly diagonally dominant.

If 𝐴 is decomposed into 𝐷 and −𝐶, 𝐴 = 𝐷 − 𝐶

2 0 0 −1
Then 𝐷 = ,𝐶 =
0 7 −5 0

78
(𝐷 − 𝐶)𝑋 = 𝐵
𝐷𝑋 = 𝐵 + 𝐶𝑋
2 0 𝑥 11 0 −1 𝑥
= +
0 7 𝑥 13 −5 0 𝑥
2𝑥 = 11 − 𝑥
7𝑥 = 13 − 5𝑥
( ) ( )
𝑥 = − 𝑥
( ) ( )
𝑥 = − 𝑥 ,𝑘 = 0,1,2, …

Choosing initial guess 𝑋 ( ) = (1,1)

𝑋 ( ) = (5, 1.1428571)

𝑋 ( ) = (4.9285714, -1.7142857)

𝑋 ( ) = (6.3571429, -1.6632653)

𝑋 ( ) = (6.3316327, -2.68367350)

𝑋 ( ) = (6.8418367, -2.6654519)

𝑋 ( ) = (6.8327259, -3.0298834)

𝑋 ( ) = (7.0149417, -3.0233757)

𝑋 ( ) = (7.0116878, -3.1535298)

Solving system of linear equations by Jacobi method using Scilab.

Scilab Code 𝟕. 𝟑

//Jacobi Method
clc
clear
A=input(" Enter the coefficient matrix of nxn:-")
B=input(" Enter the right-hand side matrix nx1:-")
X0=input(" Initial approximation nx1:-")
tol=input(" Enter the convergence tolerance ")//Infinity norm of (A*X1-B)
t1=0
t2=0
disp(" Coefficient Matrix:-")
disp(A)
disp(" Right hand side Matrix:-")
disp(B)
79
[m, n] = size( A )
if ( m ~= n )
disp( " The input matrix is not square")
t1=1
end
for i=1:m//To check if the matrix is diagonally dominant
s=0
for j=1:m
s=s+abs(A(i,j))
end
if 2*abs(A(i,i))<s then
disp(" the matrix is not strictly diagonally dominant")
t2=1
break
end
end
if (t1==1|t2==1) then
disp(" Execute again and enter square and strictly diagonally dominant
coeffiecient matrix ")
else
for i=1:m
for j=1:m
if(i==j) then
D(i,j)=A(i,j)
else
D(i,j)=0
end
end
end
disp(" Diagonal matrix of coefficient matrix:-")
disp(D)
C=D-A
//Implements the Jacobian algorithm
Error=1
while tol<Error
X1=(inv(D)*C*X0)+(inv(D)*B)//Computes the jacobian updates
Error=norm((A*X1)-B,'inf')
X0=X1
end
disp(" Solution vector is:-")
disp(X1)
end

80
Console window 𝟕. 𝟑

𝟕. 𝟒 Gauss-Seidel Method
The Gauss-Seidel Method is an iterative algorithm for determining the
solutions of a strictly diagonally dominant or positive definite system of linear
equations. (𝐴 must be strictly diagonally dominant or positive definite matrix
for system of linear equation 𝐴𝑋 = 𝐵). Let 𝐴𝑋 = 𝐵 is square system of 𝑛 linear

A symmetric 𝑛 × 𝑛, real Matrix 𝐴 is said to be positive definite if the scalar

𝑧 𝐴𝑧 is strictly positive for every non-zero column vector 𝑧 of 𝑛 real
numbers. Here 𝑧 denotes the transpose of 𝑧.

equation. 𝐴 can be decomposed such that

81
 𝐴=𝐿+𝑈
Where 𝐿 is lower triangular matrix and 𝑈 is strictly upper triangular matrix.
𝑎 𝑎 ⋯ 𝑎 𝑥 𝑏
𝑎 𝑎 ⋯ 𝑎 ⋮ ⋮
𝐴= ⋮ ⋮ ⋮ ⋮ ,𝑋= ⋮ ,𝐵 =
⋮
𝑎 𝑎 ⋯ 𝑎 𝑥 𝑏
𝑎 0 ⋯ 0 0 𝑎 ⋯ 𝑎
𝑎 𝑎 ⋯ 0 0 0 ⋯ 𝑎
𝐿= ,𝑈=
⋮ ⋮ ⋱ ⋮ ⋮ ⋮ ⋱ ⋮
𝑎 𝑎 ⋯ 𝑎 0 0 ⋯ 0
Now, 𝐴𝑋 = 𝐵 can be Re-written as
(𝐿 + 𝑈)𝑋 = 𝐵
𝐿𝑋 = 𝐵 − 𝑈𝑋
if we choose (arbitrarily) some specific value for 𝑋, say 𝑋 = 𝑋 ( ) . Like, we
have done in Jacobi method then the resulting system 𝐿𝑋 = 𝐵 − 𝑈𝑋 can be
solved using forward substitution for 𝑋 ( ) .

𝑋 ( ) = 𝐿 (𝐵 − 𝑈𝑋 ( ) ), 𝑘 = 0,1,2,··· (7.5)
The procedure (7.5) is generally continued until the changes made by an
iteration are below some tolerance. The iterative procedure (7. 5) is called the
Gauss-Seidel Method.

Solving system of linear equations by Gauss-Seidel Method using Scilab.

Scilab code 𝟕. 𝟒

//Gauss-Seidel Method
clc
clear
A=input(" Enter the coefficient matrix of nxn:-")
B=input(" Enter the right-hand side matrix nx1:-")
X0=input(" Initial approximation nx1:-")
tol=input(" Enter the convergence tolerance ")//tol>Infinity norm of (A*X1-B)
t1=0
t2=0
disp(" Coefficient Matrix:-")
disp(A)
disp(" Right hand side Matrix:-")
disp(B)
[m, n] = size( A )

82
 if ( m ~= n )
disp( " The input matrix is not square")
t1=1
end
for i=1:m//To check if the matrix is diagonally dominant
s=0
for j=1:m
s=s+abs(A(i,j))
end
if 2*abs(A(i,i))<s then
disp(" the matrix is not strictly diagonally dominant")
t2=1
break
end
end
if (t1==1|t2==1) then
disp(" Execute again and enter square and strictly diagonally dominant
coefficient matrix ")
else
for i=1:m
for j=1:m
if(i<j) then
U(i,j)=A(i,j)
else
U(i,j)=0
end
end
end
disp(" strictly upper triangular matrix of coefficient matrix:-")
disp(U)
L=A-U
disp(" lower triangular matrix of coefficient matrix:-")
disp(L)
//Implements the Gauss-Seidel algorithm
Error=1
while tol<Error
X1=inv(L)*(B-(U*X0))//Computes the Gauss-Seidel updates
Error=norm((A*X1)-B,'inf')
X0=X1
end
disp(" Solution vector is:-")
disp(X1)
end

83
Console window 7.4

84
CHAPTER
Finite Difference Method
Chapter-6 was about how to solve ODE numerically. This Chapter deals with
solving PDEs numerically. There are many methods for solving PDEs
numerically out of which only finite difference method is discussed in this
chapter.

PDE (Partial Differential Equation) is an equation that contains at least two

independent variables, a dependent variable or the unknown function and
one or more partial derivative of the unknown function.

In this chapter, a subscript notation is used for PDE such that

𝑢 = ,𝑢 = ,𝑢 = ,𝑢 = …

𝟖. 𝟏 Classification of PDE
General terminology of PDE

The order of a PDE is the order of the highest partial derivative which appears
in the equation.

E.g. 𝑢 = 𝑐 𝑢 is a second-order PDE.

The most general first-order PDE in two independent variables 𝑥 and 𝑦 is of

the form

𝐹 𝑥, 𝑦, 𝑢, 𝑢 , 𝑢 =0

Linearity of PDE

A PDE is said to be linear if it is of the first degree in the unknown function 𝑢

and its partial derivatives. Otherwise, it is called non-linear.

E.g. 𝑥𝑢 + 𝑦𝑢 = 𝑢 is linear and 𝑥 𝑢 +𝑦 𝑢 = 𝑢 is non-linear PDE.

The most general first-order linear partial differential equation in two

independent variables 𝑥 and 𝑦 has the form
𝑎(𝑥, 𝑦)𝑢 + 𝑏(𝑥, 𝑦)𝑢 + 𝑐(𝑥, 𝑦)𝑢 = 𝑑(𝑥, 𝑦)

Where,𝑎, 𝑏, 𝑐 and 𝑑 in general are function of 𝑥 and 𝑦.

85
Similarly, A second-order linear partial differential equation in two
independent variables 𝑥 and 𝑦 has the form
𝐴𝑢 + 𝐵𝑢 + 𝐶𝑢 + 𝐷𝑢 + 𝐸𝑢 + 𝐹𝑢 = 𝐺

Where 𝐴,𝐵,𝐶,𝐷,𝐸,𝐹,𝐺 are the function of 𝑥 and 𝑦.

Classification of second-order linear PDE

𝐴𝑢 + 𝐵𝑢 + 𝐶𝑢 + 𝐷𝑢 + 𝐸𝑢 + 𝐹𝑢 = 𝐺 is said to be hyperbolic, parabolic

or elliptic at point (𝑥 , 𝑦 ) according as the discriminant

𝛥(𝑥 , 𝑦 ) = 𝐵 (𝑥 , 𝑦 ) − 4𝐴(𝑥 , 𝑦 ) 𝐶(𝑥 , 𝑦 ) is positive, zero or negative

respectively. If this is true at all point then the equation is said to be
hyperbolic, parabolic or elliptic in a domain.

E.g. 1. Tricomi Equation

𝑢 + 𝑥𝑢 =0

𝐴 = 1, 𝐵 = 0, 𝐶 = 𝑥
𝛥 > 0 for 𝑥 < 0 and 𝛥 < 0 for 𝑥 > 0. So, the Tricomi equation is hyperbolic for
𝑥 < 0 and Elliptic for 𝑥 > 0.

2. 2D Laplace's Equation
𝑢 +𝑢 =0

𝐴 = 1, 𝐵 = 0, 𝐶 = 1

𝛥 < 0 for entire ℝ . So, Laplace's Equation is elliptic for entire ℝ .

3. 1D Heat Equation
𝑢 =𝛼 𝑢
𝐴 = 𝛼 , 𝐵 = 0, 𝐶 = 0
𝛥 = 0 for entire ℝ . So, 1D Heat Equation is parabolic for entire ℝ .

4. 1D Wave Equation
𝑢 =𝑐 𝑢
𝐴 = 𝑐 , 𝐵 = 0, 𝐶 = −1
𝛥 > 0 for entire ℝ . So, 1D Wave Equation is hyperbolic for entire ℝ .

86
Homogeneous and non-homogeneous linear PDE

A linear PDE is said to be homogeneous if each of its term contains either 𝑢 or

one of its partial derivatives. Otherwise, it is called non-homogeneous.

E.g. 𝑢 = 𝑐 (𝑢 + 𝑢 ) is homogeneous and 𝑢 +𝑢 = 𝑓(𝑥, 𝑦) is non-

homogeneous.

A function of independent variables of PDE that satisfy PDE in entire

space of independent variables is called Solution of PDE.

𝟖. 𝟐 Finite Difference Fundamentals

The finite difference method (FDM) works by replacing the region over which
the independent variables in the PDE are defined by a finite grid (produced
by discretization of independent variables region) of points at which the
dependent variable is approximated. The partial derivatives in the PDE at
each grid point are approximated from neighboring values by using Taylor’s
theorem.

Taylor’s Theorem

Let’s 𝑢(𝑥)have 𝑛 continuous derivative’s over the interval (𝑎, 𝑏). Then for

𝑎 < 𝑥 , 𝑥 + ℎ < 𝑏,
( ) ( )
𝑢(𝑥 + ℎ) = 𝑢(𝑥 ) + ℎ𝑢 (𝑥 ) + +⋯ ( )!
+ 𝑜(ℎ ) (8.1)
!

Where,

𝑢 = ,𝑢 = ,…,𝑢 =

The usual interpretation of Taylor’s theorem says that if we know the value of
𝑢 and the values of its derivatives at point 𝑥 then we can write down the
equation (8.1) for its value at the (nearby) point 𝑥 + ℎ. This expression
contains an unknown quantity(error) which is written in as 𝑜(ℎ ) and
pronounced ‘order ℎ to the 𝑛’. If we discard the term 𝑜(ℎ ) in (8.1) (𝑖. 𝑒.
truncate the right-hand side of (8.1)) we get an approximation to 𝑢(𝑥 + ℎ).
The error in this approximation is 𝑜(ℎ ).

87
 Let 𝑓: [𝑎, 𝑏] → ℝ, 𝑓 , 𝑓 , … , 𝑓 ( ) be continuous on [𝑎, 𝑏] and suppose 𝑓
exist on (𝑎, 𝑏). Then there exist 𝑐 ∈ (𝑎, 𝑏) such that

𝑓( ) (𝑎)(𝑏
− 𝑎) 𝒇𝒏 (𝒄)(𝒃 − 𝒂)𝒏
𝑓(𝑏) = 𝑓(𝑎) + 𝑓 (𝑎)(𝑏 − 𝑎) + ⋯ + +
(𝑛 − 1)! 𝒏!

In the (8.1) both 𝑥 and 𝑥 + ℎ are grid points for FDM. This allows us to
rearrange equation (8.1) to get so-called Finite Difference (FD)
approximations to derivatives which have 𝑜(ℎ ) errors.

Let’s see simple finite difference approximation to derivative, Truncating (8.1)

after the first derivative
𝑢(𝑥 + ℎ) = 𝑢(𝑥 ) + ℎ𝑢 (𝑥 ) + 𝑜(ℎ )
𝑢(𝑥 + ℎ) − 𝑢(𝑥 ) 𝑜(ℎ )
𝑢 (𝑥 ) = +
ℎ ℎ
𝑢(𝑥 + ℎ) − 𝑢(𝑥 )
= + 𝑜(ℎ)
ℎ
Neglecting the error term 𝑜(ℎ).
𝑢(𝑥 + ℎ) − 𝑢(𝑥 )
𝑢 (𝑥 ) ≈ (8.2)
ℎ
(8.2) is called a first order finite difference approximation to 𝑢 (𝑥 ) since the
approximation error = 𝑜(ℎ)which depends on the first power of ℎ. This
approximation is called a forward finite difference approximation since we
start at 𝑥 and step forwards to the point 𝑥 + ℎ. ℎ is called the step size (ℎ >
0).

The following simple example shows that our forward difference

approximation (8.2) works and has the stated order of accuracy. We choose a
simple function for 𝑢. Let 𝑢(𝑥) = 𝑥 . We will find the first order forward FD
approximation to 𝑢 (3) using step 𝑠𝑖𝑧𝑒 ℎ = 0.1

Using (8.2)
(3 + 0.1) − (3)
𝑢 (3) ≈ = 6.1
0.1

88
The exact answer from basic Calculus is clearly 𝑢 (3) = 6 so the error in the
approximation is 6.1 – 6 = 0.1. Repeating the problem with ℎ = 0.05 (𝑖. 𝑒.
half the step size) gives,

(3 + 0.05) − (3)
𝑢 (3) ≈ = 6.05
0.05
The error is 6.05 – 6 = 0.05. The approximation formula (8.2) is first order so
the errors should be proportional to ℎ which is seen to be the case: halving the
step size results in a halving of the error.

Constructing a Finite Difference Toolkit

We now construct common FD approximations to common partial

derivatives. For simplicity, we suppose that 𝑢 is a function of only two
variables, 𝑥 and 𝑡. We will approximate the partial derivatives of 𝑢 with
respect to 𝑥. As 𝑡 is held constant 𝑢 is effectively a function of the single
variable 𝑥 so we can use Taylor’s formula (8.1) where the ordinary derivative
terms are now partial derivatives and the arguments are (𝑥, 𝑡) instead of 𝑥.
Finally, we will replace the step size ℎ by 𝛥𝑥 so that (8.1) becomes,
( , )
𝑢(𝑥 + 𝛥𝑥, 𝑡) = 𝑢(𝑥 , 𝑡) + 𝛥𝑥𝑢 (𝑥 , 𝑡) + ⋯ ( )!
+ 𝑜(𝛥𝑥 ) (8.3)

Truncating (8.3) to 𝑜(𝛥𝑥 ) gives,

𝑢(𝑥 + 𝛥𝑥, 𝑡) = 𝑢(𝑥 , 𝑡) + 𝛥𝑥𝑢 (𝑥 , 𝑡) + 𝑜(𝛥𝑥 ) (8.4)
Now we derive some FD approximations to partial derivatives. Rearranging
(8.4) gives,
𝑢(𝑥 + 𝛥𝑥, 𝑡) − 𝑢(𝑥 , 𝑡) 𝑜(𝛥𝑥 )
𝑢 (𝑥 , 𝑡) = +
𝛥𝑥 𝛥𝑥
𝑢(𝑥 + 𝛥𝑥, 𝑡) − 𝑢(𝑥 , 𝑡)
= + 𝑜(𝛥𝑥) (8.5)
𝛥𝑥
Equation (8.5) holds at any point (𝑥 , 𝑡). In numerical schemes for solving
PDEs, we are restricted to a grid of discrete 𝑥 values, 𝑥 , 𝑥 , … , 𝑥 and discrete
𝑡 levels 𝑡 , 𝑡 , … We will assume a constant grid spacing, 𝛥𝑥 in 𝑥 so that 𝑥 =
𝑥 + 𝛥𝑥 Evaluating Equation (8.5) for a point, 𝑥 , 𝑡 on the grid gives,

𝑢 𝑥 ,𝑡 − 𝑢 𝑥 ,𝑡
𝑢 𝑥 ,𝑡 = + 𝑜(𝛥𝑥) (8.6)
𝛥𝑥

89
We will use the common subscript/superscript notation,

𝑢 = 𝑢 𝑥 ,𝑡

So that dropping 𝑜(𝛥𝑥) from (8.6) becomes

𝑢 −𝑢
𝑢 𝑥 ,𝑡 = (8.7)
𝛥𝑥
(8.7) is the first order forward difference approximation to 𝑢 𝑥 , 𝑡 that we
derive previously approximation (8.2). We now derive another FD
approximation to 𝑢 𝑥 , 𝑡 . Replacing 𝛥𝑥 by −𝛥𝑥 in (8.4) gives,

𝑢(𝑥 − 𝛥𝑥, 𝑡) = 𝑢(𝑥 , 𝑡) − 𝛥𝑥𝑢 (𝑥 , 𝑡) + 𝑜(𝛥𝑥 ) (8.8)

Evaluating (8.8) at 𝑥 , 𝑡 and rearranging as previously gives,

𝑢 −𝑢
𝑢 𝑥 ,𝑡 = (8.9)
𝛥𝑥
(8.9) is the first order backward difference approximation to 𝑢 𝑥 , 𝑡 .

Our first two FD approximations are first order in 𝑥 but we can increase the
order (and so make the approximation more accurate) by taking more terms
in the Taylor series as follows. Truncating (8.3) to 𝑜(𝛥𝑥 ),then replacing 𝛥𝑥 by
−𝛥𝑥 and subtracting this new expression from (8.3) and evaluating at 𝑥 , 𝑡 .
gives, after some algebra,

𝑢 −𝑢
𝑢 𝑥 ,𝑡 = (8.10)
2𝛥𝑥
(8.10) is called the second order central difference FD approximation to
𝑢 𝑥 ,𝑡 .

We could construct even higher order FD approximations to 𝑢 by taking

even more terms in the Taylor series but we will stop at second order
approximations to first order derivatives.

Many PDEs of interest contain second order (and higher) partial derivatives
so we need to derive approximations to them. We will restrict our attention to
second order unmixed partial derivatives 𝑖. 𝑒. 𝑢 .

Truncating (8.3) to 𝑜(𝛥𝑥 ) gives,

90
𝑢(𝑥 + 𝛥𝑥, 𝑡)
𝛥𝑥 𝑢 (𝑥 , 𝑡) 𝛥𝑥 𝑢 (𝑥 , 𝑡)
= 𝑢(𝑥 , 𝑡) + 𝛥𝑥𝑢 (𝑥 , 𝑡) + +
2! 3!
+ 𝑜(𝛥𝑥 ) (8.11)
Replacing 𝛥𝑥 by − 𝛥𝑥 in (8.11)
𝑢(𝑥 − 𝛥𝑥, 𝑡)
𝛥𝑥 𝑢 (𝑥 , 𝑡) 𝛥𝑥 𝑢 (𝑥 , 𝑡)
= 𝑢(𝑥 , 𝑡) − 𝛥𝑥𝑢 (𝑥 , 𝑡) + −
2! 3!
+ 𝑜(𝛥𝑥 ) (8.12)
Adding (8.11) and (8.12) gives
𝑢(𝑥 + 𝛥𝑥, 𝑡) + 𝑢(𝑥 − 𝛥𝑥, 𝑡) = 2𝑢(𝑥 , 𝑡) + 𝛥𝑥 𝑢 (𝑥 , 𝑡) + 𝑜(𝛥𝑥 ) (8.13)
Evaluating (8.13) at 𝑥 , 𝑡 and using our discrete notation gives,

𝑢 +𝑢 = 2𝑢 + 𝛥𝑥 𝑢 𝑥 ,𝑡 + 𝑜(𝛥𝑥 ) (8.14)

Rearranging (8.14) and dropping the 𝑜(𝛥𝑥 ) error term gives,

𝑢 − 2𝑢 + 𝑢
𝑢 𝑥 ,𝑡 = (8.15)
𝛥𝑥
(8.15) is the second order symmetric difference FD approximation to
𝑢 𝑥 , 𝑡 . These results are put into Table 8.1 to form a FD approximation
toolkit. FD approximations to partial derivatives with respect to 𝑡 are derived
in a similar manner and are included in Table 8.1.

91
partial derivative finite difference type order
approximation
𝜕𝑢 𝑢 −𝑢 forward first in 𝑥
= 𝑢 𝑥 ,𝑡
𝜕𝑥 𝛥𝑥
𝜕𝑢 𝑢 −𝑢 backward first in 𝑥
= 𝑢 𝑥 ,𝑡
𝜕𝑥 𝛥𝑥
𝜕𝑢 𝑢 −𝑢 central second in 𝑥
= 𝑢 𝑥 ,𝑡
𝜕𝑥 2𝛥𝑥
𝜕 𝑢 𝑢 − 2𝑢 + 𝑢 symmetric second in 𝑥
=𝑢 𝑥 ,𝑡
𝜕𝑥 𝛥𝑥
𝜕𝑢 𝑢 −𝑢 forward first in 𝑡
= 𝑢 𝑥 ,𝑡
𝜕𝑡 𝛥𝑡
𝜕𝑢 𝑢 −𝑢 backward first in 𝑡
= 𝑢 𝑥 ,𝑡
𝜕𝑡 𝛥𝑡
𝜕𝑢 𝑢 −𝑢 central second in 𝑡
= 𝑢 𝑥 ,𝑡
𝜕𝑡 2𝛥𝑡
𝜕 𝑢 𝑢 − 2𝑢 + 𝑢 symmetric Second in 𝑡
=𝑢 𝑥 ,𝑡
𝜕𝑡 𝛥𝑡
Table 8.1

𝟖. 𝟑 2D Laplace's Equation
𝜕 𝑢 𝜕 𝑢
+ =0 (8.16)
𝜕𝑥 𝜕𝑦
Let’s take computational domain rectangular. The computational domain is
discretized using constant grid 𝛥𝑥 and 𝛥𝑦 in 𝑥 and 𝑦 direction such that
𝑥 , 𝑥 , … , 𝑥 and 𝑦 , 𝑦 , … , 𝑦 . 𝑢(𝑥 , 𝑦 ) is denoted by 𝑢 , .

The boundary condition is given for (8.16). So, (𝑀 − 2) × (𝑁 − 2) grid points

are unknown out of 𝑀 × 𝑁 grid points.

Constructing a Finite difference scheme for (8.16) using Table 8.1

𝑢 , − 2𝑢 , + 𝑢 , 𝑢, − 2𝑢 , + 𝑢 ,
+ =0 (8.17)
𝛥𝑥 𝛥𝑦
Taking 𝛥𝑥 = 𝛥𝑦 (8.17) becomes

92
 𝑢 , +𝑢 , +𝑢, +𝑢,
𝑢, = (8.18)
4
(8.18) is our FDS for Laplace’s equation.

(8.18) can be solved using Jacobi iteration by taking some initial guess for
interior grid points such that
𝑢 , +𝑢 , +𝑢, +𝑢,
𝑢, = (8.19)
4
Superscript 𝑚 is iteration index, iterates (8.19) till |𝑢 − 𝑢 | < 𝑡𝑜𝑙.

𝑡𝑜𝑙 is pre-defined tolerance (allowed error).

|𝐴| is infinity norm of matrix 𝐴. It is equal to maximum of absolute sum

of row. 𝑖. 𝑒|𝐴| = 𝑚𝑎𝑥 𝑎 for all 𝑖.

Let’s understand (8.19) for a simple example.

The rectangular computational domain is given below

8.9 8.9 8.9 8.9 8.9
8.4 𝑢 , 𝑢 , 𝑢 , 9.2
7.2 𝑢 , 𝑢 , 𝑢 , 9.4
6.1 6.8 7.7 8.7 9.8

Setting all interior grid points to zero as the initial guess

𝑚=0
𝑢 , +𝑢 , +𝑢 , +𝑢 ,
𝑢 , =
4
7.2 + 0 + 6.8 + 0
𝑢 , = = 3.5
4
𝑢 , +𝑢 , +𝑢 , +𝑢 ,
𝑢 , =
4
0 + 0 + 7.7 + 0
𝑢 , = = 1.925
4
𝑢 , +𝑢 , +𝑢 , +𝑢 ,
𝑢 , =
4
0 + 9.4 + 8.7 + 0
𝑢 , = = 4.525
4

93
 𝑢 , +𝑢 , +𝑢 , +𝑢 ,
𝑢 , =
4
8.4 + 0 + 0 + 8.9
𝑢 , = = 4.325
4
𝑢 , +𝑢 , +𝑢 , +𝑢 ,
𝑢 , =
4
0 + 0 + 0 + 8.9
𝑢 , = = 2.225
4
𝑢 , +𝑢 , +𝑢 , +𝑢 ,
𝑢 , =
4
0 + 9.2 + 0 + 8.9
𝑢 , = = 4.525
4
Do the same for 𝑚 = 1,2,3, … 𝑀 until |𝑢 − 𝑢 | < 𝑡𝑜𝑙

The solution is 𝑢 .

Example 𝟖. 𝟏 Find the potential in square plate 0 ≤ 𝑥 ≤ 2, 0 ≤ 𝑦 ≤ 2 using

Scilab if the upper side is kept at potential 1000 𝑠𝑖𝑛 𝑥 and the other sides
are grounded.

No charge is lying on the square plate so,

+ = 0, 𝑉 is potential

Scilab code 8.1

//Laplace's equation
clc
clear
clf
dx=0.1//step size in x direction
dy=0.1//step size in y direction
nx=21
ny=21
x=[0:0.1:2]
y=[0:0.1:2]
tol=0.00000000000001
//matrix u store grid points values such that at u(i,j)=u(xi,yj)
//applying boundary condition and taking intial guess zero for interior grid
points
for i=1:nx

94
 for j=1:ny
if j==ny then
u(i,j)=1000*sin(%pi*0.5*(i-1)*dx)

else
u(i,j)=0

end
end
end
err=1//Initialization of while loop
uk=u
while err>tol//jacobi iteration
for i=2:nx-1//both for do the jacobian updates to entire grid points except for
boundary
for j=2:ny-1
uk(i,j)=(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1))*0.25

end
end
err=norm( uk - u, 'inf' )//infinity norm of |uk-u|
u=uk

end
//plotting potential of square plate
surf(x,y,uk','color_flag',0)
colorbar(0,1000)
xlabel('x','Fontsize',4)
ylabel('y','Fontsize',4)
zlabel('V','Fontsize',4)
title('POTENTIAL IN SQUARE PLATE','Fontsize',4)

95
Graphic window 8.1

𝟖. 𝟒 1D Heat Equation
Consider the following initial-boundary value problem for the heat equation
𝜕𝑢 𝜕 𝑢
=𝛼 ,0 < 𝑥 < 𝐿 (8.20)
𝜕𝑡 𝜕𝑥
BC: 𝑢(0, 𝑡) = 0 𝑢(𝐿, 𝑡) = 0

IC: 𝑢(𝑥, 0) = 𝑓(𝑥)

Discretizing space and time such that (𝑥 , 𝑥 , … , 𝑥 ), (𝑡 , 𝑡 , … , 𝑡 ) using

constant grid 𝛥𝑥 and 𝛥𝑡.

Using Table 8.1

𝑢 −𝑢
𝑢 (𝑥 , 𝑡 ) = (8.21)
𝛥𝑡
𝑢 − 2𝑢 + 𝑢
𝑢 (𝑥 , 𝑡 ) = (8.22)
𝛥𝑥
Substituting (8.21) and (8.22) in (8.20)

𝑢 −𝑢 𝑢 − 2𝑢 + 𝑢
=𝛼
𝛥𝑡 𝛥𝑥
96
 𝛼 𝛥𝑡
𝑢 =𝑢 + (𝑢 − 2𝑢 + 𝑢 ) (8.23)
𝛥𝑥
(8.23) is our Finite difference scheme for the 1D heat equation.

Stability of finite difference scheme for the heat equation

A finite difference scheme is stable if the errors made at one time step of the
calculation do not cause the errors to be magnified as the computations are
continued. A neutrally stable scheme is one in which errors remain constant
as the computations are carried forward. If the errors decay and eventually
damp out, the numerical scheme is said to be stable. If, on the contrary, the
errors grow with time the numerical scheme is said to be unstable.

Let 𝑢 = 𝜙 𝑒 then

Substituting 𝜙 𝑒 into (8.23)

𝛼 𝛥𝑡
𝜙 𝑒 =𝜙 𝑒 + 𝑒 −2+𝑒 𝜙 𝑒
𝛥𝑥
𝛼 𝛥𝑡
𝜙 𝑒 + (2 cos(𝛥𝑥𝜃) − 2)𝜙 𝑒
𝛥𝑥
𝛼 𝛥𝑡
𝜙 =𝜙 − 4 sin (𝛥𝑥𝜃 ⁄2)𝜙
𝛥𝑥
𝛼 𝛥𝑡
𝜙 = (1 − 4 sin (𝛥𝑥𝜃 ⁄2))𝜙
𝛥𝑥
Now for stability, we require that |𝜙 | ≤ |𝜙 | so that
𝛼 𝛥𝑡
|1 − 4 sin (𝛥𝑥𝜃 ⁄2) | ≤ 1
𝛥𝑥
𝛼 𝛥𝑡
−2 ≤ − 4 sin (𝛥𝑥𝜃 ⁄2) ≤ 0
𝛥𝑥
The right inequality is satisfied automatically, while the left inequality can be
re-written in the form:
𝛼 𝛥𝑡
4 sin (𝛥𝑥𝜃 ⁄2) ≤ 2
𝛥𝑥
Since sin (𝛥𝑥𝜃 ⁄2) ≤ 1this condition is satisfied for all 𝜃 provided
𝛥𝑥
𝛥𝑡 ≤
2𝛼

97
Example 8.2 Consider a laterally insulated metal bar of length 1 and such
that 𝛼 = 1 in the heat equation. Suppose that the ends of the bar are kept at
temperature 0°C and initially the temperature along the bar is given by
𝑠𝑖𝑛(𝜋𝑥) find the temperature in 0 ≤ 𝑡 ≤ 0.2 using Scilab.

Using 𝑢 =𝑢 + (𝑢 − 2𝑢 + 𝑢 )

𝛼 = 1 taking 𝛥𝑡 = 0.004 and 𝛥𝑥 = 0.1

= 0.8 ≤ 1 so (8.23) is stable

Using initial and boundary condition

𝑢(𝑥 , 𝑡 ) = 𝑢(𝑥 , 𝑡 ) = 0 for all time level

𝑢(𝑥 , 𝑡 ) = 𝑠𝑖𝑛(𝜋(𝑛 − 1)𝛥𝑥) along the entire rod.

𝑢 = 𝑢 + 0.4(𝑢 − 2𝑢 + 𝑢 )
𝑘=1
𝑢 = 𝑢 + 0.4(𝑢 − 2𝑢 + 𝑢 )
𝑢 = 𝑢 + 0.4(𝑢 − 2𝑢 + 𝑢 )
⋅
⋅
⋅
𝑢 =𝑢 + 0.4(𝑢 − 2𝑢 +𝑢 )
Similarly, for 𝑘 = 2,3, … ,50

Scilab code 8.2

//solving 1D heat equation

clc
clear
clf
x=[0:0.1:1]'
//u(n,k)
//u(space,time)
dx=0.1//step size in x
dt=0.004//step size in t
nx=11//number of space grid poits
nt=51//number of time grid points

98
//using initial and boundary condition
u(1,:)=0
u(11,:)=0
u(:,1)=sin(%pi*(0:10)*dx)
//using FDS
for n=2:10
u(n,2)=u(n,1)+0.4*(u(n+1,1)-2*u(n,1)+u(n-1,1))
end
for k=2:50
for n=2:10
u(n,k+1)=u(n,k)+0.4*(u(n+1,k)-2*u(n,k)+u(n-1,k))
end
end
plot(x,u(:,1),'m-*-')
plot(x,u(:,11),'r-+-')
plot(x,u(:,21),'b-^-')
plot(x,u(:,31),'c-')
plot(x,u(:,41),'r.')
plot(x,u(:,51),'r^')
xgrid(4)
xlabel('x','Fontsize',3)
ylabel('u(x,t)','Fontsize',3)
title('Temperature along rod for different time level','Fontsize',3)
legend('t=0 ','t=0.04 ','t=0.08 ','t=0.12 ','t=0.16 ','t=0.2 ',-4)

Graphic window 8.2

99
𝟖. 𝟓 1D Wave Equation
Consider the following initial boundary value problem for the Wave
Equation:
𝜕 𝑢 𝜕 𝑢
=𝑐 ,0 < 𝑥 < 𝐿 (8.24)
𝜕𝑡 𝜕𝑥
BC: 𝑢(0, 𝑡) = 0 𝑢(𝐿, 𝑡) = 0

IC: 𝑢(𝑥, 0) = 𝑓(𝑥)

𝜕𝑢
(𝑥, 0) = 𝑔(𝑥) (8.25)
𝜕𝑡
Discretizing space and time such that (𝑥 , 𝑥 , … , 𝑥 ), (𝑡 , 𝑡 , … , 𝑡 ) using
constant grid 𝛥𝑥 and 𝛥𝑡.

𝑢 − 2𝑢 + 𝑢
𝑢 (𝑥 , 𝑡 ) = (8.26)
𝛥𝑡
𝑢 − 2𝑢 + 𝑢
𝑢 (𝑥 , 𝑡 ) = (8.27)
𝛥𝑥
Substituting (8.25) and (8.26) in (8.24)

𝑢 − 2𝑢 + 𝑢 𝑢 − 2𝑢 + 𝑢
=𝑐
𝛥𝑡 𝛥𝑥
𝑢 = 2𝑢 − 𝑢 + (𝑐𝛥𝑡 ⁄𝛥𝑥) (𝑢 − 2𝑢 + 𝑢 )
𝑢 =𝑟 𝑢 + 2(1 − 𝑟 )𝑢 + 𝑟 𝑢 −𝑢 (8.28)

(8.28) is FDS for 1D Wave equation where 𝑟 = (𝑐∆𝑡/∆𝑥) is known as the

Courant Number. We observe that the Discrete Equation (8.28) involves three
distinct levels in which known data is transferred from steps 𝑘 − 1 and 𝑘 to
step 𝑘 + 1.

(8.28) is stable only when 𝑟 ≤ 1. Detail analysis of the stability of (8.28) is not
discussed here.

Initial Conditions - Starting the Solution

The 3-level scheme poses some challenges when imposing the initial
conditions. If we imagine a row of false mesh points at time 𝑡 = −∆𝑡 = 𝑡 , then
the initial velocity condition (8.25) can be approximated using central
differences as:

100
 𝑢 −𝑢
= 𝑔(𝑥 )
2𝛥𝑡

𝑢 = 𝑢 − 2𝛥𝑡𝑔(𝑥 ) (8.29)
(8.28) when 𝑘 = 1

𝑢 =𝑟 𝑢 + 2(1 − 𝑟 )𝑢 + 𝑟 𝑢 −𝑢 (8.30)
Substituting (8.29) in (8.30) for 2nd (𝑡 ) time level
1
𝑢 = (𝑟 𝑢 + 2(1 − 𝑟 )𝑢 + 𝑟 𝑢 ) + 𝛥𝑡𝑔(𝑥 ) (8.31)
2
BC, IC, (8.28), and (8.31) are sufficient for knowing 𝑢 at all-time levels.

Example 𝟖. 𝟑 Solve 1D wave equation such that 𝑐 = 1 . Initial and

boundary conditions are given below:

BC: 𝑢(0, 𝑡) = 0 𝑢(𝐿, 𝑡) = 0

IC: 𝑢(𝑥, 0) = 𝑠𝑖𝑛(3𝜋𝑥 ⁄10)

𝜕𝑢
(𝑥, 0) = 0
𝜕𝑡
Consider 𝛥𝑡 = 0.05 ,𝛥𝑥 = 0.1. So that 𝑟 = 0.5 ≤ 1 for the stability of (8.28).

Scilab code 𝟖. 𝟑

//1D Wave Equation

clc
clear
//u(space,time)
//u(n,k)
//u(101,1201)
x(:,1)=[0:0.1:10]
u=zeros(101,1201)
dx=0.1//step size in x
dt=0.05//step size in t
c=1
r=0.5//Courant Number
function z=f(x)
z=sin(3*%pi*x/10)
endfunction
//using initial and boundary condition
u(1,:)=0

101
u(101,:)=0
u(:,1)=f((0:100)*dx)
for n=2:100
u(n,2)=0.5*(r^2*u(n+1,1)+2*(1-r^2)*u(n,1)+r^2*u(n-1,1))
end
for k=2:1200
for n=2:100
u(n,k+1)=r^2*u(n+1,k)+2*(1-r^2)*u(n,k)+r^2*u(n-1,k)-u(n,k-1)
end
end
realtimeinit(0.05)
figure(1)
realtime(0)
plot(x,u(:,1))
xlabel("x","Fontsize",4)
ylabel("u(x,t)","Fontsize",4)
title('1D Standing wave',"Fontsize",4)
h_compound = gce()
h_axes = gca()
h_axes.data_bounds = [0,-1;10,1]
for i=1:1200
realtime(i)
drawlater()
h_compound.[Link] = [x,u(:,i+1)]
drawnow()
end

102
Graphic window 𝟖. 𝟑

Keep Learning, Keep Growing

103
 About the Author

Rohan Verma
Academic Qualifications:

• [Link]. Physics, Indian Institute of Technology (IIT) Kanpur, India

• [Link]. Physics, University of Delhi, India

Connect with Me:

• ¯ LinkedIn: rohan-verma-0534a71ab

• Email: rohanv.i341@[Link]