# A Note on Fourth Order Method for Doubly Singular Boundary Value Problems

Authors: R. K. Pandey; G. K. Gupta
Journal: Advances in Numerical Analysis
 (2012)

---

## Abstract

We present a fourth order finite difference method for doubly singular boundary value problem(p(x)y′(x))′=q(x)f(x,y),0<x≤1 with boundary conditions y(0)=A(ory'(0)=0orlimx→0p(x)y'(x)=0) and αy(1)+βy'(1)=γ, where α(>0), β(≥0), γ and A are finite constants. Here p(0)=0 and q(x) is allowed to be discontinuous at the singular point x=0. The method is based on uniform mesh. The accuracy of the method is established under quite general conditions and also corroborated through one numerical example.


---

## Body

## 1. Introduction

Consider the following class of singular two point boundary value problems:(1.1)(p(x)y′)′=q(x)f(x,y),0<x≤1,
with boundary conditions
(1.2)y(0)=A,αy(1)+βy′(1)=γ(1.3)or,y′(0)=0(orlimx→0p(x)y′(x)=0),αy(1)+βy′(1)=γ,
where α>0,β≥0 and A,γ are finite constants. Here p(0)=0 and q(x) is unbounded near x=0. The condition p(0)=0 says that the problem is singular and as q(x) is unbounded near x=0, so the problem is doubly singular [1]. Let p(x), q(x), and f(x,y) satisfy the following conditions:(A-1)
(i)
p(0)=0,p(x)>0 in (0,1],(ii)
p(x)∈C[0,1]∩C1(0,1],(iii)
p(x)=xb0g(x),(0≤b0<1),g(x)>0, and(iv)
G∈C4[0,1] and Gv(x) exists on (0,1),whereG(x)=1/g(x) on [0,1].(A-2)
(i)
q(x)>0 in (0,1],q(x) is unbounded near x=0,(ii)
q(x)∈L1(0,1),(iii)
q(x)=xa0H(x),H(0)>0 with (a0>-1),(iv)
H(x)∈C4[0,1] and Hv(x) exists on (0,1), and(v)
1+a0-b0≥0.(A-3)  f(x,y) is continuous on {[0,1]×ℝ},∂f/∂y exists, and is continuous and∂f/∂y≥0 for all 0≤x≤1 and for all real y.Thomas [2] and Fermi [3] independently derived a boundary value problem for determining the electrical potential in an atom. The analysis leads to the nonlinear singular second order problem
(1.4)y′′=x-1/2y3/2
with a set of boundary conditions. The following are of our interest:(i)
the neutral atom with Bohr radiusb given by y(0)=1,by′(b)-y(b)=0;(ii)
the ionized atom given byy(0)=1,y(b)=0.Furthermore, Chan and Hon [4] have considered the generalized Thomas-Fermi equation
(1.5)(xby′)′=cq(x)ye,y(0)=1,y(a)=0
for parameter values 0≤b<1,c>0,d>-2,e>1, and q(x)=xb+d.Such singular problems have been the concern of several researchers [5–8]. The existence-uniqueness of the solution of the boundary value problem (1.1) with boundary condition (1.2) or (1.3) is established in [1, 9–13]. Bobisud [1] has mentioned that in case limx→0+(q(x)/p′(x))≠0, the condition y′(0)=0 is quite severe, that is, it is sufficient but not necessary for forcing the solution to be differentiable at x=0. In fact if limx→0+p(x)y′(x)=0, then
(1.6)limx→0+y′(x)=limx→0+q(x)p′(x)·limx→0+f(x,y(x)).
Thus limx→0+y′(x)=0 if either limx→0+(q(x)/p′(x))=0 or f(0,y(0))=0. But limx→0+(q(x)/p′(x))≠0 if q(x) has discontinuity at x=0; thus it is natural to consider the weaker boundary condition limx→0+p(x)y′(x)=0.There is a considerable literature on numerical methods forq(x)=1 but to the best of our knowledge very few numerical methods are available to tackle doubly singular boundary value problems. Reddien [14] has considered the linear form of (1.1) and derived numerical methods for q(x)∈L2[0,1] which is stronger assumption than (A-2)(ii).Some second order methods (Chawla and Katti [15], Pandey and Singh [16, 17]) as well as fourth order methods (Chawla et al. [18–21], Pandey and Singh [22–24]) have been developed for q(x)=1,p(x). Most of the researchers have developed methods for the function p(x)=xb0,0≤b0<1 and the boundary conditions y(0)=A and y(1)=B.Chawla [19] has given fourth order method for the problem (1.1)-(1.2) with q(x)=1 and p(x)=xb0,0≤b0<1. The method is extended by Pandey and Singh [22] to a class of function p(x)=xb0g(x) satisfying (A-1)(i–iii) with G(x)=1/g(x) analytic in the neighborhood of the singular point.In this work we extend the fourth order accuracy method developed in [22] to doubly singular boundary value problems (1.1)-(1.2) and (1.1), (1.3) both. That is, we allow the data function to be discontinuous at the singular point x=0. Further, we do not require analyticity of the function G(x)=1/g(x). The convergence of the method is established under quite general conditions on the functions p(x),q(x), and f(x,y). Fourth order convergence of the method is corroborated through one example and maximum absolute errors are displayed in Table 1.Table 1
Maximum absolute errors and order of the methods for Example4.1.
N
b
0
=
0.1
Order
b
0
=
0.5
Order
b
0
=
0.1
Order
b
0
=
0.5
Order
y
(
0
)
=
ln
(
1
/
4
)
and
θ
=
0.01
y
(
0
)
=
ln
(
1
/
4
)
and
θ
=
0.98
32
8.23 (−7)a
1.10 (−6)
6.11 (−7)
8.11 (−7)
64
5.05 (−8)
4.03
6.80 (−8)
4.02
3.50 (−8)
4.13
4.60 (−8)
4.14
128
3.13 (−9)
4.22 (−9)
2.09 (−9)
2.72 (−9)
256
1.94 (−10)
4.01
2.64 (−10)
4.00
1.27 (−10)
4.04
1.67 (−10)
4.03
y
'
(
0
)
=
0
and
θ
=
0.01
y
'
(
0
)
=
0
and
θ
=
0.98
32
3.37 (−7)
4.03 (−7)
1.72 (−7)
1.95 (−7)
64
2.04 (−8)
4.05
2.46 (−8)
4.03
9.13 (−9)
4.24
9.80 (−9)
4.31
128
1.25 (−9)
1.52 (−9)
5.20 (−10)
5.39 (−10)
266
7.74 (−11)
4.02
9.56 (−11)
3.99
3.09 (−11)
4.08
3.23 (−11)
4.06
a
8.23
(
-
7
)
=
8.23
×
10
-
7.The work is organized in the following manner. In Section2, first the method is described and then its construction is explained. In Section 3, convergence of the method is established and in Section 4, the order of the method is corroborated through one example.


