EXPERIMENT NO.

NAME OF THE EXPERIMENT

Non-linear Algebraic Equations (Bisection Method)

AIM:

Write a program to find out the root of equation x3 1.8x2 10 x 17 0 and roots lie
between [1,2] by using Non-Linear Algebraic Method (Bisection Method) using Scilab.

PRE-REQUISITE:

1. Knowledge of SCILAB Programming

2. Knowledge of Numerical Methods

THEORY:

The bisection method in mathematics is a root finding method which repeatedly

bisects an interval and then selects a subinterval in which a root must lie for further
processing. It is a very simple and robust method, but it is also relatively slow. Because of
this it is often used to obtain a rough approximation to a solution which is then used as a
starting point for more rapidly converging methods. The method is also called Binary Search
method and Internal Halving method, and is similar to the concept of binary search, where
the range of possible solutions is halved for each iteration.

The method is applicable when we wish to solve the equation f ( x) 0 for the real

variable x, where f is a continuous function defined on an interval a, b while f a and

f b have opposite signs. In this case a and b are said to bracket a root since, by the

intermediate value theorem, the f must have at least one root in the interval a, b .

At each step the method divides the interval in two by computing the midpoint
c (a b) / 2 of the interval and the value of the function f c at that point. Unless c is itself

a root (which is very unlikely, but possible) there are now two possibilities: either f a and

f c have opposite signs and bracket a root or f c and f b have opposite signs and

bracket a root. The method selects a subinterval that is bracketed as new interval to be used in
the next step. In this way the interval that contains a zero of f is reduced in width by 50% at
each step. The process is continued until the interval is sufficiently small.

1
 Explicitly, if f a and f c are opposite signs, then method sets c as the new value

for b, and if f b and f c are opposite signs then the method sets c as new a. In both cases,

the new f a and f b have opposite signs, so the method is applicable to this smaller

interval.

ALGORITHM FOR BISECTION METHOD:

1. Given a continuous function f x

2. Find points a and b such that a b and f a f b 0 .

3. Take the interval a, b and find its midpoint x1 .

4. If f x1 0 then x1 is an exact root, else if f x1 f b 0 then let a x1 , else if

f a f x1 0 then let b x1 .

5. Repeat steps 2 and 3 until f xi 0 or f xi DOA ( DOA degree of accuracy).

Figure 1: Bisection Method

For the i th iteration, where, i 1, 2,..., n , the interval width is

xi 0.5

xi 1 0.5 b a
i

And the new midpoint is

xi ai 1 xi

2
PROGRAM IN SCILAB
function [ ] = bisectionMethod( a,b )
a=3;
b=5;
MaxIt=1000;
epsilon=10^-6;
c=(b+a)/2;
NumIt=0;
while NumIt<MaxIt && abs (f(c))>epsilon
if f(a)*f(c)<0
b=c;
else
a=c;
end
NumIt=NumIt+1;
c=(b+a)/2;
end
c
end
function y=f(x)
y=x^3-1.8*x^2-10*x+17;
end

OUTPUT
C = 1.6618

RESULTS:
The root of the given Non-linear Algebraic Equation by Bisection Method is: 1.6618

CONCLUSION:
Thus we successfully studied the root finding program of Non-linear Algebraic Equation by
Bisection Method.

3
INDUSTRIAL APPLICATIONS:
1. It is used in computer science research to analyze safeguard zero finding methods
2. It is simplest of other all methods
3. We can safeguard bisection to detect cases where we dont have any roots

VIVA QUESTIONS:
1. What is the difference between linear & nonlinear methods?
2. How does Bisection method differ from other numerical methods used for solving
linear algebraic equations?
3. What do you mean by Iterative process?
4. Bisection method uses Continuous functions or Non continuous functions?
5. Find root of following polynomial equation using Bisection method: f(x)=x3-x-2

4
 EXPERIMENT NO. 2

NAME OF THE EXPERIMENT

Non-linear Algebraic Equations (Secant Method)

AIM:

Write a program to find out the root of equation x3 5x 7 = 0 and roots lie between [2,3] by
using Non-Linear Algebraic Method (Secant Method) using Scilab.

PRE-REQUISITE:

3. Knowledge of SCILAB Programming

4. Knowledge of Numerical Methods

THEORY:

Secant method is an iterative tool of mathematics and numerical methods to find the
approximate root of polynomial equations. During the course of iteration, this method
assumes the function to be approximately linear in the region of interest. Although secant
method was developed independently, it is often considered to be a finite difference
approximation of Newtons method.

Mathematical Derivation of Secant Method:

Consider a curve f(x) = 0 as shown in the figure below:

5
Secant method estimates the point of intersection of the curve and the X- axis (i.e. root of the
equation that represents the curve) as exactly as possible. For that, it uses succession of roots
of secant line in the curve.

Assume x0 and x1 to be the initial guess values, and construct a secant line to the curve
through (x0, f(x0)) and (x1, f(x1)). The equation of this secant line is given by:

If x be the root of the given equation, it must satisfy: f(x) = 0 or y= 0. Substituting y = 0 in

the above equation, and solving for x, we get:

Now, considering this new x as x2, and repeating the same process for x2, x3, x4, . . . . we end
up with the following expressions:

This is the required formula which will also be used in the program for secant method in
Scilab.

Secant method falls under open bracket type. The programming effort may be a tedious to
some extent, but the secant method algorithm and flowchart is easy to understand and use for
coding in any high level programming language.

6
ALGORITHM FOR BISECTION METHOD:

Given a continuous function f x Start

1. Get values of x0, x1 and e

Here x0 and x1 are the two initial guesses
e is the stopping criteria, absolute error or the desired degree of accuracy
2. Compute f(x0) and f(x1)
3. Compute x2 = [x0*f(x1) x1*f(x0)] / [f(x1) f(x0)]
4. Test for accuracy of x2
If [ (x2 x1)/x2 ] > e, *Here [ ] is used as modulus sign*
then assign x0 = x1 and x1 = x2
goto step 4
Else,
goto step 6
5. Display the required root as x2.
6. Stop

PROGRAM IN SCILAB
function [ ] = bisectionMethod( a,b )
a=3;
b=5;
MaxIt=1000;
epsilon=10^-6;
c=(b+a)/2;
NumIt=0;
while NumIt<MaxIt && abs (f(c))>epsilon
if f(a)*f(c)<0
b=c;
else
a=c;
end
NumIt=NumIt+1;
c=(b+a)/2;
end

7
c
end
function y=f(x)
y=x^3-1.8*x^2-10*x+17;
end

OUTPUT
C = 1.6618

RESULTS:
The root of the given Non-linear Algebraic Equation by Bisection Method is: 1.6618

CONCLUSION:
Thus we successfully studied the root finding program of Non-linear Algebraic Equation by
Bisection Method.

