New Approach for Solving a Class of Doubly Singular Two-Point Boundary Value Problems Using Adomian Decomposition Method

Abstract

We propose two new modified recursive schemes for solving a class of doubly singular two-point boundary value problems. These schemes are based on Adomian decomposition method (ADM) and new proposed integral operators. We use all the boundary conditions to derive an integral equation before establishing the recursive schemes for the solution components. Thus we develop recursive schemes without any undetermined coefficients while computing successive solution components, whereas several previous recursive schemes have done so. This modification also avoids solving a sequence of nonlinear algebraic or transcendental equations for the undetermined coefficients with multiple roots, which is required to complete calculation of the solution by several earlier modified recursion schemes using the ADM. The approximate solution is computed in the form of series with easily calculable components. The effectiveness of the proposed approach is tested by considering four examples and results are compared with previous known results.

1. Introduction

Consider the following class of doubly singular two-point boundary value problems: 𝑝(𝑡)𝑦^′(𝑡)^′=𝑞(𝑡)𝑓𝑡,𝑦(𝑡),0<𝑡≤1, (1.1) with boundary conditions 𝑦(0)=𝑎₁,𝑦(1)=𝑐₁,or𝑦^′(0)=0,𝑎𝑦(1)+𝑏𝑦^′(1)=𝑐, (1.2) where 𝑎₁, 𝑐₁, 𝑎>0, 𝑏≥0, and 𝑐 are any finite constants. The condition 𝑝(0)=0 characterizes that problem (1.1) is singular and in addition to this if 𝑞(𝑡) is allowed to be discontinuous at 𝑡=0, then problem (1.1) is called doubly singular Bobisud [1]. Consider problem (1.1) with the following conditions on 𝑝(𝑡), 𝑞(𝑡), and 𝑓(𝑡,𝑦(𝑡)).

Type 1. Dirichlet boundary conditions: 𝑦(0)=𝑎₁, 𝑦(1)=𝑐₁(𝐸₁)𝑝∈𝐶[0,1], 𝑝∈𝐶¹(0,1] with 𝑝>0 in (0,1] and 1/𝑝∈𝐿¹(0,1); (𝐸₂)𝑞>0 in (0,1] and ∫¹₀(1/𝑝)(∫¹_𝑥𝑞𝑑𝑠)𝑑𝑥<∞;(𝐸₃) (i) 𝑓∈𝐶([0,1]×ℝ); (ii)  𝜕𝑓/𝜕𝑦 exists and continuous, 𝜕𝑓/𝜕𝑦≥0 for all 0≤𝑡≤1 and all real 𝑦.

Type 2. Mixed type boundary conditions: 𝑦^′(0)=0,  𝑎𝑦(1)+𝑏𝑦^′(1)=𝑐(𝐸^′₁) 𝑝∈𝐶[0,1],  𝑝∈𝐶¹(0,1] with 𝑝>0 in (0,1]; (𝐸^′₂) 𝑞>0 in (0,1], 𝑞∈𝐿¹(0,1) and ∫¹₀(1/𝑝)(∫^𝑥₀𝑞𝑑𝑠)𝑑𝑥<∞;(𝐸^′₃)  the same as (𝐸₃).

In recent years, the study of such singular boundary value problems (1.1) has attracted the attention of several researchers [2–9]. In particular, if 𝑝(𝑡)=1,  𝑞(𝑡)=𝑡^−1/2 and 𝑓(𝑡,𝑦)=𝑦^3/2 problem (1.1) with Type 1 boundary conditions is known as Thomas-Fermi equation. Thomas [10] and Fermi [11] independently derived a boundary value problem for determining the electrical potential in an atom. The analysis leads to the nonlinear SBVP 𝑦^′′=𝑡^−1/2𝑦^3/2 with boundary conditions given by 𝑦(0)=1,  𝑦(𝑏₁)=0.

Chan and Hon [12] considered the generalized Thomas-Fermi equation: (𝑡^𝑏2𝑦^′)^′=𝑐𝑞(𝑡)𝑦^𝑒, 𝑦(0)=1, 𝑦(𝑎)=0 with 𝑞(𝑡)=𝑡^𝑏2+𝑑2, where 0≤𝑏₂<1, 𝑐>0, 𝑑₂>−2, 𝑒>1, which is doubly singular problem. Problem (1.1) with 𝑞(𝑡)=𝑝(𝑡)=𝑡^𝛼 where 𝛼=2 arises in the study of the distribution of heat sources in the human head [13] with 𝑓(𝑡,𝑦)=−𝛿𝑒^−𝜃𝑦, 𝜃>0, 𝛿>0.

There is a huge literature available on numerical methods for problem (1.1) with 𝑞(𝑡)=1, but very few numerical methods are available to tackle doubly singular boundary value problems. Reddien [14] studied the linear form of problem (1.1) and derived numerical methods for 𝑞(𝑡)∈𝐿²[0,1] which is stronger assumption than (𝐸^′₂).

Chawla and Shivakumar [15] established the existence as well as uniqueness of solution for problem (1.1) where 𝑞(𝑡)=𝑝(𝑡)=𝑡^𝛼. The existence and uniqueness for problem (1.1) have also been discussed by Dunninger and Kurtz [16] and Bobisud [1]. Later Pandey and Verma [17] extended the results on the existence uniqueness for problem (1.1) with Types 1 and 2 boundary conditions.

1.1. Adomian Decomposition Method (ADM)

In this subsection, we briefly describe standard ADM for nonlinear second-order equation.

Recently, many researchers [4, 18–26] have shown interest to the study of ADM for different scientific models. Adomian [18] asserted that the ADM provides an efficient and computationally suitable method for generating approximate series solution for a large class of differential equations. Let us consider nonlinear second-order ordinary differential equation 𝐿𝑦+𝑅𝑦+𝑁𝑦=𝑔(𝑡), (1.3) where 𝐿≡𝑑²/𝑑𝑡² is the second-order linear derivative operator, 𝑅 is the linear remainder operator, 𝑁 represents the nonlinear term, and 𝑔(𝑡) is a source term. The above equation can be rewritten as 𝐿𝑦=𝑔(𝑡)−𝑅𝑦−𝑁𝑦. (1.4) The inverse operator of 𝐿 is defined as 𝐿⁻¹(⋅)=∫^𝑡₀∫^𝑡₀(⋅)𝑑𝑡𝑑𝑡.

