Variable Mesh Size Exponential Finite Difference Method for the Numerical Solutions of Two Point Boundary Value Problems

In this article, we presented a variable mesh size exponential finite difference scheme for the numerical solutions of two point boundary value problems with Dirichlet’s boundary conditions. Under appropriate condition, we have discussed the local truncation error and the convergence of the proposed method. Numerical experiments demonstrated the use and computational efficiency of the method in several model problems. Numerical results showed that the proposed method is convergent and has at least second order of accuracy which is in good agreement with the theoretically established order of the method.


Introduction
In this article we considered a method for the numerical solution of the twopoint boundary value problems of the form Two point boundary value problems are of common occurrence in many areas of sciences and engineerings.This class of problems has gained importance in the literature for the variety of their applications.In most cases it is impossible to obtain solutions of these problems using analytical methods which satisfy the given boundary conditions.In these cases we resort to approximate solution of the problems and the last few decades have seen substantial progress in the development of approximate solutions of these problems.In the literature, there are many different methods and approaches such as method of integration and discretization that are used to derive the approximate solutions of these problems [1,2,3,4].
The existence and uniqueness of the solution to problem (1) is assumed.We further assumed that problem (1) is well posed with continuous derivatives and that the solution depends differentially on the boundary conditions.The specific assumption on f (x, y) to ensure existence and uniqueness will not be considered [3,4,5].
Over the last few decades, finite difference methods [6,7,8] have generated renewed interest and in recent years, variety of specialized techniques [10,11] for the numerical solution of boundary value problems in ODEs have been reported in the literature.Recently, an exponential finite difference method with uniform step size was proposed in [12] for the numerical solution of linear two point boundary value problem.This method generated impressive numerical results for the problem (1).Hence, the purpose of this article is to propose an exponential finite difference method with variable step size for problem (1).The development of this numerical method for two-point boundary-value problems plays a paramount role in the approximate solution of boundary value problems with a small parameter affecting the highest derivative of the differential equation.Boundary value problems with such property are known as singularly perturbed two-point boundary-value problems.It is a well known fact that singularly perturbed boundary value problem possess a small interval in which the solution varies rapidly and this small interval is known as the boundary layer in the literature.The occurrence of boundarylayers creates difficulty for most standard numerical schemes with uniform mesh size in solving these problems.A variable mesh method can overcome this difficulty and is well suit for solving boundary layer problem [8,9].Our proposed variable step size exponential method for the solution of two point boundary value problems is efficient in solving such boundary layer problems without any difficulty.
We hope that others may find the proposed method an improvement and appealing to those existing finite difference methods for two-point boundary value problems.
A new method of at least quadratic order is proposed for the numerical solution of linear boundary value in [12] problems (1).Our idea is to apply the exponen-tial finite difference method to discretize equation (1) in order to get a system of algebraic equations.In addition, if we apply a linearization technique, the method results in a tridiagonal matrix for the nodal values.The elements of this matrix depend on the source function i.e. right-hand side of the ordinary differential equation as well as on its partial derivatives with respect to the dependent variable and its first-order derivative.To the best of our knowledge, no similar method for the numerical solution of problem (1) has been discussed in the literature so far.
We have presented our work in this article as follows.In the next section we derived a new variable mesh size exponential finite difference method.In Section 3, local truncation error and convergence of the new method are discussed in Section 4. The application of the proposed method to the problems in (1) has been presented and illustrative numerical results have been produced to show the efficiency of the new method in Section 5. Discussion and conclusion on the performance of the new method are presented in Section 6.

