M-Tech 1st Year

CACE Lab
COMPUTATIONAL LAB

ROLL NUMBER: 320206204002

  Newton Raphson Method:-
Definition:

In Numerical analysis, Newton’s method is known as the Newton-

Raphson method named after Isaac Newton and Jospeh Rapshon is a
method of finding successively better approximations to the roots. This
method is also called tangent method. This is a cord method in which we
approximate the curve near a root by straight line.

Derivation:

Let x0 be an initial approximation to the root of f(x)=0 then P(x o,f0) i.e.,
where f0 = f(x0) is a point on the curve, there we draw a tangent to the
curve at point P. We approximate the curve at the point M the
neighborhood of the root by the tangent to the curve at the point P.
Then point of intersection of the tangent with the x- axis is taken as the
next approximation to the root. The process is repeated until the required
accuracy is obtained. The equation of the tangent to the curve y=f(x) at
point P(xo,f0) is given by,

                                y - f(x0) = (x – x0) f1(x0)

Where f1(x0) is the slope of the tangent to the curve at P. Setting y= 0
solving for x, we get
                        X = x0 – f(x0)/ f1(x0) ,        f1(x0) ≠ 0

The next approximation to the root is given by,

                        X1 = x0 – f(x0)/ f1(x0) ,        f1(x0) ≠ 0

By repeating, we get the iteration method as,

 
Description:

Convergence of the Newton's method depends on the initial approximate

to the root. If the approximation is for away from the exact root, the
method diverges. Here we observe that the method may fail when f'(x) is
close to zero in the neighborhood of root, However if a root lies in small
interval (a,b) and if x0 € (a,b) then the method converges.

Question:

A gaseous reactant A, decomposes as follows A → 3R –r A = 0.6min-1

CA, find the conversion xA of A in a 50%A and 50% inerter feed given
mol
Ɣ0= 180 lit/min CA0 = 300 mill . lit to a lm3 mixed flow reactor.
£A A→3R
£A = 3*0.5+0.5 -1 /1 =1
For MFR
V X {A}
F { A {0 }} = τ /Ca 0 = −l { A }

C { A {0 }}( l − x )
C { A }= { A}

l+ £ { A } x { A } ¿
¿

C { A {0 }} x { A } ( 1+ x { A } ) Ɣ 1000 lit
τ=
0.6 C { A { 0}} ( 1− x
= Ɣ {0 }
= 180lit /min
{ A} )

x { A }+ x 2{ A }
= 5.56 (1-xA)
0.6

x 2+ 4.33 x { A } −3.33=0=f ( x )

Here assuming initial approximation x0 =0.5

Applying Newton's Raphson method

f ( x {0 } )
x1= x0 – f ' ( x {0 })
here f’(x) = 2xA + 4.33

−0.915
x1 = 0.5 - 6.33 = 0.67167

0.02947
x2 = 0.67167 – 6.35.673343 = 0.66647

−2.639 x 10− 6
x3 = 0.66647 - 5.66294
= 0.66647

∴ xA= 0.66647

MATLAB Code:

% Newton Raphson method

clear all
close all
clc
f=@(x) x^2+(4.33*x)-3.33
df=@(x) (2*x)+4.33
a=0.5, b=1;
x=a;
for i=1:1:100

x1=x-(f(x)/df(x));
x=x1;

end
sol=x;
fprintf('Approximate Root is %.15f',sol)
a=0.5, b=1;
x=a;
er(5)=0;
for i=1:1:5
x1=x-(f(x)/df(x));
x=x1;
er(i)=x1-sol;

end
plot(er)
xlabel('number of iterations')
ylabel('error')
title('erroe vs number of iterations')

Output:

a = 0.5
Approximate root is 0.6664704

 Method of false position:

Definition:
This method is also called linear interpolation method or Chord
method. Sometimes also called as Regula falsi method.

Derivation:

Let the root of the equation f(x)=0 lie in the interval (x k-1, xk) then P(xk-1,
fk-1), Q(xk , fk) are points on the curve f(x) = 0. Drawing the chord joining
the points P & Q taking the chord f(x)=ax+b, the next approximation is
b
given by x = − a

∵ Chord passes through P & Q

fk-1 = axk-1+b, fk = axk+b

Subtracting
fk - fk-1 = a(xk - xk-1)
a = [fk – fk-1]/[xk – xk-1]
b = fk - axk , the next approximation

b f {k } − a x { k }
xk+1 = − a = − a

x {k } − x { k− 1}
xk+1 = xk –( f {k } − f {k −1 } ¿ ) fk
¿
Description:

At the start of each iteration the required root lies in on interval, whose
length is decreasing. Hence the method converges. The disadvantage is
if the root lies initially in the interval (x 0,x1) the one of the end points is
fixed the end points is fixed for all iterations.

Question:
2 ( 1+0.2 ) s 2
R→ 0.2 –––– O ––––– 1+ 0.025
––––––– 2
s + 11s+ 1
––––––––––––> Y(s)

0.2
Find the roots of the characteristic equation,
 2 ( 1+0.2 s ) 2