Operating the inverse linear operator 𝐿⁻¹(⋅) on both the sides of (1.4) yields 𝑦=𝑦(0)+𝑡𝑦^′(0)+𝐿⁻¹𝑔(𝑡)−𝐿⁻¹𝑅𝑦−𝐿⁻¹𝑁𝑦. (1.5) Next, we decompose the solution 𝑦 and the nonlinear function 𝑁𝑦 by an infinite series as 𝑦=^∞_𝑛=0𝑦_𝑛,𝑁𝑦=^∞_𝑛=0𝐴_𝑛, (1.6) where 𝐴_𝑛 are Adomian polynomials that can be constructed for various classes of nonlinear functions with the formula given by Adomian and Rach [19] 𝐴_𝑛=1𝑑𝑛!𝑑𝜆^𝑛𝑛𝑁^∞_𝑘=0𝑦_𝑘𝜆^𝑘_𝜆=0,𝑛=0,1,2,…. (1.7) Substituting the series (1.6) into (1.5), we obtain ^∞_𝑛=0𝑦_𝑛=𝑦(0)+𝑡𝑦^′(0)+𝐿⁻¹𝑔(𝑡)−𝐿⁻¹𝑅^∞_𝑛=0𝑦_𝑛−𝐿^∞−1_𝑛=0𝐴_𝑛. (1.8) From (1.8), the various components 𝑦_𝑛 of the solution 𝑦 can be determined by using the recursive relation 𝑦₀=𝑦(0)+𝑡𝑦^′(0)+𝐿⁻¹𝑔(𝑡),𝑦_𝑘+1=−𝐿⁻¹𝑅𝑦_𝑘−𝐿⁻¹𝐴_𝑘,𝑘≥0. (1.9) For numerical purpose, the 𝑛-term truncated series 𝜓_𝑛=∑^𝑛−1_𝑚=0𝑦_𝑚 may be used to give the approximate solution.

The ADM has been used to solve nonlinear boundary value problems (BVPs) for ordinary differential equations by several researchers [4, 20–26]. Solving nonlinear BVP by standard ADM is always a computationally involved task as it requires the computation of undetermined coefficients 𝑦^′(0) in (1.9) in a sequence of nonlinear equations which increases computational complexity.

It is important to note that the standard ADM (1.9) can not be applied directly to solve two-point boundary value problem (1.1) as the component 𝑦₀ in scheme (1.9) is not independent from undetermined coefficient. Many researchers [23–25] have proposed modified ADM to overcome the difficulty by setting undermined coefficient 𝑦^′(0)=𝑐, and then it will be determined by using second boundary condition satisfying 𝜓_𝑛(1)=𝑦(1). In this case, it requires additional computational work to solve nonlinear equation 𝜓_𝑛(1)=𝑦(1) for 𝑐 and 𝑐 may not be uniquely determined.

Benabidallah and Cherruault [24] have considered the following boundary value problem: 𝑦^′′+𝑝𝑦^′𝑦(𝑎)(𝑏)=𝑓𝑡,𝑦,𝑎<𝑡<𝑏,=𝛼,𝑦=𝛽. (1.10) The inverse operator 𝐿⁻¹=∫^𝑥_𝑎∫^𝑥_𝑏𝑑𝑥𝑑𝑥 was proposed, then unidentified constant 𝑦^′(𝑏) is set, as 𝑦^′(𝑏)=∑^∞_𝑛=0𝑐_𝑛, and each solution component 𝑦_𝑛 is obtained by the following scheme: 𝑦₀=𝑦(𝑎)+(𝑥−𝑎)𝑐₀(𝑏)𝑝(𝑏,𝑦+𝑦)_𝑛+1=𝑐_𝑛+1(𝑡−𝑎)−^𝑡_𝑎𝑝𝑦_𝑛𝑑𝑡+𝐿⁻¹𝐴_𝑛+𝐿⁻¹𝑝^′𝑦_𝑛,𝑛≥0. (1.11) In order to determine the unidentified constants 𝑐_𝑛, it is also required 𝑛-term approximate solution 𝜓_𝑛=∑^𝑛−1_𝑚=0𝑦_𝑚 satisfying all the boundary conditions. It is clearly that 𝜓_𝑛(𝑎)=𝑦(𝑎). For boundary condition 𝑥=𝑏, the following scheme was proposed 𝑦₀(𝑏)=𝑦(𝑏),𝑦_𝑛(𝑏)=0,𝑛≥1. (1.12) However, since each constant 𝑐_𝑛 is determined by solving (1.12), this scheme also require additional computational work.

In order to avoid solving such nonlinear algebraic equations for a two-point boundary value problems, Jang [26] introduced extended ADM for nonsingular problems with Dirichlet boundary condition. Khuri and Sayfy [27] applied a novel approach based on the mixed decomposition spline for solving singular problems arising in physiology. The method proposed by Ebaid [4] is based on the modification of Lesnic’s work [28]. Later, Duan and Rach [29] introduced special modified inverse linear integral operators for higher order boundary value problem.

To the best of our knowledge, no one has applied ADM to solve a class of doubly singular two-point boundary value problem (1.1) with Types 1 and 2 boundary conditions.

In this work, a new approach based on ADM and new integral operator are proposed for solving doubly singular two-point boundary value problems (1.1) with Types 1 and 2. To set up the modified scheme, we first consider DSBVP (1.1) and use all the boundary conditions and Adomian decomposition method to establish the recursion scheme. Thus a modified recursion scheme is developed which does not require the computation of undermined coefficients, whereas most of previous recursive schemes do require the computation of undermined coefficients (see [23–25]). In fact, the proposed recursion scheme is useful for solving problem (1.1) whether they are linear or nonlinear singular (nonsingular). The main advantage of the method is that it provides a direct scheme for solving the singular boundary value problem, that is, without linearization and discretization. Numerical results are presented to demonstrate the effectiveness of proposed recursive scheme. The symbolic and numerical computations have been performed using “MATHEMATICA” software.

