Abstract
This paper presents a modified Galerkin method based on sinc basis functions to numerically solve nonlinear boundary value problems. The modifications allow for the accurate approximation of the solution with accurate derivatives at the endpoints. The algorithm is applied to well-known problems: Bratu and Thomas-Fermi problems. Numerical results demonstrate the clear advantage of the suggested modifications in obtaining accurate numerical solutions as well as accurate derivatives at the endpoints.
1. Introduction
In this paper, we consider the following class of nonlinear boundary value problems (NBVPs) of the following form: subject to where .
This class of NBVPs covers many interesting and important problems. Two of these important problems that we consider in this paper are as follows.
(a) Bratu equation:
(b) Thomas-Fermi equation:
Nonlinear boundary value problems are widely used in science and engineering to describe many engineering, physical, and chemical models. For instance, Bratu problem is used in a variety of applications such as the solid fuel ignition in thermal combustion theory, the Chandrasekhar model of the expansion of the universe, chemical reaction theory, and radiative heat transfer, nanotechnology. For more details see [1–9]. In addition, Thomas-Fermi equation is used to calculate the electrostatic potential in the Thomas-Fermi atom model, see [10–13].
During the past decades, a great deal of theoretical and numerical research works have been invested for the solution of nonlinear boundary value problems. Although, each nonlinear problem has its own properties, the numerical schemes used to approximate the solution were almost the same. Bratu problem has been treated numerically via several techniques such as weighted residual method [14], finite difference method [1, 6], Laplace transform decomposition method [7, 15], Adomian decomposition method [8, 9], homotopy analysis method [3], He’s variational method [16], differential transformation method [17], and B-spline method [2]. Thomas-Fermi problem has also been numerically studied using variational approach [18, 19], the -expansion method [20], the decomposition method [21, 22], and the homotopy analysis method [11, 23]. Recently, Amore et al. [24] have discussed the numerical solution of Thomas-Fermi problem using Padé-Hankel method, numerical integration, power series with Padé, and Hermite-Padé approximants and Chebyshev polynomials.
However, to the authors best knowledge, neither Bratu nor Thomas-Fermi problem has been investigated using a Galerkin method based on sinc basis functions. In addition, the many different numerical schemes mentioned above do not treat well the approximation of the derivative of the solution at the endpoints. In this work, we consider the Sinc-Galerkin method with suitable modifications to investigate the numerical solution of problems of the form (1) subject to (2). The suggested modifications allow for an accurate numerical solution with accurate derivative at the endpoints.
The Sinc-Galerkin method was developed back in 1993, see [25]. It has been used by many authors in the numerical treatment of several types of differential equations, see, for example, [26–28]. The success of the Sinc-Galerkin method led us to consider it for the numerical treatment of two important well-known problems, namely, Bratu and Thomas-Fermi problems.
The rest of the paper is organized as follows. A review of the definition and the main properties of sinc functions and the interpolation properties of the sinc basis functions are presented in Section 2. Section 3 presents the derivation of the sinc-Galerkin method for problem (1) with homogeneous boundary conditions. Modifications for the treatment of nonhomogeneous boundary condition to better approximate the derivative of the solution at the endpoints are presented at the end of Section 3. Section 4 presents the numerical results for Bratu and Thomas-Fermi problems. Finally in Section 5, we conclude with some remarks.
2. Sinc Basis Functions Properties and Quadrature Interpolations
A thorough study of sinc basis functions and their interpolation properties are presented in [27] (see also [25]). In this section, we present a brief overview of the basic definitions and the required interpolation quadratures required for our subsequent development.
2.1. The Sinc Basis Functions
The sinc function, , is defined by For , the translated sinc functions, , , are defined by The functions , in (6) are dilated and translated versions of the mother sinc function . Note that both and are defined over the whole real line. In order to make use of them in the approximation of solutions of boundary value problems (1), an appropriate one-to-one map between and is needed. To this end, we introduce the following conformal mapping: which transforms the eye-shaped domain, , in the -plane, onto the infinite strip in the -plane, , where These are shown in Figure 1.
The sinc basis functions used in our numerical approach are given by the composition of the translated sinc functions and the conformal map as follows: for . A plot of is shown in Figure 2, for .
The functions in (9) will serve as the basis functions in the numerical solution of (1) on the interval . The mesh size defines the uniform grid points, , over the real line. The grid points in are given by the inverse image of the equispaced grid points , , namely,
It is important to mention that . Moreover, the basis functions , as defined in (9), are not differentiable at the endpoints and . Therefore, if the boundary conditions in (2) are not homogeneous, that is, , and/or differentiability of the solution at the end points is required, a modification of the expansion basis functions is needed. To this end, we define new sinc basis function as It can be easily verified that are differentiable at the endpoints with
The next subsection overviews how sinc basis is used to interpolate (approximate) functions defined over an interval .
2.2. Sinc Function Interpolation and Quadratures
Sinc interpolation and quadratures have been developed for a special class of functions (defined below). A discussion of the properties of such functions in and their sinc interpolation can be found in [25, 27]. Here we present only the relevant results to our work. We begin by the definition of .
Definition 1. Let be a conformal mapping from to with inverse . Let . Then is the class of functions which is analytic in and satisfies where and is the boundary of .
The following two theorems, whose proofs can be found in [27] (see also [25]), provide interpolation formulas for functions in as well as quadrature rule for their integrals, using sinc basis functions . These formulas will prove to be essential in the development of a numerical scheme for the solution of problems (1).
Theorem 2. Let and be in . Then for all , we have where
In numerical approximations, the sum in (14) must be truncated. The following theorem, whose proof also can be found in [27], indicates the conditions under which the truncation error is of exponential order.
Theorem 3. Let , be a conformal map, and , and constants such that
where
Then
For the selections , where denotes the greatest integer less than or equal to , the exponential order of the sinc quadrature rule in (18) is .
3. The Sinc-Galerkin Approach
In this section, we present the derivation of the Sinc-Galerkin approach for the problem (1) with homogeneous boundary conditions and differentiability requirement at the endpoints. The treatment of nonhomogeneous boundary conditions and differentiability requirement at the endpoints is dealt with at the end of this section. The problem at hand is subject to An approximate solution is sought as a finite expansion using the sinc basis functions , namely, where the basis functions are defined by (9) and are expansion coefficients to be determined. Note that , since if and 0, otherwise, and . The index in is there to indicate the number of basis functions used. The coefficients in (21) are determined by orthogonalizing the residual with respect to the basis functions . This yields the following discrete system: where the inner product is taken to be where is an appropriate weight function. A complete discussion on the choice of the weight function of the form can be found in [25, 27]. For our second order problem, we choose .
Another approach to find the coefficients , , is to analyze Using the sinc quadrature rule in (18), we approximate each of the inner products in (24). First, we consider Integrating by parts the right-hand side of (25), we obtain Another integration by parts gives The first term on the right-hand side of (27) is zero because . Since are not differentiable at the end points and , the weight function is chosen so that is differentiable at and . For the mapping , a suitable weight function is . With as such and using the boundary conditions (20), , the second term on the right-hand side of (27) is zero and (27) reduces to Note that . Setting we obtain, by expanding the derivatives under the integral in (28), where Applying the sinc quadrature rule (18) to the right-hand side of (30), assuming the integrand belongs to , we obtain where , , are given by We note that , , are bounded by (see, [28])
The bounds in (36) are used in the next theorem whose proof resembles the proof of Theorem in [28], to show that the error in approximating (30) by (34) is of order for some positive constants .
Theorem 4. Let be the conformal mapping and make the selections for and as in Theorem 3. With the assumption that , and for some constant . Then there exists a constant independent of such that
The other two inner products in (24) are approximated in a similarity. Using (18) yields the following approximations:
Combining (34) and (39), substituting , we obtain, for , We introduce the following notations in order to write the system (40) in a matrix-vector form. For any function , we define the diagonal matrix, , as Also define the matrices , to be the given by Finally, denote by , , and the column vectors where the superscript denotes the transpose, and by the column vector of length , each of whose components is equal to as Then system (40) can be written in the following matrix-vector equation:
The nonlinear system (45) of equations in the unknown coefficients , can be solved using the multidimensional version of Newton’s method. Once the coefficients are solved, the approximate solution is given by (21).
Treatment of Boundary Conditions. The above system (45) was derived based on homogeneous boundary conditions and no differentiability requirement at the endpoints. However, depending on the boundary conditions special modifications to the sinc basis functions are needed.
For the problem of nonhomogeneous boundary conditions (2), , and differentiability requirement at the endpoints, we transform the problem into one with homogeneous boundary conditions by writing where and are some suitable differentiable endpoint basis functions satisfying and and are unknown parameters to be determined. The reason and are chosen as in (47)-(48) is that we want and to have contribution at the endpoint and no contribution at the other endpoint , and similarly and to have contribution at the endpoint and no contribution at the other endpoint . This way, and approximate and , respectively, and and approximate and .
Substituting (46) into (19), we find that satisfies From (46)–(48), satisfies the homogeneous boundary conditions as Moreover, we impose that .
Following the derivation presented at the beginning of this section, we seek an approximate solution to using the modified Sinc basis functions (see (11)) as We determine the coefficients , , and by orthogonalizing each term in (49) with respect to for . Note that we added two basis functions and to accommodate the extra unknowns and . The resulting system for the unknowns , , and is where and are as in (31)–(33) with replaced by . The matrices , are now rectangular of size , given by Once the coefficients , , and are solved, the approximate solution for is given by
4. Numerical Examples
In this section, we apply the Sinc-Galerkin method outlined in the previous section to solve numerically Bratu’s and Thomas-Fermi problems.
4.1. Bratu’s Problem
Consider Bratu’s boundary value problem given by subject to It is well known that for Bratu’s problem (56)-(57), there is a such that the problem has two solutions for , no solutions for , and a unique solution for , where , see [29–32].
For , the two exact solutions are given by where is the solution of . For each , there are two possible values for . For ease of reference, the corresponding values for selected values of are displayed in Table 1.
In applying the Sinc-Galerkin approach to Bratu’s problem, we are particularly interested in capturing the two solutions for . This will be accomplished by initializing Newton’s iterative scheme with two different initial solutions, thereby the iterations converge to the two different solutions.
In the numerical simulation for Bratu’s problem (56)-(57), we have used the following endpoint basis functions: which satisfy (47)-(48) for any constants and . The weight function used is . The resulting nonlinear system (52) has been solved numerically for various values of using Newton’s method with two different initializations. The two different solutions have been captured. The results of the numerical simulations are shown in Figures 3 and 4 and Tables 2 and 3.
(a)
(b)
(c)
(d)
(e)
(f)
It can be clearly seen from Figure 3 that as approaches the critical value , the obtained two numerical solutions approach the unique one. In Figure 4, the maximum value of the two solutions, that is, versus is displayed. The plot clearly shows the bifurcation point at . Table 2 displays the maximum absolute error between the exact, , and the numerical, , solutions for various values of , and Table 3 displays the absolute error in the derivative at the endpoints.
From these results, we see that the presented Galerkin method based on modified sinc basis functions with boundary treatment is successful in obtaining accurate solution to the problem. Moreover, different initialization of Newton’s method was successful in obtaining the two different solutions for . This is in contrast to other works where only one solution branch is obtained. In addition, the added endpoint basis functions and for boundary treatment were the key to obtaining an accurate numerical solution with accurate derivatives at the endpoints. This was needed because the usual sinc basis functions are not differentiable at the endpoints and the modified sinc basis functions and their derivatives vanish at the endpoints.
4.2. The Thomas-Fermi Problem
The Thomas-Fermi problem consists of the following nonhomogeneous boundary value problem: subject to
It is well known that the value is very important because it determines the energy of a neutral atom in the Thomas-Fermi model as follows: where is the unclear charge, see [33]. Therefore any numerical scheme should produce a numerical solution, , with accurate . According to Kobayashi et al. [33], the initial slope is .
Problem (60)-(61) is transformed into one with homogeneous boundary conditions by setting where is an additional endpoint basis function with unknown parameter satisfying for any . Note that we only need one additional endpoint basis function at . Then satisfies
We approximate using the modified sinc basis functions as Following the general derivations outlined in the previous section, we obtain the following nonlinear system for the unknown coefficients and . where and , , and are as before.
In our numerical simulations we adopt the following conformal mapping, which maps onto . We use a weight function . For the additional endpoint basis function chose which satisfies (64) for any . Note that one can choose other forms for as long as it satisfies (64).
A plot of the approximate solution, , is displayed in Figure 5. The approximate along with the relative absolute error when compared with the value given in [33], , for various values of is displayed in Table 4.
We can see from Table 4 that is accurately approximated. The percentage error decreases as increases. These results are comparable to those obtained by [11] using the homotopy analysis method.
5. Conclusion
In this paper, the Sinc-Galerkin method has been applied to accurately solve nonlinear second-order boundary value problems (1)-(2). The Sinc-Galerkin approach using the basic Sinc basis functions is suitable only for homogeneous boundary conditions. Also, because the basis functions are not differentiable at the end points, the numerical solution does poor approximation of the derivative at the endpoints. In this work, we have modified the basic Sinc basis functions to new ones, , which are differentiable at the endpoints. Then additional endpoint basis functions, and , with free parameters, and , have been added in order to transform nonhomogeneous boundary value problem to a homogeneous one. Moreover, the added endpoint basis functions were there to accurately approximate the solution and its derivative at the endpoints.
We have applied the presented method to two well-known problems, Bratu and Thomas-Fermi problems. For Bratu problem, we were able to capture the two solutions for with accurate derivative at the endpoints. The two solutions were captured by initializing Newton’s method used to solve the resulting nonlinear system with two different initializations. For Thomas-Fermi problem the numerical solution obtained accurately approximates the initial slope .
The results obtained using the modified Sinc-Galerkin method with added endpoint basis functions proved to be very accurate and lands itself as an excellent alternative numerical scheme for numerically solving nonlinear boundary value problems.