## 2. Finite Difference Method

This section is divided in to two parts: (i) description of the method and (ii) derivation of the method. The coefficients not specified explicitly in this section are specified in the appendices.


### 2.1. Description of the Method

In this section we first state the method; detailed derivation is given in Section2.2.For positive integerN≥2, we consider the mesh wh={xk}k=0N over [0,1]: 0=x0<x1<x2<⋯xN=1, h=1/N. For uniform mess xk=kh, where b0 is mentioned in condition (A-1). Let r(x):=f(x,y(x))  (y(x) is the solution),yk=y(xk), and so forth. Now we approximate the differential equation (1.1) with boundary condition (1.2) on the grid wh by the following difference equations:
(2.1)y~k-1Jk-1-(1Jk+1Jk-1)y~k+y~k+1Jk=a0,kr~k+a1,kr~k+1-a-1,kr~k-1,k=1(1)(N-1),(2.2)y~N-1JN-1-(1JN-1+α(βGN))y~N+γ(βGN)=a0,Nr~N-a-1,Nr~N-1+a1,Nr~N-θ,
for boundary value problem (1.1)-(1.2). In the case of boundary condition (1.3) we use the following difference equation:
(2.3)-y~1J1+y~2J1=a-1,1r~0+a0,1r~1+a1,1r~2.
Here y~=(y~k) denotes the approximate solution,y~k≈yk,Gk=G(xk), and
(2.4)Jk=∫xkxk+1‍(p(τ))-1dτ.
Note that the functions p(x) and G(x) are defined in assumption (A-1). Thus the method for boundary value problem (1.1)-(1.2) is given by difference equations (2.1)-(2.2). Further, for boundary value problem (1.1) and (1.3) the method is given by difference equations (2.1) (with k=2(1)(N-1)) and (2.2)-(2.3).


### 2.2. Construction of the Method

In this section we describe the derivation of the method. Local truncation errors are mentioned without proof.Withz(x)=p(x)y′(x) the differential equation (1.1) becomes z′(x)=q(x)f(x,y(x)). Now integrating z′(x)=q(x)f(x,y(x)) twice, first from xk to x and then from xk to xk+1 and changing the order of integration we get
(2.5)yk+1-yk=zkJk+∫xkxk+1‍(∫txk+1‍(p(τ))-1dτ)q(t)r(t)dt,
where zk=z(xk),r(x):=f(t,y(x)) and Jk=∫xkxk+1‍(p(τ))-1dτ. In an analogous manner, we get
(2.6)yk-yk-1=zkJk-1-∫xk-1xk‍(∫xk-1t‍(p(τ))-1dτ)q(t)r(t)dt.
Eliminating zk from (2.5) and (2.6) we obtain the Chawla's identity
(2.7)yk+1-ykJk-yk-yk-1Jk-1=Ik+Jk+Ik-Jk-1,k=1(1)(N-1),
where
(2.8)Ik+=∫xkxk+1‍(∫txk+1‍(p(τ))-1dτ)q(t)r(t)dt,Ik-=∫xk-1xk‍(∫xk-1t‍(p(τ))-1dτ)q(t)r(t)dt.
Via Taylor expansion of G(x), r(x) and H(x) about xk in Ik± and taking the following approximations for rk′ and rk′′(2.9)rk′=rk+1-rk-12h-16h2r′′′(ηk),xk-1<ηk<xk+1,rk′′=rk+1-2rk+rk-1h2-112h2riv(ξk),xk-1<ξk<xk+1,
we get the approximation for the smooth solution y(x) as:
(2.10)y~k-1Jk-1-(1Jk+1Jk-1)y~k+y~k+1Jk=a0,krk+a1,krk+1-a-1,krk-1+tk,k=1(1)(N-1),
where tk is local truncation error given by
(2.11)tk=(HkGk′′′b03,k+HkGiv(ξk)b04,k)fk+(Hkfk′+Hk′fk)(Gk′′b12,k+Gk′′′b13,k)+(12Hkfk′′+Hk′fk′+12Hk′′fk)(Gk′b21,k+Gk′′b22,k)+(16Hkfk′′′+12Hk′fk′′+12Hk′′fk′+16Hk′′′fk)(Gkb30,k+Gk′b31,k)+(124Hkfiv(ξk)+16Hk′fk′′′+14Hk′′fk′′+16Hk′′′fk′+124Hiv(ηk)fk)Gkb40,k-16h2f′′′(ηk)(b10,kGkHk+b11,kGk′Hk+b20,kGkHk′)-124b20,kh2GkHkfiv(ξk),k=1(1)(N-1),xk-1<ξk,ηk<xk+1.
The coefficients b03,k,b04,k,b10,k,b11,k,b12,k,b13,k,b20,k,b21,k,b22,k,b30,k,b31,k,b40,k, and so forth, are specified in the appendices and the functions G(x), H(x) are defined in the assumptions (A-1), (A-2).


