An efficient method for solving Bratu equations

Applied Mathematics and Computation 176 (2006) 704–713 www.elsevier.com/locate/amc

An efficient method for solving Bratu equations Muhammed I. Syam *, Abdelrahem Hamdan Department of Mathematical Sciences, College of Science, UAE University, P.O. Box 17551, Al-Ain, United Arab Emirates

Abstract In this paper, we present a numerical technique for solving Bratu equation. It is based on the Laplace Adomain decomposition method which produces an implicit equation in two variables. We used the predictor corrector technique to trace the solution curve generated from this equation. Numerical results and conclusions will be presented. 2005 Elsevier Inc. All rights reserved. Keywords: Laplace Adomain decomposition method; Predictor corrector; Turning and bifurcation points

1. Introduction Many problems in science and engineering require the computation of a family of solutions of a nonlinear system of the form Gðu; kÞ ¼ 0; nþ1

u ¼ uðkÞ;

ð1Þ

n

where G : R ! R is continuously differentiable function, u represents the solution and k is a real parameter, i.e., ReynoldÕs number, load, . . ., etc. To solve these problems, we need to find the solution for some k-intervals, i.e., a path solutions, (u(k), k). Equations of the form (1) are called nonlinear elliptic eigenvalue problems if the operator G with k fixed is an elliptic differential operator. As a typical example of a nonlinear elliptic eigenvalue problems, we consider the following problem: Gðu; kÞ ¼ Du þ kf ðuÞ ¼ 0; u ¼ 0;

in X;

on oX.

ð2Þ

Eq. (2) arises in many physical problems. For example, in chemical reactor theory, radiative heat transfer, combustion theory, and in modelling the expansion of the universe. In this paper, we will concentrate on one of the forms of Eq. (2), which is the Bratu equation Du þ keu ¼ 0; u ¼ 0; *

in X;

on oX.

Corresponding author. E-mail addresses: [email protected] (M.I. Syam), [email protected] (A. Hamdan).

0096-3003/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2005.10.021

ð3Þ

M.I. Syam, A. Hamdan / Applied Mathematics and Computation 176 (2006) 704–713

705

Problem (3) has some singular points which had been studied numerically and theoretically. The behavior of the solution near the singular points has been studied numerically [4,5,12] and theoretically [6,11,14,15]. In this paper, we will use the Adomain decomposition method (ADM) and the Laplace transforms to solve Eq. (3), see [1,2,18,19]. We will combine between these two techniques. Also, we will approximate the solution of this problem by a polynomial of degree n in terms of x. The boundary conditions will produce an implicit equation in two parameters. One of these parameters is k. We will use the predictor corrector technique to follow the implicit curve. Many authors used the path following method for solving problem (3) such as [4,8,9]. However, they discretize the problem first using the finite difference method or the finite element method, then they apply the path following technique. They face serious problems at the turning and bifurcation points since the Jacobian is singular at these points. For this reason, special treatments were done at these points. These treatments affect the accuracy of the solution and the complexity of the problem. However, in our approach we will produce an implicit equation so that the Jacobian of the problem is nonsingular at these points. This will give us freedom to use the predictor corrector technique. We will use the Euler method as a predictor and NewtonÕs method as a corrector. We will organize our paper in the following way. We will review the path following method in Section 2. Also, some related facts and theorems will be given there. In Section 3, we will present the Laplace Adomain decomposition method (LADM) for solving problem (3). Numerical results and conclusions will be given in Section 4. 2. Review of path following method The problem under consideration here is an implicit equation in two variables of the form H ðx; kÞ ¼ 0;

ð4Þ

where H : R · R ! R is a smooth map; that is, it has as many continuous derivatives as the discussion requires. By letting X = (x, k), Eq. (4) can be written in the form H ðX Þ ¼ 0.

ð5Þ

If the Jacobian H 0 (X) of H 0 (X) has full rank at a point X, the point X is called a regular point and Y is a regular value if H(X) = Y. While if H 0 (X) has a drop in rank at a point X, such a point will be called a singular point. If X0 is a regular point, then the Implicit Function Theorem concludes that the solution H1(0) can be locally parameterized with respect to one component of x or some parameter, say arclength s. Lemma 2.1. If H(X) is as given in (5) and H 0 (X) has full rank. Then there exist a smooth solution curve c(s) : I ! R2, where s 2 I with I some open interval such that cð0Þ ¼ X 0 ; H 0 ðcðsÞÞ_cðsÞ ¼ 0; k_cðsÞk ¼ 1; 0 H ðcðsÞÞ det > 0. c_ ðsÞ

