A new modification of the Adomian decomposition method for solving boundary value problems for higher order nonlinear differential equations

Applied Mathematics and Computation 218 (2011) 4090–4118

Contents lists available at SciVerse ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

A new modification of the Adomian decomposition method for solving boundary value problems for higher order nonlinear differential equations Jun-Sheng Duan a,⇑, Randolph Rach b a b

College of Science, Shanghai Institute of Technology, Fengxian District, Shanghai 201418, PR China 316 South Maple Street, Hartford, MI 49057-1225, USA

a r t i c l e

i n f o

Keywords: Adomian decomposition method (ADM) Adomian polynomials Nonlinear differential equations Boundary value problems (BVPs)

a b s t r a c t In this paper we propose a new modified recursion scheme for the resolution of multiorder and multi-point boundary value problems for nonlinear ordinary and partial differential equations by the Adomian decomposition method (ADM). Our new approach, including Duan’s convergence parameter, provides a significant computational advantage by allowing for the acceleration of convergence and expansion of the interval of convergence during calculations of the solution components for nonlinear boundary value problems, in particular for such cases when one of the boundary points lies outside the interval of convergence of the usual decomposition series. We utilize the boundary conditions to derive an integral equation before establishing the recursion scheme for the solution components. Thus we can derive a modified recursion scheme without any undetermined coefficients when computing successive solution components, whereas several prior recursion schemes have done so. This modification also avoids solving a sequence of nonlinear algebraic equations for the undetermined coefficients fraught with multiple roots, which is required to complete calculation of the solution by several prior modified recursion schemes using the ADM. Ó 2010 Elsevier Inc. All rights reserved.

1. Introduction We develop a new resolution method of boundary value problems (BVPs) for nonlinear ordinary and partial differential equations by the Adomian decomposition method (ADM) [1–9]. It is also well known that the ADM can provide approximate analytic solutions without using the Green function concept, which greatly facilitates analytic approximations and numerical computations. Several different resolution techniques for solving BVPs for nonlinear ordinary differential equations by using the ADM were considered by Adomian and Rach [10–12], Adomian [6], Wazwaz [8,13–19], Tatari and Dehghan [20], Dehghan and Tatari [21], Jang [22], Ebaid [23,24] and Al-Hayani [25]. In the process of resolution of a nonlinear BVP, solution of a sequence of nonlinear algebraic equations in the undetermined coefficients is often involved [8,13–21], which increases the complexity of the computations and the number of steps. One shortcoming of the method of undetermined coefficients for nonlinear BVPs is that we obtain a sequence of different sets of multiple roots and the number of roots increases more and more, because we need to solve a sequence of progressively higher order polynomials, or more complex transcendental equations, as the case may be. For a two-point BVP for second-order nonlinear differential equations, Adomian and Rach [10–12] proposed the double decomposition method in order to avoid solving such nonlinear algebraic equations, and Jang [22] and Ebaid [23]

⇑ Corresponding author. E-mail addresses: [email protected], [email protected] (J.-S. Duan), [email protected] (R. Rach). 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2011.09.037

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

4091

introduced different modified inverse linear operators. Al-Hayani [25] instead used Green’s functions to treat two-point high order BVPs. We consider multi-point nonlinear BVPs for higher order differential equations subject to mixed sets of multi-point and multi-order boundary conditions by the ADM. We utilize all of the boundary values to derive an integral equation before establishing the recursion scheme to calculate the solution components. Thus we develop a modified recursion scheme excluding all undetermined coefficients when computing successive solution components, whereas most previous recursion schemes do incorporate undetermined coefficients. Furthermore we can parametrize the new modified recursion scheme by embedding Duan’s convergence parameter [26,27] in order to achieve simple-to-integrate series, fast rate of convergence and extended region of convergence, in particular for such cases when one of the boundary points lies outside the interval of convergence of the usual decomposition series, see Example 8. We remark that the convergence of the Adomian series has already been proven by several investigators [7,28–32]. For example, Abdelrazec and Pelinovsky [32] have published a rigorous proof of convergence for the ADM under the aegis of the Cauchy–Kovalevskaya theorem for initial value problems. Furthermore prerequisites for existence and uniqueness of solutions for higher order BVPs have been provided in a comprehensive survey by Agarwal [33], and such higher order BVPs were first solved by Wazwaz [13] using the ADM. In point of fact the Adomian decomposition series is found to be a computationally advantageous rearrangement of the Banach-space analog of the Taylor expansion series about the initial solution component function, which permits solution via recursion. Furthermore convergence of the ADM is not limited to cases when only the fixed-point theorem applies, which is far too restrictive for most physical applications. 2. Description of the new approach for solving BVPs We illustrate the steps of our new method first with a second-order nonlinear differential equation

Lu ¼ Nu þ gðxÞ;

a 6 x 6 b;

ð1Þ

subject to the Dirichlet boundary conditions

uðaÞ ¼ a;

uðbÞ ¼ b;

d2 dx2

where LðÞ ¼ ðÞ is the linear differential operator to be inverted, Nu is an analytic nonlinear operator, and g(x) is the system input, which is a given continuous function. We also suppose that the solution of the BVP exists uniquely [33–36]. We take the inverse linear differential operator as

L1 ðÞ ¼

Z

x

Z

a

x

ðÞ dxdx;

ð2Þ

n

where n 2 [a, b] is a prescribed value. Then we have

L1 Lu ¼ uðxÞ  uðaÞ  u0 ðnÞðx  aÞ: Applying the operator L1 to both sides of Eq. (1) yields

uðxÞ  uðaÞ  u0 ðnÞðx  aÞ ¼ L1 Nu þ L1 gðxÞ:

ð3Þ

0

Let x = b in Eq. (3) and solve for u (n), then

u0 ðnÞ ¼

i 1 h uðbÞ  uðaÞ  ½L1 Nux¼b  ½L1 gðxÞx¼b ; ba

where

½L1 ðÞx¼b ¼

Z

b

a

Z

ð4Þ

x

ðÞ dxdx:

n

Substituting Eq. (4) into Eq. (3) yields

uðxÞ ¼ uðaÞ þ

uðbÞ  uðaÞ x  a 1 x  a 1 ðx  aÞ þ L1 g  ½L gx¼b þ L1 Nu  ½L Nux¼b : ba ba ba

ð5Þ

Thus the right hand side of Eq. (5) does not contain the undetermined coefficient u0 (n). Next, we decompose the solution and the nonlinearity

uðxÞ ¼

1 X

un ðxÞ;

Nu ¼

n¼0

1 X

An ;

ð6Þ

n¼0

where An = An(u0(x), u1(x), . . . , un(x)) are the Adomian polynomials, whose definitional formula

An ¼

n 1 X 1 d u k kk nN n! dk k¼0

!    

;

n P 0;

k¼0

was first published by Adomian and Rach [2] in 1983.

ð7Þ

4092

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

For convenient reference, we list the first five Adomian polynomials for the simple nonlinearity Nu = f(u),

A0 ¼ f ðu0 Þ; A1 ¼ f 0 ðu0 Þu1 ; A2 ¼ f 0 ðu0 Þu2 þ f 00 ðu0 Þ

u21 ; 2!

A3 ¼ f 0 ðu0 Þu3 þ f 00 ðu0 Þu1 u2 þ f 000 ðu0 Þ A4 ¼ f 0 ðu0 Þu4 þ f 00 ðu0 Þ

u31 ; 3!

 2  u2 u2 u2 u4 þ u1 u3 þ f 000 ðu0 Þ 1 þ f ð4Þ ðu0 Þ 1 : 2! 2! 4!

For a complicated nonlinearity such as Nu = f(x, u, u0 ), formula (7) still holds [27,37,38] since only the dependent variable u and its derivatives are parametrized by k, which is a grouping parameter of convenience that is later set equal to one. Sometimes k is also called the analytic parameter, but it is not a perturbation parameter and was not assumed to be small by Adomian and his collaborators. Examples 6 and 7 involve such complicated nonlinearities. Several algorithms [2,4,6–8,26,39–48] for symbolic programming have since been devised to efficiently generate the Adomian polynomials quickly and to high orders, For example, a convenient formula for the Adomian polynomials is Rach’s Rule, which reads (p. 16 in [4] and p. 51 in [6])

An ¼

n X

f ðkÞ ðu0 ÞCðk; nÞ;

n P 1;

ð8Þ

k¼1

where C(k, n) are the sums of all possible products of k components from u1, u2, . . . , unk+1, whose subscripts sum to n, divided by the factorial of the number of repeated subscripts [39]. New, efficient algorithms and subroutines in MATHEMATICA for rapid computer-generation of the Adomian polynomials to high orders are provided by Duan [26,37,46,47] and Duan and Guo [48]. From Eq. (5), the solution components are determined by the modified recursion scheme

uðbÞ  uðaÞ x  a 1 ðx  aÞ þ L1 g  ½L gx¼b ; ba ba x  a ½L1 An x¼b ; n P 0; ¼ L1 An  ba

u0 ¼ uðaÞ þ unþ1

ð9Þ ð10Þ

or explicitly,

unþ1 ¼

Z

x

a

Z

x

An dxdx 

n

xa ba

Z

Z

x

a

x

 An dxdx

n

;

n P 0:

x¼b

We emphasize that the right hand side of Eqs. (9) and (10) are independent of n, because such terms, which contain n, will cancel. In fact, if we denote ½1

h ðxÞ ¼

Z

hðxÞ dx;

½2

h ðxÞ ¼

Z

½1

h ðxÞ dx;

ð11Þ

where the right hand side denotes pure integration, i.e. without constants of integration, then we have

L1 hðxÞ 

x  a 1 ½L hðxÞx¼b ¼ ba

Z

x

a

Z

x

hðxÞ dxdx 

n

xa ba

Z a

x

Z

x n

 hðxÞ dxdx x¼b

i x  a h ½2 ½2 ½1 h ðxÞ  h ðaÞ  ðx  aÞh ðnÞ ¼ h ðxÞ  h ðaÞ  ðx  aÞh ðnÞ  x¼b ba i x  a h ½2 ½2 ½2 ½2 h ðbÞ  h ðaÞ : ¼ h ðxÞ  h ðaÞ  ba ½2

½2

½1

Thus such terms, involving the lower limit of integration n, where n 2 [a, b], cancel through addition. Therefore we have the freedom to choose n. Hence if Eq. (1) were to have a singularity at one of the end points, we would then choose the parameter n to be a different point within the specified interval. In some cases, we will instead use the following parametrized recursion scheme

u0 ¼ c; uðbÞ  uðaÞ x  a 1 x  a 1 ðx  aÞ þ L1 g  ½L gx¼b ; þL1 A0  ½L A0 x¼b ; ba ba ba x  a 1 ½L Anþ1 x¼b ; n P 0; ¼ L1 Anþ1  ba

ð12Þ

u1 ¼ c þ uðaÞ þ

ð13Þ

unþ2

ð14Þ

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

4093

where the constant c, which is Duan’s convergence parameter [26,27], is a predetermined parameter that is embedded to greatly simplify subsequent operations, i.e. always leading to simple-to-integrate series and can be adjusted in order to increase the domain of convergence [26,27], which are significant computational advantages for the ADM. In the following we consider a third-order nonlinear differential equation

Lu ¼ Nu þ gðxÞ;

ð15Þ

subject to the mixed set of Dirichlet and Neumann boundary conditions

uðx0 Þ ¼ a0 ;

u0 ðx1 Þ ¼ a1 ;

