M.Sc.-Mathematics-Numerical-Analysis-munotes

Page 1

1
1
NUMERICAL I NTEGRATIO N–I
Chapter Structure
1.1 Objective
1.2 Introduction
1.3 Newton –cote’s quadrature formula
1.4 Trapezoidal Rule
1.5 Error in Trapezoidal rule
1.6 Simpson’s one third rule
1.7 Error in Simpson’s 1/3rdrule
1.8 Simpson’s 3/8thrule
1.9 Error in Simpson’s (3/8)thrule
1.10 Review
1.11 Unit End Exercise
1.1 OBJECTIVE
*After going through this chapter you will be able to :
*Solve any definite integral in interval [a,b].
*Solve unknown function integral by numerical method .
*Learn different technique for solving definite integral by numerical
method .
*Solve definite integral very efficient way .
1.2 INTRODUCTIO N
Differentiation and integration are basic Mathematical operation
with a wide range of applications in many areas of science .
In this we are going to explose various ways for approximating the
integral of a function over a given domain .These are various reasing as of
why such approximations can be useful .
i)Not every function can be analytically integrated .
e.g.:1
3011dxxor1
0tanxdxx.
ii)Even if a closed integration formula exists it might still not be the most
efficient way of calculating the integral .munotes.in

Page 2

2


iii) It can happen that we need to integrate an unknown function in which
only some examples of the functions are known.

A simpler appro ach for approximating the value of b
afx d x would be
to complete the product of the value of the function at one of the end points of
the interval by the length of the interval .

In case we choose the end -point where the function is eva luated to be
x = a we obtained
 ~fb
afx d x aba
This approximation is called rectangular method or the rectangular
quadrature . The points 01,, . . . . . ,nxx x that are used in the quadrature formula are
called quadrature points.
















We may approximate the integrate by a linear curve i.e. ( y = a + bx) or by
a second degree curve
i.e. by ( y = a + bx + cx2) …. and soon.
Here we will learn till three degree polynomial by numerical integration
using :
i) Trapezoidal rule ( linear curved )
ii) Simpson’s one -third rule (second degree)
iii) Simpson’s third -eight rule (third degree)




y
x x0x0 + nx0 + ihx0 + nhy1 y2 yi yny=f(x)
munotes.in

Page 3

3


1.3 Newton – cote’s quadrature formula :

The newton cote’s is an extremely useful and straight forward family of
numerical integration techniques using this w e can derive the all three numerical
integration formula which we need.

Let b
aIf x d x

Let us divide the interval (a, b) into n – sub intervals of width h . as shown
in figure (1) so that x0 = a, x1 = x0 + h, x2 = x0 + 2h ….. + xn = x0 + nh = b.

Then 0
0xn h
xIf x d x
taking 0xx p h d xh d p
0
00xn h
xIf x p h d p h

= By newton’s forward interpolation formula.
2
3
00
011 2...2! 3!n
w
opp y pp phy p y y d p      

Integrating term by term we obtain
2
23
00 0 023 2.....21 2 2nn n n nnh y y y yh      

0
02
23
00 0 023 2....21 2 2xn h
xnn n n nfx d x n hy y y yn       

Which is known as the newton –cote’s quadrature formula using this we
can deduce the following rule taking n = 1, 2, 3, …..

1.4 Trapezoidal Rule :

Taking n = 1 in newton cote’s formula at the curve ( x0, y0) to ( x1, y1) as a
straight line.

i.e. polynomial of order are so that higher order difference become zero
we get,



munotes.in

Page 4

4

0001
2xh
xfx d x hy y 
012hyy
Similarly 
02122xh
xhfxd x y y
.
.
.
0
012xn hnnxhfxd x y y


Adding these integrals we get
 0
001 2 1 2. . .2xn h
nnxhfxd x y y y y y
     

This is known as the trapezoidal rule.
y
xx0 = ax0 +n x0 +2n x0 + nh y = f x ()


Working of trapezoidal rule. For b
afxd x
Let y = f (x) take n value of x in interval [ a, b]

012
012.....
.....nnxx x x xyy y y y
  01 2 1 2. . . . .2b
nnahfx d x y y y y y    
where bahn



munotes.in

Page 5

5

Example :
1) Evaluate 6
301
1dxx taking h = 1 by using trapezoidal rules.

Solution :
Divide the interval (0, 6) into six parts each of width h = 1. The value of
311fxx are given below :
01 2 3 4 5 601 2 3 4 5 610 . 5 0 . 1 1 0 . 0 3 5 7 0 . 0 1 5 4 0 . 0 0 7 9 0 . 0 0 4 6x
yyy y y y y y
by trapezoidal rule
  01 2 1 2. . . . .2b
nnahfx d x y y y y y    
 110 . 0 0 4 6 2 0 . 50 . 1 10 . 0 3 5 70 . 0 1 5 40 . 0 0 7 92    
11.0046 2 0.6692
 11.0046 1.3382
12.34262
1.1713

Ex.2 :
Evaluating 1
0sin cosxx d x taking 5 sub-intervals using trapezoidal rules.

Soln :

Here n = 5, a = 0, b = 1
100.25bahn
00 . 2 0 . 4 0 . 6 0 . 8 111 . 0 8 5 71 . 1 4 4 81 . 1 7 9 01 . 1 8 9 11 . 1 7 5 5x
y

by trapezoidal rule
  01 2 1 2. . . . .2b
nnahfx d x y y y y y     

 0.211 . 1 7 5 5 2 1 . 0 8 5 71 . 1 4 4 81 . 7 9 01 . 1 8 9 12  
0.1 2.1755 2 4.5986 munotes.in

Page 6

6

0.1 2.1755 9.1972
1.13727

Ex.3 :
Evaluate 1201dxx by trapezoidal rule .
6 – coordinate, also determine absolutely error.

Soln :
Here a = 0 , b = 1 , n = 5

100.25bahn
211fxx


00 . 2 0 . 4 0 . 6 0 . 8 110 . 9 6 1 50 . 8 6 2 10 . 7 3 5 30 . 6 0 9 80 . 5x
y

by trapezo idal rule
  01 2 1 2. . . . .2b
nnahfx d x y y y y y     

 0.210 . 5 2 0 . 9 6 1 50 . 8 6 2 10 . 7 3 5 30 . 6 0 9 82  
0.1 1.5 2 3.1687
0.1 1.5 6.3374
0.78374.

 11 111 1
2 0 001tan tan 1 tan 014fx d x d x xx    
Absolute error 0.783744
0.00165.








munotes.in

Page 7

7

1.5 Error in Trapezoidal rule :

















 The expression of y = f (x) about x = x0 by
Taylor’s theorem is given by

2
0 11 1
00 0 0.....2!xxyf x f x x x f x fx    
 11
000 11 1
00 0 1.....2!xx
xxxxfx d x fx x x f x f x   
1
01023
00 11 1
000.....26x
xxxxx xxfx d x x fx f x f x  
23
10 10 11 1
00 0 0.....26xx xxxx f x f x f x   

Put 10xxh

1
023
11 1
00.....26x
xhhfx d x h fx f x f x   

Area of trapezoidal in the interval 01,xx is
10 12hAy y
by trapezoidal rule
 10 02hAf x f x h 
by taylor’s Thm y
xx=x0x=x1x=by = f x1 ()
y = f x ()
A1A2
munotes.in

Page 8

8

2
11 1
10 0 0 0.....22 !hhAf x f x h f x f x   
23
11 1
10 0 0.....24hhAh f x f x f x  


Let 1be the error in the interval 01,xx
1
011x
xfxd x A
31 1
011.....64hf x 
31 110–112hf x
Similarly :
 2
131 112 1–112x
xfxd x A hf x 
 3
231 123 2112x
xfxd x A hf x 
.
.
.
 
131 11112n
nx
nn nxfxd x A hf x
 
The total error
1.....n  
 3
11 11 11 11
012 1 .....12nhfx fx fx fx    
Choose 11 11 11 11
10 1 2max , , , .....ff x f x f x
3
11
1112hfn
3
111112hb afh
21 11112bahf
Which is the required error in trapezoid e rule.

Practice Problem :

Evaluate the following integral by using trapezoidal rule.
i) 22
0xd x with 0n Ans. : 2.68
munotes.in

Page 9

9

ii) 5
045dxx with 10n Ans. : 0.4055

iii) 2 0.6
0xed x with 6n Ans. : 0.5357

1.6 Simpson’s one third rule :

While taking n = 2 in newton cote’s formula we get
  0
001 3 1 2 4 2 4. . . . . 2 . . . . .3nhx
nn nxhfx d x y y y y y y y y
         

This is known as the Simpson’s one third rule.

In Simpson’s rule the given interval must be divided into even number of
equal sub -interval since we find the area of two strips at 0 time.
Working of Simpson’s one third rule for b
afxd x.
Let yf x, taking n – value of x in [a, b].
0123
0123..... .....
..... .....nnxx x x x xyy y y y y
  0
001 3 1 2 4 2 4. . . . . 2 . . . . .3nhx
nn nxhfx d x y y y y y y y y
         

Where – / hba n

Ex. 1 :

Calculate upto 4 decimal places 8
2 316dxdxxx by us ing Simpson’s
(1/3)th rule taking n = 5.
Soln :
Here 21
16fxxx
 a = 3, b = 8, n = 5,
8315bahn

01234 534567 80.1601 0.1443 0.1348 0.1291 0.1259 0.125x
yyyyyy y

using Simpson’s (1/3)rd rule. munotes.in

Page 10

10


  01 13 1 24 2 4. . . . . 2 . . . . .3b
nnahfx d x y y y y y y y y         


   10.160 0.125 4 0.1443 0.1291 2 0.1348 0.12593  
 10.2851 1.0936 0.52863 
11.90733 0.6357

Ex. 2 :

Evaluate 2
0sin
54 c o sxdxx
 by taking 5 ordinates by Simpson’s (1/3)rd rule.


Soln :

We divide the interval 0,into 4 - equal part i.e. (5 ordinates).
044bahn

01 2 33 0using calculate42 4in degree mode00 . 0 6 3 9 0 . 2 0 . 2 3 0 2 0
4x
y
yyy y

by Simpson’s (1/3)rd rule.

   01 3 2 44. . . . . 2 . . . . .3b
nahfxd x y y y y y y       
 400 4 0 . 0 6 3 90 . 2 3 0 2 2 0 . 23
  
 1.1764 0.412
0.4127.




munotes.in

Page 11

11


1.7 Error in Simpson’s 1/3rd rule :


















Let yf x be a continuous functions and continuous derivation of all
orders in [a, b]. Divide the interval [a, b] into n equal sub -intervals by the points
01,, . . . . . , xnax x b & 101.2 .....xxi h i n  .
The T aylor’s expression of yf x atleast 0xx is
23
00 11 11 1
00 0 0 0.....2! 3!xx xxyf x f x x x f x fx fx    

 22
0023
00 11 11 1 1
00 0 0 0 .....2! 3!xx
xxxx xxfx d x fx x x f x f x f x d x   


2023 4
000 11 11 1 1
000 02! 3! 4!xxxx xx xxxf x f x f x f x  

23
20 20 11 1
20 0 0 0.....2! 3!xx xxxx f x f x f x   
Put 202xx h
  21 3 1 1 4 1 100 0 2422233hf x h f x h f x h f x   …………. (I)

by Simpson’s 1/3rd rule.
 10 1 243hAy y y  xixi + 1xi + 2fi + 1 fi + 2fi
munotes.in

Page 12

12

 00 0 423hfx fx h fx h   

Using taylor’s thm.

23
11 11 1
10 0 0 0 0223
11 1 1 1
00 0 04. . . . .32 6
482. . . . .26hh hAf x f x h f x f x f x
hhfx h f x f x f x       
2
11 1 3 1 1
00 0 0866 232hhfx h f x f x h f x  

  21 3 1 1 41 1 1
00 0 04222 . . . . .33hf x h f x h f x h f x    (I)


Let 1 be the error in the interval 01,xx
2
011x
xEf x d x A
05I V 45.....15 15xhf 
IV0590xhf

Similarly


 4
2
25
32
15222 90.
.
.90n
nx
x
nxIV
nnxIV
xhEf x d x A f
Ehfxd x A f x

  


 The total error
12 1.....nEE E E
 5
02 2 .....90II IV IV
nhfx f x f x  

Choose  1 2 02ax , , .....,IV
nIV IV IV
x xxfM f f f munotes.in

Page 13

13

590 2IV hnEf
590IV bahEfh
4
0180IVnbaEh f x h  
Which is the required even in compos ited Simpson’s 1/3rd rule.

1.8 Simpson’s 3/8th rule :

While taking n = 3 in newton’s take formula we set all differences uppers
then 3rd forward difference will becomes zero.
   3
001 0 2 1 0 3 2 1 09332 3 3 . . . . .28x
xfx d x h y y y y y y y y y yh       


   01 02 1 0 3 2 0381 2 6 2 33 . . . . .8hyy yy y y y y y y          
 012 33338hyyy y  
Similarly
 6
3345 03338x
xhfx d x y y y y  
.
.
.
 
33213338n
nxnnn nxhfx d x y y y y  

[where n is divisible by 3]

We add al l above integral

  36
03 3.....n
nbx x x
ax x xfx d x fx d x fx d x fx d x  

   012 3 345 6 3 2 1333 33 . . . . . 3 38nnn nhyyy y yyy y y y y y            

  01 2 4 5 2 1 3 6 9 333( ..... 2 .....8nn n nhyy yyyy y y yyy y              


Worki ng of Simpson’s (3/8)th rule for b
afx d x munotes.in

Page 14

14

Let yf x take n - value of x in [a, b]
012
012..... .....
..... .....nnxx x x xyy y y y
   01 2 4 5 3 633. . . . . 2 . . . . .8b
nahfxd x y y y y y y y y         
where bahn.

Ex 1 :

Evaluate 5.2
4logxed x by Simpson’s (3/8)th rule taking n = 6 sub interval.
Soln :
Dividing the interval [4, 5.2] in six equal part taking h = 1 we get fx
value logxfx e are given below.
5.2 4 1.20.266bahn  
44 . 2 4 . 4 4 . 6 4 . 8 5 . 0 5 . 21.3863 1.4351 1.4810 1.5260 1.5686 1.6094 1.6484x
y

by using Simpson’s (3/8)th rule.

   01 2 4 5 3 633. . . . . 2 . . . . .8b
nahfx d x y y y y y y y y         

  30 . 21.3863 1.6486 3 1.4351 1.4816 1.5886 1.6094 2 1.52608     
0.075 3.0349 3 6.1147 2 1.5260 
0.075 3.0349 18.3441 3.052 
0.075 24.431
1.8323

Ex 2 :

Evaluate 2 0.6
0xed xtaking n = 6 by Simpson’s (3/8)th rule.





munotes.in

Page 15

15

Soln :

Here 200 . 66xfx e a b n  
0.6 00.16bahn 

01 2 3 4 5 600 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 610 . 9 9 0 . 9 6 0 8 0 . 9 1 3 9 0 . 8 5 2 1 0 . 7 7 8 8 0 . 6 9 7 6x
yyy y y y y y

by using Simpson’s (3/8)th rule.
   01 2 3 4 3 633. . . . . 2 . . . . .8b
nahfx d x y y y y y y y y         


  30 . 110 . 6 9 7 6 3 0 . 9 90 . 9 6 0 80 . 8 5 2 10 . 7 7 8 8 2 0 . 9 1 3 98     
0.0375 1.6976 10.745 1.8278 
0.0375 14.2705
0.5351.

1.9 Error in Simpson’s (3/8)th rule :

Let yf x be a c ontinu ous function and hence continuous derivatives
of all under in [a, b]. Divide the interval [a, b] into n sub interval by the point
012,,, . . . . . , xnbax x x & 01, 2, .....,ixx i h i n 

The taylor’s expansion of yf x about 0xx is

2
0
00 0 0
34
00
002!
.....3! 4!II I
III IVxxyf x f x x x f x fx
xx xxfx fx    

 
33
0000 0
23 4
00 0
00[.....]2! 3! 4!xxI
xx
II III IVfx d x fx x x f x
xx xx xxfx f x f x d x  
  

3023 4 5
00 0 0
00 0 0 0.....262 4 1 2 0xI II III IVxxx xx xx xxxf x f x f x f x f x      munotes.in

Page 16

16


23
30 3030 0 045
30 30
0026.......24 120II I
a
III IVxx xxxx f x f x fxxx xxfx f x I I  

Put 303xx h

   23 4 500 0 0 092 7 8 1 2 4 3326 2 4 1 2 0I II IV IVhf x h f x h f x h f x h f x   
…. (I)
By Simpson’s 3/8th rule

 10 1 2 33338hAy y y y  
 00 0 0333 238hfx fx h fx h fx h     

Using Taylor’s thm


   
   23 4
11 11 1 1
10 0 0 0 0 0
23
2 4
00 0 0 0
23 4
00 0 0 033. . . .82 6 2 4
48 1 632 . . . . .262 4
92 78 13. . . . .262 4iv
II II III IV
II I I I I I Vhh h hAf x f x h f x f x f x f x
hhfx h f x f x f x hf x
fx h f x hf x h f x h f x            

  23 4
00 0 0 031 3 281 2 1 2 9 . . . . .82 4I II II IV hfx h f x hf x hf x h f x  
   23 4 5
00 0 3 099 2 7 3 3322 8 1 6I II III IVhf x h f x h f x h f x h f x   
…… (II)

Let 1Ebe the erro r in the interval 03,xx

 20 1 2
222
10.5 0.1667 2 0.33332
0.6666hIy y y
I 

5
50380IV hhf x


munotes.in

Page 17

17




Similarly
 6
3533380xIV
xfx d x A h f x
.
.
.
 
3533380n
nxIV
nnxfx d x A h f x  Total error =  5
03 33.....80IV IV
nhfx f x f x 
Choose 10 0 3 Max ax , .....IV IV IV
n fM f x f x f x 
5180IV hEn f
5180IV hnf
5180IV bahf
Which is required error in Simpson’s 3/8th rule.

  0
001 3 5 1 2 4 2 4. . . . . 2 . . . . .3nhv
nn nvhfv d t y y y y y y y y y
          

where 00nhbaha v b vn
In Simp son’s (1/3)rd rule the given interval must be divide into even number of
equal sub intervals.

Ex. 1 :

The velocity V km/min of a train which strats from rest is given at a fixed
interval of time (min) as follows :
:0 2 4 6 8 1 01 21 41 61 82 0:01 01 82 52 92 22 01 10 50 2 0t
v

Estimate approximately the distance covered in 20 minutes by train.

Soln :

Let S (km) be the distance covered in t (min) by train
dsvdt munotes.in

Page 18

18

ds vdt
20
0Sv d t
by Simpson’s (1/3)rd rule /hba n
20 0210
2h.
  20
01 3 2 404. . . . . 2 . . . . .3nhvd t y y y y y y       
  200 4 1 02 53 21 12 2 1 82 42 053        
 204 8 02 7 23   
309.33km.
Hence in 20 mins train covered 309 .33 kms distance.

 Simpson’s (3/8)th rule :

This is 3rd child of newton cote’s 3n we get

  0
001 2 4 5 3 633. . . . . 2 . . . . .8nhv
nvhfv d t y y y y y y y y        
where 00nhbaha v b vn
while applying the Simpson’s (3/8)th rule the number of sub interval shoul d be
taken as multiple of 3.

Ex. a) :