ð6Þ

This lemma uniquely determines the tangent vector c_ ðsÞ with specific orientation induced by the Jacobian matrix H 0 (c(s)). This solution curve c(s) is characterized as the solution of the initial value problem X_ ¼ tðH 0 ðX ÞÞ;

X ð0Þ ¼ X 0 .

ð7Þ

This fact suggests using numerical methods for solving initial value problems to numerically trace c. A typical path following method consists of a succession of two steps:

706


Predictor step: An approximation step along the curve, usually in the direction of the tangent to the curve. It is obtained as a simple numerical integration step for (7). A common method that is often used is Euler-predictor given by Z ¼ X þ htðH 0 ðX ÞÞ; where X is a point lying along the solution curve c and h > 0 is the stepsize. Corrector step: This step is used to bring the predicted point back to the curve. This requires few iterative steps for solving H(X) = 0. Gauss–Newton corrector is used since we have an overdetermined system. It is given by W ¼ Z H 0 ðZÞy H ðZÞ; where H is an 1 · 2 matrix with maximal rank and (H 0 ) = H*(HH*)1 is called the Moore– Penrose inverse and H* denotes the Hermition transpose of H. More details can be found in Allgower and Georg [3] and Syam and Siyyam [16,17]. We will classify the points into regular and singular points. At the beginning of this section, we define the meaning of regular points. In the next definition, we will retain back to Eq. (1) and we will define the singular points. ~ ðGu ðu0 ; k0 ÞÞ is Definition 2.2. A point (u0, k0) is called a singular point if the null space of Gu(u0, k0) denoted by N of dimension at least one. It is called a simple singular point, if the null space of Gu(u0, k0) is of dimension one. This means, it is spanned by a nonzero vector u0. The range in the simple singular case is then given by RangeðGu ðu0 ; k0 ÞÞ ¼ fy : wt0 y ¼ 0g for some nonzero vector w0. Now, singular points on the solution path are either turning points or bifurcation points. The determination of the solution path in a neighborhood of a turning point or bifurcation point requires special care. It is therefore important to detect singular points on the solution path. One can define the fold and the bifurcation points in the following way. Definition 2.3. A singular point (u0, k0) is said to be a fold (turning or limit) point if wt0 Gk ðu0 ; k0 Þ 6¼ 0. If in addition wt0 Guuðu0 ; k0 Þ½u0 ; u0 6¼ 0 then it is called a quadratic fold point. Definition 2.4. A point ðu0 ; k0 Þ 2 Rnþ1 is called a simple bifurcation point of the equation G(u, k) = 0 if the following conditions hold: Gðu0 ; kÞ ¼ 0;

dimðKerðGu ðu0 ; k0 ÞÞÞ ¼ 2

and e Guu ðu0 ; k0 ÞjðKerðGu ðu0 ;k0 ÞÞ2

has one positive and one negative eigenvalue, where e spans Ker Gu(u0, k0)*. For multiple bifurcation points, the dimension of the kernel of Gu(u0, k0) is greater than two. Such points frequently are of intrinsic interest, and there are special algorithms for detecting and calculating them. We refer the reader to [7,9,10]. 3. Laplace Adomain decomposition method In this section we want to describe how to use LADM for Problem (3). First, we will study the one-dimensional case when X = [0, 1]. Apply the Laplace transform to both sides of the Bratu equation to get Lðu00 Þ þ kLðeu Þ ¼ 0.


707

Using the properties of the Laplace transforms, we get s2 LðuÞ uð0Þs u0 ð0Þ þ kLðeu Þ ¼ 0.

ð8Þ

Let u 0 (0) = A. Then, using the condition u(0) = 0, Eq. (8) becomes LðuÞ ¼

1 ðA kLðeu ÞÞ. s2

ð9Þ

The Adomain decomposition method assumes that the solution u(x) of the Bratu equation and the nonlinear function eu can be written in the series form as follows: 1 1 X X un ðxÞ and eu ¼ An . ð10Þ uðxÞ ¼ n¼0

n¼0

We can derive the values of An by expanding the function eu about u0 and using Eq. (10) to get P1 2 P1 2 n¼1 un u u0 u0 u u0 u0 ðu u0 Þ u0 u0 u0 n¼1 un þe þ ¼ e þ e þe þ . e ¼ A0 þ A1 þ A2 þ ¼ e þ e 1! 2! 1! 2! Thus, the first few Adomain polynomials are given by A0 ¼ eu0 ; A1 ¼ u1 eu0 ; u21 u0 e þ u2 eu0 ; 2! u3 2u1 u2 u0 e þ u3 eu0 ; A3 ¼ 1 eu0 þ 3! 2! u4 3u2 u2 2u1 u3 þ u22 u0 e þ u4 eu0 . A4 ¼ 1 eu0 þ 1 eu0 þ 4! 3! 2!