u0 ðx2 Þ ¼ a2 ;

x1 – x2 ;

ð16Þ

d3 dx3

where L ¼ is the linear differential operator to be inverted, Nu is an analytic nonlinear operator, and g(x) is the system input, which is a given continuous function. The field of x in Eq. (15) is min{x0, x1, x2} 6 x 6 max{x0, x1, x2}, where we include two-point and three-point BVPs. We take the inverse linear operator as

L1 ðÞ ¼

Z

x x0

Z

x

x1

Z

x

ðÞ dxdxdx;

ð17Þ

n

where n is a prescribed value in the specified interval. Then we have

h i 1 L1 Lu ¼ uðxÞ  uðx0 Þ  u0 ðx1 Þðx  x0 Þ  u00 ðnÞ ðx  x1 Þ2  ðx0  x1 Þ2 : 2 Applying the inverse operator L1 to both sides of Eq. (15) yields

h i 1 uðxÞ  uðx0 Þ  u0 ðx1 Þðx  x0 Þ  u00 ðnÞ ðx  x1 Þ2  ðx0  x1 Þ2 ¼ L1 ½Nu þ g: 2

ð18Þ

Differentiate Eq. (18), then let x = x2 and solve for u00 (n), hence

u00 ðnÞ ¼

u0 ðx2 Þ  u0 ðx1 Þ 1  x2  x1 x2  x1

Z

x2

x1

Z

x

½Nu þ g dxdx:

ð19Þ

n

Substituting Eq. (19) into Eq. (18) yields

i u0 ðx Þ  u0 ðx Þ 1h 2 1 ðx  x1 Þ2  ðx0  x1 Þ2 þ L1 g 2 x2  x1 Z Z Z Z 1 ðx  x1 Þ2  ðx0  x1 Þ2 x2 x 1 ðx  x1 Þ2  ðx0  x1 Þ2 x2 x gdxdx þ L1 Nu  Nudxdx:  2 2 x2  x1 x2  x1 x1 x1 n n

uðxÞ ¼ uðx0 Þ þ u0 ðx1 Þðx  x0 Þ þ

ð20Þ

Thus in Eq. (20), the three known boundary values u(x0), u0 (x1) and u0 (x2) are included and the undetermined coefficient was replaced. P P1 Next, we decompose the solution and the nonlinearity uðxÞ ¼ 1 n¼0 un ðxÞ; Nu ¼ n¼0 An , where An = An(u0(x), u1(x), . . . , un(x)) are the Adomian polynomials. From Eq. (20), the solution components are determined by the modified recursion scheme

1 u0 ðx2 Þ  u0 ðx1 Þ u0 ¼ uðx0 Þ þ u0 ðx1 Þðx  x0 Þ þ ½ðx  x1 Þ2  ðx0  x1 Þ2  2 x2  x1 2 2 Z x2 Z x 1 ðx  x Þ  ðx  x Þ 1 0 1 þ L1 g  gdxdx; 2 x2  x1 n x1 Z Z 1 ðx  x1 Þ2  ðx0  x1 Þ2 x2 x unþ1 ¼ L1 An  An dxdx; n P 0; 2 x2  x1 x1 n

ð21Þ ð22Þ

or by the parametrized recursion scheme

u0 ¼ c; 1 u0 ðx2 Þ  u0 ðx1 Þ u1 ¼ c þ uðx0 Þ þ u0 ðx1 Þðx  x0 Þ þ ½ðx  x1 Þ2  ðx0  x1 Þ2  2 x2  x1 2 2 Z x2 Z x 1 ðx  x1 Þ  ðx0  x1 Þ þ L1 g  gdxdx 2 x2  x1 x1 n Z Z 1 ðx  x1 Þ2  ðx0  x1 Þ2 x2 x A0 dxdx; þ L1 A0  2 x2  x1 n x1 Z Z 1 ðx  x1 Þ2  ðx0  x1 Þ2 x2 x unþ1 ¼ L1 An  An dxdx; n P 1: 2 x2  x1 x1 n

ð23Þ

ð24Þ ð25Þ

4094

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

By a similar manner as used in the previous second-order differential equation, we know that the right hand side of Eqs. (21), (22), (24) and (25) are independent of n. In practice, the n-term approximation

un ðxÞ ¼

n1 X

uk ;

ð26Þ

k¼0

can serve as an accurate approximate solution. For a general kth-order nonlinear differential equation k

d

k

dx

u ¼ Nu þ gðxÞ;

k P 2;

ð27Þ

subject to k boundary conditions

uðp0 Þ ðx0 Þ ¼ a0 ;

uðp1 Þ ðx1 Þ ¼ a1 ; . . . ; uðpk1 Þ ðxk1 Þ ¼ ak1 ;

ð28Þ

where x0, x1, . . . , xk1 are not all equal, 0 = p0 6 p1 6  6 pk1 6 k  1, pj 6 j for j = 1, 2, . . . , k  1, pi – pj if xi = xj. The field of x in Eq. (27) is min{x0, x1, . . . , xk1} 6 x 6 max{x0, x1, . . . , xk1}. If there are l different points in x0, x1, . . . , xk1 then we treat an l-point BVP. By p0 = 0 we presume that at least one boundary value equation is in the Dirichlet form. If the orders of the differentiations in Eq. (28) are different, i.e. pj = j for j = 1, 2, . . . , k  1, then we take the inverse linear operator as

L1 ðÞ ¼

Z

x

Z

x

Z

Z

x

x

k

ðÞdx ;

... x0

x1

x2

ð29Þ

xk1

and we can obtain a nonlinear integral equation involving all of the boundary conditions. For the sake of clarity, we consider the case p1 = 0, p2 = p3 =    = pk1 = 2, k P 3. For this case we take the inverse linear operator as

L1 ðÞ ¼

Z

x

Z

x

Z

x

Z

x

Z

Z

x

... x0

n

x2

n

n

x

k

ðÞ dx ;

ð30Þ

n

which will involve the boundary values u(x0) and u00 (x2) after applying it to the left hand side of Eq. (27), and where n is a prescribed value in the specified interval. We note that by replacing x0 with x1 and/or x2 with xj, for 3 6 j 6 k  1, in Eq. (30), resolution of Eq. (27) is still feasible by our new approach. Operating with L1 on both sides of Eq. (27) leads to

ðx  nÞ2  ðx0  nÞ2 u  uðx0 Þ  u0 ðnÞðx  x0 Þ  u00 ðx2 Þ 2 " # lþ3 lþ3 lþ1 k4 X ðx  nÞ  ðx  nÞ ðx ðx  nÞ2  ðx0  nÞ2 0 2  nÞ ðlþ3Þ u ðnÞ   ðl þ 3Þ! ðl þ 1Þ! 2 l¼0 ¼ L1 Nu þ L1 g;

ð31Þ 0

(3)

(4)

(k1)

where there are k  2 undetermined coefficients u (n), u (n), u (n), . . . , u (n), which can be determined by the remaining k  2 boundary values u(x1), u00 (x3), u00 (x4), . . . , u00 (xk1) through Eq. (31) and the equation obtained by twice differentiating Eq. (31). We note that in this case x0 – x1 and x2,x3, . . . , xk1 are different for any two. Therefore, through decompositions of u(x) and Nu, we can set up a modified recursion scheme without incorporating any undetermined coefficients. Compared with previous techniques using the ADM for solving BVPs, we first emphasize the choice of the lower limits in the operator L1 in order to utilize the given boundary values. Then we eliminate any remaining unknown parameters arising from the lower limits of L1 for the unused boundary values by algebraic manipulation, before we design the subsequent modified recursion scheme and hence its associated decomposition series. 3. Numerical examples In the following, we give eight examples including a third-order nonlinear differential equation with three specified boundary values in the form u(a), u(c), u(b), and a fourth-order nonlinear differential equation with four specified boundary values in the form u(a), u00 (a), u(b),u00 (b). Example 1. Consider the two-point BVP for the second-order nonlinear differential equation with an exponential nonlinearity [49]

u00 ðxÞ ¼ eu ; uð0Þ ¼ 0;

0 6 x 6 1; uð1Þ ¼ 0:

ð32Þ ð33Þ

4095

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

The exact solution for this problem is

  Cð2x  1Þ u ðxÞ ¼ 2 ln C sec  lnð2Þ; 4

ð34Þ

pffiffiffi where C satisfies C sec C4 ¼ 2, hence, to 16 significant figures, C = 1.336055694906108. The plot of u⁄(x) is shown in Fig. 1 for sake of reference, since the sequence of approximate solutions almost lie on top of one another and the exact solution. Thus, for comparison between the approximate solutions, we plot instead their error functions to highlight the rapid rate of convergence. Rx Rx Applying the inverse linear operator L1 ðÞ ¼ 0 0 ðÞ dxdx to both sides of Eq. (32) and using the boundary value at x = 0 yields

uðxÞ  u0 ð0Þx ¼ L1 eu :

ð35Þ

Let x = 1, then the above equation leads to

u0 ð0Þ ¼ ½L1 eu x¼1 ;

ð36Þ

where the boundary value at x = 1 is used. Substituting Eq. (36) into (35) gives the desired nonlinear integral equation for the solution

uðxÞ ¼ L1 eu  x½L1 eu x¼1 :

ð37Þ u

Next we decompose the solution u(x) and the nonlinearity e in the usual way, then the solution components are determined by the modified recursion scheme

u0 ¼ 0; un ¼ L1 An1  x½L1 An1 x¼1 ;

n ¼ 1; 2; . . . ;

where the Adomian polynomials for f(u) = eu are

  u2 A0 ¼ eu0 ; A1 ¼ eu0 u1 ; A2 ¼ eu0 u2 þ 1 ;    : 2 By computation we have

x x2 u1 ¼  þ ; 2 2 x x3 x4  u2 ¼ þ ; 24 12 24 x x3 x4 x5 x6 þ u3 ¼  þ  þ ; 160 144 96 60 180 ...: We examine the convergence of the n-term approximation un ðxÞ ¼

Pn1

k¼0 uk .

We denote the error functions as

En ðxÞ ¼ un ðxÞ  u ðxÞ;

u x 0.2

0.4

0.6

0.02

0.04

0.06

0.08

0.10

0.12 Fig. 1. The exact solution u⁄(x).

0.8

1.0

x

4096

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

En x 0.002

0.2

0.4

0.6

0.8

1.0

x

0.002 0.004 0.006 0.008 0.010

Fig. 2. Error functions E2(x) (solid line), E3(x) (dot line) and E4(x) (dashed line).

MEn 0.01

0.001

10

4

10

5

3

4

5

6

7

n

Fig. 3. Logarithmic plots of maximal errors MEn versus n (n = 2, 3, 4, 5, 6 and 7).

and the maximal errors as

MEn ¼ max jEn ðxÞj; 06x61

which are computed by a native command in MATHEMATICA, NMaximize. In Fig. 2 we plot the curves of the error functions En(x) for n = 2, 3 and 4. In Fig. 3 we display the logarithmic plots of the maximal errors MEn versus n for n = 2,3, . . . ,7. The data points lie almost on a straight line, which means that the maximal errors decrease approximately at an exponential rate. Example 2. Consider the two-point BVP for the third-order nonlinear differential equation with a radical nonlinearity

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u000 ðxÞ ¼  1  u2 ðxÞ; 0

u ð0Þ ¼ 1;

uð0Þ ¼ 0;

0 6 x 6 p=2;

ð38Þ

uðp=2Þ ¼ 1:

ð39Þ

Applying the inverse linear operator

L1 ðÞ ¼

Z

x

Z

0

Z

x 0

x

ðÞ dxdxdx

ð40Þ

0

to both sides of Eq. (38) yields

uðxÞ  x  u00 ð0Þ

x2 ¼ L1 NuðxÞ; 2

ð41Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where NuðxÞ ¼ 1  u2 ðxÞ and the boundary values at x = 0 are used. Let x ¼ p2 in Eq. (41) then we obtain

u00 ð0Þ ¼

8

p2



4

þ

8 h

p p2

L1 NuðxÞ

i x¼p=2

:

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

4097

Inserting this into Eq. (41), we have

uðxÞ ¼ x þ

ð4  2pÞx2

p2

þ

4x2 h

p2

i L1 NuðxÞ

x¼p=2

 L1 NuðxÞ:

ð42Þ

Next we design the modified recursion scheme so as to achieve a simple-to-integrate series, where c is Duan’s convergence parameter [26,27],

u0 ¼ c; u1 ¼ c þ x þ un ¼

4x2

p2

ð4  2pÞx2

p2

þ

4x2

p2

½L1 A0 x¼p=2  L1 A0 ;

½L1 An1 x¼p=2  L1 An1 ;

n P 2;

where the Adomian polynomials for this radical nonlinearity are

A0 ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  u20 ;

u0 u1 A1 ¼  qffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 1  u20

u2 u2 u21 u0 u2 A2 ¼   0 1 3=2  qffiffiffiffiffiffiffiffiffiffiffiffiffiffi  qffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; . . . ; 2 2 2 1  u0 1  u20 2 1  u0

which were first calculated by Adomian et al. [50] in 1985. We remark that here c is not an undetermined parameter. For example, it can take its value between the two boundary values of 0 and 1. Further computation yields

1 pffiffiffiffiffiffiffiffiffiffiffiffiffi2 2 1 pffiffiffiffiffiffiffiffiffiffiffiffiffi2 3 1  c px  1c x ; 6 p p 12  c p 4 x2 1 cx6 1 1 1 2 2 u2 ¼  cpx5  c p x2 þ c px þ þ pffiffiffiffiffiffiffiffiffiffiffiffiffi  120 12 11520 720 720 1  c2  1 c2 x3 cx4 cx5 cx5  ; c p 2 x2   þ þ 160 6 24 15p2 30p  1 p 4 x2 cp4 x2 23c2 p4 x2 c3 p4 x2 p5 x2 3c2 p5 x2  þ þ   þ u3 ¼ 2 1c 53760 11520 1612800 11520 107520 358400

u1 ¼ c þ x þ

4x2 2



2x2

þ

1 c2 px5 1 3 5 c2 p2 x5 cx6 c3 x6 p x6 c p x5  c px  þ þ  þ 720 7200 720 9600 720 720 1440  x7 c 2 x7 x7 x8 c 2 x8 x8 c 2 x8 þ þ   þ þ  2 2 630 5040 630p 504p 5040p 1008p 10080p  1 p x2 1 1 2 2 p 2 x2 1 c p x2  c px  c p 2 x2  þ þ þ 3=2 2 120 24 160 840 1120 ð1  c Þ 



p 3 x2 3360



p7 x2 11612160

þ

13c2 p7 x2 c 4 p 7 x2 c2 x3 cx4 x5  þ  þ 58060800 7257600 12 24 120

cx cx5 c 2 p 4 x5 c4 p4 x5 x6 x6 4x7 4x7 þ   þ  þ þ 2 2 4 15p 30p 691200 691200 30p 60p 105p 105p3 5



x7 p2 x7 c2 p2 x7 px8 11c2 px8 c 4 p x8 x9 þ   þ  þ 105p2 60480 60480 24192 241920 241920 36288  11c2 x9 c 4 x9  þ ; 362880 362880 þ

...: For c = 0, the even-numbered solution components u2k = 0 for all k = 0, 1, 2, . . .. In Fig. 4 we plot the curves of the error functions En(x) = un(x)  u⁄(x) for n = 2, 4 and 6, where u⁄(x) = sin (x) is the exact solution. In Table 1 we list the maximal errors MEn = max06x6p/2jEn(x)j for n = 2, 4, 6, 8, 10 and 12. For c = 0.5, the maximal errors MEn = max06x6p/2jEn(x)j, for n = 2, 3, 4, 5, 6 and 7, are listed in Table 2. We see that the value of c can affect the rate of convergence for the decomposition series solution of BVPs. Thus we can use this phenomenon to design more efficient algorithms by adjusting the value of the convergence parameter c. Example 3. Consider the three-point BVP for the third-order nonlinear differential equation with a hyperbolic sine nonlinearity

4098

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

En x

0.020

0.015

0.010

0.005

0.5

1.0

x

1.5

Fig. 4. Error functions E2(x) (solid line), E4(x) (dot line) and E6(x) (dashed line).

Table 1 For c = 0, maximal errors MEn for n = 2, 4, 6, 8, 10 and 12. n

2

4

6

8

10

12

MEn

0.0240463

0.0041320

0.0018206

0.0010355

0.0006737

0.0004760

Table 2 For c = 0.5, maximal errors MEn for n = 2, 3, 4, 5, 6 and 7. n

2

3

4

5

6

7

MEn

0.0119453

0.0061582

0.0015818

0.0010641

0.0006377

0.0004828

u000 ðxÞ ¼ 1 þ x sinhðuÞ; uð0Þ ¼ 0;

uð0:25Þ ¼ 1;

0 6 x 6 1; uð1Þ ¼ 0:

ð43Þ ð44Þ

Applying the inverse linear operator in Eq. (40) to both sides of Eq. (43) gives

uðxÞ  u0 ð0Þx  u00 ð0Þ

x2 x3 ¼ þ L1 x sinhðuÞ; 2 6

ð45Þ

where the boundary value at x = 0 is used. We evaluate the above equation for x = 0.25 and 1, respectively, and then obtain through algebraic manipulation

i i 43 16 h 1 1h  L x sinhðuÞ þ L1 x sinhðuÞ ; x¼1=4 x¼1 8 3 3 i i 133 32 h 1 8h u00 ð0Þ ¼   L1 x sinhðuÞ : þ L x sinhðuÞ x¼1=4 x¼1 12 3 3

u0 ð0Þ ¼

ð46Þ ð47Þ

Substituting Eqs. (46) and (47) into (45) gives the desired nonlinear integral equation for the solution

uðxÞ ¼

h i h i 43 133 2 x3 16 2 1 x x þ þ ðx  xÞ L1 x sinhðuÞ þ ðx  4x2 Þ L1 x sinhðuÞ þ L1 x sinhðuÞ: x¼1=4 x¼1 8 24 3 3 6

We design the modified recursion scheme

u0 ¼ c; 43 133 2 x3 16 2 1 x x þ þ ðx  xÞ½L1 xA0 x¼1=4 þ ðx  4x2 Þ½L1 xA0 x¼1 þ L1 xA0 ; 8 24 3 3 6 16 2 1 1 1 2 un ¼ ðx  xÞ½L xAn1 x¼1=4 þ ðx  4x Þ½L xAn1 x¼1 þ L1 xAn1 ; n P 2; 3 3

u1 ¼ c þ

ð48Þ

4099

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

where the Adomian polynomials for the hyperbolic sine nonlinearity are

A0 ¼ sinhðu0 Þ; A1 ¼ u1 coshðu0 Þ; 1 A2 ¼ u21 sinhðu0 Þ þ u2 coshðu0 Þ; 2 1 A3 ¼ u2 u1 sinhðu0 Þ þ u31 coshðu0 Þ þ u3 coshðu0 Þ; 6 ...; which were first calculated by Adomian and Rach [51] in 1984. By computation we have

u1 ¼

1 4 ðx2  xÞ sinhðcÞ 1 x3 133x2 43x x sinhðcÞ þ þ ðx  4x2 Þ sinhðcÞ  c þ  ; þ 24 1152 72 8 6 24

u2 ¼

x8 sinhðcÞ coshðcÞ x7 coshðcÞ 133x6 coshðcÞ 1 4 þ   cx coshðcÞ 8064 1260 2880 24 7x6 sinhðcÞ coshðcÞ 43 5 x5 sinhðcÞ coshðcÞ þ x coshðcÞ þ 15360 480 4608   16 2 c coshðcÞ 2099 coshðcÞ 271 sinhðcÞ coshðcÞ þ ðx  xÞ  þ þ 3 6144 27525120 2642411520   1 1 99 coshðcÞ 37 sinhðcÞ coshðcÞ  ; þ ðx  4x2 Þ  c coshðcÞ þ 3 24 2240 322560



...: Since the exact solution u⁄(x) of this problem is unknown, we instead investigate the error remainder function, which is a measure of how well the approximation satisfies the original nonlinear differential equation,

ERn ðxÞ ¼ u000 n ðxÞ  1  x sinhðun ðxÞÞ;

0 6 x 6 1;

ð49Þ

which in general is defined as

ERn ðxÞ ¼ Lun ðxÞ þ Run ðxÞ þ Nun ðxÞ  gðxÞ;

a 6 x 6 b;

for the nonlinear differential equation

LuðxÞ þ RuðxÞ þ NuðxÞ ¼ gðxÞ; where the n-term approximation un(x) is used in place of the solution u(x) in order to examine the convergence of the decomposition series solution over the specified interval, since the limit of un(x) equals u⁄(x). For c = 1, the boundary value at x = 0.25, we plot the functions ERn(x) for n = 2, 3, 4 in Fig. 5 and the functions ERn(x) for n = 5, 6, 7 in Fig. 6, which demonstrates a remarkable accuracy for such low orders of approximation by the ADM as we expect the limit of the error function to rapidly approach zero.

ERn x

1.0

0.5

0.2

0.4

0.6

0.8

1.0

x

Fig. 5. For c = 1, curves of ER2(x) (solid line), ER3(x) (dot line) and ER4(x) (dashed line).

4100

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

ERn x 0.01

0.2

0.4

0.6

0.8

1.0

x

0.01

0.02

0.03

Fig. 6. For c = 1, curves of ER5(x) (solid line), ER6(x) (dot line) and ER7(x) (dashed line).

Example 4. Consider the four-point BVP for the fourth-order nonlinear differential equation with a product nonlinearity [20]

uð4Þ ðxÞ þ uðxÞu0 ðxÞ  4x7  24 ¼ 0; u000 ð0:25Þ ¼ 6;

uð0Þ ¼ 0;

0 6 x 6 1;

u00 ð0:5Þ ¼ 3;

ð50Þ

uð1Þ ¼ 1:

ð51Þ

The exact solution of this problem is u⁄(x) = x4. We take

L1 ðÞ ¼

Z

x

0

Z

x 0

Z

x 0:5

Z

x

ðÞ dxdxdxdx

ð52Þ

0:25

as the inverse linear operator. Then we have for the boundary values at x = 0, 0.25 and 0.5,

L1 uð4Þ ðxÞ ¼ uðxÞ  u0 ð0Þx  x3 : In addition

L1 ð4x7 þ 24Þ ¼ 

247x2 786433x3 x11  þ x4 þ : 4718592 786432 1980

Thus operating with L1 on both sides of Eq. (50) gives

uðxÞ  u0 ð0Þx þ

247x2 x3 x11 þ  x4  þ L1 uu0 ¼ 0: 4718592 786432 1980

ð53Þ

Using the boundary value at x = 1, we derive