x
y (s ) 1+ 0.02 s s2 +11 s+1
=
R( s) 0.2 x 2 ( 1+0.2 s ) x 2
1+
(1+ 0.02 s ) ( s2 +11s +1 )
Characteristic equation,

0.2 x 2 (1+0.2 s ) x 2
1+ ( 1+ 0.02 s ) ( s2 +11 s+1 )
=0

( 1+0.02 s ) ( s 2+11 s+1 ) + 0.2 x 2 ( 1+ 0.2 s ) x 2 = 0

0.02s3 + 1.22s2+12.62s+18=0

S3+16s2+631s+90=0

Let root lies between (-13,-14)

f(-13) = -1 f(-14)=468

1( − 1)
x1 = -13 - − 468+ 1 = 13.0021322
[-13.0021322, -14]
f(a) = -0.04467 f(b) = 468

−0.04467 − 0.9978678
x2 = -13.0021322 - − 468+0.0467 = -13.002227

[-13.002227, -14]
f(a) = -1.9996x10-3 f(b)= 468
−1.9996 x 10 −3(−0.997773)
x3 = -13.002227 - 468+1.9996 x 10 − 3
= -13.00223

∴ x = -13.0022

MATLAB Code:
% False position method
clc;
f=@(x) x^3+(61*(x^2))+(631*x)+90;
a=-13; b=-14;
for i=1:10
x0=a;x1=b;
fprintf('\n hence root lies between (%.4f,%.04f)',a,b);
x2(i)=x0-((x1-x0)/(f(x1)-f(x0)))*f(x0);
if f(x2(i))>0
b=x2(i);
else a=x2(i);
end
fprintf('\n therefore, x2=%.4f\n here, f(x2)=%.4f',x2(i),f(x2(i)))
p=x2(i);
end
for i=1:10
error(i)=p-x2(i);
end

Output:
here the root lies between (-13,-14)
therefore X2 = -13.00

 R-k 4th order Method :

Definition :

Runge kutta formula’s are among the oldest and best understood
schemes in numerical analysis and provide a proper way to solve the
initial value problem of a system of ODE’s.

Derivation :

If is used to solve ordinary differential equation of the form,

dy
=f ( x , y )
dx
From the derivation of R-k 2nd order method taking with respect to 4
terms, we get,

For 4th order method,

k1 = f(x,y)

ℎ 1
k2 = f(xi + 2 , yi+ 2 k1h)

ℎ 1
k3 = f(xi + 2 , yi+ 2 k2h)

k4 = f(xi + h , y+ k3h)
 1
yi+1 = yi + 6 [k1 + 2k2 + 2k3 +k4] x h

Description :

Only first order differential equations can be solved by using R-K 4 th

order method.

Question:

An ageous feed containing A (1 mol/lit) enter 2-lit PFR and reacts

away(2A->R -ra=0.05 CA2 mol/lit-s) .Find the outlet concentration
of A for a feed rate of 0.5 lit/min.
Sol: CA0 =1 mol/lit V=2lit -ra=0.05 CA2 k1=0.05 F0= 0.5 lit/min
Constant volume system
dv
dt
=F0 − F F0=F=0.5 lit/min
dV C A
F 0 C A 0 − F C A −k 1 V C A =
dt
F0 CA 0 F CA dCA
− − k CA= =f (t , x)
V V dt
0.5 ∗1 0.5 C A dCA
− − 0.05C A =
2 2 dt
dCA
0.25-0.3C A − dt =f (t ,C A )
Here, k 1=ℎ ∗ f (t i , C Ai )
ℎ k1
k 2=ℎ ∗ f (t i+ ,C Ai + )
2 2
ℎ k2
k 3=ℎ ∗ f (t i+ , C Ai + )
2 2
k 4=ℎ∗ f ( t i+ ℎ ,C Ai +k 3 )
1
C Ai+1=C Ai+ (k 1 +2 k 2 +2 k 3 +k 4 )
6
1st iteration taking h=2,
f ( 0,1 ) k 1 =−0.1
f ( 1,0.95 ) k 2 =−0.07
f ( 1 , 0.965 ) k 3=− 0.079
f ( 2,0.921 ) k 4=− 0.0526
CA1=0.9249
2nd iteration
k 1=−0.05494 ∗ f ( 2,0.9249 )
k 2=−0.038458 ∗ f (3 , 0.89743)
k 3=−0.0434 ∗ f ( 3,0.905671 )
k 4= -0.02889 * f(4, 0.8814974)
CA2=0.8836

MATLAB Code:
x0=input('enter the values of x0');
y0=input('enter the values of y0');
h=input('enter the values of h');
xf=input('enter the values of xf');
x=x0:h:xf;
y1=y0;
rkf=@(x,y) (0.25-(0.3*y));
for i=1:length(x)-1
k1=h*rkf(x(i),y(i));
k2=h*rkf((x(i)+(h/2)),(y(i)+(k1/2)));
k3=h*rkf((x(i)+(h/2)),(y(i)+(k2/2)));
k4=h*rkf((x(i)+h),(y(i)+k3));
y(i+1)=y(i)+((k1+(2*k2)+(2*k3)+k4)/6);
x(i+1)=x(i)+h;
end
disp(x);
disp(y);
disp('rk 4th mehod');
disp('iteration x y');
for i=1:length(x)
fprintf('%d\t %f\t %f\n',i,x(i),y(i));
end
Output:

enter the values of x0 0

enter the values of y0 1
enter the values of h 2
enter the values of xf 1
i2
CA 0.88365

 Gauss Elimination Method:

Definition:
This method is basically used to solve simultaneous linear
equations, done by an elimination process.

Derivation:

This method is based on the idea of reducing the given system of

equations Ax=b to an Upper triangular system of equations Ux=z using
elementary row operations.
Considering a system a
a 11 x 1 +a12 x 2 +a 13 x 3=b1
a 21 x1 + a22 x 2+ a23 x 3=b2
a 31 x1 +a 32 x 2+ a33 x 3=b 3
a11 a 12 a13 x1 b1

[ a21 a 22 a23
a31 a 32 a33 ][ ] [ ]
x2
x3
= b2
b3
By row operations.
a11 a12 a 13 x 1 b'

[ ][ ] [ ]
0 a22 a 23 x 2 = b' '
0 a32 a 33 x 3 b '' '
then finally
 a11 a12 a 13 x 1 b'

[ 0 ][ ] [ ]
0 a22 a 23 x 2 = b' '
0 a 33 x 3 b '' '

Description:
Sometimes it fails when anyone of the pivot is zero or a very small
number as the elimination progresses. The element a 332≠0 is called the
third pivot. This system uses upper triangular form.

Question:
Consider a liquid entering into a pipe of length L and diameter D at
a fixed pressure Pi , as shown in the fig. below and the flow distributes
itself into 2 pipes each of lengths L1, L2 and diameters D1, D2 with fixed
exit pressure Pa then determine the velocities of the flowing each pipe.

32 µVL
From laminar flow ∆p= D2
Q=Q1+Q2
VA¿ V 1 A 1+V 2 A2
2 2
Vπ D 2 V 1 π D 1 V 2 π D 2
= +
4 4 4
At function
 32 µVL
Pi-P¿ D2
32 µV 1 L1
P-Pa¿ D2 pipe-2
1
32 µV 2 L2
P-Pa¿ D2 pipe-3
2
32 µVL 32 µ V 1 L1
Pi-P+P-Pa¿ D2 + D2
1
32 µVL 32 µ V 1 L1
Pi-Pa¿ D2 + D2
1
2 2 2
D V − D 1 V 1 − D 2 V 2=0
D2 − D21 − D22

[ 32 µL
D2
32 µL
D2
L=5m, L1=4m,
32 µ L1
D

0
2
1

L2=3m
0

32 µ L2
D22
] V

[][ ]
V1
V2
=
0
Pi − P a
Pi − P a

D=2m D1=D2=1m µ=0.1P

2
Pi=500 N/m Pa=300 N/m2
4 −1 − 1 V 0

[ 0.4 1.28
0.4 0
0 V
][ ] [ ]
1
0.96 V 2
= 200
200
1 −0.25 −0.25 0

[ 0.4 1.28
0.4 0
0 200
0.96 200
1 −0.25 − 0.25 0
]
R2->R2-0.4R1 R3->R3-0.4R1

[ 0 1.18 0.1 200

0 0.1 1.06 200 ]
1 − 0.25 − 0.25 V 0

[ 0 1.18
0 0
0.1
][ ] [ ]
V 1 = 200
1.0515 V 2 183.05
V2=176.2 m/s V1=132.15 m/s V=77.09 m/s

MATLAB Code:

% solve Ax=b using gauss elimination

A = [1 -1/4 -1/4; 0.4 1.28 0; 0.4 0 0.96];
b = [0;200;200];
%% Gauss elimination
% Get augumented matrix
Ab = [A, b];
n=3;
% With A(1,1) as pivot element
alpha = Ab(2,1)/Ab(1,1);
Ab(2,:) = Ab(2,:)- alpha*Ab(1,:);
alpha = A(3,1)/A(1,1);
Ab(3,:) = Ab(3,:)-alpha*Ab(1,:);
% With A(2,2) as pivot element
alpha = Ab(3,2)/Ab(2,2);
Ab(3,:) = Ab(3,:)-alpha*Ab(2,:);
%% Back-substitution
x = zeros(3,1);
for i=3:-1:1
x(i) = ( Ab(i,end)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);
end

Output:

A= [1 -1/4 -1/4; 0.4 1.28 0; 0.4 0 0.96]

i=1
n=3
X=[77.0925,132.158,176.217]

 Trapezoidal method:
Definition :
The Trapezoidal rule is a numerical method that approximates the
b

value of a definite integral as ∫ f ( x )dx

a

Derivation:
The curve y=f ( x ) ; a ≤ x ≤ b be approximated by the line joining the point P,
Q

Using Newton’s forward difference formula the linear polynomial

approximation to ( x) , interpolating at the points P, Q is given by
1
f ( x )=f ( x 0 ) + ( x − x 0 ) ∆ f (x 0 )
ℎ
x 0=a , x 1=b , ℎ=b − a
b x1