The water under portion of a water tank is divided by horizontal plane’s
one meters apart into the following areas : 472, 392, 302, 198, 116, 60, 34, 12, 4
sq.m. use the Simpson’s (3/8)th rule to find the volume in cubic meter between
the two extreme ends.

Soln :

As given condition we can write :
Distance (S) : 1 2 3 4 5 6 7 8 9Area (A) : 472 398 302 198 116 60 34 12 04

Here sub interval are multiple of 3
9118bahn   munotes.in

Page 19

19

  9
01 2 4 5 3 01332 . . . . .8nhAds y y y y y y y y       
 31472 4 3 398 302 116 60 12 2 198 348       
 3476 2664 4648 
1351.5 cube meters.

Using these three numerical integration techniques. We can evaluate integration
very easily .

Example :

Evaluate the following integral by using Simpson’s 3/8th rule.

i)
6
201
1dxx with n = 6 Ans. : 0.912
ii) 
042 s i nxd x with n = 6 Ans. : 16.5679
iii) 14
0.2sin logxxexe d x with 12 Ans. : 3.2561


1.10 Review :

1) Newton’s cote’s formula
2) Trapezoidal rule & error
3) Simpson’s (1/3)rd rule & error
4) Simpson’s (3/8)th rule & error

1.11 Unit End Exercise :

1) Obtained trapezoidal rule using newton cote’s formula. Obtained error in
composite Trapezoidal rule.
2) Compute the trapezoidal approximation for 2
0xd x with n = 6. Compare
the estimate with the exact value.
Ans. : T = 1.81948, exact = 421.88563, error = 0.035076.
3) Use Simpson’s (1/3)rd rule to approximate 2
0xd x with n = 6. Compare
the estimate with the exact value.
Ans. : T = 1.8569, E.V. = 1.8856, Error = 0.01521.
4) Compare 12
01xx d x with n = 10 using trapezoidal rule. Compute
the error exactly and the error expression.
munotes.in

Page 20

20

5) Obtain Simpson’s (1/3)rd rule by from newton cote’s formula.

6) Obtain Simpson’s (3/8)th rule by from newton cote’s formula .


7) Obtain error in composite of Simpson’s (1/3)rd rule.

8) Obtain error in composite of Simpson’s (3/8)th rule.

9) Compute 33
03xx d x with n = 6 by Simpson’s (3/8)th rule.

10) Evaluate 6
01dxx with n = 6 using Simpson’s (1/3)rd rule.
Ans. : 3.2961.































munotes.in

Page 21

21


2

Numerical Integration – II


Chapter Structure
2.1 Objective
2.2 Introduction
2.3 Romberg’s Method
2.4 Gaussian Quadrature
2.4 Gaussian Quadrature
2.5 Gaussian Quadrature for 3 points formula
2.6 Numerical evaluation of double integrals
2.7 Simpson’s rule for double integrals
2.8 Review
2.9 Unit & Exercise

2.1 Objective :

After going th rough this chapter you will able to :
 Solve integral equation which hose analytical solution.
 Solve multiple integral by using numerica l method .
 Solve integral by different numerical integral methods.

2.2 Introduction :

Integral equation are of special application in applied mathematics. There
are such integral who have no analytic solution. To solve that we can use
numerical method of integration. Multiple integral equation are difficult to solve
using numerical method. It is easy to solve using this methods we can comparing
with the numerical expect value to find the error.

2.3 Romberg’s Method :

Richardson extrap ulation is used as the formation of a numerical
integration called Romberg integration.

Consider the integral
b
aIf x d x
munotes.in

Page 22

22

Let I 1, I2 be approximated values of I obtained by using trapezoidal rule
with different sub intervals of width h 1 & h 2 respectiv ely.

Let E 1 & E 2 be the corresponding errors.  The errors in trapezoidal rule is of order h2.


211 1 1II E I k h ……… (I)
222 2 2II E I k h ……… (I I)
from equations I & II
2211 2 2Ik hI k h 
2212 2 1IIk hh 
122221IIkhh


Substituting value of k in equation (I) we get
2121122
21IIII hhh

22 2 2221 1 1 2 1 21hh I II h hh  
2221 1 222
21hI h II
hh

Put 1hh & 2/2hh we get
2221
2
244hIh IIhh


2143III
This is know n as Romberg’s formula. Working of Romberg’s integration
method difference


123 12 234 12321
21 2 1 2 3 1 2 3 4
32 2 3 42 3
32 3 2 3 4
43
43 444 4/2 /233 3
44/4 /433
4/8 /83hI h III I IIIhI h I I I I
II I IhI h I I I
IIhI h I I
   
 


where hba.

munotes.in

Page 23

23

Use Romberg’s fo rmula to evaluate the integral 1
01dxx correct upto 4
decimal place.

Soln :

Given : integral is 1
01dxx
Here 0, 1 1 0 1ab h b a     
11yf xx
Taking 1h 0110 . 5x
y
 by trapezoidal rule.
 10 1110 . 5 0 . 7 522hIy y   
10.75I

Taking 1/2h 01 / 2 112 / 30 . 5x
y

By trapezoidal rule
 20 2 1 22hIy y y 
110 . 5 2 2 / 34 
20.7083I
Taking 14h

3 110142 410 . 80 . 6 6 6 70 . 5 7 1 40 . 5x
y


by trapezoidal rule
  30 1 2 3 22nhIy y y y y   
 110 . 5 2 0 . 80 . 6 6 6 70 . 5 7 1 48  
30.6970I munotes.in

Page 24

24
Taking 18h
00 . 1 2 50 . 2 50 . 3 7 5 0 . 5 0 . 6 2 5 0 . 7 5 0 . 8 7 5 110 . 8 8 8 9 0 . 8 0 . 7 2 7 30 . 6 6 6 70 . 6 1 5 40 . 5 7 1 40 . 5 3 3 30.5x
y

  40 8 1 2 3 4 5 6 7 22hIy y y y y y y y y       

 110 . 5 2 0 . 8 8 8 90 . 80 . 7 2 7 30 . 6 6 6 70 . 6 1 5 40 . 5 7 1 40 .533316    
40.6914I
 by Romberg’s integration method .

23 12
12 2 1 123
32 3 42 3
23 234
43
3440 . 6 9 810 . 7 5 4 / 3 0 . 6 9 4 4
3
4410.7083 0.6932 0.69302 33
41 0.6970 0.69314 3
1 0.69418hI
IIII I I
II I III
III  
  
234 123
123440.69313II
I

Ex. 2 :

Use Romberg’s formula to evaluate the integral 2
202dxx correct upto
4 decimal places.
Soln :

Given integral 2
202dx
x
Here 02 2 0 2abh b a    
21/ 2fx x
Taking 2h 020.5 0.1667x
y
10 12hIy y
Taking 1h
01 2
0.5 0.333 0.1667x
y
munotes.in

Page 25

25

 20 1 2 22hIy y y 
 10.5 0.1667 2 0.33332
20.6666I
 20.5 0.1667 0.66672 
10.6667I

Taking 12h 00 . 5 1 1 . 5 20.5 0.4444 0.3333 0.2353 0.1667x
y

  30 4 1 2 3 22hIy y y y y   
  10.5 0.1667 2 0.4444 0.3333 0.23534 
30.6732I

Taking 14h

00 . 2 5 0 . 5 0 . 7 5 1 1 . 2 5 1 . 5 1 . 7 5 2
0.5 0.4848 0.4444 0.3902 0.3333 0.2807 0.2353 0.1975 0.16667x
y

  40 8 1 2 3 4 5 6 7 22hIy y y y y y y y y       
 1[0 . 5 0 . 1 6 6 7 2( 0 . 4 8 4 8 0 . 4 4 4 4 0 . 3 9 0 280.3333 0.2807 0.2353 0.1975  
40.6749I

Taking 18h

00 . 1 2 50 . 2 5 0 . 3 7 5 0 . 5 0 . 6 2 5 0 . 7 5 0 . 8 7 5 1 1 . 1 2 5 1 .25 1.3750.5 0.4961 0.4848 0.4671 0.4444 0.4183 0.3902 0.3616 0.3333 0.3062 0.2807 0.2570x
y


1.5 1.625 1.875 20.2353 0.1975 0.1813 0.1667x
y munotes.in

Page 26

26


50 1 6 1 2 3 4 5 6 7
89 1 0 1 1 1 2 1 3 1 4 1 5[2 (2)]hIy y y y y y y y yyy y y y y y y            
 10.5 0.1667 2 0.4961 0.4848 0.4671 0.4444 0.418316
0.3902 0.3616 0.3333 0.3062 0.2807 0.2570
0.2353 0.2154 0.1975 0.1813)     


50.6753I
By Romberg’s integration method 21123hI II

12
123 123
34 234 1234
45 345 234520 . 6 6 6 7
10 . 6 6 6 6 0 . 6 6 6 5
1 0.6732 0.6754 0.67832
10.6749 0.6755 0.6755 0.67454
10.6753 0.6754 0.6754 0.67548hI
I
II
II I
II I




123450.6757I 

Ex. 1 :
Evaluate the fo llowing integral by Romberg’s method :
i) 2 1
0xed x Ans. : 0.7468

2.4 Gaussian Quadrature :

In general Gaussian quadrature is of the form
 
0nbkkakwx f x d x af x 
where w(x) > 0 In [a, b] is called the weight function and ak, xk are called
weights and nodes respectively.

Gaussian quadrature function for two points also known as Gauss
Legendre integration method.
Consider the integral
1
1If x d x
Let 2
1kkkIa f x munotes.in

Page 27

27

11 2 2af x a f x ……….. (I)

Here 1x & 2x are nodes & 1a, 2aare weights.
To find 4 unknowns 1a, 2a, 1x, 2x. Here we required h
conditions. The condition t hat equation I is valid where fx is a polynomial of
degree 3.
In this case equation I is true for 32,, , 1fx x x x .
where 3fx x
133 3 3 3
11 2 2 11 2 210xd x a x a x a x a x
  ….. (II)
where 2fx x
122 2 2 2
11 2 2 11 2 2132xd x a x a x a x ax
  ….. (III)
where fx x
12
11 2 2 11 2 210xd x a x a x a x ax
  ….. (IV)
where 1fx
1
12 12112dx a a a a
 ….. (V)
Multiply equation (IV) by 21x& subtract from (II) we get,
22
22 2 10ax x x
22 2 1 2 10ax x x x x  
20a or 20x or 21xx or 12xx.
The cases 220, 0ax & 21xx. Give result to invalid equations.
 we choose 12xx
 from equation IV
11 110ax ax ….. (VI)
From equation V & VI we get
11a & 21a
From equation III we get
2
1223x
113x
We take 113x & 22 11/3 ( )xx x 
Substitute the value of 12aa & 12xx in equation (I)
1
11133fx d x f f
       
This is known as Gaussian quadrature two points functions.
Note :
Above answer gives solution of 1
1fx d x munotes.in

Page 28

28

we need b
afx d x

There are, let x mt c
If 1xa t 
If 1xb t

Such that am c 
bmc

Solving this equation we get,
,22ba bamc
22ba baxt 
2badx dt

Substituting value of x & dx we get
1
1 22 2b
aba ba bafx d x f t d t       

Ex. 1 :

Evaluate 53
223xx d x by two points
Gaussian quadrature formula

Soln :
Given integral 53
223xx d x
Gaussian two points quadrature formula is given by
1
11133fxd x f f
       

Here 223 2 , 5fx x x a b  

Put 2ba t bax 
372tx
32dtdx
Where 21xt  
51xt   munotes.in

Page 29

29


2
512
2137 3 7 323 2 322 2ttxx d x d t              
21
194 2 4 9 9 2 1 322 2tt tdt 
12
1393 3 2 84ttd t 
1
1gt d t
Where  2 399 3 3 2 84tt t 
 by Gaussian two points quadrature formula
311334gg       
2
231193 3 2 833 4
1193 3 2 833             
331 314
46.5
52
223 4 6 . 5xx d x 

2.5 Gaussian Quadrature for 3 points formula :

Consider the integral
1
1If x d x
Let 3
1kkkIa f x
i.e. 11 2 2 33Ia f x a f x a f x ……….. (I)

Where the nodes 123xxx& weights 123aaa are to be determine.
To determine we use some as we had done in Gaussian two point’s
quadrature formula & we get,
 1
113 358 0 595 5fx d x f f f
          
[Few are left for exercise]


munotes.in

Page 30

30

Ex. 2 :

Evaluate 1
01dxx by Gaussian 3 points quadrature formula.

Soln :

Given integral 1
01dxx
1,0 , 11fx a bx 
Put 2ba bax
12tx
2dtdx
Where 01xt 
11xt
11 1
01 12
11312dtdx dtt xt 
1
1gt d t
13gtt

By Gaussian 3 points quadrature formula
 1
01 3 3 58 0 5559fx d x g g g 
11 1 158590 3 3353 53      
15 5 8 5 5
9333 5 33 5 
17 5 894 2 3
1
00.4947fx d x .



munotes.in

Page 31

31


Ex. 2 :

Evaluate the following example by two point Gaussian quadrature
formula or three points.
i) 12
0tte dt Ans. : 2 points = 3477.5439
3 points = 4967.1067

2.6 Numerical evaluation of double integrals :

If ,fxy is a continuous function & defined on a closed rectangle.
,/ ,Rx yaxb c yd  then
  ,,db
Rc afx yd x d y fx y d x d y  
Hence we are going to double integral using trapezoidal rule &
Simpson’ s 1/3rd rule.

Trapezoidal rule of double integral :
Consider ,iz iz
jiyx
yxIf x y d x d y
 where 1iihx x & 1jjky y.


















22,ji
jiyx
yxIf x y d x d y
Apply trapezoidal rules inner integration.
2
12 ,2 , ,2j
jy
ii iyhIf x y f x y f x y d y
  y
xixi + 1xi + 2
0yjyj + 1yj + 2R
munotes.in

Page 32

32


 again applying trapezoidal rule we get,


1211 1 1 222 1 2 2[, 2, ,4
2, 2 , ,,2 , , ]ii i j i j
ij ij ij
ij ij ijhkIf x y f x y f x y
fx y fx y fx y
fx y fx y fx y
   
    
 
 
Continuously working we can solve triple integral also
but it is more complex to solve.

Ex. 1 :
Evaluate 77
22
331dx dyxy
using trapezoidal rule with 2h & 2k.

Soln :
Given integral 77
22
331dx dyxy
Here 2: 3 , 5 , 7hx
2: 3 , 5 , 7ky

1231
2
335731 83 45 853 45 07 475 87 49 8xyxxxy
y
y

 by trapezoidal rule form multiple integral
 
  ,1 ,21, 1, 1 1, 22, 2, 1 2, 224
22
2ij i j i j
ij i j i j
ij ij ijhkIf ff
ff f
fff
   
   
  

  
  223, 3 2 3, 5 3, 7 2 8, 3 2 3, 5 5, 74
7,3 2 7,5 7,7ff f ff f
ff f     
11 8 2 3 4 5 8 2 3 4 2 5 0 7 4 5 8 2 7 4 9 8    
864.
munotes.in

Page 33

33