u0 ð0Þ ¼ 

117157 þ ½L1 uu0 x¼1 : 259522560

ð54Þ

Substituting Eq. (54) into (53) leads to the desired nonlinear integral equation for the solution

uðxÞ ¼ 

h i 117157x 247x2 x3 x11  :  þ x4 þ  L1 uðxÞu0 ðxÞ þ x L1 uðxÞu0 ðxÞ x¼1 259522560 4718592 786432 1980

We decompose the solution and the nonlinearity u ¼ als are

P1

n¼0 un ;

A0 ¼ u0 u00 ; A1 ¼ u0 u01 þ u1 u00 ; A2 ¼ u0 u02 þ u1 u01 þ u2 u00 ; A3 ¼ u0 u03 þ u1 u02 þ u2 u01 þ u3 u00 ; ...: By the modified recursion scheme

117157x 247x2 x3 x11   þ x4 þ ; 259522560 4718592 786432 1980 1 1 un ¼ L An1 þ x½L An1 x¼1 ; n ¼ 1; 2; . . . ;

u0 ¼ 

uu0 ¼

P1

n¼0 An ,

respectively, where the Adomian polynomi-

4101

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

we obtain the successive solution components, among which we list u1 and u2 as follows,

u1 ¼ 9:241884516  1012 x25  1:031557404  107 x18 þ 1:574031684  1013 x17 þ 7:86828457  1012 x16 þ 8:351515442  1011 x15  0:0005050505051x11 þ 1:766063549  109 x10 þ 1:038613547  107 x9 þ 1:343549849  106 x8  9:257546841  1012 x7  1:96922983  1010 x6  1:698263226  109 x5 þ 1:197071763  106 x3 þ 0:00005184309329x2 þ 0:0004506661406x; u2 ¼ 8:512270566820489  1020 x39 þ 2:0611808907803163  1015 x32  3:38332271178819  1021 x31  1:8299088362943092  1019 x30  2:114188144007494  1018 x29 þ 2:595882268568207  1011 x25  9:721045476235172  1017 x24  6:185141837504404  1015 x23  8:751553594634062  1014 x22 þ 6:59731186730981  1019 x21 þ 1:5698793638472888  1017 x20 þ 1:5400931909423285  1016 x19 þ 1:0315574041064238  107 x18  7:384436702019146  1013 x17  4:657207413215791  1011 x16  6:590311932698167  1010 x15 þ 7:875438264780266  1015 x14 þ 1:8306469783644867  1013 x13 þ 1:7460492288019172  1012 x12  1:6814917783267445  1017 x11  1:6625999182262368  109 x10  1:0286327883752899  107 x9  1:3412678929520298  106 x8 þ 1:8224924175027626  1011 x7 þ 3:916191589406249  1010 x6 þ 3:3907580800253106  109 x5 þ 7:322810785423231  108 x3 þ 5:013353855662456  107 x2 þ 7:649524731979151  107 x: In Figs. 7 and 8, we plot the error functions En(x) = un(x)  u⁄(x) for n = 1 and 2, respectively, with different vertical scales to better highlight the rapid rate of convergence of the decomposition series. In Fig. 9, we display the logarithmic plots of the maximal errors MEn = max06x61jEn(x)j versus n for n = 1, 2, 3, 4, where the four points lie almost in a straight line, which indicates that the maximal errors decrease approximately at an exponential rate. In contrast to our new approach, we note that the recursion scheme in [20] contains three undetermined coefficients. Example 5. Consider the two-point BVP for the fourth-order nonlinear differential equation with an exponential nonlinearity [19]:

E1 x 0.2

0.4

0.6

0.00005 0.00010 0.00015 0.00020 0.00025 0.00030 0.00035 Fig. 7. Curve of E1(x).

0.8

1.0

x

4102

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

E2 x 0.2 1. 10

7

2. 10

7

3. 10

7

4. 10

7

5. 10

7

6. 10

7

7. 10

7

0.4

0.6

0.8

1.0

x

Fig. 8. Curve of E2(x).

MEn

10

5

10

7

10

9

1.5

2.0

2.5

3.0

3.5

4.0

n

Fig. 9. Logarithmic plots of the maximal errors MEn versus n for n = 1, 2, 3 and 4.

uð4Þ ðxÞ ¼ 6e4uðxÞ ;

0 6 x 6 4  e; 1 u00 ð0Þ ¼  2 ; uð4  eÞ ¼ lnð4Þ; e

uð0Þ ¼ 1;

ð55Þ u00 ð4  eÞ ¼ 

1 : 16

ð56Þ

The exact solution for this problem is u⁄(x) = ln(e + x). We take

L1 ðÞ ¼

Z 0

x

Z

x 0

Z

x 0

Z

x

ðÞ dxdxdxdx

ð57Þ

0

as the inverse linear operator. Then we have for the boundary values at x = 0,

uðxÞ  1  u0 ð0Þx þ

x2 x3  u000 ð0Þ ¼ 6L1 e4u : 2 2e 6

ð58Þ

Let x = 4  e in Eq. (58), we obtain

ð4  eÞu0 ð0Þ þ

ð4  eÞ3 000 ð4  eÞ2 þ 6½L1 e4u x¼4e : u ð0Þ ¼ lnð4Þ  1 þ 6 2e2

ð59Þ

Differentiating Eq. (58) two times, then let x = 4  e, we obtain

ð4  eÞu000 ð0Þ ¼

1 1 þ6  e2 16

Z 0

4e

Z

x

e4u dxdx: 0

ð60Þ

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

4103

Solving for u0 (0) and u000 (0) from Eqs. (59) and (60) yields

ð4  eÞð32 þ e2 Þ lnð4Þ  1  ð4  eÞ þ 96e2 4e Z 4e Z x 4þe 6 þ e4u dxdx: u000 ð0Þ ¼ 2 16e 4e 0 0

u0 ð0Þ ¼

Z

4e 0

Z

x

e4u dxdx þ

0

6 ½L1 e4u x¼4e ; 4e

ð61Þ ð62Þ

Substituting Eqs. (61) and (62) into (58) gives the nonlinear integral equation

  Z 4e Z x ð4  eÞð32 þ e2 Þ lnð4Þ  1 x2 ð4 þ eÞx3 6x  ½L1 e4u x¼4e þ þ  ð4  eÞx e4u dxdx þ 96e2 4e 4e 2e2 96e2 0 0 Z 4e Z x x3 þ e4u dxdx  6L1 e4u ; 4e 0 0

uðxÞ ¼ 1 þ x

ð63Þ

where all four of the boundary values are used. Next we apply the decomposition of the solution and the nonlinearity in the usual way. Then we design the modified recursion scheme

u0 ¼ 1;   Z 4e Z x ð4  eÞð32 þ e2 Þ lnð4Þ  1 x2 ð4 þ eÞx3 6x  ½L1 A0 x¼4e þ þ  ð4  eÞx A0 dxdx þ u1 ¼ x 96e2 4e 4e 2e2 96e2 0 0 Z 4e Z x x3 þ A0 dxdx  6L1 A0 ; 4e 0 0 Z 4e Z x Z 4e Z x 6x x3 ½L1 An1 x¼4e þ An1 dxdx þ An1 dxdx  6L1 An1 ; n P 2; un ¼ ð4  eÞx 4e 4e 0 0 0 0 where the Adomian polynomials are

A0 ¼ e4u0 ; A1 ¼ 4e4u0 u1 ; A2 ¼ 4e4u0 ð2u21  u2 Þ:   8 3 A3 ¼ 4e4u0 u1  4u1 u2 þ u3 ; 3 ...: Then we can calculate the first several solution components,

u1 ¼ 0:3629183x  0:0676676x2 þ 0:0212088x3  0:00457891x4 ; u2 ¼ 0:00746686x  0:00658100x3 þ 0:00132942x5  0:000082625x6 þ 0:000011099x7  1:1980808  106 x8 ; u3 ¼ 0:003358787x þ 0:002575707x3 þ 0:000027352x5  0:000321646x6 þ 0:00004796097x7  0:00001045196x8 þ 1:9939469  106 x9  1:9376435  107 x10 þ 2:2175865  108 x11  1:5958993  109 x12 ; ...: In Figs. 10 and 11, we plot the error functions En(x) = un(x)  u⁄(x) for n = 2, 3, 4 and for n = 5, 6, 7, respectively, with different vertical scales to better highlight the rapid rate of convergence of the decomposition series. In Fig. 12, we display the logarithmic plots of the maximal errors MEn = max06x64  ejEn(x)j versus n for n = 2, 3, 4, 5, 6 and 7. Fig. 12 shows that the maximal errors ME2 through ME7 decrease more quickly than an exponential rate. We note that the recursion scheme in [19] contains two undetermined coefficients. In the sequel we consider two BVPs occurring in heat transfer analysis and design simulations [52,53].

4104

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

En x 0.0010 0.0005

0.2

0.4

0.6

0.8

1.0

x

1.2

0.0005 0.0010 0.0015 0.0020 Fig. 10. Error functions E2(x) (solid line), E3(x) (dot line) and E4(x) (dashed line).

En x 0.0001

0.00008

0.00006

0.00004

0.00002

0.2

0.4

0.6

0.8

1.0

x

1.2

Fig. 11. Error functions E5(x) (solid line), E6(x) (dot line) and E7(x) (dashed line).

MEn

0.001 5 10

4

1 10

4

5 10

5

1 10

5

5 10

6

3

4

5

6

7

n

Fig. 12. Logarithmic plots of the maximal errors MEn versus n for n = 2, 3, 4, 5, 6 and 7.

Example 6. Consider the nonlinear second-order homogeneous ordinary differential equation with a variable coefficient [53]

h00 ðrÞ þ

e m2 r hðrÞ ðh0 ðrÞÞ2  ¼ 0; r b 1 þ ehðrÞ 1 þ ehðrÞ

ð64Þ

4105

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

where h(r) has the domain of definition r 2 [rb, 1] and 0 < rb < 1, and subject to a mixed set of inhomogeneous Dirichlet and homogeneous Neumann boundary conditions

h0 ð1Þ ¼ 0:

hðr b Þ ¼ 1;

The physical variables and parameters are h, r, , m and rb, which represent the dimensionless temperature, radius, thermal conductivity parameter, fin parameter and fin base radius, respectively. The interested reader is referred to [53] for further details of the derivation and design limitations of this engineering model. For convenience, we rewrite this nonlinear differential equation in Adomian’s operator-theoretic form

LhðrÞ þ NhðrÞ ¼ 0;

ð65Þ

2

where LðÞ ¼ drd 2 ðÞ, and

NhðrÞ ¼

e m2 r hðrÞ ðh0 ðrÞÞ2  : r b 1 þ ehðrÞ 1 þ ehðrÞ

We define

L1 ðÞ ¼

Z

r

Z

rb

r

ðÞ drdr;

1

then we have after applying the inverse operator to Eq. (65),

hðrÞ ¼ hðr b Þ þ ðr  r b Þh0 ð1Þ  L1 NhðrÞ;

ð66Þ

which already involves both of the specified boundary conditions. Therefore we have

hðrÞ ¼ 1  L1 NhðrÞ:

ð67Þ

We decompose the solution and the nonlinearity

hðrÞ ¼

1 X

hn ðrÞ;

n¼0

NhðrÞ ¼ f ðr; hðrÞ; h0 ðrÞÞ ¼

1 X

An ðrÞ ¼

n¼0

1 X

