NUMERICAL

 METHODS  IN  
CHEMICAL  ENGINEERING  
 
CLL-­‐113  
 
Prof.  Jayati  Sarkar  (course  co-­‐ordinator)  
Feedback  
•  Email  ( jayati@chemical.iitd.ac.in)  
Books  

Text  books:  
•  S.C.  Chapra  &  R.P.  Canale,  "Numerical  Methods  for  Engineers  with  Personal  
         Computer  Applica@ons",  McGraw  Hill  Book  Company,  1985.  
•  R.L.  Burden  &  J.  D.  Faires,  "Numerical  Analysis".  
Reference  books:  
•  Atkinson,  K.E.,  "An  Introduc@on  to  Numerical  Analysis",  John  Wiley  &  Sons,  
1978.  
•  Gupta,  S.  K.,  "Numerical  Methods  for  Engineers,  New  Academic  Science,  2012.  
•   W.  H.  et  al.,  "Numerical  Recipes  in  C:  The  Art  of  Scien@fic  Compu@ng,  
             3rd  Edi@on,  Cambridge  University  Press,  2007.  
Journals:    h[ps://www.aps.org,  h[ps://www.sciencedirect.com/  
•  Interna@onal  Journal  for  Numerical  Methods  in  Engineering  
•  Physics  of  Fluid  
•  Interna@onal  Journal  of  Numerical  Methods  for  Heat  &  Fluid  Flow  
•  Journal  of  Non-­‐Newtonian  Fluid  Mechanics  
 
 
 
Policies  
•  Lectures  will  be  carried  out  in  asynchronous    
mode  and  ppts  with  video  lectures  will  be  
uploaded  in  Impartus  in  Moodle.  

•   After  every  5th  lecture  the  lecture  on  the  

stipulated  lecture  day  will  be  converted  into  a  
live  –session,  where  part  of  it  can  be  converted  
into  lectures,  part  of  it  as  help-­‐sessions  to  clear  
doubt  and/or  part  to  conduct  quizzes.  

•  Do  not  use  your  mobile  phones  in  live-­‐sessions.  

 Practicals  
• Practical    starts  from  next  week.  
• Days  of  practical  
• Reminder:  Please  bring  your  calculators  to  the  
practicals?  
• Practical  sheets  will  be  uploaded  on  Moodle.  
• Codes/Reports  to  be  submitted  on  Moodle.  
• All  announcements  related  to  the  course  will  be  via    
course  mailing  lists  
• For  coding  C/C++  language  will  be  only  admissible.  
Marking  scheme  

• Quiz – 10 %
• Minor – 20 %
• Major – 40 %
• Term Paper-15 %
• Practical Submissions-15 %
• All exams will be closed book
Course  Outline  

1.  INTRODUCTION,APPROXIMATION  AND  ROUND  OFF  

ERRORS    
2  LINEAR  ALGEBRAIC  EQUATIONS  
3.  ROOT  FINDING  AND  SOLUTION  OF  NON-­‐LINEAR  
EQUATIONS  
4.  FUNCTIONS,  INTERPOLATION,  APPROXIMATION,  
REGRESSION  
5.  NUMERICAL  INTEGRATION  AND  DIFFERENTIATION  
6.  ORDINARY  DIFFERENTIAL  EQUATIONS  
7.  PARTIAL  DIFFERENTIAL  EQUATIONS  
 Where?  

Numerical   What  
Why?  
is?  
Methods  

How?  
Where do we need numerical methods

In Engineering
Civil, Mechanical,
chemical in all
engineering
problems

Numerical methods play a very

important role. Dedicated software
like ANSYS/FLUENT/ABACUS
are used to solve computer aided
design/ computational fluid
dynamics problem.
Where do we need numerical methods

v  In life sciences and medicines

Differential modeling of cancer network
Travelling wave analysis of tumor-immune interaction with
immunotherapy
Population dynamics determining algae growth. Therapy
Planning and 3-D imaging

v  Social science

Correlating different aspects of growth and economy

v  Business and finance

Predicting stock market and stock trading
 Why do we need numerical methods  

General  Navier  Stokes  eg.  àAre  Coupled  2nd  order  

PDESàNo  analyzed  solu@on  
Numerical Methods/ Computational techniques
Use computers to solve problems by step-wise, repeated and
iterative solution methods, which otherwise would be
tedious or unsolvable by hand calculations.  
 ALGEBRIC EQUATIONS
!!! !! + !!" !! = !!
!!" !! + !!! !! = !!
(linear if !!" = !"#$%, !"!#$!%&!!"!!!" = !!" ! !! )
ODE (IVP/BVP)
!"
= ! !, !
!"
PDE
!"!,! + 2!"!,! + !!,! + ⋯ … … = 0::::::::!! = !! − !"
!!! !
! !
! !
! = ! ! … . . !"!ℎ!"#!!! .!! ! > 0 ∷ ℎ!"#$%&'()!!"#
!" !"
!" !!!
= ! ! … … ! = 0 ∷ !"#"$%&'(!!"#
!" !"
!!! !!! ! )!
+ = ! !, ! … (! < 0 ∷ !""#$%#&!!"
!" ! !" !
Example of Iterative Method
Henon Approximation !
(BABYLONIAN METHOD ~1750BC)

!!!
! [!] !
! = ! + [!]
!
Example  
!

! ! = ! [!] − ! = 0
! ! [!]
! !!! = ! [!] − !!!!!(!"#$%&!!"#ℎ!"#!!"#ℎ!")
!! ! [!]
[!] ! [!] !
! −! ! +! 1 [!] !
= ! [!] − = = (! + [!] )
2! [!] 2! [!] 2 !
Example of Iterative Method

S=3: Finding value of 3

[0]
Iteration Number ! = 0 x =0.1

! ! ! !
! [!] = !
! ! +!! = !
0.1 + !.! = 0.5× 30.1 =15.05

Iteration Number ! = 1 x[1]=15.05

! ! ! !
! [!] = ! ! ! + ! ! = ! 15.05 + !".!" = 0.5× 30.1 =7.6247

[2]
Iteration Number ! = 2 x = 7.6247
[!] ! ! ! ! !
! = !
! +!! = !
7.6247 + !.!"#$! = 0.5× 8.018 =4.009
 ! !
Henon Approximation ! ::::::! !!! = ! ![!] + ![!]

[!] ! [!] 3 ! [!] ! [!!!] =0.5*(! [!] + 3 ! [!] )

0 0.1 30 15.08
1 15.05 0.199 7.6247
2 7.6247 0.393 4.009
3 4.009 0.748 2.3787
4 2.3787 1.2612 1.81997
5 1.8199 1.6484 1.73205
6 1.7341 1.72993 1.73205
7 1.73205 1.732049 1.73205
8 1.73205 1.73205 1.73205
9 1.73205 1.73205 1.73205
#
 Accuracy and Precision
Accuracy: refers to how closely a computer or
measured value agrees with true value
Precision: refers to how closely individual computed or
measured values are to each other

(a) & (c) though tightly grouped are

both inaccurate or biased
systematic deviation from truth.

(b) & (d) are equally accurate

(centered around bull’s eye), the
latter is more precise because the
shots are tightly grouped.

Numerical methods should be sufficiently accurate or

unbiased and should be accurately precise.
Error Definitions
Numerical errors arise from the use of approximation to
represent exact mathematical operations and quantities

Truncation error Round off errors

which results when which results when
approximations are numbers having limited
used to represent exact significant figures are

m a t h e m a t i c a l used to represent exact
procedures numbers

True Value= Approximation + Error
ET = True Value – Approximation
 Shortcomings
True Value - It takes no account of the order of
magnitude of the value under consideration

.
 Shortcomings
Approximation - True or analytical values are rarely
available

Error can be normalized by using the best

available estimate of the true value to the
approximation itself  

!"##$%&!!""#$%&'!(&$) − !"#$%&'(!!""#$%&'!(&$)!
=! ×!""%
!"##$%&!!""#$%&'!(&$)
For number of terms =2 !! = !! ! − !!! = 2.718282 − 2 =0.718282

!! ! !!!! !.!"#$#$!!
For number of terms =2 !! = ! ×100 = ×100=26.424116
!! !.!"#$#$

! !
!! !!! !!!
For number of terms =2 !! = ! ×100 = ×100=50
!! !

! !
For number of terms =3 !! = !! − !! = 2.718282 − 2.5 =0.218282

! !
!! !!! !.!"#$#$!!.!
For number of terms =3 !! = ! ×100 = ×100=8.030
!! !.!"#$#$
! !
!! !!! !.!!!
For number of terms =3 !! = ! ×100 = !.!
×100=20
!!
 (1)
et =2.718282+

Terms !"#$%& !! (100%) !! (100%)

1 1 63.2120582
2 2 26.42411641 100
3 2.5 8.030145511 40
4 2.666666667 1.898821878 12.5
5 2.708333333 0.36599097 3.076923077
6 2.716666667 0.059424789 0.613496933
7 2.718055556 0.008330425 0.102197241
8 2.718253968 0.00103123 0.01459854
&
 ! !
Henon Approximation ! ::::::! !!! = ! ![!] + ![!]

[!] ! [!] 3 ! [!] ! [!!!] =0.5*(! [!] + !! [!!!] !! [!!!]

3 ! [!] ) ! [!!!] − 3 ! [!!!] − ! [!]
= =
3 ! [!!!]
×100% ×100%
0 0.1 30 15.08 770.6442059 99.33687003
1 15.05 0.199 7.6247 340.2122597 97.38481514
2 7.6247 0.393 4.009 131.4597229 90.18957346
3 4.009 0.748 2.3787 37.33430853 68.53743641
4 2.3787 1.2612 1.81997 5.076016942 30.69995659
5 1.8199 1.6484 1.73205 4.6625E-05 5.07202448
6 1.7341 1.72993 1.73205 4.6625E-05 0.11835686
7 1.73205 1.732049 1.73205 4.6625E-05 0
8 1.73205 1.73205 1.73205 4.6625E-05 0
9 1.73205 1.73205 1.73205 4.6625E-05 0
#
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING  

CLL-­‐113  
 
Prof.  Jaya*  Sarkar  
 Round off errors:

•  Originate from the fact that computers can retain only a

fixed number of significant figures during calculation.

•  Irrational numbers like √7, e, π etc. cannot be expressed by

a fixed number of significant figures.

•  Also computers use a base 2 representation so they cannot

precisely represent base 10 numbers

•  The discrepancy introduced by omission of significant

figures.
Computer  representa*on  of  numbers  
Fundamental  unit  whereby  informa3on  is  stored  consists  of  a  string  of  binary  
digits:  bits  
Base  10  systems  uses  10  digits:-­‐  1,  2,3  4,  5,  6,  7,  8,  9,0  to  represent  numbers  
Place  and  face  values  are  mul3plied  to  get  the  number  
86409!
8×10! + 6×10! + 4×10! + 0×10! + 9×10!
•  Computers represent numbers in binary base 2 system

Integer representation in a computer

•  Signed magnitude approach
First bit represents sign
0 for positive (+) number
1 for negative (-) integer
Integer representation in a computer
   For  16  bit  computer    
  0 0 0 0 0 0 0 0 1 0 1 0 1 1 0 1

             Sign        
Number  
Range  of  base  10  that  can  be  represented  on  a  16  bit  computer  
  0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
=  214  x  1  +  213  x  1  +  212  x  1  +  211  x  1  +  210  x  1+  29  x  1  +…..+21  x  1  +  20  x  1  
=  215  -­‐  1  
=  32767            range  is  from    -­‐  32767  to  32767  
Since  1000000000000000  ∼  is  redundant  as  +0  &  -­‐0  are  same  ,  so  one  more  is  
added  to  -­‐ve  integer  that  is  -­‐32768  to  32767.    
Integer  numbers  above  and  below  this  cannot  be  stored    
FLOATING POINT REPRESENTATION
Fractional values can be represented by floating point representation
m . be
where m = the mantissa, b = the base of the number system being used, and
e = the exponent
Fractional part is called mantissa or significand.
156.78 à 0.15678 x 103 For a machine that stores 7 bit
Sign of
1/34 = 0.029411765 number Magnitude
of exponent
↪︎for decimal system 10 base
0.0294 x 100
0.2941 × 10−1
​1/𝑏  ≤ m < 1
Magnitude
Sign of of mantissa
In 10 base system 0.1≤m <1 exponent
In 2 base system 0.5 ≤m <1
Smallest  Number  
 
 
  Magnitude of mantissa
= 1×2-­‐1+0×2-­‐2+0×2-­‐3
  = 0.5
   
Magnitude of exponent

= 1 × 21 + 20 = 3
 
 
 
Exponent:  1  ×  21  +  20  =  3  
Magnitude:  1  ×  2−1  +  0  ×  2−2  +  0  ×  2−3=  0.5  
Number:    0.5  x  2-­‐3  =  0.  0625  
 0111100  =  (1  ×  2−1  +  0  ×  2−2  +  0  ×  2−3)  ×  2−3  =  (0.06250)10  
0111101  =  (1  ×  2−1  +  0  ×  2−2  +  1  ×  2−3)  ×  2−3  =  (0.078125)10   ∆x1=  0.015625  
0111110  =  (1  ×  2−1  +  1  ×  2−2  +  0  ×  2−3)  ×  2−3  =  (0.093750)10  
0111111  =  (1  ×  2−1  +  1  ×  2−2  +  1  ×  2−3)  ×  2−3  =  (0.109375)10  

0110100  =  (1  ×  2−1  +  0  ×  2−2  +  0  ×  2−3)  ×  2−2  =  (0.125000)10   ∆x2=  0.03125  

0110101  =  (1  ×  2−1  +  0  ×  2−2  +  1  ×  2−3)  ×  2−2  =  (0.156250)10  
0110110  =  (1  ×  2−1  +  1  ×  2−2  +  0  ×  2−3)  ×  2−2  =  (0.187500)10  
0110111  =  (1  ×  2−1  +  1  ×  2−2  +  1  ×  2−3)  ×  2−2  =  (0.218750)10  

Max: 0011111 = (1 × 2−1 + 1 × 2−2 + 1 × 2−3) × 23 = 7

∆x1   ∆x2  

(0.06250)10   (0.125000)10  
There is a limited range of quantities that may be represented
•  Large +ve and -ve quantities cannot be represented → overflow error
•  In addition very small quantities near to zero cannot be represented →
underflow error
→There are Only a finite number of quantities that can be represented
within the range
∆𝑥  
→The degree of precision is limited
If the value of fractional number falls here, Value it takes is that of lower one
So, error is always biased → This is called “CHOPPING OFF”.
So,  maximum  Et  is  ∼(Δx)  
→ Alternative is to round off the error ∆𝑥  
Depending on the value it is represented as the nearest allowable number
So, maximum E is rounding off ∼(Δx/2)  
The  interval  between  the  numbers  ∆x  increases  as  the  number  grows  in  
magnitude  
→It  allows  floa3ng-­‐point  representa3on  to  preserve  significant  digits  
→Quan3sizing  errors  ∝  magnitude  of  number  it  represents  
         
           |∆x|/|x|      <  ε          Chopping  
           |∆x|/|x|      <    ε/2            Round  off  error  
               𝜀  à  machine  epsilon  
                       =  b(1-­‐t)  (no  of  significant  digit  in  man3ssa)  
→Though round off errors are important → testing convergence
•  No. of significant digits carried on most computers allows most engineering
problems to be performed with more than acceptable precision

Computer with IEEE format

∼24 bits from mantissa →10-38 to 1039
Double precision ∼15 to16 decimal digits of precision
∼10-308 to 10308
Common arithmetic operations
Addition
Ø  0.1557 x 101 + 0.4381 x 10-1
Ø  Hypothetical decimal computer with 4-digit mantissa and 1 digit exponent
 
The decimal of the mantissa of the smallest number is to be shifted by
(1-(-1)) = 2 places
0.1557 x 101
0.004381 x 101
0.160081 x 101
                             ↓
•  Chopped off- the last two digits shifted to the right is essentially lost
 More  Acute  when  adding  two  vastly  different  magnitude  number    (4000  &  0.001)
⇒ 0.4000 x 104
0.0000001 x 104

0.4000001 x 104 ⇒ 0.4000 x 104 (Chopped off)

↓
as if no addition has been done
Subtraction
28.86 from 36.41

⇒ 0.3641 x 102
0.2886 x 102
0.0995 x 102
0. 9950 x 101 (Addition of zero)
→ Subtracting two close numbers 0.7641 x 103 from 0.7642 x 103
⇒ 0.7642 x 103
0.7641 x 103
0.0001 x 103
0.1000 x 100
↓
Will appear as a significant digit in subsequent computation
Multiplication

0.1363 x 103 x 0.6423 x 10-1

Multiplication of two n digit mantissa will yield 2n digit results,
→ Most computers hold intermediate results in a double length register

Exponents are added, mantises are multiplied

0.08754549 x 102

⇒ 0.8754549 x 101

Chopped off
0.8754 x 102 ⇒ loosing precision
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING

CLL-­‐113  
Solution  of  Linear  Equations
Prof.  Jayati  Sarkar
 Scalar, Vector, Matrices
Scalar: T=500°C :Magnitude 3 (2,3)
2 Element of
Vector : Magnitude and direction (velocity, acceleration) 1

Represented as ordered set of scalar 1 2 3

A point in 2-dimension space

Norm

Rows Column
Matrix:

Liner operations addition / subtraction

Scalar/Matrix multiplication rules,
Eigen value/ vectors,
Rank of a matrix by Echelon method
 A Linear Equation
Coefficients
4  x  +5y+7z+3a=15
Variables
❖ 4 equations to solve this
❖ It resembles 3 dimensional hyper space::::::::::1 degree is lost

❖ Not guaranteed that will have a solution

❖ A n×n matrix will have a unique solution when

 
No. of rows : No. of Equations
No. of Columns: No. of unknowns
C
D=1×(-­‐4)-­‐(2)×(-­‐2)=0
D=1×(-­‐4)-­‐(2)×(-­‐2)=0
1.999

• In well conditioned system small change in RHS

1.0
caused small change in solutions
• In ill conditioned system small change in RHS
caused large change in solutions
λ1
λ2

2
1
Intersection of lines/planes/Hyper surfaces==➔concept becomes difficult as n increases

Linear Combination of two vectors is the resultant vector

Requirement to find the scalar that solve equation

 VECTOR SPACE
 

If V & W are linearly independent they form the basis

for that particular vector space

Rn  space  there  can  be  at  most  n  linearly  independent  vectors

 ➢As  base  vectors    
NOT  linearly  independent    
2V=m

If  the  point  is  outside  this  line  

➢ The  linear  combination  of  
these  two  vectors  wont  be  
able  to  represent  b  

In this case there are infinite combinations that give b Can only span this
line
Poorly conditioned matrix
 1

➢Direct methods
• Gauss Elimination, Gauss Jordon Elimination, Sparse method Thomas Algorithm
➢ Indirect methods/ Iterative Methods
• Jacobi, Gauss-seidel & relaxation methods
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING  

CLL-­‐113  
 
Solution  of  Linear  Equations  
Prof.  Jaya*  Sarkar  
General  𝑛×𝑛  system  

v  What  does  the  solu3on  signify?  

v  Row  format  
v  Column  format  
v  Unique  solu3on?  
Example  of  Elimina3on  
Are  opera3ons  on  original  equa3ons  jus3fied?  

L1 ≡ x1 − 2x2 = 1
L2 ≡ 2x1 + x2 = 7
 Gauss  Elimina3on  

xn=  
 Step  k=1   Pivoted  Element  
Pivo3ng  row  
k=1   #a11(0) a12(0)
......................................a1n(0) &
% (
% # (0) a21 (0)
(0)
& # (0) a21 (0) &
(0) (
%0 %a22 − (0) × a12 ( ............ %a2n − (0) × a1n ((
a11 a11
% $ ' $ '(
#b1(0) &
(1)
% # (0) a31 (0)
(0)
& # (0) a (0)
(0)
&( % (
31
A = %0 %a32 − (0) × a12 ( ............ %a3n − (0) × a1n (( (0)
$ a11 ' $ a11 '( %b(0) − a21 × b(0) (
% % 2 1 (
%. ( a11(0)
. ............ . % (
% ( % a (0)
(
% # (0) an1 (0) & # (0) an1 (0) &( b(1) = %b3(0) − 31 × b (0)

%0
(0)
%an2 − (0) × a12 ( ............ %ann − (0) × a1n ((
(0)
a11(0)
1 (
$ $ a11 ' $ a11 '' % (
%. . (
% (0) (
a
%b(0) − n1 × b(0) (
%$ n a11(0)
1
('
 At  the  end  of  Step  k=1   !b1(0) $
!a11(0) (0)
a12 (0) $
.............a1n # (1) &
# & #b2 &
(1) (1)
#0 a22 .............a2n & (1) # (1) &
# b = #b3 &
(1) (1) (1) &
A = #0 a32 .............a3n & #. &
#. . ............ . & # (1) &
# (1) (1)
& #"bn &%
#"0 an2 .............ann &%

k=1  
(0)
a
aij(1) = aij(0) − i1(0) × a1(0)j
a11
(0)
a
bi(1) = bi(0) − i1(0) × b1(0)
a11
 Step  k=2   Pivoted  Element  

Pivo3ng  row  
# &
  % (0) (
(0) (0)
a
k=2   % 11 a12 ......................................a1n (
%0 (1)
a22 (1)
......................................a2n (
% ( # &
% # (1) a32 (1) & # a (1) &(
(2)
A = %0 0 (1) (1) 32
%a33 − (1) × a23 ( ............ %a3n − (1) × a2n ((
(1) % (0) (
$ a22 ' $ a22 '( %b1 (
% %b2(1) (
%. . ............ . ( % (
% ( % (1) a32
(1)
(
% # (1) an2 (1)
(1)
& # (1) an2(1) &
(1) ( b = %b3 − (1) × b2(1)
(2)

%0 0 %an3 − (1) × a23 ( ............ %ann − (1) × a2n (( a22 (

$ $ a22 ' $ a22 '' % (
%. . (
% (1) (
(1) a
%bn − n2 × b2(1) (
(1)
%$ a22 ('
At  the  end  of  Step  k=2   !b1(0) $
!a11(0) (0)
a12 ..............a1n(0) $ # (1) &
# & #b2 &
(1) (1)
#0 a22 ..............a2n & (2) # (2) &
# b = #b3 &
(2) (2) (2) &
A = #0 0 a33 ...........a3n & #. &
#. . ............ . & # (2) &
# (2) (2)
& #"bn &%
#"0 0 an2 .........ann &%
(1)
a
k=2  
2   aij(2) = aij(1) − i2(1) × a2(1)j
3,   a22
(1)
a
3,   bi(2) = bi(1) − i2(1) × b2(1)
a22
For  general  Step  k  

k = 1, 2,....(n −1)

2  
k=k   (k+1),  

(k+1),   a (k−1)
aij(k ) = aij(k−1) − ik(k−1) × akj(k−1)
akk
(k−1)
a
bi(k ) = bi(k−1) − ik(k−1) × bk(k−1)
akk
 Back  Subs3tu3on  
At  the  end  of  Step  k=n-­‐1  
"a11(0) a12(0) ..............a1n(0) %
Back  subs3tu3on  
$ '
$0
(1)
a22 ..............a2n ' (1)
bn(n−1)
$ (2) (2) '
xn = (n−1)
(n−1) $ 0 0 a33 ...........a3n ' ann
A =
$. . ............ . ' n
$ (n−2) (n−2)
' bi(i−1) − ∑ aij(i−1) × x j
$0 0 0....an−1n−1an−1n '
j=i+1
$0
# 0 0..............ann (n−1) '
&
xi = (i−1)
a ii
"b1(0) %
$ (1) '
$b2 ' i = ( n −1), ( n − 2 ),...1
(n−1) $ (2) '
b = $b3 '
$. '
$ (n−1) '
$#bn '&
Solu3on  by  steps  of  Gauss  Elimina3on  

"3 − 0.1 0.2 % "3 − 0.1 0.2 % "3 − 0.1 0.2 %

$ ' (1) $ ' (2) $ '
(0)
A = $0.1 7 − 0.3' A = $0 7.00333 − 0.29333' A = $0 7.00333 − 0.29333'
$#0.3 − 0.2 10'& R2 = R2 + α R $0 − 0.19000 10.02000'& $#0 0 10.01204 '&
#
21 1

R3 = R3 + α 31R1 R3 = R3 + α 32 R2
0.1 −0.19
α 21 = − α 32 = −
3 7.00333
0.3
α 31 = −
3
"7.85 % "7.85 % "7.85 %
$ ' $ ' $ '
b(0) = $−19.3' ⇒ b(1) = $−19.56167' ⇒ b(2) = $−19.56167'
$#71.4 '& $#70.615 '& $#70.08929 '&
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING  

CLL-­‐113  
 
Solution  of  Linear  Equations  
Pivoting,  FLOPS,  GJ,  LU  Decomposition  
Prof.  Jaya*  Sarkar  
 Pivo%ng  
Ø  p  stores  current  pivot  rows  
Ø  Big  stores  current  pivot  element  
Ø  First  task  to  see  if  element  below  in  the  column  is  higher  than  the  pivot  
element  
Ø  New  largest  element  &  its  row  number  is  stores  in  big  &  p    
   
  If the previous search ends with 𝑝≠𝑘 then switching of row required
 
p=  k   IF  (𝑝≠𝑘)  
big=|ak,k|   FOR  (jj  =k,  n)    
FOR(  ii  =k+1,  n)   {  
{   dummy  =ap,jj  
ap,jj=ak,jj  
dummy  =|aii,k|  
 
IF  (dummy  >  big)   ak,jj=dummy  
big  =dummy   }  
p=  ii   dummy  =bp  
END  IF   bp=bk      
}   bk  =dummy  
END  IF  
 Pivo%ng  
Ø  p  stores  current  pivot  rows  
Ø  Big  stores  current  pivot  element  
Ø  First  task  to  see  if  element  below  in  the  column  is  higher  than  the  pivot  
element  
Ø  New  largest  element  &  its  row  number  is  stores  in  big  &  p    
   
  If the previous search ends with 𝑝≠𝑘 then switching of row required
 
p=  k   IF  (𝑝≠𝑘)  
big=|ak,k|   FOR  (jj  =k,  n)    
FOR(  ii  =k+1,  n)   {  
{   dummy  =ap,jj  
ap,jj=ak,jj  
dummy  =|aii,k|  
 
IF  (dummy  >  big)   ak,jj=dummy  
big  =dummy   }  
p=  ii   dummy  =bp  
END  IF   bp=bk      
}   bk  =dummy  
END  IF  
 FLOPS  involved  in  Gauss  Elimina%on    

1  
k  
 FLOPS  IN  STEP  k=1  for  Gauss  Elimina%on  
#a11(0)
a12(0)
......................................a1n(0) & #b1(0) &
% ( % (
% # (0) a21(0)
(0)
& # (0) a21(0) &
(0) (
(0)
%b(0) − a21 × b(0) (
% 0 % a 22 − × a 12 ( ............ % a 2n − × a1n ( (
$ a11(0)
' $ a11(0)
' % 2 a (0) 1 (
% ( % 11
(
% # (0) a31 (0) & # (0) a31 (0) &( (1) % (0) a31 (0)
(
A = %0 %a32 − (0) × a12 ( ............ %a3n − (0) × a1n (( b = %b3 − (0) × b1(0) (
(1) (0) (0)

% $ a11 ' $ a11 '( a11

% (
%. . ............ . ( %. . (
% ( % (
# (0) an1 (0) & # (0) & (0)
% (0) (0) a n1 (0)
( a
%b(0) − n1 × b(0) (
%0 %an2 − (0) × a12 ( ............ %ann − (0) × a1n (( % n (0) 1
('
$ $ a 11 ' $ a 11 '' $ a11

Calcula*on  per  row  

Division  
(0)
a21 1   Total  Number  of  Addi*on/Subtrac*on  
α 21 = (0)
a11 per  row=n-­‐1+1=n    
Subtrac*on   LHS=n-­‐1   RHS=1   Total  Number  of  Mul*plica*on/
  Division  
Mul*plica*on   Division   LHS=n-­‐1   RHS=1   per  row=1+n-­‐1+1=n+1  
  In  step  k=1  
n
No.  of  middle  loops:     ∑ 1 = n − 2 +1 = n −1 Total  Number  of  Addi*on/Subtrac*on=n×(n-­‐1)  
i=2 Total  Number  of  Mul*plica*on/Division=(n+1)×(n-­‐1)  
 Outer  Loop  (k)   Middle  Loop  (i)   Addi*on/ Mul*plica*on/
Subtrac*on     Division    
FLOPS   FLOPS  
 
1   2,n   (n-­‐1)×n   (n-­‐1)×(n+1)  
2   3,n   (n-­‐2)×(n-­‐1)   (n-­‐2)×(n)  
:   :   :   :  
k   k+1,n   (n-­‐k)×(n-­‐k+1)   (n-­‐k)×(n-­‐k+2)  
 
n-­‐1   n,n   (1)×(2)   (1)×(3)  
 

=  

n  

2  
 n
b (i−1)
(i-1)
i − ∑ aij(i−1) × x j
j=i+1 bn( n−1)
xi = (i−1)
, xn = ( n−1)
a ii ann
i   Back  Subs*tu*on   Mul*plica*on+   Addi*on/
  Division   Subtrac*on  

n-­‐1   xn−1 = ( bn−1

(n−2 ) (n−2 )
− an−1,n × xn ) an−1,n−1
(n−2 )
1+1   1  

n-­‐2   xn−2 = ( bn−2 × xn ) an−2,n−2 2+1   2  

(n−3) (n−3) (n−3) (n−3)
− an−2,n−1 × xn−1 − an−2,n

:   :   :  

n-­‐i   i+1   i  

:   :   :  

n-­‐(n-­‐1)   (n-­‐1)+1   (n-­‐1)  

1 1
Total   1+ ∑ (i +1) ∑i
 Total  number  of  FLOPS  in  Back  Subs%tu%on      
1 1 1
Mul%plica%on  +  Division   1+ ∑ (i +1) = 1+ ∑ 1 + ∑ i
i=n−1 i=n−1 i=n−1

n(n −1)
= 1+ ( n −1) +
2
n ( n +1)
=
2
1
Addi%on  +  Subtrac%on   ∑ (i ) =
n(n −1)
2
i=n−1

Total   n ( n +1) n(n −1)

= + = n 2 + O(n)
2 2
Elimina%on   Back  Subs%tu%on    
3
2n 3
Total  number  of  FLOPS  in  Gauss  Elimina%on   =
2n 2 2
+ n + O(n ) + O(n) = + n 2 + O(n 2 )
3 3
 3x1 − 0.1x2 + 0.2x3 = 7.85....(1)
Varia%ons  of  Gauss  Elimina%on:  Gauss  Jordon   0.1x1 + 7x2 − 0.3x3 = −19.3...(2)
0.3x1 − 0.2x2 +10x3 = 71.4...(3)

"3 − 0.1 - 0.2% "1 − 0.03333 - 0.06667% "1 − 0.03333 - 0.06667 %

$ ' (0) $ ' (1) $ '
(0)
A = $0.1 7 − 0.3 ' A * = $0.1 7 − 0.3 ' A = 0 7.00333 − 0.29333
$ '
$#0.3 − 0.2 10 '& $#0.3 − 0.2 10 '& R2 = R2 + α 21R1 $#0 − 0.19000 10.02000'&
Normalize  
Pivot  row   R3 = R3 + α 31R1
α 21 = −0.1 1
R1 = R1 / a11
α 31 = −0.3 1

"7.85 % "2.61667% "2.61667 %

(0) $ ' (0) $ ' (1) $ '
b = $−19.3' ⇒ b * = $−19.3 ' ⇒ b = $-19.56167'
$#71.4 '& $#71.4 '& $#70.615 '&
 "1 − 0.03333 - 0.06667 % "1 − 0.03333 - 0.06667 %
"1 0 − 0.0681%
(1) $ ' (1) $ '
A = $0 7.00333 − 0.29333' A* = $0 1.00000 − 0.0419 ' R1 = R1 + α12 R2 (2) $ '
A = $0 1 − 0.0419'
$#0 − 0.19000 10.02000'& Normalize   $#0 − 0.19000 10.02000'& R3 = R3 + α32 R2
$#0 0 10.0120'&
Pivot  row  
α12 = −(−0.0333 1)
R2 = R2 / a22
α 32 = −(−0.19 1)

!2.61667 $ !2.61667$ !2.5237 $

# & # & # &
b(1) = #-19.56167& ⇒ b *(1) = #-2.7932 & ⇒ b(2) = #-2.7932 &
#"70.615 &% #"70.615 &% #"70.0843&%

"1 0 − 0.0681% "1 0 − 0.0681% R1 = R1 + α13 R3 !1 0 0$

$ ' (2) $ ' R2 = R3 + α 23 R3 # &
A(2) = $0 1 − 0.0419' Normalize   A* = 0 1 − 0.0419 A(3) = #0 1 0&
$ ' α13 = −(−0.0681 1)
$#0 0 10.0120'& Pivot  row   $#0 0 1 '& #"0 0 1&%
α 23 = −(−0.0419 1)
R3 = R3 / a33
!2.5237 $ !2.5237 $ !3.0 $
# & # & # &
b(2) = #-2.7932 & ⇒ b *(2) = #-2.7932& ⇒ b(3) = #-2.5&
#"70.0843&% #"7.0 &% #"7.0 &%
 than  

-­‐1  

Steps  of  Gauss  Jordon  

[ A!I!b] ⇒ #$I!A−1 !x%& [ A] [ x ] ⇒ [b]

[ A!I ] ⇒ #$I!A−1 %& −1

"3 − 0.1 - 0.2!1.0 0.0 0.0%

[ x ] ⇒ [ A] [b]
$ '
$0.1 7 - 0.3!0.0 1.0 0.0 '
$#0.3 − 0.2 10!0.0 0.0 1.0 '&
 Varia%ons  of  Gauss  Elimina%on:  LU  Decomposi%on  
[ A] [ x ] ⇒ [b]
y1 = b1
[ L ] [U ] [ x ] ⇒ [b] y2 = b2 − α 21 y1
i−1

[y] yi = bi − ∑α ij y j
j=1
[ L ] [ y] ⇒ [b]
Forward  Subs%tu%on  
[U ] [ x ] = [ y] ⇒ [ x ]
"1 0 0 0........ %" y1 % "b1 %
$ '$ ' $ ' Back  Subs%tu%on  
$−α 21 1 0 0......... '$ y2 ' $b2 '
$−α 31 − α 32 1 0........ '$ y3 ' = $b3 '
$ '$ ' $ '
$−α 41 - α 42 - α 43 0........'$ y4 ' $b4 '
$#! ! ! ! '&$#! '& $#! '&Original  b  value  
 Varia%ons  of  Gauss  Elimina%on:  LU  Decomposi%on  
[ A] [ x ] ⇒ [b]
y1 = b1
[ L ] [U ] [ x ] ⇒ [b] y2 = b2 − α 21 y1
i−1

[y] yi = bi − ∑α ij y j
j=1
[ L ] [ y] ⇒ [b]
Forward  Subs%tu%on  
[U ] [ x ] = [ y] ⇒ [ x ]
"1 0 0 0........ %" y1 % "b1 %
$ '$ ' $ ' Back  Subs%tu%on  
$−α 21 1 0 0......... '$ y2 ' $b2 '
$−α 31 − α 32 1 0........ '$ y3 ' = $b3 '
$ '$ ' $ '
$−α 41 - α 42 - α 43 0........'$ y4 ' $b4 '
$#! ! ! ! '&$#! '& $#! '&Original  b  value  
 Varia%ons  of  Gauss  Elimina%on:  LU  Decomposi%on  
"3 − 0.1 - 0.2% "3 − 0.1 - 0.2 % "3 − 0.1 - 0.2 %
$ ' (1) $ ' (2) $ '
(0)
A = $0.1 7 − 0.3 ' A = $0 7.00333 − 0.29333' A = $0 7.00333 − 0.29333'
$#0.3 − 0.2 10 '& $#0 − 0.19000 10.02000'& $#0 0 10.01201 '&

"3 − 0.1 - 0.2 % (0)

a21
α 21 = − (0) = −
0.1
$ ' a11 3
U = $0 7.00333 − 0.29333' (0)
a31 0.3
α 31 = − (0) = −
a11 3
$#0 0 10.01201 '& (1)
a32 0.19
α 32 = − (1) =
a22 7.00333
"1 0 0 %
$ '
L = $−α 21 1 0 '
$#−α31 − α32 1 '&
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING  

CLL-­‐113  
 
Solution  of  Linear  Equations  
Thomas  Algorithm,  Iterative  Methods  (GS,  JI),  SOR  
Prof.  Jaya*  Sarkar  
 i

1 2 3 . . . i-1 i+1 … N-2 N-1 N

/K  

/K  
/K  
 -­‐ banded
diagonal   super  diagonal  

sub  diagonal   -­‐

Banded  Matrix  Structure::::Tri-­‐diagonal  Matrix  

Thomas  Algorithm  

"b1 c1 0 0 0 0 0 %" x1 % "d1 %

$ '$ ' $ '
0   $a2 b2 c2 0 0 0 0 '$ x2 ' $d2 '
$0 a3 b3 c3 0 0 0 '$ x3 ' $ d3 '
$ '$ ' $ '
$0 0 a4 b4 c4 0 0 '$ x4 ' = $d4 '
$! ! ! ! ! ! ! '$! ' $! '
$ '$ ' $ '
$0 0 0 0 aN−1 bN−1 cN−1 '$ x N−1 ' $d N−1 '
$#0 0 0 0 0 aN bN '&$# x N '& $#d '&
N
 Thomas  Algorithm  
1.  Normalize  pivot  row    
Step  k=1  
             with  pivot  element  
c1 d !1 c1% 0 0 0 0 0 # !d1% #
!b1 c1 0 0 0 0 0 # !d1 # c1! = , d1! = 1 ' (' (
& '& ' b1 b1
&a2 b2 c2 0 0 0 0 ' &d2 ' 'a2 b2 c2 0 0 0 0 ( 'd2 (
&0 a3 b3 c3 0 0 0 ' & d3 ' '0 a3 b3 c3 0 0 0 ( ' d3 (
&
!" A #$, !"d #$&0
(0) (0) '& ' !" A *(0) #$, !"d *(0) #$''0 0 a4 b4 c4 0 0
('
(, 'd4
(
(
0 a4 b4 c4 0 0 ', &d4 '
&! '! ! ! ! ! ! ! ( '! (
! ! ! ! ! ! ' &! '
' (' (
& '& ' 2  MulBplicaBons/Divisions  
&0 0 0 0 aN−1 bN−1 cN−1 ' &d N−1 ' '0 0 0 0 aN−1 bN−1 cN−1 ( 'd N−1 (
&"0 '"0 0 0 0 0 aN bN ($ '"d N ($
0 0 0 0 aN bN '$ &"d N '$

!1 c1% 0 0 0 0 0 # !d1% #
2.  EliminaBon  Step   ' (' (
'0 b2% c2 0 0 0 0 ( 'd2% ( b2! = b2 − c1!a2
             R2=R2-­‐R1*a2   '0 a3 b3 c3 0 0 0 ( ' d3 (
d2! = d2 − d1!a2
!" A (1) #$, !"d (1) #$''0 0 a4 b4 c4 0 0
('
(, 'd4
(
(
'! ! ! ! ! ! ! ( '! ( 2  MulBplicaBons/Divisions  
' (' (
'0 0 0 0 aN−1 bN−1 cN−1 ( 'd N−1 ( 2  AddiBons/SubtracBons  
'"0 0 0 0 0 aN bN ($ '"d N ($

Number  of  FLOPs  involved  in  step  k=1  is  2+2+2=6  

No  middle  loop  
 Thomas  Algorithm  
1.  Normalize  pivot  row    
Step  k=2  
             with  pivot  element  
!1 c1% 0 0 0 0 0 # !d1% # c2 d! !1 c1% 0 0 0 0 0 # !d1% #
' (' ( c!2 = , d2!! = 2 ' (' (
'0 b2% c2 0 0 0 0 ( 'd2% ( b2! b2! '0 1 c%2 0 0 0 0 ( 'd2%% (
'0 a3 b3 c3 0 0 0 ( ' d3 ( '0 a3 b3 c3 0 0 0 ( ' d3 (
!" A (1) #$, !"d (1) #$''0 0 a4 b4 c4 0 0
('
(, 'd4
(
!" A *(1) #$, !"d *(1) #$''0 0 a4 b4 c4 0 0
(' (
( (, 'd4 (
'! ! ! ! ! ! ! ( '! ( '! ! ! ! ! ! ! ( '! (
' (' ( 2  MulBplicaBons/Divisions   ' (' (
'0 0 0 0 aN−1 bN−1 cN−1 ( 'd N−1 ( '0 0 0 0 aN−1 bN−1 cN−1 ( 'd N−1 (
'"0 aN bN ($ '"d N ($ '"0
0 0 0 0 0 0 0 0 aN bN ($ '"d N ($

2.  EliminaBon  Step   !1 c1% 0 0 0 0 0 # !d1% #

b3! = b3 − c!2 a3
             R3=R3-­‐R2*a3   '
'0 1 c%2 0 0 0 0
('
( 'd2%%
(
(
'0 ( 'd3% (
d3! = d3 − d2!!a3
0 b3% c3 0 0 0
!" A *(2) #$, !"d *(2) #$''0 0 a4 b4 c4 0 0
('
(, 'd4
(
( 2  MulBplicaBons/Divisions  
'! ! ! ! ! ! ! ( '! (
' (' ( 2  AddiBons/SubtracBons  
'0 0 0 0 aN−1 bN−1 cN−1 ( 'd N−1 (
'"0 0 0 0 0 aN bN ($ '"d N ($
Number  of  FLOPs  involved  in  step  k=2  is  2+2+2=6  
No  middle  loop  
 Thomas  Algorithm  
APer  Step  k=N-­‐1  
"1 c1& 0 0 0 0 0 $ "d1& $
' (' (
'0 1 c&2 0 0 0 0 ( 'd2&& ( Number  of  FLOPs  involved  in  EliminaBon=6×(N-­‐1)  
'0 0 1 c3& 0 0 0 ( 'd3&& (
"# A ( N−1) $%, "#d ( N−1) $%''0 0 0 1 c&4 0 0
('
(, 'd4&&
(
(
'! ! ! ! ! ! ! ( '! (
' (' (
'0 0 0 0 0 1 c&N−1 ( 'd &&N−1 (
'#0 0 0 0 0 0 b& (% '#d & (%
N N

d !N Number  of  FLOPs  involved  :  1  Div  

Backward  Sweep   xN =
b!N
xi = di!!− ci! × xi+1 Number  of  FLOPs  involved    :  (  1  Mult+  1  SubtracBon)×(N-­‐1)  
i=N-­‐1,  N-­‐2,  …1  

Total  Number  of  FLOPs  involved  in  Thomas  Algorithm~6×(N-­‐1)+1+2×(N-­‐1)~8N<<<N3  

 Start  with  iniBal  guess  

x i+1 y i+1

iniBal  guess  

0  

1  

2  

3  
Start  with  iniBal  guess  

i   xi+1=x1   yi+1=x2  
=x1   =x2

iniBal  guess   0   0  

0   3.5   1.25  

1   2.875   0.9375  

x (j ) − x (j
N N−1)
2   3.03125   1.015625  
(N )
< εtolerance for ∀j
xj
3   2.992   0.996  

4   3.001   1.009  
 n
a11 ≥ ∑ a1i
i=2
i≠1
n
a22 ≥ ∑ a2i
i=1
i≠2
n
aii ≥ ∑ aij
i=1
i≠ j
a22
 )  

2
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING  

CLL-­‐113  
 
Solution  of  non-­‐Linear  Equations  
Bracketing  Methods  
Prof.  Jaya*  Sarkar  
van  der  Waals  
 T0  ,  

CA   2.00  

2.00  
 =0  
itera2on  
 fl  &  fr  
y− fl f r − fl (xr,fr)  
l
= r
x−x x − xl
xl   xr  
0− fl fr− fl (x3,f3)  
i+1 l
= r l
x −x x −x (x2,f2)  

x i+1 l
=x −f l ( x r
− x l
) (xl,fl)  

( f r
− f l
)
 l
x &x r
fl× fr <0

x i+1 l
=x −f l (x r
−x ) l

(f r
−f ) l

x l = x i+1
x r = x i+1
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING  

CLL-­‐113  
 
Solution  of  non-­‐Linear  Equations  
Open  Methods  
Prof.  Jaya*  Sarkar  
Secant  Method   f(x)=x-­‐e^x/3   x^0=1   x^1=2  

tol=1*10-­‐4  

x^(i+1)=x^i-­‐f^i*(x^i-­‐
Itn   x^i-­‐1   x^i   f^i-­‐1   f^i   x^i-­‐1)/(f^i-­‐f^i-­‐1)   f^(i+1)   |x(i-­‐1)-­‐x(i+1)|  

1   1   2   0.093906057   -­‐0.4630187   1.16861534   0.096103886   0.16861534  

2   2   1.16861534   -­‐0.4630187   0.096103886   1.311516555   0.074250358   0.688483445  
3   1.16861534   1.311516555   0.096103886   0.074250358   1.79704301   -­‐0.213552036   0.62842767  
4   1.311516555   1.79704301   0.074250358   -­‐0.213552036   1.436777893   0.03440517   0.125261338  
5   1.79704301   1.436777893   -­‐0.213552036   0.03440517   1.486766287   0.012509487   0.310276723  
6   1.436777893   1.486766287   0.03440517   0.012509487   1.515325761   -­‐0.001642036   0.078547868  
7   1.486766287   1.515325761   0.012509487   -­‐0.001642036   1.512011934   6.27852E-­‐05   0.025245647  
8   1.515325761   1.512011934   -­‐0.001642036   6.27852E-­‐05   1.512133976   2.94813E-­‐07   0.003191785  
9   1.512011934   1.512133976   6.27852E-­‐05   2.94813E-­‐07   1.512134552   -­‐5.33749E-­‐11   0.000122617  
10   1.512133976   1.512134552   2.94813E-­‐07   -­‐5.33749E-­‐11   1.512134552   0   5.75656E-­‐07  
11   1.512134552   1.512134552   -­‐5.33749E-­‐11   0   1.512134552   0   1.0422E-­‐10  
Fixed Point Iteration

guess  

x i+1 = g ( x i )
x i+1 − x i < ε tol
Fixed  Point  Itera*on  Method   f(x)=x-­‐e^x/3   x^0=1  
tol=1*10-­‐4  
g(x)=Exp(x)/3   g(x)=LN(3x)  

itera*on  Number   x^i   g(x^i)   Error=|x^i1-­‐x^i-­‐1|   itera*on  Number   x^i   g(x^i)   Error=|x^i1-­‐x^i-­‐1|  
0   1   0.906093943   0   1   1.098612289  
1   0.906093943   0.824879185   0.093906057   1   1.098612289   1.192660116   0.098612289  
2   0.824879185   0.760535032   0.081214758   2   1.192660116   1.274798493   0.094047828  
3   0.760535032   0.713140191   0.064344153   3   1.274798493   1.34140041   0.082138377  
4   0.713140191   0.680129473   0.047394841   4   1.34140041   1.392326439   0.066601917  
5   0.680129473   0.658044437   0.033010718   5   1.392326439   1.429588334   0.050926029  
6   0.658044437   0.643670808   0.022085035   6   1.429588334   1.455998813   0.037261895  
7   0.643670808   0.634485097   0.014373629   7   1.455998813   1.474304423   0.026410479  
8   0.634485097   0.628683586   0.009185712   8   1.474304423   1.48679859   0.01830561  
9   0.628683586   0.625046831   0.005801511   9   1.48679859   1.4952375   0.012494167  
10   0.625046831   0.622777817   0.003636755   10   1.4952375   1.500897346   0.00843891  
11   0.622777817   0.621366327   0.002269014   11   1.500897346   1.504675448   0.005659846  
12   0.621366327   0.620489894   0.00141149   12   1.504675448   1.507189515   0.003778103  
13   0.620489894   0.619946314   0.000876433   13   1.507189515   1.508858957   0.002514066  
14   0.619946314   0.619609416   0.00054358   14   1.508858957   1.509965996   0.001669442  
15   0.619609416   0.619400705   0.000336899   15   1.509965996   1.51069942   0.001107039  
16   0.619400705   0.619271443   0.00020871   16   1.51069942   1.511185024   0.000733424  
17   0.619271443   0.6191914   0.000129262   17   1.511185024   1.511506416   0.000485604  
18   0.6191914   0.61914184   8.0043E-­‐05   18   1.511506416   1.511719069   0.000321392  
19   0.61914184   0.619111156   4.956E-­‐05   19   1.511719069   1.511859748   0.000212653  
20   0.619111156   0.61909216   3.06839E-­‐05   20   1.511859748   1.511952803   0.000140679  
21   0.61909216   0.619080399   1.89965E-­‐05   21   1.511952803   1.512014351   9.30548E-­‐05  
x i+1 − x i < ε tol
Newton  Raphson  Itera*on  Method   f(x)=x-­‐e^x/3   f`(x)=1-­‐e^x/3   tol  1*10^-­‐4  

itera*on  Number   x^i   f(x^i)   f`(x^i)   x^i+1   Error=|x^i1-­‐x^i-­‐1|  

0   2   -­‐0.4630187   -­‐1.4630187   1.683518263   0.316481737  
1   1.683518263   -­‐0.111303955   -­‐0.794822218   1.543481972   0.140036291  
2   1.543481972   -­‐0.016804879   -­‐0.560286851   1.513488623   0.029993349  
3   1.513488623   -­‐0.000694853   -­‐0.514183476   1.51213725   0.001351373  
4   1.51213725   -­‐1.38198E-­‐06   -­‐0.512138632   1.512134552   2.69846E-­‐06  
5   1.512134552   -­‐5.5056E-­‐12   -­‐0.512134552   1.512134552   1.07503E-­‐11  
6   1.512134552   0   -­‐0.512134552   1.512134552   0  
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING  

CLL-­‐113  
 
Solution  of  non-­‐Linear  Equations  
Error  Estimation/Convergence  
Prof.  Jaya*  Sarkar  
 /guesses

i l
3. Stopping criteria x − x < εtol or, x − x < εtol
3.  Si+1
t   i+1
Convergence

i i i−1
E = x −x
Bisection Method

itera*on:0  

itera*on:1  

itera*on:2  
Bisection Method

itera*on:0  

itera*on:1  

itera*on:2  
Rate of Convergence

Bisec*on  Method:     i+1 Ei 1

E = α = ,η = 1 Linear  Rate  of  Convergence  
2 2
Quadra3c  Rate  of  Convergence  
Newton  Raphson  Method   η=2
Non-­‐Linear  Rate  of  Convergence  
Regula  Falsi/Secant  Method   1<η < 2
Bisec*on  Method:    

E i < ε tol
L
Ei = i
< ε tol
2
! L $
ln # &
ε
" tol %
i=
ln 2
Bisec*on  Method   f(x)=x-­‐e^x/3   x^l=1   x^r=2   tol=1*10-­‐4  

x^(i+1)=0.5(x^l
Itn   x^l   x^r   f^l   f^r   +x^r)   f^(i+1)   f^l  f^(i+1)   Error  
1   1   2   0.093906057   -­‐0.4630187   1.5   0.006103643   0.000573169   0.5  
2   1.5   2   0.006103643   -­‐0.4630187   1.75   -­‐0.168200892   -­‐0.001026638   0.25  
3   1.5   1.75   0.006103643   -­‐0.168200892   1.625   -­‐0.067806346   -­‐0.000413866   0.125  
4   1.5   1.625   0.006103643   -­‐0.067806346   1.5625   -­‐0.027744394   -­‐0.000169342   0.0625  
5   1.5   1.5625   0.006103643   -­‐0.027744394   1.53125   -­‐0.010067718   -­‐6.14498E-­‐05   0.03125  
6   1.5   1.53125   0.006103643   -­‐0.010067718   1.515625   -­‐0.001796801   -­‐1.0967E-­‐05   0.015625  
7   1.5   1.515625   0.006103643   -­‐0.001796801   1.5078125   0.002199369   1.34242E-­‐05   0.0078125  
8   1.5078125   1.515625   0.002199369   -­‐0.001796801   1.51171875   0.000212816   4.6806E-­‐07   0.00390625  
9   1.51171875   1.515625   0.000212816   -­‐0.001796801   1.513671875   -­‐0.000789104   -­‐1.67934E-­‐07   0.001953125  
10   1.51171875   1.513671875   0.000212816   -­‐0.000789104   1.512695313   -­‐0.000287423   -­‐6.11681E-­‐08   0.000976563  
11   1.51171875   1.512695313   0.000212816   -­‐0.000287423   1.512207031   -­‐3.71233E-­‐05   -­‐7.90042E-­‐09   0.000488281  
12   1.51171875   1.512207031   0.000212816   -­‐3.71233E-­‐05   1.511962891   8.78913E-­‐05   1.87046E-­‐08   0.000244141  
13   1.511962891   1.512207031   8.78913E-­‐05   -­‐3.71233E-­‐05   1.512084961   2.53953E-­‐05   2.23202E-­‐09   0.00012207  
14   1.512084961   1.512207031   2.53953E-­‐05   -­‐3.71233E-­‐05   1.512145996   -­‐5.86119E-­‐06   -­‐1.48846E-­‐10   6.10352E-­‐05  

1   itera*on  (i)   Number  of  itera*ons  required  to  reach  desired  accuracy  
0   2   4   6   8   10   12   14   16  
0.1  

y  =  e-­‐0.693x   # 1 &
Error  

0.01   ln % (
$ 1×10 −4 '
0.001   i= = 13.2877 ≈ 14
ln(2)
0.0001  

0.00001  
 g ( x ) − g ( xi )
g! (ξ ) =
( x − x i
)
( )
x i

(x )
i
xi
(i+1)
= g! (ξ ) E ( ) = α E ( ) = α 2 E i−1 = ..... = α i E 1
i i
E
 0.2  

0.1  

0  
0   0.5   1   1.5   2   2.5  
-­‐0.1  
f(x)  -­‐0.2   x  
-­‐0.3  

-­‐0.4  

-­‐0.5  

itera*on  
Number   x^i   g(x^i)   Error=|x^i1-­‐x^i-­‐1|  
itera*on  Number   x^i   g(x^i)   Error=|x^i1-­‐x^i-­‐1|   0   1.52   1.524075065  
0   1.0986   0.999987711  
1   0.999987711   0.906082808   0.098612289   1   1.524075065   1.530298442   0.004075065  
2   0.906082808   0.82487   0.093904903   2   1.530298442   1.539851762   0.006223377  
3   0.82487   0.760528047   0.081212808  
4   0.760528047   0.71313521   0.064341953  
3   1.539851762   1.55463295   0.00955332  
5   0.71313521   0.680126085   0.047392837   4   1.55463295   1.577782944   0.014781189  
6   0.680126085   0.658042208   0.033009125  
7   0.658042208   0.643669373   0.022083877  
5   1.577782944   1.614734675   0.023149994  
8   0.643669373   0.634484186   0.014372835   6   1.614734675   1.675518026   0.036951731  
9   0.634484186   0.628683013   0.009185187  
7   1.675518026   1.7805205   0.060783351  
10   0.628683013   0.625046473   0.005801173  
11   0.625046473   0.622777594   0.00363654   8   1.7805205   1.977647904   0.105002474  
12   0.622777594   0.621366189   0.002268879   9   1.977647904   2.408575792   0.197127404  
13   0.621366189   0.620489808   0.001411405  
14   0.620489808   0.619946261   0.000876381   10   2.408575792   3.706038451   0.430927888  
15   0.619946261   0.619609383   0.000543547   11   3.706038451   13.5640941   1.297462659  
16   0.619609383   0.619400685   0.000336878  
17   0.619400685   0.61927143   0.000208698   12   13.5640941   259232.8096   9.858055654  
18   0.61927143   0.619191392   0.000129254   13   259232.8096   #NUM!   259219.2455  
19   0.619191392   0.619141835   8.00382E-­‐05  
 0.2  

0.1  

0  
0   0.5   1   1.5   2   2.5  
-­‐0.1  
f(x)  -­‐0.2   x  
-­‐0.3  

-­‐0.4  
itera*on  Number   x^i   g(x^i)   Error=|x^i1-­‐x^i-­‐1|  
-­‐0.5  
0   1   1.098612289  
1   1.098612289   1.192660116   0.098612289   itera*on  Number   x^i   g(x^i)   Error=|x^i1-­‐x^i-­‐1|  
2   1.192660116   1.274798493   0.094047828   0   0.619   0.618962282  
3   1.274798493   1.34140041   0.082138377   1   0.618962282   0.618901347   3.77176E-­‐05  
4   1.34140041   1.392326439   0.066601917   2   0.618901347   0.618802895   6.0935E-­‐05  
5   1.392326439   1.429588334   0.050926029   3   0.618802895   0.618643807   9.84519E-­‐05  
6   1.429588334   1.455998813   0.037261895   4   0.618643807   0.618386685   0.000159088  
7   1.455998813   1.474304423   0.026410479   5   0.618386685   0.617970975   0.000257123  
8   1.474304423   1.48679859   0.01830561   6   0.617970975   0.617298499   0.00041571  
9   1.48679859   1.4952375   0.012494167   7   0.617298499   0.616209708   0.000672475  
10   1.4952375   1.500897346   0.00843891  
8   0.616209708   0.614444351   0.001088791  
11   1.500897346   1.504675448   0.005659846  
9   0.614444351   0.611575375   0.001765357  
12   1.504675448   1.507189515   0.003778103  
10   0.611575375   0.606895219   0.002868976  
13   1.507189515   1.508858957   0.002514066  
14   1.508858957   1.509965996   0.001669442   11   0.606895219   0.599213165   0.004680156  
15   1.509965996   1.51069942   0.001107039   12   0.599213165   0.586474412   0.007682054  
16   1.51069942   1.511185024   0.000733424   13   0.586474412   0.564986049   0.012738753  
17   1.511185024   1.511506416   0.000485604   14   0.564986049   0.527658048   0.021488363  
18   1.511506416   1.511719069   0.000321392   15   0.527658048   0.459305447   0.037328001  
19   1.511719069   1.511859748   0.000212653   16   0.459305447   0.32057246   0.068352601  
20   1.511859748   1.511952803   0.000140679   17   0.32057246   -­‐0.039034656   0.138732987  
21   1.511952803   1.512014351   9.30548E-­‐05   18   -­‐0.039034656   #NUM!   0.359607116  
Fixed Point Iteration: Rate of Convergence
1  
0   2   4   6   8   10   12   14   16   18   20  

0.1  
itera*on  number  (i)  
g(x)=Exp(x)/3  

0.01  

Error(i+1)  
itera*on  Number   x^i   g(x^i)   Error=|x^i1-­‐x^i-­‐1|  
0   0.7   0.671250902  
0.001  
1   0.671250902   0.652227804   0.028749098  
2   0.652227804   0.639937678   0.019023099  
0.0001  
3   0.639937678   0.632120897   0.012290125  
4   0.632120897   0.627199008   0.007816781  
0.00001  
5   0.627199008   0.624119588   0.004921889  
6   0.624119588   0.622200619   0.003079419  
0.000001  
y  =  0.0508e-­‐0.473x  
7   0.622200619   0.621007779   0.00191897  
 
8   0.621007779   0.620267458   0.001192839  
9   0.620267458   0.619808431   0.000740321  
ln[E (
i+1)
10   0.619808431   0.619523988   0.000459027   ] = ln[E 1 ] + i × ln α
11   0.619523988   0.619347793   0.000284444  
12   0.619347793   0.619238677   0.000176195   α = g! (ξ )
13   0.619238677   0.619171112   0.000109116  
14   0.619171112   0.619129279   6.75652E-­‐05  
≈ g! x ( )
≈ 0.6192
ln ( 0.6192 ) ≈ -0.4793
Derivation of Newton Raphson Method
i
Current guess: x

Taylor Series Expansion around

xi
2

f ( x) = f ( x
( x−x ) i
i
) + (x − x i
) f "( x i
) +
2!
f "" ( x i ) +...

Assumptions: f ( x i+1 ) = 0, x i+1 ≈ x i

f (x )i

0 = f (x i
) + (x i+1
−x i
) f "( x ) i
x i+1
=x − i

f "( x )i
Error Estimation: Newton Raphson Method
Taylor Series of true value x 2

0 = f (x
( x−x ) i
i
) + ( x − x ) f "( x )
i i
+
2!
f "" ( x i ) +...
2
E (x i
) Mean Value Theorem
0= f ( x ) + E ( x ) f !( x ) +
i i i

2!
f !! ( x i
) + ...
2
E (x i
) xi < ξ < x
0 = f ( xi ) + E ( xi ) f !( xi ) + f !! (ξ ) (1)
2!
Newton Raphson Method
0 = f ( x i ) + ( x i+1 − x i ) f " ( x i )

0 = f ( x i ) + ( x i+1 − x + x − x i ) f " ( x i )

0 = f ( x i ) + (−E(x i+1 ) + E(x i )) f " ( x i ) …..(2)

Error Estimation: Newton Raphson Method
Taylor Series of true value x 2

0 = f (x
( x−x ) i
i
) + ( x − x ) f "( x )
i i
+
2!
f "" ( x i ) +...
2
E (x i
) Mean Value Theorem
0= f ( x ) + E ( x ) f !( x ) +
i i i

2!
f !! ( x i
) + ...
2
E (x i
) xi < ξ < x
0 = f ( xi ) + E ( xi ) f !( xi ) + f !! (ξ ) (1)
2!
Newton Raphson Method
0 = f ( x i ) + ( x i+1 − x i ) f " ( x i )

0 = f ( x i ) + ( x i+1 − x + x − x i ) f " ( x i )

0 = f ( x i ) + (−E(x i+1 ) + E(x i )) f " ( x i ) …..(2)

Error Estimation: Newton Raphson Method
Subtracting Eqn. 1 and Eqn. 2 we get:

2
E (x i
) f "" (ξ )
E (x i+1
)=− 2! f "( xi ) xi < ξ < x
i 2
E ( x ) f ## ( x i )
E ( x i+1 ) ≈ −
2! f #( xi )
$ f ## ( x i
) '
i 2
η  
E (x ) ≈ −
i+1
& ) E (x )
% 2! f # ( x ) )
i
& (

α  
Rate of convergence is QUADRATIC
 0.2  
0.1  
0  
itera*on   -­‐0.1   0   0.5   1   1.5   2   2.5  

Number   x^i   f(x^i)   f`(x^i)   x^i+1   Error=|x^i+1-­‐x^i|   f(x)  

-­‐0.2   x  
0   0.9   0.080132296   0.180132296   0.455147534   0.444852466   -­‐0.3  
-­‐0.4  
1   0.455147534   -­‐0.070321113   0.474531354   0.60333819   0.148190656  
-­‐0.5  
2   0.60333819   -­‐0.006065658   0.390596153   0.61886742   0.015529231  
3   0.61886742   -­‐7.38629E-­‐05   0.381058717   0.619061256   0.000193836  
4   0.619061256   -­‐1.16283E-­‐08   0.380938732   0.619061287   3.05254E-­‐08  

1  
0.1   0   1   2   3   4   5  
itera*on  no.  

Error(i+1)  
0.01  

f(x)=x-­‐e^x/ 0.001  
0.0001  
Newton  Raphson  Method   3   f`(x)=1-­‐e^x/3   tol  1*10^-­‐4   0.00001  
0.000001  
0.0000001  
1E-­‐08  
Itn.   Error=|x^i+1-­‐
No   x^i   f(x^i)   f`(x^i)   x^i+1   x^i|  
0   1.2   0.093294359   -­‐0.106705641   2.074315156   0.874315156  
1   2.074315156   -­‐0.578716129   -­‐1.653031285   1.724221281   0.350093876  
2   1.724221281   -­‐0.14516277   -­‐0.869384051   1.557249308   0.166971973  
3   1.557249308   -­‐0.024667085   -­‐0.581916393   1.51485991   0.042389398  
4   1.51485991   -­‐0.001401371   -­‐0.516261281   1.512145449   0.002714461  
5   1.512145449   -­‐5.58108E-­‐06   -­‐0.51215103   1.512134552   1.08973E-­‐05  
Convergence  Rate  of  different  Methods  
1  
0   2   4   6   8   10   12   14   16  
Itera*on  Number  
0.1  

0.01  
Error  (i+1)  

Bisec*on  
0.001  
FPI  
NR  
0.0001  
Secant  
Regula  Falsi  
0.00001  

0.000001  

0.0000001  
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING  

CLL-­‐113  
 
Solution  of  non-­‐Linear  Equations  
NR/Multivariable  case  
Prof.  Jaya*  Sarkar  
 Problem Areas: Newton Raphson Method
i+1 i
f (x ) i

x =x −
f "( x ) i 0.2  
0.1  
0  
-­‐0.1   0   0.5   1   1.5   2   2.5  
f(x)  
-­‐0.2   x  
-­‐0.3  
f(x)=x-­‐e^x/ -­‐0.4  
-­‐0.5  
Newton  Raphson  Itera*on  Method   3   f`(x)=1-­‐e^x/3   tol  1*10^-­‐4  

100  
itera*on   Error=|x^i+1-­‐
Number   x^i   f(x^i)   f`(x^i)   x^i+1   x^t|   itera*on  no.  

Errortrue(i+1)  
0   1.1   0.098611325   -­‐0.001388675   72.11110792   70.59897337  
1   72.11110792   -­‐6.92365E+30   -­‐6.92365E+30   71.11110792   69.59897337   10  

2   71.11110792   -­‐2.54707E+30   -­‐2.54707E+30   70.11110792   68.59897337  

3   70.11110792   -­‐9.37014E+29   -­‐9.37014E+29   69.11110792   67.59897337  
4   69.11110792   -­‐3.44708E+29   -­‐3.44708E+29   68.11110792   66.59897337  
1  
5   68.11110792   -­‐1.26811E+29   -­‐1.26811E+29   67.11110792   65.59897337   0   2   4   6  
Problem Areas: Newton Raphson Method
f ( x) = 0
if f ! ( x (i) ) = 0
f ( x (i ) )
x (i+1) = x (i ) −
f " ( x (i ) )
Newton Raphson Method: Poor Convergence

Loop  and  Divergence  

Oscilla;on  around    
Minima  or  Maxima    

An  ini;al  guess  near  one  root    

lands  near  a  loca;on  several    
roots  away  
itn  no   x(i)   f(i)   f`(i)   x(i+1)=   Error   f ( x ) = x10 −1 = 0
1   0.5   -­‐0.999023438   0.01953125   51.65   51.15  
2   51.65   1.35115E+17   2.61597E+16   46.485   5.165  
3   46.485   4.71117E+16   1.01348E+16   41.8365   4.6485   48  

f(x)  
4   41.8365   1.64268E+16   3.92643E+15   37.65285   4.18365  
5   37.65285   5.72768E+15   1.52118E+15   33.887565   3.765285  
6   33.887565   1.99712E+15   5.89336E+14   30.4988085   3.3887565   38  
7   30.4988085   6.96352E+14   2.28321E+14   27.44892765   3.04988085  
8   27.44892765   2.42803E+14   8.84562E+13   24.70403489   2.744892765  
9   24.70403489   8.46601E+13   3.42698E+13   22.2336314   2.470403488   28  
10   22.2336314   2.95192E+13   1.32768E+13   20.01026826   2.22336314  
11   20.01026826   1.02927E+13   5.14371E+12   18.00924143   2.001026826  
12   18.00924143   3.58884E+12   1.99278E+12   16.20831729   1.800924143   18  
13   16.20831729   1.25135E+12   7.72043E+11   14.58748556   1.620831729  
14   14.58748556   4.36319E+11   2.99105E+11   13.128737   1.458748556  
8  
15   13.128737   1.52135E+11   1.15879E+11   11.8158633   1.3128737  
16   11.8158633   53046236849   44894084748   10.63427697   1.18158633  
17   10.63427697   18496079117   17392888267   9.570849276   1.063427697  
-­‐2  
18   9.570849276   6449184014   6738361278   8.613764348   0.957084927   0   0.5   1   1.5   2  
19  
20  
8.613764348  
7.752387914  
2248691422  
784070216.9  
2610579222  
1011391879  
7.752387914  
6.977149123  
0.861376434  
0.77523879   x  
21   6.977149123   273388379.9   391833936.9   6.279434214   0.69771491  
22   6.279434214   95324633.58   151804496   5.651490799   0.627943415   100  
23   5.651490799   33237644.28   58812172.68   5.086341736   0.565149063  
24   5.086341736   11589249.7   22785041.39   4.577707606   0.50863413   10  
25   4.577707606   4040921.242   8827392.637   4.119936959   0.457770647  
26   4.119936959   1408981.851   3419913.618   3.707943555   0.411993403   1  
27   3.707943555   491281.3301   1324945.547   3.337149955   0.370793601   0   10   20   30   40   50  

Error  
28  
29  
3.337149955  
3.003436907  
171298.9439  
59727.98466  
513312.0965  
198868.7844  
3.003436907  
2.703098245  
0.333713047  
0.300338662  
0.1  
itn  
30   2.703098245   20825.59662   77047.13161   2.4328014   0.270296845  
31   2.4328014   7261.172653   29851.07068   2.189554759   0.24324664   0.01  
32   2.189554759   2531.55048   11566.50899   1.97068574   0.218869019  
33   1.97068574   882.4332477   4482.872281   1.773840237   0.196845503   0.001  
34   1.773840237   307.4217666   1738.723478   1.597031348   0.176808889  
35   1.597031348   106.9280696   675.8043276   1.438807931   0.158223417   0.0001  
36   1.438807931   37.02141119   264.2563358   1.298711343   0.140096589  
37   1.298711343   12.64980151   105.1026588   1.178354716   0.120356627   0.00001  
38   1.178354716   4.161315905   43.80103747   1.083349754   0.095004962  
39   1.083349754   1.226829103   20.55503401   1.023664661   0.059685092  
40   1.023664661   0.263505419   12.34296217   1.002316024   0.021348637  
41   1.002316024   0.023403117   10.21038368   1.000023934   0.00229209  
Improvements
f ( x (i ) )
x (i+1) = x (i ) −
f " ( x (i ) ) + β

f ( x (i) )
x (i+1) = x (i) − m
f " ( x (i) )
Improvements
f ( x (i ) )
x (i+1) = x (i ) −
f " ( x (i ) ) + β

f ( x (i) )
x (i+1) = x (i) − m
f " ( x (i) )
 f ( x) = 0 f ( x)
Improvements f !( x)
=0

Both  have  the  same  solu;on   x

f ( x)
1.  Modified  Equa;on  to  solve:   u ( x) = =0
f !( x)

Objec;ve  is  to  find  the  root  of  u(x)=0     u ( x (i) )

x (i+1) = x (i) −
u" ( x (i) )
f ( x (i) ) / f " ( x (i) )
= x (i) − 2
# f " ( x (i) ) f " ( x (i) ) − f ( x (i) ) f "" ( x (i) )% / # f " ( x (i) )%
$ & $ &
f ( x (i) ) f " ( x (i) )
= x (i) −
# f " ( x (i) ) f " ( x (i) ) − f ( x (i) ) f "" ( x (i) )%
$ &

2.  Modified  Equa;on  to  solve:   F ( x) = f 2 ( x) = 0

fn ( x1, x2 , x3 .......xn ) = 0
 E1(
100 )
E (j
i+1)
= x (j
i+1)
− x (j )
i
= 0.01

E2(
100 )
= 0.02

E3(
100 )
tol   = 0.002
 ∂g1 ∂g1 ∂g1
g1 ( x ) = g1 ( x ) +(i )

∂x1
(
(i )
|x ( i ) x1 − x1 +
∂x2
) (
(i )
|x ( i ) x2 − x2 +.... + )
∂xn
(
|x ( i ) xn − xn( )
i
)
∂g1 ∂g1 ∂g1
t
x =x
1
(i+1)
1 +
∂x1
(t (i )
|x ( i ) x1 − x1 + )
∂x2
(
t (i )
)
|x ( i ) x2 − x2 +.... +
∂xn
(
|x ( i ) xnt − xn( )
i
)
∂g1 ∂g1 ∂g1
E1(
i+1)
|x ( i ) E1( ) + |x ( i ) E2( ) +.... + |x ( i ) En( )
i i i
=
∂x1 ∂x2 ∂xn
 (i+1) ∂g j (i ) ∂g j (i ) ∂g j
|x ( i ) En( )
i
Ej = |x ( i ) E1 + |x ( i ) E2 +.... +
∂x1 ∂x2 ∂xn

" ∂g j ∂g j ∂g j %
$ + +... + ' |x ( i ) ≤ 1
# ∂x1 ∂x2 ∂xn &
for  j=1  to  n  

Very  Stringent  Condi;on  

Newton Raphson Method: Multivariable
f j ( x) = 0
(i)
x⇒x
∂f j ∂f j ∂f j
(
fj x
(i+1)
) ( )
= fj x +
(i)

∂x1
(
(i+1)
|x ( i ) x1 − x1 +(i )
)
∂x2
( (i+1)
|x ( i ) x2 − x2 +..... + )
(i )
∂xn
(
|x ( i ) xn( ) − xn( )
i+1 i
)
0   "Δx (i+1) %
$ 1 '
In  Matrix  form:   $Δx2(i+1) '
# x (i+1) − x (i) &
(
% 1 1) (
$
$ !'
'
% (i+1) i) ( $Δx (i+1) '
− fj x = %( )
(i)
# ∂f j ∂f j
.......
∂f j
& (
%
( |x ( i ) %
x 2 − x )
(
2 ( # n &
$ ∂x1 ∂x2 ∂xn ' (
% ! (
% x (i+1) − x (i) (
(
$ n )
n '

for  j=1  to  n  

 Jacobian  Matrix  (J)  

x
(i+1)
= x
(i)
( )
− J −1 f x
(i)

if  Singular  
" T −1 T %
x
(i+1)
= x
(i)
− $( J J ) J ' f x
# &
(i)
( )
Improvements  
" T −1 T %
x
(i+1)
= x
(i)
− $( J J + β I ) J ' f x
# &
(i)
( )
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING  

CLL-­‐113  
 
Solution  of  non-­‐Linear  Equations  
Root  Finding  of  Polynomials  
Prof.  Jaya*  Sarkar  
 Roots of Polynomial
15   30  

f(x)  
10  

f(x)  

f(x)  
25   5  
10  
20   0  
5   15   -­‐4   -­‐2   0   2   4  
-­‐5  
10   x  
-­‐10  
0   5  
-­‐4   -­‐2   0   2   4   -­‐15  
x   0   x  
-­‐5   -­‐20  
-­‐4   -­‐2   -­‐5   0   2   4  
-­‐10   -­‐25  

f (x) = x − 7x + 63 3
f (x) = x + x − 5x + 3 2
f (x) = x 3 − 3x 2 + x − 3
2
= ( x −1) ( x − 2 ) ( x + 3) = 0 = ( x −1) ( x + 3) = 0 = ( x 2 +1) ( x − 3) = 0

x = 1, 2, −3 Real Unique Roots x = 1,1, −3 Repeated Roots x = −i, i, 3 Complex Roots

 Roots of Polynomial: Newton Raphson Method
P3 (x) = x 3 − 7x + 6
= ( x −1) ( x − 2 ) ( x + 3) = ( x −1) P2 (x) itn   x(i)   f(x)   f`(x)  
x(i+1)   Error  
1   0   6   -­‐7  0.857142857  0.857142857
15   2   0.857142857   0.629737609   -­‐4.795918367  0.988449848  0.131306991

P3(x)  
3   0.988449848   0.046599285   -­‐4.068900694  0.999902397  0.011452549
10   4   0.999902397   0.00039044   -­‐4.000585589  0.999999993  9.75957E-­‐05

5   P2 (x) = x 2 + x − 6 = (x − 2)P1 (x)

itn x(i) f(x) f`(x) x(i+1) Error
20  

P2(x)  
0   1 1 "4 3 2.33333333 1.33333333
2 2.33333333 1.77777778 7 2.07936508 0.25396825
-­‐4   -­‐2   0   2   4   15   3 2.07936508 0.40312421 6.23809524 2.01474211 0.06462297
-­‐5   x   10   4 2.01474211 0.0739279 6.04422634 2.00251095 0.01223116
5 2.00251095 0.01256107 6.00753286 2.00042007 0.00209089
5   6 2.00042007 0.00210051 6.0012602 2.00007006 0.00035001
-­‐10  
0  
x   7 2.00007006 0.00035028 6.00021017 2.00001168 5.8378E"05

P1(x)  
10  
-­‐5   -­‐5   0   5  
-­‐10   5  

0   x  
-­‐5   0   5  
-­‐5  
x = 1, 2, −3 Real Unique Roots P1 (x) = x + 3
Roots of Polynomial: Newton Raphson Method

i+1 i
f ( xi )
x = x −m
f "( xi )

2
( x − (a + ib))( x − (a − ib)) = x − rx − s

Deflate the polynomial by 2

Pn (x) ⇒ Pn−2 (x) ⇒ Pn−4 (x) ⇒ ......... ⇒ P2 (x) if n is Even
Pn (x) ⇒ Pn−2 (x) ⇒ Pn−4 (x) ⇒ ......... ⇒ P1 (x) if n is Odd
Bairstow's Method
Pn ( x ) = a0 + a1 x + a2 x 2 + a3 x 3 +......... + an x n

Q ( x ) = x 2 − rx − s
Quadra*c  
Is  this  a  factor?  

Pn ( x ) = Pn−2 ( x ) Q(x) + R(x)

Pn−2 ( x ) = b2 + b3 x + b4 x 2 ........ + bn x n−2 R ( x ) = b0 + b1 ( x − r )

Bairstow's Method
Pn ( x ) = Pn−2 ( x ) Q(x) + R(x)

a0 + a1 x + a2 x 2 + a3 x 3 +......... + an x n
= ( b2 + b3 x + b4 x 2 ........ + bn x n−2 ) × ( x 2 − rx − s ) + #$b0 + b1 ( x − r )%&

a0 + a1 x + a2 x 2 + a3 x 3 +............. + an−2 x n−2 + an−1 x n−1 + an x n

= b2 x 2 + b3 x 3 + b4 x 4 ...... + bn−2 x n−2 + bn−1 x n−1 + bn x n
− b2 rx − b3rx 2 − b4 rx 3 − b5rx 4 ..... − bn−1rx n−2 − bn rx n−1 + 0
− b2 s − b3sx − b4 sx 2 − b5 sx 3 − b6 sx 4 ...... − bn sx n−2 + 0 +0
+ ( b0 − b1r ) + b1 x
Bairstow's Method
n
x : an = bn n
x : bn = an x n : bn = an
n−1 n−1 x n−1 : bn−1 = an−1 + bn r
x : an−1 = bn−1 − bn r x : bn−1 = an−1 + bn r
x i : bi = ai + bi+1r + bi+2 s
x n−2 : an−2 = bn−2 − bn−1r − bn s x n−2 : bn−2 = an−2 + bn−1r + bn s ......i = n - 2..to..0

! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !

x2 : a2 = b2 − b3r − b4 s x2 : b2 = a2 + b3r + b4 s
x1 : a1 = b1 − b2 r − b3s x1 : b1 = a1 + b2 r + b3s
x0 : a0 = b0 − b1r − b2 s x0 : b0 = a0 + b1r + b2 s
Bairstow's Method
x n : bn = an R ( x ) = b0 + b1 ( x − r )
x n−1 : bn−1 = an−1 + bn r
i
b0 ( r, s ) ⇒ 0
x : bi = ai + bi+1r + bi+2 s
......i = n - 2..to..0 b1 ( r, s ) ⇒ 0

0  
∂b1 ∂b1
b1 ( r + Δr, s + Δs ) = b1 ( r, s ) + Δr + Δs
∂r ∂s
0  
∂b0 ∂b0
b0 ( r + Δr, s + Δs ) = b0 ( r, s ) + Δr + Δs
∂r ∂s
" ∂b1 ∂b1 %
$ '"Δr % "−b1 ( r, s ) %
$ ∂r ∂s '
$ '=$ ' "Δr % "−b1 ( r, s ) % "Δr % "
−1 −b1 ( r, s )
%
$ ∂b0 ∂b0 '#Δs & $#−b0 ( r, s )'& [ M ]$ ' = $ ' $ ' = [M ] $ '
$# ∂r #Δs & $#−b0 ( r, s )'& #Δs & $#−b0 ( r, s )'&
∂s '&
Bairstow's Method
x n : cn = bn " ∂b1 ∂b1 %
x n−1 : cn−1 = bn−1 + cn r $ ' " %
∂s '"Δr % = $−b1 ( r, s ) '
!c2 c3 $!Δr $ !−b1 ( r, s ) $
$ ∂r $ ' # &# & = # &
x i : ci = bi + ci+1r + ci+2 s $ ∂b0 ∂b0 '#Δs & $#−b0 ( r, s )'& "c1 c2 %"Δs % #"−b0 ( r, s )&%
......i = n - 2..to..1 $# ∂r ∂s '&

c1
−b0 + b1
c2 −b1 − c3Δs r (i+1) = r (i) + Δr
Δs = , Δr =
c c2 s (i+1) = s (i) + Δs
c2 − 1 c3
c2

ε a,r =
Δr
100% Roots  
r
Δs r ± r 2 + 4s
ε a,s = 100% x=
s 2
Bairstow's Method
1. Pn-­‐2   is   a   third   order   polynomial   or   greater.   Bairstow's  
method   is   repeated   for   the   quoEent   for   finding   new  
values  of  r  and  s.  
 
2.   The   quoEent   is   quadraEc.   Roots   can   be   evaluated    
directly.   x=
2
r ± r + 4s
  2
3.  The  quoEent  is  first  order  polynomial,  the  remaining    
single  root  can  be  evaluated  as     x = −r / s
 P5 (x) = x 5 − 3.5x 4 + 2.75x 3 + 2.125x 2 − 3.875x +1.25

P5(x)   a5   a4   a3   a2   a1   a0  
1   -­‐3.5   2.75   2.125   -­‐3.875   1.25  

8  

P5(x)  
6  
4  
2   Real  Roots:  
0   x1=-­‐1,  
-­‐1.5   -­‐1   -­‐0.5  
-­‐2  
0   0.5   1   1.5   2   2.5  
x  
3  
x2=0.5  
-­‐4   x3=2  
-­‐6  
-­‐8  
 b5 = a5
P5 (x) = x 5 − 3.5x 4 + 2.75x 3 + 2.125x 2 − 3.875x +1.25 b4 = a4 + b5r
P5(x)   a5   a4   a3   a2   a1   a0  
1   -­‐3.5   2.75   2.125   -­‐3.875   1.25  
b3 = a3 + b4 r + b5 s
b2 = a2 + b3r + b4 s
itn   r   s   b5   b4   b3   b2   b1   b0  
1   -­‐1.00   -­‐1.00   1.00   -­‐4.50   6.25   0.38   -­‐10.50   11.38  
b1 = a1 + b2 r + b3s
2   -­‐0.64   0.14   1.00   -­‐4.14   5.56   -­‐2.03   -­‐1.80   2.13  
3   -­‐0.51   0.47   1.00   -­‐4.01   5.27   -­‐2.45   -­‐0.15   0.17  
b0 = a0 + b1r + b2 s
4   -­‐0.50   0.50   1.00   -­‐4.00   5.25   -­‐2.50   0.00   0.00  
5   -­‐0.50   0.50   1.00   -­‐4.00   5.25   -­‐2.50   0.00   0.00  
c5 = b5
6   -­‐0.50   0.50   1.00   -­‐4.00   5.25   -­‐2.50   0.00   0.00  
c4 = b4 + c5r
itn   c5   c4   c3   c2   c1   c3 = b3 + c4 r + c5 s
1   1.00   -­‐5.50   10.75   -­‐4.88   -­‐16.38  
2   1.00   -­‐4.79   8.78   -­‐8.34   4.79   c2 = b2 + c3r + c4 s
3   1.00   -­‐4.52   8.05   -­‐8.69   8.08  
4   1.00   -­‐4.50   8.00   -­‐8.75   8.37   c1 = b1 + c2 r + c3s
5   1.00   -­‐4.50   8.00   -­‐8.75   8.38  
6   1.00   -­‐4.50   8.00   -­‐8.75   8.38   c1
dels   delr   Errs=   Err_r  
−b0 + b1
c2 −b1 − c3Δs
1.138109017   0.35583014   824.0656852   55.23855773   Δs = , Δr =
0.331624608   0.133056764   70.5984393   26.03274392   c1 c2
0.030468695   0.011426623   6.091274153   2.286758475   c2 − c3
-­‐0.00020233   -­‐0.000313592   0.040466059   0.062718311   c2
1.03876E-­‐08   6.52645E-­‐08   2.07753E-­‐06   1.30529E-­‐05  
-­‐1.67709E-­‐14   -­‐1.08671E-­‐14   #DIV/0!   #DIV/0!  

(
x = −0.5 ± −0.5 + 4 × 0.5 2
) 2 x1 = 0.5, x2 = −1 (
x = r ± r 2 + 4s ) 2
 2 3 b5   b4   b3   b2  
P3 ( x ) = b2 + b3 x + b4 x + b5 x 1.00   -­‐4.00   5.25   -­‐2.50  

P3(x)   a3   a2   a1   a0  
1   -­‐4   5.25   -­‐2.5  

P3(x)  
2  
0  
-­‐2   -­‐1   0   1   2   3  
-­‐2  
x  
Real  Roots:  
-­‐4  
x1=-­‐1,  
-­‐6  
x2=0.5  
-­‐8  
-­‐10  
x3=2  
-­‐12  
-­‐14  
-­‐16  
-­‐18  
 P3 ( x ) = b2 + b3 x + b4 x 2 + b5 x 3 P3 (x) = x 3 − 4x 2 + 5.25x − 2.5
b3 = a3
P3(x)   a3   a2   a1   a0  
1   -­‐4   5.25   -­‐2.5  
b2 = a2 + b3r
itn   r   s   b3   b2   b1   b0   c3   c2   c1   dels   delr   Errs=   Err_r  
b1 = a1 + b2 r + b3s
1   0.00   5.00   1.00   -­‐4.00   10.25   -­‐22.50   1.00   -­‐4.00   15.25   88.42   24.67   3.584459459   0.264049955  

2   24.67   93.42   1.00   20.67   608.44   16936.41   1.00   45.33   1820.08   1445.09   -­‐45.30   -­‐70.04163723   -­‐0.029443207  
b0 = a0 + b1r + b2 s
3   -­‐20.63   1538.50   1.00   -­‐24.63   2051.95  -­‐80234.30   1.00   -­‐45.26   4524.33   -­‐2283.16   -­‐5.11   88.70166669   0.006859424  

4   -­‐25.74   -­‐744.65   1.00   -­‐29.74   26.09   21471.67   1.00   -­‐55.48   709.46   510.76   9.68   -­‐31.79737754   -­‐0.041372997  
c3 = b3
5   -­‐16.06   -­‐233.89   1.00   -­‐20.06   93.64   3185.91   1.00   -­‐36.13   440.05   180.68   7.59   -­‐21.33269308   -­‐0.142716826   c2 = b2 + c3r
6   -­‐8.47   -­‐53.21   1.00   -­‐12.47   57.66   172.60   1.00   -­‐20.94   181.80   54.93   5.38   -­‐17.75781197   3.126486167  
c1 = b1 + c2 r + c3s
7   -­‐3.09   1.72   1.00   -­‐7.09   28.91   -­‐104.11   1.00   -­‐10.19   62.13   17.68   4.57   11.94047118   0.235791492  

8   1.48   19.40   1.00   -­‐2.52   20.92   -­‐20.41   1.00   -­‐1.04   38.77   -­‐20.95   -­‐0.04   -­‐14.49929367   0.022664987  

9   1.45   -­‐1.56   1.00   -­‐2.55   0.00   1.48   1.00   -­‐1.11   -­‐3.16   0.37   0.34   0.209251221   -­‐0.284808763  

10   1.78   -­‐1.18   1.00   -­‐2.22   0.11   0.33   1.00   -­‐0.44   -­‐1.85   -­‐0.03   0.19   -­‐0.016737978   -­‐0.152346768  

11   1.97   -­‐1.22   1.00   -­‐2.03   0.03   0.04   1.00   -­‐0.06   -­‐1.31   -­‐0.03   0.03   -­‐0.01612642   -­‐0.025990154   c1
−b0 + b1
12   2.00   -­‐1.25   1.00   -­‐2.00   0.00   0.00   1.00   0.00   -­‐1.25   0.00   0.00   -­‐0.000526834   3.52232E-­‐05  
c2 −b1 − c3Δs
13   2.00   -­‐1.25   1.00   -­‐2.00   0.00   0.00   1.00   0.00   -­‐1.25   0.00   0.00   -­‐9.69272E-­‐10   -­‐6.18632E-­‐08   Δs = , Δr =
c1 c2
c2 − c3
c2
(
x = 2 ± 2 2 − 4 ×1.25 ) 2 x3 = (1+ 0.5i), x4 = (1− 0.5i)

P1 ( x ) = b2 + b3 x x5 = −b2 b3 = 2
(
x = r ± r 2 + 4s ) 2
 NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING  

CLL-­‐113  
 
Regression  Analysis  and  Curve  Fitting  
Least  Square  Method  
 
Prof.  Jaya*  Sarkar  
 Curve Fitting
ü  Data is given at discrete values along a continuum
•  Growth rate as a function of concentration:
C   0.5   0.8   1.5   2.5   4  
X   1.1   2.4   5.3   7.6   8.9  

𝑓(𝑥)   𝑓(𝑥)  
𝑓(𝑥)  

𝑥   𝑥   𝑥  
Did not attempt to connect Used straight line segments or
the points but characterizes linear interpolation to connect the Captures meandering
general upward trend . points. Fails when the underlying suggested by data
Curve fitting which does this curve is highly curve linear or data
is called ‘Regression’ is widely placed
 𝑓(𝑥)  
𝑓(𝑥)  

𝑥   𝑥  

Ø  When data is precise as obtained from steam tables.

Ø  Density ρ , viscosity η as a function of temperature.
Ø  Heat capacity Cp as a function of Temperature and Pressure
↓
The approach is to fit a curve that passes directly through each point.
Ø  Converts discrete data into analytical form such that any in between points within the range can
be interpolated.
Linear Regression
Data: (x1,y1) (x2,y2) (x3,y3) (x4,y4) ……(xN,yN)

Model: y=ao+a1 x

Predicted:  
Inadequate Criteria
Best fit line
p a s s e s
through two
points

Any straight line passing through the mid point also

minimizes the error so it gets cancelled.

Another criteria can be to minimize the sum of

absolute values of the discrepancies
1

2 Any of the straight line lying between straight lines 1 & 2

passing through the 4 points will satisfy this, we will not
end up with a UNIQUE solution.
Error Minimization Criteria

This gives undue advantage to a point

which is an outlier with a large error.
1
 No  of  
Data  
Point   x  i   y  i   x2  i   x  I  y  i   ypredicted  
1   1   0.5   1.00   0.50   0.91  
2   2   2.5   4.00   5.00   1.75  
3   3   2.0   9.00   6.00   2.59  
4   4   4.0   16.00   16.00   3.43  
5   5   3.5   25.00   17.50   4.27   7.0  
6   6   6.0   36.00   36.00   5.11   6.0  

7   7   5.5   49.00   38.50   5.95   5.0  

    28   24   140   119.5     4.0  

y  
3.0  
2.0   Data  

a1   a0   1.0   Predicted  
0.0  
0.8393   0.0714   0   2   4   6   8  

x  
Quantification of Error of Linear Regression
Standard Variation for the Regression Line

Standard Error of the Estimate

Standard Derivation for the data

2

Sy =
St
=
( y − y)
i

N −1 N −1

r: Correlation Coefficient
Qua*fica*on  of  Errors  

No  of  Data  
Point   x  i   y  i   ypredicted   (yi-­‐ybar)2   (yi-­‐ypredicted)2   Standard Error of the Estimate

1   1   0.5   0.91   8.58   0.17  

2   2   2.5   1.75   0.86   0.56  
3   3   2.0   2.59   2.04   0.35  
4   4   4.0   3.43   0.33   0.33  
Standard Derivation for the data
5   5   3.5   4.27   0.01   0.59  
6   6   6.0   5.11   6.61   0.80   2

7   7   5.5   5.95   4.29   0.20   Sy =

St
=
( yi − y)
    28   24   24   22.71   2.99   N −1 N −1
St   Sr   r2  
22.71   2.99   0.87  
           
Sy/x   Sy   r  
0.77   1.95   0.93   r2: Coefficient of Determination
Least  Square   x = b0 + b1 y
Method   n∑ xi yi − ∑ xi ∑ yi b0 = x − a1 y b1   b0  
b1 = 2
1.0346   0.4528  
No  of  
Data  
n∑ y −
6.5  
2
i (∑ y ) i

Point   x  i   y  i   y2  i   x  I  y  i   xpredicted  
1   1   0.5   0.25   0.50   0.97   5.5  
2   2   2.5   6.25   5.00   3.04   4.5  
3   3   2.0   4.00   6.00   2.52   3.5  

y  
4   4   4.0   16.00   16.00   4.59   2.5   Data  
5   5   3.5   12.25   17.50   4.07   1.5   y-­‐Predicted  
6   6   6.0   36.00   36.00   6.66   x-­‐Predicted  
0.5  
7   7   5.5   30.25   38.50   6.14  
    28   24   105   119.5     -­‐0.5  
0   2   4   6   8  
6.5   x  
6.5  
5.5  
5.5  
4.5  
4.5  
3.5  
3.5  
y  

y  
2.5  
Data   2.5  
1.5  
1.5   Data  
y-­‐Predicted  
0.5   0.5   x-­‐Predicted  
-­‐0.5   0   1   2   3   4   5   6   7   8   -­‐0.5  
0   2   4   6   8  
x   x  
 NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING  

CLL-­‐113  
 
Regression  Analysis  and  Curve  Fitting  
Multi-­‐linear  Regression  

Prof.  Jaya*  Sarkar  

 No  of  
Data  
Point   x  i   y  i   x2  i   y2  i   x  I  y  i   ypredicted  xpredicted  
1   1   0.5   1.00   0.25   0.50   0.91   0.97  
2   2   2.5   4.00   6.25   5.00   1.75   3.04  
3   3   2.0   9.00   4.00   6.00   2.59   2.52   x = b + b y
0 1
4   4   4.0   16.00   16.00   16.00   3.43   4.59   n∑ xi yi − ∑ xi ∑ yi
5   5   3.5   25.00   12.25   17.50   4.27   4.07   b1 = 2
n∑ yi − (∑ yi )
2

6   6   6.0   36.00   36.00   36.00   5.11   6.66  

7   7   5.5   49.00   30.25   38.50   5.95   6.14   b0 = x − a1 y
    28   24   140   105   119.5      
a1   a0   b1   b0   a0(P)=-­‐b0/b1   a1(P)=1/b1  
0.8393   0.0714   1.0346   0.4528   -­‐0.4377   0.9666  
 6.5  
5.5  
4.5  
3.5  

y  
2.5   Data  
1.5   y-­‐Predicted  
0.5   x-­‐Predicted  
-­‐0.5  
0   2   4   6   8  
x  

6.5  
6.5  
5.5  
5.5  
4.5  
4.5  
3.5  
3.5  
y  

y  
2.5  
Data   2.5  
1.5  
1.5   Data  
y-­‐Predicted  
0.5   0.5   x-­‐Predicted  
-­‐0.5   0   1   2   3   4   5   6   7   8   -­‐0.5  
0   2   4   6   8  
x   x  
 Polymer( Water( Compressive, Flexural,
cement,ratio, cement,ratio, Strength, Strength, 40.00-­‐60.00  
(x1) (x2) (MPa),(y1) (MPa)(y2)
20.00-­‐40.00  
Pure,cemet,mortar 0 0.4 50.63 7.39
50/50(5 5 0.4 45.31 8.59 60.00   0.00-­‐20.00  
50/50(10 10 0.4 35.63 9.30
50/50(15 15 0.4 22.21 7.30 Compress
50/50(20 20 0.4 27.18 6.40 ive   40.00  
Pure,cemet,mortar 0 0.5 42.40 7.11 Strength  
50/50(5 5 0.5 27.87 7.37 (MPa)   20.00  
50/50(10 10 0.5 26.30 6.90
50/50(15 15 0.5 29.60 6.70
50/50(20 20 0.5 24.27 5.10 0.00  0  
5  
Flexural  
Strength   8.00-­‐10.00   10  
10.00   0.5  
(MPa)   6.00-­‐8.00  
8.00  
4.00-­‐6.00   15  
6.00  
4.00  
20   0.4  
2.00  
0.00  0  

5  

Polymer  
10  
0.5  
Cement   15  
Water-­‐
Ra*o   20   0.4   Cement  Ra*o  
Solve  by  Gauss  EliminaHon/Gauss  Seidel  Method  

=  
 Polymer-­‐ Water-­‐ Compressiv "N
cement  raHo   cement   e  Strength   $
∑ x ∑u %'"a % "$∑ y
i i
0
i
%
'
$ '
  (x)   raHo  (u)   (MPa)  (y)   x2   u2   xu   yx   yu   $∑ x ∑ x x ∑u x ''$a ' = $$∑ y x '
$ i i i i i 1 i i
'
Pure  cemet  mortar   0   0.4   50.63   0   0.16   0   0  20.252   $ '
$∑ u
# i ∑ x u ∑u u '&#a & $#∑ y u
i i i i
2
i i
'
&
50/50-­‐5   5   0.4   45.31   25   0.16   2   226.55  18.124  
50/50-­‐10   10   0.4   35.63   100   0.16   4   356.3  14.252  
50/50-­‐15   15   0.4   22.21   225   0.16   6   333.15   8.884   !10 100 4.5$!a0 $ !331.4 $
20   400   0.16   8   543.6  10.872   # &# & # &
50/50-­‐20   0.4   27.18   # 100 1500 45 a
&# 1 & #= 2791.35 &
Pure  cemet  mortar   0   0.5   42.40   0   0.25   0   0   21.2   #"4.5 45 2.05 &%#"a2 &% #"147.6 &%
50/50-­‐5   5   0.5   27.87   25   0.25   2.5   139.35  13.935  
50/50-­‐10   10   0.5   26.30   100   0.25   5   263   13.15  
!a0 $ !71.1333 $
50/50-­‐15   15   0.5   29.60   225   0.25   7.5   444   14.8   # & # &
50/50-­‐20   20   0.5   24.27   400   0.25   10   485.4  12.135   #a1 & = #−1.0453 &
#"a2 &% #"−61.200 &%
N=10   100   4.5   331.4   1500   2.05   45.00   2791.35   147.60  

y = 71.1333−1.0453x − 61.200u
Actual   40.00-­‐60.00  
60.00   20.00-­‐40.00  
u\x   0   5   10   15   20   40.00  
y  
20.00   0.00-­‐20.00  
0.4  50.63   45.31   35.63   22.21   27.18   0.00  0  
0.5  42.40   27.87   26.30   29.60   24.27   5  
10  
0.5  
Predicted   x   15   u  
20   0.4  
u\x   0   5   10   15   20  
0.4   46.65   41.43   36.20   30.97   25.75   60.00  
0.5   40.53   35.31   30.08   24.85   19.63   40.00  
y  
20.00   40.00-­‐60.00  
0.00  0   20.00-­‐40.00  
5   0.00-­‐20.00  
y = 71.1333−1.0453x − 61.200u
x  10   0.5  
15   u  
20   0.4  
 No  of   !e1 $ !1 1 $ !0.5 $
Data   # & # & # &
Point   x  i   y  i  
#e2 & #1 2 & #2.5 &
#e3 & #1 3 & #2.0 &
1   1   0.5   # & # &!a0 $ # &
2   2   2.5   #e4 & + #1 4 &# & = # 4.0
a &
#e & #1 &" 1% #3.5 &
3   3   2.0   5 5
4   4   4.0  
# & # & # &
#e6 & #1 6 & #6.0 &
5   5   3.5   #"e &% #"1 7 &% #"5.5 &%
6   6   6.0   7

7   7   5.5   e + XΦ = Y
    28   24  
a1   a0  
0.8393   0.0714  
 !e1 $
# &
#e2 &
+ #e3 &
# &
#! &
#
"e N &
%
Y = XΦ + e
 NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING  

CLL-­‐113  
 
Regression  Analysis  and  Curve  Fitting  
Functional  and  non-­‐linear  Regression  

Prof.  Jaya*  Sarkar  

 !e1 $
# &
#e2 &
+ #e3 &
# &
#! &
#
"e N &
%

Y = XΦ + e

Φ
N  
Standard  Error  
used  
Functional Regression
 Exponential Model
Lineariza1on  

ln ( y) = ln (α1 ) + β1 x
Y = a0 + a1 x
Example:  Popula1on  growth,  Radioac1vity  Decrease,  Arrhenius  Law  
 Power Law Model

Lineariza1on  

ln ( y) = ln (α 2 ) + β2 ln ( x )
Y = a0 + a1 X
n−1
µ = mγ
Example:  non-­‐Newtonian  fluid  flow   Y = a0 + a1 X
 Saturation Growth Rate
Lineariza1on  

1 1 1
= + β3
y α3 x
Y = a0 + a1 X

Example:  Enzyma1c  kine1cs  

NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING  

CLL-­‐113  
 
Curve  Fitting  
Newtons  Divided  Difference  Formula  

Prof.  Jaya*  Sarkar  

 Curve Fitting and Interpolation
•   To  es&mate  intermediate  values  between  precise  data  points.    
( x1, y1 ) ( x2 , y2 )........ ( xn , yn )
• The  most  common  method  used  for  this  purpose  is  polynomial  
interpola&on  by  fi:ng  a  (n-­‐1)th  order  polynomial.  
2 n−1
Y = a0 + a1 x + a2 x +... + an−1 x
 Curve Fitting and Interpolation
2 n−1
( x1, y1 ) ( x2 , y2 )........ ( xn , yn ) Y = a0 + a1 x + a2 x +... + an−1 x
! y1 $ !1 x1 x12 .........x1n−1 $
# & # 2 n−1
&
# y2 & #1 x2 x2 .........x2 &
Y = #! & #
X= ! ! !......... ! &
# & # &
#! & #! ! !......... ! &
#" y &% # 2 n−1 &
n "1 xn xn .........xn %
n×n

−1
Φ = (X X) X Y
T T

T −1 Inverting X becomes difficult as n goes beyond 6

=X −1
(X ) T
X Y
= X −1Y
 Curve Fitting and Interpolation
S.no   Temp   P   Y   X                      
1   0   0.0002   0.0002   1   0   0   0   0   0  
2   20   0.0012   0.0012   1   20   400   8000   160000   3200000  
3   40   0.0060   0.0060   1   40   1600   64000   2560000   102400000  
4   60   0.0300  
5   80   0.0900   0.0300   1   60   3600   216000   12960000   777600000  
6   100   0.2700   0.0900   1   80   6400   512000   40960000   3276800000  
Phi   0.2700   1   100   10000   1000000   100000000   10000000000  
a0=0.0002   0.3000  
Data  Points  
a1=0.000852167   0.2500  
FiQed  Curve  

Pressure  
a2=-­‐8.14375E-­‐05   0.2000  
−1
a3=2.67604E-­‐06   Φ=X Y 0.1500  
a4=-­‐3.39063E-­‐08   0.1000  
a5=1.71354E-­‐10   0.0500  
0.0000  
Y = a0 + a1 x + a2 x 2 + a3 x 3 + a4 x 4 + a5 x 5 0   20   40   60  
Temp  
80   100   120  
 Newton`s Divided-Difference
Interpolating Polynomials
Linear Interpolation  The  simplest  form  of  interpola&on  is  to  connect  two  data  
points  with  a  straight  line.  Using  similar  triangles    

f1 ( x ) − f ( x0 ) f ( x1 ) − f ( x0 )
=
x − x0 x1 − x0

 finite-­‐divided-­‐difference  

f ( x1 ) − f ( x0 )
f1 ( x ) = f ( x0 ) + ( x − x0 )
( x1 − x0 )
 Newton`s Divided-Difference
Interpolating Polynomials
Linear Interpolation f ( x1 ) − f ( x0 )
f1 ( x ) = f ( x0 ) + ( x − x0 )
( x1 − x0 )
X   y=ln(x)   Series1   Series1:   Series2:  
X   f1(x)   f1(x)   2.00  
1   0.00   x0   x1  
1   0.00   0.00  
2   0.69   1   6.00  
2   0.36   0.46   1.50  
3   1.10   f(x0)   f(x1)   3   0.72   0.92  
0   1.79   Series-­‐1  

f(x)  
4   1.39   4   1.08   1.39   1.00  
Series2  
5   1.61   5   1.43   1.85   Data  
x0   x1  
6   1.79   1   6   1.79   2.31   0.50  
4.00   Series-­‐2  
f(x0)   f(x1)   Error(2)  
   
0.00  
0   1.39   Series1   48.30   0   5   10  
Series2   33.33   x  
 Newton`s Divided-Difference
Interpolating Polynomials
Quadratic    If  three  data  points  are  available,  this  can  be  accomplished  
Interpolation with  a  second-­‐order  polynomial  (also  called  a  quadra&c  
polynomial  or  a  parabola  ).  

f2 ( x ) = b0 + b1 ( x − x0 ) + b2 ( x − x0 ) ( x − x1 ) a0 = b0 − b1 x0 + b2 x0 x1
2
= b0 + b1 x − b1 x0 + b2 x + b2 x0 x1 − b2 xx0 − b2 xx1 a1 = b1 − b2 x0 − b2 x1
a2 = b2
f2 ( x ) = a0 + a1 x + a2 x 2
 Newton`s Divided-Difference
Interpolating Polynomials
Quadratic
Interpolation f2 ( x ) = b0 + b1 ( x − x0 ) + b2 ( x − x0 ) ( x − x1 )

b0 = f ( x0 )
Linear Interpolation
f ( x1 ) − f ( x0 )
b1 =
( x1 − x0 )
f ( x2 ) − f ( x1 ) f ( x1 ) − f ( x0 )
−
b2 =
( x2 − x1 ) ( x1 − x0 )
Second Order Curvature
( x2 − x0 )
 Newton`s Divided-Difference Interpolating Polynomials
Quadratic f2 ( x ) = b0 + b1 ( x − x0 ) + b2 ( x − x0 ) ( x − x1 )
Interpolation
x0   x1   x2  
1   4   6   b0   0.00  
b0 = f ( x0 )     Error(2)  
f(x0)   f(x1)   f(x2)   b1   0.46   Series1   18.32  
0   1.39   1.79   b2   -­‐0.05  
f ( x1 ) − f ( x0 )
b1 =
( x1 − x0 ) 2.00  
Series1:                
f2(x)  
f ( x1 ) − f ( x0 ) X   y=ln(x)   1.50  
f ( x2 ) − f ( x1 )
− 1   0.00   0.00  

f(x)  
( x2 − x1 ) ( x1 − x0 ) 1.00   Data  
b2 = 2   0.69   0.57  
( x2 − x0 ) 3   1.10   1.03   0.50  
Series-­‐1  

4   1.39   1.39  
0.00  
5   1.61   1.64   0   2   4   6   8  
6   1.79   1.79   x  
 Newton`s Divided-Difference
Interpolating Polynomials
General
Interpolating fn ( x ) = b0 + b1 ( x − x0 ) +... + bn ( x − x0 ) ( x − x1 ).... ( x − xn−1 )
Polynomial f (x ) − f x
 First  finite  divided  difference   f !" xi , x j #$ =
i ( ) j
b0 = f ( x0 )
(x − x )
i j
b1 = f [ x1, x0 ]  Second  finite  divided  difference  
f !" xi , x j #$ − f !" x j , xk #$
b2 = f [ x2 , x1, x0 ] f !" xi , x j , xk #$ =
( xi − x k )
!  nth  finite  divided  difference  
bn = f [ xn , xn−1,...., x1, x0 ] f [ xn , xn−1,...., x1 ] − f [ xn−1,...., x1, x0 ]
f [ xn , xn−1,...., x1, x0 ] =
( xn − x0 )
 Newton`s Divided-Difference Interpolating Polynomials
fn ( x ) = b0 + b1 ( x − x0 ) + b2 ( x − x0 ) ( x − x1 ) + b3 ( x − x0 ) ( x − x1 ) ( x − x2 )
Cubic
Interpolating i   xi   f(xi)   First   Second   Third  
Polynomial 0   x0   f(x0)   f[x1,x0]   f[x2,x1,x0]   f[x3,x2,x1,x0]  
1   x1   f(x1)   f[x2,x1]   f[x3,x2,x1]      
2   x2   f(x2)   f[x3,x2]          

3   x3   f(x3)              

i   xi   f(xi)   First   Second   Third   b0   0.00  

0   1   0.00   0.46   -­‐0.06   0.01   b1   0.46  
1   4   1.39   0.22   -­‐0.02      
b2   -­‐0.06  
2   5   1.61   0.18          
b3   0.01  
3   6   1.79              
X   y=ln(x)   Series1:f3(x)  
1   0.00   0.00  
2.0  
2   0.69   0.63  
3   1.10   1.08   1.5  
4   1.39   1.39  

f(x)  
5   1.61   1.61   1.0  
Series1  
6   1.79   1.79  
0.5  
Series2  
0.0  
    Error(2)   0   2   4   6   8  
Series1   9.29   x  
 NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING  

CLL-­‐113  
 
Curve  Fitting  
Lagrange  Polynomial  and  Newtons  Divided  Difference  Formula  

Prof.  Jaya*  Sarkar  

Newton`s Divided-Difference Error Analysis
fn ( x ) = b0 + b1 ( x − x0 ) +... + bn ( x − x0 ) ( x − x1 ).... ( x − xn−1 )
Also, as was the case with the Taylor series, a formulation for the truncation
error can be obtained.
(n+1)
f (ξ) n+1
Rn = ( xi+1 − xi )
(n +1)!
(n+1)
f (ξ)
= ( x − x0 ) ( x − x1 ).... ( x − xn )
(n +1)!
= f [ x, xn , xn−1,...., x1, x0 ] ( x − x0 ) ( x − x1 ).... ( x − xn )
≅ f [ xn+1, xn , xn−1,...., x1, x0 ] ( x − x0 ) ( x − x1 ).... ( x − xn )
 Newton`s Divided-Difference Error Analysis
fn ( x ) = b0 + b1 ( x − x0 ) +... + bn ( x − x0 ) ( x − x1 ).... ( x − xn−1 )

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

fn+1 ( x ) = b0 + b1 ( x − x0 ) +... + bn ( x − x0 ) ( x − x1 ).... ( x − xn−1 ) + bn+1 ( x − x0 ) ( x − x1 ).... ( x − xn )

fn+1 ( x ) = fn ( x ) + bn+1 ( x − x0 ) ( x − x1 ).... ( x − xn )

fn+1 ( x ) = fn ( x ) + Rn
 Ea ( 2 ) = 0.69 − 0.57 = 0.12
X   y=ln(x)   f2(x)  
1   0.00   0.00   R ( 2 ) = 0.63− 0.57 = 0.06
2   0.69   0.57  
3   1.10   1.03  
4   1.39   1.39  
5   1.61   1.64  
6   1.79   1.79  

X   y=ln(x)   f3(x)  
1   0.00   0.00  
2   0.69   0.63  
3   1.10   1.08  
4   1.39   1.39  
5   1.61   1.61  
6   1.79   1.79  
2  
n n n

n n
n-­‐1  
1  
 Lagrange Interpolation
(1, 0) ( 4,1.386294)

f1 ( x ) = y1P1 (x) + y2 P2 (x) P1 (x) =

( x − x2 )
, P2 (x) =
( x − x1 )
( x1 − x2 ) ( x2 − x1 )
P1 (2) =
( 2 − 4)
=, P2 (2) =
( 2 −1)
f1 ( 2 ) = y1P1 (2) + y2 P2 (2)
(1− 4) ( 4 −1)

f1 ( 2 ) = 0 × 0.667 +1.386295 × 0.333 = 0.462098

Lagrange Interpolation (1, 0) ( 4,1.386294) (6,1.791760)
f2 ( x ) = y1P1 (x) + y2 P2 (x) + y3 P3 (x)

P1 (x) =
( x − x 2 ) ( x − x3 )
, P2 (x) =
( x − x1 ) ( x − x3 )
, P3 (x) =
( x − x1 ) ( x − x2 )
( x1 − x2 ) ( x1 − x3 ) ( x2 − x1 ) ( x2 − x3 ) ( x3 − x1 ) ( x3 − x2 )
f2 ( 2 ) = y1P1 (2) + y2 P2 (2) + y3 P3 (2)

P1 (2) =
( 2 − 4) ( 2 − 6 )
, P2 (2) =
( 2 −1) ( 2 − 6 )
, P3 (2) =
( 2 −1) ( 2 − 4)
(1− 4) (1− 6) ( 4 −1) ( 4 − 6) (6 −1) (6 − 4)
f2 ( 2 ) = 0 × 0.533+1.386295 × 0.666 +1.79176 × (−0.2 ) = 0.565
 Derivation of Lagrange Polynomial from
Newton Divided Difference Formula
( x1, y1 ) ( x2 , y2 ) f ( x1 ) − f ( x2 ) f ( x1 ) f ( x2 )
f [ x1, x2 ] = = +
( x1 − x2 ) ( x1 − x2 ) ( x2 − x1 )
f1 ( x ) = f ( x1 ) + f [ x1, x2 ] ( x − x1 )
" f (x ) f ( x ) %
1 2
= f ( x1 ) + $$ + '' ( x − x1 )
# ( x1 − x2 ) ( x2 − x1 ) &

= f ( x1 ) ×
( x − x2 )
+ f ( x2 ) ×
( x − x1 )
( x1 − x2 ) ( x2 − x1 )
= y1 × P1 ( x ) + y2 × P2 ( x )
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING

CLL-­‐113  
Curve  Fitting  
Numerical  Differentiation

Prof.  Jayati  Sarkar

Differentiation

xi x i+1
Differentiation FD

xi x i+1

BD

x i-1 x i

CD

x i-1 x i x i+1
+
(Backward Difference)
=
undetermined  coefficients:
3  point  difference  formula
 Truncation Error
2 3

f ! (xi ) = f (xi )[a1 + a2 + a3 ] + Δxf ! (xi )[a1 − a3 ] +

(Δx) f !! (xi )[a1 + a3 ] +
(Δx) f !!! (xi )[a1 − a3 ]
2! 3!
1 1

Δx

Error 3

=
(Δx)
f """ (xi )[a1 − a3 ]
3!
2

≈
(Δx)
f ### (xi )
3!
+

Error
 Error
−1 2 3
a3 = a2 = a1 = −
2Δx Δx 2Δx
 Comparisons
fi+1 − fi Δx ##
Forward  Difference: f !(xi ) = − fi (ξ )
Δx 2!
2
Central  Difference: fi+1 − fi−1 (Δx) !!!
f !(xi ) = − fi (ξ )
2Δx 3!
3  Pt  FWD  Difference: 2
−3 fi + 4 fi+1 − fi+2 (Δx) !!!
f !(xi ) = − fi (ξ )
2Δx 3
4 3 2
f (x) = −0.1x − 0.15x − 0.5x − 0.25x +1.2
f !(x) = −0.4x3 − 0.45x2 − x − 0.25 = −0.9125@x = 0.5
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING

CLL-­‐113  
Numerical  Differentiation  

Prof.  Jayati  Sarkar

 Comparisons
fi+1 − fi Δx !!
Forward  Difference: f !(xi ) = − fi (ξ )
Δx 2!
2
Central  Difference: fi+1 − fi−1 (Δx) !!!
f !(xi ) = − fi (ξ )
2Δx 3!
3  Pt  FWD  Difference: 2
−3 fi + 4 fi+1 − fi+2 (Δx) !!!
f !(xi ) = − fi (ξ )
2Δx 3
4 3 2
f (x) = −0.1x − 0.15x − 0.5x − 0.25x +1.2
3 2
f !(x) = −0.4x − 0.45x − x − 0.25 = −0.9125@x = 0.5
delx x(i-1) x(i) x(i+1) x(i+2) f(i-1) f(i) f(i+1) f(i+2)
1 -­‐5.000000000000E-­‐01 5.00E-­‐01 1.500000000000E+00 2.500000000000E+00 1.212500000000E+00 9.25E-­‐01 -­‐1.312500000000E+00 -­‐8.800000000000E+00
0.1 4.000000000000E-­‐01 5.00E-­‐01 6.000000000000E-­‐01 7.000000000000E-­‐01 1.007840000000E+00 9.25E-­‐01 8.246400000000E-­‐01 7.045400000000E-­‐01
0.01 4.900000000000E-­‐01 5.00E-­‐01 5.100000000000E-­‐01 5.200000000000E-­‐01 9.340378490000E-­‐01 9.25E-­‐01 9.157871490000E-­‐01 9.063971840000E-­‐01
0.001 4.990000000000E-­‐01 5.00E-­‐01 5.010000000000E-­‐01 5.020000000000E-­‐01 9.259116253499E-­‐01 9.25E-­‐01 9.240866246499E-­‐01 9.231714971984E-­‐01
0.0001 4.999000000000E-­‐01 5.00E-­‐01 5.001000000000E-­‐01 5.002000000000E-­‐01 9.250912412504E-­‐01 9.25E-­‐01 9.249087412497E-­‐01 9.248174649972E-­‐01
0.0000 4.999900000000E-­‐01 5.00E-­‐01 5.000100000000E-­‐01 5.000200000000E-­‐01 9.250091249125E-­‐01 9.25E-­‐01 9.249908749125E-­‐01 9.249817496500E-­‐01
0.0000
01 4.999990000000E-­‐01 5.00E-­‐01 5.000010000000E-­‐01 5.000020000000E-­‐01 9.250009124991E-­‐01 9.25E-­‐01 9.249990874991E-­‐01 9.249981749965E-­‐01
0.0000
001 4.999999000000E-­‐01 5.00E-­‐01 5.000001000000E-­‐01 5.000002000000E-­‐01 9.250000912500E-­‐01 9.25E-­‐01 9.249999087500E-­‐01 9.249998175000E-­‐01
0.0000
0001 4.999999900000E-­‐01 5.00E-­‐01 5.000000100000E-­‐01 5.000000200000E-­‐01 9.250000091250E-­‐01 9.25E-­‐01 9.249999908750E-­‐01 9.249999817500E-­‐01
0.0000
00001 4.999999990000E-­‐01 5.00E-­‐01 5.000000010000E-­‐01
Forward 5.000000020000E-­‐01
Central 9.250000009125E-­‐01
3PtFWD 9.25E-­‐01 9.249999990875E-­‐01 9.249999981750E-­‐01
1E-­‐10 4.999999999000E-­‐01 5.00E-­‐01 5.000000001000E-­‐01 5.000000002000E-­‐01 9.250000000913E-­‐01 9.25E-­‐01 9.249999999088E-­‐01 9.249999998175E-­‐01
1E-­‐11 4.999999999900E-­‐01 Difference Difference
5.00E-­‐01 5.000000000100E-­‐01 5.000000000200E-­‐01 Difference
9.250000000091E-­‐01 9.25E-­‐01
  Error
9.249999999909E-­‐01  
9.249999999818E-­‐01
1E-­‐12 4.999999999990E-­‐01 f`(i) f`(i)
5.00E-­‐01 5.000000000010E-­‐01 5.000000000020E-­‐01 f`(i)
9.250000000009E-­‐01 9.25E-­‐01 E_FD
9.249999999991E-­‐01 E_CD
9.249999999982E-­‐01 E_3Pt
1E-­‐13 4.999999999999E-­‐01 5.00E-­‐01 5.000000000001E-­‐01 5.000000000002E-­‐01 9.250000000001E-­‐01 9.25E-­‐01 9.249999999999E-­‐01 9.249999999998E-­‐01
-­‐2.237500000000E+00 -­‐1.262500000000E+00 3.875000000000E-­‐01 1.325000000000 0.350000000000 1.300000000000
1E-­‐14 5.000000000000E-­‐01 5.00E-­‐01 5.000000000000E-­‐01 5.000000000000E-­‐01 9.250000000000E-­‐01 9.25E-­‐01 9.250000000000E-­‐01 9.250000000000E-­‐01
1E-­‐15 5.000000000000E-­‐01
-­‐1.003600000000E+00 -­‐9.160000000000E-­‐01
5.00E-­‐01 5.000000000000E-­‐01 5.000000000000E-­‐01
-­‐9.049000000000E-­‐01 0.091100000000
9.250000000000E-­‐01 9.25E-­‐01 9.250000000000E-­‐01
0.003500000000
9.250000000000E-­‐01
0.007600000000
-­‐9.212851000000E-­‐01 -­‐9.125350000000E-­‐01 -­‐9.124294000000E-­‐01 0.008785100000 0.000035000000 0.000070600000
-­‐9.133753500999E-­‐01 -­‐9.125003500000E-­‐01 -­‐9.124992994000E-­‐01 0.000875350100 0.000000350000 0.000000700600
-­‐9.125875034999E-­‐01 -­‐9.125000034998E-­‐01 -­‐9.124999929999E-­‐01 0.000087503500 0.000000003500 0.000000007000
f (x)= −0.1x4 −0.15x3 −0.5x2 −0.25x+1.2 -­‐9.125087500284E-­‐01 -­‐9.125000000332E-­‐01 -­‐9.124999999222E-­‐01 0.000008750028 0.000000000033 0.000000000078
-­‐9.125008749722E-­‐01 -­‐9.125000000054E-­‐01 -­‐9.125000000054E-­‐01 0.000000874972 0.000000000005 0.000000000005
-­‐9.125000866028E-­‐01 -­‐9.124999994503E-­‐01 -­‐9.124999988952E-­‐01 0.000000086603 0.000000000550 0.000000001105
-­‐9.125000088872E-­‐01 -­‐9.125000033361E-­‐01 -­‐9.125000088872E-­‐01 0.000000008887 0.000000003336 0.000000008887
f !(x) = −0.4x3 − 0.45x2 − x− 0.25 -­‐9.124999644783E-­‐01 -­‐9.125000199894E-­‐01 -­‐9.124999644783E-­‐01 0.000000035522 0.000000019989 0.000000035522
-­‐9.125000755006E-­‐01 -­‐9.125000755006E-­‐01 -­‐9.125006306121E-­‐01 0.000000075501 0.000000075501 0.000000630612
= −0.9125@x = 0.5 -­‐9.125034061697E-­‐01 -­‐9.125034061697E-­‐01 -­‐9.125145083999E-­‐01 0.000003406170 0.000003406170 0.000014508400
-­‐9.124923039394E-­‐01 -­‐9.124923039394E-­‐01 -­‐9.126033262419E-­‐01 0.000007696061 0.000007696061 0.000103326242
-­‐9.126033262419E-­‐01 -­‐9.126033262419E-­‐01 -­‐9.137135492665E-­‐01 0.000103326242 0.000103326242 0.001213549267
-­‐9.103828801926E-­‐01 -­‐9.103828801926E-­‐01 -­‐9.159339953158E-­‐01 0.002117119807 0.002117119807 0.003433995316
-­‐8.881784197001E-­‐01 -­‐9.436895709314E-­‐01 -­‐9.436895709314E-­‐01 0.024321580300 0.031189570931 0.031189570931
Error Comparisons
fi+1 − fi Δx !!
Forward  Difference: f !(xi ) =
Δx
−
2!
fi (ξ )
2
(Δx) f !!! ξ
Central  Difference:
f −f
f !(xi ) = i+1 i−1 − i ( )
2Δx 3!
1.00E+01
3  Pt  FWD  Diff: 2
1.00E+00

−3f + 4 f − f (Δx) 1.00E-­‐01

f !(xi ) = i i+1 i+2 − fi!!! (ξ ) 1.00E-­‐02
2Δx 3
1.00E-­‐03
1.00E-­‐04
Error

1.00E-­‐05
1.00E-­‐06 FWD
1.00E-­‐07
CD
3Pt
1.00E-­‐08
1.00E-­‐09
1.00E-­‐10

1.0E-­‐10 1.0E-­‐08 1.0E-­‐06 1.0E-­‐04 1.0E-­‐02 1.0E+00

 4

4
 1.00E+01
4
1.00E-­‐01
1.00E-­‐03

Error
4
1.00E-­‐05
1.00E-­‐07
1.00E-­‐09
1.0E-­‐10 1.0E-­‐07 1.0E-­‐04

delx
 1.00E-­‐01
1.00E-­‐02
1.00E-­‐03
1.00E-­‐04
1.00E-­‐05

Error
1.00E-­‐06
1.00E-­‐07
1.00E-­‐08
1.00E-­‐09
1.00E-­‐10

1.0E-­‐10 1.0E-­‐07 1.0E-­‐04 1.0E-­‐01

delx
8
Summary

minimum
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING

CLL-­‐113  
Numerical  Integration  
Newton  Cotes  Formula

Prof.  Jayati  Sarkar

Newton's Forward Difference Method
(Curve Fitting)

xN
Newton's Forward Difference Method (Curve
Fitting)

-­‐
 Newton's Forward Difference Method (Curve
Fitting)

2
4. α = 3
Δ 3 y1 = Δ (Δ 2 y1 )
3
y4 = b0 + 3b1 + 3.2b2 + 3.2.1b3 Δ y1
1
b3 = = Δ (Δ (Δy1 ))
b3 = (y4 − b0 − 3b1 − 6b2 ) 3! = Δ (Δ (y2 − y1 ))
6
1% 1 2 ( = Δ (Δy2 − Δy1 )
= ' y4 − y1 − 3 (Δy1 ) − 6 × (Δ y1 )*
6& 2! )
= Δ ((y3 − y2 ) − (y2 − y1 ))
1% 1 (
= ' y4 − y1 − 3 (y2 − y1 ) − 6 × (y3 − 2y2 + y1 )* = Δy3 − 2Δy2 + Δy1
3! & 2! )
= (y4 − y3 ) − 2 (y3 − y2 ) + (y2 − y1 )
1
= (y4 − 3y3 + 3y2 − y1 ) = y4 − 3y3 + 3y2 − y1
 Newton's Forward Difference Method (Curve
Fitting)

2
2 3 i
b1 = y1 Δ y1 Δ y1 Δ y1
b2 = b3 = bi =
2! 3! i!

Δ 2 y1 Δ N−1 y1
P (α ) = y1 + Δy1α + α (α −1) +..... + α (α −1).. (α − (N − 2 ))
2! (N −1)!
Integration-Applications
1 b
• To calculate the mean: = ∫ f (x)dx
b− a a

• Mass flow:

• Net heat loss:

Numerical Integration
Geometric Interpretation of Integration

• Newton Cotes Integration Formula:

Euler Backward: I = f (x1 )(x2 − x1 )

truncation error
x1 x2
h
Geometric Interpretation of Integration
• Newton Cotes Integration Formula:

Euler Forward: I = f (x2 )(x2 − x1 )

truncation error

x1 x2
h
Trapezoidal Rule:

1
I = ( f (x1 )+ f (x2 ))(x2 − x1 )
2
truncation error
x1 x2
h
Trapezoidal Rule:

x1 x2
h
Newton's Forward Difference Method
Δ 2 y1 Δ N−1 y1
y (α ) = y1 + Δy1α + α (α −1) +..... + α (α −1).. (α − (N − 2 ))
2! (N −1)!

2!

Δ 2 y1 = Δ (Δy1 )
= Δ (y2 − y1 )
= (Δy2 − Δy1 )
= ((y3 − y2 ) − (y2 − y1 ))
= y3 − 2y2 + y1
x1 x2
h
Trapezoidal  Rule Error
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING

CLL-­‐113  
Numerical  Integration  
Newton  Cotes  Formula

Prof.  Jayati  Sarkar

Trapezoidal Rule:

1
I = ( f (x1 )+ f (x2 ))(x2 − x1 )
2
truncation error
x1 x2
h
Simpsons 1/3rd and 3/8th Rule

a b a b
x1 x2 x3 x1 x2 x3 x4
h = (b − a) 2
h = (b − a) 3
2
Δ y1 Δ 2 y1 Δ 3 y1
y (α ) = y1 + Δy1α + α (α −1) y (α ) = y1 + Δy1α + α (α −1) + α (α −1)(α − 2 )
2! 2! 3!
h3 f """ (ξ ) h4 f iv (ξ )
E (α ) = α (α −1)(α − 2 ) E (α ) = α (α −1)(α − 2 )(α − 3)
3! 4!
 Simpsons 1/3rd Rule
h=
(b − a)
2

x : x1 → x3 x3 2

α :0 → 2
∫ ydx → ∫ yhdα
x1 0

a b
x1 h x2 h x3

I
 Simpsons 3/8th Rule
h=
(b − a)
3
x4 3
x : x1 → x4
∫ ydx → ∫ yhdα
α :0 →3 x1 0

a b
x1 x2 x3 x4

I      =

Simpsons 3/8th Rule

 Error Analysis of Simpsons 3/8th Rule
h=
(b − a)
3
x4 3
x : x1 → x4
∫ ydx → ∫ yhdα
α :0 →3 x1 0
3
h5 f iv (ξ ) 27h5 f iv (ξ )# 9 9 11 &
E (α ) = ∫ α (α −1)(α − 2 )(α − 3)dα = %$ − + −1('
0 4! 24 5 2 3
27h5 f iv (ξ ) a b
= [54 −135 +110 − 30 ]
=
h 5 f iv (ξ ) 3

∫ (α 4
− 6α 3 +11α 2 − 6α )dα
24 × 30 x1 x2 x3 x4
4! 0
3h5 f iv (ξ ) 3h5 f iv (ξ )
=− −
3 80
=
h5 f iv (ξ )$α 5 3α 4 11α 3
& − + − 3α 2
'
)
E    = 80
4! % 5 2 3 (0 5
(b − a) f iv (ξ )
−
=
h5 f iv (ξ )$ 35 35 11× 33
− + − 33
' 6480
& )
4! % 5 2 3 ( Simpsons 3/8th Rule
 Error in Simpsons 1/3 rd Rule
h=
(b − a)
2

x : x1 → x3 x3 2

α :0 → 2
∫ ydx → ∫ yhdα
x1 0
a b
2
h4 f """ (ξ ) 2
h5 f iv (ξ ) x1 x2 x3
E (α ) = ∫ α (α −1)(α − 2 )dα + ∫ α (α −1)(α − 2 )(α − 3)dα
0 3! 0 4!
h h
h4 f !!! (ξ ) 2
3 h5 f iv (ξ ) 2
= ∫ (α − 3α 2 + 2α )dα + ∫ (α 4
− 6α 3 +11α 2 − 6α )dα
3! 0 4! 0
h5 f iv (ξ )
4 2 5 2 E =−
90
4 iv 5
h f !!! (ξ )% α 3 2
( h f (ξ ) % α 3 4 11 3 2
(
= ' −α +α * + ' − α + α − 3α *
3! & 4 )0 4! & 5 2 3 )0
5iv

=
4
h f !!! (ξ )$ 2
− 2 3
+
4
2 2
' h f (ξ )$ 2
+ −
3
5
× 2
iv
4
+
11
× 2 3
− 3×
5
2 2
'
=−
(b − a) f (ξ)
& ) & )
3! % 4 ( 4! % 5 2 3 ( 2880
h5 f iv (ξ ) # 32 88 & h5 f iv (ξ )
=0+ % − 24 + −12 ( = −
4! $ 5 3 ' 90
 Repeated use of Trapezoidal Rule
xN+1
I= ∫ ydx
x1

h!" f (x1 ) + f (x2 )#$ h!" f (x2 ) + f (x3 )#$ h!" f (x3 ) + f (x4 )#$
= + + +...
2 2 2
h!" f (xN−1 ) + f (xN )#$ h!" f (xN ) + f (xN+1 )#$
+ +
2 2
 Repeated use of Trapezoidal Rule
xN+1
I= ∫ ydx =−
3
(Δx) f # (ξ )
x1 12
h!" f (x1 ) + f (x2 )#$ h!" f (x2 ) + f (x3 )#$ h!" f (x3 ) + f (x4 )#$
= + + +...
2 2 2
h!" f (xN−1 ) + f (xN )#$ h!" f (xN ) + f (xN+1 )#$
+ +
2 2
" N %
# f (x1 ) + 2∑ f (xi ) + f (xN+1 )&
$ i=2 '
=h
2

# N &
$ f (x1 ) + 2∑ f (xi ) + f (xN+1 )'
% i=2 ( 1
= (b − a) 3
2N (b − a) 1
=− N × f ## ∝ 2
WIDTH AVERAGE  HEIGHT 12N 3 N
 Repeated use of Simpson`s 1/3rd Rule
xN+1 x3 x5 xN+1
I= ∫ ydx = ∫ ydx + ∫ ydx +... + ∫ ydx
x1 x1 x3 xN−1

h!" f (x1 ) + 4 f (x2 ) + f (x3 )#$ h!" f (x3 ) + 4 f (x4 ) + f (x5 )#$ h!" f (x5 ) + 4 f (x6 ) + f (x7 )#$
= + + +... 5
3 3 3 (Δx) f iv (ξ )
=−
h!" f (xN−3 ) + 4 f (xN−2 ) + f (xN−1 )#$ h!" f (xN−1 ) + 4 f (xN ) + f (xN+1 )#$ 90
+ +
3 3

h# N N−1 &
= % f (x1 ) + 4 ∑ f (xi ) + 2 ∑ f (xi ) + f (xN+1 )(
3$ i=2,4,6 i=1,3,5 '
5 N /2
( b − a)
# N N−1 & Error = −
90N 5
∑ (ξi )
f iV

% f (x1 ) + 4 ∑ f (xi ) + 2 ∑ f (xi ) + f (xN+1 )( i=1

$ i=2,4,6 i=1,3,5 ' 5
5
= (b− a)
3N =−
(b − a)
N× f iV
iV 1
∝ 44
55
WIDTH AVERAGE  HEIGHT 90N
180N N
N
 Repeated use of Simpson`s 3/8 th Rule
x N+1 x4 x7 x N+1
I= ∫ y dx = ∫ y dx + ∫ y dx +... + ∫ y dx
=−
5
3(Δx ) f iv (ξ )
x1 x1 x4 x N−2 80

3h !" f ( x1 ) + 3 f ( x2 ) + 3 f ( x3 ) + f ( x4 )#$ 3h !" f ( x4 ) + 3 f ( x5 ) + 3 f ( x6 ) + f ( x7 )#$ 3h !" f ( x7 ) + 3 f ( x8 ) + 3 f ( x9 ) + f ( x10 )#$

= + + +..
8 8 8
3h !" f ( x10 ) + 3 f ( x11 ) + 3 f ( x12 ) + f ( x13 )#$ 3h !" f ( x N −2 ) + 3 f ( x N −1 ) + 3 f ( x N ) + f ( x N+1 )#$
+ +.. +
8 8
5 N /3
3(b − a)
3h # N −1 N N −2 &
= % f ( x1 ) + 3 ∑ f ( xi ) + 3 ∑ f ( xi ) + 2 ∑ f ( xi ) + f ( x N+1 )(
Error = −
80N 5
∑ f (ξ )
iV
i
i=1
8$ i=2,5,8,.. i=3,6,9.. i=4,7,10... '
5

=−
( b − a) N × f iV ∝ 1
80N 5 N4
3 $ N −1 N N −2 '
= (b − a) × & f ( x1 ) + 3 ∑ f ( xi ) + 3 ∑ f ( xi ) + 2 ∑ f ( xi ) + f ( x N+1 ))
8N % i=2,5,8,.. i=3,6,9.. i=4,7,10... ( 11
AVERAGE&HEIGHT&
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING

CLL-­‐113  
Numerical  Integration  
Newton  Cotes  Formula:  Worked  Out  Examples

Prof.  Jayati  Sarkar

 Integration-Example
2 3 4 5
f (x) = 0.2 + 25x − 200x + 675x − 900x + 400x
4

3
0.8
I= ∫ f (x) dx

f(x)
2
0
1

0
0 0.2 0.4 0.6 0.8
Analytical Solution: x

0.8
25 2 200 3 675 4 900 5 400 6
I = 0.2x + x − x + x − x + x
2 3 4 5 6 0

I = 1.6405333
 Trapezoidal Rule
Analytical Solution: 1.6405333
Single Domain, N=1

h = 0.8
h
I = ( f (x1 ) + f (x2 ))
2
h
= ( f (0 ) + f (0.8))
2
0.8
= (0.2 + 0.232)
2
= 0.1728
 Trapezoidal Rule
Analytical Solution: 1.6405333
2 Domains, N=2:

h = 0.4

h
I = ( f (x1 ) + 2 × f (x2 ) + f (x3 ))
2
h
= ( f (0 ) + 2 f (0.4) + f (0.8))
2
0.4
= (0.2 + 2 × 2.456 + 0.232)
2
= 1.0688 x1 x2 x3
 Trapezoidal Rule
Analytical Solution: 1.6405333
3 Domains, N=3:

h = 0.8 3

h
I = ( f (x1 )+ 2 × f (x2 )+ 2 × f (x3 )+ f (x4 ))
2
h
= ( f (0)+ 2 f (0.2667)+ 2 f (0.5334)+ f (0.8))
2
0.2667
= (0.2 + 2 ×1.4329 + 2 × 3.4872 + 0.232)
2
= 1.3698 x1 x2 x3 x4
 Trapezoidal Rule
Analytical Solution: 1.6405333
4 Domains, N=4:

h = 0.8 4

h
I = ( f (x1 )+ 2 × f (x2 )+ 2 × f (x3 )+ 2 × f (x4 )+ f (x5 ))
2
h
= ( f (0)+ 2 f (0.2)+ 2 f (0.4)+ 2 f (0.6) + f (0.8))
2
0.2
= (0.2 + 2 ×1.288+ 2 × 2.456 + +2 × 3.464 + 0.232)
2
x1 x2 x3 x4 x5
= 1.4848
 Error Estimate in Trapezoidal Rule
3
2 3
f (x) = 0.2 + 25x − 200x + 675x − 900x + 400x 4 5
=−
(b− a)
× f ##
2
12N
f ! (x) = 25 − 400x + 2025x2 − 3600x3 + 2000x4

f !! (x) = −400 + 4050x −10800x2 + 8000x3

0.8

∫ f !!(x)dx 1 0.8 2 3
f !! (x) = 0

(0.8 − 0)
= ∫
0.8 0
(−400 + 4050x −10800x + 8000x )dx

0.8
2
1 " x 10800 3 8000 4 %
= $ −400x + 4050 − x + x '
0.8 # 2 3 4 &0
= −60
Error in Trapezoidal Rule =−
(b− a)
3

× f ##
2
Analytical Solution: 1.6405333 12N
Number  of   Actual Estimated  
Domains

1 =1.6405333-­‐0.1728   =-­‐(0.8)3×(-­‐60)/12×  (1)2  

=1.4677333   =1.31072  
Error  actual=1.4677333/1.6405333=89.47%   Errorestimatel=1.31072  /1.6405333=79.90%  

2 =1.6405333-­‐1.0688   =-­‐(0.8)3×(-­‐60)/12×  (2)2  

=0.5717333   =0.64  
Error  actual=0.5717333/1.6405333=34.85%   Errorestimatel=0.64  /1.6405333=39.01%  

3 =1.6405333-­‐1.3698   =-­‐(0.8)3×(-­‐60)/12×  (3)2  

=0.2707333   =0.28444  
Error  actual=0.2707333/1.6405333=16.50%   Errorestimatel=0.28444  /1.6405333=17.34%  

4 =1.6405333-­‐1.4848   =-­‐(0.8)3×(-­‐60)/12×  (4)2  

=0.1557333   =0.16  
Error  actual=0.1557333/1.6405333=9.49%   Errorestimatel=0.16  /1.6405333=9.75%  
 Error in Trapezoidal Rule
1.6

1.2

0.8
Error_actual
Error_Es3mated

0.4

0
0 0.2 0.4 0.6 0.8

h
Simpson`s 1/3rd Rule
h = 0.4
h
I = ( f (x1 ) + 4 × f (x2 ) + f (x3 ))
3
h
= ( f (0 ) + 4 × f (0.4) + f (0.8))
3
0.4 x1 x2 x3
= (0.2 + 4 × 2.456 + 0.232)
3
= 1.367
Et = 1.6405333 −1.367 = 0.2735333
ε t = 16.67%
Et,2T = 1.6405333 −1.0688 = 0.5717333
ε t,2T = 34.85% x1 x2 x3
 Error Estimate in Simpsons 1/3rd Rule
2 3 4 5 5
f (x) = 0.2 + 25x − 200x + 675x − 900x + 400x (b− a)
=− × f iv
f ! (x) = 25 − 400x + 2025x2 − 3600x3 + 2000x4 180N 4
f !! (x) = −400 + 4050x −10800x2 + 8000x3
f iii (x) = 4050 − 21600x + 24000x2 f iv (x) = −21600 + 48000x
0.8
5
∫ f iv (x)dx 0.8 (b− a)
iv 1 =− × f iv
f (x) = 0

(0.8 − 0)
=
0.8
∫ (−21600 + 48000x)dx 180N 4
0
0.85
2
0.8 =− 4
× (−2400 )
1 " x % 180 × 2
= $ −21600x + 48000 '
0.8 # 2 &0 = 0.27307
Et = 0.2735333
= −2400
Simpson`s 3/8th Rule
h = 0.2667
3h
I = ( f (x1 ) + 3 × f (x2 ) + 3 × f (x3 ) + f (x4 ))
8
3 0.8
= × ( f (0)+ 3× f (0.2667)+ 3× f (0.5334)+ f (0.8))
8 3
= 0.1(0.2 + 3×1.43287 + 3× 3.472 + 0.232 ) x1 x2 x3 x4

= 1.519221
Et = 1.6405333 −1.519221 = 0.1213123
ε t = 7.395%
Et,3T = 1.6405333 −1.3698 = 0.2707333
ε t,3T = 16.5%
x1 x2 x3 x4
 Error Estimate in Simpsons 3/8th Rule
2 3 4 5 5
f (x) = 0.2 + 25x − 200x + 675x − 900x + 400x (b− a)
=− × f iv
f ! (x) = 25 − 400x + 2025x2 − 3600x3 + 2000x4 80N 4
f !! (x) = −400 + 4050x −10800x2 + 8000x3
f iii (x) = 4050 − 21600x + 24000x2 f iv (x) = −21600 + 48000x
0.8
5
∫ f iv (x)dx 0.8 (b− a)
iv 1 =− × f iv
f (x) = 0

(0.8 − 0)
=
0.8
∫ (−21600 + 48000x)dx 80N 4
0
0.85
2
0.8 =− 4
× (−2400 )
1 " x % 80 × 3
= $ −21600x + 48000 '
0.8 # 2 &0 = 0.121363
Et = 0.1213123
= −2400
 Plug Flow Reactor
FA0 − FA
XA =
FA0

dFA = −FA0 dXA

−rA = kCAn
n n

FA (z) − FA (z+ dz) + rAdV = 0

dFA −rA = kC A0 (1− xA )
V= ∫ rA
XA_ EXIT
dxA FA0 dxA
−dFA + rAdV = 0 V = FA0 ∫ V= n ∫ n
−rA kCA0 0 (1− xA )
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING  

CLL-­‐113  
 
Numerical  Integration  
Newton  Cotes  Formula:  Worked  Out  Examples  

Prof.  Jaya*  Sarkar  

+  
 n   Order  of  Error  
 

Trapezoidal  Rule   3   2 !h$

2 I # & − I(h) 4  
"2%
2 2 −1
1/3rd  Simpsons   5   4 !h$
2 I # & − I(h) 6  
"2%
Rule  
2 4 −1
3/8th  Simpsons   5   6  
!h$
Rule   2 4 I # & − I(h)
"2%
 
2 4 −1
 True  value:  1.640533  

Segments   h   Integral   Richardson   Error(%)  

Extrapola*on  
(4/3)I(h/2)-­‐I(h)/3  
 
1   0.8   0.1728   1.3675   16.65%  
2   0.4   1.0688   1.6235   1.038%  
4   0.2   1.4848  
 Open Integration

BeDer  esEmate  by  passing  the  straight    

Trapezoidal  Rule   line  through  two  intermediate  points    
rather  than  the  two  extreme  points  
b  

a  
x  
Open Integration

∫ f (z)dz
−1
Gauss Legendre Quadrature Formula
Gauss Legendre Quadrature Formula

1  
Gauss Legendre Quadrature Formula

1  

" −1 % " 1 %
−2 $ ' 2$ '
−2u2 # 3& 2u1 # 3&
a1 = = =1 a2 = = =1
(u1 − u2 ) "$ 1 − "$ −1 %'%' (u1 − u2 ) "$ 1 − "$ −1 %'%'
# 3 # 3 && # 3 # 3 &&
Gauss Legendre Quadrature Formula
 Gauss Legendre Quadrature -Example
0.8

∫ f (x)dx
0

x=
( 0.8 − 0 )
Z+
" 0.8 + 0 %
x = 0.4Z + 0.4
0.8 1 1

2
$
# 2
'
&
I= ∫ f (x)dx = ∫ 0.4 × f (z)dz = ∫ g(z)dz
0 −1 −1
2
f (z) = 0.2 + 25 ( 0.4 + 0.4Z ) − 200 ( 0.4 + 0.4Z )
3 4 5
+675 ( 0.4 + 0.4Z ) − 900 ( 0.4 + 0.4Z ) + 400 ( 0.4 + 0.4Z )

g(−0.57735) = 0.516740523 g(0.57735) = 1.3058376

# 1 & # 1 & =1.6405333-­‐1.822578123=-­‐0.182044823  

I GQ = 1× g % − ( +1× g % ( = 1.822578123 Error  actual=-­‐0.182044823/1.6405333=-­‐11.1%  
$ 3' $ 3'
 NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING

CLL-­‐113  
Ordinary  Differential  Equations,  Initial  Value  
Problem  
Runge  Kutta  Methods,  Euler  Method
Prof.  Jayati  Sarkar
 Reactions  in  Series:

Model  Equations
These  Set  of  Model  Equations  can  be  written  in  standard  form  for  1st  order  ODE

dyi
= fi (y1, y2 ,....yN )
dt
for  i=1,2,…N

The  initial  conditions  of  the  above  equations  are:

t = 0, yi (0) = αi .......(known)
 How  to  recast  in  standard  form?
We  introduce  a  new  dependent  variable

+1
Q(t)
Convert  to  Standard  Form

dy2
=1
dt

dy1 Ah Q (y2 )
= (y1 − T0 )−
dt VρCp ρCp
 How  to  recast  in  standard  form?

They  will  be  recast  in  the  form

yi+1 = yi + φh
Euler Explicit Method

yi+1 = yi + hf (xi , yi )
Euler Implicit Method

yi+1 = yi + hf (xi+1, yi+1 )

xi xi+1 x
Semi Implicit (Crank Nicholson) Method

h!
yi+1 = yi + " f (xi , yi )+ f (xi+1, yi+1 )#$
2
y

xi xi+1 x
 7"

6"

5"

4"

y"
3"

Explicit  Euler: 2"

1"
ytrue"
0"
0" 0.5" 1" 1.5" 2"yEuler" 2.5"

True  Solution  @  x=0.5 x" y"Euler_local"

Error@  x=0.5
Second  Step: True  Solution  @  x=1.0

Explicit  Euler: y(1) = ytrue (0.5)+ hf (0.5, ytrue (0.5))

= 3.21875 + 0.5× f (0.5, 3.211875)
= 3.21875 + 0.5× (1.25)
= 3.84375

Local  Error@  x=1.0

Global  Error@  x=1.0

5.875-­‐3=2.875 3.84375-­‐3=0.84375

εt = 0.84375×100% 3 = 95.8% = 28.125%

εt = 2.875×100% 3 = 95.8% 7"

6"

5"

4"
y"

3"

2"

1"
ytrue"
0"
0" 0.5" 1" 1.5" 2"yEuler" 2.5"

x" y"Euler_local"
 7"

yi+1 = yi + hf (xi , yi )
6"

5"

4"

y"
3"

2"

1"
ytrue"
0"
0" 0.5" 1" 1.5" 2"yEuler" 2.5"

x" y"Euler_local"

h 0.5          

f(x)   yEul
x =dy/dx ytrue er y  Euler_local Global_Error Local_error
0 8.5 1 1      

0.5 1.25 3.21875 5.25 5.25 -­‐63.10679612 -­‐63.10679612

1 -­‐1.5 3 5.87 3.84375 -­‐95.83333333 -­‐28.125
5
5.12
1.5 -­‐1.25 2.21875 5 2.25 -­‐130.9859155 -­‐1.408450704
2 0.5 2 4.5 1.59375 -­‐125 20.3125
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING

CLL-­‐113  
Ordinary  Differential  Equations  
Runge  Kutta  Methods,  Euler  Method

Prof.  Jayati  Sarkar

Round  off  Errors:

Truncation  /Discretization  Errors:

❖ Local  Truncation  Error

❖ Propagated  Truncation  Error
The  sum  of  the  two  is  the  Global  Truncation  Error
Global  Truncation  Error  ~  O(h)
 dy
=
dx
Error@  x=0.5
 Heun's  Method

x f(x) ytrue k1 k2 S yH Error

_t
0 8.50 1.00 1.00  
yHeun (i +1) = y (i )+ hs!" y (i ), x (i )#$ 0.5 1.25 3.22 8.50 1.25 4.88 3.44 -­‐6.80
yHeun (i +1) = y (i )+ h [w1k1 + w2 k2 ] 1 -­‐1.50 3.00 1.25 -­‐1.50 -­‐0.13 3.38 -­‐12.50
1.5 -­‐1.25 2.22 -­‐1.50 -­‐1.25 -­‐1.38 2.69 -­‐21.13
k1 = f !" y (i ), x (i )#$ 2 0.50 2.00 -­‐1.25 0.50 -­‐0.38 2.50 -­‐25.00
k2 = f !" y (i )+ hk1, x (i )+ h#$
 Mid  Point  Method

x f(x yt k1 k2 S yMi Error

) ru d _t
yMid (i +1) = y(i )+ hs!" y (i ), x (i )#$ 0 8.50 1.00 1.00  

y (i +1) = y (i )+ h [w1k1 + w2 k2 ] 0.5 1.25 3.22 8.50 4.22 4.22 3.11 3.40
Mid
1 -­‐1.50 3.00 1.25 -­‐0.59 -­‐0.59 2.81 6.25
k1 = f !" y (i ), x (i )#$ 1.5 -­‐1.25 2.22 -­‐1.50 -­‐1.66 -­‐1.66 1.98 10.56
k2 = f !" y (i )+ 0.5hk1, x (i )+ 0.5h#$ 2 0.50 2.00 -­‐1.25 -­‐0.47 -­‐0.47 1.75 12.50
 Ralston  Method
-­‐

¾  h

x f(x) yt k k2 S yR Err
ru 1 or_
yR (i +1) = y (i )+ hs!" y (i ), x (i )#$ 0 8.50 1.00 1.00  

yR (i +1) = y(i )+ h [w1k1 + w2 k2 ] 0.5 1.25 3.22 8.50 2.58 4.53 3.27 -­‐1.51
w1 = 0.33, w2 = 0.67 1 -­‐1.50 3.00 1.25 -­‐1.15 -­‐0.36 3.09 -­‐2.92
k1 = f !" y (i ), x (i )#$ 1.5 -­‐1.25 2.22 -­‐1.50 -­‐1.51 -­‐1.51 2.33 -­‐5.18
k2 = f !" y (i )+ 0.75hk1, x (i )+ 0.75h#$ 2 0.50 2.00 -­‐1.25 0.00 -­‐0.41 2.13 -­‐6.44
 Comparison

6.00

4.50

ytrue
3.00
yH
y

yMid
yR
1.50 YEuler

0.00
0 1 1 2 2

x
 ]

yi+1 = yi + h [w1 fi ] + h!"w2 fi + w2 qhk1 f yi + w2 phfti + O (h2 )#$

 hf (yi , ti )
qhf
ph

h
t1 t2
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING

CLL-­‐113  
Ordinary  Differential  Equations  
Runge  Kutta  Methods,  Euler  Method

Prof.  Jayati  Sarkar

RK-n Methods
n
yi+1 = yi + h∑ wmkm
m=1

k1 = f (yi , ti )
km = f (yi + h"#qm,1k1 + qm,2 k2 +...... + qm,m−1km−1 $%, (ti + pmh))

Accuracies

Eulers RK-­‐2 RK-­‐3   RK-­‐4   RK-­‐5

O(h2) O(h3)   O(h4)   O(h5)   O(h5)  

RK-4 Methods
yi+1 = yi + hS(yi , ti )

S(yi , ti ) = w1k1 + w2 k2 + w3k3 + w4 k4

p2 q21      

k1 = f (yi , ti ) p3 q31 q32    

p4 q41 q42 q43  

k2 = f (yi + hq2,1k1, (ti + p2 h)) w1 w2 w3 w4

k3 = f (yi + hq3,1k1 + hq3,2 k2 , (ti + p3h))

k4 = f (yi + hq4,1k1 + hq4,2 k2 + hq4,3k3, (ti + p4 h))

Classical RK-4 Methods
Popular RK-4 Methods p2 q21      

p3 q31 q32  
Classic RK-4  

p4 q41 q42 q43  

0.5 0.5      
w1 w2 w3 w4
0.5 0 0.5    

1 0 0 1  
 1/6  1/3  1/3  1/6 RK-Fehlberg Method

RK-Gill 1/4 1/4      

3/8 3/32 9/32    

0.5 0.5       12/13 1932/2197 -­‐7200/2197 7296/219  

0.5 (√2-­‐1)/2 (2-­‐√2)/2   7  
 
 25/216  1408/2565  -­‐1/5
1 0 -­‐1/√2 (2+√2)/2   2197/410
 1/6  (2-­‐√2)/6 (2+√2)/6    1/6
Multi Step Method
Heun's Method
k1 = f (yi , xi )
Predictor: Euler Method
0 2
y = yi + hk1
i+1 +O(h )
Corrector: Trapezoidal Rule Corrector repeated multiple times
0 0
k = f (y , xi+1 ) ( )
m
2 i+1 k = f yi+1, xi+1
m
2

1 h m+1 h
yi+1 = yi + (k1 + k20 ) yi+1 = yi + (k1 + k2m ) m= 1 to M
2 m+1 m m 2
yi+1 = yi+1 = yi
h!
yi+1 = yi + " f (yi , ti ) + f yi+1, ti+1 #$
( ) CRANK NICHOLSON (SEMI-IMPLICIT METHOD)
+O(h3 )
2
Improved Multi Step Method
non-starting Heun's Method

3
+O(h )

3
+O(h )

use original Heun's Method at i=0

Stability of Euler Method
−λt
Analytical solution y=e

0.75

0.5
y

0.25

0
0 1.25 2.5 3.75 5

t
Stability of Euler Method
−λt
Analytical solution y=e

Euler Method

A=
To enforce the stability we want the error at tn+1 to be smaller than tn
Stability  Envelope
1.5

0.75
Globally  Stable
0
A= =-­‐ve Converges  but    
oscillatory  after  
-­‐0.75
A  this

-­‐1.5

-­‐2.25 Converges    
Euler  (-­‐inf) Implicit  (0) before  this  
Crank  Nicholson  (-­‐1)
Oscillatory    
-­‐3
and  diverges
0 1 2 3 4

hlambda
 Richardson Extrapolation

# h&
Δ = yi+1 − yi+1 % (
same $2'
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING

CLL-­‐113  
Ordinary  Differential  Equations  
Runge  Kutta  Methods,  Euler  Method

Prof.  Jayati  Sarkar

 Richardson Extrapolation

# h&
Δ = yi+1 − yi+1 % (
same $2'
 3

3
 Concept of Adaptive Step-Size
Step-Size? Value of λ?

Euler Stability Envelope

Linearize the Equation
 Adaptive Step-Size

If Δnew~εtol 1/n
# ε tol &
hnew = % ( h
$ Δ '
If εtol~10-2, n=2 0.5 0.5
! 0.01 $ ! 0.01 $
Δ1=0.0625 hnew = # & h = 0.4h Δ2=0.00625 hnew = # & h = 1.265h
" 0.0625 % " 0.00625 %
a0 = 0 Explicit Implicit
Backward Difference Formula

Implicit Euler
 α (α +1) 2 α (α +1).. (α + n− 2)
Pn−1 (i ) = fi +α∇fi + ∇ fi +... + ∇ n−1 f
2! 2!
t − ti
α=
h
ti+1 1 1
yi+1 − yi = ∫ Pn−1 (i )dt = ∫ Pn−1 (i )hdα = h ∫ Pn−1 (i )dα
ti 0 0

1
yi+1 = yi + h ∫ Pn−1 (i )dα
0
O(hn)

LTE:::O(hn+1)
GTE:::O(hn)
 α (α +1) 2 α (α +1).. (α + n− 2)
Pn−1 (i ) = fi +α∇fi + ∇ fi +... + ∇ n−1 f
2! 2!
t − ti
α=
h ti+1 1 1
yi+1 − yi = ∫ Pn−1 (i )dt = ∫ Pn−1 (i )hdα = h ∫ Pn−1 (i )dα
ti 0 0

1 1
❖ n=1
yi+1 = yi + h ∫ P0 (i )dα = yi + h ∫ fi dα = yi + hfi
0 0
LTE:::O(h2) GTE:::O(h1)
Adams  Bashforth  1st  Order  M
1
ethod
1 1 $ 2 &
❖ n=2 yi+1 = yi + h ∫ P1 (i )dα = yi + h ∫ $% f (i )+α∇fi &'dα = yi + h)α fi + ( fi − fi−1 )* = h 3 fi − 1 fi−1
α $ &
)% *'
0 0 % 2 '0 2 2
LTE:::O(h3) GTE:::O(h2) Adams  Bashforth  2ndOrder  Method
1 1$ α (α +1) 2 ' $ 23 4 5 '
❖ n=3 yi+1 = yi + h ∫ P2 (i )dα = yi + h ∫ & f (i )+α∇fi + ∇ f )dα = yi + h& fi − fi−1 + fi−2 )
0 0 % 2! ( %12 3 12 (
4 3 rd
@i = 0 y1 = y0 + a1 f0 + a2 f−1 + a3 f−2
@i =1 y2 = y1 + a1 f1 + a2 f0 + a3 f−1
@i = 2 y3 = y2 + a1 f2 + a2 f1 + a3 f0
AM    n  method
NUMERICAL  METHODS  IN  
CHEMICAL  ENGINEERING

CLL-­‐113  
Ordinary  Differential  Equations  
Boundary  Value  Problem  
Shooting  method,  Finite  Difference  Method
Prof.  Jayati  Sarkar
ODE-IVP : Auxiliary condition at 1 point y

y0

ODE-BVP : Auxiliary condition at 2 points

x
CA (x = L ) = 0
Assume:
dy1
= y2
dx
dy2
+ py2 + q = 0
dx
dy2
= −py2 − q
dx

d ! y1 $ ! y2 $
# &=# &
dx " y2 % "−py2 − q%
 Ta=30°C  
d 2T
2
= 0.01(T − Ta )
dx
T(0)=30°C T(L)=80°C

y1 = T

Actual BC:
y1 = 30@x = 0
d ! y1 $ ! y2 $
# &=# & y1 = 80@x = L
dx " y2 % "0.01(y1 − 30)%
Modified BC:
y1 = 30@x = 0
y2 = 0@x = 0
  Fundamentals   of  CFD  (CLL:113)
NUMERICAL*METHODS*IN*
CHEMICAL*ENGINEERING*
Dr.  Jayati   Sarkar  
CLL#113&
Email:  jayati@chemical.iitd.ac.in
&
Lecture 8
Ordinary&Differential&Equations&
PDE&
Department of Chemical Engineering, IIT Delhi,
NewProf.&Jaya*&Sarkar&
Delhi, India
 Classification  of  Partial  Differential  Equations

• A general second order PDE

aφxx + bφxy + cφyy + dφx + eφy + fφ + g = 0
where
a,b,c,d ,e, f , g = F(x, y) ≠ F(φ)

• The behavior of the PDE depends on the sign of

the discriminant
if D< 0 D =b2 −4ac if D >0

if D =0 PDE is called hyperbolic

PDE is called elliptic
PDE is called parabolic
 Elliptic  Partial  Differential  Equation
Steady state Heat Conduction in a
1D slab
∂ " ∂T %
$k ' = 0
∂x # ∂x &
with T (0 ) = T 0 , T ( L ) = T L

For const k, the soln is given by

(TL − T0 )
T (x) = T0 + x
T
L
Important lessons to learn!
T0
1.The temperature at any point x is influenced by
TL temperatures at both boundaries.
2.In the absence of source terms, T (x) is
x bounded
by the boundary temperatures.
 Parabolic  Partial  Differential  Equation

t >0
T0 T0

0 L/2 x

T⎯⎯→∝

Unsteady state Heat Conduction in a 1D slab

∂T ∂2T
=α 2 α = k / ρCp
with ∂t ∂x
i
T (x,0) = T (x)
T (0,t )= T0 T
(L,t )= T0
 Parabolic  Partial  Differential  Equation  (Contd.)
∞ −αn2π 2
t
Solution is given by $ nπ x ' L2
T (x, t ) = T0 + ∑ Bn sin & )e
n=1
% L (

2 L $ nπ x '
where Bn = ∫ (Ti (x) − To )sin & ) dx
L0 % L (
Important lessons to learn!
The boundary temperature T0 influences the temperature
T(x,t) at every point in the domain, just as with elliptic
PDE’s
The initial conditions only affect future temperatures, not past
temperatures.
 Parabolic  Partial  Differential  Equation  (Contd.)
∞ −αn2π 2
t
Solution is given by $ nπ x ' L2
T (x, t ) = T0 + ∑ Bn sin & )e
n=1
% L (
2 L $ nπ x '
where Bn = ∫ (Ti (x) − To )sin & ) dx
Important lessons to learn! L0 % L (
The initial conditions influence the temperature at every point
in the domain for all future times. The amount of influence
decreases with time, and may affect different spatial points to
different degrees.
A steady state is reached for t ∞. Here, the solution becomes
independent of Ti
(x ,0) . It also recovers its elliptic spatial behavior.
The temperature is bounded by its initial and boundary
conditions in the absence of source term
 Hyperbolic  Partial  Differential  Equations

• 1D flow of fluid in a
channel

∂ (ρCPT ) ∂ (ρCPUT )
+ =0
∂t ∂t
Temperature Variation with time
with
Important lessons to learn! T (x,0) = Ti (x)
1.The upstream boundary conditions and
T (x ≤ 0,t )= T0
not the downstream condition affects the
solution of the domain. Solution is given by
2. The inlet boundary condition T ( x , t ) = T ((x −Ut ),0 )
propagates with a finite speed U.
or x
3. The inlet boundary condition T (x,t ) = Ti for t <
U
is not felt at point x until t=x/U. x
=T 0 for t ≥
U
 Laplace  Equation  (Elliptic)
• Heat  Transfer  in  2D

8
Laplace  Equation  (Elliptic)-­‐Dirichlet  Boundary  Condition
• Heat  Transfer  in  2D

• Gauss  Seidel

9
Laplace  Equation  (Elliptic)

10
• Solution  after  9th  iteration
Laplace  Equation  (Elliptic)
• Post  Processing(Heat  Fluxes)

11
• Solution  after  9th  iteration Laplace  Equation  (Elliptic)
• Post  Processing(Heat  Fluxes)

12
Laplace  Equation  (Elliptic)-­‐Neumann  Boundary  Condition

13
Laplace  Equation  (Elliptic)-­‐Neumann  Boundary  Condition

14
Unsteady  Heat  Conduction  Equation  (Parabolic)

15
Unsteady  Heat  Conduction  Equation  (Parabolic)

16
Unsteady  Heat  Conduction  Equation  (Parabolic)

17
NUMERICAL*METHODS*IN*
CHEMICAL*ENGINEERING*

CLL#113&
&
Ordinary&Differential&Equations&
PDE&

Prof.&Jaya*&Sarkar&

18
19