Solving Nonlinear Two Point Boundary Value Problems Using Exponential Finite Difference Method

In this article, we presented an exponential finite difference scheme for solving nonlinear 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 demonstrate the use and computational efficiency of the method. Numerical results show that this method is at least fourth order accurate, which is good agreement with the theoretically established order of the method.


Introduction
Two point boundary value problems for ordinary differential equations arise in many branches of sciences and engineering.The existence of the solutions of the two point boundary value problems, either associated with system of linear or nonlinear ordinary differential equations and boundary conditions are specified at two points of the domain, depends on the domain considered for the solution of the problems.In most case it is impossible to obtain solutions of these problems using analytical methods which satisfy the given specified boundary conditions.In these cases we resort to approximate solutions and the last few decades have seen substantial progress in the development of approximate solutions of these problems.

P.K.Pandey
In the literature, there are many different methods and approaches such as method of integration and discretization which be used to derive the approximate solutions in the domain of these problems [1,2,3,4].
In this article we proposed a method for the numerical solution of the boundary value problems of the form subject to the boundary conditions y(a) = α and y(b) = β, where α and β are real constants and f is continuous on (x, y, y ′ ) for all x ∈ [a, b] y, y ′ ∈ ℜ.
The existence and uniqueness of the solution to problem (1) is assumed.Further we 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, y ′ ) to ensure existence and uniqueness will not be considered [3,4,5].
Over the last few decades, high order finite difference method [6,7,8] have generated renewed interest and in recent years, variety of specialized techniques [9,10] for the numerical solution of boundary value problems in ODEs have been reported in the literature.Recently, an exponential finite difference method was proposed in [11] for the numerical solution of linear two point boundary value problem.This method generated impressive numerical results for the linear problem in (1).Hence, the purpose of this article is to propose an exponential finite difference method to nonlinear problems (1), in the hope that others may find the proposed method an improvement to those existing finite difference method.
A new method of at least order four is an extension of the method which was developed for the numerical solution of linear problems based on local assumption.Our idea is to apply the exponential finite difference method to discretize equation (1) in order to get a nonlinear system of algebraic equations.To the best of our knowledge, no similar method for the numerical solution of problem (1) has been discussed in literature so far.
In the next section we discussed our exponential finite difference method.In section 3, we derived our method ; local truncation error and convergence of the method are discussed in Section 4 & 5 respectively.The application of the developed method to the problems (1) has been presented and illustrative numerical results have been produced to show the efficiency of the new method in Section 6. Discussion and conclusion on the performance of the method are presented in Section 7.

The Exponential Difference Method
We defined N + 1 finite numbers of nodal points of the domain [a,b], in which the solution of the problem (1) is desired, as x i = a + ih, i = 0, 1, 2, ........, N using uniform step length where h = b−a N , x 0 = a and x N = b.Suppose we wish to determine numerical approximation of the theoretical solution y(x) of the problem (1) at the nodal point x i , i = 1, 2, ....., N − 1 and denote as y i .Let f i denotes the approximation of the theoretical value of the source function f (x, y(x), y ′ (x)) at node x = x i , i = 0, 1, 2, ....., N .We can define other notations f i±1 , y i±1 , in the similar way used in this article.To develop the exponential difference method for the numerical solution of the problem (1), we need the following definitions: Define and We note that c and d from equations ( 7) and ( 8) respectively, are finite parameters to be determined.We proposed the exponential difference method for solving problem (1) numerically as,

Derivation of the Method
By the Taylor series expansion about node x = x i , from (3) we have Let us define G 1 i±1 = ( ∂f ∂y ′ ) i±1 , so from (5)we have

P.K.Pandey
Similarly from ( 4) and ( 6), we have By the Taylor series expansion of G 1 i±1 about node x = x i and from ( 13) and ( 14), we have On expanding (1) in Taylor series about x = x i , then substitute in (7) together with ( 15), we have y ′ i will provide fourth order approximation for y ′ if we choose parameter c in ( 16) Thus from ( 16) and ( 17) we have find y ′ i , a fourth order approximation for y ′ i i.e So from ( 9) and (18), we have Let us define Using the approximations defined above, we can prove that f i+1 + f i−1 − 2 f i will provide a fourth order approximation for Finally, following the idea in [11] for the source function f (x, y), from (11), we proposed our fourth order exponential difference method for solving problem (1) numerically as, For each nodal point x = x i , i = 1, 2, ...., N − 1, we will obtain a system of nonlinear equations given by (23).