2.7 Simpson’s rule for double integrals :

It is same as trapezoidal we had done we need to apply Simpson’s 1/3rd
rule twice we get,



  11 1 2 2 2 1 2 212 1 ,
,, y , , ,4, , 4 , y9
44ij
ij ij i j i j i jij ij i j xy
xy x x y x y x yhkIff x y f x y f x
ff f ff
        

Ex. :
Evaluate 12
00xy dx dy by Simpson’s rule taking 0.5, 0.5hk.

Soln :
Given integral 12
00Ix y d x d y
Here 0.5 0, 0.5, 1hx
10 , 0 . 5 , 1 , 1 . 5 , 2ky
12 31
2
3
4
500 . 5 100 0 00.5 0 0.25 510 0 . 511.5 0 0.75 1.520 1 2xxx xy
y
y
y
y
y

By Simpson’s 1/3rd rule for double integration integral.
    
  
  11 15 12 14 13
21 25 22 24 23
31 35 32 3 4 33,, 4 ,, 2 ,9
4, , 4, , 2 ,
,, 4 , , 2 ,hkIf x y f x y f x y f x y f x yfxy fxy fxy fxy fxyfxy fxy fxy fx y fxy       

    
0.5 0.5004 0 0 0 4 014 0 . 2 50 . 7 52 0 . 59
02 4 0 . 5 1 . 5 2 1 1         
  
 0.2502 41 29  munotes.in

Page 34

34

360.259
1.

2.8 Review :

 Romberg’s method of integration.
 Gaussian quadrature formula for 2 & 3 points.
 Multiple integral by trapezoidal rule & Simpson’s 1/3rd rule.

2.9 Unit & Exercise :

1) Evaluate
0sinxd x 2h by Romberg’s integration method
Ans. : 1.99857
2) Evaluate 5
1dxx by Romberg’s integration method
Ans. : 1.6289
3) /2
2
01c o sxx x d x
 by Romb erg’s integration method
Ans. : 2.032
4) For an integral b
afxd x, derive the two point
Gaussian quadrature formula.
5) For an integral b
afxd x, derive the Romberg’s integration method
using trapezoidal rule.
6) For an integral b
afxd x, derive the three point Gaussian
quadrature formula.
7) Evaluate 22
00xy d x d y by trapezoidal rule by taken 1hk.
8) Evaluate 23
2
001dx dyxy by Simpson’s 1/3rd rule by taking 1hk.
9) Evaluate 1
023dxx by two points Gaussian quadrature.
10) Evaluate 3
2cos 2
1 sinxdx dyx Gaussian 2 & 3 points formula.
11) Evaluate 5
41 sindxxby Gaussian 2 & 3 points quadrature formula.
12) Evaluate 0.6 0.4
00sin cosxy d x d yby Trapezoidal & Simpson’s 1/3rd rule, by
taking 2hk. munotes.in

Page 35

35

13) Evaluate 55
22
11.1ddx dyxyusing T rapezoidal rule, by taking 2hk.
Ans. : 4.1345
14) Evaluate /2 /2
00sinxy d x d y
by Simpson’s 1/3rd rule, by taking /4hk. Ans. : 2.0091
15) Evaluate 1.5 1
00xyed x d yby Trapezoidal & Simpson’s 1/3rd rule, by taking 0.5hk
Ans. : Trapezoidal = 6.2334
Simpson’s = 5.0171.



































munotes.in

Page 36

36





3

Approximation of function

Chapter Structure
3.1 Objective
3.2 Introduction
3.3 General least squares method
3.4 Fitting of a straight line
3.5 Fitting of a curve of second degree polynomial to fit the second
degree polynomial OR (parabola)
3.6 Fitting of a polynomial of degree M
3.7 Fitting of curve for exponential function
3.8 Review
3.9 Unit End Exercise

3.1 Objective :

After going through this chapter you will be able to :
 Fit a straight line by least square.
 Fit a second degree polynomial by least square method.
 Fitting a polynomial of degree M by least square method.
 Fit a curve for exponential function.

3.2 Introduction :

The function approximations needs for many branches of applie d
mathematics and computer science . In mathematics, least squares can be applied
to approximating a given functions. The best approximation can be defined as
that which minimizes the difference between the original function and the
approximation we also go ing to used weighted approximation.

The process of finding the equations of curve of best fit which may be
most useful for predicting the unknown n values is known as curve fitting. Here
we get two value are by observation and other by the predicating va lue. The
difference between this values is called residual an error.

Thus, the principle of least square’s status that “The sum of the residuals
squares of is minimum”.
munotes.in

Page 37

37






3.3 General least squares method :

Let ,iixy where 1, 2, .....,in be given n -points.
To fit the curve
012,, ,, . . . . . ,nyp x a a a a ……. (I)

To the given n -points where 012,, , . . . . . ,naaa aare unknowns parameters
whose values are to be determined. Taking ixx the value of y obtained from
equation (I) we can write as 01,,, . . . . . ,ii kyP x a a a

The quantity iy is called the predicated or expected values of y at point 0xx and iy is called the observed values of y.

The difference iiyy is called the residual or error corresponding to
point ixx.

Let 01,, . . . . . ,kEE a a a be the sum of square of the residual then
21n
ii
iEy Y
 201
1;,, . . . . . . ,n
ii k
iyP x a a a


by the principle of least square the value of E will be minimum i.e.


010, 0, ....., 0kEE E
aa a   

These equation are calle d normal equations.
To solve these normal equations with help of given points and we obtain
the values of ia where 1, 2, .....,ik.
Suppose the values of ia are ** *00 1 1 ,, . . . . . ,kkaa aa aa 
substitute all these values in equation I .
** *01,,, . . . . . ,ii kyP x a a a

Which is required best fitting curve. munotes.in

Page 38

38




3.4 Fitting of a straight line :

Given the general form of a straight line ya xb
Where a, b are unknown parameters w hose values are to be determine .
Let 11 1 2 22,, ,, . . . . . , , , . . . . . . , ,kkk nn nPxy Px y P x y P x y
are the given n -points.



Taking kxx the observed value of y is ky and corresponding value on
then fitting of a straight line is kkya xb.
If kE is the residual of approximation at kxx then kk kEy Y
Let ,EE a b be the sum of squares of the residuals then
21n
ii
iEy Y
21n
ii
iya x b  ………. (I I)

by principle of least square the values E will be minimum if

2
10, 0, .....,............... 0n
ii
iEE
ab
ya x ba 
 yn
xkxhy0
P1P3
P2P5Pnnn
(, )xyyk k (a ) x+b
munotes.in

Page 39

39



120n
ii i
iya x b x
 
2
1
10n
ii i
ixy a x b x
 

2
1
11nnii iiixy a x b x

2
11 1nn nii i iii ixy a x b x    ………. (III)

Similarly
2
10n
ii
iya x bb 

121 0n
ii
iya x b
 

10n
ii
iyb a x
 

11nn
ii
iiya x b

11nn
ii
iiya xn b ………. (IV)

From the equation (III) & (IV) are called normal equation of a straight
line.

Solving these normal equations with the help of given points and we
obtained the values of *aa & *bb
**ya xb 

Which is required best fit of a straight line.


Ex :

Use least squares method to fit a straight line for the following data :
32 1 . 5 2 0 . 5210 . 5 12x
y

Also estimate the total error. munotes.in

Page 40

40


Soln :

Equation of straight line
ya xb ……. (I)

Normal equation of a straight line are
ya xb n ..…. (II)

2xy a x b x   ..…. (III)

2
232 9 6
21 4 21.5 0.5 2.25 0.7521 4 2
0.5 2 0.25 153 . 5 1 9 . 5 2 . 2 5xy x x y
xy x x y


     


We get two equation
3.5 5 5ab …… (IV)

2.25 19.5 5ab .…. (V)

Solving equation (IV) and (V) we get
0.0862 0.7862ab

Put value of ‘a’ and ‘b’ in equation (I)
0.0862 0.7862yx 

This is equation straight lin e
2
obseved expected
value value32 0 . 5 2 7 6 1 . 4 7 2 4 2 . 1 6 7 9210 . 6 1 3 6 1 . 6 1 3 8 2 . 6 0 4 31.5 0.5 0.6569 1.1569 1.338421 0 . 9 5 8 6 0 . 0 4 1 4 0 . 0 0 0 7 10.5 2 0.7431 1.2569 1.57987.6921ii i i i xy Y E y Y E 
  
    



N = 5 munotes.in

Page 41

41


 The total error
52
17.6921ii
iEy Y
 

Case i) :

Fitting a straight line for odd numbers i.e. (n = odd number)

The procedure of fitting a straight line by the method of least square for
odd numbers may be simplified as follows :
xMxh

Where M is the middle term h is difference between any two successive
value.
Here 0x normal equation can be written as
yyn b bn  
22xyxy x ax  

Case ii) :

Fittin g a straight line for even numbers of term if n is even them

2xMxh

Where M is the arithmetic mean of two middle values.
h is difference between any two successive value.

Here 0x normal equation given by
yyn b bn  
22xyxy a x ax  

Ex :
Fit a straight line ya xb to the following data : .
2010 2011 2012 2013 2014 2015 201612 16 21 24 28 32 38x
y

munotes.in

Page 42

42



Soln :

The value of x is to high  we us e change of origin.

Here n = 7 which odd number
middle value M = 2013, h = 1.

2
22010 12 3 9 362011 16 2 4 322012 21 1 1 212013 24 0 0 0
2014 28 1 1 28
2015 32 2 4 642016 38 3 9 114171 0 28 117xMxy x x yh
yx y x y 



     

Normal equation are

21174.178628xyax 

17124.42867ybh

 The remained best

 The fit of straight line is

ya xb

4.1786 24.4286yx

20134.1786 24.42861xy

4.1786 8411.5218 24.4286yx 

4.1786 8387.0932yx

munotes.in

Page 43

43

Ex. :

Fit a straight line by method of least squares also find an
estimated value for the year 2020.

Year 2011 2012 2013 2014 2015 2016Production30 35 42 48 53 60termunits

Soln :

Magnitude of years are high
 we use charge of origin
Here 6nwhich even
M = Arithmetic mean o f middle year

2013 20142013.52
1yearh
22 0 1 3 . 5 / 1xx 
2
2Year Production 2 2013.52011 30 5 25 1502012 35 3 9 1052013 42 1 1 422014 48 1 1 48
2015 53 3 9 159
2016 60 5 25 300268 0 90 210xx x x y
xy
yx x x y



     
Normal equation for even numbers of terms
2210370yax 
26844.676ybh 
 The best fit of a straight line is
ya xb
34 4 . 6 7yx
32 2 0 1 3 . 5 4 4 . 6 7yx 
61 2 0 8 1 4 4 . 6 7yx  munotes.in

Page 44

44

61 2 0 3 6 . 3 3yx
Where 2020x then
62 0 2 0 1 2 0 3 6 . 3 3y
83.67y tones .

3.5 Fitting of a curve of second degree polynomial to fit the
second degree polynomial OR (parabola) :

2Px a x b x c  …. (I)

Where a, b, c are u nknown parameters where values are to be determine.

Taking kxx the observed value of y is ky and corresponding value on
the fitting of a second degree polynomial is

2
kk kya x b xc 

If kE is the error of approximation at the points kxx then, kk kEy Y

This error may be negative or +ve zero.

Let ,,EE a b c be the sum of squares of re siduals then
21n
ii
iEy Y
 22
1n
ii i
iya x b x c  

By the principle of least squares. The value of E will be minimum.

i.e. 0, 0EEab & 0Ec
Now,
0Ea
 22
120n
ii i i
iya x b x c x
  
24 3 21,nii i i iixy a x b x c x   ………. (II)



munotes.in

Page 45

45


0Eb
 2
120n
ii i i
iya x b x c x
  
32
11 1 1nn n nii i i iii i ixy a x b x c x     ……. (II I)

0Ec

 2
121 0n
ii i
iya x b x c
  
2
11nn
ii i
iiya x b xn c   …… . (IV)

The equation (II) (III) (IV) are called normal equation as solve the normal
equations with the helps of given data and we obtain **,aa bb & *cc

Put these values in equation ………. (I)
*2 * *yP x a x b xc  

Which is required best fit for second degree polynomial .


Ex. :

Fit a second degree parabola using the method of least square.
31 0 2 483 3 8 1 1x
y

Soln :

Let equation of 2nd degree polynomial.
2ya x b xc  ………. (I)

and normal equation are
2ya x b xn c  ………. (II)

32xy a x b x c x    ………. (II I)

24 3 2xy ax bx cx      ………. (I V) munotes.in

Page 46

46



23 4 2
234 238 9 2 7 3 1 2 4 7 2
13 1 1 0 1 0 3 0 3
03 0 0 0 0 0 0 0 0
28 4 8 1 6 1 6 2 4
41 1 1 6 6 4 2 5 6 4 4 1 7 621 3 0 4 4 3 5 4 8 7 1 2 5xy x x x x y x y
xy x x x x y x y  
  
          
 Normal equ ation are
30 2 5 11abc ……… . (V)

44 30 2 87ab c  ……… . (VI)

354 44 30 125abc ……… . (VII )

Solving equation (V), (VI), (VII) we get 0.2269, 3.0774, 2.3303abc  

Put value of a, b, c in equation (I)
20.2269 3.0774 2.3303yx x  

Which is the required best fit of second degree polynomial.

3.6 Fitting of a polynomial of degree M :

To fit the mth degree polynomial. 1212 1 0.....mm m
mm myP x a x ax a x a xa
    

Where 12 1 0,,, . . . . . , ,mm maa a a a unknown parameters are whose val ues to
be determine. Taking kxx the observed value of y is ky and corresponding
values as the fitting of a mth degree polynomial is

1212 1 0.....mm m
km k m k m k kya x a x a x a x a 
   

If kE is the residual of approximation at the point kxx then, kk kEy Y
This error may be negative an +ve an zero.



munotes.in

Page 47

47


Let 11 2 0,,, . . . . . ,mm mEE a a a a  be the sum of squares residual then,
21n
ii
iEy Y
 21
11 0
1.....n
mm
im im i
iya xa x a x a

   

By principle of least squares the value of E will be minimum,

If
100, 0, ....., 0mmEE E
aa a   

Now  1
11 0
1 02. . . . . 1 0n
mm
im im i
iEya xa x a x aa

     

2
01 2
11.....nnmii i m iiiyn a a xa x a x    

Similarly


10E
a
2101
1.....nmii i i m i
ixy a x a x a x 
.
.
.
0mE
a we get
1201
1.....nmm m mii i i miixy a x ax a x
  

The above equation is called normal equa tion.

Solve the normal equation with the help of given data and we obtain ** * * *11 2 2 1 0 0,, , . . . . . . , ,mm m m m maa a a a a a a a a     .

Put all these values in equation …...... (I)
** 1 * *11 0.....mm
mmyP x a x a x a xa
   

Which is required best fit for mth degree polynomial.
munotes.in

Page 48

48



3.7 Fitting of curve for exponential function :

Let the exponential curve
0bxya e
taking l og on both side we get
log logya b x
We can write this as yAB x ……. (I)

Where 0log , A logyy a and Ba.

Thus equation is straight line.

 we can use fitting of curve of straight line by least square method.
Similarly we can curve for any exponential curve like
,xbya bya x .

Ex. :

Fit a curve of the type bxya e to the following data :
26 1 0 1 42.5 3.2 4.7 5.9x
y

Soln :

Given curve xya b
taking log both sides we get log log logya x b

Put log , log , log ,yY A = a x x Y Ab x   ……. (II)

This is straight line least square method.

 Normal equation a re
yn Ab x 
2xy A x b x    munotes.in

Page 49

49

2log log y22 . 5 0 . 6 9 3 1 0 . 9 1 6 3 0 . 4 8 0 4 0 . 6 3 5 163 . 2 1 . 7 9 1 8 1 . 1 6 3 2 3 . 2 1 0 5 2 . 0 8 9 210 4.7 2.3026 1.5476 5.3020 3.563514 5.9 2.6391 1.7750 6.9648 4.65447.4266 5.4021 15.9577 10.9672xy x x y x x yxy
 
 Normal equation

5.4021 4 7.4266Ab
10.9672 7.4266 15.9577Ab

0.5482 0.4322Ab

Equation become

0.5482 0.4322YX
But logaAe

0.5482 logae
0.5482ae
1.7301a

Equation (I) becomes
1.7301; 0.4322yx 
Which is the required least fit of curve.

3.8 Review :

 We have learn to fit a straight line by least squares method.
 We have fit quadratic curve fit by least squares method (second degree).
 We have fitting of a polynomial of degree M.
 We have fitting of a curve for exponential function.









munotes.in

Page 50

50



3.9 Unit End Exercise :

1) Fit a straight time yab x for the following data by using least square
method.

i) 4387911 10 14 12 18x
y

ii) Year 2009 2010 2011 2012 2013 2014 2015Production 2.5 3 4.2 4.8 5.3 6.4 7.3x

iii) 0.2 0.4 0.8 1 1.2 1.4346 9 1 3 1 0x
y

iv) 1121 . 563521263x
y


2) Fitting of a second degree curve using the method of least square :

i) 53 0 8 511 8 2 0 4x
y 

ii) 1941 1951 1961 1971 1981 1991 20011.1 1.3 1.6 2 2.7 3.4 4.1x
y

