# Comparing Numerical Solution to a Second Order Boundary Value Problem by Variational Techniques

DOI : 10.17577/IJERTV10IS090175

Text Only Version

#### Comparing Numerical Solution to a Second Order Boundary Value Problem by Variational Techniques

Shubham V. Deshmukh

Research Scholar Department of Mechanical Engineering, Shri Ramdeobaba College of Engineering and Management, Nagpur, Maharashtra, India 440013.

Vishal V. Shukla

Associate Professor Department of Mechanical Engineering, Shri Ramdeobaba College of Engineering and Management, Nagpur, Maharashtra, India 440013.

Abstract This Research paper presents the techniques to determine the solution for a second-order Boundary Value Problem (BVP) that is frequently witnessed in engineering phenomena, by multiple variational methods. Further, the solutions obtained by these methods are compared with the exact solution of the second-order differential equation. The methods that are illustrated here are Subdomain Method, Collocation Method, Least Square Method, Petrov Galerkin Method, Galerkin Method, Rayleigh-Ritz Method and the Finite Difference Method. The solutions obtained by numerical methods seem to be fairly in good agreement with the exact solutions. Further, the comparison and average percentage error between solutions obtained by Exact Method and above- stated Numerical Methods have been carried out in MS-Excel. The interval or spacing of 0.1 is considered as the step size within the specified domain of 0 to 1. It is observed and confirmed that to achieve higher accuracy, the intervals or step size must be reduced. This analysis reveals that the Finite Difference Method (FDM) gives the most accurate solution among the other Numerical Methods considered in this study. And therefore, most of the Finite Element Analysis (FEA) and Computational Fluid Dynamics (CFD) software includes FDM. The accuracy level attained by various methods from high to low is Finite Difference Method, Rayleigh-Ritz Method, Galerkin Method, Petrov Galerkin Method, Collocation Method, Subdomain Method, Least Square Method. Galerkin Method and Rayleigh-Ritz Method give identical solutions. The study also concludes that the Finite Difference Method is high in accuracy and Least Square Method appears to be the poorest in achieving the accuracy.

Keywords Boundary Value Problem, Approximate Method, Weighted Residue Methods.

1. INTRODUCTION

It is said that the backbone of any physical system is its differential equation. Ordinary differential or partial differential equation is used to model the physical system [1]. The ultimate aim of numerical analysis is to find the approximate solution to the differential equation when the exact solution is difficult to obtain. Both finite difference method and variational methods of approximation attempts are performed to find the approximate solution in the form of a linear combination of suitable approximation function

difference methods with the exact method by comparing their solutions to the second ordered differential boundary value problem.

2. SOLVING SECOND ORDER DIFFERENTIAL BOUNDARY VALUE PROBLEM BY VARIOUS

2

2

NUMERICAL METHODS

Given Differential Equation is: y + y + 1 = 0 ; Domain:

x2

0 x 1

Boundary Conditions: i. Y(0) = 0 ; ii. Y(1) = 0. By Exact method:

By Exact Method, the solution is given by:

Yexact = C1 cos(x) + C2 sin(x) 1 (1) Applying boundary condition i to equation 1:

0 = C1 1

C1 = 1

Applying boundary condition ii to equation 1:

0 = (1) cos(1) + C2 sin(1) 1

C2 = 0.5463

Putting values of C1 and C2 in equation 1:

Yexact = cos(x) + (0.5463)sin(x) 1 (2) By Variational methods of approximations:

Given Differential Equation is:

y

y

2

x2 + y + 1 = 0. (3)

Let Y Yexact = Yapproximate.

Let C1, C2, C3 and C4 be the constants. Let us assume a polynomial:

Yapproximate = C1 + C2x + C3x2 + C4x3 (4) Applying boundary condition i to equation 4:

C1 = 0

Applying boundary condition ii to equation 4:

0 = 0 + C2 + C3 + C4 C2 = C3 C4

Putting values of C1 and C2 in equation 4:

Yapproximate = (C3 C4)x + C3x2 + C4x3 (5) Yapproximate = C3x C4x + C3x2 + C4x3 (6) Yapproximate = C3(x2 x) + C4(x3 x) (7) Differentiating equation 7 with respect to x:

Yapproximate = C3(2x 1) + C4(3×2 1) (8)

