Advance_numerical_method.pptx

Advanced Numerical Methods (S. A. Sahu) Code: AMC 51151

Syllabus We can divide our syllabus in following four major sections:

1. Method of solution of system of equations Solution of Non-linear Simultaneous Equations (Newton-Raphson Method)

Tridiagonal system of equations (Thomas algorithm)

2. Numerical Integration Evaluation of Double & Triple Integrals

(a) Constant Limit

Simpson’s

Rule

(b) Variable Limit

&

Trapezoidal

Rule

3. Methods of Interpolation  Newton’s F/W and B/W formula  Gauss’s F/W and B/W formula

 Bessel’s Formula  Stirling Formula  Lagrange’s Interpolation Formula  Newton’s Divided Difference Formula

4. Numerical Solution of Differential Equations O.D.E. •First Order & Higher Order Ordinary Diffe. Eqs. •Initial Value & Boundary Value Problems

P.D.E. •Laplace & Poisson Equations •Heat Conduction • Wave Equations

Solution of Non-linear Simultaneous Equations

(Newton-Raphson Method)

This Method is extension of N-R Method for Non-Linear Equation

• Method for Simultaneous equations: Suppose we are given the following Simultaneous equations

f ( x, y)  0

&

g ( x, y)  0

Why simultaneous..?

…(1)

We consider that ( x0 , y0 ) be the first approximation

And ( x0  h, y0  k ) be the real roots It means, we must have :

f ( x0  h, y0  k )  0   and   g ( x0  h, y0  k )  0 

….(2)

Now we expand the relation, we had in eq. (2) Using Taylor Series Expansion h2 f ( x0  h)  f ( x0 )  h f ( x0 )  f ( x0 )  ............. 2

 f   f  f ( x0  h, y0  k )  f ( x0 , y0 )  h    k    ....... (3)  x 0  y 0

(similarly for g ( x0  h, y0  k ) )

Now since h and k are very small, we can neglect the terms containing h2 (or k 2)and higher order terms…

• Now, eq.(2), eq.(3) and assumption suggests  f   f  f ( x0  h, y0  k )  0  f ( x0 , y0 )  h    k    x 0  y 0  g   g  g ( x0  h, y0  k )  0  g ( x0 , y0 )  h    k    x 0  y 0

..(4a)

..(4b)

• Or simply:  f   f  h    k     f ( x0 , y0 )  x 0  y 0

…5(a)

 g   g  h    k     g ( x0 , y0 )  x 0  y 0

..5(b)

Let us say

f ( x0 , y0 )  f 0 and

g ( x0 , y0 )  g0

Cramer’s Rule for or the solution of a system of linear equation

• We have from Eqs. 5(a) & 5(b)  f   f  h    k     f0  x 0  y 0  g   g   k     g0   x 0  y 0

h

Or in matrix form  f     x 0    g    x 0 

 f       y 0    g       y 0 

h   f0                   k    g0 

Now let us call

  f0   g  0 

 f     x 0    g    x 0   f       y 0    Dx  g       y 0 

 f       y 0    g       y 0   f     x 0   g      x 0

 D,

 f0  g0

   D y   

And hence we can have (from Cramer’s rule) Dy Dx h and k  D D

• This gives the next approximation, as

x1  x0  h

and

y1  y0  k

Using these values we can proceed further for next approximation (better value) This completes the method.

• Solve the given simultaneous equations (up to two iterations)

f ( x, y )  2 x  y  1  0 3

2

g ( x, y )  xy  y  4  0 3

take x0  1.2, y0 1.7

Q.2

f ( x, y )  x  3log x  y  0 2

g ( x, y )  2 x  xy  5 x  1  0 2

take x0  3.4 , y0  2.2

Solution of Q.2 f 0  3.4  3* log(3.4)  (2.2)  0.1545 2

g 0  2 * (3.4) 2  (3.4) * (2.2)  5* (3.4)  1  0.36 3 3  f   1.88235   1 1 x 3.4  x 0  g     4 x  y  5  6.4  x 0

 f     2 y  4.4  y 0  g      x  3.4  y 0

1.88235 D1   6.4

 0.1545 Dx1   0.36

1.88235 Dy1   6.4

 4.4   21.76001   3.4  h

Dx1 2.1093   0.09693 D1 21.76001

