# Numerical Scheme for Solving Singular Two-Point Boundary Value Problems

Authors: N. Ratib Anakira; A. K. Alomari; I. Hashim
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

Singular two-point boundary value problems (BVPs) are investigated using a new technique, namely, optimal homotopy asymptotic method (OHAM). OHAM provides a convenient way of controlling the convergence region and it does not need to identify an auxiliary parameter. The effectiveness of the method is investigated by comparing the results obtained with the exact solution, which proves the reliability of the method.


---

## Body

## 1. Introduction

Consider singular two-point boundary value problems (BVPs) of the form(1)1pu′′(x)+1q(x)u′(x)+1r(x)=g(x),0≤t≤1,
subject to the boundary conditions
(2)u(0)=α1,u(1)=β,oru′(0)=α2,u(1)=β,
where p,q,r, and g are continuous functions on (0,1] and the parameters α1,α2, and β are real constants.Problems of the forms (1) and (2) are encountered in the fields of fluid mechanics, reaction-diffusion processes, chemical kinetics, optimal control, and other branches of applied mathematics [1, 2]. Many different numerical methods have been proposed by various authors regarding singular two-point boundary value problems such as variational iteration method (VIM) [3], cubic splines [4], differential transformation method (DTM) [5], the Adomian decomposition method (ADM) [6], continuous genetic algorithm (CGA) [7], sinc-Galerkin method and homotopy perturbation method (HPM) [8–13], and homotopy analysis method (HAM) [14–23]. The existence of a unique solution of (1) and (2) was discussed in [24, 25].Recently, Marinca et al. introduced and developed optimal homotopy asymptotic method (OHAM) in a series of papers [26–29] for the approximate solutions of nonlinear problems. OHAM does not depend on the presence of a small parameter. An advantage of OHAM is that it does not need to identify the ℏ-curve as in HAM for convergence region. In OHAM, the control and adjustment of the convergence region is provided in a convenient way. Furthermore, it has a built-in convergence criteria similar to HAM but with a greater degree of flexibility.In this paper, OHAM is presented to create an approximate analytic solution of singular two-point boundary value problems. The method is directly applied without any linearization and discretizations and without splitting the nonhomogeneous term. The structure of this paper is organized as follows: Section2 is devoted to the analysis of the proposed method, in Section 3, three examples are employed to illustrate the accuracy and computational efficiency of this approach, and lastly, conclusions are given in the last section.


## 2. Analysis of the Method