I =∫ f ( x ) dx=∫ f ( x) dx
a x0
x1 x1

¿ f ( x 0 ) ∫ dx +
0x
1
ℎ [∫
x0
]
( x − x 0 ) dx ∆ f 0

1 1 x
ℎ 2 x0 |{
¿ ( x 1 − x0 ) f ( x0 ) + ∗ ( x − x 0 )2 1 ∆ f 0

ℎ
¿ ℎ ∗ f ( x0 )+ [ f ( x1 ) − f x0¿ ]
2
ℎ
¿ [ f (x 1)+ f ( x 0) ]
2
b

∫ f ( x ) dx= b −2 a [ f ( b ) +f (a)]
a

Description :
Basically, the right hand side of the trapezium rule is the area of
the trapezoid with width b-a, and ordinates f ( a ) , f ( b ) which is an
approximation to the area under the curve y=f ( x )
Question:
A gaseous reactor A, decomposes as A->3R and -r a=0.6CA in
batch reactor with conversion XA=0.925 feed as 50% feed and 50% inert
feed and CA0¿ 300 m. mol /lit to a lm3 reactor . Calculate the volumetric flow
rate of the initial feed .
XA
V T d XA
= =∫ - l A =0.6 C A
F A 0 C A 0 0 −l A
0.925
1000lit ( 1+ X A )
=∫ d XA
γ0 0 ( 0.6 ( 1− X A ) )
0.925
0.6 1+ X A
1000 ∗ =∫ d XA
γ0 0 1− X A

X 0 0.1825 0.37 0.555 0.75 0.925

F(X) 1 1.45396 2.175 3.4944 7 25.6667
Here h=0.1825
0.1825
¿ [ 1+ 1.45396+2.175+3.4944+7 +25.6667 ]
2
¿ 3.773
γ 0.6
0=¿1000 ∗ =159 lit / min¿
3.7773

MATLAB Code:

clc
clear all
close all
f=@(x) ((1+x)/(1-x));
x0=input('enter x0');
xn=input('enter final value xn');
N=input('enter total num of steps');
h=((xn-x0)/N);
area=0;
while (x0<xn)
area=area+(h/2)*(f(x0)+f(x0+h));
x0=x0+h;
end
fprintf('area is %f',area);

Output:

Enter x0 0
Enter final value xn 0.925
Enter no of steps 5
Area is 5.02249

 Gauss Seidel Iteration Method:

Definition:

Sometimes, this method is a base modification of Jacobi

Methods ,We assume that the pivots a ii ≠0 , for all i .

Derivation:
Taking the equations as,
a 11 x 1=b 1 − ( a12 x 2 +a13 x3 )
a 22 x2 =b2 − ( a 21 x 1+ a23 x 3 )
a 33 x 3=b3 − ( a31 x 1+ a32 x 2 )
From this we take the Jacobi Iteration as
1
x k+1
1 = [ b1 − ( a12 x2k + a13 x k3 ) ]
a 11
 1
x k+1
2 = [b2 − ( a21 x k1 +a23 x3k )
]
a 22
1
x k+1
3 = [b − ( a31 x k1 +a 32 x2k ) where
] k=0,1,2,…..
a 33 3

Description:
Here, we replace the complete vector xk, in the right hand side of
the above equation at the end of each iteration. This method is also
called the method of simultaneous displacement.

Question:
An ageous solution of pyridine containing 27% by wt of pyridine and
75% water is to be extracted with chlorobenzene. The feed and solvent
are mixed well in a batch extractor then allowed for phase separation.
The extract phase contains 11% pyridine 88% chlorobenzene 1% water
by weight raffinate phase contains 6% pyridine, 94% water by weight
.Calculate the quantities of 2 phases.
Feed +solvent= raffinate phase +extract phase
100+z=x +y (1)
Material balance for pyridine
100*0.27+0=x*0.11+y*0.06
0.11x+0.06y=27 (2)
Material balance for water
100*0.73+0=x*0.01+y*0.94
0.01x+0.94y=73 (3)
The simultaneous linear equations are
x +y-z=100
0.11x+0.06y=27
0.01x+0.94y=73
From this z=x+y-100
0.06 x
x=27 −
0.11
0.01 x
y=73 −
0.94
1st Iteration:
Initially x0=0(say)
X1=27/0.11 =245.45
Y1= (73-0.01(245.45))/0.94=75.048
Z1=245.45+75-100=220.49
2nd Iteration
X2 = (27-0.06(75.048))/0.11=204.519
Y2= (73-0.01(201.519))/0.94=75.48
Z2=204.519+75.48-100=179.99

MATLAB Code:

% Solving Ax = b using Gauss Siedel method

A = [1 1 -1; 0.11 0.06 0; 0.01 0.94 0];
b = [100;27;73];
Ab = [A b];
% Rearrange to get diagonally dominant matrix
As = [Ab(2,:); Ab(3,:); Ab(1,:)];
% Initializing
n = 3;
x = zeros(n,1);
err = zeros(n,1);
%% G-S Ierations
for iter = 1:3
for k = 1:n
xold = x(k);
num = As(k,end) - As(k,1:k-1)*x(1:k-1) - As(k,k+1:n)*x(k+1:n);
x(k) = num/As(k,k);
err(k) = abs(x(k)-xold);
end
disp(['Iter',num2str(iter),'; Error = ', num2str(max(err))]);
end
Output:

iter 3
k3
n3
X [204.2815;75.486;179]

 Simpson’s 1/3rd rule:

Definition:
Simpson’s 1/3rd rule is an extension of trapezoidal rule when the
integrand is approximated by a 2nd order polynomial

Derivation:
Approximating f (x) by a 2nd order polynomial as
( a+b )
(
f 2 ( x)=b0 +b1 ( x − a ) +b2 ( x − a ) x −
2 )
Where b0= f ( a )
( a+b )
b=1
( f
2 )
− f (a )

a+b
−a
2
b 2=¿)/(b-a)
b b