INDUSTRIAL APPLICATIONS:
1. The secant method is a root-finding algorithm that uses a succession of roots of secant
lines to better approximate a root of a function f.
2. The secant method does not require that the root remain bracketed like the bisection
method does, and hence it does not always converge.
3. The secant method has also been used for capturing the peak response of complex
multi-story buildings

VIVA QUESTIONS:
1. How does Secant method differ from Bisection method?
2. How does Secant method differ from False position method?
3. How is Secant method a form of Newton Raphson method?
4. Define Order of Convergence.
5. Define Broyden's method in relation to secant method.

8
 EXPERIMENT NO. 3

NAME OF THE EXPERIMENT

Non-linear Algebraic Equations (Newton Raphson Method)

AIM:

Write a program to find out the root of equation x3 - 5x + 3 = 0 and initial guess of zero by
using Non-Linear Algebraic Method (Newton Raphson Method) using Scilab.

PRE-REQUISITE:

5. Knowledge of SCILAB Programming

6. Knowledge of Numerical Methods

THEORY:

The idea of the method is as follows: one starts with an initial guess which is reasonably
close to the true root, then the function is approximated by its tangent line (which can be
computed using the tools of calculus), and one computes the x-intercept of this tangent line
(which is easily done with elementary algebra). This x-intercept will typically be a better
approximation to the function's root than the original guess, and the method can be iterated.

Suppose : [a, b] R is a differentiable function defined on the interval [a, b] with values in
the real numbers R. The formula for converging on the root can be easily derived. Suppose
we have some current approximation xn. Then we can derive the formula for a better
approximation, xn+1 by referring to the diagram on the right. The equation of the tangent line
to the curve y= (x) at the point x = xn is:
y = f (xn) (x-xn) + f (xn)
Where, f denotes the derivative of the function f.

The x-intercept of this line (the value of x such that y = 0) is then used as the next
approximation to the root, xn+1. In other words, setting y to zero and x to xn+1 gives:
0 = f (xn) (xn+1 - xn) + f (xn)
Solving for xn+1 gives:

xn+1 = xn [f(xn) / (f(xn)]

We start the process off with some arbitrary initial value x0. (The closer to the zero, the better.
But, in the absence of any intuition about where the zero might lie, a "guess and check"

9
method might narrow the possibilities to a reasonably small interval by appealing to the
intermediate value theorem.) The method will usually converge, provided this initial guess is
close enough to the unknown zero, and that '(x0) 0. Furthermore, for a zero of multiplicity
1, the convergence is at least quadratic (see rate of convergence) in a neighbourhood of the
zero, which intuitively means that the number of correct digits roughly at least doubles in
every step. More details can be found in the analysis section below.

The Householder's methods are similar but have higher order for even faster convergence.
However, the extra computations required for each step can slow down the overall
performance relative to Newton's method, particularly if f or its derivatives are
computationally expensive to evaluate.

ALGORITHM FOR NEWTON RAPHSON METHOD:

1. Start
2. Read x, e, n, d
x is the initial guess
e is the absolute error i.e. the desired degree of accuracy
n is for operating loop
d is for checking slope*
3. Do for i =1 to n in step of 2
4. f = f(x)
5. f1 = f'(x)

10
 6. If ( [f1] < d), then display too small slope and goto 11.
*[ ] is used as modulus sign*
7. x1 = x f/f1
8. If ( [(x1 x)/x1] < e ), the display the root as x1 and goto 11.
*[ ] is used as modulus sign*
9. x = x1 and end loop
10. Display method does not converge due to oscillation.
11. Stop

PROGRAM IN SCILAB

f=@(x) x^3-5*x+3;
df=@(x) 3*(x^2)-5;
a=0; b=1;
x=a;
for i=1:1:100
x1=x-(f(x)/df(x));
x=x1;
end
sol=x;
fprintf('\n Approximate Root is %.15f',sol)
a=0; b=1;
x=a;
er(5)=0;
for i=1:1:5
x1=x-(f(x)/df(x));
x=x1;
er(i)=(sol-x1)/sol;
end
plot(er)
xlabel('Number of iterations')
ylabel('Error')
title('Error Vs. Number of iterations')

11
OUTPUT
Approximate Root is 0.656620431047110

RESULTS:
The root of the given Non-linear Algebraic Equation by Newton Raphson Method is: ______

CONCLUSION:
Thus we successfully studied the root finding program of Non-linear Algebraic Equation by
Newton Raphson Method.

INDUSTRIAL APPLICATIONS:
1. These methods are used to approximate the solution of differential equations.
2. Bisection method and Newton-Raphson methods are used to find the roots and fixed
points of equations.
3. Civil Engineering: "Solution of the characteristic equation in the dynamic analysis of
a plane frame using the newtonraphson method".
4. Newton method for nonlinear systems: Application in Electronic Engineering.
5. Newtons method is a basic tool in numerical analysis and numerous applications,
including operations research and data mining.

12
VIVA QUESTIONS:
1. Newton Raphson method is used to solve linear or nonlinear equations?
2. Give example of Newton Raphson method being used in Optimization.
3. What do you mean by fixed point of an equation?
4. What is Fourier transforms of nonlinear differential equations?
5. Write down Colebrook equation. What is it used for?

13
 EXPERIMENT NO. 4

NAME OF THE EXPERIMENT

Non-linear Algebraic Equations (False Position Method/ Regula Falsi Method)

AIM:

Write a program to find out the root of equation sin x x 2 1 (roots lie between [0, 1]) by

using False Position Method with the help of SCILAB.

PRE-REQUISITE:

7. Knowledge of SCILAB Programming

8. Knowledge of Numerical Methods

THEORY

In numerical analysis, double false position became a root-finding algorithm that combines
features from the bisection method and the secant method.

Like the bisection method, the false position method starts with two points a0 and b0 such
that f a0 and f b0 are of opposite signs, which implies by the intermediate value
theorem that the function f has a root in the interval a0 , b0 , assuming continuity of the
function f. The method proceeds by producing a sequence of shrinking intervals ak , bk that
all contain a root of f.
At iteration number k, the number

14
 ck bk f (bk )
bk ak
f bk f ak
is computed. As explained below, ck is the root of the secant line through ak , f ak and
b , f b .
k k If f ak and f ck have the same sign, then we set ak 1 ck and bk 1 bk ,
otherwise we set ak 1 ak and bk 1 ck . This process is repeated until the root is approximated
sufficiently well.
The above formula is also used in the secant method, but the secant method always retains
the last two computed points, while the false position method retains two points which
certainly bracket a root. On the other hand, the only difference between the false position
method and the bisection method is that the latter uses ck ak bk / 2 .

f bk f ak
y f bk x bk
bk ak
We now choose ck to be the root of this line (substituting for x), and setting y 0 and see that

f bk f ak
f bk ck bk 0
bk ak
Solving this equation gives the above equation for ck .

ALGORITHM:

1. Start
2. Read values of x0, x1 and e
*Here x0 and x1 are the two initial guesses
e is the degree of accuracy or the absolute error i.e. the stopping criteria*
3. Computer function values f(x0) and f(x1)
4. Check whether the product of f(x0) and f(x1) is negative or not.
If it is positive take another initial guesses.
If it is negative then goto step 5.
5. Determine:
x = [x0*f(x1) x1*f(x0)] / (f(x1) f(x0))
6. Check whether the product of f(x1) and f(x) is negative or not.
If it is negative, then assign x0 = x;
If it is positive, assign x1 = x;
7. Check whether the value of f(x) is greater than 0.00001 or not.
If yes, goto step 5.
If no, goto step 8.
*Here the value 0.00001 is the desired degree of accuracy, and hence the stopping
criteria.*
8. Display the root as x.
9. Stop

15
PROGRAM IN SCILAB

format long;
f=@(x) sin(x)+x/2-1;
a=0;b=1;
stepsize =inf;
estep=1e-15;
n=0;
while (stepsize>=estep)
c=b-(f(b)*(b-a))/(f(b)-f(a));
if f(a)*f(c)<0
stepsize= b-c;
b=c;
else stepsize=c-a;
a=c;
end
n=n+1;
end
fprintf ('Approximate root after %d iterations is
%0.15f',n,c);

OUTPUT
Approximate root after 17 iterations is 0.704576912921746

RESULT:
The root of the given Non-linear Algebraic Equation by False Position Method is:
0.704576912921746 after 17 iterations.

CONCLUSION:
Thus we successfully studied the root finding program of Non-linear Algebraic Equation by
False Position Method.

16
INDUSTRIAL APPLICATIONS:
1. Modification of the Regula Falsi method to accelerate system convergence in the
prediction of trace quantities of atmospheric pollutants.
2. The Regula Falsi method is applied to prediction of trace quantities of air pollutants
produced by combustion reactions such as those found in industrial point sources.
3. The application of various classical root-finding methods is in digital maximum
power point tracking.

VIVA QUESTIONS:
1. Find root of following equation using Regula Falsi method: x + x/4 = 15
2. Explain the Two point Bracketing Method.
3. Explain the successive under-interpolation (SUI).
4. Explain the Illinois Algorithm.
5. Explain the Anderson & Bjrk's algorithm.

17
 EXPERIMENT NO. 5

NAME OF THE EXPERIMENT

Non-linear Algebraic Equations (Successive Substitution Method)

AIM:

PRE-REQUISITE:

9. Knowledge of Scilab Programming

10. Knowledge of Numerical Methods

THEORY:

In modular arithmetic, the method of successive substitution is a method of solving problems

of simultaneous congruence by using the definition of the congruence equation. It is
commonly applied in cases where the conditions of the Chinese remainder theorem are not
satisfied.
There is also an unrelated numerical-analysis method of successive substitution,
a randomized algorithm used for root finding, an application of fixed-point iteration.
The method of successive substitution is also known as back substitution.
Consider the following sets of liner equations:

a11x1 + a12x2 + a13x3 + a14x4 . . . . . . . . . . .. . . . + a1nxn = b1

a21x1 + a22x2 + a23x3 + a24x4 . . . . . . . . . . .. . . . + a2nxn = b2

a31x1 + a32x2 + a33x3 + a34x4 . . . . . . . . . . .. . . . + a3nxn = b3

.. ..... ...... ........ ........ ...................

.. ..... ...... ........ ........ ...................

an1x1 + an2x2 + an3x3 + an4x4 . . . . . . . . . . .. . . . + annxn = bn

Express these equations in the form: Ax = b where,

18
Now, decompose matrix A into Diagonal component matrix D, Lower triangular matrix L,
and upper triangular matrix U, such that:

A = D + L + U where,

Rewrite the system linear equation as:

( D + L ) x = b [ U + ( 1 ) D ] x

Where, is the relaxation factor and > 1

During the iteration process, the left hand side of the equation is solved by using the previous
value for x on right hand side.

The iteration equation is expressed analytically as:

x(k+1) = ( D + L ) -1 (b [U + ( 1 ) D ] x (k) ) = L (k) + c

Where, x(k) is the kth approximation for x and x(k+1) is (k+1)th approximation for x.

Taking the advantage of triangular matrix form of ( D + L), the (k+1)th can be evaluated
sequentially by using forward substitution. The expression obtained for x(k+1) is:

ALGORITHM:

write the first equation in its equivalent form

substitute it into the next
simplify, use the modular multiplicative inverse if necessary
continue until the last equation
back substitute, then simplify
rewrite back in the congruence form

19
PROGRAM IN SCILAB:

function[]=substitution(x)
x=0.5;
for i=1:100
x1=((x^3+(0.61506*x)-0.032823)/1.97844)^(1/2);
error=x1-x;
if abs(error)<0.5*10^-5
char 'iteration,root,relative error'
i,x1,error;
break
end
x=x1;
end

OUTPUT

ans =

iteration,root,relative error

i = 32

x1 = 0.3021

error = -4.9255e-006

RESULTS:
The root of the given Non-linear Algebraic Equation by Successive Substitution Method is:
0.3021 after 32 iterations.

CONCLUSION:
Thus we successfully studied the root finding program of Non-linear Algebraic Equation by
using Successive Substitution Method.

APPLICATIONS:
1. Used to solve Equations of State.
2. It is commonly applied in cases where the conditions of the Chinese remainder
theorem are not satisfied.
3. Numerical method for solving a nonlinear equation for the unknown.

20
VIVA QUESTIONS:
1. Successive substitution method is used to solve linear or nonlinear equations?
2. Explain the Successive over Relaxation method.
3. What do you mean by simultaneous methods?
4. Find the positive root of x3-2x-8=0 by method of successive substitution correct up to
two places of decimal.
5. Why is initial approximation necessary?

21
 EXPERIMENT NO. 6

NAME OF THE EXPERIMENT

Simulation of batch reactor
AIM:
To simulate the performance of 1) Variable pressure (Constant Volume) Batch Reactor, and
2) Variable volume (Constant Pressure) Batch Reactor using SimzLab software.
THEORY:
There are two types of batch reactors: 1) Variable pressure (Constant Volume) Batch Reactor,
and 2) Variable volume (Constant Pressure) Batch Reactor.
Variable pressure (Constant Volume) Batch Reactor is a constant density reaction system.
Most liquid-phase reactions as well as all gas phase reactions occurring in a constant volume
bomb fall in this class
In a constant volume system, the measure of reaction rate of component i becomes:

1 dNi d Ni / V dCi
ri (1)
V dt dt dt
P
Or for ideal gases, where C
RT
1 dPi
ri (2)
RT dt
For constant volume batch reactor the following equation:

a
PA C A RT PA0 P P0 (3)
n
Gives the concentration or partial pressure of reactant A as a function of the total pressure P
at time t, initial partial pressure of A, PA0 and initial total pressure of the system, P0, and for
the reaction:
aA bB .... rR sS ....
We have n r s .... a b .... (4)
Similarly for any product R we can find:
r
PR CR RT PR 0 P P0 (5)
n
Equations (4) and (5) are the desired relationships between total pressure of the system and
the partial pressure of reacting materials and products.

22
Let us introduce one other useful term, the fractional conversion, or the fraction of any
reactant, say A, converted to something else, or the fraction of A reacted away. We call this,
simply, the conversion of A, with symbol XA.
Suppose that NA0 is the initial amount of A in the reactor at time t = 0, and that NA is the
amount present at time t. Then the conversion of A in the constant volume system is given by

CA
X A 1 (6)
C A0

Differentiating we get,

dC A
dX A (7)
C A0

For first order irreversible unimolecular type reaction A B

dC A
rA kC A (8)
dt
For above relation, rearranging and integrating we obtain

CA
CA t
dC A
CA 0
CA 0
k dt OR ln
C A0
kt

(9)

Variable volume batch reactors used for isothermal constant pressure operations, or for
reactions having a single stoichiometry. For such systems the volume is linearly related to the
conversion, or

V V0
V V0 1 A X A or X A (10)
V0 A

Differentiating we get
dV
dX A (11)
V0 A

Here A is the fractional change in volume of the system between no conversion and
complete conversion of reactant A. Thus

VX A 1 VX A 0
A (12)
VX A 0

For variable volume system,

CA 1 X A
(13)
C A0 1 A X A

23
The rate of reaction (disappearance of component A) can be written for a variable volume
system as:

C A0 dX A
rA (14)
1 A dt
Or in terms of volume:

C dV C A0 d ln V
rA A0 (15)
VA dt A dt
SIMULATION USING SIMZLAB
1) Variable Pressure Batch Reactor
For the first order reaction A 2B carried out in a variable pressure (constant volume) batch
reactor the following data are given:
Reaction rate constant, k = 2x10-3 s-1 (value at 300 K)
Activation energy, Ea = 40.0 kJ/mol
Initial pressure = 100.0 kPa
Reaction temperature = 300 K
Reaction time = 500 s
Select the pressure range from 0 200 kPa and simulate the variable pressure (constant
volume) batch reactor using SimzLab by clicking on the Run Reaction button.
Generate the pressure vs. time and the fractional conversion and the final pressure after a
reaction time of 500 seconds.
2) Variable Volume Batch Reactor
For the first order reaction 2A B carried out in a variable volume (constant pressure) batch
reactor the following data are given
Reaction rate constant, k = 0.100 s-1(value at 300 K)
Activation energy, Ea = 20.0 kJ/mol
Initial volume = 5 x 10-3 m3
Reaction temperature = 300 K
Reaction time = 40 s
Simulate the variable volume (constant pressure) batch reactor using SimzLab by clicking on
the Run Reaction button.
Generate the volume vs. time, the concentration of A vs. time and the fractional conversion
vs. time plots for the variable pressure (constant volume) batch reactor.

24
Determine the final volume, the final concentration of A and the fractional conversion after
reaction time of 40s.
OUTPUT
1) Variable Volume Batch Reactor
Table 1: Concentration, volume and conversion variation for Variable Volume (Constant
Pressure) batch reactor
Conct. Final 12.5 11.5 0.714 3.22E-3
Time (s) Conv.
(mol/m3) Volume 13.0 10.9 0.728 3.18E-3

0 40.1 0 5.00E-3 13.5 10.4 0.741 3.15E-3

0.500 38.1 4.88E-2 4.88E-3 14.0 9.89 0.753 3.12E-3

1.00 36.3 9.52E-2 4.76E-3 14.5 9.40 0.765 3.09E-3

1.50 34.5 0.139 4.65E-3 15.0 8.95 0.777 3.06E-3

2.00 32.8 0.181 4.55E-3 15.5 8.51 0.788 3.03E-3

2.50 31.2 0.221 4.45E-3 16.0 8.09 0.798 3.01E-3

3.00 29.7 0.259 4.35E-3 16.5 7.70 0.808 2.98E-3

3.50 28.3 0.295 4.26E-3 17.0 7.32 0.817 2.96E-3

4.00 26.9 0.330 4.18E-3 17.5 6.97 0.826 2.93E-3

4.50 25.6 0.362 4.09E-3 18.0 6.63 0.835 2.91E-3

5.00 24.3 0.394 4.02E-3 18.5 6.30 0.843 2.89E-3

5.50 23.1 0.423 3.94E-3 19.0 6.00 0.850 2.87E-3

6.00 22.0 0.451 3.87E-3 19.5 5.70 0.858 2.86E-3

6.50 20.9 0.478 3.81E-3 20.0 5.43 0.865 2.84E-3

7.00 19.9 0.503 3.74E-3 20.5 5.16 0.871 2.82E-3

7.50 18.9 0.528 3.68E-3 21.0 4.91 0.877 2.81E-3

8.00 18.0 0.551 3.62E-3 21.5 4.67 0.883 2.79E-3

8.50 17.1 0.573 3.57E-3 22.0 4.44 0.889 2.78E-3

9.00 16.3 0.593 3.52E-3 22.5 4.23 0.895 2.76E-3

9.50 15.5 0.613 3.47E-3 23.0 4.02 0.900 2.75E-3

10.0 14.8 0.632 3.42E-3 23.5 3.82 0.905 2.74E-3

10.5 14.0 0.650 3.38E-3 24.0 3.64 0.909 2.73E-3

11.0 13.4 0.667 3.33E-3 24.5 3.46 0.914 2.72E-3

11.5 12.7 0.683 3.29E-3 25.0 3.29 0.918 2.71E-3

12.0 12.1 0.699 3.25E-3 25.5 3.13 0.922 2.70E-3

25
26.0 2.98 0.926 2.69E-3 33.5 1.41 0.965 2.59E-3

26.5 2.83 0.929 2.68E-3 34.0 1.34 0.967 2.58E-3

27.0 2.69 0.933 2.67E-3 34.5 1.27 0.968 2.58E-3

27.5 2.56 0.936 2.66E-3 35.0 1.21 0.970 2.58E-3

28.0 2.44 0.939 2.65E-3 35.5 1.15 0.971 2.57E-3

28.5 2.32 0.942 2.65E-3 36.0 1.10 0.973 2.57E-3

29.0 2.21 0.945 2.64E-3 36.5 1.04 0.974 2.56E-3

29.5 2.10 0.948 2.63E-3 37.0 0.991 0.975 2.56E-3

30.0 2.00 0.950 2.62E-3 37.5 0.943 0.977 2.56E-3