#### 2.2.1. Discretization of the Boundary Condition atx=1

We write (2.6) for k=N(2.12)yNJN-1-yN-1JN-1=pNyN′-IN-JN-1.Now, we use the boundary condition atx=1, Taylor expansion of r(x) and following approximations for rN′ and rN′′(2.13)rN′=1h[(1+θ)θrN+θ(1-θ)rN-1+1θ(θ-1)rN-θ]-h2θ6r′′′(ηN),xN-1<ηN<xN,rN′′=2h2[1θrN+1(1-θ)rN-1+1θ(θ-1)rN-θ]-h(1+θ)3r′′′(ηN′),xN-1<ηN′<xN;0<θ<1
in (2.12) and get the discretization at x=1 as follows:
(2.14)y~N-1JN-1-(1JN-1+α(βGN))y~N+γ(βGN)=a0,NrN-a-1,NrN-1+a1,NrN-θ+tN,
where the local truncation error tN is given by
(2.15)tN=1JN-1[h(1+θ)3HNfNGN′′′a03,N-+HNfNG(ξN-)iva04,N-+(HNfN′+HN′fN){GN′′a12,N-+GN′′′a13,N-}1JN-1+(12HNfN′′+HN′fN′+12HN′′fN)(GN′a21,N-+GN′′a22,N-)1JN-1+(16HNfN′′′+12HN′fN′′+12HN′′fN′+16HN′′′fN)(GNa30,N-+GN′a31,N-)1JN-1+(124HNf(ξN)iv+16HN′fN′′′+14HN′′fN′′+16HN′′′fN′+124H(ηN)ivfN)GNa40,k-1JN-1-θ6h2f′′′(ηN-)(a10,N-GNHN+a11,N-GN′HN+a20,N-GNHN′)1JN-1-h(1+θ)3a20,N-h2GNHNf′′′(ηN-)],xN-1<ξN-,ηN-<xN.
The coefficients aij,N- are specified in the appendices.The discretization (2.14) involves unknown yN-θ, which we approximate in the following way:
(2.16)yN-θ=y¯N-θ-h3θ2(1-θ)6yN′′′,
where
(2.17)y¯N-θ=θ(1-θ)[θ(1-θ)yN-1+((1+θ)θ+hαβ)-hγβ].
Now, let r¯N-θ=f(xN-θ,y¯N-θ) then replacing rN-θ by r¯N-θ in (2.14) we obtain the following discretization for a smooth solution y(x) at k=N(2.18)y~N-1JN-1-(1JN-1+α(βGN))y~N+γ(βGN)=a0,NrN-a-1,NrN-1+a1,Nr¯N-θ+tN¯,
where
(2.19)tN¯=tN(1)-h36θ2(1-θ)a1,NyN′′′∂f∂y(xN-θ,yN-θ*),yN-θ*∈(min{yN-θ,y¯N-θ},max{yN-θ,y¯N-θ}).


#### 2.2.2. Discretization of the Boundary Condition atx=0