and undetermined coefficients [2]. The given differential

x

Subs ting

and

in equation 5:

2y

2y

titu

C3 = A

C4 = B

equation is + y + 1 = 0, comparing with A +

2 2

2 2

x2 x2

B y + C y + D y + E y + Fy + G = 0. Here, A = 1,

Yapproximate = (A + B)x Ax2 Bx3 (9)

Yapproximate = A(x x2) + B(x x3) (10)

xu

u2

x u

Differentiating equation 9 with respect to x:

B = C = 0. So, B2 4AC = 0. Hence, the given

differential equation is parabolic [3]. This research compares

Yapproximate x

= (A + B) 2Ax 3Bx2 (11)

the accuracy of various weighted residue methods and finite

Differentiating equation 11 with respect to x:

2Yapproximate = 2A 6Bx (12)

x2

Substituting equations 9 and 12 in equation 3:

2A 6Bx + (A + B)x Ax2 Bx3 + 1 = R

Where R = Residue.

(2A + 1) + (A 5B)x Ax2 Bx3 = R (13)

A(x x2 2) + B(x3 5x) + 1 = R (14)

Now, we have Weighted Residue Statement for the given differential equation is:

Wi R dx = 0 ; 0 x 1 (15)

A. Subdomain method

1. Least square method

In the least square method, differentiate the equation of residue i.e. equation (13) with respect to A and B to get weight function w2 and w1 respectively.

w1 = R = x x2 2

A

w2 = R = x3 5x

B

Solving both cases as given below:

Wi R dx = 0 ; 0 x 1 ; i.e. i = 1

1

1

1

0 ( w1)Rdx = 0

0

0

(x x2 2)[(2A + 1) + (A 5B)x Ax2 Bx3]dx = 0

As the given domain is 0 x 1. In the subdomain

101 A 119 B = 11

(22)

method, we divide the domain into two parts and we will solve each part individually.

In the subdomain, we take the weight function as a unity.

i.e. w1 = w2 = 1

1

30 20 6

Wi R dx = 0 ; 0 x 1 ; i.e. i = 2

1

1

1

0 ( w2)Rdx = 0

0

0

(5x x3)[(2A + 1) + (A 5B)x Ax2 Bx3]dx = 0

Solving for the first part of domain = 0 x

; i.e. i = 1

101 A 220 B = 11

(23)

1 2 20 21 4

2(w1)R dx = 0

Solving equations 22 and 23 simultaneously:

55

0

0

1

2(1) ((2A + 1) + (A 5B)x Ax2 Bx3)dx = 0

A =

101

and B = 0.

0 Substituting the values of A and B in equation 9:

11 A 41 B = 1

(16)

55 55 2

12 64 2

Y least square method = x x

(24)

Solving for second part of domain = 1 x 1 ; i.e. i = 2

2

1

1 ( w2)Rdx = 0

101

2. Petrov Galerkin method

101

1

1

2

1 (1) ((2A + 1) + (A 5B)x Ax2 Bx3)dx = 0

In the Petrov Galerkin Method we take weight function (w1 and w2) in combination such as (1 and x) or (x and x2) or (1

2 2 3 3 9

11 A 135 B = 1

(17)

and x ) or (1 and x ) or (x and x ) respectively.

12 64 2

Solving equations 16 and 17 simultaneously:

Here, we will solve by using x and x2 as a weight function.

w1 = x

A = 6

11

and B = 0.

w2 = x2

Substituting the values of A and B in equation 9:

Ysubdoain method = 6 x2 + 6 x (18)

11 11

B. Collocation method

Solving both cases as given below:

Wi R dx = 0 ; 0 x 1 ; i.e. i = 1

1

1

1

0 ( w1)R dx = 0

0

0

(x)[(2A + 1) + (A 5B)x Ax2 Bx3]dx = 0

In the collocation method, we take residue at gauss point 0.

11 A 28 B = 1

(25)

i.e. R(x) = 0.

A domain is 0 x 1. We will take two gauss points x =

(1) and x = (2), where the value of residue is 0.

12 15 2

Wi R dx = 0 ; 0 x 1 ; i.e. i = 2

1

1

1

0 ( w2)R dx = 0