The organization of the paper is as follows. Section 2 presents the inverse integral operator with modified Adomian decomposition method for doubly singular boundary value problems of (1.1) with Types 1 and 2 boundary conditions. In Section 3, the convergence of the method is discussed. In Section 4, we illustrate our method with numerical results along with graphical representation.

2. Modified Adomian Decomposition Method

In this section, we establish two modified recursive schemes for solving doubly singular boundary value problems (1.1) with Types 1 and 2 boundary conditions.

We again write (1.1) with Type 1 boundary conditions as: 𝑝(𝑡)𝑦^′(𝑡)^′(𝑡)𝑓(𝑡)𝑦(0)=𝑞𝑡,𝑦,0<𝑡≤1,=𝑎₁,𝑦(1)=𝑐₁, (2.1) which can be rewritten as 𝐿𝑦(𝑡)=𝑞(𝑡)𝑁𝑦(𝑡), (2.2) where 𝐿𝑦(𝑡)=(𝑝(𝑡)𝑦^′(𝑡))^′ is the linear differential operator to be inverted and 𝑁𝑦(𝑡)=𝑓(𝑡,𝑦(𝑡)) is an analytic nonlinear operator. It is also assumed that the solution of the problem (2.1) exits and unique. A twofold integral operator 𝐿⁻¹(⋅), regarded as the inverse operator of 𝐿(⋅), is proposed as 𝐿⁻¹(⋅)=^𝑡₀1𝑝(𝑠)¹_𝑠(⋅)𝑑𝑥𝑑𝑠. (2.3) To establish the recursive scheme, we operate 𝐿⁻¹(⋅) on the left hand side of (2.1) and use initial condition 𝑦(0)=𝑎₁; we have 𝐿⁻¹𝑝(𝑡)𝑦^′(𝑡)^′=^𝑡₀1𝑝(𝑠)¹_𝑠𝑝(𝑥)𝑦^′(𝑥)^′𝑑𝑥𝑑𝑠=𝑐₂^𝑡₀1𝑝(𝑠)𝑑𝑠−𝑦(𝑡)+𝑎₁,𝐿⁻¹𝑝(𝑡)𝑦^′(𝑡)^′=𝑐₂^𝑡₀1𝑝(𝑠)𝑑𝑠−𝑦(𝑡)+𝑎₁, (2.4)(2.5) where 𝑐₂=𝑝(1)𝑦^′(1).

We again operate the inverse operator 𝐿⁻¹(⋅) on both sides of (2.1) and use (2.5), yielding 𝑦(𝑡)=𝑎₁+𝑐₂^𝑡₀1𝑝(𝑠)𝑑𝑠−^𝑡₀1𝑝(𝑠)¹_𝑠𝑞(𝑡)𝑁𝑦(𝑡)𝑑𝑥𝑑𝑠. (2.6) For simplicity, we set ℎ(𝑡)=∫^𝑡₀(1/𝑝(𝑠))𝑑𝑠 and [𝐿⁻¹(⋅)]_𝑡=1=∫¹₀(1/𝑝(𝑠))∫¹_𝑠(⋅)𝑑𝑥𝑑𝑠.

Then (2.6) can be written as 𝑦(𝑡)=𝑎₁+𝑐₂ℎ(𝑡)−𝐿⁻¹𝑞(𝑡)𝑁𝑦(𝑡). (2.7) Using the boundary condition 𝑦(1)=𝑐₁ in (2.7), we get 𝑐₂=𝑐₁−𝑎₁ℎ(1)+1ℎ(1)𝐿⁻¹𝑞(𝑡)𝑁𝑦(𝑡)_𝑡=1. (2.8) By using 𝑐₂ into (2.7), we obtain 𝑦(𝑡)=𝑎₁+𝑐₁−𝑎₁ℎ(1)ℎ(𝑡)+ℎ(𝑡)ℎ(1)𝐿⁻¹𝑞(𝑡)𝑁𝑦(𝑡)_𝑡=1−𝐿⁻¹𝑞(𝑡)𝑁𝑦(𝑡). (2.9) It is important to note that the right hand side of (2.9) does not contain any unknown constants.

The equation (2.9) can be rewritten in the following operator equation form as 𝑦=𝑓+𝑀𝑦, (2.10) where 𝑓=𝑎₁+𝑐₁−𝑎₁ℎ(1)ℎ(𝑡)𝑦=ℎ(𝑡),𝑀ℎ(1)𝐿⁻¹𝑞(𝑡)𝑁𝑦(𝑡)_𝑡=1−𝐿⁻¹𝑞(𝑡)𝑁𝑦(𝑡). (2.11) Next, the solution 𝑦(𝑡) and the nonlinear function 𝑁𝑦(𝑡) are decomposed by infinite series of the form 𝑦(𝑡)=^∞_𝑛=0𝑦_𝑛(𝑡),𝑁𝑦(𝑡)=^∞_𝑛=0𝐴_𝑛𝑦₀(𝑡),𝑦₁(𝑡),…,𝑦_𝑛(𝑡), (2.12) respectively, where 𝐴_𝑛(𝑦₀(𝑡),𝑦₁(𝑡),…,𝑦_𝑛(𝑡)) are Adomian polynomials [19].