An ðr; h0 ðrÞ; h1 ðrÞ; . . . ; hn ðrÞ; h00 ðrÞ; h01 ðrÞ; . . . ; h0n ðrÞÞ;

ð68Þ

n¼0

where An(r) are the multivariable Adomian polynomials. The definitional formula for the Adomian polynomials is

An ¼

  1 @n 0  ; n f ðr; hðr; kÞ; h ðr; kÞÞ n! @k k¼0

where the solution and its ordinary derivatives are parametrized by the grouping parameter k, i.e. hðr; kÞ ¼ P n Nhðr; kÞ ¼ 1 n¼0 k An ðrÞ. Upon substitution, we have 1 X

hn ðrÞ ¼ 1  L1

n¼0

1 X

An ðrÞ:

n¼0

Thus we have derived the recursion scheme for this nonlinear BVP

h0 ðrÞ ¼ 1; hnþ1 ðrÞ ¼ L1 An ðrÞ;

n P 0;

or

hnþ1 ðrÞ ¼ 

Z

r

rb

Z

r

An ðrÞ drdr;

1

where the first three Adomian polynomials are

ð69Þ P1

n n¼0 k hn ðrÞ,

and

4106

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118





e h00 ðrÞ 2 m2 rh0 ðrÞ A0 ¼  ; eh0 ðrÞ þ 1 rb ðeh0 ðrÞ þ 1Þ  2 e2 h1 ðrÞ h00 ðrÞ 2eh00 ðrÞh01 ðrÞ m2 r eh0 ðrÞh1 ðrÞ m2 rh1 ðrÞ þ  A1 ¼  þ ; 2 2 eh0 ðrÞ þ 1 rb ðeh0 ðrÞ þ 1Þ rb ðeh0 ðrÞ þ 1Þ ðeh0 ðrÞ þ 1Þ 

A2 ¼







e3 h1 ðrÞ2 h00 ðrÞ 2 e2 h2 ðrÞ h00 ðrÞ 2 2e2 h1 ðrÞh00 ðrÞh01 ðrÞ   ðeh0 ðrÞ þ 1Þ3 ðeh0 ðrÞ þ 1Þ2 ðeh0 ðrÞ þ 1Þ2  0 2 e h1 ðrÞ 2eh00 ðrÞh02 ðrÞ m2 r e2 h0 ðrÞh1 ðrÞ2 þ þ  eh0 ðrÞ þ 1 eh0 ðrÞ þ 1 r b ðeh0 ðrÞ þ 1Þ3 2 m2 r eh1 ðrÞ m2 reh0 ðrÞh2 ðrÞ m2 rh2 ðrÞ þ þ  : 2 2 r ð b eh0 ðrÞ þ 1Þ r b ðeh0 ðrÞ þ 1Þ rb ðeh0 ðrÞ þ 1Þ

By computation we obtain

 m2 3r þ r 3 þ 3rb  r3b ; 6ð1 þ eÞr b  1 h2 ðrÞ ¼ m4 15r 4 þ 2r 6  10r3 r b ð3 þ r2b Þ 360ð1 þ eÞ3 r 2b

h1 ðrÞ ¼

þ6rð8  15rb þ 5r 3b Þ þ r b ð48 þ 90r b  45r3b þ 8r 5b Þ ;

1 m4 ðr  r b Þ 6ð1 þ rÞ3 ð8 þ 9r þ 3r2 Þeð1 þ eÞ2 r b h3 ðrÞ ¼  5 3 360ð1 þ eÞ r b 1  þ m2 r6 ð10  40eÞ  55ð1 þ eÞ þ 90r4 e þ r8 ð1 þ 5eÞ 4   8r5 ð1 þ 2eÞrb ð3 þ r 2b Þ þ 8r3 8 þ 5ð1 þ 2eÞr b ð3 þ r2b Þ

 þ 2r2 r b 48  90rb þ 45r 3b  8r5b þ 10erb ð3 þ r 2b Þ2   þ2r b 96 þ 90rb þ 16r 2b  45r3b þ 8r 5b  2eð3 þ r 2b Þ 16 þ 5r b ð3 þ r2b Þ :

If e = 0, Eq. (64) will become a linear differential equation and its exact solution can be readily obtained by MATHEMATICA 8.0 in terms of the Airy functions and their derivatives,



h ðrÞje¼0

qffiffiffiffiffi qffiffiffiffiffi

qffiffiffiffiffi qffiffiffiffiffi 2 2 2 2 0 0 Ai 3 mr Bi r 3 mr  Ai r 3 mr Bi 3 mr b b b

qffiffiffiffiffib  qffiffiffiffi ffi : ¼ qffiffiffiffiffi qffiffiffiffiffi  2 2 2 2 0 0 Ai 3 mr Bi 3 mr r b  Ai 3 mr r b Bi 3 mr b

b

b

ð70Þ

b

We note that by the following relations [54]

pffiffiffi z ðI1=3 ðfÞ  I1=3 ðfÞÞ; 3 rffiffiffi z ðI1=3 ðfÞ þ I1=3 ðfÞÞ; BiðzÞ ¼ 3

AiðzÞ ¼

z ðI2=3 ðfÞ  I2=3 ðfÞÞ; 3 z 0 Bi ðzÞ ¼ pffiffiffi ðI2=3 ðfÞ þ I2=3 ðfÞÞ; 3 0

Ai ðzÞ ¼

ð71Þ ð72Þ

where f ¼ 23 z3=2 , we can then obtain the very same expression for the exact solution, for e = 0, as reported in [53] in terms of the modified Bessel functions. 1.00

0.95

0.90

0.85

0.80

0.2

0.4

0.6

0.8

1.0

r

Fig. 13. The exact solution (70) (solid line), and the approximations u2(r) (dot line), u3(r) (dashed line) and u4(r) (dot-side line) as e = 0.

4107

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

ERn r

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

r

0.05 0.10 0.15 0.20 0.25 0.30 Fig. 14. Curves of ER1(r) (solid line), ER2(r) (dot line) and ER3(r) (dashed line).

ERn r 0.0030 0.0025 0.0020 0.0015 0.0010 0.0005

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

r

0.0005 Fig. 15. Curves of ER3(r) (solid line), ER4(r) (dot line) and ER5(r) (dashed line).

ERn r

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

r

0.00001

0.00002

0.00003 Fig. 16. Curves of ER5(r) (solid line), ER6(r) (dot line) and ER7(r) (dashed line).

P For rb = 0.25, m = 0.5 and e = 0, we plot the exact solution (70), and the n-term approximations un ðrÞ ¼ n1 k¼0 hk ðrÞ for n = 2, 3, 4 in Fig. 13, where u4(r) and the exact solution almost overlap. Since the exact solution u⁄(x) of this problem as e – 0 is unknown, we instead investigate the error remainder function

ERn ðrÞ ¼ u00n ðrÞ þ

e 1 þ eun ðrÞ

ðu0n ðrÞÞ2 

m2 r un ðrÞ ; r b 1 þ eun ðrÞ

4108

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118 1.00

0.95

0.90

0.85

0.80

0.2

0.4

0.6

0.8

1.0

r

Fig. 17. The curves of the exact solution (70) (solid line), 7-term approximations u7(r) for e = 0.3 (dashed line), 0.1 (dot line), 0.1 (dot-side line) and  0.3 (dot-dot-side line).

for the n-term approximation un(r) in order to examine the convergence of the decomposition series solution since the limit of un(r) equals h⁄(r). For rb = 0.25, m = 0.3 and e = 0.2, we plot the functions ERn(r) for n = 1, 2, 3 in Fig. 14, the functions ERn(r) for n = 3, 4, 5 in Fig. 15 with different vertical scales, and the functions ERn(r) for n = 5, 6, 7 in Fig. 16 with different vertical scales, which demonstrates a remarkable accuracy for such low orders of approximation by the ADM as we expect the limit of the error function to approach zero. Next we examine the effect of e on the solution. For rb = 0.25 and m = 0.5, we plot the curves of the exact solution (70), 7-term approximations u7(r) for e = 0.3, 0.1, 0.1 and 0.3 in Fig. 17. The trend of the solution is consistent with that in [53] by the double decomposition method. Example 7. Consider the nonlinear second-order homogeneous partial differential equation [52]

e

hxx ðx; tÞ þ

1 þ ehðx; tÞ

ðhx ðx; tÞÞ2  K 2

hðx; tÞ ht ðx; tÞ  ¼ 0; 1 þ ehðx; tÞ 1 þ ehðx; tÞ

ð73Þ

where K depends on the physical properties and design parameters, and where h(x, t) has the domain of definition x 2 [0, 1], t 2 [0, 1) and subject to a mixed set of homogeneous Neumann and inhomogeneous Dirichlet boundary conditions, which includes a sinusoidally varying boundary value,

hx ð0; tÞ ¼ 0;

hð1; tÞ ¼ 1 þ S cosðBtÞ:

ð74Þ

The physical variable and parameters are h, x, t, e, K, S and B, which represent the dimensionless temperature, distance, time, thermal conductivity parameter, fin parameter, amplitude of oscillation and frequency of oscillation, respectively. The interested reader is referred to [52] for further details in regard to the derivation and design limitations of this engineering model. For convenience, we rewrite this nonlinear partial differential equation in Adomian’s operator-theoretic form

Lx hðx; tÞ þ Nhðx; tÞ ¼ 0; where Lx ðÞ ¼

@2 @x2

Nhðx; tÞ ¼

ð75Þ

ðÞ and

e 1 þ ehðx; tÞ

ðhx ðx; tÞÞ2  K 2

hðx; tÞ ht ðx; tÞ  : 1 þ ehðx; tÞ 1 þ ehðx; tÞ

We will use the partial solution concept of the ADM to solve the nonlinear partial differential equation in this example, which states that the temporal partial solution for its set of initial conditions, or t-partial solution, and the various spatial partial solutions for their respective sets of boundary conditions, such as the x-partial solution, are all identically equal, which was discovered by Adomian and Rach [55,56] and proven by Wazwaz [57]. Thus, for this particular BVP, we shall solve for approximations of the x-partial solution. Usually we solve for the partial solution in the coordinate, for which we most accurately know the values of its auxiliary conditions. But we also reserve the choice to average partial solutions for matching at the boundaries, when the values of more than one set of boundary conditions are known accurately, for sake of computational advantages. For sake of the form of boundary conditions, we define the inverse operator as

L1 x ðÞ ¼

Z 1

x

Z

x

ðÞ dxdx: 0

Applying L1 x to both sides of Eq. (75), we have

hðx; tÞ ¼ 1 þ S cosðBtÞ  L1 x Nhðx; tÞ;

4109

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

where the boundary conditions are used. Decomposing the solution and the nonlinearity

hðx; tÞ ¼

1 X

hn ðx; tÞ;

Nhðx; tÞ ¼

n¼0

1 X

An ðx; tÞ;

n¼0

where

An ðx; tÞ ¼

  1 @n  ; f ð x; t; hðx; t; kÞ; h ðx; t; kÞ; h ðx; t; kÞ Þ x t n  n! @k k¼0

for the nonlinearity

Nhðx; tÞ ¼ f ðx; t; hðx; tÞ; hx ðx; tÞ; ht ðx; tÞÞ: We list the first two Adomian polynomials as follows:

2 ð0;1Þ e hð1;0Þ ðx; tÞ 0 K 2 h0 ðx; tÞ h ðx; tÞ A0 ¼  þ  0 ; eh0 ðx; tÞ þ 1 eh0 ðx; tÞ þ 1 eh0 ðx; tÞ þ 1