A2 ¼

From the previous discussion, we see that Ak can be given by the following formula: " !# 1 X 1 dk Ak ¼ f aj uj ; j P 0. k! dak j¼0 a¼0

Substituting (10) in Eq. (9) to get ! !! 1 1 X X 1 L un ðxÞ ¼ 2 A kL An . s n¼0 n¼0 Using the linearity property of the Laplace transform, Eq. (11) becomes ! 1 1 X X 1 Lðun ðxÞÞ ¼ 2 A k LðAn Þ . s n¼0 n¼0

ð11Þ

ð12Þ

Equating the terms yield A ; s2 k Lðu1 ðxÞÞ ¼ 2 LðA0 Þ; s k Lðu2 ðxÞÞ ¼ 2 LðA1 Þ; s .. . Lðu0 ðxÞÞ ¼

Lðun ðxÞÞ ¼

k LðAn1 Þ. s2

ð13Þ

708


Applying the inverse Laplace transform to system (13), we get u0 ðxÞ ¼ Ax; 1 x eAx u1 ðxÞ ¼ k 2 þ 2 ; A A A 5 x eAx xeAx e2Ax þ þ u2 ðxÞ ¼ k2 ; A3 4A4 4A4 2A3 A4 .. .

ð14Þ

1

From (14) the sequence fun gn¼0 has only two unknown which are A and k, so the solution of the Bratu equation is given by uðxÞ ¼ u0 ðxÞ þ u1 ðxÞ þ u2 ðxÞ þ .

ð15Þ

Using the boundary condition u(1) = 0, we get the following implicit equation: H ðA; kÞ ¼ uð1Þ ¼ u0 ð1Þ þ u1 ð1Þ þ u2 ð1Þ þ ¼ g0 ðAÞ þ kg1 ðAÞ k2 g2 ðAÞ þ k3 g3 ðAÞ þ ¼ 0;

ð16Þ

where g0 ðAÞ ¼ A; 1 1 eA g1 ðAÞ ¼ þ ; A2 A A2 5 1 eA1 eA e2A þ þ g2 ðxÞ ¼ ; 4A4 2A3 A4 A3 4A4 .. . As we described in Section 2, we can apply the predictor corrector method to trace the implicit equation H(A, k). Numerical results will presented in the next section. Now, Consider the n-dimensional Bratu problem Du þ keu ¼ 0; in X; u ¼ 0; on oX; where X = B1 is the unit ball in Rn ; n > 1. This problem is equivalent to the one-dimensional boundary value problem n1 0 u ðrÞ þ keuðrÞ ¼ 0; r u0 ð0Þ ¼ uð1Þ ¼ 0;

u00 ðrÞ þ

r 2 ð0; 1Þ;

ð17Þ

for the profile u(r) = u(jxj). For more details about the relation between these two problems, see [13]. Multiply both sides of Eq. (17) by r to get ru00 ðrÞ þ ðn 1Þu0 ðrÞ þ kreuðrÞ ¼ 0;

r 2 ð0; 1Þ;

0

u ð0Þ ¼ uð1Þ ¼ 0. Apply the Laplace transform to both sides of the last equation to get

d 2 dLðeu Þ ½s LðuÞ uð0Þs u0 ð0Þ þ ðn 1Þ½sLðuÞ uð0Þ k ¼ 0. ds ds

Substituting (10) in the last equation to get " ! # " ! # ! 1 1 1 X X d 2 X d sL un ðxÞ Bs þ ðn 1Þ sL un ðxÞ B k L An ¼ 0; ds ds n¼0 n¼0 n¼0


709

where B = u(0). Simple calculations yields to 1 1 1 X X X d d Lðun ðxÞÞ þ ðn 3Þs LðAn Þ ¼ 0. Lðun ðxÞÞ ðn 2ÞB k s2 ds ds n¼0 n¼0 n¼0 Equating the terms yield d Lðu0 ðxÞÞ þ ðn 3ÞsLðu0 ðxÞÞ ðn 2ÞB ¼ 0; ds d d s2 Lðu1 ðxÞÞ þ ðn 3ÞsLðu1 ðxÞÞ ðn 2ÞB k LðA0 Þ ¼ 0; ds ds d d s2 Lðu2 ðxÞÞ þ ðn 3ÞsLðu2 ðxÞÞ ðn 2ÞB k LðA1 Þ ¼ 0; ds ds d 2 d Lðu3 ðxÞÞ þ ðn 3ÞsLðu3 ðxÞÞ ðn 2ÞB k LðA2 Þ ¼ 0; s ds ds .. .