∫ f ( x ) dx =∫ f 2 ( x ) dx
a a
b
a+ b
¿ ∫ b0 + b1 ( x −a ) + b2 ( x −a ) x −
a
( 2
dx )
2
¿ b0 x +b1 ([ x2
2
−ax ) +b3
3] [(
x 3 ( 3 a+ b ) x ( a ( a+b ) x )
−
4
+
2 )] b
a
b2 −a2
¿ b0 ( b − a ) +b1
b
[2
− a( b− a) +b 2 ¿ ]
∫ f ( x ) dx= b −6 a f ( a ) + 4 f a+2 b + f (b)
[ ( ) ]
a

ℎ a+ b
[3
¿ f ( a) + 4 f (
2 )
+ f (b)
]
Description:
It is a process used for integration. It has a mole accuracy
than trapezoidal rule, Since trapezoidal rule was based on approximating
the integrand by a first order polynomial and then integrating the
polynomial over interval of integration.

Question:
A gaseous reactant A decomposes as A3R, -rA=0.6 CA/min
with CA0=300m.mol/lit to a lm3 PFR. where conversion is XA=0.88 find
the volumetric flow rate.
For PFR
XA
V T
= ( d X A )/− r A ¿
F A0 C A0 ∫
=
0 ¿

A3R 50% feed 50% inerts

EA=3*0.5+0.5-1
0.88
1000 lit d X A ( 1+ x )
10 −3 mol
= ∫ 0.6∗ C A 0 ∗(1 − x )
300∗ ∗ γ0 0
lit
0.88
1000 lit 1+ x
0.6 ∗ =∫ dx
γ0 0 1−x
 1+ x
Here f ( x )= 1 − x
Taking step size=h=0.2
X 0 0.22 0.44 0.66 0.88
F(X) 1 1.5641 2.57143 4.8823 15.667

0.22
¿ [ 1+ 15.667+4 (1.5641+2.51743+4.8823) ]
3
=3.86748
1000
γ 0=0.6∗
3.86748
γ 0=156 lit /min

MATLAB Code:
clc
clear all;
close all;
f=@(x) 1+x/1-x;
a=0; b=0.88;
n=b-a;
h=(b-a)/n;
p=0;
for i=a:b
p=p+1;
x(p)=i;
y(p)=i+1/i-1;
end
l=length(x);
x
y
answer=(h/3)*((y(1)+y(l))+2*(y(3)+y(5))+4*(y(2)+y(4)+y(6)))

Output:

h 0.22
x 3.47
  Secant Method:

Definition:
One of the drawback of newton’s method is it is necessary to
evaluate f’(x) at various points, which may not be practical for some
choices of f. The secant method avoids this issue by using a finite
difference to approximate the derivative.

Derivation:
It can be noted that f(x) is approximated by a secant line through 2
points on the graph of f, rather than a tangent line through one point on
the graph.
Since a secant line is defined using 2 points on the graph of f(x), it
is necessary to choose 2 initial iterates.
From the newton Raphson method of solving a nonlinear equation
f(x)=0, we have
Xi+1 = Xi – (f(Xi)/f’(Xi))
To overcome the drawbacks the derivative of the function f(X) is
taken as
f’(Xi) = (f(Xi) – f(Xi-1))/(Xi – Xi-1)
Xi+1 = Xi – ( (f(Xi)* (Xi – Xi-1))/ f(Xi) – f(Xi-1))
Which is secant method.

Description:
This method uses 2 initial guesses but unlike the bisection method,
the 2 initial guesses do not need to bracketed the root of the equation.
This is an open method which may or may not converge, but if secant
method converges, it will typically converge faster than bisection
method, but it typically converges slower than Newton Raphson method.

Question:
 The floating ball has a specific gravity of 0.6 and radius 5.5cm, depth X