2 ð1;0Þ 2 eK 2 h0 ðx; tÞh1 ðx; tÞ K 2 h1 ðx; tÞ e h1 ðx; tÞ h0 ðx; tÞ   A1 ¼ eh0 ðx; tÞ þ 1 ðeh0 ðx; tÞ þ 1Þ2 ðeh0 ðx; tÞ þ 1Þ2 þ

ð1;0Þ ð1;0Þ ð0;1Þ eh1 ðx; tÞhð0;1Þ ðx; tÞ 2eh0 ðx; tÞh1 ðx; tÞ h ðx; tÞ 0 þ  1 ; 2 eh0 ðx; tÞ þ 1 eh0 ðx; tÞ þ 1 ðeh0 ðx; tÞ þ 1Þ

where the superscript (1, 0) denotes the differentiation with respect to x, while the superscript (0, 1) denotes the differentiation with respect to t, etc. By the recursion scheme

h0 ðx; tÞ ¼ 1 þ S cosðBtÞ; hnþ1 ðx; tÞ ¼

L1 x An ;

ð76Þ

n P 0;

ð77Þ

we obtain

0.06

0.06

0.04

0.04

0.02 0.00 0.0

10 5

t

0.02 0.00 0.0

0.5 x

10 5

t

0.5 x 1.0

0

1.0

(a)

0

(b)

0.06

0.06

0.04

0.04

0.02 0.00 0.0

10 5 0.5 x

t

0.02 0.00 0.0

10 5

t

0.5 x 1.0

(c)

0

1.0

0

(d)

Fig. 18. Surfaces of the absolute error remainder functions (a) jER2(x, t)j, (b) jER3(x, t)j, (c) jER4(x, t)j and (d) jER5(x, t)j.

4110

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

h1 ðx; tÞ ¼

 ðx2  1Þ K 2 S cosðBtÞ  BS sinðBtÞ þ K 2

; 2ðeS cosðBtÞ þ e þ 1Þ



 ðx4  6x2 þ 5Þ 2S B2 ðe þ 1Þ  K 4 cosðBtÞ h2 ðx; tÞ ¼ 3 48ðeS cosðBtÞ þ e þ 1Þ  3B2 eS2 þ B2 eS2 cosð2BtÞ þ BeK 2 S2 sinð2BtÞ  þ2BeK 2 S sinðBtÞ  4BK 2 S sinðBtÞ þ 2K 4 ; ...:

Since the exact solution h⁄(x, t) of this problem is unknown, we instead investigate the error remainder function

ERn ðx; tÞ ¼ Lx un ðx; tÞ þ Nun ðx; tÞ; P for the n-term approximation un ðx; tÞ ¼ n1 k¼0 hk ðx; tÞ in order to examine the convergence of the decomposition series solu⁄ tion since the limit of un(x, t) equals h (x, t). We take S = 0.1, B = 1, K = 0.5, e = 0.2, and plot the error surfaces jERn(x, t)j for 0 6 x 6 1, 0 6 t 6 4p and n = 2, 3, 4 and 5, respectively, in Fig. 18, where the errors decrease significantly as n increases. In Fig. 19 we plot the surfaces of the approximations u5(x, t) in 0 6 x 6 1 and 0 6 t 6 4p for S = 0.1, B = 1, K = 0.5 and different e. We note that Yang et al. [52,53] solved the nonlinear BVPs in Examples 6 and 7 by the Adomian-Rach double decomposition method. Chiu and Chen [58] used the double decomposition method to design heat exchangers by solving nonlinear conduction–convection–radiation heat transfer equations subject to nonlinear boundary conditions. Serrano [59] investigated several nonlinear spatio-temporal hydrologic models by decomposition including moving boundary conditions, and Patel and Serrano [60] solved nonlinear boundary value problems of the multidimensional groundwater transport equations for irregularly shaped aquifer domains by the double decomposition method. Qin and Sun [61,62] solved nonlinear boundary value problems for nonlinear partial differential equations to design packed-bed catalytic reactors by decomposition. In [26,27] several different kinds of Duan’s parametrized recursion scheme were introduced, which leads to expansion of the domain of convergence for the decomposition solution series of initial value problems. We now demonstrate the effect of Duan’s convergence parameter for the following expository example of a nonlinear BVP, in which the right hand boundary

1.1

1.1

1.0

1.0

0.9

0.9

10

0.8 0.0

5

10

0.8 0.0

t

5

0.5 x

0.5 x 1.0

0

1.0

(a)

(b)

1.1 1.0 0.9

10

0.8 0.0

5

t

0.5 x 1.0

0

(c) Fig. 19. Surfaces of u5(x, t) (a) e = 0.3, (b) e = 0 and (c) e = 0.3.

0

t

4111

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

point lies just outside the interval of convergence of the unparametrized decomposition series. The intended difficulty is that the approximate solutions are similar to a truncated geometric series, so that a divergent sequence of approximations will result when the specified interval includes the point x = 1. Consider the Taylor expansion series of this particular exact solution, but where both boundary points lie within the interval  1 < x < 1. Example 8. Consider the nonlinear BVP with a quadratic nonlinearity

u00 ¼ 6u2 ;

0 6 x 6 1;

uð0Þ ¼ 1;

uð1Þ ¼ 1=4:

ð78Þ ð79Þ 2

⁄

The exact solution for this BVP is u (x) = (1 + x)

. According to Eq. (5), we have

3 u ¼ 1  x þ 6L1 u2  6x½L1 u2 x¼1 ; 4

ð80Þ

Rx Rx where L1 ðÞ ¼ 0 0 ðÞ dxdx. P The Adomian polynomials for the nonlinearity f(u) = u2 with the decomposition u ¼ 1 n¼0 un are

A0 ¼ u20 ;

A2 ¼ 2u0 u2 þ u21 ; . . . ; An ¼

A1 ¼ 2u0 u1 ;

n X

uk unk :

k¼0

If we use the following unparametrized recursion scheme

3 u0 ¼ 1  x; 4 unþ1 ¼ 6L1 An  6x½L1 An x¼1 ;

n P 0;

we obtain

57x 3x3 9x4 þ 3x2  þ ; 32 2 32 873x 57x3 555x4 9x5 9x6 27x7  u2 ¼ þ  þ  ; 896 16 128 4 16 448 ...:

u1 ¼ 

P ⁄ This particular sequence of partial sums un ðxÞ ¼ n1 k¼0 uk ðxÞ does not converge to the exact solution u (x) on the interval 0 6 x 6 1, but instead exhibits divergence. This can be displayed by the plots of the absolute errors jEn(x)j = jun(x)  u⁄(x)j. In Fig. 20 we show the jEn(x)j for n = 3, 4, 5 and 6. The fundamental difficulty in this expository example is that the right hand boundary point lies outside the interval of convergence of the unparametrized decomposition series, where the approximate solutions are similar to a truncated geometric series, which diverges at x = 1, e.g. consider the Taylor expansion series of this particular exact solution. We remark that these solution components cannot be improved, for example, by the Padé approximants technique, as applied to the sequence of n-term approximations, whereas Duan’s parametrized recursion scheme [26,27] can overcome the observed divergence, by pre-adjusting the unparametrized decomposition components as they are individually computed. Another way of looking at this difficulty is that the original constant term is in effect too large, which limits the interval of

En x

0.15

0.10

0.05

0.2

0.4

0.6

0.8

1.0

x

Fig. 20. Functions jE3(x)j (solid line), jE4(x)j (dot line), jE5(x)j (dashed line) and jE6(x)j (dot-side line).

4112

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

convergence, thus we diminish its contribution to the initial component and delay the remainder of its contribution to later solution components. The resulting parametrized decomposition series extends the interval of convergence to include the boundary point lying outside the previous interval of convergence for this BVP. Then we show how to easily solve such nonlinear BVPs by using the parametrized recursion scheme

3 u0 ¼ c þ 1  x; 4 u1 ¼ c þ 6L1 A0  6x½L1 A0 x¼1 ; unþ1 ¼ 6L1 An  6x½L1 An x¼1 ;

n P 1:

We have checked that, for c = 0.2,0.3 and  0.4, respectively, the maximal errors MEn = max06x61j un(x)  u⁄(x)j become gradually smaller and approach zero for n from 1 to 30. For example, as c = 0.2 we obtain

4 3x  ; 5 4 1 801x 48x2 6x3 9x4 u1 ¼  þ  þ ; 5 800 25 5 32 2 3 33003x 24x 951x 36591x4 36x5 9x6 27x7 þ u2 ¼   þ  þ  ; 112000 25 500 16000 25 20 448 ...:

u0 ¼

The absolute values of the error functions jEn(x)j, for n = 2, 3, 4 and 5, are plotted in Fig. 21. The maximal errors MEn, for n = 1 through 16, are listed in Table 3. Alternatively, we can also use the parametrized recursion scheme

u0 ¼ c; 3 u1 ¼ c þ 1  x þ 6L1 A0  6x½L1 A0 x¼1 ; 4 unþ1 ¼ 6L1 An  6x½L1 An x¼1 ; n P 1: We have checked that for c = 0.1, 0.2, 0.3, 0.4 and 0.5, respectively, the maximal errors MEn become gradually smaller and approach zero. For example, for c = 0.4, we have

En x

0.025

0.020

0.015

0.010

0.005

0.2

0.4

0.6

0.8

1.0

x

Fig. 21. Functions jE2(x)j (solid line), jE3(x)j (dot line), jE4(x)j (dashed line) and jE5(x)j (dot-side line).

Table 3 For c = 0.2, maximal errors MEn for n = 1,2, . . . ,16. n MEn

1 0.2

n MEn

9 0.000199

2 0.027527

3 0.013583

4 0.005345

5 0.001447

6 0.000180

7 0.000464

10 0.000065

11 0.000009

12 0.000026

13 0.000022

14 0.000013

15 0.000004

8 0.000375 16 0.0000007

4113

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118 Table 4 For c = 0.4, maximal errors MEn for n = 1,2, . . . ,16. n MEn

1 0.6

n MEn

9 0.000660

2 0.082902

3 0.014544

4 0.010142

5 0.007595

6 0.002760

7 0.000979

8 0.001312

10 0.000115

11 0.000285

12 0.000185

13 0.000033

14 0.000067

15 0.000055

16 0.000015

6 0.0125

7 0.00625

14 0.000049

15 0.000024

Table 5 For c = 0.4, maximal errors MEn for n = 1,2, . . . ,16. n MEn

1 0.4

n MEn

9 0.001563

2 0.2 10 0.000781

3 0.1

4 0.05

11 0.000391

12 0.000195

5 0.025 13 0.000098

8 0.003125 16 0.000012

2 ; 5 3 123x 12x2 u1 ¼  þ ; 5 100 25 81x 36x2 123x3 24x4 þ u2 ¼   þ ; 125 25 125 125 ...; u0 ¼

and the maximal errors MEn = max06x61 jun(x)  u⁄(x)j, for n = 1,2, . . . ,16, are listed in Table 4. Next, we consider the parametrized recursion scheme with the infinite decomposition of the convergence parameter c [26,27], where we have used a specified rule to decompose c,

3 u0 ¼ c þ 1  x; 4 c unþ1 ¼  nþ1 þ 6L1 An  6x½L1 An x¼1 ; 2

n P 0:

For c = 0.4 we have