Substituting the series (2.12) into (2.9) gives ^∞_𝑛=0𝑦_𝑛(𝑡)=𝑎₁+𝑐₁−𝑎₁ℎ(1)ℎ(𝑡)+ℎ(𝑡)ℎ(1)𝐿⁻¹𝑞(𝑡)^∞_𝑛=0𝐴_𝑛𝑡=1−𝐿⁻¹𝑞(𝑡)^∞_𝑛=0𝐴_𝑛. (2.13) Using the above equation (2.13), the solution components 𝑦₀,𝑦₁,…,𝑦_𝑛 are determined by the following recursive scheme: 𝑦₀(𝑡)=𝑎₁+𝑐₁−𝑎₁ℎ(1)ℎ(𝑡),𝑦_𝑛+1(𝑡)=ℎ(𝑡)ℎ(1)𝐿⁻¹𝑞(𝑡)𝐴_𝑛𝑡=1−𝐿⁻¹𝑞(𝑡)𝐴_𝑛𝑦,𝑛≥0,or_𝑛+1(𝑡)=ℎ(𝑡)ℎ(1)¹₀1𝑝(𝑠)¹_𝑠𝑞(𝑡)𝐴_𝑛𝑑𝑥𝑑𝑠−^𝑡₀1𝑝(𝑠)¹_𝑠𝑞(𝑡)𝐴_𝑛𝑑𝑥𝑑𝑠,𝑛≥0. (2.14) We further modify the above algorithm (2.14) to get a more efficient and economic algorithm. To obtain modified algorithm, the zeroth component 𝑦₀ is divided into the sum of two parts, namely, 𝑓₀+𝑓₁, where 𝑓₀=𝑎₁,  𝑓₁=((𝑐₁−𝑎₁)/ℎ(1))ℎ(𝑡). The first part, 𝑓₀, is kept in 𝑦₀ and the rest part, 𝑓₁, is added to 𝑦₁. Thus, the zeroth component will be 𝑦₀=𝑎₁ in modified algorithm as 𝑦₀(𝑡)=𝑎₁,𝑦₁(𝑡)=𝑐₁−𝑎₁ℎ(1)ℎ(𝑡)+ℎ(𝑡)ℎ(1)𝐿⁻¹𝑞(𝑡)𝐴_0𝑡=1−𝐿⁻¹𝑞(𝑡)𝐴₀,⋮𝑦_𝑛+1(𝑡)=ℎ(𝑡)ℎ(1)𝐿⁻¹𝑞(𝑡)𝐴_𝑛𝑡=1−𝐿⁻¹𝑞(𝑡)𝐴_𝑛,𝑛≥0. (2.15) Note that new modified recursion scheme (2.15) does not require the computation of undermined coefficients, since this algorithm does not contain any undermined coefficient. This slight change plays a major role for minimizing the size of calculations.

Now, we again consider the same problem (1.1) with Type 2 boundary conditions: 𝑝(𝑡)𝑦^′(𝑡)^′(𝑡)𝑓(𝑡)𝑦=𝑞𝑡,𝑦,0<𝑡≤1,_′(0)=0,𝑎𝑦(1)+𝑏𝑦^′(1)=𝑐, (2.16) which can be rewritten as 𝐿𝑦(𝑡)=𝑞(𝑡)𝑁𝑦(𝑡), (2.17) where 𝐿𝑦(𝑡)=(𝑝(𝑡)𝑦^′(𝑡))^′ is the linear differential operator to be inverted and 𝑁𝑦(𝑡)=𝑓(𝑡,𝑦(𝑡)) is an analytic nonlinear operator. It is also assumed that the solution of problem (2.16) exits and is unique. In this case, we propose the following inverse integral operator: 𝐿⁻¹(⋅)=¹_𝑡1𝑝(𝑠)^𝑠₀(⋅)𝑑𝑥𝑑𝑠. (2.18) To setup an algorithm, we operate the inverse operator 𝐿⁻¹(⋅) on the left hand side of (2.16) and using the boundary condition 𝑦^′(0)=0, we have 𝐿⁻¹𝑝(𝑡)𝑦^′(𝑡)^′=¹_𝑡1𝑝(𝑠)^𝑠₀𝑝(𝑥)𝑦^′(𝑥)^′=𝑑𝑥𝑑𝑠_1𝑡𝑦^′(𝑠)𝐿𝑑𝑠,₋₁𝑝(𝑡)𝑦^′(𝑡)^′=𝑐₀−𝑦(𝑡), (2.19)(2.20) where 𝑐₀=𝑦(1).

By operating the inverse operator 𝐿⁻¹(⋅) on both sides of (2.16) and using (2.20) we have 𝑦(𝑡)=𝑐₀−¹_𝑡1𝑝(𝑠)^𝑠₀𝑞(𝑡)𝑁𝑦(𝑡)𝑑𝑥𝑑𝑠. (2.21) Thus (2.21) can be written as 𝑦(𝑡)=𝑐₀−𝐿⁻¹𝑞(𝑡)𝑁𝑦(𝑡). (2.22) Using boundary condition 𝑎𝑦(1)+𝑏𝑦^′(1)=𝑐 into (2.22), we obtain 𝑦(𝑡)=𝑐𝑎−𝑏𝑎𝑝(1)¹₀𝑞(𝑡)𝑁𝑦(𝑡)𝑑𝑡−𝐿⁻¹𝑞(𝑡)𝑁𝑦(𝑡). (2.23) The equation (2.23) may be rewritten in operator equation form 𝑦=𝑓+𝑀𝑦, (2.24) where 𝑐𝑓=𝑎𝑦𝑏,𝑀=−𝑎𝑝(1)¹₀𝑞(𝑡)𝑁𝑦(𝑡)𝑑𝑥−𝐿⁻¹𝑞(𝑡)𝑁𝑦(𝑡). (2.25) Substituting the series (2.12) into (2.23) gives ^∞_𝑛=0𝑦_𝑛(𝑡)=𝑐𝑎−𝑏𝑎𝑝(1)¹₀𝑞(𝑡)^∞_𝑛=0𝐴_𝑛𝑑𝑡−𝐿⁻¹𝑞(𝑡)^∞_𝑛=0𝐴_𝑛. (2.26) Upon matching both sides of (2.26), the solution components 𝑦₀,𝑦₁,…,𝑦_𝑛 can be determined by the following modified recursive scheme: 𝑦₀(𝑡)=𝑐𝑎,𝑦_𝑛+1(𝑡)𝑏=−𝑎𝑝(1)¹₀𝑞(𝑡)𝐴_𝑛𝑑𝑡−𝐿⁻¹𝑞(𝑡)𝐴_𝑛,𝑛≥0. (2.27) Thus, we have established two modified recursive schemes which give the complete determination of solution components 𝑦_𝑛, and hence the approximate series solution 𝜓_𝑛(𝑡) can be obtained for the singular boundary value problems (1.1) with Types 1 and 2 boundary conditions. For numerical purpose, the truncated 𝑛-term approximate series solution is given by 𝜓_𝑛(𝑡)=^𝑛−1_𝑚=0𝑦_𝑚(𝑡). (2.28) Note that the truncated 𝑛-term approximate series solution is obtained by adding 𝑛 solutions components 𝑦₀,𝑦₁,…,𝑦_𝑛−1, that is, 𝜓_𝑛=∑^𝑛−1_𝑚=0𝑦_𝑚.