k

Dy1 1.6664   0.0766 D1 21.76001

 4.4   2.1093   3.4 

 0.1545  1.6664  0.36 

• And hence we have the second approximation as

x1  x0  h  3.4  0.09693  3.4969 y1  y0  k  2.2  0.0766  2.2766 • We repeat the process for next approximation

Tridiagonal system of equations (Thomas algorithm)

• In applications, many times we can have system of equations which gives a coefficient matrix of special structure, majority of zeros . Sparse matrices [Upper Triangular & Lower Triangular Matrices]

What, If we have a matrix such that except the principal diagonal and its upper and lower diagonal, all the elements are zero….???

Such a matrix is called TRIDIAGONAL MATRIX

Something which looks like…….

Tridiagonal Systems Tridiagonal Systems: • The non-zero elements are in the main diagonal, super diagonal and subdiagonal.

• aij=0 if |i-j| > 1

CISE301_Topic3

5 3  0  0 0

1 0 0 0  x1   b1       4 1 0 0  x2  b2  2 6 2 0  x3   b3      0 1 4 1  x4  b4  0 0 1 6  x5  b5 

25

Mathematically…. • An n×n matrix A is called a tridiagonal matrix if ai , j  0 whenever i  1  j or j  1  i

• The system of equations which gives rise to a tridiagonal coefficient matrix is called ….. Tridiagonal System of Equations

A tridiagonal system for n unknowns may be written as

ai xi 1  bi xi  ci xi 1  di i  1, 2,3...................n and

…(1)

(a1  cn  0)

B.V.P. with second order ODE, Heat conduction equations…. Generate equations like (1)

Now if we put i = 1, 2…n in eq. (1), we get

a1 x0  b1 x1  c1 x2  d1

..2(a)

a2 x1  b2 x2  c2 x3  d 2

..2(b)

a3 x2  b3 x3  c3 x4  d3

..2(c)

….. + …… +…….

= ……

Which in matrix form can be written as

b1 c1 0 0 0.....0   a b c 0 0 ....0   2 2 2  0 a3 b3 c3 0....0    ...............................  .............................c  n   0 0 0 0 0 an bn 

 x1   d1  x  d   2  2  x3   d3       ..  ..  ..  ..       xn   d n 

• Thomas algorithm (named after L. H. Thomas), is a simplified form of Gauss elimination method • Suppose we have a system like

b1 x1  c1 x2  d1 a2 x1  b2 x2  c2 x3  d 2 a3 x2  b3 x3  c3 x4  d3 And we have to find all the unknowns…

x1 , x2 , x3 ,...........xn

Then Thomas Algo. Gives the values of unknown in the order last to first..

It means first we shall get xn and then by back substitution, we will get xn 1 , xn  2 , xn 3 ,........ x1

 n  xn and

xi  i   i ( xi 1 )

i  n  1, n  2, .........3, 2,1

Where

di  ai i 1 d1 1  and i  b1 bi  aii 1

( for i  2,3,4....n)

and

ci c1 1  and i  b1 bi  aii 1

( for i  2,3,4....n)

Q.1 Solve the following system of equations

5 x1  2 x2  12 x1  5 x2  2 x3  9 x2  5 x3  2 x4  8 x3  5 x4  6

In matrix form

5 1  0  0 CISE301_Topic3

2 5 1 0

0 2 5 1

0   x1  12       0   x2   9   2   x3   8      5   x4   6  34

2 1   0.4 5

2  2   0.43478  0.4348 and so on .......... 4.6

1  2.4,  2  1.4348, 3  1.4381,  4  1.  4  x4  1 Now by back substitution we get x3  1, x2  1, x1  2

Section: 2 Finite Difference

Interpolation

Numerical Integration

Finite Differences

Let we have

y  f ( x)

..(1)

which is tabulated for equally spaced values..

x  x0 , x0  h, ...... x0  nh

And returns the following values

y0 , y1 , ........... yn It means

f ( x0 )  y0 f ( x0  h)  f ( x1 )  y1 . . . f ( x0  nh)  f ( xn )  yn

Then we can define the following Differences

• Forward Difference



DELTA

• Backward Difference



Nabla

• Central Difference



delta (lower case)