Integrating the differential equation (1.1) twice, first from x1 to x; then from x1 to x2 and interchanging the order of integration we get
(2.20)y2-y1=J1∫0x1‍q(t)r(t)dt+∫x1x2‍(∫tx2‍(p(τ))-1dτ)q(t)r(t)dt.
Via Taylor expansion of r(x),H(x), and G(x) about x=x1 we get the following discretization for k=1:
(2.21)-y~1J1+y~2J1=a-1,1r0+a0,1r1+a1,1r2+t1,
where the local truncation error t1 is given by
(2.22)t1=(16H1f1′′′+12H1′f1′′+12H1′′f1′+16H1′′′f1)(-6x14+a0ψ(4))+(124H1f1iv+16H1′f1′′′+14H1′′f1′′+16H1′′′f1′+124H1ivf1)(24x14+a0ψ(5))+1J1{H1f1(16a03,1+G1′′′+124a04,1+G1iv)+1J1+(H1f1′+H1′f1)(12a12,1+G1′′+16a13,1+G1′′′)+1J1+(12H1f1′′+H1′f1′+12H1′′f1)(a21,1+G1′+12a22,1+G1′′)+1J1+(16H1f1′′′+12H1′f1′′+12H1′′f1′+16H1′′′f1)(a30,1+G1+a31,1+G1′)+1J1+(124H1f1iv+16H1′f1′′′+14H1′′f1′′+16H1′′′f1′+124H1ivf1)a40,1+G1}+-h26f′′′(η1)[-x12+a0ψ(2)H1+2x13+a0ψ(3)H1′+-h26f′′′(η1)+1J1{(a10,1+G1+a11,1+G1′)H1+a20,1+G1H1′}x12+a0ψ(2)]-h212fiv(ξ1)[x13+a0ψ(3)H1+12J1a20,1+G1H1],0<ξ1<x2.
The coefficients ai,j are specified in the appendices and the functions G(x),H(x) are defined in the assumptions (A-1) and (A-2).The discretization (2.22) involves y0; so for y0 we use the following approximation:
(2.23)y¯0=y1+x12+a0-b0(1+a0)(2+a0-b0)H1G1f1+x13+a0-b0{(-1(1+a0)(2-b0)+1(2+a0-b0)(3+a0-b0)+1(2-b0)(3+a0-b0))G1′H1f1+x13+a0-b0+(-1(1+a0)(2+a0)+1(3+a0-b0))G1H1′f1}+x12+a0-b02(-1(1+a0)(2+a0)+1(3+a0-b0))(f2-f0),
where
(2.24)y0=y¯0+τ0,τ0=c*h4+a0-b0,c*issomeconstant.
Let r¯0=f(x0,y¯0); then in (2.21) replacing r0 by r¯0 we obtain the following discretization for k=1:
(2.25)1J1y1-1J1y2+b-1,1r¯0+b0,1r1+b1,1r2+t¯1(1)=0,
where t¯1=t1+b-1,1τ0(∂f/∂y)(x0,y0*),y0*∈(min{y0,y¯0},max{y0,y¯0}), and coefficients bi,j, are specified in appendices.


## 2.1. Description of the Method

In this section we first state the method; detailed derivation is given in Section2.2.For positive integerN≥2, we consider the mesh wh={xk}k=0N over [0,1]: 0=x0<x1<x2<⋯xN=1, h=1/N. For uniform mess xk=kh, where b0 is mentioned in condition (A-1). Let r(x):=f(x,y(x))  (y(x) is the solution),yk=y(xk), and so forth. Now we approximate the differential equation (1.1) with boundary condition (1.2) on the grid wh by the following difference equations:
(2.1)y~k-1Jk-1-(1Jk+1Jk-1)y~k+y~k+1Jk=a0,kr~k+a1,kr~k+1-a-1,kr~k-1,k=1(1)(N-1),(2.2)y~N-1JN-1-(1JN-1+α(βGN))y~N+γ(βGN)=a0,Nr~N-a-1,Nr~N-1+a1,Nr~N-θ,
for boundary value problem (1.1)-(1.2). In the case of boundary condition (1.3) we use the following difference equation:
(2.3)-y~1J1+y~2J1=a-1,1r~0+a0,1r~1+a1,1r~2.
Here y~=(y~k) denotes the approximate solution,y~k≈yk,Gk=G(xk), and
(2.4)Jk=∫xkxk+1‍(p(τ))-1dτ.
Note that the functions p(x) and G(x) are defined in assumption (A-1). Thus the method for boundary value problem (1.1)-(1.2) is given by difference equations (2.1)-(2.2). Further, for boundary value problem (1.1) and (1.3) the method is given by difference equations (2.1) (with k=2(1)(N-1)) and (2.2)-(2.3).


## 2.2. Construction of the Method

In this section we describe the derivation of the method. Local truncation errors are mentioned without proof.Withz(x)=p(x)y′(x) the differential equation (1.1) becomes z′(x)=q(x)f(x,y(x)). Now integrating z′(x)=q(x)f(x,y(x)) twice, first from xk to x and then from xk to xk+1 and changing the order of integration we get
(2.5)yk+1-yk=zkJk+∫xkxk+1‍(∫txk+1‍(p(τ))-1dτ)q(t)r(t)dt,
where zk=z(xk),r(x):=f(t,y(x)) and Jk=∫xkxk+1‍(p(τ))-1dτ. In an analogous manner, we get
(2.6)yk-yk-1=zkJk-1-∫xk-1xk‍(∫xk-1t‍(p(τ))-1dτ)q(t)r(t)dt.
Eliminating zk from (2.5) and (2.6) we obtain the Chawla's identity
(2.7)yk+1-ykJk-yk-yk-1Jk-1=Ik+Jk+Ik-Jk-1,k=1(1)(N-1),
where
(2.8)Ik+=∫xkxk+1‍(∫txk+1‍(p(τ))-1dτ)q(t)r(t)dt,Ik-=∫xk-1xk‍(∫xk-1t‍(p(τ))-1dτ)q(t)r(t)dt.
Via Taylor expansion of G(x), r(x) and H(x) about xk in Ik± and taking the following approximations for rk′ and rk′′(2.9)rk′=rk+1-rk-12h-16h2r′′′(ηk),xk-1<ηk<xk+1,rk′′=rk+1-2rk+rk-1h2-112h2riv(ξk),xk-1<ξk<xk+1,
we get the approximation for the smooth solution y(x) as:
(2.10)y~k-1Jk-1-(1Jk+1Jk-1)y~k+y~k+1Jk=a0,krk+a1,krk+1-a-1,krk-1+tk,k=1(1)(N-1),
where tk is local truncation error given by
(2.11)tk=(HkGk′′′b03,k+HkGiv(ξk)b04,k)fk+(Hkfk′+Hk′fk)(Gk′′b12,k+Gk′′′b13,k)+(12Hkfk′′+Hk′fk′+12Hk′′fk)(Gk′b21,k+Gk′′b22,k)+(16Hkfk′′′+12Hk′fk′′+12Hk′′fk′+16Hk′′′fk)(Gkb30,k+Gk′b31,k)+(124Hkfiv(ξk)+16Hk′fk′′′+14Hk′′fk′′+16Hk′′′fk′+124Hiv(ηk)fk)Gkb40,k-16h2f′′′(ηk)(b10,kGkHk+b11,kGk′Hk+b20,kGkHk′)-124b20,kh2GkHkfiv(ξk),k=1(1)(N-1),xk-1<ξk,ηk<xk+1.
The coefficients b03,k,b04,k,b10,k,b11,k,b12,k,b13,k,b20,k,b21,k,b22,k,b30,k,b31,k,b40,k, and so forth, are specified in the appendices and the functions G(x), H(x) are defined in the assumptions (A-1), (A-2).