3. Convergence of Method

In this section, we discuss the convergence analysis of modified ADM for doubly singular boundary value problem (1.1).

Many authors [30–33] established the convergence of ADM for differential and integral equations. The first proof of convergence of ADM for the general functional equation was given by Cherruault [31]. Cherruault and Adomian [32] proposed a new convergence proof of Adomian decomposition method for the general nonlinear functional equation based on the properties of convergent series. Recently, Hosseini and Nasabzadeh [33] introduced a simple technique to determine the rate of convergence of ADM for initial value problem.

Now, we discuss the convergence of doubly singular two-point boundary value problems (1.1) with Types 1 and 2 boundary conditions. To do so, note that (2.9) and (2.23) may be written in the operator equation form 𝑦=𝑓+𝑀𝑦, (3.1) where 𝑀 is an operator from a Banach space 𝐶[0,1] to 𝐶[0,1], 𝑓 is a given function in 𝐶[0,1], and we are looking for 𝑦∈𝐶[0,1] satisfying (3.1).

For Type 1 boundary conditions, we have 𝑓=𝑎₁+𝑐₁−𝑎₁ℎ(1)ℎ(𝑡),𝑀𝑦=ℎ(𝑡)ℎ(1)𝐿⁻¹𝑞(𝑡)𝑁𝑦(𝑡)_𝑡=1−𝐿⁻¹𝑞(𝑡)𝑁𝑦(𝑡). (3.2)(3.3) For Type 2 boundary conditions, 𝑓 and 𝑀 are given as 𝑐𝑓=𝑎,𝑀𝑦𝑏=−𝑎𝑝(1)¹₀𝑞(𝑡)𝑁𝑦(𝑡)𝑑𝑡−𝐿⁻¹𝑞(𝑡)𝑁𝑦(𝑡). (3.4)(3.5) Let the solution 𝑦 and the nonlinear function 𝑁𝑦 be considered as sum of series 𝑦=^∞_𝑛=0𝑦_𝑛,𝑁𝑦=^∞_𝑛=0𝐴_𝑛𝑦₀,𝑦₁,…,𝑦_𝑛, (3.6) where 𝐴_𝑛(𝑦₀,𝑦₁,…,𝑦_𝑛) are Adomian polynomials [19].

Substituting the series (3.6) into the operator equation (3.1) yields ^∞_𝑛=0𝑦_𝑛=𝑓+^∞_𝑛=0𝐵_𝑛𝑦₀,…,𝑦_𝑛,𝑦₀+𝑦₁+⋯+𝑦_𝑛+⋯=𝑓+𝐵₀𝑦₀+⋯+𝐵_𝑛−1𝑦₀,…,𝑦_𝑛−1+⋯, (3.7) where ^∞_𝑛=0𝐵_𝑛𝑦₀,…,𝑦_𝑛=ℎ(𝑡)ℎ(1)𝐿⁻¹𝑞(𝑡)^∞_𝑛=0𝐴_𝑛𝑡=1−𝐿⁻¹𝑞(𝑡)^∞_𝑛=0𝐴_𝑛,_∞𝑛=0𝐵_𝑛𝑦₀,…,𝑦_𝑛𝑏=−𝑎𝑝(1)¹₀𝑞(𝑡)^∞_𝑛=0𝐴_𝑛𝑑𝑡−𝐿⁻¹𝑞(𝑡)^∞_𝑛=0𝐴_𝑛, (3.8) for Types 1 and 2 boundary conditions, respectively.

Upon matching both sides of (3.7), we obtain the following scheme: 𝑦₀𝑦=𝑓,₁=𝐵₀𝑦₀,⋮𝑦_𝑛+1=𝐵_𝑛𝑦₀,…,𝑦_𝑛,⋮ (3.9)

Let 𝜓_𝑛=∑^𝑛−1_𝑚=0𝑦_𝑚 be approximate solution, then the modified ADM for (3.1) is equivalent to the following problem: 𝑦₀𝜓=𝑓,_𝑛+1=𝑦₀+𝑀𝜓_𝑛,𝑛≥0. (3.10) Finding the solution of (3.1) is equivalent to finding the sequence 𝜓_𝑛 such that 𝜓_𝑛=𝑦₀+𝑦₁+⋯+𝑦_𝑛−1 satisfies (3.10).

Theorem 3.1. Let 𝑀(𝑦) be the nonlinear operator defined by (3.3) or (3.5) which satisfies the Lipschitz condition ‖𝑀(𝜑)−𝑀(𝜉)‖≤𝛿‖𝜑−𝜉‖, for all 𝜑,𝜉∈𝐶[0,1] with Lipschitz constant 𝛿, 0≤𝛿<1. If ‖𝑦₀‖<∞, then there holds ‖𝑦_𝑘+1‖≤𝛿‖𝑦_𝑘‖, for all 𝑘∈ℕ∪{0} and the sequence {𝜓_𝑛} defined by (3.10) converges to 𝑦=∑^∞_𝑛=0𝑦_𝑛=lim_𝑛→∞𝜓_𝑛=𝜓.