To illustrate the basic idea of OHAM [30], we consider the following differential equation:(3)L(u(x))+g(x)+N(u(x))=0,B(u,dudx)=0,
where L is the chosen linear operator, N is the linear or nonlinear operator, u(x) is an unknown function, x denotes an independent variable, g(x) is a known function, and B is a boundary operator.According to the basic idea OHAM we construct a homotopyh(v(x,p),p):R×[0,1]→R which satisfies
(4)(1-p)[L(v(x,p))-u0(x)]=H(p)[L(v(x,p))+g(x)+N(v(x,p))],B(v(x,p),∂v(x,p)∂x)=0,
where x∈R and p∈[0,1] is an embedding parameter, H(p) is a nonzero auxiliary function for p≠0, H(0)=0, and v(x,p) is an unknown function. Obviously, when p=0 and p=1 it holds that v(x,0) = u0(x) and v(x,1)=u(x), respectively. Thus, as p varies from 0 to 1, the solution v(x,p) approaches from u0(x) to u(x) where u0(x) is the initial guess that satisfies the linear operator and the boundary conditions
(5)L(u0(x))=0,B(u0,du0dx)=0.
Next, we choose the auxiliary function H(p) in the form
(6)H(p)=pC1+p2C2+p3C3+⋯,
where C1,C2,C3,… are constants which can be determined later. H(p) can be expressed in many forms as reported by Marinca et al. [26–29].To get an approximate solution, we expandv(x,p,Ci) in Taylor’s series about p in the following manner:
(7)v(x,p,Ci)=u0(x)+∑k=1∞‍uk(x,C1,C2,…,Ck)pk.
Substituting (7) into (4) and equating the coefficient of like powers of p, we obtain the following linear equations. The zeroth-order problem is given by (5); the first- and second-order problems are given as
(8)L(u1(x))+g(x)=C1N0(u0(x)),B(u1,du1dx)=0,(9)L(u2(x))-L(u1(x))=C2N0(u0(x))+C1[L(u1(x)+N1(u0(x),u1(x))],B(u2,du2dx)=0.
The general governing equations for uk(x) are
(10)L(uk(x))-L(uk-1(x))=CkN0(u0(x))+∑i=1k-1‍Ci[L(uk-i(x))+Nk-i×(u0(x),u1(x),…,uk-1(x))],B(uk,dukdx)=0,
where k=2,3,… and Nm(u0(x),u1(x),…,um(x)) is the coefficient of pm in the expansion of N(v(x,p)) about the embedding parameter p:
(11)N(v(x,p,Ci))=N0(u0(x))+∑m=1∞‍Nm(u0(x),u1(x),…,um(x))pm.
It has been observed that the convergence of the series (7) depends upon the auxiliary constants C1,C2,C3,…. If it is convergent at p=1, one has
(12)v(x,Ci)=u0(x)+∑k=1∞‍uk(x,C1,C2,…,Ck).
The result of the mth-order approximation is given
(13)u~(x,C1,C2,C3,…,Cm)=u0(x)+∑i=1m‍ui(x,C1,C2,…,Ci).
Substituting (12) into (3) yields the following residual:
(14)R(x,C1,C2,C3,…,Cm)=L(u~(x,C1,C2,C3,…,Cm))+g(x)+N(u~(x,C1,C2,C3,…,Cm)).
If R=0, then u~ will be the exact solution. Generally such a case will not arise for nonlinear problems, but we can minimize the functional
(15)J(C1,C2,C3,…,Cm)=∫ab‍R2(x,C1,C2,C3,…,Cm)dx,
where a and b are the endpoints of the given problem. The unknown constants Ci(i=1,2,3,…,m) can be identified from the conditions
(16)∂J∂C1=∂J∂C2=⋯=∂J∂Cm=0.
With these constants known, the approximate solution (of order m) is well determined.


## 3. Numerical Examples

To illustrate the effectiveness of the OHAM we will consider three examples of singular two-point BVPs.Example 1.
Consider the following singular two-point BVP [2, 31]:
(17)u′′(x)+1xu′(x)+u(x)=g(x),0≤x≤1,
subject to the boundary conditions
(18)u(0)=0,u(1)=0;
the exact solution of this problem in case of
(19)g(x)=4-9x+x2-x3
is given by
(20)u(x)=x2-x3.
According to the OHAM formulation described in the above section, we start with
(21)L[v(x,p)]=xd2v(x,p)dx2+1xdv(x,p)dx,N[v(x,p)]=xd2v(x,p)dx2+dv(x,p)dx+xv(x,p)-xg(x).
Now, apply (4) at p=0 to give the zeroth-order problem as follows:
(22)u0′′(x)=0;
with conditions
(23)u0(0)=0,u0(1)=0,
it gives us
(24)u0(x)=0.
Now, apply (8) to give the first-order problem as follows:
(25)u1′′(x,C1)=-4C1+9xC1-x2C1+x3C1+c1u0(x)+u0′(x)x+C1u0′(x)x-u1′(x)x+C1(u0)′′(x)
subject to the boundary conditions
(26)u1(0)=0,u1(1)=0
and having the solution
(27)u1(x,C1)=1400(9-400x2+400x3-25x4+16x5)C1.
The second-order problem can be defined by (9):
(28)u2′′(x,C1,C2)=-4C2+9xC2-x2C2+x3C2+C2u0(x)+C1u1(x)+C2u0′(x)x+u1′(x)x+C1u1′(x)x-u2′(x)x+C2u0′′(x)+u1′′(x)+C1u1′′(x)
subject to the boundary problem
(29)u2(0)=0,u2(1)=0
and has the solution
(30)u2(x,C1,C2)=(9400-x2+x3-x416+x525)C1+(177744100-1591x21600+x3-x48+2x525-x6576+x71225)C12+1400(9-400x2+400x3-25x4+16x5)C2.
By applying (10), the third-order problem is defined as follows:
(31)u3′′(x,C1,C2,C3)=-4C3+9xC3-x2C3+x3C3+C3u0(x)+C2u1(x)+C1u2(x)+C3u0′(x)x+C2u1′(x)x+u2′(x)x+C1u2′(x)x-u3′(x)x+C3u0′′(x)+C2u1′′(x)+u2′′(x)+C1u2′′(x)
subject to the boundary problem
(32)u3(0)=0,u3(1)=0
and has the solution
(33)u3(x,C1,C2,C3)=9C1400-x2C1+x3C1-x4C116+x5C125+1777C1222050-1591800x2C12+2x3C12-14x4C12+425x5C12-1288x6C12+2x7C121225+22038893C13406425600-694523x2C13705600+x3C13-4791x4C1325600+325x5C13-1192x6C13+3x7C131225-x8C1336864+x9C1399225+9C2400-x2C2+x3C2-x4C216+x5C225+1777C1C222050-1591800x2C1C2+2x3C1C2-14x4C1C2+425x5C1C2-1288x6C1C2+2x7C1C21225+9C3400-x2C3+x3C3-x4C316+x5C325.
Using (24), (27), (30), and (33), the third-order approximate solution by OHAM for p=1 is as follows:
(34)u~(x,C1,C2,C3)=u0(x)+u1(x,C1)+u2(x,C1,C2)+u3(x,C1,C2,C3).
Following the procedure described in Section 2 on the domain between a=0 and b=1, using the residual error,
(35)R=xu~′′(x,C1,C2,C3)+u~′(x,C1,C2,C3)+xu~(x,C1,C2,C3)-x(4-9x+x2-x3).
The less square method can be applied as
(36)J(C1,C2,C3)=∫01‍R2dx,dJdC1=dJdC2=dJdC3.
Thus, the following optimal values of Ci’s are obtained:
(37)C1=-1.076627460,C2=-0.005293863,C3=-0.000381048.
By considering these values our approximate solution becomes
(38)u~(x,C1,C2,C3)=-0.0000109358+0·x+1.00005x2-1.00002x3+7.25718×10-6x4-0.000285432x5+0.000442819x6-0.000208215x7+0.0000338527x8-0.0000125769x9.
It is clear that the given solution is very close to the exact one since the other terms’ approach zero which leads to that the solution is converge. Moreover, Table 1 exhibits the approximate solution obtained by using the OHAM and CGA [7]. It is clear that the obtained results in a our method are in very good agreement with the exact solution, which proves the reliability of the method. In Figure 1 we plot the approximate solution and the exact solution.Table 1
Comparison of the exact solution and the OHAM solution for Example1.
x
Exact solution
OHAM solution (3 terms)
CGA solution [7]
Absolute error
0.1
0.0090
0.0089895718978
0.0089999999973
1.04281000
×
1
0
-
5
0.2
0.0320
0.0319909691543
0.0319999999954
9.03084567
×
1
0
-
6
0.3
0.0630
0.0629929391328
0.0629999999949
7.06086713
×
1
0
-
6
0.4
0.0960
0.0959950247523
0.0959999999950
4.97524765
×
1
0
-
6
0.5
0.1250
0.1249967655986
0.1249999999952
3.23440137
×
1
0
-
6
0.6
0.1440
0.1439978797237
0.1439999999957
2.12027629
×
1
0
-
6
0.7
0.1470
0.1469984038747
0.1469999999965
1.59612524
×
1
0
-
6
0.8
0.1280
0.1279987010935
0.1279999999976
1.29890642
×
1
0
-
6
0.9
0.0810
0.0809992354577
0.0809999999988
7.64542249
×
1
0
-
7Figure 1
Exact and approximate solution using OHAM for Example1.Example 2.
Let us consider the singular two-point BVP [2, 3]:
(39)u′′(x)+1xu′(x)+u(x)=g(x),0≤x≤1,u′(0)=0,u(1)=1716.
The exact solution of this problem in the case of
(40)g(x)=54+x216
is given by
(41)u(x)=1+x216.
According to the OHAM formulation described in the previous section, we start with
(42)L[v(x,p)]=xd2v(x,p)dx2+1xdv(x,p)dxN[v(x,p)]=xd2v(x,p)dx2+dv(x,p)dx+xv(x,p)-xg(x).
Applying OHAM, we have the following zeroth-, first-, second-, and the third-order problem solutions:(43)u0(x)=1716,u1(x,C1)=13C1256-3x2C1264-x4C1256,u2(x,C1,C2)=13C1256-3x2C164-x4C1256+379C129216-35x2C121024-7x4C121024-x6C129216+13C2256-3x2C264,u3(x,C1,C2,C3)=13C1256-3x2C164-x4C1256+379C124608-35x2C12512-7x4C12512-x6C124608+19565C13589824-881x2C13589824-147x4C1316384-11x6C1336864-x8C13589824-13C2256-3x2C264-x4C2256+379C1C24608-35x2C1C2512-7x4C1C2512-x6C1C24608+13C3256-3x2C364-x4C3256.
Using (43) we obtain the following third-order approximate solution by OHAM:
(44)u~(x,C1,C2,C3)=u0(x)+u1(x,C1)+u2(x,C1,C2)+u3(x,C1,C2,C3).
Following the procedure described in Section 2 on the domain between a=0 and b=1, using the residual
(45)R=xu~′′(x,C1,C2,C3)+u~′(x,C1,C2,C3)+xu~(x,C1,C2,C3)-x(54+x216).
The following optimal values of Ci’s are obtained:
(46)C1=-1.0769,C2=-0.00539311,C3=-0.000394546.
By considering these values, our approximate solution becomes
(47)u~(x,C1,C2,C3)=1.0+0.0624971x2+6.3086×10-6x4-6.10932×10-6x6+2.11738×10-6x8.
Therefore, we have the third-order approximate solution of Example 2.Table2 shows a comparison between the OHAM solution and the solutions of VIM [3] both with the exact solution. As it is an evident from the compared results, it was found that OHAM gives better results. In Figure 2, both the approximate solution by using OHAM and the exact solution have been plotted.Table 2
Comparison of the exact solution and the OHAM solution for Example2.
x
Exact solution
OHAM solution (3 terms)
VIM solution [3]
Absolute error
0.1
1.0006
1.0006255449104
0.8665109375000
5.44910481
×
1
0
-
7
0.2
1.0025
1.0025004672985
0.8723000000000
4.67298560
×
1
0
-
7
0.3
1.0056
1.0056253598830
0.8819692708333
3.59883028
×
1
0
-
7
0.4
1.0100
1.0100002486730
0.895550000000
2.48673017
×
1
0
-
7
0.5
1.0156
1.0156251578232
0.9130859375000
1.57823275
×
1
0
-
7
0.6
1.0225
1.0225001009617
0.9346333333333
1.00961799
×
1
0
-
7
0.7
1.0306
1.0306250751620
0.9602609375000
7.51620191
×
1
0
-
8
0.8
1.0400
1.0400000614013
0.9900500000000
6.14013284
×
1
0
-
8
0.9
1.0506
1.0506250362014
1.0240942708333
3.62014673
×
1
0
-
8Figure 2
Exact and approximate solution using OHAM for Example2.Example 3.
Consider the singular two-point BVP [2, 8]
(48)(1-x2)u′′(x)+32(1x-1)u′(x)+(x2-1)u(x)=g(x),0≤x≤1,u′(0)=0,u(1)=0.
The exact solution of this problem in the case of
(49)g(x)=5-29x2+13x22+3x32-x42
is given by
(50)u(x)=x2-x3.
According to the OHAM formulation described in the above section, we start with(51)L[v(x,p)]=d2v(x,p)dx2+dv(x,p)dxN[v(x,p)]=(x-12)d2v(x,p)dx2+32(1-x)dv(x,p)dx+x2-2x2v(x,p)-xg(x).
Applying OHAM, we have the following zeroth-, first-, and second-order solutions:
(52)u0(x)=0,u1(x,C1)=1741C1120+23C1e-23C1e-23xC1+23x2C12-14x3C13+19x4C18-4x5C15+x6C112,u2(x,C1,C2)=1741C1120+23C1e-23C1ex-23xC1+23x2C12-14x3C13+19x4C18-4x5C15+x6C112-4057507C126048-46C12e2-11141C1210e+46C12e(-1-x)+64201C1260e-x+64201xC1260+46xC12e-66271x2C12120-23x2C12e+69x2C124e-x+141841x3C12720+23x3C126e-2623x4C1248+41x5C123-1301x6C12360+1369x7C121680-49x8C12480+x9C12216+1741C2120+23C2e-23C2ex-23xC2+23x2C22-14x3C23+19x4C28-4x5C25+x6C212.
Using (52) and also by adding the solutions of the third-order problems, we obtain the following third-order approximate solution by OHAM:
(53)u~(x,C1,C2,C3)=u0(x)+u1(x,C1)+u2(x,C1,C2)+u3(x,C1,C2,C3).
Following the procedure described in Section 2 on the domain between a=0 and b=1, using the residual
(54)R=(x-x22)u~′′(x,C1,C2,C3)+32(1-x)u~′(x,C1,C2,C3)+(x22-x)u~(x,C1,C2,C3)-x(5-29x2+13x22+3x32-x42).
The following optimal values of Ci’s are obtained:
(55)C1=-2.5182,C2=-0.0304705,C3=-0.00712009.
The approximate solution now becomes
(56)u~(x,C1,C2,C3)=-1.40959×106+2.81716×107e-3-x+51363.2e-2-x+57.9185e-x+1.40766×106x+38729.3e-3-xx-717179·x2+285281·e-3-xx2+808.272e-2-xx2+248195·x3+1844.25e-3-xx3-65742·x4+4841.17e-3-xx4+14295.2x5-2702.04x6+477.666x7-84.9559x8+14.2263x9-1.7433x10+0.117783x11-0.00308038x12.
Table3 exhibits the approximate solution obtained by using the OHAM and He’s HPM [8]. It can be seen that the solution obtained by our procedure is nearly identical with that given by the exact solution, which proves the reliability of the method. In Figure 3 we compare the exact solution and the approximate solution obtained by OHAM.Table 3
Comparison of the exact solution and the OHAM solution for Example3.
x
Exact solution
OHAM solution (3 terms)
OHAM error
He's HPM error [8]
0.1
0.0090
0.0099841073222
9.841
×
1
0
-
4
8.936
×
1
0
-
4
0.2
0.0320
0.0291990486850
2.801
×
1
0
-
3
8.631
×
1
0
-
4
0.3
0.0630
0.0600046093219
2.995
×
1
0
-
3
8.120
×
1
0
-
4
0.4
0.0960
0.0941688334742
1.831
×
1
0
-
3
7.408
×
1
0
-
4
0.5
0.1250
0.1241809960356
8.190
×
1
0
-
4
6.507
×
1
0
-
4
0.6
0.1440
0.1437292946566
2.707
×
1
0
-
4
5.435
×
1
0
-
4
0.7
0.1470
0.14693478562367
6.521
×
1
0
-
5
4.215
×
1
0
-
4
0.8
0.1280
0.1279073412152
9.266
×
1
0
-
5
2.878
×
1
0
-
4
0.9
0.0810
0.0808082770963
1.917
×
1
0
-
4
1.467
×
1
0
-
9Figure 3
Exact and approximate solution using OHAM for Example3.


## 4. Conclusions

In this work, OHAM has been applied successfully to solve singular two-point BVPs. The results which are obtained by using OHAM are in a good agreement with the exact solution as well as the results which are already presented in the literature like CGA, VIM, and HPM. This shows that the method is efficient and reliable for the solution of singular two-point boundary value problems.



---

*Source: 468909-2013-04-15.xml*