• Forward Difference First Forward Difference We define the FORWARD DIFFERENCE OPERATOR, denoted by  as

The expression f ( x  h)  f ( x) gives the FIRST FORWARD DIFFERENCE of and the operator  is called the FIRST FORWARD DIFFERENCE OPERATOR.

• Notation (First Forward Differences)

y0  y1  y0 y1  y2  y1 . . yr  yr 1  yr

Second Forward Differences are given by

 y0  y1  y0 2

 y1  y2  y1 2

. .  yr  yr 1  yr 2

• And hence in general the pth Forward Diff.

 yr   p

p 1

yr 1  

p 1

yr

..(2)

• Forward Difference Table

• Forward Diff. Terms

x  Argument y  function / Entry y0  Leading term y0 ,  2 y0 ,  3 y0 ....  Leading differences

• Backward differences are defined and denoted by following

y1  y1  y0 y2  y2  y1 . . yr  yr  yr 1

Following the concept of Forward Diff.

Make a table for Backward Difference up to 5th Difference

• Central Differences • The central Diff. operator following relation



is defined by the

 y1/ 2  y1  y0  y3/ 2  y2  y1 Define second third and other central differences And make a central diff. table

• Write the forward Diff. table if…

x : 10 20 30 40 y : 1.1 2.0 4.4 7.9

• Other operators

Shift Operator

E

Average Operator



Shift Operator is the operation of increasing of argument x by h and hence defined by

Ef ( x)  f ( x  h) E 2 f ( x )  f ( x  2h) E 3 f ( x)  f ( x  3h) .......... The inverse shift operator is defined by

1

E f ( x)  f ( x  h)

We have defined the following Differences • Forward Difference



DELTA

y0  y1  y0

• Backward Difference



Nabla

y1  y1  y0

• Central Difference



delta

 y1/ 2  y1  y0

• Shift Operator

E



• Average Operator





E 1 f ( x)  f ( x  h)  yx  ( y

h x 2

y

h x 2

)

Relations among the operators

(i )   E  1

or

(ii )   1  E 1 or (iii )   E

1/ 2

E

1/ 2

E 1  E 1  1   1 1/ 2 (iv)   ( E  E 1/ 2 ) 2

• Proof of (i)

yx  yx  h  yx  Eyx  yx  ( E  1) yx    E 1

( y0  f ( x0 )  yx  h  f ( x  h)  Eyx )

Interpolation

Suppose we have following values for y=f(x)

x : x0

x1 x2 ... xn

y : y0

y1

y2 ... yn

• Interpolation is the technique of estimating the value of a function for any intermediate value of the independent variable • The process of computing the value of the function outside the range is called Extrapolation

Interpolation Equal spaced Unequal Differences Forward Interpolation Newton’s Divided Diff. (Newton’s Fwd. Int.) Backward Interpolation Lagrange’s Interpolation (Newton’s Bwd. Int.) Central Diff. Interpolation • Gauss Forward • Gauss Backward • Stirling Formula • Bessel’s Formula • Everett’s Formula

Newton’s Forward Interpolation

Suppose the function

y  f ( x)

is tabulated for equally spaced values..

x  x0 , x0  h, ...... x0  nh And it takes the following values

y0 , y1 , ........... yn

• Suppose we have to calculate f(x) at x=x0+ph (p is any real number), then we have

y p  f ( x0  ph)  E p f ( x0 )  (1  ) p y0 Now we expand the term (1  ) by Binomial theorem p

  p ( p  1) 2 p( p  1)( p  2) 3 y p  (1  ) p y0  1  p      .... y0 2 3   or p ( p  1) 2 p( p  1)( p  2) 3 p( p  1)...[ p  (n  1)] n y p  y0  py0   y0   y0  ....   y0 2 3 n ........(1)

(sinceif y  f ( x) is polynomial of nth degree then  n 1 y0 and other terms will be zero )

The above blue colored equation (1) is Newton’s Forward Interpolation Formula Because it contains y0 and forward differences of y0

Example Estimate f (3.17)from the data using Newton Forward Interpolation. x: 3.1 f(x): 0

3.2 0.6

3.3 1.0

3.4 1.2

3.5 1.3

Re member : x ( x0  ph)  p 61

x  x0 h

Solution First let us form the difference table x 3.1

y 0