Proof. Since 𝜓₁=𝑦₀+𝑦₁,𝜓₂=𝑦₀+𝑦₁+𝑦₂,…,𝜓_𝑛=𝑦₀+𝑦₁+𝑦₂+⋯+𝑦_𝑛,…, (3.11) we have 𝑦_𝑘+1=𝜓_𝑘+1−𝜓_𝑘, 𝑘=1,2,….
We now show that the sequence {𝜓_𝑛} is convergent. To prove this, it is sufficient to show that {𝜓_𝑛} is cauchy sequence in Banach space 𝐶[0,1]. Now using the nonlinearity of 𝑀(𝑦), we have 𝑦_𝑛+1=𝜓_𝑛+1−𝜓_𝑛=𝑀𝜓_𝑛−𝑀𝜓_𝑛−1≤𝛿𝜓_𝑛−𝜓_𝑛−1=𝛿𝑦_𝑛. (3.12) Thus we obtain 𝑦_𝑛+1≤𝛿𝑦_𝑛≤𝛿²𝑦_𝑛−1⋯≤𝛿^𝑛+1𝑦₀. (3.13) Now for every 𝑛,𝑚∈ℕ, and 𝑛≥𝑚, we have 𝜓_𝑛−𝜓_𝑚=𝜓_𝑛−𝜓_𝑛−1+𝜓_𝑛−1−𝜓_𝑛−2+⋯+𝜓_𝑚+1−𝜓_𝑚≤𝜓_𝑛−𝜓_𝑛−1+𝜓_𝑛−1−𝜓_𝑛−2+⋯+𝜓_𝑚+1−𝜓_𝑚≤𝛿^𝑛𝑦₀+𝛿^𝑛−1𝑦₀+⋯+𝛿^𝑚+1𝑦₀≤𝛿^𝑚+11+𝛿+𝛿²+⋯+𝛿^{𝑛−𝑚−1}𝑦₀. (3.14) As we know ∑^∞_𝑗=0𝛿^𝑗 is geometric series with common ration 𝛿, 0≤𝛿<1. Therefore, we have 𝜓_𝑛−𝜓_𝑚≤𝛿^𝑚+1𝑦1−𝛿⁰, (3.15) which converges to zero, that is, ‖𝜓_𝑛−𝜓_𝑚‖→0, as 𝑚→∞. This implies that there exits 𝜓 such that lim_𝑛→∞𝜓_𝑛=𝜓. But, we have 𝑦=∑^∞_𝑛=0𝑦_𝑛=lim_𝑛→∞𝜓_𝑛, that is, 𝑦=𝜓 which is solution of (3.1). Hence, we have convergence in norm.

4. Numerical Examples and Discussion

In this section, we consider four examples to demonstrate the effectiveness of proposed recursive schemes (2.15) and (2.27). The numerical results are compared with known results and maximum absolute error is also calculated. We also plot approximate and exact solutions to show how approximate solutions converge to exact solution.

Example 4.1. Consider linear two-point boundary value problem [6] 𝑡^𝛼𝑦^′′=𝛽𝑡^{𝛼+𝛽−2}𝛽𝑡^𝛽𝑦(0)(1)+𝛼+𝛽−1𝑦,0<𝑡≤1,=1,𝑦=𝑒, (4.1) with exact solution 𝑦(𝑡)=𝑒^𝑡𝛽. For every 𝛼>0 and 𝛽>0, (4.1) is called singular.

In particular, if 𝛼=0.5,  𝛽=1, that is, 𝑝(𝑡)=𝑡^0.5 and 𝑞(𝑡)=𝑡^−0.5 then problem (4.1) is called doubly singular boundary value problem. Applying modified recursive scheme (2.15) to (4.1), the solution components are obtained as 𝑦₀𝑦=1,₁=0.384948𝑡^0.5+𝑡+0.333333𝑡²,𝑦₂=−0.551338𝑡^0.5+0.128316𝑡^1.5+0.166667𝑡²+0.076989𝑡^2.5+0.155556𝑡³+0.023809𝑡⁴,𝑦₃=0.224316𝑡^0.5−0.183779𝑡^1.5−0.097435𝑡^2.5+0.011111𝑡³+0.015886𝑡^3.5+0.017460𝑡⁴+0.004277𝑡^4.5+0.007442𝑡⁵+0.000721𝑡⁶,𝑦₄=−0.077011𝑡^0.5+0.074771𝑡^1.5+0.026485𝑡^2.5−0.022142𝑡^3.5+0.000396𝑡⁴−0.004971𝑡^4.5+0.000881𝑡⁵+0.000655𝑡^5.5+0.0006418𝑡⁶+0.0001096𝑡^6.5+0.0001715𝑡⁷+0.000012𝑡⁸,⋮ (4.2) For numerical purpose, the approximate series solution is 𝜓₆=1−0.002003𝑡^0.5+𝑡+0.002065𝑡^1.5+0.5𝑡²+0.000603𝑡^2.5+0.166667𝑡³−0.0006957𝑡^3.5+0.041666𝑡⁴−0.0000457𝑡^4.5+0.0083333𝑡⁵+0.0000802𝑡^5.5+0.001388𝑡⁶+1.0684407×10⁻⁶𝑡^6.5+0.000198𝑡⁷−4.663181×10⁻⁶𝑡^7.5+0.000024𝑡⁸−3.885381×10⁻⁸𝑡^8.5+2.745287×10⁻⁶𝑡⁹+1.676726×10⁻⁷𝑡^9.5+2.664448×10⁻⁷𝑡¹⁰+1.5360215×10⁻⁸𝑡^10.5+2.063867×10⁻⁸𝑡¹¹+9.172396×10⁻¹⁰𝑡¹². (4.3) Now, we define error function as 𝐸_𝑛(𝑡)=𝜓_𝑛(𝑡)−𝑦(𝑡) and the maximum absolute errors as 𝐸_𝑛∞=max_0<𝑡≤1𝐸_𝑛(𝑡). (4.4) For different values of 𝛼 and 𝛽, the maximum absolute error ‖𝐸_𝑛‖_∞, for 𝑛=4,6, and 8, is shown in Tables 1 and 2, and we observe that when 𝑛 increases the error decreases. On the other hand, Chawla and Katti [6] applied 𝑀₃ method based on three evaluations of 𝑓 to solve the same problem in which a huge amount of computation work is needed to obtain the numerical solution. The results are compared with our approximate series solution in Table 3. Furthermore, the approximate series solutions 𝜓₂,  𝜓₄ and exact solution 𝑦 are plotted in Figures 1 and 2. It can be seen from the figures that series solution is very close to exact solution.