30.5 1.90 0.953 2.62E-3 38.0 0.897 0.978 2.56E-3

31.0 1.81 0.955 2.61E-3 38.5 0.853 0.979 2.55E-3

31.5 1.72 0.957 2.61E-3 39.0 0.812 0.980 2.55E-3

32.0 1.63 0.959 2.60E-3 39.5 0.772 0.981 2.55E-3

32.5 1.55 0.961 2.60E-3 40.0 0.734 0.982 2.55E-3

33.0 1.48 0.963 2.59E-3

Figure 1: Volume vs Time plot for variable pressure reactor

26
 Figure 2: Concentration vs Time plot for variable pressure reactor

Figure 3: Conversion vs Time plot for variable pressure reactor

2) Variable Pressure Batch Reactor
Table 2: Conversion and pressure variation for Variable pressure (Constant volume) batch
reactor

27
Time (s) Conv. Pressure (kPa)
0 0 100.0

25.0 4.88E-2 104.9

50.0 9.52E-2 109.5

75.0 0.139 113.9

100.0 0.181 118.1

125.0 0.221 122.1

150.0 0.259 125.9

175.0 0.295 129.5

200.0 0.330 133.0

225.0 0.362 136.2

250.0 0.394 139.3

275.0 0.423 142.3

300.0 0.451 145.1

325.0 0.478 147.8

350.0 0.503 150.3

375.0 0.528 152.8

400.0 0.551 155.1

425.0 0.573 157.3

450.0 0.593 159.3

475.0 0.613 161.3

500 0.632 163.2

28
 Figure 4: Pressure vs Time plot for Variable pressure batch reactor

Figure 5: Conversion vs Time plot for Variable pressure batch reactor

RESULTS
1) Variable Volume (Constant Pressure) Batch Reactor
Final concentration after 40 s = 0.734 mol/m3
Fractional conversion after 40 s = 98.2 %
Final volume after 40 s = 2.55 x 10-3 m3
2) Variable Pressure (Constant Volume) Batch Reactor
Fractional conversion after 500 s = 63.2 %
Final pressure after 500 s = 163.2 kPa

29
CONCLUSION
The simulation of a variable pressure (constant volume) batch reactor and a variable volume
(constant pressure) batch reactor was successfully carried out using SimzLab software.

INDUSTRIAL APPLICATIONS:
1. SimzLab is a desktop app (Mac and Win) with web access for access to simulations of
interest to chemists and chemical engineers.
2. Acceleration of scale-up projects, from laboratory bench to pilot and full-scale plants
3. Safety studies through simulation of breakdowns and of strongly exothermic reactions
4. Evaluation of VOC emissions of production
5. Cost reduction through optimization of operating conditions
6. Perpetuation and diffusion of process knowledge through acquisition and storage of
data
7. Feasibility studies, in particular when considering an existing reactor for a new
process
8. Investment risk reduction through beforehand representation of new equipment
configurations

VIVA QUESTIONS:
1. How is a chemical reaction defined?

2. What is the level of detail in the reactor description?

3. Does BatchReactor take into account all type of reactions?

4. Do you have any examples of industrial processes created with BatchReactor?

5. What are the hardware requirements for BatchReactor?

30
 EXPERIMENT NO. 7
NAME OF THE EXPERIMENT:
Simulation of Plug Flow Reactor
AIM:
To simulate the performance of Variable Flow Rate Plug Flow Reactor using SimzLab and
validate the predicted SimzLab results with the results obtained by hand calculations.
THEORY:
In a plug flow reactor the composition of the fluid varies from point to point along a flow
path; consequently, the material balance for a reaction component must be made for a
differential element of volume dV and can be written as:

FA FA dFA rA dV (1)

Since dFA d FA0 1 xA FA0 dX A (2)

We obtain on replacement

FA0 dX A rA dV (3)

This, then is the equation which accounts for A in the differential section of reactor of
volume dV. For the reactor as a whole the expression must be integrated. Now FA0 is the feed
rate which is constant but the rate of reaction certainly depends on the concentration or
conversion of materials. Grouping the terms accordingly and integrating them between
appropriate limits, we obtain:
V X Af
dV dX A
0 FA0 r
0 A
(4)

Thus,
X Af
V V dX A

FA0 v0 C A0 C A0
r
0 A
(5)

Here is the space time.

X Af
dX A
i.e. C A0 r
0 A
(6)

Above equation is the generalised performance equation for a PFR.

For constant volume system where a first order reaction takes place, we have:

CA CA0 1 X A (7)

rA kCA (8)

31
rA kCA0 1 X A (9)

Substituting above equation in (6)

X Af
dX A
C A0
0
kC A0 1 X A
(10)

Integrating above equation we get,

k ln 1 X A (11)

Equation (10) is the equation for a constant volume PFR when a first order reaction takes
place.
For variable volume PFR where a first order reaction takes place, we have

V V0 1 A X A (12)

And

C A0 1 X A
CA (13)
1 A X A
Substituting above equation in (6)
X Af
1 A X A dX A
C A0
0
kC A0 1 X A
(14)

Solving above equation we get

1
k ln 1 A A X A (15)
1 X A

SIMULATION USING SIMZLAB:

Variable flow rate plug flow reactor
For the first order reaction A 4B carried out in a variable plug flow reactor the following
data are given:
Reaction rate constant, k = 1 s-1
Activation energy, Ea = 100 kJ/mol
Initial Flow Rate = 2 x 10-2 m3/s
Reaction Temperature = 300 K
Simulate the variable flow rate plug flow reactor using SimzLab by clicking on the Run
Reaction button.Generate the outlet flow rate vs reactor volume and the fractional conversion
vs. reactor volume plots for the variable flow rate plug flow reactor.