y

2y

3y 4y

0.6 3.2 0.6

- 0.2 0.4

3.3 1.0

0 - 0.2

0.2 3.4 1.2 3.5 1.3 Here x0 = 3.1, x = 3.17, h = 0.1.

62

0.1 -0.1

0.1

0.1

Solution p=

x  x0 h

=

0.07 0 .1

= 0.7

Newton forward formula is: P(x)

= y0 + py0+ p( p2! 1) 2y0+ p(p  13)(! p  2) 3y0+

p(p  1)( p  2)( p  3) 4!

P(3.17)=0+0.70.6+ 0.7(02.7  1) (-0.2)+ 0.7(0.7 61)(0.7  2) 0+ 0.7(0.7  1)(0.7  2)(0.7  3) 24

 0.1

= 0.4384 Thus f(3.17) = 0.4384.

63

4y0

Newton’s Backward Interpolation

Suppose the function

y  f ( x)

is tabulated for equally spaced values..

x  x0 , x0  h, ...... x0  nh And it takes the following values

y0 , y1 , ........... yn

• Suppose we have to calculate f(x) at x=xn+ph (p is any real number, may be -ve), then we have

y p  f ( xn  ph)  E f ( xn )  E p

( 1)(  p )

 (1  )

p

we have (1  ) for E 1 Now we expand the term (1  ) p by Binomial theorem

Re member : x ( xn  ph)  p

x  xn h

yn

NEWTON GREGORY BACKWARD INTERPOLATION FORMULA Taking p = x hx , we get the interpolation formula as: n

P(xn+ph) = y0 + pyn + p(p  1)(p  2) 3!

p(p  1) 2!

 yn +

3yn +  + p (p  1) (p  2)  (p  n  1) n!

66

2

nyn

Example Estimate f(42) from the following data using newton backward interpolation. x:

20

25

30

35

40

45

f(x):

354

332

291

260

231

204

67

Solution The difference table is: x f 20 354

f

2f

3f 4f 5f

- 22 25 332

and p = - 0.6 - 19

- 41 30 291

29 10

- 31 35 260

-37 -8

2 - 29

40 231

Here xn = 45, h = 5, x = 42

45 8

0 2

- 27 45 204

68