Table 1 

Maximum error of Example 4.1, when 𝛽=1.

Table 2 

Maximum error of Example 4.1, when 𝛽=5.

Table 3 

Maximum error of Example 4.1.

Figure 1 

Comparison of exact 𝑦 and approximate 𝜓2,  𝜓4 solutions of Example 4.1, 𝛼=0.50,  𝛽=1.

Figure 2 

Comparison of exact 𝑦 and approximate 𝜓2,  𝜓4 solutions of Example 4.1, 𝛼=0.50,  𝛽=5.

Example 4.2. Consider nonlinear singular two-point boundary value problem 𝑡^𝛼𝑦^′′=𝑡^{𝛼+𝛽−2}𝑒^𝑦𝑡^𝛽𝑒^𝑦𝑦(0)𝛽−𝛼−𝛽+1,0<𝑡≤1,=ln4(1)𝛽,𝑦=ln5, (4.5) with exact solution 𝑦(𝑡)=ln(𝛽/(4+𝑡^𝛽)) and for any real 𝛼>0 and 𝛽>0, this equation is called singular.

If 𝛼=0.25,  𝛽=1, that is, 𝑝(𝑡)=𝑡^0.25, 𝑞(𝑡)=𝑡^−0.75, then problem (4.5) is called doubly singular boundary value problem.

Now we apply recursive scheme (2.15) to obtain the solution components as 𝑦₀𝑦=−1.38629,₁=0.001856𝑡^0.75−0.25𝑡+0.025𝑡²,𝑦₂=−0.001605𝑡^0.75−0.000066𝑡^1.75+0.00625𝑡²+0.000042𝑡^2.75−0.004861𝑡³+0.000240𝑡⁴,𝑦₃=−0.000283𝑡^0.75+0.000057𝑡^1.75−2.461715×10⁻⁸𝑡^2.5−0.000030𝑡^2.75−0.000347𝑡³+4.475846×10⁻⁸𝑡^3.5−0.000011𝑡^3.75+0.000714𝑡⁴+8.882529×10⁻⁷𝑡^4.75−0.000103𝑡⁵+3.434065×10⁻⁶𝑡⁶,⋮ (4.6) To compare the results numerically, the 6-term truncated series solution is given by 𝜓₆=−1.38629−3.372935×10⁻⁷𝑥^0.75−0.25𝑡−1.108619×10⁻⁷𝑡^1.75+0.03125𝑡²+7.059037×10⁻⁹𝑡^2.5−3.299920×10⁻⁸𝑡^2.75−0.005208𝑡³+1.307442×10⁻¹¹𝑡^3.25−1.557418×10⁻⁸𝑡^3.5+1.289457×10⁻⁷𝑡^3.75+0.000976𝑡⁴−5.530529×10⁻¹¹𝑡^4.25+1.078097×10⁻⁸𝑡^4.5+3.090869×10⁻⁷𝑡^4.75−0.000195𝑡⁵−1.480660×10⁻¹¹𝑡^5.25+1.980188×10⁻⁹𝑡^5.5−3.668865×10⁻⁷𝑡^5.75+0.0000405𝑡⁶+2.567699×10⁻¹²𝑡^6.25−1.009792×10⁻⁹𝑡^6.5+1.315036×10⁻⁷𝑡^6.75−7.607217×10⁻⁶𝑡⁷+5.654798×10⁻¹¹𝑡^7.5−1.378289×10⁻⁸𝑡^7.75+8.830360×10⁻⁷𝑡⁸+4.154631×10⁻¹⁰𝑡^8.75−4.867214×10⁻⁸𝑡⁹+9.674748×10⁻¹⁰𝑡¹⁰. (4.7) Similarly, for 𝛼=0.5,  𝛽=1, that is, 𝑝(𝑡)=𝑡^0.5 and 𝑞(𝑡)=𝑡^−0.5, the solution components are obtained using the scheme (2.15) as 𝑦₀𝑦=−1.38629,₁=0.006023𝑡^0.5−0.25𝑡+0.020833𝑡²,𝑦₂=−0.005737𝑡^0.5−0.000501𝑡^1.5+0.010416𝑡²+0.000150𝑡^2.5−0.004513𝑡³+0.000186𝑡⁴,𝑦₃=−0.000302𝑡^0.5+0.000478𝑡^1.5−7.557899×10⁻⁷𝑡²−0.000093𝑡^2.5−0.000693𝑡³−0.000045𝑡^3.5+0.0007378𝑡⁴+2.788479×10⁻⁶𝑡^4.5−0.0000851𝑡⁵+2.348635×10⁻⁶𝑡⁶,⋮ (4.8) The maximum absolute error ‖𝐸_𝑛‖_∞, for 𝑛=4,6, and 8, is listed in Tables 4 and 5. As the value of 𝑛 increases, the maximum absolute error decreases. In addition, to verify how close is our approximate series solutions 𝜓₂,  𝜓₄ to exact solution 𝑦, we have plotted approximate solution and exact solution in Figures 3 and 4. It is shown from figures that the approximate series solution 𝜓₄ is very near to the exact solution 𝑦.

Table 4 

Maximum error for Example 4.2, when 𝛽=1.

Table 5 

Maximum error for Example 4.2, when 𝛽=5.

Figure 3 