32
Determine the outlet flow rate and fractional conversion after a space time of 5 seconds.
OUTPUT
Table 1: Variation of conversion and outlet volumetric flow rate with changing volume
Vol Flow 2.50E-2 0.495 4.97E-2 5.10E-2 0.683 6.10E-2
Conv
(m3) (m3/s)
2.60E-2 0.505 5.03E-2 5.20E-2 0.689 6.13E-2
0 0 2.00E-2
2.70E-2 0.515 5.09E-2 5.30E-2 0.694 6.16E-2
1.00E-3 0.045 2.27E-2
2.80E-2 0.524 5.15E-2 5.40E-2 0.699 6.19E-2
2.00E-3 0.084 2.51E-2
2.90E-2 0.534 5.20E-2 5.50E-2 0.703 6.22E-2
3.00E-3 0.119 2.72E-2
3.00E-2 0.542 5.26E-2 5.60E-2 0.708 6.25E-2
4.00E-3 0.150 2.90E-2
3.10E-2 0.551 5.31E-2 5.70E-2 0.713 6.28E-2
5.00E-3 0.178 3.07E-2
3.20E-2 0.559 5.36E-2 5.80E-2 0.717 6.30E-2
6.00E-3 0.204 3.22E-2
3.30E-2 0.568 5.41E-2 5.90E-2 0.722 6.33E-2
7.00E-3 0.228 3.37E-2
3.40E-2 0.575 5.45E-2 6.00E-2 0.726 6.36E-2
8.00E-3 0.250 3.50E-2
3.50E-2 0.583 5.50E-2 6.10E-2 0.730 6.38E-2
9.00E-3 0.271 3.62E-2
3.60E-2 0.591 5.54E-2 6.20E-2 0.734 6.41E-2
1.00E-2 0.290 3.74E-2
3.70E-2 0.598 5.59E-2 6.30E-2 0.739 6.43E-2
1.10E-2 0.309 3.85E-2
3.80E-2 0.605 5.63E-2 6.40E-2 0.743 6.45E-2
1.20E-2 0.326 3.96E-2
3.90E-2 0.612 5.67E-2 6.50E-2 0.747 6.48E-2
1.30E-2 0.343 4.06E-2
4.00E-2 0.619 5.71E-2 6.60E-2 0.750 6.50E-2
1.40E-2 0.358 4.15E-2
4.10E-2 0.625 5.75E-2 6.70E-2 0.754 6.52E-2
1.50E-2 0.374 4.24E-2
4.20E-2 0.632 5.79E-2 6.80E-2 0.758 6.55E-2
1.60E-2 0.388 4.33E-2
4.30E-2 0.638 5.83E-2 6.90E-2 0.762 6.57E-2
1.70E-2 0.402 4.41E-2
4.40E-2 0.644 5.87E-2 7.00E-2 0.765 6.59E-2
1.80E-2 0.415 4.49E-2
4.50E-2 0.650 5.90E-2 7.10E-2 0.769 6.61E-2
1.90E-2 0.428 4.57E-2
4.60E-2 0.656 5.94E-2 7.20E-2 0.772 6.63E-2
2.00E-2 0.440 4.64E-2
4.70E-2 0.662 5.97E-2 7.30E-2 0.776 6.65E-2
2.10E-2 0.452 4.71E-2
4.80E-2 0.667 6.00E-2 7.40E-2 0.779 6.67E-2
2.20E-2 0.463 4.78E-2
4.90E-2 0.673 6.04E-2 7.50E-2 0.782 6.69E-2
2.30E-2 0.474 4.85E-2
5.00E-2 0.678 6.07E-2 7.60E-2 0.785 6.71E-2
2.40E-2 0.485 4.91E-2

33
7.70E-2 0.789 6.73E-2 8.50E-2 0.812 6.87E-2 9.30E-2 0.833 7.00E-2
7.80E-2 0.792 6.75E-2 8.60E-2 0.815 6.89E-2 9.40E-2 0.835 7.01E-2
7.90E-2 0.795 6.77E-2 8.70E-2 0.817 6.90E-2 9.50E-2 0.837 7.02E-2
8.00E-2 0.798 6.79E-2 8.80E-2 0.820 6.92E-2 9.60E-2 0.840 7.04E-2
8.10E-2 0.801 6.80E-2 8.90E-2 0.823 6.94E-2 9.70E-2 0.842 7.05E-2
8.20E-2 0.804 6.82E-2 9.00E-2 0.825 6.95E-2 9.80E-2 0.844 7.06E-2
8.30E-2 0.806 6.84E-2 9.10E-2 0.828 6.97E-2 9.90E-2 0.846 7.08E-2
8.40E-2 0.809 6.86E-2 9.20E-2 0.830 6.98E-2 0.100 0.848 7.09E-2

34
 Figure 1: Concentration vs volume plot for variable volume PFR

Figure 2: Outlet volumetric flow rate vs volume plot for variable volume PFR
RESULTS:
1) Variable Flow Rate Plug Flow Reactor using SimzLab:
Outlet flow rate after space time of 5 sec = 7.09 x 10-2 m3/s
Fractional conversion after space time of 5 sec = 0.848
2) Variable Flow Rate Plug Flow Reactor using hand calculations:
Outlet flow rate after space time of 5 sec = 0.07088 m3/s
Fractional conversion after space time of 5 sec = 0.848

CONCLUSION
The simulation of a variable volume plug flow reactor was successfully carried out using
SimzLab software. The values of conversion and outlet flow rate predicted from simulation
was validated against results obtained from hand calculations and were found to be in
complete agreement with each other.
INDUSTRIAL APPLICATIONS:

35
Plug flow reactors are used for some of the following applications:

1. Large-scale production
2. Fast reactions
3. Homogeneous or heterogeneous reactions
4. Continuous production
5. High-temperature reactions

VIVA QUESTIONS:
1. Name the softwares that can be used for simulation & modelling of a chemical
reactor.
2. What is difference between Batch, CSTR and Plug Flow reactor?
3. Differentiate between steady state and unsteady state reactors.
4. Give mass balance and energy balance reactions for a PFR.
5. What are the input parameters given while simulation of a PFR using SimzLab?

36
 EXPERIMENT NO. 8
NAME OF THE EXPERIMENT
Linear regression using Least Squares Method and SCILAB
AIM:
To use MATLAB for determining the activation energy E in the equation:

k Ae E / RT where T is measured in Kelvin.

The reaction rate constant for the decomposition of a particular acid has been determined at
various temperatures as given in the following table:
t in Celsius 50 70.1 89.4 101
k in hr-1 1.08 x 10-4 7.34 x 10-4 45.4 x 10-4 138 x 10-4

THEORY

Least Squares Method is a mathematical procedure for finding the best-fitting curve to
a given set of points by minimizing the sum of the squares of the offsets ("the residuals") of
the points from the curve. The sum of the squares of the offsets is used instead of the offset
absolute values because this allows the residuals to be treated as a continuous differentiable
quantity. However, because squares of the offsets are used, outlying points can have a
disproportionate effect on the fit, a property which may or may not be desirable depending on
the problem at hand.

Figure 1: Vertical and perpendicular offsets for Least Squares Method

In practice, the vertical offsets from a line (polynomial, surface, hyperplane, etc.) are
almost always minimized instead of the perpendicular offsets. This provides a fitting function
for the independent variable x that estimates y for a given x (most often what an
experimenter wants), allows uncertainties of the data points along the x and y -axes to be
incorporated simply, and also provides a much simpler analytic form for the fitting
parameters than would be obtained using a fit based on perpendicular offsets. In addition, the
fitting technique can be easily generalized from a best-fit line to a best-fit polynomial when
sums of vertical distances are used. In any case, for a reasonable number of noisy data points,
the difference between vertical and perpendicular fits is quite small.

37
 The linear least squares fitting technique is the simplest and most commonly applied