Solution Newton backward formula is: P(x) = 3 4 p(p + 1)(p + 2)(p + 3) + 2) yn+pyn+ p(p2!+ 1) 2yn+ p(p + 1)(p  y +  yn + n 3! 4! p( p  1)( p  2)( p  3)( p  4) 5!

5yn

(-0.6)(0.40(1.4) P(42)=204+(-0.6)(-27)+ (-0.6)(0.4) 2+ 0+ 2 6 (-0.6)(0.4)(1.4)(2.4) 24

Thus,

8+

(-0.6)(0.4)(1.4)(2.4)(3.4) 120

 45 = 219.1430

f(42) = 219.143 69

Central Difference Formula • Gauss Forward • Gauss Backward • Stirling Formula • Bessel’s Formula • Everett’s Formula

• If y=f(x) is any function, which takes values

x0  2h, x0  h, x0 , x0  h, x0  2h and returns corresponding values

y2 , y1 , y0 , y1 , y2 • Then the difference table can be written as

Central Difference table

x x0  2h

y y2

1st Diff.

2nd Diff.

3rd Diff.

4th Diff.

y2 (  y3/2 ) x0  h

 2 y2 (  2 y1 )

y1 y1 (  y1/2 )

x0

 3 y2 (  3 y1/2 )  2 y1 (  2 y0 )

y0 y0 (  y1/2 )

x0  h

 3 y1 (  3 y1/2 )  2 y0 (  y1 )

y1 y1 (  y3/2 )

x0  2h

 4 y2 (  4 y0 )

y2

Table-1

1. GAUSS FORWARD INTERPOLATION FORMULA p p  p  1 3  p  1 4   y-1+    y-2 + P(x) = y0+   y0+   2y-1+  1 2 3 4        

p=

 p  2 5    y-2 +  where 5  

x - x0 p = p( p -1) ( p - 2) ( p -r +1) and r h r!

()

OR y p  y0  py0 

p ( p  1) 2 ( p  1) p( p  1) 3 ( p  1) p( p  1)( p  2) 4  y1   y1   y2  ..... 2! 3! 4!

…..(A) 74

Derivation of Gauss Fwd. Int. Formula We have Newton’s Forward Int. Formula y0  py0 

p( p  1) 2 p( p  1)( p  2) 3 p( p  1)...[ p  (n  1)] n  y0   y0  ....   y0 ...(1) 2 3 n

We have following relations (from table)

• Note: (1) Gauss forward formula contains odd differences just below the central line and even differences on the central line

(2) Formula is useful to interpolate the value of y for 0
Example Find f(30) from the following table values using Gauss forward difference formula: x:

21

25

29

33

37

F(x):

18.4708

17.8144

17.1070

16.3432

15.5154

78

Solution The difference table is f

x

f

21

18.4708

 2f

 3f

 4f

-0.6564 25

17.8144

- 0.0510 -0.7074

29

17.1070

- 0.0054 0.0564

-0.7638 33

16.3432

- 0.0022 - 0.0076

- 0.0640 -0.8278

37

15.5154

p P(x) = y0 +   y0 + 1

p 2    y-1 +  2

 p  1 3    y-1 + 3  

 p  1 4    y-2 +  4  

f (30) = 16.9217 79

2. Gauss Backward interpolation formula Again starting with Newton’s Forward Int. Formula y0  py0 

p( p  1) 2 p( p  1)( p  2) 3 p( p  1)...[ p  (n  1)] n  y0   y0  ....   y0 ...(1) 2 3 n

we have y0  y1   2 y1  y0   2 y1  y1 similarly  2 y0  3 y1   2 y1 and 3 y0   4 y1   3 y1

Also we can have

3 y1  3 y2   4 y2  3 y1   4 y2   3 y2 And similarly other differences.. Putting these values we finally get

y p  y0  py1 

p( p  1) 2 ( p  1)( p  1) 3 ( p  1)( p  1)( p  2) 4  y1   y2   y2  ..... 2 3 4

or in central difference notation y p  y0  p y1/ 2 

p( p  1) 2 ( p  1)( p  1) 3 ( p  1)( p  1)( p  2) 4  y0   y1/ 2   y0  ..... 2 3 4

this is Gauss Backward Formula

• Note: (1) Gauss Backward formula contains odd differences just above the central line and even differences on the central line

• Formula is useful to interpolate the value of y for p lying between -1 & 0

Example Estimate cos 51 42' by Gauss backward interpolation from the following data: x: cos x:

50 0.6428

51

52

53

54

0.6293

0.6157

0.6018

0.5878

83

Solution x0 = 52, x = 51 42' = 51.7, h = 1 x  x0 h

P=

=

51.7  52 1

= -0.3

The difference table is X

y

50

0.6428

y

2y

3y

4y

- 0.0135 51

0.6293

- 0.0001 - 0.0136

x0 = 52

0.6157

- 0.0002 - 0.0003

- 0.0139 53

0.6018

- 0.0002 - 0.0001

- 0.0140 54

0.0004

0.5878

84

Solution The Gauss backward formula is: P(x) = y0 +

p   y-1 1

+

 p  1 2    y-1  2 

P(51.7) = 0.6198

85

+

 p  1 3    y-2  3 

+

 p  2 4    y-2  3 

3. STIRLING’S FORMULA • This formula gives the average of the values obtained by Gauss forward and backward interpolation formulae. For using this formula we should have – ½ < p< ½. We can get very good estimates if - ¼ < p < ¼. The formula is:

P(x) =

 y 0  y 1  p 2 y0+p  + 2   2!

 y-1+ 2

86

p(p 2  1) 3!

 3 y 1  3 y  2  p 2 (p  1) 4   +  y-2 2 4 !  

+ ...

• Stirling Formula (Derivation) [Mean of Gauss Fwd. & Gauss Bkwd]

• Eq. (3) is Stirling Formula

• Note: • Stirling Formula contains means of odd differences just above and below the central line and even differences on central line

Example • Using Sterling Formula estimate f (1.63) from the following table: x: 1.50 1.60 1.70 1.80 1.90 f(x):17.609 20.412 23.045 25.527 27.875 Solution X0 = 1.60, x = 1.63, h = 0.1 p=

1.63 – 1.60 = 0.3 0.1

89

Solution • Difference table x y 1.50 17.609

y

 2y

 3y

 2y

2.803 x0= 1.60 20.412 1.70 23.045 1.80 25.527

2.633 2.482 2.348

-0.170 -0.151 -0.134

1.90 27.875

90

0.019 0.017

-0.02

4. BESSELS’ INTERPOLATION FORMULA We have Gauss Fwd. Interpolation Formula as y p  y p  py0 

p ( p  1) 2 ( p  1) p( p  1) 3 ( p  1) p( p  1)( p  2) 4  y1   y1   y2  ..... 2! 3! 4!

Now we use following relations  2 y0   2 y1   3 y1



 2 y1   2 y0   3 y1

and similarly  4 y2   4 y1   5 y2 and other even differences

• We get the following Bessel’s Formula 1 ( p  ) p( p  1)  p ( p  1)  y1   y0  ( p  1) p ( p  1)( p  2)  4 y2   4 y1 3 2 y p  y0  py0    y1   ...  2! 2 3! 4! 2   2

2

….(B)

BESSELS’ INTERPOLATION FORMULA/other way • This formula involves means of even difference on and below the central line and odd difference below the line. • The formula is P(x) =

+

y 0 + y1 æ 1 ö p (p – 1) + çp – ÷ Dy 0 + è 2 2ø 2!

1ö æ p ç ÷ p (p – 1) è 2ø 3!

3y-1 +

æD2 y- 1 + D2 y0 ö ç ÷ 2 è ø

(p + 1) (p – 1) (p – 2) 4!

93

æD4 y –2 + D4 y –1 ö ç ÷ 2 è ø

BESSELS’ INTERPOLATION FORMULA Note • •

When p = 1/2 , the terms containing odd differences vanish. Then we get the formula in a more simple form:

y 0 + y1 p (p – 1) æD2 y- 1 + D2 y0 ö P(x) = + ç ÷ 2 2! 2 è ø

(p + 1) p(p – 1)(p – 2) æD4 y –2 + D4 y –1 ö + + ç ÷ 4! 2 è ø

94

Example Use Bessels’ Formula to find (46.24)1/3 from the following table of x1/3.

X: X1/3:

41 45 49 53 3.4482 3.5569 3.6593 3.7563

95

Solution x0 = 45, x = 46.24, Difference table is: x

Y

y

h = 4 p = 0.31  2y

 3y

41 3.4482 0.1087 x0 =

45 3.5569

-0.0063 0.1024

49 3.6593

-0.00091 -0.0054

0.0970 53 3.7563

Applying Bessels’ formula, we get (46.24)1/3 = 3.5893 96

5. Laplace-Everett’s Formula Gauss’s forward interpolation formula is y p  y0  py0 

p  p  1 2  p  1 p  p  1 3 y   p  1 p  p  1 p  2 4 y   p  2 p  1 p  p  1 p  2 5 y  ......  y1  1 2 2 2! 3! 4! 5!

We eliminate the odd differences in above equation by using the relations y0  y1  y2 , 3 y1   2 y0   2 y1 , 5 y2   4 y1   4 y2

And then to change the terms with negative sign, putting p  1  q , we obtain y p  qy0 

q  q 2  12  3!

 y1  2

q  q 2  12  q 2  22  5!

 y2  ......  py1  4

p  p 2  12  3!

 y0  2

p  p 2  12  p 2  22  5!

…(E) This is known as Laplace-Everett’s formula.

 4 y1  .....

Choice of an interpolation formula The following rules will be found useful:  To find a tabulated value near the beginning of the table, use Newton’s forward formula.  To find a value near the end of the table, use Newton’s backward formula.  To find an interpolated value near the centre of the table, use either Stirling’s or Bessel’s or Everett’s formula.

Choice of central difference  If interpolation is required for p lying between 1 1  4 and 4 , prefer Stirling’s formula.  If interpolation is desired for p lying between 1 3 4 and 4 , use Bessel’s or Everett’s formula.  Gauss’s Backward interpolation formula is used to interpolate the value of “y” for a negative value of p lying between -1 and 0.

 Gauss’s forward interpolation formula is used to interpolate the value of “y” for p (0
Bessel’s and Everett’s central difference interpolation formula can be obtained by each other by suitable rearrangements.