### 2.2.1. Discretization of the Boundary Condition atx=1

We write (2.6) for k=N(2.12)yNJN-1-yN-1JN-1=pNyN′-IN-JN-1.Now, we use the boundary condition atx=1, Taylor expansion of r(x) and following approximations for rN′ and rN′′(2.13)rN′=1h[(1+θ)θrN+θ(1-θ)rN-1+1θ(θ-1)rN-θ]-h2θ6r′′′(ηN),xN-1<ηN<xN,rN′′=2h2[1θrN+1(1-θ)rN-1+1θ(θ-1)rN-θ]-h(1+θ)3r′′′(ηN′),xN-1<ηN′<xN;0<θ<1
in (2.12) and get the discretization at x=1 as follows:
(2.14)y~N-1JN-1-(1JN-1+α(βGN))y~N+γ(βGN)=a0,NrN-a-1,NrN-1+a1,NrN-θ+tN,
where the local truncation error tN is given by
(2.15)tN=1JN-1[h(1+θ)3HNfNGN′′′a03,N-+HNfNG(ξN-)iva04,N-+(HNfN′+HN′fN){GN′′a12,N-+GN′′′a13,N-}1JN-1+(12HNfN′′+HN′fN′+12HN′′fN)(GN′a21,N-+GN′′a22,N-)1JN-1+(16HNfN′′′+12HN′fN′′+12HN′′fN′+16HN′′′fN)(GNa30,N-+GN′a31,N-)1JN-1+(124HNf(ξN)iv+16HN′fN′′′+14HN′′fN′′+16HN′′′fN′+124H(ηN)ivfN)GNa40,k-1JN-1-θ6h2f′′′(ηN-)(a10,N-GNHN+a11,N-GN′HN+a20,N-GNHN′)1JN-1-h(1+θ)3a20,N-h2GNHNf′′′(ηN-)],xN-1<ξN-,ηN-<xN.
The coefficients aij,N- are specified in the appendices.The discretization (2.14) involves unknown yN-θ, which we approximate in the following way:
(2.16)yN-θ=y¯N-θ-h3θ2(1-θ)6yN′′′,
where
(2.17)y¯N-θ=θ(1-θ)[θ(1-θ)yN-1+((1+θ)θ+hαβ)-hγβ].
Now, let r¯N-θ=f(xN-θ,y¯N-θ) then replacing rN-θ by r¯N-θ in (2.14) we obtain the following discretization for a smooth solution y(x) at k=N(2.18)y~N-1JN-1-(1JN-1+α(βGN))y~N+γ(βGN)=a0,NrN-a-1,NrN-1+a1,Nr¯N-θ+tN¯,
where
(2.19)tN¯=tN(1)-h36θ2(1-θ)a1,NyN′′′∂f∂y(xN-θ,yN-θ*),yN-θ*∈(min{yN-θ,y¯N-θ},max{yN-θ,y¯N-θ}).


### 2.2.2. Discretization of the Boundary Condition atx=0