form of linear regression and provides a solution to the problem of finding the best fitting
straight line through a set of points. In fact, if the functional relationship between the two
quantities being graphed is known to within additive or multiplicative constants, it is common
practice to transform the data in such a way that the resulting line is a straight line, say by
plotting T vs. l instead of T vs. l in the case of analysing the period T of a pendulum as a
function of its length l . For this reason, standard forms for exponential, logarithmic, and
power laws are often explicitly computed. The formulas for linear least squares fitting were
independently derived by Gauss and Legendre.

For nonlinear least squares fitting to a number of unknown parameters, linear least
squares fitting may be applied iteratively to a linearized form of the function until
convergence is achieved. However, it is often also possible to linearize a nonlinear function at
the outset and still use linear methods for determining fit parameters without resorting to
iterative procedures. This approach does commonly violate the implicit assumption that the
distribution of errors is normal, but often still gives acceptable results using normal
equations, a pseudoinverse, etc. Depending on the type of fit and initial parameters chosen,
the nonlinear fit may have good or poor convergence properties. If uncertainties (in the most
general case, error ellipses) are given for the points, points can be weighted differently in
order to give the high-quality points more weight.

Vertical least squares fitting proceeds by finding the sum of the squares of the vertical
deviations R 2 of a set of n data points

(1)

from a function . Note that this procedure does not minimize the actual deviations from the
line (which would be measured perpendicular to the given function). In addition, although the
unsquared sum of distances might seem a more appropriate quantity to minimize, use of the
absolute value results in discontinuous derivatives which cannot be treated analytically. The
square deviations from each point are therefore summed, and the resulting residual is then
minimized to find the best fit line. This procedure results in outlying points being given
disproportionately large weighting.

The condition for R 2 to be a minimum is that

(2)

For i 1, 2,...., n . For a linear fit,

(3)

so

(4)

38
 (5)

(6)

These lead to the equations

(7)

(8)

PROGRAM IN SCILAB

>> t=[50,70.1,89.4,101]

t=

50.0000 70.1000 89.4000 101.0000

>> T=t+273.15

T=

323.1500 343.2500 362.5500 374.1500

>> T1=1./T

T1 =

0.0031 0.0029 0.0028 0.0027

39
>> k=[1.08E-4,7.34E-4,45.4E-4,138E-4]

k=

0.0001 0.0007 0.0045 0.0138

>> polyfit (T1,log(k),1)

OUTPUT
ans =

1.0e+004 *

-1.1459 0.0026

RESULTS:
The slope (-E/R) is found to be -11459 and the intercept log(A) is found to be 26. The
activation energy is found to be 95270.126 kJ.

CONCLUSION
The method of least squares was studied successfully using Scilab.

APPLICATIONS:
Least square method has following applications:
1. Model fitting
2. Multi-objective least squares
3. Control
4. Estimation
5. Statistics

40
VIVA QUESTIONS:
1. What is regression analysis? Define linear & nonlinear regression.
2. What is Activation energy?
3. What is the significance of Arrhenius equation?
4. Differentiate between vertical and perpendicular offsets with respect to least square
method.
5. Define best-fit line and best-fit polynomial.

41
 EXPERIMENT NO. 9

NAME OF THE EXPERIMENT

4th Order Runge Kutta Method using SCILAB
AIM:
To find the value of y art x 2 , given that h 0.5, x0 0, y 0 0.5 . The Ordinary Differential
Equation (ODE) given is
dy
x y
dx
THEORY:
The 4th order RungeKutta method is often referred as "RK4", "Classical RungeKutta
method" or simply as "The RungeKutta method".
Let an initial value problem be specified as follows.
dy
y f (t , y), y(t 0 ) y 0 (1)
dt
Here, y is an unknown function (scalar or vector) of time t which we would like to
approximate; we are told that y the rate at which y changes, is a function of t and of y itself.
At the initial time t 0 the corresponding y-value is y0 . The function f and the data t0 , y0 are
given.
Now pick a step-size h>0 and define

y n 1 y n
h
k1 2k 2 2k 3 k 4 (2)
6
t n 1 t n h (3)
For n = 0, 1, 2, 3,.., using
k1 hf t n , y n (4)

h k
k 2 hf t n , y n 1 (5)
2 2

h k
k 3 hf t n , y n 2 (6)
2 2
k 4 hf (t n h, y n k 3 ) (7)

Here y n 1 is the RK4 approximation of yt n 1 , and the next value y n 1 is determined by the
present value ( y n ) plus the weighted average of four increments, where each increment is the
product of the size of the interval, h, and an estimated slope specified by function f on the
right-hand side of the differential equation.

k 1 is the increment based on the slope at the beginning of the interval, using y , (Euler's
method) ;

42
 k1
k 2 is the increment based on the slope at the midpoint of the interval, using y ;
2
k
k 3 is again the increment based on the slope at the midpoint, but now using y 2 ;
2
k 4 is the increment based on the slope at the end of the interval, using y k 3 .
In averaging the four increments, greater weight is given to the increments at the midpoint.
If is independent of, so that the differential equation is equivalent to a simple integral, then
RK4 is Simpson's rule.
The RK4 method is a fourth-order method, meaning that the local truncation error is on the
order of , while the total accumulated error is order.

PROGRAM IN SCILAB

function rungekutta
h=0.5;
x=0;
y=0.5;
fprintf ('Step 0: x=%12.8f, y=%12.8f\n',x,y);
for i=1:4
k1 = h*f(x,y);
k2 = h*f(x+h/2, y+k1/2);
k3 = h*f(x+h/2, y+k2/2);
k4 = h*f(x+h, y+k3);
y = y + (k1+2*k2+2*k3+k4)/6;
x = x + h;
fprintf('Step %d: x = %6.4f, y = %18.15f\n', i, x, y);
end
function v = f(x,y)
v=x+y;

OUTPUT

Step 0: x 0.00000000, y = 0.50000000

Step 1: x = 0.5000, y = 0.972656250000000
Step 2: x = 1.0000, y = 2.076019287109375
Step 3: x = 1.5000, y = 4.219063043594360
Step 4: x = 2.0000, y = 8.075955485925078

RESULT
The value of y at x=2 using the Runge Kutta method in SCILAB was found to be
8.075955485925078

43
CONCLUSION
Runge Kutta method was successfully studied using SCILAB.

APPLICATIONS:
1. To solve the Schrodinger equation for hydrogen and positronium atoms
2. To solve nonlinear partial differential equations.
3. To solve the problem of initial and boundary conditions.
4. It is helpful in estimating the Force-Deformation behavior in the non-linear range in
the field of structural engineering.
5. Applications in stability of structures and earthquake engineering and structural
dynamics.

VIVA QUESTIONS:
1. Define The Euler and Improved Euler Methods.
2. What are Semi-implicit additive RK methods?
3. What is the difference between 2nd order and 4th order RK method?
4. Write down equations for the 4th order RK method.
5. What is the difference between Eulers method and RK method?