Comparison of exact 𝑦 and approximate 𝜓2,  𝜓4 solutions of Example 4.2, when 𝛼=0.50,  𝛽=1.

Figure 4 

Comparison of exact 𝑦 and approximate 𝜓2,  𝜓4 solutions of Example 4.2, when 𝛼=0.50,  𝛽=5.

Example 4.3. The heat conduction model in human head problem Duggan and Goodman [13] is 𝑡²𝑦^′′=−𝑡²𝑒^−𝑦𝑦,0<𝑡≤1,_′(0)=0,2𝑦(1)+𝑦^′(1)=0. (4.9) We apply the recursive scheme (2.27) to (4.9), where 𝑝(𝑡)=𝑞(𝑡)=𝑡², 𝑓=𝑒^−𝑦, 𝑎=2, 𝑏=1, and 𝑐=0 then scheme (2.27) for (4.9) becomes as 𝑦₀𝑦=0,_𝑛+1=12¹₀𝑡²𝐴_𝑛𝑑𝑡+¹_𝑡𝑠^−2𝑠₀𝑥²𝐴_𝑛𝑑𝑥𝑑𝑠,𝑛≥0. (4.10)

We can easily determine the solution components from above scheme as 𝑦₀𝑦=0,₁=0.333333−0.166667𝑡²,𝑦₂=−0.086111+0.055555𝑡²−0.008333𝑡⁴,𝑦₃=0.032672−0.023611𝑡²+0.005555𝑡⁴−0.000529𝑡⁶,𝑦₄=−0.014584+0.011258𝑡²−0.003287𝑡⁴+0.000529𝑡⁶−0.0000373𝑡⁸,⋮ (4.11)

The truncated series solution 𝜓₈ is obtained by adding 𝑦₀,𝑦₁,…,𝑦₆, given by 𝜓₈=0.270259−0.127429𝑡²−0.0047677𝑡^4.−0.0002650𝑡⁶−4.7439600×10⁻⁶𝑡⁸−2.4024518×10⁻⁶𝑡¹⁰+2.2847637×10⁻⁷𝑡¹²−3.485344×10⁻⁸𝑡¹⁴+2.384089×10⁻⁹𝑡¹⁶−1.1874334×10⁻¹⁰𝑡¹⁸. (4.12) The comparison of numerical results of solutions with known solution is given in Table 6. It can be noted that the approximate solution 𝜓₈ is very close to known solution [13]. However, Duggan and Goodman [13] utilized the maximum principle based on lower and an upper solutions technique in which a lot of calculation is required as we need to solve sequence of linear differential equation analytically.

Table 6 

Comparison of numerical results with known results of Example 4.3.

Example 4.4. The nonlinear singular boundary problem considered by Chawla et al. [34] is 𝑦^′′+1𝑡𝑦^′=−𝑣𝑒^𝑦𝑦,0<𝑡<1,_′(0)=0,𝑦(1)=0, (4.13) with exact solution 𝑦(𝑡)=2ln((𝐵+1)/(𝐵𝑡²+1)), where 𝐵=(4−𝑣−2(4−2𝑣)^1/2)/𝑣, 𝑣=1.

In the above problem, we have 𝑝(𝑡)=𝑞(𝑡)=𝑡 and 𝑓=−𝑒^𝑦,  𝑎=1,  𝑏=0,  𝑐=0, by applying recursive scheme (2.27), we have following recursive scheme: 𝑦₀𝑦=0,_𝑛+1=01¹₀𝑡𝐴_𝑛𝑑𝑡+¹_𝑡𝑠^−1𝑠₀𝑥𝐴_𝑛𝑑𝑥𝑑𝑠,𝑛≥0. (4.14) We have successive solution components 𝑦_𝑛 as 𝑦₀𝑦=0,₁=0.25−0.25𝑡²,𝑦₂=0.046875−0.0625𝑡²+0.015625𝑡⁴,𝑦₃=0.0130208−0.0195313𝑡²+0.0078125𝑡⁴−0.00130208𝑡⁶,⋮ (4.15) and then the truncated series is obtained as 𝜓₆=0.316294−0.342438𝑡²+0.0289497𝑡⁴−0.003107𝑡⁶+0.0003280𝑥⁸−0.0000274𝑡¹⁰+1.2715657×10⁻⁶𝑡¹². (4.16) In order to check whether the approximate solution converges to exact solution 𝑦, we plot approximate solution 𝜓_𝑛, for 𝑛=2,4 and exact solution 𝑦 in Figure 5. It is clear from the figure that the approximate solution 𝜓₄ converges rapidly to exact solution. The maximum absolute error ‖𝐸_𝑛‖_∞, for 𝑛=5,8, and 10, is presented in Table 7. It can be noted that Chawla et al. [34] employed fourth-order finite difference method to obtain the numerical solution in which a lot of computation work is required. However, the maximum absolute error obtained was 1.28×10⁻⁵ using 𝑁=64.

Table 7 

Numerical results and maximum error of Example 4.4.

Figure 5 

Comparison of exact 𝑦 and approximate 𝜓2,  𝜓4 solutions of Example 4.4.

5. Conclusion

In this paper, we we have illustrated how the modified recursive schemes (2.15) and (2.27) can be used to solve a class of doubly singular two-point boundary value problems (1.1) with Types 1 and 2 boundary conditions. The accuracy of the numerical results indicates that the method is well suited for the solution of this type problem. The main advantage of this approach is that it provides a direct scheme to obtain approximate solutions, and we have also shown graphically that these approximate solutions are almost identical to the analytic solution. Another advantage of a modified recursion scheme is that this does not require the computation of undermined coefficients, whereas most of previous recursive schemes do require the computation of undermined coefficients (see [23–25]). The method provides a reliable technique which requires less work compared to the traditional techniques such as finite difference method, Cubic spline method, and standard ADM. The numerical results of the examples are presented and only a few terms are required to obtain accurate solutions. By comparing the results with other existing methods, it has been proved that proposed modified ADM is a more powerful method for solving the singular problems (1.1).