s2

ð18Þ

Following the same procedure as described for the one-dimensional Bratu equation, we will find an implicit equation between k and B H ðB; kÞ ¼ 0. Then, we apply the predictor corrector technique to trace the implicit solution curve generated from the equation H(B, k) = 0. More details will be presented in the next section. 4. Numerical results In this section we will present two numerical examples. In the first one, we will study the one-dimensional Bratu problem on the interval (0, 1) while we will discuss the problem in the n-dimensional domain in Example 4.2. Conclusions will be drawn. Example 4.1. Consider the following Bratu problem: u00 þ keu ¼ 0; 0 < x < 1; uð0Þ ¼ uð1Þ ¼ 0. Following the discussion in Section 3 and using the system (13) we get: A or s2 A0 ¼ eu0 ¼ eAx ; Lðu0 ðxÞÞ ¼

u0 ðxÞ ¼ Ax;

k 1 x eAx Ax Lðu1 ðxÞÞ ¼ 2 Lðe Þ or u1 ðxÞ ¼ k 2 þ 2 ; s A A A Ax 1 x e A1 ¼ u1 eu0 ¼ k 2 þ 2 eAx ; A A A k 5 x eAx xeAx e2Ax 2 Lðu2 ðxÞÞ ¼ 2 LðA1 Þ or u2 ðxÞ ¼ k þ þ 3 4 ; s A 4A 4A4 2A3 A4 2 u1 þ u2 e u 0 A2 ¼ 2 2 1 x eAx 5 x eAx xeAx e2Ax þ 2 þ þ 2 2 4 3 4 3 4 A A 4A 2A A A A 4A eAx ¼ k2 2 .. .

710


where u 0 (0) = A. The rest of components of the decomposition series were obtained using Mathematica. We used 14 terms in evaluating the approximate solution 14 X uk ðxÞ. uapp ðxÞ ¼ k¼0

Using the second boundary condition we get uapp ð1Þ ¼

14 X

uk ð1Þ ¼

14 X

k¼0

kk gk ðAÞ ¼ 0;

k¼0

where gk(A), k = 0, 1, . . . are given in Eq. (16). Let H : R2 ! R be a map which be given by H ðA; kÞ ¼

14 X

kk gk ðAÞ.

k¼0

Let k0 = 0.01. We apply the fixed point method to find an approximate value for A0 so that H(A0, k0) = 0. Direct calculations give us Aðmþ1Þ ¼

14 X

kk0 gk ðAðmÞ Þ.

k¼1

Applying the previous iterations to get A0 ’ A(9) = 0.00500417. Moreover, H is a smooth map with H(A0, k0) = 0. This means, we can use (A0, k0) as an initial point for the predictor corrector technique which is described in Section 2. In this example we will use the stepsize h = 0.05. The graphs of the exact solution and uapp(x) at k* = 3.513830718 are given in Fig. 1. We see that the technique, used in this paper, works efficiently and accurately with this problem. In Fig. 2, we sketch the graph of the maximum norm of uapp(x) against k. From this graph, we see that our technique gives us two solution when k < k*, one solution when k = k*, and no solution when k > k*. We want to mention that the approximated value for k* obtained by our technique is 3.51383077 which is correct up to eight digits. Example 4.2. Consider the two-dimensional Bratu problem Du þ keu ¼ 0; in X; u ¼ 0; on oX; where X = B1 is the unit ball in R2 . This problem is equivalent to the one-dimensional boundary value problem 1 u00 ðrÞ þ u0 ðrÞ þ keuðrÞ ¼ 0; r u0 ð0Þ ¼ uð1Þ ¼ 0;

r 2 ð0; 1Þ;

Fig. 1. The graphs of the exact solution and uapp(x) at k*.


711

Fig. 2. The graph of the maximum norm of uapp(x) against k.

for the profile u(r) = u(jxj). We will use similar procedure as we did in the last example. For n = 2, system (18) becomes d Lðu0 ðxÞÞ þ sLðu0 ðxÞÞ ¼ 0; ds d d s2 Lðu1 ðxÞÞ þ sLðu1 ðxÞÞ k LðA0 Þ ¼ 0; ds ds d 2 d Lðu2 ðxÞÞ þ sLðu2 ðxÞÞ k LðA1 Þ ¼ 0; s ds ds d d s2 Lðu3 ðxÞÞ þ sLðu3 ðxÞÞ k LðA2 Þ ¼ 0; ds ds .. . s2