3) Fit the curve 2Xab x to the following data :

10 20 30 40 5081 0 1 52 1 3 0x
y

4) Fit a curve of the type bYa x to the following data :

10 20 30 40 503.5 6.2 9.5 15.3 20.4x
y

5) The pressure & volume of the gas are relation by the equation pv k
where & k are constant fit this relation for the following data using
principle of least square :
munotes.in

Page 51

51

0.5 10 1.5 2.0 2.5 3.01.62 1 0.75 0.62 0.52 0.46x
y

6) Using ‘Principle of least squares’ fit a straight line yab x for n -points.

Using principle of least squares fit a second degree polynomial
2ya x b xc  for n -points.

7) Using principle of least squares fit a polynomial of degree M for n -points.

8) Using principle of least squares fit a straight line for relation xya b.

9) Using principle of least squares for a straight line for relation bxya e.






























munotes.in

Page 52

52




4

Least squares approximation

Chapter Structure
4.1 Objective
4.2 Introduction
4.3 Orthogonal Function
4.4 Gram – Schmidt or thogonalizing process
4.5 Chebyshev polynomials
4.5 Chebyshev polynomials
4.6 Discrete Fourier Transforms (DFT)
4.7 Fast Fourier Transforms (FFT)
4.8 Review
4.9 Unit End Exercise


4.1 Objective :

After going through this chapter you will be able to :
 Fit polynomial for continuous function by least squares method.
 Orthogonalization approximation by Gram – Schmidt process.
 Fit a curve for cheby shev polynomial by least square method.
 Compute different point of sequence by DFT & FFT.

4.2 Introduction :

A more general least squares problems is the weighted least squares
approximation problem. We consider a weight function w (x) to be continuous an
(a, b) with positive mass.

i.e. 0b
awx

Given set of points 11 2 2,, . . . . . ,nnxy xy xy to given importance to the
points we assigning weights. If all the data points have amount of weights then
the weights as 1wx.

Here wx and y are known functions.
munotes.in

Page 53

53
 The above 1k equations is a system of line equation with 1k
unknowns.
 This system has a unique solution.
Let the solution is

** *00 1 1 ,, . . . . . ,kkaa aa aa 

Put all values in equation (I) we get
** 1* 2 * *
01 2 ..... .....kkkkPx a a x ax ax ax    

Which is required best approximation.

Ex. :

Constr uct a least squares linear approximation to the function 3yx
on the interval [0, 1] with respect to the weight function 1wx.

Soln :

Let the linear approximation be 01px a a x ………. (I)

Where 01,aa are unknown parameter which are to be determine.
Let 01,EE a a be the sum of squares of residuals, then,

2 b
aEw x y p x d x 
1 22
0101xa a x d x 

By principle of least squares the values of E will be minimum.


00E
a and
10E
a

 For
00E
a

13
01021 0xa a x d x  

113
0100aa x d x x d x 
munotes.in

Page 54

54

1124
010024xxaa      
1
0124aa

0142 1aa ………. (II)

For
10E
a
13
0102xa a x x d x 

1124
0100ax ax d x x d x

1123 5
010023 5xx xaa   

0 1123 5a a

0115 10 1aa …. (III)

Solving equation (II) & (III) we get
010.8 1.1aa 

Put 00.8a & 11.1a in I we get
0.8 1.1px x

Which is the remind least square linear or approxim ation.


4.3 Orthogonal Function :

The set of function 01 ,, . . . . . ,nfxfx fxin [a, b] is called
a set of orthogonal functions, with respect to a weight function wx, if
,0b
jia
iwx f x f xd x i f i jc if i j
Where jc is a real positive number further more,
if 1, 0, 1, .....,jcj n than the orthogonal set is called an orthonormal set.

munotes.in

Page 55

55





Least squares approximation of a function using orthogonal polynomial :

Let yf x be cont inuous function an [a, b] and it is approximated by
kth degree polynomial given by,
00 1 1 .....kkpx a p x a p x ap x  ………. (I)

Which orthogonal polynomial and 01,, . . . . . ,kaa a are unknown
parameters which are to be determine.

Let 01,, . . . . . ,kEE a a a by the sum of squares of residual then,

2 b
aEw x y P x d x 
 2
00 1 1 .....b
kkaEw x y a P x a P x a P x d x    
 By principle of least squares the value of E will be minimum if,

1
010, 0, ....., 0kE EE
aa a  
Now,

 00 0 1 1
02, . . . . . , 0b
kkaEwx P x y a P x a P x a P x d xa  
 00 0 0 .....bb
kkaaPxy d x w x a Px a Px Pxd x   

Since 0nknkPx is an orthogonal set we have
200b
awx P xd x c
& 000b
iaPxPx d x i
Applying the above orthogonal property we get
00b
aawxP x y d x c a 
 00
01.b
aaw x P x Y d xc
Similarly

 10 0 1 1
12. . . . . 0b
kkaEwx P x y a P x a P x a P x d xa     munotes.in

Page 56

56

 11
11 b
aaw x P x Y d xc
In general we have
10, 1, .....,b
nna
naw x P x y d x n kc 
Where 2b
nnacw x P x d x

Let the values of na be ** *00 1 1 ,, . . . . . ,kkaa aa aa  substitute all these
value in equation I, we get ** *
00 1 1,. . . . .kkpx a p x ap x ap x  

Which is the required least squares approximation.

4.4 Gram – Schmidt or thogonalizing process :

The classical Gram – Schmidt process can be used for this p starting with
the set of monomials purpose ,, . . . . . ,ncx x we can generate the orthonormal set
of polynomials kPx with respect to weight function wx can be generated
from the relation.
 1
*
0Pr 1, 2, .....k
k
kk r
rPx x a xk
 
where *
2kbr
kr ka
rwx xg xd xawx xg xd x & *
01gx

Ex. :

Using Gram – Schmidt orthogonal rizing process first three orthogonal
polynomial 01,Px Px and 2Px on [0, 1] with respe ct to weight function 1wx using this polynomial obtain the least squares approximation second
degree for the function yx on [0, 1].

Soln :

Use Gram – Schmidt orthogonalizing process to determine 012,,Px Px Px.
 1
**
01, 2, .....k
k
kk r r
rfx x a fx k
 
where ** 2bbk
kr r raaaw x x f x d x w x f x d x & *
01fx
here *
00 1fx P x munotes.in

Page 57

57

now **
11 0 0fx x afx
where  121112
10000
011 1 12xax d x d x x 
1012a
 *
1111122fx x x P x     
Now,
** *
22 0 0 2 1 1fx x a fx a fx 
Where 
3 21
20 203 1113! 11x xd xax dx 


 132
20 1
21 210 2
01 11
11 / 2 2 461111 1 1/41324 2xx d xxx d xa
xx d xxd x       


211a

*2
21/3 1 1 1/2fx x x   

21/3 1/2xx  

*2
22 1/ 6fx x x P x 

Let 00 1 1 22Px a P x a Px a P x ………. (I)

Where 012,,aaa are unknown parameter whose values are to be
determine.

Using the condition of orthogonality


1
0
2i
i b
iawx y P x d xawx P x d x





 321
1
00
0 1 2
011 322311b
axxd x
a
x dx
 
 munotes.in

Page 58

58



15/2 3/2
1
00
2132
021
11 / 2 5/2 2 3/2 53451/3 1/ 2 1/ 4 11 / 2
32 4xx
xx d x
xd x xxx        

145a

 12
0
2 2211 / 6
11 / 6b
axx x d xaxx d x




15/2 3/2
0
143 2 2
0/6
21 1233 3 6xx x d xxx x x x d x
  


17/2 5/2 3/2
0154 3 3 201
7/2 5/2 6 3/2
21 1 1
543 3 3 3 2 3 6xx x
xx x x xx 

221
75911111 1523963 6

1/315 42.9 23745 36 

247a

Put all values in equation (I) we get
 2 24 411 / 2 1 / 635 7px x x x   
20.1714 1.3714 0.5714px x x 



munotes.in

Page 59

59

4.5 Chebyshev polynomials :

The set of polynomials defined by 1cos cos 0nTx n x n  on
[-1, 1] are called the chebyshev polynomial of degrees.

It can be written as
1cos cosnTx n x
cosθn = 0 , 1 , . . . . .n
Where 1θ=c o sx on =c o sθx
we derive a recursive relation by noting that

taking 001nT x
taking 1 1nT x x
1 cos 1θnTx n
cosθc o sθ s i n θs i nθnn 1 cos 1θ = cos θcosθ sin θ sinθnTx n n n 
11 2co sθc o sθnnTxTx n
But we know that cosθ= ,x so, 11 2cos θc o sθnnTx Tx 1121nn nTx x T xTx n  
This above is a three term recurrence relation to generate the chebyshev
polynomials.


11.5.1 The orthogonal polynomial of the chebyshev polynomial :

1) nTx is a chebyshev polynomial of degree n, if n is even than nTx is
even. If n is odd then nTx is odd.

2) 1nTx for 1, 1x.

3) The chebyshev polynomial nTx are orthogonal with respect to the weight
function 21
1wxx
 & in the interval 1, 1.
1
2
100210mnif m nTx T xdx if m n
x
if n m

 
 
  munotes.in

Page 60

60

The first seven chebyshev polynomial are :

01Tx
1Tx x
2
221Tx x
3
343Tx x x
42
4881Tx x x
53
516 20 5Tx x x x
64 2
632 48 18 1Tx x x x

It is also possible to express powers of x in terms of chebyshev
polynomials we find
0 1Tx
1xT x
2
2012xT x T x
3
31134xT x T x
4
42 01438xT x T x T x 
5
53 1151 016xT x T x T x 
6
64 1 0161 5 1 032xT x T x T x T x  
and so on …..

Ex. :

Use chebyshev polynomial to obtain the least squares approximation of
second degree for the function 3fx x on the interval 1, 1 with respect to
the weight function 21
1wxx
.


Soln :

The least square approximation of second degree be
00 1 1 22px a T x a T x a T x ……. (I)

Where 012,,TTTand chebyshev polynomial & 012,,aaa are unknown
parameter whose values are to be determine. munotes.in

Page 61

61

Using the condition of orthogonalizing we obtain

20, 1, 2, .....b
ia
i b
iawx y T x d x
ai
wxT xd x



13
2 1
0 2
1
2 111
1
1
1xd xxa
dx
x






1231 11
100
2 11314
1Tx TxTx d x xTx Txdx
x




10004

[ By property of chebyshev polynomial]
0
00a.


3 11
1
31 1 2
11
1 211
1
22111 34 1
11xT xdxTx TxTx d xxaTxTx Txdx dxxx

 
 
13104 2342a


1322
1
2 1
2
2
11
1xT x wx d x xaTxdx
x





1
31 2
1
2
1
2
22
11 34
1
1Tx Tx Txd xxTx Tx xd x



 munotes.in

Page 62

62

10042
0  Put 01,aa and 2a in equation (I)
 01 23004Px T x T x T x
134Px T x
Which is the required least square approximation.

4.6 Discrete Fourier Transforms (DFT) :

A finite sequence of n real a number is given by 0, 1, . . . . . , 1xx x x n
The n -points discrete Fourier transfor m of sequence x is defined by
 12/01njk nU
jFxk X k x j e
n

Where 0, 1, 2, ....., 1kn  The discrete Fourier transformation of n -points sequence ‘r’ of complex
numbers. X0 , 1 , . . . . . , 1kxx x n

The n -points inverse d iscrete Fourier transform of sequence X is defined
by
 112 / n01Xn
ijk
N
jFk x k x j e
n



Ex. :

Compute the 4 points DFT of sequence3, 4, 5, 6x .

Soln :

Given real sequence
03 , 14 , 25 , 36xxxx   
 The 4 point s DFT of sequences is given by
 32/ 44
01
4ijkjFx k x k x j e

where k = 0, 1, 2, 3

0/ 2 3 / 24100 1 2 32ik ik ikFx x e x e x e x e    munotes.in

Page 63

63

/2 3 /24103 4 5 62ik ik ikFx e e e  
 000
41103 4 5 6 1 8 922Fx e e e    
409Fx
/2 3 /24113 45 62ii iFx e e e   


   411 3 4 cos sin 5 cos sin 6 cos3 / 2 sin 3 / 222 2Fx i i i        
  41() 1 3 4 5 1 62Fx i i    
 4112 2 12Fx i i  
41Fx i j 
3
4123 4 5 62iz i iFx e e e


    412 3 4 cos sin 5 cos 2 sin 2 6 cos3 sin 32Fx i i i         
 4123 4 1 5 1 6 12Fx 
421Fx
3/ 2 3 9/ 245133 4 6 32ii iFx e e   
  4133 45 1 62Fx i i     
 4132 22Fx i 
431Fx i
The 4 points DFT of sequence ‘ x’ is

Ex. :2

Compute the 4 points inverse DFT of sequences
111 1,, ,22 22xi i 




munotes.in

Page 64

64

Soln :

Given complex sequence is
111 1,, ,22 22xi i 

11 1 10, 1 , 2, 322 2 2xx i xx i      

 The 4 points inverse DFT of sequence ‘x’ is defined by
 312 / 44
01
4ji kiFx k x k xje 


10 / 2 3 / 24101 2 32ik i k kFx k x e x e x e x e    


1/ 2 3 / 2411 1 1 1
22 2 2 2ik i k kFx k i e e i e           

where k = 0 , 1, 2, 3

 10 0 0411 101 / 2 1 / 222 2Fx i O e i e    

1
411 1 1 1022 2 2 2Fx i i    

1
400Fx

1
400Fx

 1/ 2 3 / 2411 111 / 2 1 / 222 2iiFx i e e i e       

  1
411 1 1 11122 2 2 2Fx i i i i              

1
411Fx

1
411Fx
munotes.in

Page 65

65

  12 3
411 121 / 2 1 / 222 2ii iFx i e e i e         

  1
411 1 1 121 1 122 2 2 2Fx i i            

1
411 1 1 1222 2 2 2Fx i i    

1
421Fx

1
421Fx

 13 / 2 3 3 / 2411 131 / 2 1 / 222 2ii iFx i e e i e      

  1
411 1 1 131 122 2 2 2Fx i i i             

1
411 1 1 131122 2 2 2Fx    

1
431Fx

1
431Fx

 The 4 points inverse DFT of sequence ‘x’ is 0, 1, 1, 1xk.

4.7 Fast Fourier Transforms (FFT) :

The Fast Fourier Transforms (FFT) is the e fficient implementation of the
Discrete Fourier Transforms (DFT).
The FFT was discovered by curvely & Tukar in 1965. The FFT is based
on the following observation.

 Let data 1,, , 0 , 1 , . . . . . , 2kk k nkxy x kn be given and that p & q are
exponential polynomial of deg ree at must 1n which interpole part of the
data according to 22 22,1 , 0 , 1 , . . . . . , 1jj jjpx f x qx f x j n    then
the exponential polynomial of degree at most 21n which interpolates all the
given data is
munotes.in

Page 66

66

kkPx f x 0,1......2 1Kn
By assumption we have
111122inx inkpx e px e qxn  
Proof : Since ninx ixee has degree n & p and q have degree atmost 1n. It is
clear that phas e degree at most 21n.
We need to verify the interpolation
111122kkinx inx
kkkpx e px e qxn  
Where /1kkkinx ink nee e 
kkPx Px if k is even
/kqx n  if k is odd

Let k be even i.e. k = 2j then by assumption on P.
222jjjPx Px f x
121 22 2 jjjjxxnnn n   
  

 To compute for curve vectors.
012 1,,, . . . . . ,nxx x x x of discrete data by DFT ˆa of a is
compute component wise as

1
2/
01ˆ0, 1, ....., 1n
ik j N
kj
jaa e k NN

 

of course. The FFT c an be used to compute these co -efficient . There is an
analogs formula for the inverse DFT

1
2/
0ˆJ= 0 , 1 ,. . . . . ,N - 1n
ik j N
jk
kaa e


4.8 Review :

 Least square’s method for continuous functions.
 Approximation using orthogonal function using Gram – Schmidt
orthogona lizing process.
 Chebyshev polynomials & it’s properties.
 Discrete Fourier Transform & Fast Fourier Transform.
 Computing sequences of points by DFT & FFT.

munotes.in

Page 67

67

4.9 Unit End Exercise :

1) Construct a least squares linear approximation to the function 2yxon the
interval [0, 1] with respect to the weight function 1wx.

2) Construct a least squares quadratic approximation to the function xyeon
the interval [0, 1] with respect to the weight function 1wx.

3) Obtain the least squares polynomial approximation of degree 2 for the
function yx on it interval [0, 2] with respect to the weight function 1wx.

4) Using the Gram – Schmidt ortho gonalizing process compute the first three
orthogonal polynomial 01,px px & 01,px px which an orthogonal
an the interval 1, 1 with respect to the weight function 21
1wxx
.

5) Use chebyshev polynomial to obtain the least square approximation of
second degree for the function. 4Fx x on the interval 1, 1 with
respect to the weight function 21
1wxx
.

6) Use chebyshev polyno mial to obtain the least squares approximation for the
function 4332 2fx x x x  function 21
1wxx
.

7) Compute the 4 points DFT of sequences 1, 2, 3, 4x .