Integrating the differential equation (1.1) twice, first from x1 to x; then from x1 to x2 and interchanging the order of integration we get
(2.20)y2-y1=J1∫0x1‍q(t)r(t)dt+∫x1x2‍(∫tx2‍(p(τ))-1dτ)q(t)r(t)dt.
Via Taylor expansion of r(x),H(x), and G(x) about x=x1 we get the following discretization for k=1:
(2.21)-y~1J1+y~2J1=a-1,1r0+a0,1r1+a1,1r2+t1,
where the local truncation error t1 is given by
(2.22)t1=(16H1f1′′′+12H1′f1′′+12H1′′f1′+16H1′′′f1)(-6x14+a0ψ(4))+(124H1f1iv+16H1′f1′′′+14H1′′f1′′+16H1′′′f1′+124H1ivf1)(24x14+a0ψ(5))+1J1{H1f1(16a03,1+G1′′′+124a04,1+G1iv)+1J1+(H1f1′+H1′f1)(12a12,1+G1′′+16a13,1+G1′′′)+1J1+(12H1f1′′+H1′f1′+12H1′′f1)(a21,1+G1′+12a22,1+G1′′)+1J1+(16H1f1′′′+12H1′f1′′+12H1′′f1′+16H1′′′f1)(a30,1+G1+a31,1+G1′)+1J1+(124H1f1iv+16H1′f1′′′+14H1′′f1′′+16H1′′′f1′+124H1ivf1)a40,1+G1}+-h26f′′′(η1)[-x12+a0ψ(2)H1+2x13+a0ψ(3)H1′+-h26f′′′(η1)+1J1{(a10,1+G1+a11,1+G1′)H1+a20,1+G1H1′}x12+a0ψ(2)]-h212fiv(ξ1)[x13+a0ψ(3)H1+12J1a20,1+G1H1],0<ξ1<x2.
The coefficients ai,j are specified in the appendices and the functions G(x),H(x) are defined in the assumptions (A-1) and (A-2).The discretization (2.22) involves y0; so for y0 we use the following approximation:
(2.23)y¯0=y1+x12+a0-b0(1+a0)(2+a0-b0)H1G1f1+x13+a0-b0{(-1(1+a0)(2-b0)+1(2+a0-b0)(3+a0-b0)+1(2-b0)(3+a0-b0))G1′H1f1+x13+a0-b0+(-1(1+a0)(2+a0)+1(3+a0-b0))G1H1′f1}+x12+a0-b02(-1(1+a0)(2+a0)+1(3+a0-b0))(f2-f0),
where
(2.24)y0=y¯0+τ0,τ0=c*h4+a0-b0,c*issomeconstant.
Let r¯0=f(x0,y¯0); then in (2.21) replacing r0 by r¯0 we obtain the following discretization for k=1:
(2.25)1J1y1-1J1y2+b-1,1r¯0+b0,1r1+b1,1r2+t¯1(1)=0,
where t¯1=t1+b-1,1τ0(∂f/∂y)(x0,y0*),y0*∈(min{y0,y¯0},max{y0,y¯0}), and coefficients bi,j, are specified in appendices.


## 2.2.1. Discretization of the Boundary Condition atx=1

We write (2.6) for k=N(2.12)yNJN-1-yN-1JN-1=pNyN′-IN-JN-1.Now, we use the boundary condition atx=1, Taylor expansion of r(x) and following approximations for rN′ and rN′′(2.13)rN′=1h[(1+θ)θrN+θ(1-θ)rN-1+1θ(θ-1)rN-θ]-h2θ6r′′′(ηN),xN-1<ηN<xN,rN′′=2h2[1θrN+1(1-θ)rN-1+1θ(θ-1)rN-θ]-h(1+θ)3r′′′(ηN′),xN-1<ηN′<xN;0<θ<1
in (2.12) and get the discretization at x=1 as follows:
(2.14)y~N-1JN-1-(1JN-1+α(βGN))y~N+γ(βGN)=a0,NrN-a-1,NrN-1+a1,NrN-θ+tN,
where the local truncation error tN is given by
(2.15)tN=1JN-1[h(1+θ)3HNfNGN′′′a03,N-+HNfNG(ξN-)iva04,N-+(HNfN′+HN′fN){GN′′a12,N-+GN′′′a13,N-}1JN-1+(12HNfN′′+HN′fN′+12HN′′fN)(GN′a21,N-+GN′′a22,N-)1JN-1+(16HNfN′′′+12HN′fN′′+12HN′′fN′+16HN′′′fN)(GNa30,N-+GN′a31,N-)1JN-1+(124HNf(ξN)iv+16HN′fN′′′+14HN′′fN′′+16HN′′′fN′+124H(ηN)ivfN)GNa40,k-1JN-1-θ6h2f′′′(ηN-)(a10,N-GNHN+a11,N-GN′HN+a20,N-GNHN′)1JN-1-h(1+θ)3a20,N-h2GNHNf′′′(ηN-)],xN-1<ξN-,ηN-<xN.
The coefficients aij,N- are specified in the appendices.The discretization (2.14) involves unknown yN-θ, which we approximate in the following way:
(2.16)yN-θ=y¯N-θ-h3θ2(1-θ)6yN′′′,
where
(2.17)y¯N-θ=θ(1-θ)[θ(1-θ)yN-1+((1+θ)θ+hαβ)-hγβ].
Now, let r¯N-θ=f(xN-θ,y¯N-θ) then replacing rN-θ by r¯N-θ in (2.14) we obtain the following discretization for a smooth solution y(x) at k=N(2.18)y~N-1JN-1-(1JN-1+α(βGN))y~N+γ(βGN)=a0,NrN-a-1,NrN-1+a1,Nr¯N-θ+tN¯,
where
(2.19)tN¯=tN(1)-h36θ2(1-θ)a1,NyN′′′∂f∂y(xN-θ,yN-θ*),yN-θ*∈(min{yN-θ,y¯N-θ},max{yN-θ,y¯N-θ}).


## 2.2.2. Discretization of the Boundary Condition atx=0