to which the ball is submerged under water is given by
X3 - 0.165*X2 + 3.993*10-4 = 0
Find the depth X to which the ball is submerged under water
Given
f(X) = X3 - 0.165*X2 + 3.993*10-4 = 0
let us assume initial guesses of the root as
X-1 = 0.02, X0 = 0.05
Iteration – 1
X1 = X0 – ( (f(X0)* (X0 – X-1))/ f(X0) – f(X-1))
Substituting we get
X1 = 0.06461
|Ea| = |(X1 – X0)/X1|*100%
=22.62%
Iteration – 2
X2 = X1 – ( (f(X1)* (X1 – X0))/ f(X1) – f(X0))
X2 = 0.06241
|Ea| = |(X2 – X1)/X2|*100%
= 3.525%
Iteration – 3
X3 = X2 – ( (f(X2)* (X2 – X1))/ f(X2) – f(X1))
X3 = 0.06238
|Ea| = |(X3 – X2)/X3|*100%
= 0.0595%
Here the root converges at X = 0.06238

MATLAB Code:

f=@(x) ((x^3)-(0.165*(x^2))+(3.993*0.0001));
x(0)=0.02;
x(1)=0.05;
n=0.0005
iteration=0;
for i=2:10
x(i)=x(i-1)-(f(x(i-1)))*((x(i-1)-x(i-2))/(f(x(i-1))-f(x(i-2))));
iteration=iteration+1;
if abs((x(i)-x(i-1))/x(i))*100<n
root=x(i);
iteration=iteration;
break
end
end

Output:
i3
x 0.06238

 Bisection Method:

Definition:

One of the first numerical methods developed to find the root of a

nonlinear equation f (x )=0 was the bisection method (also called
binary-search method). The method is based on the following theorem.
An equation f (x )=0 , where f (x ) is a real continuous function, has at
least one root between x ℓ and x u if f (x ℓ )f ( x u )< 0
Note that if f (x ℓ )f ( x u )> 0 , there may or may not be any root between
x ℓ and x u . If f (x ℓ )f ( x u )< 0 , then there may be more than one root

between x ℓ and x u . So the theorem only guarantees one root between

x ℓ and x u .

Derivation:
The steps to apply the bisection method to find the root of the equation
f (x )=0 are
1. Choose x ℓ and x u as two guesses for the root such that
f (x ℓ )f ( x u )<0 , or in other words, f (x ) changes sign between x ℓ

and x u .

2. Estimate the root, x m , of the equation f (x )=0 as the mid-point

between x ℓ and x u as
x ℓ +x u
xm =
2
3. Now check the following
a) If f (x ℓ )f ( x m )<0 , then the root lies between x ℓ and x m ; then
x ℓ=x ℓ and x u=x m .

b) If f (x ℓ )f ( x m )>0 , then the root lies between x m and x u ; then

x ℓ =x m and x u=x u .

c) If f (x ℓ )f ( x m )=0 ; then the root is x m . Stop the algorithm if this

is true.

4. Find the new estimate of the root

xℓ + xu
xm =
2
Find the absolute relative approximate error as
x new
m - x old
m
|∈a| = | |× 100
x mnew
where
new
xm = estimated root from present iteration
old
xm = estimated root from previous iteration
5. Compare the absolute relative approximate error |∈a| with the pre-
specified relative error tolerance ∈s . If |∈a|>∈s , then go to Step
3, else stop . Note one should also check whether the number of
iterations is more than the maximum number of iterations allowed.
 If so, one needs to terminate the algorithm and notify the user
about it.

Description:

 At least one root exists between the two points if the function is
real, continuous, and changes sign

f (x)

xℓ
x
xu

 If the function f (x ) does not change sign between the two points,
roots of the equation f (x )=0 may still exist between the two
points
f (x)

x
xℓ xu
  If the function f (x ) does not change sign between two points,
there may not be any roots for the equation f (x )=0 between the
two points

f (x) f (x)

xℓ xu
x x
xℓ xu

 If the function f (x ) changes sign between the two points, more

than one root for the equation f (x )=0 may exist between the two
points f (x)

xu
x
xℓ

Question:

You have a spherical storage tank containing oil. The tank has a
diameter of 6 ft. You are asked to calculate the height h to which a
dipstick 8 ft long would be wet with oil when immersed in the tank
3
when it contains 4 ft of oil. The equation that gives the height, h , of
the liquid in the spherical tank for the given volume and radius is given
by :
3 2
f ( h )=h −9 h +3 .8197=0

Solution
From the physics of the problem, the dipstick would be wet between
h=0 and h=2r , where
r= radius of the tank,
that is
0≤h≤2 r
0≤h≤2(3)
0≤h≤6
Let us assume
hℓ =0 , hu =6
Check if the function changes sign between hℓ and hu .
3 2
f (hℓ )=f ( 0 )=( 0 ) −9 ( 0 ) +3. 8197=3 . 8197
3 2
f (hu )=f (6)=(6 ) −9(6 ) +3 . 8197=−104 .18
Hence
f ( hℓ ) f ( hu ) =f ( 0 ) f ( 6 )=( 3 . 8197 ) (−104 .18 )< 0

So there is at least one root between hℓ and hu that is between 0 and

6.

Iteration 1
The estimate of the root is
h ℓ +hu
hm =
2
0+6
=
2
=3
f ( hm )=f ( 3 )=( 3 )3 −9 ( 3 )2 +3 .1897=−50 .180
f ( hℓ ) f ( hm )=f ( 0 ) f ( 3 )= ( 3. 1897 ) (−50 .180 ) < 0