8) Use chebyhev polynomial to obtain the least squares approximation for the
function 32565 3fx x x x  on the interval 1, 1 with respect to
the weight function 21
1wxx
.

9) Compute the 4 points DFT of sequenc e 0, 1, 1, 1x.
munotes.in

Page 68

68
10) Compute the 4 points inverse 0 DFT of sequence 5, 1 , 1, 1xi i .

11) Determine the normal equation if the cubie polynomial 2302 3 xya a a x a x  is fitted to the data ,, 0 , 1 , 2 , . . . . . ,iixy i m .

12) Prove that 212x 02Tx Tx.

13) Prove that nTx is a polynomial in x of degree A.

14) Express the following polynomial as sums of chebyshev polynomials.
a) 231xx x 
b) 2412xx

15) Determine the least square method for continuous function in ,ab.




munotes.in

Page 69

69
munotes.in

Page 70

70
munotes.in

Page 71

71



Soln :

Given equation 23dyxydx with initial condition 01y.
201 1dydx 
2
232dy d yydx dx put 01y
32 1 1 5 

232
3222dy dy d yydx dx dx

221 5 21 1 2 

43 2 243 2 222 4dy dy d y dy d ydyydx dx dx dx dx dx  


21 1 2 21 5 41 5 5 4  

 By Talyor’s series method.


23 4
00 0 II II I II V
00 0 0 0 0.....2! 3! 4!xx xx xxyx y x x y y y y    

23 400 010 1 5 1 2 5 4 . . . . .2! 3! 4!xx xyx x      

2
34 5912 . . . . .24xyx x x x    

For 0.1y put 0.1x

 23 4 590.1 1 0.1 0.1 2 0.1 0.1 .....24y     

10 . 10 . 0 2 50 . 0 0 20 . 0 0 0 2 2 5   

0.1 1.127225y  .

munotes.in

Page 72

72



5.3.1 Improving Accuracy for Taylor Series Method :

The error in Taylor method is in the order of !0nxx. The accuracy can
be improved by dividing the entire interval into subintervals 01 12 23,, , . . . . .xx x x xx of equal length and computing ,1 , 2 , . . . . . ,iyx i n
successively using the Taylor series expansion using iyx as initial condition we
compute 1iyx it is given by,

II
2
11 1 1. . . . .1! 2! !mmii i
ii i i i i i iyx y x y xyx yx x x x x x xm         

In above expression put 1iixx h where 0, 1, ....., 1in
 we get
 2
II I
1 .....1! 2! !m
mii i i ihh hyx yx y x y x y xm   

Now denote 11ii i i iyx yx y yx y , above expression
becomes
2
II I
1 .....1! 2! !mmii i i ihh hyy y y ym    . This formula can be used recursively
to obtain iy values.


5.4 Picards Method :

Consider the differential equation ,dyfx ydx ………. (I)
with initial condition 00,xx yy.

Integrating the above equation in the interval 0,xx we get,

00,xx
xxdy f x y dx

00 ,x
xyx yx f xy d x

00 ,x
xyx yx f xy d x ………. (II)

Since y appears under the integral. Since on the eight the integration can
not formed. Thus we shall replace the variable y by constant or a function of ‘x’.
since we know that the initial value of y at 0xx.
munotes.in

Page 73

73

This value we assue as first approximation to the solution.

01
00 ,x
xyy f x y d x
To obtain second approximation 0xx put on right ha nd side of equation
(II) we get,
 
021
0 ,x
xyy f x y d x

In this we get 34,, . . . . .yy
 In general

01
0 ,x
n n
xyy f x y d x

This equation is known as Picard’s method.

Ex. :

Solve the differential equation ydyxedx with 00y by Picard method
find 0.2y . Check the error with analytically.

Soln. :

Given differential equation
ydyxedx …. (I)
with initial condition 00y i.e. 000, 0xy
by Picard’s method

0II
00 ,x
xyy f x y d x

0
00x
y
xyx e d x

20
0002xxxxe dx xdx   

 
021
0 ,x
xyy f x y d x

22
22
000xxxxyxe dx xe dx 
munotes.in

Page 74

74

Put 22xtx d x d t 

222
00 012x xx xtt x ed t e e e  

2
221xye

 
032
0 ,x
xyy f x y d x
2 2
2
0001xxx yxe dx x e dx  .

Now this integration is more difficult therefore we stop at 2y
2
221xyx y e 

 To find 0.2y taking 0.2x we get,
20.2
20.2 1 0.0202013ye   .

We shall solve the given equation analytically
ydyxedx
yed y x d x

Integrating both side we get,
yed y x d x
22yxec

Put 0, 0xy we get, 1c
212yxe 
2log 12x y 

Which is particular solution of d.E.
20.20.2 log 12y  
0.0202027
Error 0.0202027 0.0202013
61.4 10

munotes.in

Page 75

75



5.5 Euler’s Method :

Euler’s method is the one step method and has limited applicat ion because
of it’s law accuracy. Consider the differential equation ,dyfx ydx
………. (I) with initial value 00yx y.

Let, 01, 2, 3, .....ixxi h i 
Integration equation (I) by the limit 0x to 1x, we get,
11
00,xx
xxdy f x y dx
1
010 ,x
xyx yx f xy d x
1
010 ,x
xyx yx f xy d x
1
010 ,x
xyy f x y d x
Assuming th at 00,,fx y fx y in 01xx x
1
010 0 0 ,x
xyy f x y d x
10 0 0 10,yyf x xxx 
10 0 0,yyh f x y
Again integrating equation (I) between 1x and 2x, we get,
22
11,xx
xxdy f x y dx
2
121 ,x
xyx yx f xy d x
2
121 ,x
xyy f x y d x
21 1 1,yy h f x y

Similarly we can obtain 34,, . . . . .yy
In general 1,nn n nyy h f x y
This is known as Euler’s formula.
munotes.in

Page 76

76

Ex. : Using Euler’s method find an approximate value y corresponding to 2x
given that 231dyxdx with 12y taking interval 0.2h also find error in it.

Soln. :

Given differentiable equation is
231dyxdx
2,3 1fx y x  with initial 001, 2xy
By Euler’s method
1,nn n nyy h f x y ………. (I)

Taking 10 0 00, ,ny y h f x y 
Put 001, 2xy
2
00,1 , 2 3 1 1 4fx y f  
10 0 0,yyh f x y
120 . 2 42 . 8y 
Taking 1n in equation (1) we get 110 . 21 . 2x 

21 1 1,yy h f x y
2.8 0.2 1.2, 2.8f 
22.8 0.2 3 1.2 1 
3.864.

Taking 2n in eq uation (1) wi th 21.2 0.2 1.4x
32 2 2,yyh f x y
3.864 0.2 1.4, 3.864f 
35.24y.

Taking 3n in equation (1) with 31.4 0.2 1.6x
43 3 3,yy h f x y
45.24 0.2 1.6, 5.24yf
6.976.

Taking 4n in equation (1) with 51.6 0.2 1.8x
54 4 4,yyh f x y
6.976 0.2 1.8, 6.976f  .
9.12.


munotes.in

Page 77

77





 Thus the required approximate value of 21 1 . 7 2y .
We shall solve the given equation analytically
231dyxdx
231dy x dx

Integrating both side
231dy x dx
333xyx c 
3yx xc .

Taking initial value 001, 2xy we get value of c
321 1c 
22c
0c.
3yx x
3222y
82
10

Error = Exact Value – Approximate Value
10 9.12
0.88.

5.5.1 Accuracy of Euler’s Method :

Since Euler’s method use Taylor’s series iteratively, the truncation error
causes loss of accuracy. The truncation introduced by the step itself is known as
the local truncation error and the sum of the propagated error and local error is
called Global Truncation Error.
Consider Taylor’s expansion
23
II I I I I
1.....2! 3!nn n n nhhyy h y y y   

Since only first two terms are used in Euler’s formula the local truncation
error is given by,
23
II III
,1.....2! 3!tn n nhhEy y

The local truncation error of Eu ler’s method is of the order 2h. If the final
estimation required n steps, the global truncation error at the target point b will be munotes.in

Page 78

78

 2212
11. . . . .n
tg i n
iEl c h c c c h  
21tgEl n c h

Where 12 .....nccc cn

But 0 bxnh

0 1tgEl b x c h  .

5.6 Euler’s Modified Formula :

Consider the differential equation ,dyfx ydx ………. (I)
with initial condition 00yx y.
Integrating both side bet ween limit 0xto we get,
11
00,xx
xxdy f x y dx
1
010 ,x
xyy f x y d x
 By trapezoidal rule, we get,
 10 0 0 1 1 ,,2hyy f x y f x y  ………. (II)
But 11,fx y which occ urs on the right hand side of equation (II), cannot
be calculate since 1y is unknown so first we calculate 1y from Euler’s formula.

 By Eul er’s formula
10 0 0,yyh f x x
Put this value in equation (II) we set first approximation11y.
   1
10 0 0 1 0 0 0 ,, ,2hyy f x y f x y h f x y  
To obtain 2nd approximation to 1y i. e. 21y put 111yy in right hand side
of equation (II)
   2
10 0 0 1 0 0 0 ,, ,2hyy f x y f x y h f x y  
Continuous in this process until 1111kyy| the 1thkapproximation to 1y is
1
10 0 0 1 1 ,,2kk hyy f x y f x y 
In general munotes.in

Page 79

79

1
11 1 ,,2k
nn n n n nhyy f x y f x y
  
Where 1 ,knn n nyy h f x y [by Euler formula]

Ex. :

Using Euler’s modified formula. Find approximation value of y when 0.3xgiven that dyxydx and 01y with 0.1h.

Soln. :

Given differential equation is dyxydx with initial value condition
00x and 01y & 00.1h.
,fx y x y 

For 1st approximation :
00,0 1 1fx y
10 0 0,yyh f x y
10 . 1 1 1 . 1 
11,0 . 1 , 1 . 1 1 . 1 0 . 1 1 . 2fx y f   
1
10 0 0 1 1 ,,2hyy f x y f x y 
10 . 0 5 11 . 2 
1.11.

For 2nd approximation :
11,0 . 1 , 1 . 1 1 1 . 2 1fx y f 
21 1 1,yy h f x y
1.11 0.1 1.21 1.231 
22,0 . 2 , 1 . 2 3 1 0 . 2 1 . 2 3 1 1 . 4 3 1fx y f   
1
21 1 1 2 2 ,,2hyy f x y f x y 
1.11 0.05 1.21 1.431 
1
21.242y.

For 3rd approxi mation :
22,0 . 2 , 1 . 2 4 2 0 . 2 1 . 2 4 2 1 . 4 4 2fx y f   
32 2 2,yyh f x y
1.242 0.1 1.442
31.3862y.
munotes.in

Page 80

80
33,0 . 3 , 1 . 3 8 6 2fx y 
0.3 1.3862
1.6862.

1
32 2 2 3 3 ,,2hyy f x y f x y 
1.242 0.05 1.442 1.6862 
1.3984.

 The approximate value of 0.3 1.3984y .

5.7 Runge – Kutta Method :

Runge-Kutta method is also called as RK -method it is the generali zation
of the concept used in modified Euler’s method.

The R unge-Kutta meth od do on required the calculation of higher order
derivatives their designer to give greater accuracy.

 First order Runge -Kutta Method :

Consider a differential equation ,dyfx ydx …. (I)
with initial condition 00yx y  By Euler’s formula
10 0 0,yyh f x y
11 0 0 0,yx y h fx y

Taking 10xxh we get,
By Taylor’s series
2
II I
10 0 0.....2!hyx y h y y     Euler’s method agrees with Taylor’s series upto the first 2 term’s.
Hence Euler’s formul a is the first order Runge -Kutta method.

 Second Order Runge -Kutta Method :

Consider a differential equation ,dyfx ydx …. (I)
with initial condition 00yx y.
 By Euler’s modified formula

102hyy
 11 1 ,,2nn n n n nhyy f x yf x y   munotes.in

Page 81

81

i.e.     1 ,, ,2nn n n n n n nhyy f x yf x h y h f  
Let ,nnkh f x y
2 ,,nh nh n nkh f x yf x y
21 1,nhkh f x yk

Putting the values of 1k and 2k in (II) we get,
11 212nnyy k k 
This is known as Runge -Kutta formula of order.

 3rd Order Runge -Kutta Method :

Consider a differential equation ,dyfx ydx …. (I)
with initial condition 00yx y.
To determine 1y the 3rd order R -K is given by
 10 1 43146yy k kk  
Where 10 0,kh f x y
120 0 ,22khkh f x y 
30 0 2,kh f xh yk 
Where 20 0 1,kh f xh y k  .

5.6.1 4th Order Runge -Kutta Method or Runge -Kutta Method :

Consider a differentiable equation
,dyfx ydx with 00yx y.

To determine 1y the 4th order R -K formula is given by
 10 1 2 2 3412226yy kk k kk    

Where 10 0,kh f x y

1200 ,22kh kh f x y

2
300 ,22k h kh f x y

40 0 3,kh f xh yk  . munotes.in

Page 82

82





Ex. :

Using R -K method of order 2 approximate value of y where 1x
given that 23dyxydx with initial condition 11 . 2y .


Soln. :

Given differential equatio n 23dyxydx …… (I)
Where 2,3fx y x y with 001, 1.2xy and 0.1h.
2
00,1 , 1 . 2 3 1 1 . 2 4 . 4 4fx y f   10 0,0 . 1 4 . 4 4 0 . 4 4 4kh f x y    00 1h, 1.1, 1.2 0.444fx y k f   
1.1, 1.644f
231 . 1 1 . 6 4 4
6.
20 0 1h, 0.1 6 0.6kh f x y k    

 By R -K method of 2nd order
11 212nnyy k k 
 111.2 0.444 0.62y 
1.722

Ex. 2 :

Using R -K method of order 3, approximate value of y when 0.2x given
that 2 dyxydx with initial condition 11y and 00.1h.

Soln. :

Given diffe rential equation
2 dyxydx 2,fx y x y  with initial 001, 1xy

1st approximation and 00.1h munotes.in

Page 83

83

2
00,1 , 1 1 1 0fx y f   10 0,0 . 1 0 0kh f x y  
2
00,1 . 1 1 1 . 1 1 0 . 2 1hk fx y f   20 0 2h, 1.1, 1.021kh f x y k f  
21.1 1.02
0.189. 30 0 2h, 0.1 0.189 0.0189kh f x y k    
21
00,1 . 0 5 , 1 1 . 0 5 1 0 . 1 . 2 52 2kh fx y f   
1
40 0,0 . 1 0 . 1 0 2 5 0 . 0 1 0 2 522hkkh f x y   

 R-K method of 3rd order
 10 1 43146yy k kk  
 0.1110 4 0 . 0 1 0 2 5 0 . 0 1 8 96y   
1.00998.

2nd approximation
111.1, 1.00998, 0.1xy h 
2
11,1 . 1 1 . 0 0 9 9 8 0 . 2fxy  
10 0,0 . 1 0 . 2 0 . 0 2kh f x y  
11 1h, 1.2, 1.02998 0.41002fx y k f   
21 1 1h, 0.1 0.41002 0.041kh f x y k    
11 2h, 1.2, 1.05098fx y k f  
0.38902
31 1 2h, 0.1 0.38902 0.0389kh f x y k    
 1
11,1 . 5 , 1 . 0 1 9 9 8 0 . 3 0 2 5 222k hfx y f  
1
41 1,0 . 1 0 . 3 0 2 5 2 0 . 0 3 0 2 5 222k hkh f x y    
 11 4 310.26yy k h k k   
 11.00998 0.02 4 0.030252 0.03896  
1.03996
0.2 1.03996y 
munotes.in

Page 84

84

Ex.3:
Using R-K method of 4th order + 0 Find approximate value of y when
x = 0.2 given that 3 dyxydx with initial condition
y(0) = 1 $ h = 0.1


Soln: Given differential equation
3 dyxydx
3
00,0 , 1 , 0 . 1o f x y x y with x y h    
1st approximation
3
00,0 , 1 0 1 1fxy f  
10 0,0 . 1 1 0 . 1kh f x y  
 1
20 0/2 , 0 . 1 0 . 0 5 1 . 0 52kkh f x h y f   
30.1 0.05 1.05
0.105
 2
30 0/2 , 0 . 1 0 . 0 5 1 . 0 5 2 52kkh f xh y f   
30.1 0.05 1.0525
0.10526

40 0 3,0 . 1 0 . 1 1 . 0 5 2 6kh f x h y k f   
30.1 0.1 1.10526
0.110626

According to 4th order R -K method

 10 1 2 341226yy k k kk   
 110 . 1 2 0 . 1 0 5 2 0 . 1 0 5 2 6 0 . 1 1 0 6 2 66   
1.10591.

2nd approximation
110.1, 1.105191, 0.1xy h 
11 1,0 . 1 0 . 1 , 1 . 1 0 5 9 1kh f x y f
30.1 0.1 1.105191
0.1106191.

 1
21 1,0 . 1 5 , 1 . 1 6 0 52 2kh kh f x y h f  munotes.in

Page 85

85
munotes.in

Page 86

86
munotes.in

Page 87

87

6

Numerical Solution of Differential
Equation – II