3 3 (x2) [(2A + 1) + (A 5B)x Ax2 Bx3]dx = 0

Solving both cases as given below:

1

0

37 A 17 B = 1

(26)

R (x = ) = 0

60 12 3

3 Solving equations 25 and 26 simultaneously:

Putting x = 1 in equation 14:

3

16 46

A = 310

531

and B = 10

531

A B = 1 (19) Substituting the values of A and B in equation 9:

9 27

100

310 2

10 3

R (x = 2) = 0

3

Putting (x = 2) in equation 14:

3

Ypetrov galerkin method =

E. Galerkin method

177

x x

531

+ x

531

(27)

16 A 98 B = 1 (20)

In Galerkin Method, differentiate equation of Y

i.e.

9 27

approximate

Solving equations 19 and 20 simultaneously:

equation (10) with respect to A and B to get weight function

A = 9

16

Sub

and B = 0.

ting the values of A and B in equation 9:

w1 and w2 respectively.

w1 = yapproximate = (x x2)

stitu

Ycollocation method = 9 x 9 x2 (21)

A

w2 = yapproximate = (x x3)

16 16

B

Solving both cases as given below:

Wi R dx = 0 ; 0 x 1 ; i.e. i = 1

1

0 ( w1 ) R dx = 0

1

1

(x x2)[(2A + 1) + (A 5B)x Ax2 Bx3]dx = 0

Wi R dx = 0 ; 0 x 1 ; i.e. i = 2

0

3 A 9 B = 1

(28)

(wi) [

2y

+ y + 1 ] dx = 0

1

1

10 20 6

0 x2

Wi R dx = 0 ; 0 x 1 ; i.e. i = 2

1 2y 1 1

(wi) dx + (wi) ydx + (wi) dx = 0

1 0 x2 0 0

1

1

0 ( w2)R dx = 0

(x x3)[(2A + 1) + (A 5B)x Ax2 Bx3]dx = 0

1 2y 1 1

0 (wi) x2 dx = 0 (wi) ydx 0 (wi) dx

2

0 Applying integration by parts to term (1(wi) y dx ),

2

2

9 A 76 B = 1

(29)

0 x2

20 105 4

considering u = wi and v = y , where the formula is

Solving equations 28 and 29 simultaneously:

A = 5 and B = 0.

9

x2

u

u v dx = u v dx [ v dx ] dx.

x

Substituting the values of A and B in equation 9:

dy 1 dwi dy

Ygalerkin method = 5

9

x 5 x2

9

(30)

[ wi

]1 (

0

0

dx 0 dx

Ã— ) dx dx

1 1

1. Rayleigh-Ritz method

= (wi) y dx (wi) dx

In Rayleigh-Ritz Method we take Residue (R) as a given

dy 1

0 0

dwi dy 1

differential equation.

[ wi

]1 = (

0

0

dx dx

Ã— ) dx (wi) ydx dx

So,

R = 2y + y + 1. x2

0 0

1

(wi) dx

In Rayleigh-Ritz Method, we take the coefficient of C3 and

C4 from equation 7 as weight function W1 and W2 respectively.

0

0

Putting i = 2: dy

0

1 dw2 dy 1

W1 =

yapproximate

(C3)

= (x2

x)

[ w2

]1 = (

dx 0 dx

Ã— ) dx ( w2)ydx dx 0

W2 = yapproximate = (x3 x)

(C4)

d(w1)

= (2x 1)