3 3x  ; 5 4 1 369x 27x2 9x3 9x4 þ u1 ¼   þ ; 5 800 25 10 32 2 1 36699x 18x 1707x3 15903x4 81x5 27x6 27x7  þ u2 ¼  þ  þ  : 10 112000 25 2000 16000 100 80 448

u0 ¼

The maximal errors MEn = max06x61 jun(x)  u⁄(x)j, for n = 1 through 16, are listed in Table 5, where this sequence {MEn} decreases strictly monotonically as compared with those in Tables 3 and 4. We note that the class of second-order differential equations subject to a set of only Neumann boundary conditions currently are not amenable by our new approach, including such problems as higher order differential equations for the case of sets of (k  1)-order generalized Neumann boundary conditions occurring at all boundary points, where all boundary points are distinct. We thus plan further research for extending our new approach to include more sophisticated applications. 4. Conclusion We have presented a new modification of the ADM to conveniently solve a wide class of multi-order and multi-point nonlinear BVPs. Our new approach yields a rapidly convergent sequence of analytic functions for the solution without resort to incorporating undetermined coefficients within the recursion scheme for computing successive solution components as practiced in prior art. These undetermined coefficients are inserted due to the nature of the boundary conditions and the choice of the inverse linear operator. Thus we avoid the complications resulting from the necessity of evaluating such undetermined coefficients at each stage of approximation. As this resulting sequence of nonlinear algebraic equations produces more than one root for the possible value of each undetermined coefficient, we are then confronted with the ambiguity of multiple roots. Usually this difficulty is handled by discarding unphysical roots, which does not lend itself to automation. Furthermore we can also parametrize the new modified recursion scheme [26,27] by inserting a predetermined parameter in order to achieve simple-to-integrate series without regard to the possible complexity of the original initial solution component, which is comprised of the boundary values and the input function. We have verified that the value of the

4114

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

predetermined parameter c can favorably affect the rate of convergence and the domain of convergence of the resulting parametrized decomposition series solution, thus permitting another broad avenue to adjust the convergence properties for an even more rapid rate of convergence. We have demonstrated the practicality and efficiency of our new modification of the ADM by eight numerical examples for nonlinear BVPs. In point of fact, the convergence was already sufficiently rapid that we plotted the error of the low-order approximations and also presented the logarithmic plots of the maximal errors instead of the approximations to demonstrate the rapid convergence of our new modified decomposition series. The ADM permits the convenient solution of nonlinear ordinary and partial differential equations for either initial value problems or boundary value problems for any analytic nonlinearity f(u), f(t, u(t),u0 (t), . . . , u(k1)(t)), . . . , and f ðt; ~ x; uðt; ~ xÞ; ut ðt; ~ xÞ; . . . ; uxy ðt; ~ xÞ; . . . ; uzkz 1 ðt; ~ xÞÞ, including mixed partial derivatives, where ~ x ¼ ðx; y; zÞ, and so on. We have solved these eight nonlinear multi-order and multi-point BVPs for an exponential, a radical, a hyperbolic sine, a quadratic and a product nonlinearity in a straightforward procedure using the Adomian polynomials, including cases of two-, three- and four-point as well as second-, third- and fourth-order BVPs. The overall efficiency of our new modification of the ADM is further enhanced by the new algorithms and subroutines crafted by Duan [26,46,47,37] for generating the Adomian polynomials quickly and to high orders. We remark that, because of the usual rapid convergence of the Adomian decomposition series, we often only require a low-stage approximation of the solution for high accuracy. Furthermore the rigorous mathematical convergence of the Adomian decomposition series to the exact solution has already been established by many researchers. In many physical models, including several examples within this paper, a closed-form solution is not always possible. However our expository examples have once again demonstrated that only a few terms of the truncated decomposition series does provide an excellent approximation even in the case of nonlinear BVPs. Acknowledgement This work was supported in part by the Innovation Program of the Shanghai Municipal Education Commission (No. 11YZ225). Appendix A. Undetermined coefficient method in the ADM We demonstrate this method by the nonlinear BVP in Example 1,

u00 ðxÞ ¼ eu ; uð0Þ ¼ 0;

0 6 x 6 1; uð1Þ ¼ 0:

Applying the inverse linear operator L1 ðÞ ¼ ary value at x = 0 yields

Rx Rx 0

0

ðÞ dxdx to both sides of the differential equation and using the bound-

uðxÞ ¼ Bx þ L1 eu ;

ðA:1Þ

where B = u0 (0) is an undetermined coefficient. Next we decompose the solution u(x) and the nonlinearity eu,

uðxÞ ¼

1 X

un ðxÞ;

eu ¼

n¼0

1 X

An ;

ðA:2Þ

n¼0

where An are the Adomian polynomials for the nonlinearity Nu = eu. If we use the recursion scheme

u0 ¼ Bx; un ¼ L1 An1 ;

n ¼ 1; 2; . . . ;

we have

u1 ¼ u2 ¼

1 þ eBx  Bx

; B2 5 þ e2Bx  2Bx þ eBx ð4  4BxÞ 4B4

;

... and

u1 ðxÞ ¼ Bx; u2 ðxÞ ¼ Bx þ u3 ðxÞ ¼ Bx þ ...:

1 þ eBx  Bx B2 1 þ eBx  Bx 2

B

; þ

5 þ e2Bx  2Bx þ eBx ð4  4BxÞ 4B4

;

4115

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

For the boundary value at x = 1, if we match u1(x), then B = 0, and u1(x) = 0. B If we match u2(x), then we need to solve the equation B þ 1Bþe ¼ 0. By the native command ‘‘FindRoot’’ in MATHEMB2 ATICA, we obtain B = 0.4347754841, and then

u2 ðxÞ ¼ 5:29017 þ 5:29017e0:434775x þ 1:86526x: If we match u3(x), then we need to solve the equation

Bþ

1  B þ eB 2

B

þ

5  2B þ ð4  4BÞeB þ e2B 4B4

¼ 0:

Solving the transcendental equation yields B = 0.4604486353 and then

u3 ðxÞ ¼ 32:5257 þ 5:5618e0:920897x þ 26:9639e0:460449x þ 6:83319x þ 10:2437e0:460449x x: If we match u4(x), we obtain by a similar manner, B = 0.4632279437 and

u4 ðxÞ ¼ 217:363 þ 8:43432e1:38968x þ 56:0355e0:926456x þ 152:893e0:463228x þ 30:1678x þ 23:4421e0:926456x x þ 80:3867e0:463228x x þ 10:859e0:463228x x2 : By the exact solution (34) we calculate that u0 (0) = 0.4636325917 . . ., the true value of B. In Fig. 22 we plot the curves of the exact solution u⁄(x), and the approximations u2(x), u3(x) and u4(x), where u4(x) and u⁄(x) almost overlap one another. If we use the alternate recursion scheme

u0 ¼ 0; u1 ¼ Bx þ L1 A0 ; un ¼ L1 An1 ;

n ¼ 2; 3; . . . ;

we have instead

u1 ¼ Bx þ

x2 ; 2

u2 ¼

Bx3 x4 þ ; 6 24

u3 ¼

B2 x4 Bx5 x6 þ þ ; 24 30 180

u4 ¼

B3 x5 11B2 x6 17Bx7 17x8 þ þ þ ; 120 720 2520 20160 ...:

and

0.2

0.4

0.6

0.8

1.0

x

0.02

0.04

0.06

0.08

0.10

Fig. 22. Exact solution u⁄(x) (solid line), approximations u2(x) (dot line), u3(x) (dashed line) and u4(x) (dot-side line).

4116

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118 Table 6 The roots of matching equations in B. n

Real roots for B

2 3 4 5 6 7 8 9

0.5 0.4642857143 0.463477239; 28.33652276 0.4636228337 0.4636337428; 11.90426815 0.4636327238 0.4636325836; 14.71306454 0.4636325900; 9.315847233; 23.85884627

u1 ðxÞ ¼ 0; u2 ðxÞ ¼ Bx þ

x2 ; 2

u3 ðxÞ ¼ Bx þ

x2 Bx3 x4 þ þ ; 2 6 24

u4 ðxÞ ¼ Bx þ

x2 Bx3 x4 B2 x4 Bx5 x6 þ þ þ þ þ ; 2 6 24 24 30 180

u5 ðxÞ ¼ Bx þ

x2 Bx3 x4 B2 x4 Bx5 B3 x5 x6 11B2 x6 17Bx7 17x8 þ þ þ þ þ þ þ þ þ ; 2 6 24 24 30 120 180 720 2520 20160

...: Matching the boundary value at x = 1, by u2(x) we have B ¼  12 and

x 2

u2 ðxÞ ¼  þ by u3(x) we have

u3 ðxÞ ¼ 

13 24

x2 ; 2 þ 7B ¼ 0, hence B ¼  13 and 6 28

13x x2 13x3 x4 þ  þ ; 28 2 168 24 2

B by u4(x) we solve 197 þ 6B þ 24 ¼ 0, which yields two roots 0.463477239 and 28.33652276 for B. Matching u5(x) we solve 360 5 3683 6720

2

3

B þ 3041B þ 41B þ 120 ¼ 0, which yields one real root 0.4636228337 for B. 2520 720 The real roots of the matching nonlinear algebraic equations in B for un(x), n = 2, 3, . . . , 9, are listed in Table 6. As n increases, we observe that one root is approaching 0.46363259. . .. If such a recursion scheme were to contain more than one undetermined parameter, then the matching equations would become a system of multivariable nonlinear equations. We note that the Padé approximants technique [6,63,64] is sometimes used in this method of undetermined coefficients [8,15].

Appendix B. Adomian-Rach double decomposition method Consider the nonlinear BVP

Lu ¼ Nu þ gðxÞ; uðaÞ ¼ a;

a 6 x 6 b;

uðbÞ ¼ b:

In the double decomposition method, the inverse linear operator L1 is taken as a twofold indefinite integration for second-order differential equations [10,11], i.e.

L1 ðÞ ¼ C 0 þ C 1 x þ

ZZ

ðÞ dxdx;

ðB:1Þ

RR where C0 and C1 are constants of integration, which are called the matching coefficients, and where I2x ðÞ ¼ ðÞ dxdx denotes 1 pure integrations. Applying the operator L to both sides of the differential equation yields the nonlinear integral equation

uðxÞ ¼ C 0 þ C 1 x þ I2x Nu þ I2x gðxÞ;

ðB:2Þ

where C0 and C1 are arbitrary constants of integration to be determined by decomposition for each stage of approximation. The double decomposition method decomposes the solution u(x), the nonlinearity Nu, and the constants C0 and C1:

4117

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

uðxÞ ¼

1 X

un ðxÞ;

Nu ¼

n¼0

1 X

An ;

C0 ¼

n¼0

1 X

C 0;n ;

n¼0

C1 ¼

1 X

C 1;n :

n¼0

Upon substitution of these series into Eq. (B.2), we design the recursion scheme

u0 ¼ C 0;0 þ C 1;0 x þ I2x gðxÞ; un ¼ C 0;n þ C 1;n x þ I2x An1 ;

n ¼ 1; 2; . . . :

The constants C0,k and C1,k are determined by matching each partial sum un(x) to the boundary values for the case of Dirichlet boundary conditions. That is, matching u1(x) = u0 to the boundary values determines the values of C0,0 and C1,0, matching P u2(x) = u0 + u1 to the boundary values determines the values of C0,1 and C1,1, . . ., and matching unþ1 ðxÞ ¼ nk¼0 uk to the boundary values determines the values of C0,n and C1,n. For example, for the nonlinear BVP in Example 1,