Integrating the differential equation (1.1) twice, first from x1 to x; then from x1 to x2 and interchanging the order of integration we get
(2.20)y2-y1=J1∫0x1‍q(t)r(t)dt+∫x1x2‍(∫tx2‍(p(τ))-1dτ)q(t)r(t)dt.
Via Taylor expansion of r(x),H(x), and G(x) about x=x1 we get the following discretization for k=1:
(2.21)-y~1J1+y~2J1=a-1,1r0+a0,1r1+a1,1r2+t1,
where the local truncation error t1 is given by
(2.22)t1=(16H1f1′′′+12H1′f1′′+12H1′′f1′+16H1′′′f1)(-6x14+a0ψ(4))+(124H1f1iv+16H1′f1′′′+14H1′′f1′′+16H1′′′f1′+124H1ivf1)(24x14+a0ψ(5))+1J1{H1f1(16a03,1+G1′′′+124a04,1+G1iv)+1J1+(H1f1′+H1′f1)(12a12,1+G1′′+16a13,1+G1′′′)+1J1+(12H1f1′′+H1′f1′+12H1′′f1)(a21,1+G1′+12a22,1+G1′′)+1J1+(16H1f1′′′+12H1′f1′′+12H1′′f1′+16H1′′′f1)(a30,1+G1+a31,1+G1′)+1J1+(124H1f1iv+16H1′f1′′′+14H1′′f1′′+16H1′′′f1′+124H1ivf1)a40,1+G1}+-h26f′′′(η1)[-x12+a0ψ(2)H1+2x13+a0ψ(3)H1′+-h26f′′′(η1)+1J1{(a10,1+G1+a11,1+G1′)H1+a20,1+G1H1′}x12+a0ψ(2)]-h212fiv(ξ1)[x13+a0ψ(3)H1+12J1a20,1+G1H1],0<ξ1<x2.
The coefficients ai,j are specified in the appendices and the functions G(x),H(x) are defined in the assumptions (A-1) and (A-2).The discretization (2.22) involves y0; so for y0 we use the following approximation:
(2.23)y¯0=y1+x12+a0-b0(1+a0)(2+a0-b0)H1G1f1+x13+a0-b0{(-1(1+a0)(2-b0)+1(2+a0-b0)(3+a0-b0)+1(2-b0)(3+a0-b0))G1′H1f1+x13+a0-b0+(-1(1+a0)(2+a0)+1(3+a0-b0))G1H1′f1}+x12+a0-b02(-1(1+a0)(2+a0)+1(3+a0-b0))(f2-f0),
where
(2.24)y0=y¯0+τ0,τ0=c*h4+a0-b0,c*issomeconstant.
Let r¯0=f(x0,y¯0); then in (2.21) replacing r0 by r¯0 we obtain the following discretization for k=1:
(2.25)1J1y1-1J1y2+b-1,1r¯0+b0,1r1+b1,1r2+t¯1(1)=0,
where t¯1=t1+b-1,1τ0(∂f/∂y)(x0,y0*),y0*∈(min{y0,y¯0},max{y0,y¯0}), and coefficients bi,j, are specified in appendices.


## 3. Convergence of the Method

