Finite Difference Method for Boundary Value Problems in Ordinary Differential Equations

In this paper, a second order numerical method based on finite difference method with uniform mesh is presented for solving boundary value problems. Numerical solution is found for the boundary value problem using finite difference method and the results are tabulated and compared with analytical solution. Graphs are also depicted in support of the numerical results. The present method is simple and it approximates the exact solution very well.


Introduction
Ordinary and partial differential equation have long played important roles in the history of modelling physical phenomena; they continue to serve as indispensable tools in future investigations. In the theory and application of ordinary and partial differential equations, one may be interested to find a solution to a differential equation satisfying certain defined conditions. If the conditions are given at only one point of the independent variable, we have initial conditions; where as if the conditions are given at more than one point of the independent variable, they are called boundary conditions (BC). The problem of finding a solution of an ℎ order differential equation together with n-initial conditions is called an initial value problem (IVB), however if the n-boundary conditions are considered, the problem is called boundary value problem (BVP).
Boundary value problems for linear second order equations are particularly important because of numerous applications in science and technology. In this paper as in most physical applications, boundary conditions are always imposed at end points of an interval.
While there are many numerical methods for solving such boundary value problems, the method of finite difference is most commonly used. The basis of the finite difference method is the replacement of all derivatives occurring by the corresponding difference quotients; this is applicable to any problem in differential equations and the resulting linear system of equation is solved by any standard procedure. These roots are the value of the required solution at the pivotal points.
Some researchers have studied and numerically solved second order boundary-value problems using different method with different boundary conditions, for instance, the collocation method [2], finite difference method [9], shooting method [1], quantic spline [14], non-polynomial spline [6] and references therein. Hence, the purpose of this paper is to develop a numerical method for solving second order boundary-value problem (1).
This paper is arranged in the following manner. In section 2, some definitions and notions are presented. In section 3, we formulate and described the proposed method. We have discussed convergence of the proposed method under appropriate condition in Section 4. The Numerical examples and applications of linear boundary value problems of the second order are discussed in Section 5 and Section 6, respectively. Finally, conclusion drawn is presented.
finite difference in second order linear boundary value problem of ordinary differential equation.
Definition 2.1. An equation containing the derivative of one or more dependent variables with respect to one or more independent variables is said to be differential equation.

Definition 2.2.
A differential equation is said to be an ordinary differential equation if it contains only ordinary derivatives of one or more dependent variables with respect to a single independent variable.
3. An ℎ order ordinary differential equation is said to be linear in if it can be written in the form Where ( ) are continuous function 2.1. Linear boundary value problem The differential equation and the boundary condition are linear (these are called BVP) and for which the boundary conditions are based on only two points say 1 = and 2 = where < . The differential equation may be written in the form ( ) ( ) + −1 ( ) ( −1) + ⋯ + 0 ( ) = ( )

With boundary conditions
( 1 ) = , ( 2 ) = The differential equation is homogeneous when ( ) = 0 in [ , ] otherwise non-homogeneous. The boundary value problem is called homogeneous when the differential and all the boundary condition are homogeneous. That is with boundary conditions ( 1 ) = ( 2 ) = .

Linear system of equations
A number of problems in numerical analysis can be reduced to system of linear equations. Let us consider a system of real linear equations, which may be written in matrix form = Where, = ( ) is the × matrix of coefficients, = ( ) the column matrix of the unknowns and, = ( ) the column matrix of the right-hand sides. With inequality holding for at least one row , then the Jacobi iterations converge for each , when applied to the equation = .

The Finite Difference Method for Linear BVP in ODEs
The solution of BVP by finite difference method is accomplished by the following steps.
1. Discretizing the continuous solution domain into a discrete finite difference grid.
2. Approximating the exact derivatives in the Ordinary Differential Equation by finite difference approximation.
3. Substituting finite difference approximation into the Ordinary Differential Equation to obtain algebraic finite difference equation.
4. Solving the resulting system of equations by any standard procedure.