[(x3 x)

1

( w2)dx

0

1

1

y ]1 = ( 3×2 1)(C3 (2x 1) + C4(3×2

dx x 0 0

1

1

d(w2)

dx

= (3×2

1)

1))dx ( x3 x)(C3( x2 x) + C4(x3 x)dx

0

0

1

1

( x3 x) dx

Solving both cases as given below:

Wi R dx = 0 ; 0 x 1 ; i.e. i = 1

0

9 C3 + 76 C4 = 1

(32)

1 2y

20 105 4

0 (wi) [

x2

+ y + 1 ] dx = 0

Solving equations 31 and 32 simultaneously:

5

1 2y 1 1

C3 =

and C4 = 0.

0 (wi) x2 dx + 0 (wi) ydx + 0 (wi) dx = 0

9

Subs

g the values of C

and C

in equation 7:

2 titutin 3 4

( ) ( )

( ) ( )

1 y 1 1

wi 2 dx = wi ydx (wi) dx

Y rayleigh ritz method = 5 (x2 x) (33)

0 x 0 0

2

9

9

Applying integration by parts to term (1(wi) y dx ) ,

0 x2

2

2

considering u = wi and v = y , where the formula is

x2

u

u v dx = u v dx [ v dx ] dx.

2. Finite difference method (FDM)

In the Finite Difference Method, we have to define the nodes at a point where the value is to be determined. As Domain is

1 x 0 x 1, we will take a step size of 0.1 i.e. h = 0.1. As

dy dwi dy

[ wi ]1 ( Ã— ) dx 2 1

dx 0

0 dx dx

h = 0.1 ; Therefore h

= . Let Y0, Y1, Y2, Y3, Y4, Y5,

100

1 1 Y6, Y7, Y8, Y9, Y10 be the values of Y at a particular node

= (wi) y dx (wi) dx

when value of x is 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,

1 0 1 0 1 respectively.

dy dwi dy

[ wi ]1 = ( Ã— ) dx (wi) ydx

According to given boundary conditions i and ii:

dx 0

0 dx dx 0

1

Y0 = Y(0) = 0, Y10 = Y(1) = 0.

For the value of an unknown node we have a formula:

(wi) dx

0

Putting :

[yi+1]+[yi1]2yi + Y + 1 = 0.

2 i

2 i

h

i = 1

dy 1 dw1 dy 1

For the value of i equals to 1, 2, 3, 4, 5, 6, 7, 8, 9 equations

[ w1

]1 = ( Ã—

0

0

dx 0 dx

1

) dx ( w1)ydx dx 0

are:

199Y1 100Y2 = 1 (34)

100Y1 199Y2 + 100Y3 = 1 (35)

( w1)dx

0

0 0

0 0

[(x2 x) y]1 = 1(2x 1)(C3(2x 1) + C4(3×2

x

1))dx 1(x2 x)ydx 1(x2 x)dx

100Y2 199Y3 + 100Y4 = 1 (36)

100Y3 199Y4 + 100Y5 = 1 (37)

100Y4 199Y5 + 100Y6 = 1 (38)

100Y5 199Y6 + 100Y7 = 1 (39)

0 0

3 C3 + 9 C4 = 1

(31)

100Y6 199Y7 + 100Y8 = 1 (40)

10 20 6

100Y7 199Y8 + 100Y9 = 1 (41)

100Y8 199Y9 = 1 (42)

3. SOLUTION OF SECOND-ORDER DIFFERENTIAL BOUNDARY VALUE PROBLEM BY

VARIOUS NUMERICAL METHODS

Substituting the values of x in equations 2, 18, 21, 24, 27,

30, 33 and solving equations 34, 35, 36, 37, 38, 39, 40, 41, 42 simultaneously we get the results shown in table 1. The solutions of second order boundary value problem by various numerical methods at 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6,

0.7, 0.8, 0.9, 1.0 is mentioner under table 1.

The graphical representation of solutions of second-order differential boundary value problem by various numerical methods has been shown in figure 1 in scatter with smooth

line form. From figure 1, it can be concluded that the considered differential equation is parabolic. The parabolic nature of the considered second-order differential equation has been explained in the introduction portion of this research article. But, as the solutions by various numerical methods are in good agreement with each other the curves in figure 1 is overlapping on each other and the vision of graphs of each numerical method is not much clear. To overcome this issue in figure 2 bar graph representation is shown.

Table 1. Solution of second-order differential BVP by various numerical methods

 Value of y Value of x Exact Subdomain Collocation Least Square Petrov Galerkin Galerkin FDM Rayleigh Ritz 0.0 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.1 0.04954 0.04909 0.05063 0.04901 0.05068 0.05000 0.04959 0.05000 0.2 0.08860 0.08727 0.09000 0.08713 0.08979 0.08889 0.08868 0.08889 0.3 0.11678 0.11455 0.11813 0.11436 0.11746 0.11667 0.11689 0.11667 0.4 0.13380 0.13091 0.13500 0.13069 0.13379 0.13333 0.13393 0.13333 0.5 0.13949 0.13636 0.14063 0.13614 0.13889 0.13889 0.13962 0.13889 0.6 0.13380 0.13091 0.13500 0.13069 0.13288 0.13333 0.13393 0.13333 0.7 0.11678 0.11455 0.11813 0.11436 0.11588 0.11667 0.11689 0.11667 0.8 0.08860 0.08727 0.09000 0.08713 0.08798 0.08889 0.08868 0.08889 0.9 0.04954 0.04909 0.05063 0.04901 0.04932 0.05000 0.04959 0.05000 1.0 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000

#### Fig. 2 Bar graph representation of the solution of second-order differential boundary value problem by various numerical methods

4. ERROR ANALYSIS FOR THE SOLUTIONS OF VARIOUS NUMERICAL METHODS WITH RESPECT

TO EXACT METHOD

For error analysis, the solutions by the exact method are considered to be real without any bias i.e. with zero error. For calculating the error in the solutions of numerical methods, the solution difference of the exact method and numerical method is estimated. If in some cases difference comes negative then modulus is applied. Table 2 indicates

the error in the solution of various numerical methods with respect to exact method solutions. Figure 3 indicates a graphical representation of error analysis for the solutions of various numerical methods with respect to the exact method.

The assumptions considered during estimating error analysis is as follows:

1. Consider Exact Method has zero error.

2. Error in values of Yapproximate is estimated by formula: Error in values of Yapproximate = | (Yexact method Ynumerical method)

|

Table 2. Error analysis for the solutions of various numerical methods with respect to exact method

 Difference in value of y x Exact Subdomain Collocation Least Square Petrov Galerkin Galerkin FDM Rayleigh Ritz 0.0 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.1 0.00000 0.00045 0.00109 0.00053 0.00114 0.00046 0.00005 0.00046 0.2 0.00000 0.00133 0.00140 0.00147 0.00119 0.00029 0.00008 0.00029 0.3 0.00000 0.00223 0.00135 0.00242 0.00068 0.00011 0.00011 0.00011 0.4 0.00000 0.00289 0.00120 0.00311 0.00001 0.00047 0.00013 0.00047 0.5 0.00000 0.00313 0.00114 0.00335 0.00060 0.00060 0.00013 0.00060 0.6 0.00000 0.00289 0.00120 0.00311 0.00092 0.00047 0.00013 0.00047 0.7 0.00000 0.00223 0.00135 0.00242 0.00090 0.00011 0.00011 0.00011 0.8 0.00000 0.00133 0.00140 0.00147 0.00062 0.00029 0.00008 0.00029 0.9 0.00000 0.00045 0.00109 0.00053 0.00022 0.00046 0.00005 0.00046 1.0 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000

Fig. 3 Graphical representation of error analysis for the solutions of various numerical methods with respect to the exact method

5. PERCENTAGE ERROR ANALYSIS FOR THE SOLUTIONS OF VARIOUS NUMERICAL METHODS

WITH RESPECT TO EXACT METHOD

In this section, the percentage error for solutions of numerical methods with respect to the exact method at eleven points within the limit of 0 to 1 is estimated. Along with this, the average percentage error in the solution of each numerical method is also calculated. It has been assumed that the solution obtained by the exact method is a true solution with zero percentage of error. Table 3 indicates percentage error analysis for solutions of numerical methods with respect to the exact method. The average percentage error with respect to exact method in solutions of subdomain method, collocation method, least square method, petrov galerkin method, galerkin method, finite difference method, rayleigh ritz method is 1.38199 %, 1.13488 %, 1.51389 %, 0.66057 %, 0.34850 %, 0.07803 %, 0.34850 % respectively.

Table 3 Percentage error analysis for solutions of numerical methods with respect to the exact method

 Percentage error in value of y with respect to exact method (%) x Exact Subdomain Collocation Least Square Petrov Galerkin Galerkin FDM Rayleigh Ritz 0.0 0 0 0 0 0 0 0 0 0.1 0 0.90835 2.20024 1.06984 2.30117 0.92854 0.10092 0.92854 0.2 0 1.50112 1.58013 1.65914 1.34311 0.32731 0.09029 0.32731 0.3 0 1.90957 1.15601 2.07227 0.58229 0.09491 0.09419 0.09491 0.4 0 2.15994 0.89686 2.32436 0.00747 0.35127 0.09715 0.35127 0.5 0 2.24388 0.81726 2.40160 0.43013 0.43013 0.09319 0.43013 0.6 0 2.15994 0.89686 2.32436 0.68759 0.35127 0.09715 0.35127 0.7 0 1.90957 1.15601 2.07227 0.77067 0.09419 0.09419 0.09419 0.8 0 1.50112 1.58010 1.65914 0.69977 0.32731 0.09029 0.32731 0.9 0 0.90835 2.20024 1.06984 0.44408 0.92854 0.10092 0.92854 1.0 0 0 0 0 0 0 0 0 Average Percentage Error (%) 0 1.38199 1.13488 1.51389 0.66057 0.34850 0.07803 0.34850

Figure 4 and figure 5 indicate scatter with smooth lines and bar graph representation of percentage solution error by numerical methods with respect to exact method respectively.

#### Fig. 5 Bar graph representation of percentage solution error by numerical methods with respect to the exact method

6. RESULT AND DISCUSSION

The chosen second order differential equation has been solved by Subdomain Method, Collocation Method, Least Square Method, Petrov Galerkin Method, Galerkin Method, Rayleigh-Ritz Method and Finite Difference Method. The solutions obtained by these numerical methods are compared with the exact solutions. The average percentage error has been estimated between solutions of Numerical Methods and Exact Methods. Figure 6 indicates a bar graph of average percentage error by numerical methods with respect to the exact method and table 4 indicates the ranking of various numerical methods on the basis of average percentage error with the exact method.

#### Fig. 6 Average percentage error by numerical methods with respect to the exact method

Table 4. Ranking of various numerical methods on the basis of average percentage error with the exact method.

 Numerical Methods Average Percentage Error (%) Rank Finite Difference Method 0.07803 1 Rayleigh Ritz Method 0.34850 2 Galerkin Method 0.34850 2 Petrov Galerkin Method 0.66057 3 Collocation Method 1.13488 4 Subdomain Method 1.38199 5 Least Square Method 1.51389 6

From table 4, the numerical methods on the basis of minimum average percentage error with respect to exact solution and accuracy in their solution is listed as below:

1. Finite Difference Method (FDM)

2. Rayleigh-Ritz Method and Galerkin Method

3. Petrov Galerkin Method

4. Collocation Method

5. Subdomain Method

6. Least Square Method

7. CONCLUSION

1. In this research work, an attempt has been made to solve second-order differential equations by various methods of approximations and the Finite Difference Method. It has been found that from table 4, the Finite Difference Method gives closest solutions to exact solutions. Hence, the Finite Difference Method can be concluded as the most accurate method for solving a differential equation.

2. From Table 1, It can be concluded that the Galerkin and Rayleigh-Ritz Method gives identical solutions. Hence, Galerkin and Rayleigh-Ritz methods are almost in agreement with each other.

3. From Table 1, it is observed that at x = 0.5 (i.e. centre point of the domain) we get the peak value for all numerical methods as the equation is parabolic.

4. Further from Table 1, at points (0.4 and 0.6), (0.3 and 0.7), (0.2 and 0.8) and (0.1 and 0.9), the identical solution

obtained in each respective method, but only in the Petrov Galerkin Method, this trend cant be observed.

ACKNOWLEDGEMENTS

The authors would like to thanks the editor and reviewers for giving their valuable time.

REFERENCES

1. N. Ahmad and S. Charan, A comparative study of numerical solutions of second order ordinary differential equations with boundary value problems by shooting method and finite difference method, International Journal of Statistics and Applied Mathematics., vol. 4, no. 1, (January, 2019), pp. 18-22.

2. A. A. Anulo, A. S. Kibret, G. G. Gofna, A. D. Negassa, Numerical Solution of Linear Second Order Ordinary Differential Equations with Mixed Boundary Conditions by Galerkin Method, Mathematics and Computer Science., vol. 2, no. 5, (January, 2017), pp. 66-78.

3. H. K. Versteeg and W. Malalasekera, An Introduction to Computational Fluid Dynamics The Finite Volume Method, Logman Scientific and Technical, England, 1995.

4. J. N. Reddy, An Introduction to the Finite Element Method, McGraw-Hill Education, New York, 1993.