Local Truncation Error
from equations ( 19),( 20) and ( 22), by Taylor series expansion of f on each node x = x i ,we have From ( 23) and ( 24), the truncation error T i at the nodal point x = x i may be written as [8,12,13], By the Taylor series expansion of y at nodal point x = x i and second order expansion of exponential function, we have

Convergence of the Method
Let us write the second order expansion of the exponential function in (23) and then simplify, we have Let us define Let us define column matrix S N ×1 as where [.....] t is transpose of column matrix.The difference method (26) represents a system of nonlinear equations in unknown y i , i = 1, 2, ..., N .Let us replace above defined S by φ &y, so we can write (26) in matrix form as, where is tridiagonal matrix.Let Y be the exact solution of (26) and replace above defined S by Y &T , so we can write (26) in matrix form as where T is truncation error matrix in which each element has O(h 6 ).Let us define , After linearization of f i+1 , we have where G i+1 = ( ∂f ∂Y ) i+1 and H i+1 = ( ∂f ∂Y ′ ) i+1 .Thus Similarly, we can linearize f i−1 , f i , f i and obtained the following results : where By Taylor series expansion of G i±1 & H i±1 about x = x i , and from (29)-(32), we can write where P = (P lm ) N ×N is a tridiagonal matrix defined as , where E i = (y i − Y i ), i = 1, 2, ...., N .Let assume that the solution of difference equation (26) has no roundoff error.So from ( 27),( 28) and (33) we have Let us define We further define H 0 = {H i , H 1 i , H 2 i , ( ∂H ∂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, 1×N , denotes the row sum of the matrix J = (J lm ) N ×N where It is easy to see that J is irreducible [12].By row sum criterion and for sufficiently small h, J is monotone [14].Thus J −1 exist and J −1 ≥ 0. For the bound of J, we define [15,16] where Let d l (J) ≥ 0, ∀ l and d * (J) = min 1≤l≤N d l (J).
Thus from ( 34) and ( 35), we have It follows from ( 25) and ( 36) that E → 0 as h → 0. Thus we conclude that method (23) converges and the order of convergence is at least four.

Numerical Results
To illustrate our method and demonstrate its computationally efficiency, we consider some model problems.In each case, we took uniform step size h.
In Table 1 -Table 5, we have shown the maximum absolute error (MAY), computed for different values of N and is defined as 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 -1 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 nonlinear problem given by The analytical solution is y(x) = log( 1.0 1+x ).For comparison purpose, we computed the MAY by the method in [17].The MAY computed by both methods for different values of N are presented in Table 1.Problem 2. The second model problem is a nonlinear problem The analytical solution is y(x) = 1 1+x .For comparison purpose, we computed the MAY by the method in [17].The MAY computed by both methods for different values of N are presented in Table 2. Problem 3. The third model problem is a nonlinear problem given by The analytical solution is y(x) = 1 √ 1+x .For comparison purpose, we computed the MAY by the method in [17].The MAY computed by both methods for different values of N are presented in Table 3. Problem 4. The fourth model problem is a nonlinear problem given by where f (x) is calculated so that y(x) = 1 − (x 2 + 1) 2 is analytical solution.For comparison purpose, we also computed the MAY by the method in [17].The MAY computed by both methods for different values of N are presented in Table 4.The analytical solution is y(x) = exp(x).Solving this model problem by method in [17], for each nodal point we obtained a system of linear equations.We applied Gauss-Seidel iterative for the solution of resulting system of linear equations.For comparison purpose, we computed the MAY by the method in [17].The MAY computed by both methods for different values of N are presented in Table 5.   1-4, show that the exponential finite difference method has less discretization error and is definitely better than the finite difference method [17].Table 5 shows the accuracy and efficiency of the exponential finite difference method for solving linear problem numerically.Note that for small N, method [17] yield good results except in model problems 2 and 4.However, as the N becomes larger, the exponential finite difference method shows less error than the method [17].It is an advantage of the exponential finite difference method over existing method [17].

Conclusion
A new approach to obtain the numerical solution of general two point boundary value problems has been developed.The new scheme 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 nonlinear equation given by (23), which is always difficult to be solved.On the other hand, the new method has a high order of convergence which yield smaller discretization error.The decision to use a certain difference scheme does not only depend on the given order of the method but also on its computational efficiency.The numerical results of model problems showed that the new method is computationally efficient.It is also observed from the results that method has high accuracy i.e. small discretization error.In the present article high order finite difference method has been derived based on exponential function.This new method leads to the possibility to develop new difference methods to solve third order and fourth order boundary value problems in ordinary differential equations.Works in this direction is in progress.

Table 5 :
Maximum absolute error in y(x) = exp(x) for Problem 5.