The Variable Mesh Size Exponential Difference Method
We define N finite numbers of nodal points of the domain [a,b], in which the solution of the problem (1) is desired, as a ≤ x 0 < x 1 < x 2 < ...... < x N < x N +1 = b, using nonuniform step length h such that x i+1 = x i + h i+1 , i = 0, 1, 2, ....., N and r i = hi+1 hi .Suppose that we wish to determine the numerical approximation of the theoretical solution y(x) of the problem (1) at the nodal point x i , i = 1, 2, ....., N .We denote the numerical approximation of y(x) at node x = x i as y i .Let us denote f i as the approximation of the theoretical value of the source function f (x, y(x)) at node x = x i , i = 0, 1, 2, ....., N + 1.We can define other notations used in this article i.e. f i±1 , and y i±1 , in the similar way.Following the ideas in [11,12], we propose an approximation to the theoretical solution y(x i ) of the problem (1) by the exponential finite difference scheme as, where a 0 , a 1 , a 2 and b 0 are unknown functions of the argument r i and φ(x i ) is an unknown sufficiently differentiable function of x.Let us define a function F i (h, y) and associate it with (2) as, Assume that φ(x i ) can be expanded in Taylor series about the point x = x i−1 .Hence we write φ(x i ) in Taylor series , The application of (4) in the expansion of exp(φ(x i )) will provide an O(h 2 i ) approximation of the form as, Expand F i (h, y) in Taylor series about mesh point x = x i and using (5) in it, we have On comparing the coefficients of h p i , p = 0, 1, 2, 3 both sides in ( 6), we get the following system of nonlinear equations To determine the unknown functions a 0 , a 1 , b 0 , φ(x i−1 ) and φ ′ (x i−1 ) in ( 7), we have to assign arbitrary value to some unknown functions.To simplify the system of equations in (7), we have considered the following assumption: Using ( 8) in (7) and solved the reduced system of equations, we obtained Write f ′ i for y (3) i in (9) and substituting the values of φ(x i−1 ) and φ ′ (x i−1 ) from ( 8) and ( 9) in ( 4), we have Finally substitute the values of a 0 , a 1 , b 0 and φ(x i ) from ( 9) and ( 10) in ( 2), we obtain our proposed exponential finite difference method as For each nodal point, we will obtain the nonlinear system of equations given by ( 11) or a linear system of equations if the source function is f(x).In the derived numerical method (11), the exponential function exp( . If f i in the denominator of the argumnt becomes zero in the domain of the solution, we take the series expansion of the function exp( ) and neglecting the second and higher order terms.Therefore method (11) becomes For the computational purpose in Section 4, we have used the following second order finite difference approximation in place of h i f ′ i in (11) and in ( 12):

Local Truncation Error
We can write the following expression for the term in (11) with the help of (13): Write the expansion for the exponential function in the ( 14) by neglecting the third and higher order terms, so we will obtain, (15) From ( 11) and ( 15), the truncation error T i at the nodal point x = x i may be written as [8,13,14], By the Taylor series expansion of y at nodal point x = x i and using i and etc., we have (16) can be simplified as follows : Thus we have obtained a truncation error at each node of O(h 4 i ).

Convergence of the Method
Let us substitute ( 15) into (11) and then simplify (11), we have Thus where Let us define The difference method (18) represents a system of nonlinear equations in unknown y i , i = 1, 2, ..., N .Let us write (18) in matrix form as, where is a tridiagonal matrix.Let Y be the exact solution of (18), so it will satisfy the matrix equation where Y is a column matrix of order N × 1 which can be obtained by replacing y with Y in matrix y and T is a truncation error matrix in which each element has O(h 4 i ).
Variable Mesh Size Exponential Finite Difference Method

15
Let us define , After linearization of f i+1 , we have where G i+1 = ( ∂f ∂Y ) i+1 .Thus Similarly, we can linearize f i−1 ,and f i , to obtain the following results : By taking the Taylor series expansion of G i±1 about x = x i , and from the difference of ( 19) and ( 20), we can write where P = (P lm ) N ×N is a tri-diagonal matrix defined as and , where Let us assume that the solution of difference equation (11) has no roundoff error.So from ( 18), ( 19) and (20) we have Let us define G 0 = {G i : i = 1, 2, ..., N }, We further define H 0 = {( ∂G ∂x ) i , i = 1, 2, ....., N }.Let there exist some positive constant W such that t 0 ≤ W, ∀ t 0 ∈ H 0 .So it is possible for very small h i , ∀i = 1, 2, ..., N , 1×N , denotes the row sum of the matrix J = (J lm ) N ×N where On neglecting the higher order terms i.e.O(h 3 i ) in R i then it is easy to see that J is irreducible [13].By the row sum criterion and for sufficiently small h i , ∀i = 1, 2, ..., N , J is monotone [15].Thus J −1 exist and J −1 ≥ 0. For the bound of J, we define [16,17] where We note that higher order terms i.e.O(h 3 i ) in the above expressions are neglected.Let d l (J) ≥ 0, ∀ l and Thus from ( 25) and (26), we have It follows from ( 16) and ( 27) that E → 0 as h i → 0. Thus we conclude that method (11) converges and the order of the convergence of method ( 11) is at least quadratic.

Numerical Results
To illustrate our method and demonstrate its computational efficiency, we considered some model problems.In each model problem, we took non uniform step size h i .In Table 1 -Table 8, we have shown the maximum absolute error (MAY), computed for different values of N and is defined as The starting value of the step length h 1 is calculated by formula where r = r i , ∀ i = 1, 2, ...., N in computation.In case of uniform mesh r = 1 , the above formula for computation of step length becomes h = b−a N .The order of the convergence (O N ) of the method ( 11) is estimated by the formula where m can be estimated by considering the ratio of N ′ s.
We have used Newton-Raphson iteration method to solve the system of nonlinear equations arised from equation (23).All computations were performed on a MS Window 2007 professional operating system in the GNU FORTRAN environment version 99 compiler (2.95 of gcc) on Intel Duo Core 2.20 Ghz PC.The solutions are computed on N nodes and iteration is continued until either the maximum difference between two successive iterates is less than 10 (−10) or the number of iteration reached 10 3 .
Problem 1.The first model problem is a linear problem [18] given by The analytical solution is y( . The MAY computed by method (11) for different values of N and ǫ are presented in Table 1.The analytical solution is y(x) = ) 2 .The MAY computed by method ( 11) for different values of N are presented in Table 2. Problem 3. The third model problem is a linear problem [19] given by where f (x) is calculated so that y is the analytical solution.The MAY computed by method (11)for different values of N and ǫ are presented in Table 3 .Problem 4. The fourth model problem is a linear problem [20] given by where f (x) is calculated so that y( ) is analytical solution.The MAY computed by method (11) for different values of N and ǫ are presented in Table 4 and Table 5. Problem 5.The fifth model problem is a linear problem [21] given by where f (x) is calculated so that y is the analytical solution.The MAY computed by method (11) for different values of N and ǫ are presented in Table 6, Table 7 and Table 8.

25
We have described a new method for numerically solving two-point boundary value problems and several model problems considered to demonstrate the performance of the proposed method.Numerical result for examples 1 which is presented in table 1, for different values of r i show as ǫ decreases N increases, with uniform step size maximum absolute errors in our method increases.On the other hand for the variable mesh, maximum absolute error decreases with increase in N. The numerical results for examples 2 and 4 are accurate in both for uniform and non-uniform mesh sizes.The results for examples 3 and 5 are same as for example 1 with uniform mesh size.The maximum absolute error decreases with increase in N and decrease in ǫ except for N = 100 in tables 3 & 6 with non-uniform mesh size.For the comparison purpose, A numerical results by existing finite difference method in [22] is considered for uniform mesh in examples 4 and 5.This comparison show that our method has less discretization error and perform better for considered model problems.Note that for small N, method [22] yield good results in model problems 4 and 5.However, as the N becomes larger, the exponential finite difference method shows less error than the method [22].Over all method (11) is convergent and convergence of the method depends on choice of mesh ratio r i .

Conclusion
A new method to find the numerical solution of two point boundary value problems has been developed.The new method has advantages and disadvantages when considered individually.For example, at each nodal point x = x i , i = 1, 2, ....., N, we will obtain a system of algebraic equations given by (11).If the source function is f(x) then the system of equations from (11) is linear otherwise we will obtain nonlinear system of equations, which is always difficult to be solved, disadvantage of the new method.The decision to use a certain difference scheme depends on its computational efficiency for accurate solution.It is obvious that special method required for some special problem where the solution is not regular and varies rapidly.But on the other hand, the new method produces good numerical approximate solutions for variety of model problems without any modification either in method or in problem and its rate of convergence is quadratic.It may be an advantage of the new method.The numerical results of the model problems showed that the new method is computationally efficient and plays an important role to obtain accurate numerical solutions.The idea presented in this article leads to the possibility to develop difference methods to solve third order and fourth order boundary value problems in ordinary differential equations.Works in these direction are in progress.

10 P
subject to the boundary conditions y(a) = α and y(b) = β, where α and β are real constants and f is continuous on (x, y) for all x ∈ [a, b] y ∈ ℜ. .K. Pandey and B. D. Pandey