Hence the root is bracketed between hℓ and hm , that is, between 0 and
3. So, the lower and upper limits of the new bracket are
hℓ =0 , hu =3
At this point, the absolute relative approximate error |∈a| cannot be
calculated, as we do not have a previous approximation.

Iteration 2
The estimate of the root is
h ℓ +hu
hm =
2
0+3
=
2
=1.5
f ( hm )=f ( 1. 5 )=( 1 .5 )3−9 (1 . 5 )2 +3 . 8197=−13. 055
f ( hℓ ) f ( hm )=f ( 0 ) f ( 1 .5 ) =( 3. 8197 )(−13 . 055 ) <0

Hence, the root is bracketed between hℓ and hm , that is, between 0

and 1.5. So the lower and upper limits of the new bracket are
hℓ =0 , hu =1 .5
The absolute relative approximate error |∈a| at the end of Iteration 2 is
new old
hm −h m
|∈a|=| new
|×100
hm
1.5−3
=| |×100
1.5
=100 %
None of the significant digits are at least correct in the estimated root
hm=1 . 5
as the absolute relative approximate error is greater that 5% .

Iteration 3
The estimate of the root is
h ℓ +hu
hm =
2
0+1 .5
=
2
=0. 75
f ( hm )=f ( 0 .75 ) =( 0 .75 )3 −9 ( 0 . 75 )2 +3 . 8197=−0 .82093
f ( hℓ ) f ( hm )=f ( 0 ) f ( 0 . 75 )=( 3 . 8197 ) (−0 .82093 )< 0
Hence, the root is bracketed between hℓ and hm , that is, between 0 and
0.75. So the lower and upper limits of the new bracket are
hℓ =0 , hu =0 . 75
The absolute relative approximate error |∈a| at the end of Iteration 3 is
new old
hm −h m
|∈a|=| |×100
hnew
m
0.75−1.5
=| |×100
0.75
=100 %
Still none of the significant digits are at least correct in the estimated
root of the equation as the absolute relative approximate error is greater
than 5% .
The height of the liquid is estimated as 0.75 ft at the end of the third
iteration.
Seven more iterations were conducted and these iterations are shown in
the table:
Iteratio hℓ hu hm |∈a|% f ( hm )
n
3
0.00
1.5
0.00 ---------
0.75
1 0.00 6 - −50.180
0.375
2 0.00 3 100 −13.055
0.5625
3 0.375 1.5 100 −0.82093
0.6562
4 0.5625 0.75 100 2.6068
5
5 0.6562 0.75 33.333 1.1500
0.7031
6 5 0.75 14.286 0.22635
3
7 0.6562 0.75 6.6667 −0.28215
0.6796
8 5 0.70313 3.4483 −0.024077
9
9 0.6562 0.67969 1.7544 0.10210
0.6679
10 5 0.67969 0.8695 0.039249
7
0.6679 7
0.6738
7
3
 th
At the end of the 10 iteration,
|∈a|=0.86957%
The estimated root is 0.67383.

MATLAB Code:

clc;
clear all;
close all;
a=0;
b=6;
bif=@(x) (x^3-(9*(x^2))+3.8197);
for i=1:100
c=(a+b)/2;
if bif(a)*bif(c)<0
b=c;
else
a=c;
end
end

Output:

i 100
c 0.6775

 Gauss Jordan Method:

Definition:

The Gauss-Jordan method is also known as Gauss-Jordan

elimination method. It was introduced by the mathematicians Carl
Friedrich Gauss and Wilhelm Jordan, after their name it is called so.
This method is very useful in solving a linear system of equations.
It is a technique in which a system of linear equations is resolved by the
means of matrices. This method allows the isolation of the coefficients
of a system of linear equations. In Gauss-Jordan method, given matrix
can be fetched to row echelon form and made simpler.

Derivation:
There are 3 kinds of elementary row operations which are discussed
below:

i) Switching of Rows:
A row can be interchanged with another row.
Ri↔Rj

ii) Row multiplication:

A non-zero number can be multiplied to every element in a row.
aRi→Ri , where a ≠ 0

iii) Addition of Rows:

We may replace a row by the sum of elements of a row and a multiple of
corresponding elements of another row.
Ri+aRj→Ri , where i ≠ j

Row Reduced Echelon Form

A matrix can be called to be in reduced row echelon form (RREF) if the
following conditions are satisfied:
1) The first nonzero element of a row must be 1 which is known as
leading 1.
2) The leading 1 of a row must be positioned at the right side of the
leading 1 of its previous row.
3) In case of a column containing a leading 1, all the remaining elements
in that column must be 0.
4) At the bottom of the matrix, there is a rows having only zeros as the
elements
Description:
Step 1: Write the augmented matrix for the given system of linear
equations.
Step 2: Use a set of elementary row operations in order to perform row
reduction on this matrix until we obtain a unique reduced row echelon
form.
Step 3: Once RREF is obtained, write down the system of linear
equations from this form.
Step 4: In this way, we get a value corresponding to each unknown
variable which would be the solution. If simple equations are obtained,
then they should be solved by the methods of solving equation such as
substitution method or elimination method.