44
 EXPERIMENT NO. 10

NAME OF THE EXPERIMENT

Numerical Integration in SCILAB
AIM:
To compute the heat required to raise the temperature of 1kg of a given material from -100K
to 200K using the Trapezoidal rule, Simpsons 1/3rd rule, Simpsons 3/8th rule and then by
using the quadrature function in SCILAB
The given integral function is
T2

H m c T dT (kJ/kg)
T1

Here, c T 0.132 1.56 104 T 2.64 107 T 2 (T is in K)

THEORY
Numerical Integration is nothing but finding an approximate value of
b
I f x dx (1)
a

There are two different strategies to develop numerical integration formulae. One is similar to
what we have adopted in numerical differentiation. That is, we approximate a polynomial for
the given function and integrate that polynomial within the limits of the integration. This
restricts us to integrate a function known at discrete tabular points. If these points are
uniformly spaced then the corresponding integration formulae are called as Newton - Cotes
formulae for numerical integration. On the other hand if we know the function explicitly but
could not integrate in the usual means because of the nature of the function then we can use
the concept called quadrature rule to find an approximate value to the integration.
In this method the function is evaluated at some predetermined abscissa (nodal) points and
then these values are added after multiplying them with some weights which are again
predetermined, to find an approximate value to the given integral.
Newton Cotes formulae
If the function values are known at uniformly spaced-points then as we have already
discussed, by interpolation we can find a polynomial of degree n for the given n 1 data-
points, which approximates the function at non tabular points and passes through the tabular
points. However, if the degree of this polynomial is too high then the approximated integral
value is prone to large errors due to round-off and local irregularities present in the function.
To avoid these large errors, the interval of integration is divided into sub intervals and a
polynomial of a lower degree is obtained in these sub intervals and integrated in the end
points of the subinterval. Then the required result can be obtained by adding these values
obtained in each of these sub intervals. In the following discussion let us consider three
widely used Newton-Cotes formulae obtained from polynomials of degree one (Trapezoidal
rule), two (Simpsons 1/3rd rule) and three (Simpsons 3/8th rule), to approximate I.
Trapezoidal Rule:

45
If f x is a function which is known at discrete points x1 , x2 , x3 ,..., xn . The plot for f x is
divided into number of uniformly spaced trapezoids using straight lines to approximate the
area under the curve.

The integral of f x is then given as

b
h
f x dx 2 f x 2 f x 2 f x .....2 f x f x
a
0 1 2 n 1 n (2)

ba
Here, h (n is the number of equally spaced panels in which the domain [a,b] is
n
discretized into).
Simpsons rule:
Simpsons 1/3rd rule:

The Simpsons 1/3rd rule is similar to Trapezoidal rule, with the only difference that instead
of using straight lines, in this method we approximate each part of the curve using parabolas.
This method can be used only when the domain is discretized into even number of intervals.

46
The integral of function f x whose values at discrete points x1 , x2 , x3 ,..., xn are known then
the integral can be evaluated using Simpsons 1/3rd rule, formula for which is given as
follows:
b
h
f x dx 3 f x f x 2 f x f x .... 4 f x f x ....
a
0 n 2 4 1 3 (3)

ba
Here, h (n is the number of equally spaced panels in which the domain [a,b] is
n
discretized into).
Simpsons 3/8th Rule
The concept for this rule is similar to the 1/3rd rule with the only difference that this rule is
applicable only if the domain is discretized into number of intervals that are divisible by 3.

The integral of function f x whose values at discrete points x1 , x2 , x3 ,..., xn are known then
the integral can be evaluated using Simpsons 3/8th rule, formula for which is given as
follows:
b
3h
f x dx
a
f x0 f x n 3 f x1 f x2 f x4 .... 2 f x3 ....
8
(4)

ba
Here, h (n is the number of equally spaced panels in which the domain [a,b] is
n
discretized into).
Adaptive Quadrature Method
Adaptive Simpson's method, also called adaptive Simpson's rule, is a method of numerical
integration proposed by G.F. Kuncir in 1962. It uses an estimate of the error we get from
calculating a definite integral using Simpson's rule. If the error exceeds a user-specified
tolerance, the algorithm calls for subdividing the interval of integration in two and applying
adaptive Simpson's method to each subinterval in a recursive manner. The technique is
usually much more efficient than composite Simpson's rule since it uses fewer function
evaluations in places where the function is well-approximated by a cubic function.

Adaptive Quadrature in SCILAB

Quadrature is a numerical method used to find the area under the graph of a function, that is,
to compute a definite integral.
b
q f x dx (5)
a

q = quad(fun,a,b) tries to approximate the integral of function fun from a to b to

within an error of 1e-6 using recursive adaptive Simpson quadrature. fun is a function
handle. Limits a and b must be finite. The function y = fun(x) should accept a vector
argument x and return a vector result y, the integrand evaluated at each element of x.

47
q = quad(fun,a,b,tol) uses an absolute error tolerance tol instead of the default
which is 1.0e-6. Larger values of tol result in fewer function evaluations and faster
computation, but less accurate results.

SCILAB COMMANDS
>> c = @(t)(0.132+1.5e-4*t+2.64e-7*t.^2);
>> qs = quad (c,-100,200,1e-6)

OUTPUT
qs = 42.642

RESULTS
The value of H obtained by: Trapezoidal rule is 41.848 kJ/kg (dividing the domain in six
equal parts), Simpsons 1/3rd rule is 41.848 kJ/kg (dividing the domain in six equal parts),
Simpsons 3/8th rule is 41.848 kJ/kg (dividing the domain in six equal parts) and using quad
in SCILAB is 42.642 kJ/kg.

CONCLUSION
The Trapezoidal rule, Simpsons 1/3rd rule, Simpsons 3/8th rule and the adaptive quadrature
method in SCILAB for numerical integration were studied successfully.

APPLICATIONS:

1. The integrand f(x) may be known only at certain points, such as obtained by sampling.
Some embedded systems and other computer applications may need numerical
integration for this reason.
2. A formula for the integrand may be known, but it may be difficult or impossible to
find an antiderivative that is an elementary function. An example of such an
integrand is f(x) = exp(x2), the antiderivative of which (the error function, times a
constant) cannot be written in elementary form.
3. It may be possible to find an antiderivative symbolically, but it may be easier to
compute a numerical approximation than to compute the antiderivative. That may be

48
 the case if the antiderivative is given as an infinite series or product, or if its
evaluation requires a special function that is not available.

VIVA QUESTIONS:
1. Define extrapolation methods.
2. Define Newton Cotes formulae.
3. Define numerical quadrature.
4. Differentiate between Simpsons 1/3rd rule and 3/8th rule.
5. What is the need of numerical integration?

49