# Modified Sinc-Galerkin Method for Nonlinear Boundary Value Problems

Authors: M. A. Hajji; Q. M. Al-Mdallal
Journal: Journal of Mathematics
 (2013)

---

## 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.


---

## Body

## 1. Introduction

In this paper, we consider the following class of nonlinear boundary value problems (NBVPs) of the following form:(1)y′′(x)+λf(x,y(x))=g(x),x∈(a,b),
subject to
(2)y(a)=α0,y(b)=α1,
where b∈ℝ.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:(3)y′′+λey=0,0<x<1,y(0)=0,y(1)=0.(b) Thomas-Fermi equation:(4)y′′-y3/2x=0,0<x<∞,y(0)=1,y(∞)=0.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 Section2. 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,sinc(x), is defined by
(5)sinc(x)={sin(πx)πx,x≠0,1,x=0.
For h>0, the translated sinc functions, S(k,h)(x), k∈Z, are defined by
(6)S(k,h)(x)≡sinc(x-khh)={sin[(π/h)(x-kh)][(π/h)(x-kh)],x≠kh,1,x=kh.
The functions S(k,h)(x),k∈ℤ, in (6) are dilated and translated versions of the mother sinc function sinc(x). Note that both sinc(x) and S(k,h)(x) 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 (a,b) and (-∞,∞) is needed. To this end, we introduce the following conformal mapping:
(7)w=ϕ(z)=ln(z-ab-z)
which transforms the eye-shaped domain, 𝒟E, in the z-plane, onto the infinite strip in the w-plane, 𝒟S, where
(8)𝒟E={z=x+iy:|arg(z-ab-z)|<d≤π2},𝒟S={w=u+iv:|v|<d≤π2}.
These are shown in Figure 1.Figure 1
Eye-shaped domain and the infinite strip.The sinc basis functions used in our numerical approach are given by the composition of the translated sinc functionsS(k,h) and the conformal map ϕ as follows:
(9)Sk(z)=S(k,h)∘ϕ(z)=sinc(ϕ(z)-khh),
for z∈𝒟E. A plot of Sk(x) is shown in Figure 2, for (a,b)=(0,1).Figure 2
The graph ofSk(x)=S(k,h)∘ϕ(x) for k=-2,0,2,h=π/12 and b=1.The functionsSk(z) in (9) will serve as the basis functions in the numerical solution of (1) on the interval (a,b). The mesh size h defines the uniform grid points, {kh},k∈Z, over the real line. The grid points {xk} in 𝒟E=(a,b) are given by the inverse image of the equispaced grid points {kh}, k∈Z, namely,
(10)xk=ϕ-1(kh)=a+bekh1+ekh.It is important to mention thatlimx→aSk(x)=limx→bSk(x)=0. Moreover, the basis functions Sk,k∈Z, as defined in (9), are not differentiable at the endpoints x=a and x=b. Therefore, if the boundary conditions in (2) are not homogeneous, that is, α0α1≠0, 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
(11)ξk(z)=Sk(z)ϕ′(z).
It can be easily verified that ξk(z) are differentiable at the endpoints with
(12)limz→aξk(z)=limz→bξk(z)=limz→aξk′(z)=limz→bξk′(z)=0.The next subsection overviews how sinc basisSk is used to interpolate (approximate) functions defined over an interval (a,b).


### 2.2. Sinc Function Interpolation and Quadratures

Sinc interpolation and quadratures have been developed for a special class of functionsB(𝒟E) (defined below). A discussion of the properties of such functions in B(𝒟E) 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 B(𝒟E).Definition 1.
Letϕ:𝒟E→𝒟S be a conformal mapping from 𝒟E to 𝒟S with inverse ψ≡ϕ-1. Let Γ={ψ(u)∈𝒟E:-∞<u<∞}. Then B(𝒟E) is the class of functions f which is analytic in B(𝒟E) and satisfies
(13)(1)∫ψ(t+L)‍|f(z)|dz→0,t→±∞,(2)∫∂𝒟E‍|f(z)|dz<∞,
where L={iv:|v|<d≤π/2} and ∂𝒟E is the boundary of 𝒟E.The following two theorems, whose proofs can be found in [27] (see also [25]), provide interpolation formulas for functions in B(𝒟E) as well as quadrature rule for their integrals, using sinc basis functions Sk. These formulas will prove to be essential in the development of a numerical scheme for the solution of problems (1).Theorem 2.
Letϕ′ and f be in B(𝒟E). Then for all z∈Γ, we have(14)f(z)=∑j=-∞∞‍f(zj)Sj(z)+Ef,
where
(15)Ef=sin(πϕ(z)/h)2πi∫∂𝒟E‍ϕ′(w)f(w)dw(ϕ(w)-ϕ(z))sin(πϕ(w)/h).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.
Letf∈B(𝒟E), ϕ be a conformal map, and α,β, and C constants such that
(16)|F(z)ϕ′(z)|≤C{exp(-α|ϕ(z)|),z∈Γ-,exp(-β|ϕ(z)|),z∈Γ+,
where
(17)Γ-={z∈Γ:ϕ(z)=x∈(-∞,0)}=(a,(a+b)2)Γ+={z∈Γ:ϕ(z)=x∈[0,-∞)}=((a+b)2,b).
Then
(18)∫0b‍F(z)dz=h∑j=-MN‍F(zj)ϕ′(zj)+𝒪(e-αMh)+𝒪(e-βNh)+𝒪(e-2πd/h).
For the selectionsN=⌊|(α/β)M+1|⌋,h=2πd/(αM), where ⌊x⌋ denotes the greatest integer less than or equal to x, the exponential order of the sinc quadrature rule in (18) is 𝒪(exp(-2πdαM)).


## 2.1. The Sinc Basis Functions

The sinc function,sinc(x), is defined by
(5)sinc(x)={sin(πx)πx,x≠0,1,x=0.
For h>0, the translated sinc functions, S(k,h)(x), k∈Z, are defined by
(6)S(k,h)(x)≡sinc(x-khh)={sin[(π/h)(x-kh)][(π/h)(x-kh)],x≠kh,1,x=kh.
The functions S(k,h)(x),k∈ℤ, in (6) are dilated and translated versions of the mother sinc function sinc(x). Note that both sinc(x) and S(k,h)(x) 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 (a,b) and (-∞,∞) is needed. To this end, we introduce the following conformal mapping:
(7)w=ϕ(z)=ln(z-ab-z)
which transforms the eye-shaped domain, 𝒟E, in the z-plane, onto the infinite strip in the w-plane, 𝒟S, where
(8)𝒟E={z=x+iy:|arg(z-ab-z)|<d≤π2},𝒟S={w=u+iv:|v|<d≤π2}.
These are shown in Figure 1.Figure 1
Eye-shaped domain and the infinite strip.The sinc basis functions used in our numerical approach are given by the composition of the translated sinc functionsS(k,h) and the conformal map ϕ as follows:
(9)Sk(z)=S(k,h)∘ϕ(z)=sinc(ϕ(z)-khh),
for z∈𝒟E. A plot of Sk(x) is shown in Figure 2, for (a,b)=(0,1).Figure 2
The graph ofSk(x)=S(k,h)∘ϕ(x) for k=-2,0,2,h=π/12 and b=1.The functionsSk(z) in (9) will serve as the basis functions in the numerical solution of (1) on the interval (a,b). The mesh size h defines the uniform grid points, {kh},k∈Z, over the real line. The grid points {xk} in 𝒟E=(a,b) are given by the inverse image of the equispaced grid points {kh}, k∈Z, namely,
(10)xk=ϕ-1(kh)=a+bekh1+ekh.It is important to mention thatlimx→aSk(x)=limx→bSk(x)=0. Moreover, the basis functions Sk,k∈Z, as defined in (9), are not differentiable at the endpoints x=a and x=b. Therefore, if the boundary conditions in (2) are not homogeneous, that is, α0α1≠0, 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
(11)ξk(z)=Sk(z)ϕ′(z).
It can be easily verified that ξk(z) are differentiable at the endpoints with
(12)limz→aξk(z)=limz→bξk(z)=limz→aξk′(z)=limz→bξk′(z)=0.The next subsection overviews how sinc basisSk is used to interpolate (approximate) functions defined over an interval (a,b).


## 2.2. Sinc Function Interpolation and Quadratures

Sinc interpolation and quadratures have been developed for a special class of functionsB(𝒟E) (defined below). A discussion of the properties of such functions in B(𝒟E) 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 B(𝒟E).Definition 1.
Letϕ:𝒟E→𝒟S be a conformal mapping from 𝒟E to 𝒟S with inverse ψ≡ϕ-1. Let Γ={ψ(u)∈𝒟E:-∞<u<∞}. Then B(𝒟E) is the class of functions f which is analytic in B(𝒟E) and satisfies
(13)(1)∫ψ(t+L)‍|f(z)|dz→0,t→±∞,(2)∫∂𝒟E‍|f(z)|dz<∞,
where L={iv:|v|<d≤π/2} and ∂𝒟E is the boundary of 𝒟E.The following two theorems, whose proofs can be found in [27] (see also [25]), provide interpolation formulas for functions in B(𝒟E) as well as quadrature rule for their integrals, using sinc basis functions Sk. These formulas will prove to be essential in the development of a numerical scheme for the solution of problems (1).Theorem 2.
Letϕ′ and f be in B(𝒟E). Then for all z∈Γ, we have(14)f(z)=∑j=-∞∞‍f(zj)Sj(z)+Ef,
where
(15)Ef=sin(πϕ(z)/h)2πi∫∂𝒟E‍ϕ′(w)f(w)dw(ϕ(w)-ϕ(z))sin(πϕ(w)/h).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.
Letf∈B(𝒟E), ϕ be a conformal map, and α,β, and C constants such that
(16)|F(z)ϕ′(z)|≤C{exp(-α|ϕ(z)|),z∈Γ-,exp(-β|ϕ(z)|),z∈Γ+,
where
(17)Γ-={z∈Γ:ϕ(z)=x∈(-∞,0)}=(a,(a+b)2)Γ+={z∈Γ:ϕ(z)=x∈[0,-∞)}=((a+b)2,b).
Then
(18)∫0b‍F(z)dz=h∑j=-MN‍F(zj)ϕ′(zj)+𝒪(e-αMh)+𝒪(e-βNh)+𝒪(e-2πd/h).
For the selectionsN=⌊|(α/β)M+1|⌋,h=2πd/(αM), where ⌊x⌋ denotes the greatest integer less than or equal to x, the exponential order of the sinc quadrature rule in (18) is 𝒪(exp(-2πdαM)).


## 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 y(a)=y(b)=0 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
(19)𝒩(y)=y′′(x)+λf(x,y(x))=g(x),x∈(a,b),
subject to
(20)y(a)=0,y(b)=0.
An approximate solution is sought as a finite expansion using the sinc basis functions Sk, namely,
(21)y(x)≈ym(x)=∑k=-MN‍ckSk(x),m=M+N+1,
where the basis functions Sk(x) are defined by (9) and ck are expansion coefficients to be determined. Note that ym(xj)≈cj, since Sk(xj)=δkj=1 if k=j and 0, otherwise, and xj=ϕ-1(jh). The index m in ym is there to indicate the number of basis functions used. The coefficients ck in (21) are determined by orthogonalizing the residual 𝒩(ym)-g with respect to the basis functions {Sj}j=-MN. This yields the following discrete system:
(22)〈𝒩(um)-g,Sj〉=0,-M≤j≤N,
where the inner product is taken to be
(23)〈h(x),g(x)〉=∫ab‍h(x)g(x)w(x)dx,
where w(x) is an appropriate weight function. A complete discussion on the choice of the weight function of the form w(x)=1/(ϕ′(x))n can be found in [25, 27]. For our second order problem, we choose w(x)=1/ϕ′(x).Another approach to find the coefficientsyk, -M≤k≤N, is to analyze
(24)〈y′′,Sj〉+λ〈f(x,y(x)),Sj〉-〈g(x),Sj〉=0,-M≤j≤N.
Using the sinc quadrature rule in (18), we approximate each of the inner products in (24). First, we consider
(25)〈y′′(x),Sj〉=∫ab‍y′′(x)Sj(x)w(x)dx.
Integrating by parts the right-hand side of (25), we obtain
(26)〈y′′(x),Sj〉=y′(x)Sj(x)w(x)|ab-∫ab‍y′(x)ddx[Sj(x)w(x)]dx.
Another integration by parts gives
(27)〈y′′(x),Sj〉=y′(x)Sj(x)w(x)|ab-y(x)ddx[Sj(x)w(x)]|ab+∫ab‍y(x)d2dx2[Sj(x)w(x)]dx.
The first term on the right-hand side of (27) is zero because Sj(a)=Sj(b)=0. Since Sj are not differentiable at the end points a=a and x=b, the weight function w(x) is chosen so that (d/dx)[Sj(x)w(x)] is differentiable at x=a and x=b. For the mapping ϕ(x)=ln((x-a)/(b-x)), a suitable weight function is w(x)=1/ϕ′(x). With w(x) as such and using the boundary conditions (20), y(a)=y(b)=0, the second term on the right-hand side of (27) is zero and (27) reduces to
(28)〈y′′(x),Sj〉=∫ab‍y(x)d2dx2[Sj(x)w(x)]dx.
Note that (d/dx)[Sj(x)]=(d/dϕ)[Sj]ϕ′(x). Setting
(29)dndϕn[Sj(x)]=Sj(n)(x),
we obtain, by expanding the derivatives under the integral in (28),
(30)〈y′′(x),Sj〉=∫ab‍y(x)×[Sj(x)g20(x)+Sj(1)(x)g21(x)+Sj(2)(x)g22(x)]dx,
where
(31)g20(x)=w′′(x),(32)g21(x)=ϕ′′(x)w(x)+2ϕ′(x)w′(x),(33)g22(x)=(ϕ′(x))2w(x).
Applying the sinc quadrature rule (18) to the right-hand side of (30), assuming the integrand belongs to B(𝒟E), we obtain
(34)〈y′′(x),Sj〉≈h∑k=-MN‍y(xk)ϕ′(xk)×[1h2δjk(2)g22(xk)+1hδjk(1)g21(xk)+δjk(0)g20(xk)],
where δjk(i), i=0,1,2, are given by
(35)δjk(0)=ddϕ[S(k,h)∘ϕ(x)]|xk={1,j=k,0,j≠k,δjk(1)=hddϕ[S(k,h)∘ϕ(x)]|xk={0,j=k,(-1)k-j(k-j),j≠k,δjk(2)=h2d2dϕ2[S(j,h)∘ϕ(x)]|xk={-π23h2,j=k,-2(-1)k-j(k-j)2,j≠k.
We note that |δkj(i)|, i=0,1,2, are bounded by (see, [28])
(36)|δkj(i)|≤1,i=0,1,|δkj(2)|≤π23.The bounds in (36) are used in the next theorem whose proof resembles the proof of Theorem  3.1 in [28], to show that the error in approximating (30) by (34) is of order 𝒪(exp(-c/h)) for some positive constants c,h.Theorem 4.
Letϕ be the conformal mapping and make the selections for h and N as in Theorem 3. With the assumption that y(x)g2i(x)∈B(𝒟E),i=0,1,2, and
(37)|y(z)g2i(z)ϕ′(z)|≤C1{exp(-α|ϕ(z)|),z∈Γ-,exp(-β|ϕ(z)|),z∈Γ+,
for some constant C1. Then there exists a constant C2 independent of M such that
(38)|〈y′′(x),Sj〉-h∑k=-MN‍∑i=02‍y(xk)ϕ′(xk)hiδjk(i)g2i(xk)|≤C2M2exp(-πdαM).The other two inner products in (24) are approximated in a similarity. Using (18) yields the following approximations:
(39)〈f(x,y(x)),Sj〉≈h∑k=-MN‍δjk(0)f(xk,y(xk))w(xk)ϕ′(xk)=hf(xj,y(xj))w(xj)ϕ′(xj),〈g(x),Sj〉≈h∑k=-MN‍δjk(0)g(xk)w(xk)ϕ′(xk)=hg(xj)w(xj)ϕ′(xj).Combining (34) and (39), substituting ck=y(xk), we obtain, for -M≤j≤N,
(40)∑k=-MN‍∑i=02‍ckϕ′(xk)hiδjk(i)g2i(xk)+λf(xj,cj)w(xj)ϕ′(xj)-g(xj)w(xj)ϕ′(xj)=0.
We introduce the following notations in order to write the system (40) in a matrix-vector form. For any function g(x), we define the m×m diagonal matrix, m=M+N+1, as
(41)𝒟(g)=diag[g(x-M),…,g(x0),…,g(xN)].
Also define the matrices I(i),i=0,1,2, to be the m×m given by
(42)I(i)=[δjk(i)],j,k=-M,…,N.
Finally, denote by x→, y→, and f→(x,c) the column vectors
(43)x→=[x-M,x-M+1,…,xN-1,xN]t,c→=[c-M,c-M+1,…,cN-1,cN]t,f→(x,c)=[f(x-M,c-M),f(x-M+1,c-M+1),…,f(xN,cN)]t,
where the superscript t denotes the transpose, and by 1→ the column vector of length m, each of whose components is equal to 1 as
(44)1→=[1,1,…,1]t.
Then system (40) can be written in the following matrix-vector equation:
(45)[∑i=02‍1hiI(i)𝒟(g2iϕ′)]c→+λ𝒟(wϕ′)f→(x,c)=𝒟(gwϕ′)1→.The nonlinear system (45) of m=M+N+1 equations in the m unknown coefficients ck,-M≤k≤N, can be solved using the multidimensional version of Newton’s method. Once the coefficients ck are solved, the approximate solution ym(x) 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), y(a)=α0,y(b)=α1, and differentiability requirement at the endpoints, we transform the problem into one with homogeneous boundary conditions by writing
(46)z(x)=y(x)-𝒫a(x;ca)-𝒫b(x;cb),
where 𝒫a(x;ca) and 𝒫b(x;cb) are some suitable differentiable endpoint basis functions satisfying
(47)𝒫a(a;ca)=α0,𝒫a(b;ca)=0,𝒫a′(b;ca)=0,(48)𝒫b(a;ca)=0,𝒫b(b;ca)=α1,𝒫b′(a;cb)=0,
and ca and cb are unknown parameters to be determined. The reason 𝒫a(x) and 𝒫b(x) are chosen as in (47)-(48) is that we want 𝒫a(x) and 𝒫a′(x) to have contribution at the endpoint x=a and no contribution at the other endpoint x=b, and similarly 𝒫b(x) and 𝒫b′(x) to have contribution at the endpoint x=b and no contribution at the other endpoint x=a. This way, 𝒫a(a) and 𝒫a′(a) approximate y(a) and y′(a), respectively, and 𝒫b(b) and 𝒫b′(b) approximate y(b) and y′(b).Substituting (46) into (19), we find that z(x) satisfies
(49)z′′+λf(x,𝒫a(x;ca)+z(x)+𝒫b(x;cb))=g(x)-𝒫a′′(x;ca)-𝒫b′′(x;cb).
From (46)–(48), z(x) satisfies the homogeneous boundary conditions as
(50)z(a)=z(b)=0.
Moreover, we impose that z′(a)=z′(b)=0.Following the derivation presented at the beginning of this section, we seek an approximate solution toz(x) using the modified Sinc basis functions ξk(x) (see (11)) as
(51)z(x)=∑k=-MN‍ckξk(x).
We determine the coefficients ck,-M≤k≤N, ca, and cb by orthogonalizing each term in (49) with respect to ξj(x) for -M-1≤j≤N+1. Note that we added two basis functions ξ-M-1(x) and ξN+1(x) to accommodate the extra unknowns ca and cb. The resulting system for the unknowns ck,-M≤k≤N, ca, and cb is
(52)[∑i=02‍1hiI~(i)𝒟(g~2i(ϕ′)2)]c→+λ𝒟(w~ϕ′)f~→(x,c)=𝒟(g~w~ϕ′)1→,
where
(53)f~→(x,c)=f→(x,cϕ′+𝒫a(x;ca)+𝒫b(x;cb)),g~=g-𝒫a′′(x;ca)-𝒫b′′(x;cb),w~(x)=w(x)ϕ′(x),
and   g~2i are as in (31)–(33) with w replaced by w~. The matrices I~(i),i=0,1,2, are now rectangular of size (m+2)×m, given by
(54)I~(i)=[δjk(i)],-M-1≤j≤N+1,-M≤k≤N.Once the coefficients ck,-M≤k≤N, ca, and cb are solved, the approximate solution for y(x) is given by
(55)y(x)≈𝒫a(x;ca)+∑k=-MN‍ckξk(x)+𝒫b(x;cb).


## 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(56)y′′(x)+λey=0.x∈(0,1),
subject to
(57)y(0)=0,y(1)=0.
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 λ*=3.513830719, see [29–32].Forλ<λ*, the two exact solutions are given by
(58)yex(x)=-2ln[cosh[(x-1/2)(θ/2)]cosh(θ/2)],
where θ is the solution of θ=2λcosh(θ/4). For each λ<λ*=3.513830719, there are two possible values for θ. For ease of reference, the corresponding θ values for selected values of λ are displayed in Table 1.Table 1
The two valuesθ1 and θ2 for selected values of λ.
λ
0.5
1.0
1.5
2.0
2.5
3.0
3.5
3.513830719
θ
1
1.03357
1.51716
1.93973
2.35755
2.81155
3.37351
4.55185
4.79871
θ
2
13.03824
10.93870
9.58170
8.50720
7.54810
6.57657
5.05434
4.79871In 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:
(59)𝒫a(x;ca)=cax(1-x)2,𝒫b(x;cb)=cbx2(1-x)
which satisfy (47)-(48) for any constants ca and cb. The weight function used is w(x)=1. 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.Table 2
max
0
≤
x
≤
1
|
y
ex
(
x
)
-
y
ap
(
x
)
| for selected values of λ.
λ
Lower solution
Upper solution
0.1
4.07348
×
1
0
-
14
7.02362
×
1
0
-
8
0.5
1.1485
×
1
0
-
12
1.25449
×
1
0
-
9
1.0
5.39141
×
1
0
-
12
9.6844
×
1
0
-
10
1.5
1.44806
×
1
0
-
11
8.11187
×
1
0
-
10
2.0
3.37709
×
1
0
-
11
6.79725
×
1
0
-
10
2.5
7.40526
×
1
0
-
11
6.93944
×
1
0
-
10
3.0
1.55144
×
1
0
-
10
8.17894
×
1
0
-
10
3.5
2.46664
×
1
0
-
9
3.10425
×
1
0
-
9Table 3
|
y
ap
′
(
0
)
-
y
ex
′
(
0
)
| and |yap′(1)-yex′(1)| for selected values of λ.
λ
Lower solution
Upper solution
|
(
y
ap
′
(
0
)
-
y
ex
′
(
0
)
|
|
y
ap
′
(
1
)
-
y
ex
′
(
1
)
|
|
y
ap
′
(
0
)
-
y
ex
′
(
0
)
|
|
y
ap
′
(
1
)
-
y
ex
′
(
1
)
|
0.1
8.18169
×
1
0
-
9
8.18169
×
1
0
-
9
5.35475
×
1
0
-
4
5.35475
×
1
0
-
4
0.5
2.1811
×
1
0
-
7
2.1811
×
1
0
-
7
2.48943
×
1
0
-
4
2.48943
×
1
0
-
4
1.0
9.54514
×
1
0
-
7
9.54514
×
1
0
-
7
2.00222
×
1
0
-
4
2.00222
×
1
0
-
4
1.5
2.38345
×
1
0
-
6
2.38345
×
1
0
-
6
1.67767
×
1
0
-
4
1.67767
×
1
0
-
4
2.0
4.80072
×
1
0
-
6
4.80072
×
1
0
-
6
1.40583
×
1
0
-
4
1.40583
×
1
0
-
4
2.5
8.78773
×
1
0
-
6
8.78773
×
1
0
-
6
1.15255
×
1
0
-
4
1.15255
×
1
0
-
4
3.0
1.5852
×
1
0
-
5
1.5852
×
1
0
-
5
8.90305
×
1
0
-
4
8.90305
×
1
0
-
4
3.5
3.7798
×
1
0
-
5
3.7798
×
1
0
-
5
4.94364
×
1
0
-
4
4.94364
×
1
0
-
4Numerical solutions of Bratu’s problem (56)-(57) for λ=1.0, 2.0, 3.0, 3.5, 3.513, 3.5138, with M=N=80,α=1/2,d=π/12, and h=πd/(αN).
(a)
(b)
(c)
(d)
(e)
(f)Figure 4
A plot ofmax0≤x≤1|y(x)| versus λ for 0≤λ≤λ*.It can be clearly seen from Figure3 that as λ approaches the critical value λ*=3.513830719, the obtained two numerical solutions approach the unique one. In Figure 4, the maximum value of the two solutions, that is, max0≤x≤1|y(x)| versus 0≤λ≤λ* is displayed. The plot clearly shows the bifurcation point at λ*. Table 2 displays the maximum absolute error between the exact, yex(x), and the numerical, yap(x), 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 𝒫a(x) and 𝒫b(x) 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 Sk(x) are not differentiable at the endpoints and the modified sinc basis functions ξk(x) and their derivatives ξk′(x) vanish at the endpoints.


### 4.2. The Thomas-Fermi Problem

The Thomas-Fermi problem consists of the following nonhomogeneous boundary value problem:(60)y′′(x)-y3x=0,x∈(0,∞),
subject to
(61)y(0)=1,y(∞)=0.It is well known that the valuey′(0) is very important because it determines the energy of a neutral atom in the Thomas-Fermi model as follows:
(62)E=67(4π3)2/3Z7/3y′(0),
where Z is the unclear charge, see [33]. Therefore any numerical scheme should produce a numerical solution, yap(x), with accurate yap′(0). According to Kobayashi et al. [33], the initial slope y′(0) is y′(0)=-1.588071.Problem (60)-(61) is transformed into one with homogeneous boundary conditions by setting
(63)z(x)=y(x)-𝒫a(x;ca),
where 𝒫a(x;ca) is an additional endpoint basis function with unknown parameter ca satisfying
(64)𝒫a(0;ca)=1,𝒫a(∞;ca)=0.
for any ca. Note that we only need one additional endpoint basis function at x=0. Then z(x) satisfies
(65)z′′-(z+𝒫a(x;ca))3/2x=-𝒫a′′(x;ca),z(0)=0,z(∞)=0.We approximatez(x) using the modified sinc basis functions as
(66)z(x)≈∑k=-MN‍ckξk(x).
Following the general derivations outlined in the previous section, we obtain the following nonlinear system for the unknown coefficients ca and ck,-M≤k≤N.
(67)[∑i=02‍1hiI~(i)𝒟(g~2i(ϕ′)2)]c→-𝒟(w~ϕ′)f~→(x,c)=𝒟(g~w~ϕ′)1→,
where
(68)f~→(x,c)=f→(x,cϕ′+𝒫a(x;ca)),f(x,y)=y3/2x,g~=-𝒫a′′(x;ca),
and w~, g~2i, and   I~(i) are as before.In our numerical simulations we adopt the following conformal mapping,(69)ϕ(x)=ln(sinh(x)),
which maps (0,∞) onto (-∞,∞). We use a weight function w(x)=1. For the additional endpoint basis function 𝒫a(x;ca) chose
(70)𝒫a(x;ca)=(cax+1)e-x,
which satisfies (64) for any ca. Note that one can choose other forms for 𝒫a(x;ca) as long as it satisfies (64).A plot of the approximate solution,yap(x), is displayed in Figure 5. The approximate y′(0)≈yap′(0) along with the relative absolute error when compared with the value given in [33], y′(0)=-1.588071, for various values of N is displayed in Table 4.Table 4
Approximations of the initial slopey'(0) and the corresponding errors to the result in [33].
N  (=M)
y
ap
'
(
0
)
Error (%)
30
-
1.53333
3.44718
40
-
1.55499
2.08332
50
-
1.56659
1.35237
60
-
1.57311
0.941921
70
-
1.57678
0.711127
80
-
1.57867
0.591687
90
-
1.57934
0.549815Figure 5
A plot of the numerical solution of the Thomas-Fermi problem (56)-(57), using M=N=80, α=1/2, and d=π/12.We can see from Table4 that y′(0) is accurately approximated. The percentage error decreases as N increases. These results are comparable to those obtained by [11] using the homotopy analysis method.


## 4.1. Bratu’s Problem

Consider Bratu’s boundary value problem given by(56)y′′(x)+λey=0.x∈(0,1),
subject to
(57)y(0)=0,y(1)=0.
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 λ*=3.513830719, see [29–32].Forλ<λ*, the two exact solutions are given by
(58)yex(x)=-2ln[cosh[(x-1/2)(θ/2)]cosh(θ/2)],
where θ is the solution of θ=2λcosh(θ/4). For each λ<λ*=3.513830719, there are two possible values for θ. For ease of reference, the corresponding θ values for selected values of λ are displayed in Table 1.Table 1
The two valuesθ1 and θ2 for selected values of λ.
λ
0.5
1.0
1.5
2.0
2.5
3.0
3.5
3.513830719
θ
1
1.03357
1.51716
1.93973
2.35755
2.81155
3.37351
4.55185
4.79871
θ
2
13.03824
10.93870
9.58170
8.50720
7.54810
6.57657
5.05434
4.79871In 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:
(59)𝒫a(x;ca)=cax(1-x)2,𝒫b(x;cb)=cbx2(1-x)
which satisfy (47)-(48) for any constants ca and cb. The weight function used is w(x)=1. 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.Table 2
max
0
≤
x
≤
1
|
y
ex
(
x
)
-
y
ap
(
x
)
| for selected values of λ.
λ
Lower solution
Upper solution
0.1
4.07348
×
1
0
-
14
7.02362
×
1
0
-
8
0.5
1.1485
×
1
0
-
12
1.25449
×
1
0
-
9
1.0
5.39141
×
1
0
-
12
9.6844
×
1
0
-
10
1.5
1.44806
×
1
0
-
11
8.11187
×
1
0
-
10
2.0
3.37709
×
1
0
-
11
6.79725
×
1
0
-
10
2.5
7.40526
×
1
0
-
11
6.93944
×
1
0
-
10
3.0
1.55144
×
1
0
-
10
8.17894
×
1
0
-
10
3.5
2.46664
×
1
0
-
9
3.10425
×
1
0
-
9Table 3
|
y
ap
′
(
0
)
-
y
ex
′
(
0
)
| and |yap′(1)-yex′(1)| for selected values of λ.
λ
Lower solution
Upper solution
|
(
y
ap
′
(
0
)
-
y
ex
′
(
0
)
|
|
y
ap
′
(
1
)
-
y
ex
′
(
1
)
|
|
y
ap
′
(
0
)
-
y
ex
′
(
0
)
|
|
y
ap
′
(
1
)
-
y
ex
′
(
1
)
|
0.1
8.18169
×
1
0
-
9
8.18169
×
1
0
-
9
5.35475
×
1
0
-
4
5.35475
×
1
0
-
4
0.5
2.1811
×
1
0
-
7
2.1811
×
1
0
-
7
2.48943
×
1
0
-
4
2.48943
×
1
0
-
4
1.0
9.54514
×
1
0
-
7
9.54514
×
1
0
-
7
2.00222
×
1
0
-
4
2.00222
×
1
0
-
4
1.5
2.38345
×
1
0
-
6
2.38345
×
1
0
-
6
1.67767
×
1
0
-
4
1.67767
×
1
0
-
4
2.0
4.80072
×
1
0
-
6
4.80072
×
1
0
-
6
1.40583
×
1
0
-
4
1.40583
×
1
0
-
4
2.5
8.78773
×
1
0
-
6
8.78773
×
1
0
-
6
1.15255
×
1
0
-
4
1.15255
×
1
0
-
4
3.0
1.5852
×
1
0
-
5
1.5852
×
1
0
-
5
8.90305
×
1
0
-
4
8.90305
×
1
0
-
4
3.5
3.7798
×
1
0
-
5
3.7798
×
1
0
-
5
4.94364
×
1
0
-
4
4.94364
×
1
0
-
4Numerical solutions of Bratu’s problem (56)-(57) for λ=1.0, 2.0, 3.0, 3.5, 3.513, 3.5138, with M=N=80,α=1/2,d=π/12, and h=πd/(αN).
(a)
(b)
(c)
(d)
(e)
(f)Figure 4
A plot ofmax0≤x≤1|y(x)| versus λ for 0≤λ≤λ*.It can be clearly seen from Figure3 that as λ approaches the critical value λ*=3.513830719, the obtained two numerical solutions approach the unique one. In Figure 4, the maximum value of the two solutions, that is, max0≤x≤1|y(x)| versus 0≤λ≤λ* is displayed. The plot clearly shows the bifurcation point at λ*. Table 2 displays the maximum absolute error between the exact, yex(x), and the numerical, yap(x), 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 𝒫a(x) and 𝒫b(x) 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 Sk(x) are not differentiable at the endpoints and the modified sinc basis functions ξk(x) and their derivatives ξk′(x) vanish at the endpoints.


## 4.2. The Thomas-Fermi Problem

The Thomas-Fermi problem consists of the following nonhomogeneous boundary value problem:(60)y′′(x)-y3x=0,x∈(0,∞),
subject to
(61)y(0)=1,y(∞)=0.It is well known that the valuey′(0) is very important because it determines the energy of a neutral atom in the Thomas-Fermi model as follows:
(62)E=67(4π3)2/3Z7/3y′(0),
where Z is the unclear charge, see [33]. Therefore any numerical scheme should produce a numerical solution, yap(x), with accurate yap′(0). According to Kobayashi et al. [33], the initial slope y′(0) is y′(0)=-1.588071.Problem (60)-(61) is transformed into one with homogeneous boundary conditions by setting
(63)z(x)=y(x)-𝒫a(x;ca),
where 𝒫a(x;ca) is an additional endpoint basis function with unknown parameter ca satisfying
(64)𝒫a(0;ca)=1,𝒫a(∞;ca)=0.
for any ca. Note that we only need one additional endpoint basis function at x=0. Then z(x) satisfies
(65)z′′-(z+𝒫a(x;ca))3/2x=-𝒫a′′(x;ca),z(0)=0,z(∞)=0.We approximatez(x) using the modified sinc basis functions as
(66)z(x)≈∑k=-MN‍ckξk(x).
Following the general derivations outlined in the previous section, we obtain the following nonlinear system for the unknown coefficients ca and ck,-M≤k≤N.
(67)[∑i=02‍1hiI~(i)𝒟(g~2i(ϕ′)2)]c→-𝒟(w~ϕ′)f~→(x,c)=𝒟(g~w~ϕ′)1→,
where
(68)f~→(x,c)=f→(x,cϕ′+𝒫a(x;ca)),f(x,y)=y3/2x,g~=-𝒫a′′(x;ca),
and w~, g~2i, and   I~(i) are as before.In our numerical simulations we adopt the following conformal mapping,(69)ϕ(x)=ln(sinh(x)),
which maps (0,∞) onto (-∞,∞). We use a weight function w(x)=1. For the additional endpoint basis function 𝒫a(x;ca) chose
(70)𝒫a(x;ca)=(cax+1)e-x,
which satisfies (64) for any ca. Note that one can choose other forms for 𝒫a(x;ca) as long as it satisfies (64).A plot of the approximate solution,yap(x), is displayed in Figure 5. The approximate y′(0)≈yap′(0) along with the relative absolute error when compared with the value given in [33], y′(0)=-1.588071, for various values of N is displayed in Table 4.Table 4
Approximations of the initial slopey'(0) and the corresponding errors to the result in [33].
N  (=M)
y
ap
'
(
0
)
Error (%)
30
-
1.53333
3.44718
40
-
1.55499
2.08332
50
-
1.56659
1.35237
60
-
1.57311
0.941921
70
-
1.57678
0.711127
80
-
1.57867
0.591687
90
-
1.57934
0.549815Figure 5
A plot of the numerical solution of the Thomas-Fermi problem (56)-(57), using M=N=80, α=1/2, and d=π/12.We can see from Table4 that y′(0) is accurately approximated. The percentage error decreases as N 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 Sk is suitable only for homogeneous boundary conditions. Also, because the basis functions Sk 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, ξk, which are differentiable at the endpoints. Then additional endpoint basis functions, wa(x;ca) and wb(x;cb), with free parameters, ca and cb, 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 y′(0).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.



---

*Source: 875930-2013-01-08.xml*