Unit Structure
6.1 Objective
6.2 Introduction
6.3 Simultaneous First Order Differential Equations
Both Euler’s Method
6.4 Second Order Differential Equation
6.5 Multi-Step Methods (Prediction – Correction Method)
6.6 Adams Baeshforth Moulton Method
6.7 Accuracy of Multi-step Method
6.8 Model Differential Equation
6.9 Model Difference Problem
6.10 Stability of Euler Method
6.11 Review
6.11 Review
6.12 Unit End Exercise

6.1 Objective :

After studying this chapter you will be able to :
 Solve simultaneous first order differential equations.
 Solve higher order differential equation .
 Find solution of initial value problem of ordinary first order differential
equation by multi -step method by :
1) Milne – Simpson method
2) Adams Bashfarth maulton method
 Accuracy of multi -step method.
 Stability of numerical solution.

6.2 Introduction :

In previous chapter we have solve 1st order differential equation by
different method. Here we are going to s olve simultaneous differential equation
with some method also solve higher under different equation take munotes.in

Page 88

88

2
2,,dyfx y zdx with given initial condition to prove the efficiency of
estimations in one step method we need to use multi -step method by
predication and correction formula.

6.3 Simultaneous First Order Differential Equations :

Both Euler’s Method :

Taylor’s method, Picard’s method, Euler method and Runge -Kutta
method can be used to find the approximate solution to the system of first or der
differential equations.
Consider the first order differential equation s.
,,dxftxydt
,,dygtxydt
with initial condition 00,xx yy when 0tt or 00xt x and 00yt y.
Taking small change, assuming that ,thxk  and yA.
The 4th order R -K method given by
10 0 0,,kh f t x y
10 0 0,,hg t x yA
1 1
20 0 0 ,,222kh kh f t x y  A
1 1
20 0 0 ,,222kh hf t x y  AA

2 2
30 0 0 ,,222kh kh f t x y A
2 2
30 0 0 ,,222kh hf t x y AA

40 0 3 0 3,, kh f t h x k y   A
40 0 3 0 3,, hf t h x k y   AA
 10 1 2 341226xx k k kk    munotes.in

Page 89

89

 10 1 2 341226yy   AAA A
The extension of the R -K method to a system of n equation is quite
straight forward.
Similarly we solving 1st order differential equation by Taylor’s series
method.

6.3.1 Taylor’s method for simultaneous 1st order differential equation :

Consider the 1st order differential equation
,,dxftxydt
,,dygtx ydt
with initial condition 00,xx yy and 0tt. Let h be the small
change then by Taylor’s series expansion,
11 0xx t x th 
23
II II I I
00 0 0.....2! 3!hhxt h x t x t x t   
23
I II III
00 0 0.....2! 3!hhxh x x x   

11 0yy t y th 
23
I II III
00 0 0.....2! 3!hhyh y y y   

 Using Taylor series method evaluate 0.3x and 0.3y given that logdxytdt and cosdytxdt with i nitial condition 12x and 11y.

Soln. :

By Taylor series method
2
II I
11 00 0.....2!hxx t xh x x    ………. (I)
2
II I
11 00 0.....2!hyy t yh y y    ………. (II)

munotes.in

Page 90

90

Given differential equation
log , cosdx dyyt t xdt dt  
,, l o g , ,, c o sftxy y tg txy t x  
initial value 00 02, 1, 1xy t 


I
0log 1 1 log1 0dx
dxdt yt xttdt   
I
01x.


0cos cos 1 2dx
dydt txttdt  
I
01.4597y.

2
22I
2011.4597 1dx
dxdt yttdt t  
II
02.4597x.

2
221
20-sin sin1 1dy
dydt txttdt   
II
01.8415y.


3
33II II
0 2 320111.8415 11dx
dxdt yytt dt t    
0.8415.

3
33II
30cos sin 1 2.4592dx
dydt txttdt   
3.3006.
munotes.in

Page 91

91

4
44III343011 3.30061dx
dxdt ytt dt t   
2.3006.

4
44IV
40sin sin1 0.8415dy
dydt txttdt 
1.68297.

Put all these value in equation (I) and (II)
231
40.2 0.20.2 2 0.2 1 2.45972! 3!0.21.8415 2.30064!xx   

2.2466.


23
1
40.2 0.20.2 1 0.2 1.4597 1.8415 3.30062! 3!
0.21.682974!yy   

1.25082.

 Using R -K method of 4th order find approximate value of x & y at 0.1t the
following system 22,2dx dyxy xydt dt   with initial condition
00 01, 1, 0xy t .

Soln. :

To compute 10.1xx and 10.1yy use R -K method of 4th order.
 10 1 2 341226xx k k kk    ………. (I)
 10 1 2 341226yy   AAA A ………. (II)
Given 22,dy dxxy x ydt dt  
22000 000,, , ,, ,ft x y x y y t x y x y   munotes.in

Page 92

92

with 00 01, 1, 0xyt

2 2
10 0 0 0 0,, 0 . 1 0 . 1 1 1 0 . 2kh f t x y x y      
Im ,h

1
10 0 02 ,,22 2hk kh f t x y A
  20.1 0.05,1.118,1.176 0.1 2 1.18 1.761 0.3619g    
0.1 0.05, 1.1, 1.15f 
20.1 1.1 1.15
0.236.

0.1 0.05;1.118,1.176 0.1 2 1.118 1.1781 0.3619g    
0.1 0.05, 1.1, 1.15g 
20.1 2 1.1 1.15
0.35225.

2 2
30 0 0 ,,222kh kh f t x yA
  20.1 9 0.05;1.118,1.176 0.1 2 1.118 1.1781 0.3619   
0.1 0.05, 1.118, 1.761f 
20.1 1.118 1.1761
0.2426.

12
30 0 0 ,,22 2k h hg t x y AA

40 0 3 0 3,, kh f th xk y A
0.1 0.1, 1.2426, 1.3619f 
20.1 1.2426 1.3619
0.2906.
munotes.in

Page 93

93

40 0 3 0 3,, hg t h x k y AA
0.1 0.1, 1.2426, 1.3619g 
20.1 2 1.2426 1.3619g
0.4339.

Put all value in equation (I) and (II) we get,
 1110 . 2 2 0 . 2 3 6 2 0 . 2 4 2 6 0 . 2 9 2 66x   
1.2416.

1110 . 3 2 0 . 3 5 2 2 5 2 0 . 3 6 1 9 0 . 4 3 3 96y   
1.3604.
0.1 1.2416, y 0.1 1.3604x  .

6.4 Second Order Differential Equation :

Consider the second order differential equation
2
2,,dy d yfx ydt dx

000 0,Ixxdyyx y ydx
Substituting dyzdx, we get,
,,dyfx y zdx
with initial condition I00 0,yx yzx y.
These constitute system of simultaneous equations.

Ex. :

Use Runge -Kutta method to find 0.2y for the equation
2
2dy d yxydx dx given that 1, 0dyydx when 0x with 0.2h.


munotes.in

Page 94

94

Soln. :

Given 2nd order differential equation
2
2dy d yxydx dx
Taking ,,dyzf x y zdx we get,
2
2,,dyxz y g x y zdx 
with initial condition 0, 1, 0xy z with 0.2h.


1,, 0 . 2 0 0kh f x y z h z   
1,, , 0 . 2 0 0 1 0 . 2hg x y z h x z y     A



11 1
21,, 0 . 2 0 0 . 1 0 . 0 2222 2kM h kh f x y z h z           A
11 1
2 ,,222 2 2 20.2 0 0.1 0 0.01 1 0.202kM hh h hg x y z h x z y           AA


22
30.202,, 0 . 2 0 0 . 2 0 222 2 2k h kh f x y z      A
 
 22
3,, 0 . 2 0 . 1 0 . 1 0 1 , 0 . 9 922 20.2 0.0101 0.99 0.20002k h hg x y z  
   AA



40 0 0 3 3,9 , 0 . 2 0 0 . 2 0 0 0 2 0 . 0 4 0 0 0 4kh f xh y z h z M        A

 40 0 3 3,, 0 . 2 0 . 2 0 . 2 0 0 0 2 1 0 . 2 0 20.2 0.040004 0.9798 0.2039608hg x h y k z         AA
munotes.in

Page 95

95

Put all these value in R -K method of 4th order
 01 2 3 410.2 2 26yy k k k k   
110 2 0 . 0 2 2 0 . 0 2 0 20 . 0 4 0 0 46     
 110 . 1 2 0 4 0 46 
0.9799.

 01 2 3 410.2 2 26zz    AAA A
110 . 2 2 0 . 0 2 0 2 2 0 . 2 0 0 0 2 0 . 2 0 3 9 66      
 11.20800086
0.20133.

6.5 Multi-Step Methods (Prediction – Correction Method):

A pair of multi -step methods are used in conjunction with each other,
are for predicting the value of 1iy and the other for correcting the predicted
value of 1iy such method are called prediction – correction method.
6.5.1 Milne – Simpson Method :

Consider a differential equ ation
,dyfx ydx ………. (I)
with initial value 00yx y
Integrating equation I with the limit 0x and nx
44
00,xx
xxdy f x y dy
4
00 ,x
n
xyy f x y d x ………. (II)
We use newton forward difference interpolation formula in the form munotes.in

Page 96

96



200 0 0 0 03
001,, , ,2!
12,. . . . .nnfx y fx y n fx y fx ynn nfx yh   
  
Taking 00 0,fx y f from equation (II)

4
023
00 0 0 011 2.....2! 3!x
n
xnn nn nyy fn f f f d x     
where xn
0xx n h  00x
dx hdn nxn

4 23 2
23
40 0 0 0 0
033.....26nn n n nyyh fn f f f d n      

423 2 4
23 2 3
40 0 0 0 0011.....22 3 2 6 4nn n nyyh n f f f nn f            
neglecting 4th and higher power of 
2340 0 0 0 020 84833yyh f f f f    
Put 1E 

2340 0 0 0 012 24 1 20 1 8 13hyy f E f E f E f      
 
 2
40 0 0
320012 24 1 20 2 13
83 2 1hyy f E f E EfE E E f       
We know that Efx fx h
01 ,,Efx fx E x y f


   40 0 10 2 10 3 2 1 0 12 2 20 2 8 2 33hyy f h ff f ff f f ff          munotes.in

Page 97

97

 40 1 2 38483hyy f f f  
 40 12 34223hyy ff f 
11 140 1 2 34223phyy y y y 
In general
11 113 2 14223nn n n nphyy y y y    
This is known as Milne’s predication method.

6.5.2 Simpson’s Method :

Consider a differential equation
,dyfx ydx ………..
(I)
with initial condition 00yx y
Integrating equation (I) between the limit 0x to 2x we get,
2
020 ,x
xyy f x y d x
By Newton’s forward differenc e interpolation formula.

23
00 0 011 2,. . . . .2! 3!nn nn nfxy f n f f f     

 2
02
23 2 3
20 0 0 0 0132 . . . . .26x
xnnyy f f f n n n f d x       
0xx n h d xn d h limit change
0
202xhx
x

32 2 2
23
20 0 0 0 0
032
.....26nn n nnyyh fn f f f d n       munotes.in

Page 98

98


223 2 4
23 2 3
00 0 0 0011.....22 3 2 6 4nn n nyh n f f f nn f            
Neglecting 3rd and higher orde r difference
Put 1E 
240 0 0 0122 1 13yyh f E f E f    
 200 0 066 1 2 13hyf E f E E f    
 00 1 243hyf E f  
11 120 0 1243hyy y yy  
11 142 2 3 443hyy y y y  
In general
11 111 1 143nn n n nchyy y y y   
This is known as Simpson’s correction formula.

Ex. :

Given differential equation 21dyydx
Where 00y estimate 0.8y using the Milne Simpson predi cation
correction method taking0.2h.

Soln.:

Given differential equation
21dyydx
2,1 0 . 2fx y y h  
012300 . 2 0 . 4 0 . 6xxxx
By R -K method of 4th order
012300 . 2 0 2 7 0 . 4 2 2 8 0 . 6 8 4 1yyyy munotes.in

Page 99

99


 21
11 1,0 . 2 , 0 . 2 0 2 7 1 0 . 4 2 0 2 7yf x y f  
1.0411.

 21
22 2,0 . 4 , 0 . 4 2 2 8 1 0 . 4 2 2 8yf x y f  
1.1787.

 21
33 3,0 . 6 , 0 . 6 8 4 1 1 0 . 6 8 4 1yf x y f  
1.4679.

Milne predication formula is
11 140 1 2 34223hyy y y y 
 40 . 202 1 . 0 4 1 1 1 . 1 7 8 7 2 1 . 4 6 7 93  
1.0238.
Using Simpson’s correction formulas is given by
 11 142 2 3 443Chyy yy y   ………. (II)
1
44 4,0 . 8 , 1 . 0 2 3 8yf x y f
211 . 0 2 3 8
2.0482.

 40.21.1787 1.787 4 1.4679 2.04823Cy   
1.7853.



We can again use the correction formula II to refine the estimate.
1
44 4,0 . 8 , 1 . 7 8 5 3yf x y f
211 . 7 8 5 3
4.1873.

  40.21.787 1.1787 4 1.4679 4.18733Cy   
1.9278. munotes.in

Page 100

100

Which is not same again we use
1
44 4,0 . 8 , 1 . 9 2 7 8yf x y f
211 . 9 2 7 8
4.7164.

  40.21.787 1.1787 4 1.4679 4.71643Cy   
1.9631.

Note :

Milne formula is used to predict the value of 1iy, evaluation of 1if,
correction of 1iy, then to improved value of 1if.
It is also possible to use the correction formula repeatedly redefine the
estimate value of 1iy, before moving on to the next stage.

6.6 Adams Baeshforth Moulton Method :

Consider a differential equation
,dyfx y Idx
with initial condition 00yx y
Integrating equation ‘I’ between the limit 0x to 1x we get,
1
210 ,x
xyy f x y d x ……… (II)
Using Newton Backward difference interpolation formula in the term.

2
00 00 001,, , , . . . . .2!nnfx y fx y nfx y fx y ’  ’

1
223 2
22
10 0 0 0 032.....2! 6x
xnn n n nyy fn f f f d x   ’  ’ ’
0xx n hdx hdn
 where 0
101xhx
x munotes.in

Page 101

101

1 2
2
10 0 0 0
0.....2nnyy h fn f f d x  ’  ’

123 2 4
23 2 3
10 0 0 0 001.....22 3 2 4nn h nyyh n f f f nn f  ’   ’  ’  
Neglecting 4th and higher order and put 11E’  we get,
2111
10 0 0 0 015 3111 . . . . .21 2 8yyh f Ef E f Ef       

 11 2 1 2 300 0 015 3 311 2 1 3 321 2 8 8yh f Ef E E E EE f            

   00 0 1 0 1 2 0 1 2 3 24 12 10 2 9 3 324hyf f f f f f f f f f           
 10 0 1 2 355 59 37 924hyy f f f f    
111 110 0 1 2 355 59 37 924hyy y y y y
    
I II I 143 3 2 1 055 59 37 924phyy yy y y   
In general
11 1 111 2 355 59 37 924nn n n n nphyy y y y y     
This is known as adam’s basehforths predication formula for correction
formula
Given differential equation
,dyfx ydx
with initial condition 0yx y ………. (I)
Integration between limit 01xx we get,
1
010 ,x
xyy f x y d x ………. (II)
Using Newton Backward difference inter pollution formula munotes.in

Page 102

102


23
11 1 111 2,. . . . .2! 3!nn nn nfx y f nf f f  ’  ’ ’
where 1xnxh

1
022
10 1 1 1 111 2.....2! 3!x
xnn nn nyy fn f f f d x   ’  ’ ’
1xx n h where xn
dx hdn 01x
10x


 0
22 3 2 3
10 1 1 1 1
11132 . . . . .26yyhfn f nn f nn n f d n  ’  ’  ’
Neglecting 4th and higher order, we get and put 11E’ 


  11 1 2 3
10 1 111 111 2 2 1 3 322 2 4yyh f Ef E E E E E            



   10 1 10 1 01 1 0 12 24 12 2 2 1 3 324hyy f f f f f f f f f f         
 10 1 0 1291 9 524hyy f f f f   
11 1 110 1 0 1 291 9524hyy y y yy   
11 1 134 3 2 191 9524nChyy yy y y   
In general
11 1 111 1 291 9 524nn nn n nChyy yy y y     
This is known as Adams – Moultan correction formula.


munotes.in

Page 103

103

Ex. :

Using Adam – Baeshforth Moultan method
2dy ydx x
with initial 12y estimate 2y assuming 0.25h.

Soln. :

Given differential equations
2dy ydx x
2,yfxyx
with initial condition 12y
001; 2xy 
By using R -K method of 4th order we get,
01 2 3
01 2 31, 1.25, 1.5, 1.752, 3.13, 4.5, 6.13xx x x
yy y y  
  

1
00 022,41yf x y 
1
11 123 . 1 3,1 . 2 5 , 3 . 1 3 5 . 0 0 81.25yf x y f  
1
22 224 . 5,1 . 5 , 4 . 5 61.5yf x y f  
1
33 326 . 1 3,1 . 7 5 , 6 . 1 3 7 . 0 0 5 71.75yf x y f  
By Adams – Bushforth predication formu la
111 143 3 2 1 055 59 37 924phyy yy y y   
  40.256.13 55 1.0057 59 6 37 5.008 9 424py     