u00 ðxÞ ¼ eu ; uð0Þ ¼ 0;

0 6 x 6 1; uð1Þ ¼ 0;

the Adomian-Rach modified recursion scheme is

u0 ¼ C 0;0 þ C 1;0 x; un ¼ C 0;n þ C 1;n x þ I2x An1 ;

n P 1;

where the An are the Adomian polynomials for the nonlinearity Nu = eu. Matching u1(x) = u0 to the boundary values determines that C0,0 = C1,0 = 0. Thus u0 = 0, and

u1 ¼ C 0;1 þ C 1;1 x þ I2x eu0 ¼ C 0;1 þ C 1;1 x þ

x2 : 2

Matching u2(x) = u0 + u1 to the boundary values determines C0,1 = 0 and C1,1 = 1/2. Thus u1 = (x2  x)/2 and

u2 ¼ C 0;2 þ C 1;2 x þ I2x u1 eu0 ¼ C 0;2 þ C 1;2 x 

x3 x4 þ : 12 24 3

4

x x x Matching u3(x) = u0 + u1 + u2 to the boundary values determines that C0,2 = 0 and C1,2 = 1/24. Thus u2 ¼ 24  12 þ 24 , and so on. We observe that these results are consistent with the solution components as computed in Example 1.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

G. Adomian, Stochastic Systems, Academic, New York, 1983. G. Adomian, R. Rach, Inversion of nonlinear stochastic operators, J. Math. Anal. Appl. 91 (1983) 39–46. G. Adomian, Nonlinear Stochastic Operator Equations, Academic, Orlando, 1986. G. Adomian, Nonlinear Stochastic Systems Theory and Applications to Physics, Kluwer Academic, Dordrecht, 1989. G. Adomian, R. Rach, R. Meyers, An efficient methodology for the physical sciences, Kybernetes 20 (1991) 24–34. G. Adomian, Solving Frontier Problems of Physics: The Decomposition Method, Kluwer Academic, Dordrecht, 1994. R. Rach, A new definition of the Adomian polynomials, Kybernetes 37 (2008) 910–955. A.M. Wazwaz, Partial Differential Equations and Solitary Waves Theory, Higher Education, Beijing, Springer, Berlin, 2009. A.M. Wazwaz, Linear and Nonlinear Integral Equations: Methods and Applications, Higher Education, Beijing, Springer, Berlin, 2011. G. Adomian, R. Rach, Analytic solution of nonlinear boundary-value problems in several dimensions by decomposition, J. Math. Anal. Appl. 174 (1993) 118–137. G. Adomian, R. Rach, A new algorithm for matching boundary conditions in decomposition solutions, Appl. Math. Comput. 58 (1993) 61–68. G. Adomian, R. Rach, Modified decomposition solution of linear and nonlinear boundary-value problems, Nonlinear Anal. 23 (1994) 615–619. A.M. Wazwaz, Approximate solutions to boundary value problems of higher order by the modified decomposition method, Comput. Math. Appl. 40 (2000) 679–691. A.M. Wazwaz, The modified Adomian decomposition method for solving linear and nonlinear boundary value problems of 10th-order and 12th-order, Int. J. Nonlinear Sci. Numer. Simul. 1 (2000) 17–24. A.M. Wazwaz, A reliable algorithm for obtaining positive solutions for nonlinear boundary value problems, Comput. Math. Appl. 41 (2001) 1237–1244. A.M. Wazwaz, The numerical solution of fifth-order boundary value problems by the decomposition method, J. Comput. Appl. Math. 136 (2001) 259– 270. A.M. Wazwaz, The numerical solution of sixth-order boundary value problems by the modified decomposition method, Appl. Math. Comput. 118 (2001) 311–325. A.M. Wazwaz, A reliable algorithm for solving boundary value problems for higher-order integro-differential equations, Appl. Math. Comput. 118 (2001) 327–342. A.M. Wazwaz, The numerical solution of special fourth-order boundary value problems by the modified decomposition method, Int. J. Comput. Math. 79 (2002) 345–356. M. Tatari, M. Dehghan, The use of the Adomian decomposition method for solving multipoint boundary value problems, Phys. Scripta 73 (2006) 672– 676. M. Dehghan, M. Tatari, Finding approximate solutions for a class of third-order non-linear boundary value problems via the decomposition method of Adomian, Int. J. Comput. Math. 87 (2010) 1256–1263. B. Jang, Two-point boundary value problems by the extended Adomian decomposition method, J. Comput. Appl. Math. 219 (2008) 253–262.

4118

J.-S. Duan, R. Rach / Applied Mathematics and Computation 218 (2011) 4090–4118

[23] A.E. Ebaid, Exact solutions for a class of nonlinear singular two-point boundary value problems: the decomposition method, Z. Naturforsch. 65a (2010) 1–6. [24] A. Ebaid, A new analytical and numerical treatment for singular two-point boundary value problems via the Adomian decomposition method, J. Comput. Appl. Math. 235 (2011) 1914–1924. [25] W. Al-Hayani, Adomian decomposition method with Green’s function for sixth-order boundary value problems, Comp. Math. Appl. 61 (2011) 1567– 1575. [26] J.S. Duan, Recurrence triangle for Adomian polynomials, Appl. Math. Comput. 216 (2010) 1235–1241. [27] J.S. Duan, R. Rach, On the domain of convergence of the Adomian series solution, Appl. Math. Comput., in press. [28] Y. Cherruault, Convergence of Adomian’s method, Kybernetes 18 (1989) 31–38. [29] A. Rèpaci, Nonlinear dynamical systems: on the accuracy of Adomian’s decomposition method, Appl. Math. Lett. 3 (1990) 35–39. [30] K. Abbaoui, Y. Cherruault, Convergence of Adomian’s method applied to differential equations, Comput. Math. Appl. 28 (1994) 103–109. [31] L. Gabet, The theoretical foundation of the Adomian method, Comput. Math. Appl. 27 (1994) 41–52. [32] A. Abdelrazec, D. Pelinovsky, Convergence of the Adomian decomposition method for initial-value problems, Numer. Methods Partial Differ. Eqs. 27 (2011) 749–766. [33] R.P. Agarwal, Boundary Value Problems for Higher Order Differential Equations, World Scientific, Singapore, 1986. [34] D.R.K.S. Rao, K.N. Murthy, A.S. Rao, On three-point boundary value problems associated with third order differential equations, Nonlinear Anal. 5 (1981) 669–673. [35] K.N. Murty, B.D.C.N. Prasad, Three-point boundary value problems – existence and uniqueness, Yokohama Math. J. 29 (1981) 101–105. [36] A. Boucherif, N. Al-Malki, Nonlinear three-point third-order boundary value problems, Appl. Math. Comput. 190 (2007) 1168–1177. [37] J.S. Duan, New recurrence algorithms for the nonclassic Adomian polynomials, Comput. Math. Appl. 62 (2011) 2961–2977. [38] J.S. Duan, New ideas for decomposing nonlinearities in differential equations, Appl. Math. Comput. 218 (2011) 1774–1784. [39] R. Rach, A convenient computational form for the Adomian polynomials, J. Math. Anal. Appl. 102 (1984) 415–419. [40] G. Adomian, R. Rach, On composite nonlinearities and the decomposition method, J. Math. Anal. Appl. 113 (1986) 504–509. [41] A.M. Wazwaz, A new algorithm for calculating Adomian polynomials for nonlinear operators, Appl. Math. Comput. 111 (2000) 53–69. [42] F. Abdelwahid, A mathematical model of Adomian polynomials, Appl. Math. Comput. 141 (2003) 447–453. [43] J. Biazar, M. Ilie, A. Khoshkenar, An improvement to an alternate algorithm for computing Adomian polynomials in special cases, Appl. Math. Comput. 173 (2006) 582–592. [44] Y. Zhu, Q. Chang, S. Wu, A new algorithm for calculating Adomian polynomials, Appl. Math. Comput. 169 (2005) 402–416. [45] M. Azreg-Aı¨nou, A developed new algorithm for evaluating Adomian polynomials, CMES-Comput. Model. Eng. Sci. 42 (2009) 1–18. [46] J.S. Duan, An efficient algorithm for the multivariable Adomian polynomials, Appl. Math. Comput. 217 (2010) 2456–2467. [47] J.S. Duan, Convenient analytic recurrence algorithms for the Adomian polynomials, Appl. Math. Comput. 217 (2011) 6337–6348. [48] J.S. Duan, A.P. Guo, Reduced polynomials and their generation in Adomian decomposition methods, CMES-Comput. Model. Eng. Sci. 60 (2010) 139– 150. [49] M.R. Scott, W.H. Vandevender, A comparison of several invariant imbedding algorithms for the solution of two-point boundary-value problems, Appl. Math. Comput. 1 (1975) 187–218. [50] G. Adomian, R. Rach, D. Sarafyan, On the solution of equations containing radicals by the decomposition method, J. Math. Anal. Appl. 111 (1985) 423– 426. [51] G. Adomian, R. Rach, Light scattering in crystals, J. Appl. Phys. 56 (1984) 2592–2594. [52] Y.T. Yang, S.K. Chien, C.K. Chen, A double decomposition method for solving the periodic base temperature in convective longitudinal fins, Energy Convers. Manage. 49 (2008) 2910–2916. [53] Y.T. Yang, C.C. Chang, C.K. Chen, A double decomposition method for solving the annular hyperbolic profile fins with variable thermal conductivity, Heat Transfer Eng. 31 (2010) 1165–1172. [54] M. Abramowitz, I.A. Stegun (Eds.), Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (10th printing), Dover, New York, 1972. [55] G. Adomian, R. Rach, Equality of partial solutions in the decomposition method for linear or nonlinear partial differential equations, Comput. Math. Appl. 19 (1990) 9–12. [56] G. Adomian, R. Rach, A further consideration of partial solutions in the decomposition method, Comput. Math. Appl. 23 (1992) 51–64. [57] A.M. Wazwaz, Equality of partial solutions in the decomposition mehtod for partial differential equations, Int. J. Comput. Math. 65 (1997) 293–308. [58] C.H. Chiu, C.K. Chen, Application of the decomposition method to thermal stresses in isotropic circular fins with temperature-dependent thermal conductivity, Acta Mech. 157 (2002) 147–158. [59] S.E. Serrano, Hydrology for Engineers, Geologists, and Environmental Professionals: An Integrated Treatment of Surface, Subsurface, and Contaminant Hydrology, second revised ed., HydroScience, Ambler, PA, 2010. [60] A. Patel, S.E. Serrano, Decomposition solution of multidimensional groundwater equations, J. Hydrol. 397 (2011) 202–209. [61] X.Y. Qin, Y.P. Sun, Approximate analytic solutions for a two-dimensional mathematical model of a packed-bed electrode using the Adomian decomposition method, Appl. Math. Comput. 215 (2009) 270–275. [62] X.Y. Qin, Y.P. Sun, Approximate analytical solutions for a mathematical model of a tubular packed-bed catalytic reactor using an Adomian decomposition method, Appl. Math. Comput. 218 (2011) 1990–1996. [63] R. Rach, J.S. Duan, Near-field and far-field approximations by the Adomian and asymptotic decomposition methods, Appl. Math. Comput. 217 (2011) 5910–5922. [64] H. Chu, Y. Zhao, Y. Liu, A MAPLE package of new ADM-Padé approximate solution for nonlinear problems, Appl. Math. Comput. 217 (2011) 7074–7091.