Description of the Method
We first consider the case of a linear two-point boundary value problem of the form Where we assume that ( ) ≥ 0 and ( ) ≥ 0. In order to develop the numerical method for finding the solution of differential equation (6), the interval [ , ] is divided into equal subintervals using the mesh points = + ℎ, = 0, 1, … , , 0 = , = , and ℎ = ( − )/ .
Consider the special case when ( ) is constant. Applying Taylor series expansions on ( − ℎ) and ( + ℎ) respectively and adding the two, we obtain second order central difference If we apply the above equation (7) at the point = and insert it into the differential equation (6), we get We then define an approximation to ( ) by omitting the (ℎ 2 ) error terms. We let , = 1, 2, … , − 1 satisfy: Since the boundary conditions are given, we set 0 = and = . Hence, we have − 1 linear equations for the − 1 unknowns 1 , … , −1 .
For later use, it is useful to define the local truncation error of the method from Eqs. (8) and (9), given by This is basically how well the true solution satisfies the discrete equations. Using the Taylor expansions defined above, we see that We next consider the approximate discrete system in more detail. Multiplying our discrete equation (9) by ℎ 2 and collecting terms, we get − +1 + [2 + ℎ 2 ( )] − −1 = ℎ 2 ( ) Setting = 2 + ℎ 2 ( ), we may write this system of linear equations in matrix form as: Thus, we need to solve the linear system = , where is a symmetric tri-diagonal matrix. Note that if min [ , ] ( ) ≥ * > 0, then is strictly diagonally dominant and hence invertible.
Hence the Jacobi iteration of linear system of equation converges and by theorem (2.2) converges for each , when applied to = .
Now let us return to the case when depends on . One possibility is to write and approximate the first term at = (as above) by To get an approximation to ′ that is also second order, we could use the Taylor series to write In that case, our discrete equation (10) would become There are two drawbacks to this approach. First, the matrix we get will not be symmetric, and thus will involve more work to solve. To see this, let = 2 ( ) + ℎ 2 ( ), = ( ) − ℎ ′ ( ) 2 , = ( ) − ℎ ′ ( ) 2 Then, we may write this system of linear equations in matrix form as: When is not constant, we do not expect that +1 = . Another drawback to this approach is that we would also have to compute ′.
Instead, we consider another way to approximate the derivatives in the equation. Using Taylor series expansions, we get either  Inserting these quantities, we obtain: Although the middle terms seems to only (ℎ), we get by applying the Mean Value Theorem, that And so Denoting ( ± ℎ 2 ) by ±1/2 , and omitting the (ℎ 2 ) terms, we are led to the approximation: with homogeneous boundary conditions ( ) = ( ) = 0 has a unique solution ∈ 2 [ , ].
Now for the numerical solution of the differential equation (11) we choose an equidistance grid = + ℎ, = 1, 2, … , + 1 with the step size ℎ = ( − )/ + 1 and ∈ at the internal grid points , = 1, 2, … , , we For approximate values to the exact solution ( ). Here we have set = ( ) and = ( ). The system has to be complemented by the two boundary conditions For an abbreviated notation we introduce the × tri-diagonal matrix. Let And the vectors = ( 1 , … , ) and = ( 1 , … , ) . Then our system of equations including the boundary conditions = This tri-diagonal linear system is non-singular and symmetric can be easily solved for from any right-hand side of .

Stability and Convergence Analysis
The error and convergence analysis are initiated by first establishing the following two lemmas.

Numerical Examples
To validate the applicability of the method, we have considered the following model problems with exact solution.    Table 2.  An energy balance on a differential element of the rod yields the following second order boundary value problem of ordinary differential equation Where is the temperature of the rod (c) and 2 = ℎ where ℎ is heat transfer coefficient (J/s-m 2 -k), is perimeter of the rod ( ), is the thermal conductivity, is surface area of the rod ( 2 ) and is temperature of surrounding (℃).

Conclusion
The accuracy of finite difference method for second order linear Boundary value problem in ODEs depends on the size of the sub-interval ℎ and also on the order of approximation. As we reduce ℎ, the accuracy improves but the number of equations to be solved also increases. When a large number of subdivisions are used, we do not treat a difference equation as a set of simultaneous equations. As a result, the amount of computation time required may become excessive, and good accuracy may be difficult to achieve. The use of higher-order approximations will yield greater accuracy for the same mesh size but results in considerable complication, especially near the end points of the interval. In practice, it is advisable to solve the linear system for several different values of ℎ. A comparison of the solutions at the same mesh points indicate the more accuracy being obtained.