In this section we show that under suitable conditions the method (2.1)-(2.2) for the boundary value problems (1.1)-(1.2) described in Section 2 is of fourth order accuracy.LetY~=(y~1,…,y~N)T,R(Y~)=(r~1,…,r~N)T,Q=(q1,…,qN)T, and T=(t1,…,tN)T, then the finite difference method given by (2.1)-(2.2) for the boundary value problems (1.1)-(1.2) can be expressed in matrix form as
(3.1)BY~+PR(Y~)+GY~=Q,
and Y=(y1,…,yN)T corresponding to smooth solution y(x) satisfies the perturbed system
(3.2)BY+PR(Y)+GY+T=Q,
where B=(bij) and P=(pij) are (N×N) tridiagonal matrices.LetE=Y~-Y=(e1,…,eN)T, then from Mean Value Theorem R(Y~)-R(Y)=ME,M=diag{U1,…,UN}(Uk=∂fk/∂yk≥0), then from (3.1) and (3.2) we get the error equation as
(3.3)(B+PM+W)E=T.We recall that by the notationZ≥0 we mean that all components zi of the vector Z satisfy zi≥0. Similarly by B≥0 we mean all the elements bij of the matrix B satisfy bij≥0.From Corollary of Theorems 7.2 and 7.4 in [25] an (irreducible) tridiagonal matrix B is monotone if
(3.4)bk,k+1≤0,bk,k-1≤0,∑j=1Nbk,j{≥0,k=1,2,…,N,>0,foratleastonei.
Furthermore, from Theorem 7.3 in [25], B-1 exists and the elements of B-1 are nonnegative. Now, if B+PM+W is monotone and B+PM+W≥B, then from Theorem 7.5 in [25]
(3.5)(B+PM+W)-1≤B-1.
Let vector norm ∥E∥ and matrix norm ∥B∥ be defined by
(3.6)∥E∥=max1≤k≤N|ek|,∥B∥=max1≤k≤N∑j=1N‍|bkj|,
then from (3.3) we get
(3.7)∥E∥=∥(B+PM+W)-1|T|∥,
where |T|≡(|t1|,|t2|,…,|tN|).Furthermore, ifB+PM+W≥B and B is also monotone matrix, then from (3.5) we get
(3.8)∥E∥≤∥B-1|T|∥.
Now from [22],B-1=(bi,j-1)  of matrix B is given by
(3.9)bi,j-1={P(xi)[P(xN)-P(xj)+βGN/α][P(xN)+βGN/α],i≤j,P(xj)[P(xN)-P(xi)+βGN/α][P(xN)+βGN/α],i≥j,
where
(3.10)P(x)=∫0x‍(p(τ))-1dτ.
Let xriv,r(i) for i=0(1)3, and y′′′ be bounded on (0,1]. Now, as b0 is fixed in [0,1) and a0 is chosen such that 1+a0-b0>0, for sufficiently small h we get
(3.11)|tk|≤Ch5xk-1+a0,|t¯N|≤C~h4,
for suitable constants C and C~.Now from (3.8)–(3.11) we get
(3.12)|ei|≤C*h5[{1-P(xi)[P(xN)+βGN/α]}∑j=1i‍xj-1+a0P(xj)xj-1+a0C*h5+P(xi)[P(xN)+βGN/α]{∑j=i+1N-1βGNhα‍[P(xN)-P(xj)+βGNα]xj-1+a0+βGNhα}],
where C*=max{C,C-}. It is easy to see that
(3.13)h∑j=1i‍∫0xi‍P(xj)xj-1+a0<∫0xi‍P(x)x-1+a0dx,h∑j=i+1N-1‍βGNα[P(xN)-P(xj)+βGNα]xj-1+a0<∫xixNβGNα‍[P(xN)-P(x)+βGNα]x-1+a0dx.
From (3.12)-(3.13) we obtain
(3.14)|ei|≤C*h4(1-b0){(1-b0)xi1-b0-a0(1-b0+a0)[xia0-α+(1-b0+a0)(1-a0)βα+(1-b0)β]×sup[0,1]|G(x)|+(1-b0)d1,a0≠0xi1-b0[ln(1xi)+(2-b0)βα+(1-b0)β]sup[0,1]|G(x)|+d2,a0=0,
where
(3.15)d1=α(1-b0+a0)(1-b0+a0)/-a0(α+(1-b0)β)(1-b0)(1-b0)/-a0×[(2-b0+a0)(1-b0+a0)/-a0(2-b0)(1-b0)/-a0sup[0,1]|G′(x)|+(3-b0+a0)(1-b0+a0)/-a02(3-b0)(1-b0)/-a0sup[0,1]|G′′(x)|+(4-b0+a0)(1-b0+a0)/-a03!(4-b0)(1-b0)/-a0sup[0,1]|G′′′(x)|+(5-b0+a0)(1-b0+a0)/-a04!(5-b0)(1-b0)/-a0sup[0,1]|G'v(x)|]d2=α(α+(1-b0)β)e2[1(2-b0)e1/(2-b0)sup[0,1]|G′(x)|+12(3-b0)e2/(3-b0)sup[0,1]|G′′(x)|α(α+(1-b0)β)e2+13!(4-b0)e3/(4-b0)sup[0,1]|G′′′(x)|+14!(5-b0)e4/(5-b0)sup[0,1]|G'v(x)|].
In view of the following inequalities
(3.16)x1-b0(xa0-α+(1-b0+a0)(1-a0)βα+(1-b0)β)≤(α+(1-b0)βα+(1-b0+a0)(1-a0)β)(1-b0+a0)/-a0,x1-b0(ln(1x)+βα+(1-b0)β)≤1(1-b0)e-[α-(1-b0)β]/[α+(1-b0)β]
for x∈(0,1],(B+PM+W)-1≤B-1 and from (3.14) it is easy to establish that
(3.17)∥E∥∞=O(h4).
Thus we have established the following result.Theorem 3.1.
Assume thatp(x), q(x), and f(x,y) satisfy assumptions given in (A-1),(A-2), and (A-3), respectively. Let r(x):=f(x,y(x)), then the method is of fourth order accuracy for sufficiently small mesh size h provided, xr(iv)(x),r(i),i=0(1)3, and y′′′ are bounded on (0,1].Remark 3.2.
We have developed the numerical method for the boundary conditiony′(0)=0 but could not establish the fourth order convergence although the order of the accuracy is verified through one example in the next section.


## 4. Numerical Illustrations

To illustrate the convergence of the method and to corroborate their order of the accuracy, we apply the method to following example. The maximum absolute errors are displayed in Table1.Example 4.1.
Consider(4.1)(xb011+x3.5y′)′=(1.9-b0)x-0.1((1.9-b0)x1.9-b0ey-0.9+3.5x3.5/(1+x3.5))(4+x1.9-b0)(1+x3.5),y(0)=ln(14)or(y′(0)=0),y(1)+5y′(1)=ln(15)-(1.9-b0),
with exact solutiony(x)=ln(1/(4+x1.9-b0)).The method is applied on the example forb0=0.1,0.5. Maximum absolute errors and order of accuracy for Example 4.1 are displayed in Table 1.We use the following approximation ofJk in our numerical program:
(4.2)Jk=xk+11-b0-xk1-b01-b0Gk+(xk+12-b0-xk2-b02-b0-Gk(xk+11-b0-xk1-b0)1-b0)Gk′+12Gk′′(xk+13-b0-xk3-b03-b0+xk2xk+11-b0-xk1-b01-b0-2xkxk+12-b0-xk2-b02-b0)+16Gk′′′(xk+14-b0-xk4-b04-b0-3xkxk+13-b0-xk3-b03-b0+3xk2xk+12-b0-xk2-b02-b0-xk3xk+11-b0-xk1-b01-b0).


### 4.1. Numerical Results for Uniform Mesh Case

Remark 4.2.
In the discretization (2.18), θ may take any value in (0,1) but from the Table 1, maximum absolute errors are less for the value of θ close to one. So it is better to take value of θ close to one.


## 4.1. Numerical Results for Uniform Mesh Case

Remark 4.2.
In the discretization (2.18), θ may take any value in (0,1) but from the Table 1, maximum absolute errors are less for the value of θ close to one. So it is better to take value of θ close to one.



---

*Source: 349618-2012-12-05.xml*