8.0113.
1
44 428 . 0 1 1 3,2 , 8 . 0 1 1 3 8 . 0 1 1 32yf x y f  
Adoms – Bashfourth correction formula
11 1 143 43 2 191 9524Chyy yy y y    munotes.in

Page 104

104

 40.256.13 9 8.0113 9 7.0057 5 6 5.00824Cy   
8.0073.
Hence, 28 . 0 0 7 3y .

6.7 Accuracy of Multi-step Method :

We know that for each differential equation there is an optimum steps
size h. if his also large, accuracy diminished and if it too small round off error
would dominate and reduce accuracy.
By computing the predicated used corrected values of1iy, we can
estimates the size and sign of error.
Let’s denote the predicated value 1npy and corrected value 1ncy.
Similarly denote the truncation error in predicated value by tpE and corrected
value by tcE.
1 tp npEy y  
1 tc ncEyy
Where y denotes the exact value o f y of 1nx then differences between
the error is
11 tp tc n ncpEE y y  
A large difference indicates that step size is too large. In such cases we
must reduce the size of h.
Both the Milne and Simpson formula are of order 44 and their error
terms are of order 5h.
The truncation error in Mile’s formula is
55128θ90tpEy h
The truncation error in Simpson’s formula’s
5521θ10tcEy h
Let’s assume that
5512θθyy
28tp
tcE
E
28tp tcEE munotes.in

Page 105

105

11 tp tc n ncpEE y y   becomes
11 28tc tc n ncpEE y y   
11 29tc n ncpEy y  
1129nncptcyy
E

If the answer is required to a precision of a decimal digit s then
110.5 1029nn cpdtcyy
E
 
1129 0.5 10 15 10ddnn cpyy   | 
Similarly for Adams – Bashfourth method we get,
 112700.5 10 7 1019ddnn cpyy   | 





6.7 Stability :

Stability of a numerical method ensures that small changes in the initial
conditions s hould not lead to large changes in the solution. This is particularly
important as the initial conditions. May not be given exactly. The approximate
solution computed with error in initial conditions is further used as the initial
condition for computing s olution at the next grid point. This accounts for large
deviation in the solution started with small initial errors also round off error’s
in computations may also affect the accuracy of the solutions at a grid point.

Euler method is found to stable :

Stability is the necessary and sufficient condition for convergence.

6.7.1 Stability of Numerical Solutions :

Consider a differential equation
,dyfx ydx
with initial condition 00yx y. munotes.in

Page 106

106

The solution with nearby initial values are close 0yx is called the
stable. If the solution with near by initial values diverse from yx then the
solution is unstable.
For e.g. 1dyydx with initial condition 01y
Exact solution of given DE is
1dydxy
log 1 logyx c
11xyc e
By initial condition 01y
001yc e
0c
1yx
Which is exact solution of given differential equation if 101 . 0 0 0 1y 
then find for the same 01 . 0 0 0 1y .
001yc e
01.0001 1ec
0.0001 0c
1.0001yx .

The exact solution of this differential equation is 1yx where 0c.
However if we use other initial values 0c and the so lution with
diverge from the 1yx.
It is difficult to obtain accurate numerical solution to an unstable initial
value problem. If a small numerical error occur solution diverge from the two
solution.



6.8 Model Differential Equation :

Consider the model first order differential equation
dyydx ………. (I) munotes.in

Page 107

107

with initial condition 00yy where  is constant and it may be oral
as complex number
dyydx
dydxy
Integrating both sides
dydxy
log logyx c
xyx c e ………. (II)
To find ‘c’ use 01 0y
00yc e
0yc.
Put 0yc in equation (II) we get,
0xyx y e
Which is exact solution of differential equation.
 If 0, a small change in the initial condition causes only small
change in the solution and therefore the problem is a stable problem.
 If 0 large changes in solutions will occurs and the problem is
unstable.

6.9 Model Difference Problem :

Consider the model difference problem
1,n 1 , 2 , 3 , . . . . .nnyk y
Where the initial value 0y is given as 6 is complex number
0nnEy k y
0n Ek y
Auxiliary equatio n is
0Ek
Ek
nCF A k and 0PF
The complete solution is
nyC F P F
nnyA k ………. (I) munotes.in

Page 108

108

To find A put 0A
00yA k
0yA
Put this value in equation I
0n
nyyV.
The solution is bounded if 61 which is the solution of difference
equation.
The connec tion between the exact solution and the difference solution is
clear if we evaluate the exact solution at the points nxn h where 1, 2, 3, .....n and 0h.
0xyy e
0nx
nyx e y
0nh
nyey
0n
nyk y
hke
If exact solution is bounded then 1heV. this is possible if 1Re 0hh.
i.e. in the 12hh plane the region of stability of the exact solution is
the left half plane as shown in fig.














Stability region of exact solution.
A single step method is called
i) Also olutely stable if 1k 0Oh
Ohmunotes.in

Page 109

109

ii) Relatively stable if hke
iii) Periodically stable if 1k and  is purely imaginary.

6.10 Stability of Euler Method :

Given differential equation is
dyydx
with initial condition 00yy
,fx y y
By Euler’s method state that
1 ,nn n nyy h f x y
nnyh y
1nyh d
1nnyk y
Where 1kh.  The solution of this differential equation is 0n
nyk y.
Since the exact problem has an exponentially decaying solution for 0, a stable numerical method should exhibit the same b ehavior.  In order to ensure stability of Euler’s method we need that the so called
growth factor 1h 0nyo as nof if 11kh  
Here we disc uss the following cases.

Case I :

If  is real and 0 then
11h V 
11 1h   
20h  
2n .
Thus Euler’s method i s only conditionally stable. i.e. the step size has to
be choosen sufficiently small to ensure stability. The set of h function the
growth factor is less than are is called the linear stability domain D.

munotes.in

Page 110

110

Case II :

 is purely imaginary
2
2211 1 1kh i h h      
Euler method’s unstable where  is ure imaginary.

Case III :

 is a complex number
Let zc
11x 
1211ih   
1211hi h  
22
1211hh
22
1211hh 
A rather small circular subset of the left half of the complex plane.


Diagram

Im(O,h)Rc(Oh)
(O)-2

Stability region of Euler method.
where 1
11
11 1h   
20h 
02h
 Euler method is stable inside the circle.
munotes.in

Page 111

111

6.11 Review

In this chapter we have learn

 Simultaneous first under differential equation solution by numerical
method.
 Solution of second order differential equations by numerical method.
 Multi-step method : (Predication – correction method)

i) Milne – Simpson Method
ii) Adam – Bashfourth Maulton Method

 Accuracy of multi -step method.
 Stability.


6.12 Unit End Exercise :

1) Using Taylor’s series method evaluate 0.2x and 0.2y given that sint dx dyye x tdt dt  with 01x and 00y.

2) Using Taylor’s series method compute 0.1x and 0.1y correct upto
h-decimal places given that dxytdtand dyxtdt with 01x and 01y.

3) Using R -K method of 4th order find the approximate value of x and y at
t = 0.1 the following system 2, 3dx dyxy x ydt dt   with 001, y 0x.

4) Using R -K method solve the differential equation 2
2d x tdxxdt dtwith 101 , 00xx taking 0.1h to find 0.2x and 10.2x .
munotes.in

Page 112

112

5) Use Taylor’s series method to find 0.2x and 10.2x given that
2 222dx d xtxdt dt with 01x and 100x with 0.2h.


6) Use Taylor series method to find the value of y at t = 0.1 correct upto
decimal place if y satisfies the equations 2
2dytydt given that 1, 1dyydt
when 2t with 0.1h.

7) Given that 222dyxydxwith 1.09, 0.1 1, 0.2 0.89yy yf   and 0.3 0.7605y  use Milne’s Simpson’s method to determine 0.4y
correct to 4 decimal places.

8) Given 2dyxydx with 00y evaluate 0.8y using Milne –
Simpson’s method obtain the starting values from Euler’s method.

9) Using Adam – Bashforth method determine 4.4y given that
1252xy y with 41y 4.1 1.0049, 4.2 1.0097, 4.3 1.0143yyy  .


10) Given 2dyyxdx with 01 ,0 . 21 . 2 1 8 6 ,0 . 41 . 4 6 8 1 ,0 . 61 . 7 3 7 8yy y y    compute 0.8y
and 1.0y by Adam – Bashforth method correction upto 4th decimal
places.

11) Determine Milne Simpson’s method to solve ordinary differential
equation with initial condition 00yx y.
munotes.in

Page 113

113

12) Determine Adams – Moultan correction formula solve 1st order
differential equation with initial condition 00yx y.

13) Determine Adams – Bashforth predication formula to solve 1st order
differential equation with initial condition 00yx y.


14) Determine accuracy of multi -step method of Milne’s Simpson’s method.


______________________
munotes.in

Page 114

114

7

Numerical Solution s of Partial
Differential Equation s

Unit Structure

7.0 Objective
7.1 Introduction
7.2 Finite Difference Approximations to Derivatives
7.3 Laplace Equation of Two Dimension
7.4 Jacobi’s Iteration Formula
7.5 One Dimensional Heat Equation (Parabolic Equation)
7.6 Crank – Nicholson Difference Method
7.7 Given One Dimensional Heat Equation 2
2UUxt under the conditions
7.8 One Dimensional Wave Equation (Hyperbolic Equation)
7.9 Review
7.10 Unit End Exercise

7.0 Objective :

After studying this chapter you will be able to :
 Solve partial differential equation by numerical method.
 Solve Laplace equation using numerical method .
 Sole heat equation of one dimension by numerical method.
 Solve wave equation of are dimension by numerical method.

7.1 Introduction :

Partial differential equations are used in a number of physical problems
such as fluid flow heat transfer, solid mechanics and biological process . There
are three types of equations. Hyp erbolic equations are most commonly
associated with advection, parabolic equations are most commonly associated
with diffusion and elliptic equation are most commonly associated with steady
states of either parabolic o r hyperbolic parabolic problems classi fication of
general linear partial differential equations.
munotes.in

Page 115

115

General partial differential equation is of the form.
 
222
22A, B, C, D,
E, F, G, 0yyy yxy xy xy xyxx y yxdyxy xyu x ydx      .… (I)
Where U is or known function of x and y and A, B, C, D, E, F
and G are also functions.
This equati ons is called
1) Elliptic : I f 240BA C
For eg. 22
22UU0xy Laplace equations
22
22UU,fx yxy Poisson’s equations

2) Parabolic : If 240BA C
For eg. 2
22UUCyx are dimensional he at conduction equation.

3) Hyperbolic : If 240BA C
For eg. 22
222UUCyx the wave equation.

7.2 Finite Difference Approximations to Derivatives :

Let the Parabolic : If ,xy plane be divided into a network of
rectangles of sizes xh and yk by drawing the sets of lines ,0 , 1 , 2 , . . . . .xi h i
0, 1, 2, .....yj k j


Grid point





0yxmunotes.in

Page 116

116

The point of intersection of these families of lines are called mesh
points, lattice points on grid points.
Let U, yx Bethe function of the variable x and y
U, y U , y OUxxh x hh 
Expand U,xh y in Taylor’s series expansion
2U, U , U U . . . . .2!xhxh y x y h x x   
2U, U ,UU . . . . .2!xxh y x y hxxh  
U, U ,UOxxh y x yhh 
This is forward difference approximation for Ux. Similarly we have the
approximation.
U, U , OUxxy x hy hh
This is backward difference approximation for Ux.
2U,U,O
U2xxh y xh y hh

This is central difference approximation for Ux.
Let U, U , U ,xy i h i j ij
U1 , U , OUij i j hxh  for forward
U, U 1 , OUij i j hxh for backward
2U1 , U1 , O
U2ij ij hxh
 for central
22U1 , 2 U , U1 ,UOij i j i jxx hh 
Similarly w e have the approximation
U, 1 U,UOij i jykh for forward
U, U, 1UOij i jykh for back ward
U, 1 U, 1UO2ij ijykk  for central munotes.in

Page 117

117

22U, 1 2 U, U, 1UOij i j ijyy kk  

7.3 Laplace Equation of Two Dimension :

Consider the Laplace equation of two diameter are given by, UU0xx yy ….. (I)
Consider a square region R for which U,xy is known at the
boundary and divide R into small squares of side h as shown in fig. where 12 1 6,, . . . . . ,bb bare boundary values.


















We know that
1,21, ,UUU 2 Uij
ij i j xxh
 …... (II)
,1 , ,12U2 U UUij ij ijyyk …... (III)
Replace (II) and (III) in equation (I) we get,
  ,1 , , , 1 1, , 1
2U2 U U U2 U U0ij i j ij ij ij i jhk   ….. (IV)



b13b12b11b10b9
b14
b15
b16
b1b8b7b6b5U7 U8 U9
b2 b3 b4U4 U5 U6
U6 U2 U3
munotes.in

Page 118

118

For value of hk i.e. fo r square grid of the mesh size h equation
can be written as
1, , 1, , 1 , , 1U2 U UU2 U U0ij i j ij i j i j i j 
,1 , 1 , , 1 , 11UU U U U4ij i j i j ij ij      …... (V)
These shows that the values of U,xy is the average of its four
neighbors to the East, West, North, South is called standard five points formula
[S, F, P, F].
This formula is also known as Liebman’s averaging procedure.

Note :

The Laplace equation remains unchanged when coordinate are rotat ed
through 45°.
A formula similar to the (VI) is sometimes used with convenience it is
given as
 1, 1 1, 1 1, 1 1, 11UU U U U4ij i j i j i j i j       …... (VI)
This is known as diagonal five point formula as these points Lies on the
diagonals (DFPE). But it is less accur ate than standard five point’s formula.
















We use the following five point formula to set the initial value of U at
the centre.
 51 5 9 1 31U4bbbb   E(i,j)(i+1, j+1)
(i+1, j-1)(i+1, j+1)
(i-1, j-1)munotes.in

Page 119

119

Then the approximate values of 1379U,U,U,U are calculated by the
diagonal five point formula.
 11 3 5 1 51U4bbbb  
 33 5 7 51U4bbbU  
 71 5 5 1 1 1 31U4bUbb  
 95 7 9 1 11U4Ubb b  
The values of the remaining interion point i.e. 24 16U, U , U and 8U are
obtai ned by the standard five point formula.
 23 3 5 11U4bU U U  
 41 5 7 1 51U4UUU b  
 63 7 9 51U4Ub UU  
 85 9 1 1 71U4UUbU  
Thus we obtain all initial values 123,,. . . . .gUUU U once their accuracy
can be improved by th e repeate d application of either Jacobi iteration formula
as Gauss – Seidel iteration formula.

7.4 Jacobi’s Iteration Formula :

Let ,Unij be the nth iterative value of ,Uij then Jacobi’s iterative
procedur e is given below.
1,1 , 1 , , 1 , 11UU U U U4n nnnn
ij i j i j ij ij
   

 Gauss – Seidel Method :

This method utilizes the latest iterative value available and scans
the mesh points symmetrically from left to right along successive rows.
The formula is given below.
11 11, 1, , 1 , 1,1UU U U U4nn n n n
ij ij i j i j ij 
   
munotes.in

Page 120

120

Ex. :

Solve Laplace equation UU0xx yy in the domain of the figure given
below by Gauss – Seidel method.











Soln. :

  1
12 11U1 0 1 0 U U4nn n   
  11
21 31U2 0 2 0 U U4nn n 
 11
31 31U1 0 1 0 U U4nn 
  11 1
41 31U2 0 2 0 U U4nn n  
We use SFPF
1, 1, , 1 , 11UU U U U4ij i j i j i j i j    
and DFPF as
1, 1 1, 1 1, 1 1, 11UUUUU4ij i j i j i j i j      
Now initially 1234U0 , U 0 , U0 , U0 ,

First Iteration :
 1
11U1 0 0 1 0 0 54   
 2
21U2 0 0 2 0 5 1 1 . 2 54   
 1
31U1 0 1 1 . 2 5 1 0 0 7 . 8 1 2 54   102010
2020102010munotes.in

Page 121

121

 1
41U2 0 7 . 8 1 2 5 2 0 5 1 3 . 2 04   


Second Iteration :
 2
11U1 0 1 0 1 1 . 2 5 1 3 . 2 0 1 1 . 1 1 2 54  
 2
21U2 0 2 0 1 1 . 1 1 2 5 7 . 8 1 2 5 1 4 . 7 34   
 2
31U1 0 1 0 1 4 . 7 3 1 3 . 2 0 1 1 . 9 84  
 2
41U2 0 2 0 1 1 . 1 1 2 5 1 1 . 9 8 1 5 . 7 74  

Third Iteration :
 3
11U1 0 1 0 1 4 . 7 3 1 5 . 7 7 1 2 . 6 34  
 3
21U2 0 2 0 1 2 . 6 3 1 1 . 9 8 1 6 . 1 54 
 3
31U1 0 1 0 1 6 . 1 5 1 5 . 7 7 1 2 . 9 84  
 3
41U2 0 2 0 1 2 . 6 3 1 2 . 9 8 1 6 . 4 04 

Fourth Iteration :
 4
11U1 0 1 0 1 6 . 1 5 1 6 . 4 0 1 3 . 1 44  
 4
21U2 0 2 0 1 3 . 1 4 1 2 . 9 8 1 6 . 5 34 
 4