Question:

The waste acid from a nitrating process contains 21% HNO 3, 55% H2SO4
and 24% H2O by weight. This acid is concentrated to contain 28%
HNO3, 62% H2SO4 by the addition of concentrated sulphuric acid
containing 93% H2SO4 and concentrated nitric acid containing 90%
HNO3. Calculate the weights of waste and concentrated acids that must
be combined to obtain 1000 lb of the desired mixture.

Solution:
93% H2SO4 , 90% HNO3

21% HNO3, 28% HNO3

55% H2SO4
Mixing 62% H2SO4
24% H2O
10% H2O
Let,
X be the spent acid
Y be the conc. H2SO4
Z be the conc. HNO3
X+Y+Z = 1000
X(0.21) + Z(0.9) = 1000 (0.28) HNO3
X(0.55) +Y(0.93) = 1000 (0.62) H2SO4
X(0.24) + Y(0.07) + Z(0.1) = 1000 (0.10) H2O

These above equations are represented in the form of a matrix as below:

0.24 0.07 0.1 X 100

0.55 0.93 0 Y = 620
0.21 0 0.9 Z 280

0.24 0.07 0.1 100

0.55 0.93 0 620 R1 / 0.24
0.21 0 0.9 280

1 0.291667 0.4167 416.67

0.55 0.93 0 620 R1 x 0.55
0.21 0 0.9 280

Solving similarly to reduce into Echelon form

1 0 0 126.789
0 1 0 591.68
0 0 1 291.52

Therefore,
X = 126.789;
Y = 591.68;
Z= 291.52

MATLAB Code:

a=[ 0.24 0.07 0.1 100

0.55 0.93 0 620
0.21 0 0.9 280];
%Gauss jordhan method
[m,n]=size(a);
for j=1:m-1
for z=2:m
if a(j,j)==0
t=a(1,:);a(1,:)=a(z,:);
a(z,:)=t;
end
end
for i=j+1:m
a(i,:)=a(i,:)-a(j,:)*(a(i,j)/a(j,j));
end
end
for j=m:-1:2
for i=j-1:-1:1
a(i,:)=a(i,:)-a(j,:)*(a(i,j)/a(j,j));
end
end
for s=1:m
a(s,:)=a(s,:)/a(s,s);
x(s)=a(s,n);
end
disp('gauss jordhan method');
a
x

Output:
X= [126.789 591.68 281.527]

 Implicit and Explicit methods:

Implicit Method:

Definition:

In numerical analysis and scientific computing, the backward

Euler method (or implicit Euler method) is one of the most
basic numerical methods for the solution of ordinary differential
equations.

Derivation:

Consider the ordinary differential equation

dy/dx = f(t,y)

integrating the differential equation from tn to tn+1 yields a differential

form and approximating the integral on right by rectangle method finally
we get
y(tn+1)-y(tn) = h*f(tn+1,y(tn+1))
here, yn is supposed to approximate y(tn) and is as given above.

Description:

For the given differential equation, a numerical method produces a

sequence y0,y1,y2,……… such that yk approximates y(t+kh), here h is
the step size. Generally the backward euler method is an implicit
method, where the new approximation yk+1 appears on both sides of the
equation and thus the method needs to solve an algebraic equation for
the unknown yk+1.

Explicit Method:

Definition:
This explicit method calculate the state of a system at a later time from
the state of the system at the current time.

Derivation:

Consider an ordinary differential equation as

dy/dt = -y2
with initial condition y(0)= 1
consider a grid tk=a*(k/n)
for 0<k<n;
where time step delta t = a/n
yk=y(tk)
for explicit method we have

yk+1 = yk – deltat*yk^2

Description:

Generally, when a direct computation of the dependent variables can be

made interms of known quantites, the computation is said to be explicit.

Question:

Astell ball of 50mm diameter is cooled by exposing it to an air stream at

320k,under these conditions the convective heat transfer coefficient is
100w/m^2-k. Initially the ball is at a temperature of 1120k what would
be the temperature of the ball after 415.95sec

Solving with respect to heat balance we get

dT/dt = -3.33*10^-3*T + 0.106668

for implicit method

x0=0,y0=1120, h=207.975
solving we get
y1= 792.46
y2=599

for Explicit method

Solving as mentioned above we get

y1=366
y2=135

MATLAB Code:
clc
clear all
close all
% y' =-2x-y
% initial condition
x(1) = 0;
y(1) = -2;
h = 0.1;% step size

%Euler Forward
%%%%%%%%%%%%%%
for i = 1:10
y(i+1) = y(i)+h*(-2*x(i)-y(i));
x(i+1) = x(i)+h;
end

xlabel('x')
ylabel('y')

%Euler Backward
%%%%%%%%%%%%%%
% yn+1 = (yn -2hxn+1)/(1+h)
for i = 1:10
x(i+1) = x(i)+h;
y(i+1) = (y(i)-2*h*x(i+1))/(1+h);
end
xlabel('x')
ylabel('y')

Output:

%Euler forward
X 415.95
Y 599

%Euler backward
X 415.95
Y 135.23