ð19Þ

Using Mathematica, we can generate the terms of the approximate solution. The solution of the first equation in system (19) is Lðu0 ðxÞÞ ¼ BS which implies that u0(x) = B. Thus, A0 = eB. The second equation of system (19) gives us Lðu1 ðxÞÞ ¼ k

eB ; 4s3

which implies that u1 ðxÞ ¼ k

eB 2 x. 8

Similarly, we see that e2B 2 x; 8 e2B Lðu2 ðxÞÞ ¼ k2 ; 24s5 e2B 4 x; u2 ðxÞ ¼ k2 576 .. .

A1 ðxÞ ¼ k

712


Fig. 3. The graphs of the exact solution and uapp(x) at k* = 2.

Fig. 4. The graph of the maximum norm of uapp(x) against k.

Also, we will use 14 terms to approximate our solution. The approximate value generated by our procedure for the turning point is k = 1.99999993. The graphs of the exact solution and uapp(x) at k* = 2 are given in Fig. 3, while Fig. 4 represents graph of the maximum norm of uapp(x) against k. We note that the maximum norm of uapp(x) is uapp(0). Similar study can be done for other dimensions. From these examples, we can see that this approach works accurately and efficiently. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]

G. Adomain, Solving Frontier Problems on Physics: The Decomposition Method, Kluwer Academic Publisher, Boston, 1994. G. Adomain, A review of the decomposition method in applied mathematics, J. Math. Anal. Appl. 135 (1988) 501–544. E.L. Allgower, K. Georg, Continuation and path following, Acta Numerica, Cambridge University Press, Cambridge, 1993, pp. 1–64. E.L. Allgower, C.S. Chien, K. Georg, C.-F. Wang, Conjugate gradient methods for continuation problems, J. Comput. Appl. Math. 38 (1991) 1–16. J.P. Abbott, An efficient algorithm for the determination of certain bifurcation points, J. Comput. Appl. Math. 4 (1978) 19–27. U.M. Ascher, R.M.M. Matheij, R.D. Russell, Numerical Solution of Boundary Value Problems for Ordinary Differential Equations, Society for Industrial and Applied Mathematics, Philadelphia, PA, 1995. B.S. Attili, Efficient finite differences for parameter-dependent singularities in two-point boundary value problems, Comput. Methods Appl. Mech. Engrg. 190 (2001) 5543–5549. R.E. Bank, T.F. Chan, PLTMGC: a multi-grid continuation program for parametrized nonlinear elliptic systems, SIAM J. Sci. Statist. Comput. 7 (1986) 540–559. J.H. Bolstad, H.B. Keller, A multigrid continuation method for elliptic problems with folds, SIAM J. Sci. Statist. Comput. 7 (1986) 1081–1104. T.F. Chan, Newton-like pseudo-arclength methods for computing simple turning points, SIAM J. Sci. Statist. Comput. 5 (1984) 135– 148. E. Deeba, S.A. Khuri, S. Xie, An algorithm for solving boundary value problems, J. Comput. Phys. 159 (2000) 125–138. A.I. Fedoseyev, M.J. Friedman, E.J. Kansa, Continuation for nonlinear elliptic partial differential equations discretized by the multiquardic method, Int. J. Bifurcat. Chaos 10 (2) (2000) 481–492.


713

[13] B. Gidas, W. Ni, L. Nirenberg, Symmetry and related properties via the maximum principle, Comm. Math. Phys. 68 (1979) 209–243. [14] H.B. Keller, D.S. Cohen, Some position problems suggested by nonlinear heat generation, J. Math. Mech. 16 (1967) 1361–1376. [15] J.B. Rosen, Approximate solution and error bounds for quasilinear elliptic boundary value problems, SIAM J. Numer. Anal. 8 (1970) 80–103. [16] H. Siyyam, M. Syam, The modified trapezoidal rule for line integrals, J. Computat. Appl. Math. 84 (1997) 1–14. [17] M. Syam, H. Siyyam, Numerical differentiation of implicitly defined curves, J. Computat. Appl. Math. 108 (1999) 131–144. [18] A.M. Wazwaz, A new modification of the Adomian decomposition method for linear and nonlinear operators, Appl. Math. Comput. 122 (2001) 393–405. [19] A.M. Wazwaz, A comparison between Adomain decomposition method and Taylor series method in the series solutions, Appl. Math. Comput. 97 (1998) 37–44.