31U1 0 1 0 1 6 . 5 3 1 6 . 4 0 1 3 . 2 34  
 4
41U2 0 2 0 1 3 . 1 4 1 3 . 2 3 1 6 . 5 94 






munotes.in

Page 122

122

Fifth Iteration :
 5
11U1 0 1 0 1 6 . 5 3 1 6 . 5 9 1 3 . 2 84  
 5
21U2 0 2 0 1 3 . 2 8 1 3 . 2 3 1 6 . 6 34 
 5
31U1 0 1 0 1 6 . 6 3 1 6 . 5 9 1 3 . 3 14  
 5
41U2 0 2 0 1 3 . 2 8 1 3 . 3 1 1 6 . 6 54 

Sixth Iteration :
 6
11U1 0 1 0 1 6 . 5 3 1 6 . 5 5 1 3 . 3 34  
 6
21U2 0 2 0 1 3 . 3 2 1 3 . 3 1 1 6 . 6 64 
 6
31U1 0 1 0 1 6 . 6 7 1 6 . 6 5 1 3 . 3 34  
 6
41U2 0 2 0 1 3 . 3 2 1 3 . 3 3 1 6 . 6 64 
12 34U1 3 . 3 3 , U1 6 . 6 6 , U1 3 . 3 3 , U1 6 . 6 6    .
Ex. 2 :

Solve UU0xx yy by h i.e. bman iteration process for the domain of
the figure given below :















050100500100200100U1 U2 U3
U4 U5 U6
U7 U8 U9100
200
100050100500munotes.in

Page 123

123




Soln. :

We use standard five point formula
1, 1, , 1 , 11UU U U U4ij i j i j i j i j    
and diagonal five point formula
1, 1, 1 1, 1 1, 1 1, 11UU UUUU4ij i j i j i j i j i j           
Values given on the figure are symmetrical about middle line
1237UUUU and 28UU & 46UU.
  51U2 0 0 1 0 0 2 0 0 1 0 0 1 5 04    [Standard formula]
  41U0 2 0 0 1 5 0 1 0 0 1 1 2 . 54  [Diagonal formula]
 Similarly 1357UUUU1 1 2 . 5
  21U1 0 0 1 1 2 . 5 1 5 0 1 2 2 . 54  
118.75.
82UU1 1 8 . 7 5
  41U1 1 2 . 5 2 0 0 1 1 2 . 5 1 5 04  
143.75.
46UU 1 4 3 . 7 5
12 34 5
67 89U1 1 2 . 5 , U1 1 8 . 7 5 , U1 1 2 . 5 , U1 4 3 . 7 5 , U1 5 0 ,U1 4 3 . 7 5 , U1 1 2 . 5 , U1 1 8 . 7 5 , U1 1 2 . 5 .     

Now by Gauss – Seidel Method :
  11 11, 1, , 1 , 11UU U U U4nn nn n
ij i j i j i j i j 
   
  1
12 41U1 0 0 U 5 0 U4nn n   
11 135 7=U U Unn n 
   11 1 1
21 3 51UU U 1 0 0 U4nn n n     
18=Un munotes.in

Page 124

124

  11 1
45 1 71U2 0 0 U U U4nn n n    

16=Un
 11 1 1 154 6 2 81UU U U U4nn n n n     

First Iteration we get,
 1
11U1 0 0 1 1 8 . 7 5 5 0 1 4 3 . 7 54 
=1 0 3 . 1 2 5.
 1
21U1 0 3 . 1 2 5 1 0 3 . 1 2 5 1 0 0 1 5 04  
=1 1 4 . 0 6.
11
28UU 1 1 4 . 0 6
 1
41U2 0 0 1 5 0 1 0 3 . 1 2 5 1 0 3 . 1 2 54  
=1 3 9 . 0 6.
11
46UU 1 3 9 . 0 6 .
 1
51U1 3 9 . 0 6 1 3 9 . 0 6 1 1 4 . 0 6 1 1 4 . 0 64  
=1 2 6 . 5 6.
Similarly the Liebman’s iteration are given by,
Iteration 1397UUUU 28UU 46UU 5U

2nd 100.8 106.9 132.1 119.5 3rd 97.3 103.5 128.8 116.2 4th 95.6 101.9 126.9 114.4 5th 94.7 101.0 726 135.5 6th 94.2 100.5 125.5 113 7th 94 100.3 125.3 112.8 8th 93.9 100.2 125.2 112.7 9th 93.9 100.1 125.1 112.6
1379UUUU9 3 . 9
28UU 1 0 0
46UU 1 2 5 . 1
5U1 1 2 . 6. munotes.in

Page 125

125




7.5 One Dimensional Heat Equation (Parabolic Equation) :

Consider a one dimensional heat equation
2
2yUxtD …… (I)
with initial condition U, 0xf x and boundary conditions 0U0 ,tT and 1U1 ,tT.
Divide xt – plane into small rectangles of size h and k in x and t directions respectively
Let U, U ,x t ih jk where ,0 , 1 , 2 , . . . . .ij t t
We write partial derivative in equation (I) we get,
21, 1,222ij i j ijUU Uyxh 
and ,1ij i jUUUtk 
 From equation (I)
22UUtxD 1, , 1,1
2U2 U U UUij i j ijij ij
khD    
  ,1 1 , , 1 , 2UU U2 U Uij i j i j ij i jkhD   
Put 2/khD ….. (II)
We get,
,1 , j 1 , j , 1 ,UU U 2 U Uij i i ij i j  
,1 1 , , j 1 ,UU 1 2 U Uij i j i i j     ….. (III)
It gives formula for unknown temperature ,1Uij at ,1ij when
reaming values are known. Hence the method is called explicit method. These
method is valid 1U<2 choose k in such a way that co -efficient of ,Uij in
equation (III) will becom e zero.
i.e. 112 02  
put 12 in equation (III) we get, munotes.in

Page 126

126

212khD
22hkD
 Equation (III) reduce to the form
,1 1 , 1 ,1UU U2ij i j i j  ……….
(IV)
This is called the Bendre Schmidt or Schmidt recurrence relation.
,1 1 , 1 ,1UU U2ij i j i j 
Given the values of U at the interval points when the boundary
conditions are known.

7.6 Crank – Nicholson Difference Method :

Crank – Nichol son proposed a method in which 22yx is repl aced
by the average of its finite difference approximation on the thjand 1thj
rows thus.
2
1, , 1, 1 1 , 1 1, 1
22 2U2 U U U , U 2 UU1
2ij i j ij i j i j ij y
xh h         heat equation
22UytxD can be written as ,1 , 1 , , 1 , 1 ,1 ,1 1 ,122UU U2 U U U , 2 UU 1
2ij ij i j ij i j i j ij i j
kh hD               1, , 1, 1, 1 , 1 1, 1,1 , 22 2U2 U U U , 2 UUUU2ij i j ij ij i j ij
ij ijk
hh hD     
  
Put 2khD ,1 , 1 , , 1 , 1 ,1 ,1 1 ,12U U U 2 U U U , 2 U Uij ij i j ij i j i j ij i j              ,1 1 ,1 1 ,1 , 1 , 1 , 21 U U U 21 U U Uij i j i j ij i j i j                 ….. (I)

This is known as “Crank – Nicolson” Difference formula.
This formula is convergent for all values of . munotes.in

Page 127

127

To choose the value of k in such away that the co -efficient of ,Uij in
equation (I) will become zero i.e. 1.
21khD
2khD
For 2khD equation (I) become reduce to
,1 1 , 1 , 1 ,1 1 ,11UU U U U4ij i j i j i j i j       .

Ex. :

Solve one dimensional heat equation 2
22U Uxt
When 01 . 5t and 0xh with initial condition U, 0 5 0 0xh x x h   and the boundary conditions.
U0 , 0 0 1 . 5tt 
U, 0 0 1 . 5ht t  
Using Schmidt method.

Soln. :

Given equation
2
22U Uxt
with heat equation we get,
12, 1,2hk D
2210.2522 2htk 
The bender Schmidt equation
,1 1 , 1 ,1UU U2ij i j i j 
U, 0 5 0 0 . 2 5 , 1xh xt x t    
Where 0, 1, 2, 3, 4i and 0.25, 0.50, 0.75, 1.00, 1.25, 1.5j.

,0 U, 0 5 0 U 5 0ixh x h i  
munotes.in

Page 128

128

U0 , 0 0U1 , 0 5 04 1 1 5 0  
U2 , 0 5 04 1 1 0 0 

U3 , 0 5 04 3 U4 , 0 5 04 4 0  

Thus we can generate successfully U,xt

j 0 1 2 3 4
0.00 0 150 100 50 0
0.25 0 50 100 50 0
0.50 0 50 50 50 0
0.75 0 25 25 25 0
0.10 0 12.5 25 12.5 0
1.25 0 12.5 12.5 12.5 0
1.3 0 6.25 12.5 6.25 0

1, 1.5 2, 1.5 3, 1.5U6 . 2 5 , U1 2 . 5 , U6 . 2 5   .

7.7 Given One Dimensional Heat Equation 2
2UUxt under
the conditions.

U, 0 2 0x  for 03x and U0 , 0t
U3 , 3 0t for 0t taking 1, 1hk and using Cra nk –
Nicolson method to compute ‘U’ for one formed step on.

Soln.:

Crank – Nicolson difference formula is given by

,1 1 ,1 1 ,1 , 1 , 1 , 21 U U U 21 U U Uij i j i j ij i j i j               ………. (I)
Given that 1, 1, 1hkD

2 21111k
hD 
Put 1in equation (I) we get, munotes.in

Page 129

129

,1 1 ,1 ,1 1 , 1 ,4U U U U Uij i j i j i j i j          ,1 1 , 1 , 1 ,1 1 ,11UU U U U4ij i j i j i j i j       
………. (II)
,0U, 0 2 0 U , 0 2 0 U 2 0 0 , 1 , 2 , 3i x ih i    
0, 0 1, 0 2,0 3,0U2 0 , U2 0 , U2 0 , U2 0
0,U0 , 0 U0 , 0 U 0j t jk  
0,1U0.
U3 , 3 0 U3, 3 3 0th k 
 3, 1 3, 1U3 0 U 3 0 .















J = 1
J = 00
20U1, i
20U2,1
2030
20123tyi = 0munotes.in

Page 130

130


The mesh point 11U and 21U we need to find using equation (II)
1, 1 0, 0 2, 0.5 0, 1 2, 11UU U U U4  
1, 1 2, 11U2 0 2 0 0 U4  
1, 1 2, 14U 4 0 U
1, 1 2, 14U U 4 0 ………. (III)

For 2,1U
2,1 1, 0 3, 0 1,1 3,11UU U U U4  
2,1 1,11U2 0 2 0 U 3 04  
2,1 1,1U7 0 U
2,1 1,1UU 7 0 ………. (IV)

Solving equation (III) and (IV) we get,
1, 1U1 5 . 3 3
2,1U2 1 . 3 3.

7.8 One Dimensional Wave Equation (Hyperbolic
Equation) :

Consider a one dimensional wave equation 22 222UUxt  under the
boundary condition U0 , U , U1 , 0tt and initial conditions U, 0 , U , 0txf xxg x .
Divide xt – plane with small rectangles of sides h and k in x and t
directions resp.
Let U, U ,x t ih ik where ,0 1 2ij. First we unite initial and
boundary conditions in difference notations. 0,U0 , 0 U0 , 0 U 0j t jk   U, 0 U , 0tn k j k  where nhA
,,U0ni j munotes.in

Page 131

131
U, 0 U , 0x f x ih f ih  ,0Uif ih
1 1UU
U, 02i i
txg x g i hk  
,1 , 1UU U,,2ij ijxttk
Partial derivative of second ord er is given by
2
1,21, 2U UU2 Uij
ij i jht
 
2
i, 12,1 , 2U UU2 Uj
ij ijkt

Put 22Ux and 22Ut in equation (I) we get, 1, j 1, , 1 , 12
22U2 U U U 2 U Uii j i j i j i j i j
hkD     
22
,j 1 , 1 1 ,j 1 , 2U2 U U U2 U Uii j i j ii j i jkhD
     
Put khD 22,j 1 , 1 , 1 , , 1U2 1 U UUUii j i j i j i j       ……….
(III)
This scheme is called the explicit scheme for one dimensional wave
equation.
Equation (III) is valid where 01 and invalid where 1.
In case of 1 we get equation (III)
,1 1 , 1 , , 1UUU Uij i j i j ij  

Ex. :

Solve UUxx tt with conditions
1U, 0 0 , U 1 , 0 , U, 02xxxt x  and U, 0 0x taking 0.1h and 0.1k for 00 . 2t.



munotes.in

Page 132

132

Soln. :

Given equation UUxx tt
i.e. 2222UUxt
2 10 . 11, 2 10.1khD   
Where 1 solution of value equation is
,1 1 , 1 , ,1ij i j i j ijUUUU   ………. (II)
U0 , U1 , 0 , U , 0tt o j  and U, 0oj
1U, 02xxx
1Ui , 02ii
0.1 1 0.1U0 . 1 , 0 0 . 0 4 52
0.2 1 0.2U0 . 2 , 0 0 . 0 82
Similarly we get
U0 . 3 , 0 0 . 1 0 5 , 40 . 4 , 0 0 . 1 2
U0 . 5 , 0 0 . 1 2 5 , U0 . 6 , 0 0 . 1 2
U0 . 7 , 0 0 . 1 0 5
Now, U, 00xx
,1UU 0ij i jk for ,1 ,00 U Uii jjt 
Putting 0j in equation (I) we get,
,1 1, 0 1, 0 1UU U Uii i i
,1 1, 0 1, 0 ,1 12U U U U Uii i ii   
,1 1,0 1,01UU U2ii i  
Now, for 1i 1,1 0,0 2,01UU U2
100 . 8 02
0.040.
munotes.in

Page 133

133

For 2i 2,1 1,0 3,01UU U2
 10.045 0.1052
0.075.

For 3i 3,1 2,0 4,01UU U2
 10.08 0.1202
0.1.

For 4i 4,1 3,0 5,01UU U2
 10.105 0.1252
0.115.



For 5i 5,1 4,0 6,01UU U2
 10.12 0.122
0.12.


For 6i 6,1 5,0 7,01UU U2
 10.125 0.1052
0.115.

Putting 2j in equation (II) we get,
,2 1 ,1 1 ,0UU UUii ii 
For 1i 1, 2 0,1 2,1 1,0UUUU
00 . 0 7 50 . 0 4 5 
0.03.
munotes.in

Page 134

134

For 2i 2,2 1,1 3,1 2,0UU U U
0.040 0.100 0.08
0.060.

For 3i 3,2 2,1 4,1 3,0UUUU
0.075 0.115 0.105
0.085.

For 4i 4,2 3,1 5,1 4,0UUU U
0.1 0.12 0.12 
0.1.

For 5i 5,2 4,1 6,1 5,0UUUU
0.115 0.115 0.125
0.105.

1, 2 2, 2 3, 2 4, 2 5, 2U0 . 0 3 , U0 . 0 6 , U0 . 9 8 5 , U0 . 1 , U0 . 1 0 5     .

7.9 Review :

In this chapter we have learn ,
 Classification of partial differential equation .
 Numerical methods of solv ing elliptic partial differential equation .
 Numerical methods by solving parabolic partial differential equations.
 Numerical methods of solving hyperbolic partial differential equations.

7.10 Unit End Exercise :

1) Solve the Laplace equations 22
22UU0xy at the interior mesh
p + s of the square region with boundary p + s shown in fig.






010203040506070U7 U8 U9
U4 U5 U6
U1 U2 U320
30
405060708090munotes.in

Page 135

135

2) The temperature U in the steady heat below in a square plate
bounded by 0, 0, 4xyy x  satisfies Laplace equation
22
22UU0xy.

3) Given one – dimensional heat 2
2UUxt  under the condition U0 , 0 U1 ,tt for 0t take 13h, 136k and use
explicit method to compute U for one time step only.

4) Using Schmidt method find the values of U,xt satisfying the
parabolic equation 2
24U Uxt subject to the conditions.

5) Given one – dimensional heat equation 2
2UUxt under the
conditions U, 0 4 0x  for 03x and U0 , 0 , U3 , 6 0tt for 0t take 1, 1hk and use Crank
Nicolson method to compute x for one time step only.

6) Solve by Crank Nicolson method 2
2UUxt 01 a n d 0xt 
given that U, 0 1 0 0xx for 01x and U0 , 0 U1 ,tt
for 0t take 11,24hk and compute U for one time step
only.

7) Given one dimensional wave equation 22226U Uxt under the
conditions U0 , 0 U5 ,tt and 32Ux , 0 5xx for 05x for 0tand U, 0 0tx for 05x.
munotes.in

Page 136

136

8) Given hyperbolic equation 2222UUxt subject to the condition U0 , 0 U1 ,tt for 0t and 3U , 0 sintxx for 01x
take 0.25h for 0.2k and use explicit method to compute U
for two time steps.
9) Derive the five point formula for Laplace’s equation.

10) What is Crank – Nicholson Method? Why is it known as implicit
method?

11) What is Bender – Schmidt recurrence equation? Derive the
formula.

12) Discuss the impact of size of the incremental width t for the
time variable on the solution of a heat flow equation.



munotes.in