On Numerov's method for a class of strongly nonlinear two-point boundary value problems

Applied Numerical Mathematics 61 (2011) 38–52

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

On Numerov’s method for a class of strongly nonlinear two-point boundary value problems ✩ Yuan-Ming Wang a,b,∗ a

Department of Mathematics, East China Normal University, Shanghai 200241, People’s Republic of China Scientific Computing Key Laboratory of Shanghai Universities, Division of Computational Science, E-Institute of Shanghai Universities, Shanghai Normal University, Shanghai 200234, People’s Republic of China b

a r t i c l e

i n f o

Article history: Received 15 September 2009 Received in revised form 15 June 2010 Accepted 14 August 2010 Available online 18 August 2010 Keywords: Strongly nonlinear two-point boundary value problem Numerov’s method Fourth-order accuracy Monotone iterations Upper and lower solutions

a b s t r a c t The purpose of this paper is to give a numerical treatment for a class of strongly nonlinear two-point boundary value problems. The problems are discretized by fourthorder Numerov’s method, and a linear monotone iterative algorithm is presented to compute the solutions of the resulting discrete problems. All processes avoid constructing explicitly an inverse function as is often needed in the known treatments. Consequently, the full potential of Numerov’s method for strongly nonlinear two-point boundary value problems is realized. Some applications and numerical results are given to demonstrate the high efficiency of the approach. © 2010 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction Nonlinear two-point boundary value problems arise from many fields of applied sciences and have been investigated extensively in the literature both analytically and numerically (cf. [1,3,7,31–33]). In this paper, we seek high-accuracy numerical solution of the following strongly nonlinear two-point boundary value problem:

⎧ ⎨ ⎩

−

d dx

 k (u )

u (0) = α ,

du dx

 = f (x, u ),

0  x  1,

(1.1)

u (1) = β,

where α , β are given constants, and the functions f (x, u ) and k(u ) (which, in general, are nonlinear in u) are prescribed smooth functions of their respective arguments. The consideration of problem (1.1) is motivated by some heat-conduction problems and diffusion problems. For example, if the ends of a rod are kept at given temperature and the thermal conductivity is temperature dependent then the steadystate temperature distribution u (x) in the rod is governed by the above problem (1.1), and in this case, k(u ) is the thermal conductivity and f (x, u ) is the internal source that may also be temperature dependent. Various aspects of such heatconduction problems, such as the qualitative analysis of the equations and the computation of the solutions, have been

✩ This work was supported in part by the National Natural Science Foundation of China, No. 10571059, E-Institutes of Shanghai Municipal Education Commission, No. E03004, the Natural Science Foundation of Shanghai, No. 10ZR1409300 and Shanghai Leading Academic Discipline Project No. B407. Address for correspondence: Department of Mathematics, East China Normal University, Shanghai 200241, People’s Republic of China. E-mail address: [email protected].

*

0168-9274/$30.00 © 2010 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2010.08.003

Y.-M. Wang / Applied Numerical Mathematics 61 (2011) 38–52

39

investigated in the literature (cf. [6,20,25,26,29–33]). For some other applications of problem (1.1) in the theory of diffusion we see [1,10–12,20,31–34] and the references therein. As we know, it is difficult to give the analytical solution of problem (1.1), even if the function f (·, u ) is linear in u. Therefore, numerical methods are of considerable practical interest. If k(u ) ≡ 1, problem (1.1) is reduced to a semilinear problem. For solving such a semilinear problem, various numerical methods have been developed in the literature. These methods include, for example, finite difference methods in [4,7,15,18,23,43–46], Petrov-Galerkin methods in [21], variation methods in [14], shooting methods in [22,36], spline methods in [5,16,24] and multiderivative methods in [41]. In the context of finite difference discretizations, one of the well-known methods is Numerov’s method (cf. [3,23,28]). Because Numerov’s method possesses fourth-order accuracy and a compact property, it has attracted considerable attention and has been extensively used in practical computations (cf. [2,3,13,15,18,19,37,43–46]). The compact property of Numerov’s method means that the difference stencil in this method only utilizes mesh points directly adjacent to the mesh point at which a difference approximation is being made. This property makes that the boundary conditions can be easily treated in the same manner as in the standard second-order method. In other words, Numerov’s method requires only a regular threepoint difference stencil similar to that used in the standard second-order method but still possesses fourth-order accuracy. For the developments of Numerov’s method, we refer to the survey paper [4]. If k(u ) depends on u, problem (1.1) becomes more complicated due to the nonlinearity of the function k(u ). In this case, Numerov’s method cannot be directly applied, even if k(u ) is linear in u. However, in the article [44], the author showed that after a proper transformation of (1.1), one can use Numerov’s method to find indirectly an accurate numerical solution of (1.1). The numerical method proposed there consists of three steps: Step 1. Using a transformation T (see (2.1)) to transform (1.1) into a semilinear problem; Step 2. Applying Numerov’s method to solve the resulting semilinear problem; Step 3. Making use of the inverse T −1 to obtain the numerical solution of problem (1.1). Although the above method possesses the same accuracy as Numerov’s method regardless of the strong nonlinearity of problem (1.1), it is only valid for the situation when the inverse T −1 can be constructed explicitly. This limits its applications since the construction of the inverse T −1 is not always easy to do. The purpose of this paper is to search a new technique for developing a direct method that avoids constructing the inverse T −1 and directly offers the numerical solution of (1.1) but still maintains the same accuracy as Numerov’s method. Specifically, our approach here is to formulate (1.1) as a coupled system of a semilinear two-point boundary value problem and a nonlinear functional equation and then to discretize the coupled system by Numerov’s method. To solve the resulting discrete problem we develop a linear monotone iterative algorithm by the method of upper and lower solutions and its associated monotone iterations. This approach makes it possible to compute accurately the numerical solution of (1.1) by using Numerov’s method without finding the inverse T −1 . Consequently, the full potential of Numerov’s method for strongly nonlinear two-point boundary value problem (1.1) is realized. The outline of the paper is as follows: In the next section, we formulate the strongly nonlinear problem (1.1) as a coupled system of a semilinear two-point boundary value problem and a nonlinear functional equation and then apply Numerov’s method to discretize the coupled system. Section 3 is devoted to a linear monotone iterative algorithm for the resulting discrete problem. It is shown that using an upper solution and a lower solution as initial iterations the iterative algorithm yields two sequences which converge monotonically from above and below, respectively, to a unique solution of the resulting discrete problem. The convergence of the method is investigated in Section 4 where the fourth-order accuracy of the method is proven. In Section 5, we give some applications to several model problems arising from heat-conduction, population dynamics and chemical engineering. Some numerical results are presented to demonstrate the monotone property of the iterative algorithm and the high accuracy of the numerical solution. The final section is for some concluding remarks. 2. Numerov’s method Let I be an interval in R such that α , β ∈ I . Throughout this paper, we assume that the nonlinear functions f (·, u ) and k(u ) satisfy the following basic hypothesis:

(H )

f (·, u ) and k(u ) are C 1 -functions of u , and there exists a positive constant k0 such that k(u )  k0 > 0 for all u ∈ I .

In order to apply Numerov’s method to (1.1), we form the problem as a coupled system of a semilinear two-point boundary value problem and a nonlinear functional equation. Let δ ∈ I . Define

u v = T (u ) =

k(s) ds,

∀u ∈ I .

δ

Using the transformation T we transform problem (1.1) into the following coupled system:

(2.1)

40

Y.-M. Wang / Applied Numerical Mathematics 61 (2011) 38–52

⎧ ⎨ ⎩

−

d2 v dx2



= f (x, u ),

v (0) = T (α ),



T u (x) = v (x),

0  x  1,

(2.2)

v (1) = T (β).

It is clear that u (x) is a solution of (1.1) in I (i.e., u (x) ∈ I for each x) then (u (x), v (x)) = (u (x), T (u (x))) is a solution of the coupled system (2.2) and vice versa. As compared with the treatment in [44], a new feature of the above system is that it avoids the inverse T −1 . This makes it possible to compute accurately the numerical solution of (1.1) without finding the inverse T −1 . Let γ (x) be any nonnegative function satisfying

γ (x)  − f u (x, u )/k0 for all 0  x  1 and u ∈ I ,

(2.3)

where f u = d f /du. Then we have

γ (x)k(u ) + f u (x, u )  0 for all 0  x  1 and u ∈ I .

(2.4)

In other words, the function

g (x, u ) = γ (x) T (u ) + f (x, u )

(0  x  1)

(2.5)

is monotone nondecreasing in u for u ∈ I . By adding γ (x) v (x) on both sides of the first equation in (2.2) and using the relation (2.1), we obtain the following equivalent system:

⎧ ⎨ ⎩

−

d2 v dx2



+ γ (x) v (x) = g (x, u ),

v (0) = T (α ),



T u (x) = v (x),

0  x  1,

(2.6)

v (1) = T (β).

Our Numerov’s scheme for the numerical solution of (2.2) is based on the above equivalent form (2.6) (for the reason why we don’t discretize (2.2) directly, see Remark 2.1). Let h = 1/ N be the mesh size, and let xi = ih (0  i  N ) be the mesh points. For convenience, we introduce the finite difference operators δh2 and Ph as follows:

δh2 u (xi ) = u (xi −1 ) − 2u (xi ) + u (xi +1 ), h  2

Ph u (xi ) =

12

1  i  N − 1,



u (xi −1 ) + 10u (xi ) + u (xi +1 ) ,

1  i  N − 1.

(2.7)

Using the following Numerov’s formula (cf. [4,28])

 δh2 u (xi ) = Ph u (2) (xi ) + O h6 , where u (2) (xi ) =



2

d u ( xi ) , dx2

1  i  N − 1,

(2.8)

we have from (2.6) that



  −δh2 v (xi ) + Ph γ (xi ) v (xi ) = Ph g xi , u (xi ) + O h6 , v (x0 ) = T (α ),





T u (xi ) = v (xi ),

v (x N ) = T (β).

1  i  N − 1,

(2.9)

After dropping the O (h6 ) term, we derive a finite difference scheme of (2.6) as follows:



−δh2 v i + Ph (γi v i ) = Ph g i (u i ), v 0 = T (α ),

T (u i ) = v i ,

1  i  N − 1,

v N = T (β),

(2.10)

where u i and v i represent the approximations of u (xi ) and v (xi ), respectively, and

γi = γ (xi ),

g i (u i ) = g (xi , u i ).

To write (2.10) in a vector form, we let J = ( J i , j ) and B = ( B i , j ) be the ( N − 1)-order symmetric tridiagonal matrices with the following elements:

J i , i = 2, B i , i = 5 /6 ,

J i ,i −1 = J i ,i +1 = −1, B i ,i −1 = B i ,i +1 = 1/12,

1  i  N − 1,

and let Γ be the diagonal matrix given by

Γ = diag(γ1 , γ2 , . . . , γ N −1 ). Also we define the following ( N − 1)-dimensional column vectors:

(2.11)

Y.-M. Wang / Applied Numerical Mathematics 61 (2011) 38–52

U = (u 1 , u 2 , . . . , u N −1 )T ,



41

V = ( v 1 , v 2 , . . . , v N −1 )T ,

T

G (U ) = g 1 (u 1 ), g 2 (u 2 ), . . . , g N −1 (u N −1 ) ,



T

T (U ) = T (u 1 ), T (u 2 ), . . . , T (u N −1 ) ,











T

H = T (α ) + h2 g (0, α ) − γ0 T (α ) /12, 0, . . . , 0, T (β) + h2 g (1, β) − γ N T (β) /12 .

(2.12)

Then system (2.10) can be written in the vector form:





J + h 2 B Γ V = h 2 B G (U ) + H ,

(2.13)

T (U ) = V .

This is a coupled system of a nonlinear algebraic equation and a nonlinear functional equation. To analyze it and develop a linear monotone iterative algorithm for its solution, we review some well-known properties about the matrices J and B. For two constants M and M satisfying M  M > −π 2 we define

⎧ 12 ⎪ ⎪ ⎪ M, ⎪ ⎪ ⎪ ⎪ ⎨ 1, h( M , M ) = 12 12 ⎪ min { , (1 + πM2 ) }, ⎪ ⎪ π2 M ⎪ ⎪ ⎪ ⎪ ⎩ 12 (1 + M ), π2 π2

M > −8, M > 0, M > −8, M  0, (2.14)

M  −8, M > 0, M  −8, M  0.

Lemma 2.1. (See Lemma 2.3 of [45] or Lemma 3.1 of [46].) Let M = diag( M 1 , . . . , M N −1 ) be a diagonal matrix and let M and M be two constants such that

−π 2 < M  M i  M ,

i = 1, 2, . . . , N − 1.

(2.15)

Then the inverse ( J + h2 B M )−1 exists and is nonnegative provided that h < h( M , M ). The following lemma will be used in Section 4 to prove the convergence of the resulting finite difference scheme (2.10) (or (2.13)). Lemma 2.2. (See Lemma 2.4 of [45].) Let M, M and M i (i = 1, 2, . . . , N − 1) be the given constants satisfying (2.15), and let M = diag( M 1 , . . . , M N −1 ). Also let Z and R be the vectors in R N −1 such that





J + h2 B M Z = R .

(2.16)

Then when h < h( M , M ),

⎧ −2 ⎪ M  0, ⎨ ( R ∞ /8)h ,  Z ∞  ( R ∞ /(8 + M ))h−2 , −8 < M < 0, ⎪ ⎩ ( R ∞ /(2π (1 − θ)))h−2 , −π 2 < M  −8, h  H ,

(2.17)

where H < h( M , M ) and

θ =−

M

π 2 (1 − π 2 H 2 /12)

< 1.

Remark 2.1. The formulation of the system (2.13) is based on the equivalent form (2.6). This system can be clearly written as



J V = h2 BF (U ) + H ,

(2.18)

T (U ) = V , where



T

F (U ) = f 1 (u 1 ), f 2 (u 2 ), . . . , f N −1 (u N −1 ) ,

f i (u i ) = f ( xi , u i )

(i = 1, 2, . . . , N − 1).

(2.19)

System (2.18) is the discretization of the original problem (2.2) by Numerov’s method. As compared with (2.18), an important property of (2.13) is that the function G (U ) in this system is monotone nondecreasing in U (see (2.4)). This property, as we shall see later, makes it more convenient to design a linear monotone iterative algorithm for the computation of the solutions. This is our main reason for considering (2.13) instead of (2.18).

42

Y.-M. Wang / Applied Numerical Mathematics 61 (2011) 38–52

Remark 2.2. Much research has been done for the numerical solutions of the initial value problems (IVPs) related to the second-order ordinary differential equations. Various numerical methods have appeared in the literature, including the linear multistep methods and hybrid methods which are often used in practical applications (see [8,9,17,27,35,37–40,42]). A detailed survey for the recent developments of these methods is given in [9]. It is interesting to examine if the techniques in these methods can be applied to present two-point boundary value problems. Further investigations will be made in the future. 3. Linear monotone iterative algorithm To develop a linear monotone iterative algorithm for the solutions of (2.13) we use the method of upper and lower solutions and its associated monotone iterations. A pair of ordered upper and lower solutions of (2.13) are defined as follows: Definition 3.1. A vector (U˜ , V˜ ) ∈ R N −1 × R N −1 is called an upper solution of (2.13) if





J + h2 B Γ V˜  h2 B G (U˜ ) + H ,

(3.1)

T (U˜ )  V˜ .

Similarly, (Uˆ , Vˆ ) ∈ R N −1 × R N −1 is called a lower solution of (2.13) if it satisfies the above inequalities in the reversed order. A pair of upper and lower solutions (U˜ , V˜ ) and (Uˆ , Vˆ ) are said to be ordered if (U˜ , V˜ )  (Uˆ , Vˆ ). In the above definition, inequalities between vectors are in the sense of componentwise. It is clear that every solution of (2.13) is an upper solution as well as a lower solution. For any W in R N −1 , we denote w i the i-th component of W . Given a pair of ordered upper and lower solutions (U˜ , V˜ ) and (Uˆ , Vˆ ), we define the sectors

  S = (U , V ) ∈ R N −1 × R N −1 ; (Uˆ , Vˆ )  (U , V )  (U˜ , V˜ ) , uˆ i , u˜ i = {u i ∈ R; uˆ i  u i  u˜ i },

(3.2)

and define the interval I0 by

  I0 = min uˆ i , max u˜ i . i

i

We assume that I0 ⊂ I throughout the paper, where I is defined in hypothesis ( H ). To compute the solutions of (2.13) we use the following linear iterative scheme:









J + h2 B Γ V (m+1) = h2 B G U (m) + H ,

 Γ (m) U (m+1) = Γ (m) U (m) − T U (m) + V (m+1) ,

(3.3)

m = 0, 1 , 2 , . . . ,

where the initial iteration (U (0) , V (0) ) is either (U˜ , V˜ ) or (Uˆ , Vˆ ), and

 (m) (m) Γ (m) = diag γ1 , . . . , γ N −1 , (m)

The functions u i





γi(m) = max k(s); u (i m)  s  u (i m) .

(3.4)

s

, u (i m) in the definition of γi(m) are the respective components of U (m) and U (m) which are obtained from

(3.3) with (U (0) , V (0) ) = (U˜ , V˜ ) and (U (0) , V (0) ) = (Uˆ , Vˆ ), respectively.

To show that the sequences given by (3.3) are well-defined it is crucial that the sequences {U (m) } and {U (m) } possess the property U (m)  U (m) for every m. It is clear from (3.4) that

  Γ (m) U − U  − T (U ) + T U   0 whenever U (m)  U   U  U (m) .

(3.5)

Moreover, by the monotone nondecreasing property of g (x, u ) in u,



G (U )  G U 



whenever Uˆ  U   U  U˜ .

(3.6) that the inverse ( J + h B Γ )−1 exists and is a nonnegative 2

We observe from Lemma 2.1 and the nonnegative property of Γ matrix provided that h < h(0, Γ ∞ ). These properties lead to the following well-defined and monotone properties of the sequences from (3.3). Lemma 3.1. Let hypothesis ( H ) be satisfied and let h < h(0, Γ ∞ ). Also let (U˜ , V˜ ) and (Uˆ , Vˆ ) be a pair of ordered upper and lower solutions of (2.13). Then the sequences {(U (m) , V (m) )}, {(U (m) , V (m) )} and {Γ (m) } given by (3.3) and (3.4) with (U (0) , V (0) ) = (U˜ , V˜ )

Y.-M. Wang / Applied Numerical Mathematics 61 (2011) 38–52

43

and (U (0) , V (0) ) = (Uˆ , Vˆ ) are all well-defined and possess the monotone property

    (Uˆ , Vˆ )  U (m) , V (m)  U (m+1) , V (m+1)  U (m+1) , V (m+1)  U (m) , V (m)  (U˜ , V˜ )

(3.7)

for every m = 0, 1, 2, . . . . Proof. Let m = 0 in (3.3). Since U (0) = U˜ , U (0) = Uˆ and U˜  Uˆ , the elements

γi(0) of the diagonal matrix Γ (0) are well-

defined. By hypothesis ( H ), the inverse (Γ (0) )−1 exists and is a nonnegative matrix. Hence the first iterations and (U (1) , V (1) ) exist and satisfy the relation

(U (1) , V (1) )

⎧   ⎪ J + h2 B Γ V (0) − V (1) = J + h2 B Γ V˜ − h2 B G (U˜ ) − H , ⎪ ⎨   J + h2 B Γ V (1) − V (0) = h2 B G (Uˆ ) + H − J + h2 B Γ Vˆ , ⎪ ⎪   ⎩ J + h2 B Γ V (1) − V (1) = h2 B G (U˜ ) − G (Uˆ ) .

(3.8)

By (3.1) and (3.6), the right-hand side of the above relation is nonnegative, and thus by the nonnegative property of ( J + h2 B Γ )−1 , V (0)  V (1)  V (1)  V (0) . Knowing this relation we have from (3.1), (3.3) and (3.5) that

 Γ (0) U (0) − U (1) = T (U˜ ) − V (1)  T (U˜ ) − V˜  0,  Γ (0) U (1) − U (0) = − T (Uˆ ) + V (1)  − T (Uˆ ) + Vˆ  0,  Γ (0) U (1) − U (1) = Γ (0) (U˜ − Uˆ ) − T (U˜ ) + T (Uˆ ) + V˜ − Vˆ  0. In view of (Γ (0) )−1  0, U (0)  U (1)  U (1)  U (0) . This proves that (3.7) holds for m = 0. Assume, by induction, that for some m  1, (U (m) , V (m) ), (U (m) , V (m) ), (U (m−1) , V (m−1) ) and (U (m−1) , V (m−1) ) exist and (3.7) holds when m is replaced by m − 1. It follows from Uˆ  U (m)  U (m)  U˜ that Γ (m) is well-defined and its inverse (Γ (m) )−1 is nonnegative. This ensures that (U (m+1) , V (m+1) ) and (U (m+1) , V (m+1) ) exist. By (3.3) and (3.6),

  

J + h2 B Γ



2



2



J + h BΓ J + h BΓ



 







(m−1)



V (m) − V (m+1) = h2 B G U (m−1) − G U (m) V

(m+1)

V

(m+1)



−V

(m)

−V

(m+1)

2

 

(m)





 0,

=h B G U −G U  0,   (m)  (m) 2 =h B G U −G U  0.

The nonnegative property of ( J + h2 B Γ )−1 leads to that V (m)  V (m+1)  V (m+1)  V (m) . This relation together with (3.3) and (3.5) yields

    Γ (m) U (m) − U (m+1)  Γ (m−1) U (m−1) − U (m) − T U (m−1) + T U (m)  0,     Γ (m) U (m+1) − U (m)  Γ (m−1) U (m) − U (m−1) − T U (m) + T U (m−1)  0,     Γ (m) U (m+1) − U (m+1)  Γ (m) U (m) − U (m) − T U (m) + T U (m)  0 which implies that U (m)  U (m+1)  U (m+1)  U (m) . The conclusion of the lemma follows from the principle of induction.

2

In view of the monotone property (3.7), the limits





lim U (m) , V (m) = (U , V ),

m→∞





lim U (m) , V (m) = (U , V )

m→∞

(3.9)

exist and satisfy

    (Uˆ , Vˆ )  U (m) , V (m)  U (m+1) , V (m+1)  (U , V )  (U , V )  U (m+1) , V (m+1)  U (m) , V (m)  (U˜ , V˜ ), m = 0, 1, 2, . . . .

(3.10)

The following theorem shows that the limits (U , V ) and (U , V ) are the maximal and minimal solutions of (2.13) in S , respectively, in the sense that if (U , V ) is any solution of (2.13) in S then (U , V )  (U , V )  (U , V ). Theorem 3.1. Let the conditions in Lemma 3.1 hold. Then the sequences {(U (m) , V (m) )} and {(U (m) , V (m) )} converge monotonically to the maximal solution (U , V ) and the minimal solution (U , V ) of (2.13) in S , respectively. Moreover, the relation (3.10) holds. Proof. Let

γ i = maxs {k(s); u i  s  u i } and Γ = diag(γ 1 , . . . , γ N −1 ), where u i and u i are the respective components of

the limits U and U in (3.9). Then by the definition of Γ (m) , Γ (m)  Γ (m+1)  Γ for every m, and so the sequence {Γ (m) } converges as m → ∞. Letting m → ∞ in (3.3) shows that both (U , V ) and (U , V ) are solutions of (2.13) in S . We next show

44

Y.-M. Wang / Applied Numerical Mathematics 61 (2011) 38–52

that these solutions are the maximal and minimal solutions of (2.13) in S . Let (U , V ) be any solution of (2.13) in S . Then by (2.13) and (3.3),



J + h2 B Γ





 





V (m+1) − V = h2 B G U (m) − G (U ) ,

   Γ (m) U (m+1) − U = Γ (m) U (m) − U − T U (m) + T (U ) + V (m+1) − V , Since (U , V ) ∈ S , the above relation for m = 0 gives



J + h2 B Γ







m = 0, 1 , 2 , . . . .

(3.11)



V (1) − V = h2 B G (U˜ ) − G (U )  0,

 Γ (0) U (1) − U = Γ (0) (U˜ − U ) − T (U˜ ) + T (U ) + V (1) − V  V (1) − V .

(3.12)

This yields (U (1) , V (1) )  (U , V ) by using the nonnegativity of ( J + h2 B Γ )−1 and (Γ (0) )−1 . A similar argument gives (U (1) , V (1) )  (U , V ). It follows from an induction argument that (U (m) , V (m) )  (U , V )  (U (m) , V (m) ) for every m. This shows (U , V )  (U , V )  (U , V ), which proves the maximal and minimal property of (U , V ) and (U , V ). 2 It is obvious from the maximal and minimal property of (U , V ) and (U , V ) that if (U , V ) = (U , V ) (≡ (U ∗ , V ∗ )) then S . To ensure (U , V ) = (U , V ), it is necessary to impose some additional conditions on f (·, u ). To give such a sufficient condition, we define

(U ∗ , V ∗ ) is the unique solution of (2.13) in

  f u ( xi , ξ ) , − k(η) ξ,η∈uˆ i ,u˜ i

σ = max max i

σ = min i

  f u ( xi , ξ ) − k(η) ξ,η∈uˆ i ,u˜ i min

(3.13)

and assume

σ > −π 2 ,

h < h(σ , σ ).

(3.14)

Theorem 3.2. Let the conditions in Lemma 3.1 hold. If, in addition, condition (3.14) holds, then the sequences {(U (m) , V (m) )} and {(U (m) , V (m) )} converge monotonically from above and below, respectively, to the unique solution (U ∗ , V ∗ ) of (2.13) in S . Proof. It suffices to show (U , V ) = (U , V ), where (U , V ) and (U , V ) are the respective limits of the sequences {(U (m) , V (m) )} and {(U (m) , V (m) )}. By (2.13) (or (2.18)),





J ( V − V ) = h 2 B F (U ) − F (U ) ,

(3.15)

T (U ) − T (U ) = V − V . Applying the mean-value theorem we obtain





T (U ) − T (U ) = T u Θ  (U − U ),

F (U ) − F (U ) = F u (Θ)(U − U ),

(3.16)

where





F u (Θ) = diag f u (x1 , θ1 ), . . . , f u (x N −1 , θ N −1 ) ,



(θ1 , θ2 , . . . , θ N −1 ) ∈ Uˆ , U˜ ,





 





T u Θ  = diag k θ1 , . . . , k θ N −1 ,

θ1 , θ2 , . . . , θ N −1 ∈ Uˆ , U˜ .

(3.17)

Using the relation (3.16) in (3.15) leads to







J − h2 BF u (Θ) T u Θ 

−1

( V − V ) = 0.

(3.18)

The condition (3.14) implies that the matrix − F u (Θ)( T u (Θ  ))−1 satisfies the condition of Lemma 2.1, and therefore by this lemma, the inverse ( J − h2 BF (Θ)( T (Θ  ))−1 )−1 exists. It follows from (3.18) that V = V . This implies T (U ) − T (U ) = 0 u

u

which ensures U = U . This proves (U , V ) = (U , V ).

2

Remark 3.1. (a) If k(u ) ≡ 1 then all the conclusions in Theorems 3.1 and 3.2 hold true for the semilinear problem (1.1). In this situation, the iteration process (3.3) is reduced to that in [43] and the uniqueness condition (3.14) becomes the one in [46]. (b) For each m, the iterations V (m) and U (m) from (3.3) can be computed, respectively, by solving a tridiagonal system of linear algebraic equations and a diagonal system of linear algebraic equations, and thus (3.3) provides a linear monotone iterative algorithm. Remark 3.2. We see from the monotone property (3.10) that for each m, (U (m) , V (m) ) is an upper bound of the maximal solution (U , V ), while (U (m) , V (m) ) gives a lower bound of the minimal solution (U , V ). Moreover, these bounds are improved, step-by-step, as m increases. If (U , V ) = (U , V ) (≡ (U ∗ , V ∗ )) then they become improved upper and lower bounds of (U ∗ , V ∗ ). This demonstrates the superiority of the monotone convergence over the ordinary convergence.

Y.-M. Wang / Applied Numerical Mathematics 61 (2011) 38–52

45

The linear iterative scheme (3.3) can be replaced by the following form









J + h2 B Γ V (m+1) = h2 B G U (m) + H ,

 Γ (0) U (m+1) = Γ (0) U (m) − T U (m) + V (m) ,

(3.19)

m = 0, 1 , 2 , . . .

which results in less computations in the process of iterations but still maintains the monotone convergence (3.10) of the sequences. The latter result follows the same line as in the proofs of Lemma 3.1 and Theorem 3.1. However, the sequence given by (3.19) converges slower than the one by (3.3). Specifically, we have the following comparison result. Theorem 3.3. Let the conditions in Lemma 3.1 hold, and let (U˜ , V˜ ) and (Uˆ , Vˆ ) be a pair of ordered upper and lower solutions of (2.13). Denote by {(U (m) , V (m) )} and {(U (m) , V (m) )} the sequences from (3.3) with (U (0) , V (0) ) = (U˜ , V˜ ) and (U (0) , V (0) ) = (Uˆ , Vˆ ), and by

{(U ∗ , V ∗ )} and {(U ∗ , V ∗ )} the sequences from (3.19) with (U ∗ , V ∗ ) = (U˜ , V˜ ) and (U ∗ , V ∗ ) = (Uˆ , Vˆ ). Then we have (m)

(m)



(m)



(m)

(0)







U (m) , V (m)  U (∗m) , V (∗m) ,





(0)



U (m) , V (m)  U (∗m) , V (∗m) ,

(0)

(0)

m = 0, 1, 2, . . . .

(3.20)

Proof. By (3.3) and (3.19),



J + h2 B Γ





 







V (∗m+1) − V (m+1) = h2 B G U (∗m) − G U (m) ,

(3.21)

and

    Γ (0) U (∗m+1) − U (m+1) = Γ (0) U (∗m) − U (m) + Γ (0) − Γ (m) U (m) − U (m+1)   − T U (∗m) + T U (m) + V (∗m) − V (m+1) .

(3.22)

Since Γ (0)  Γ (m) , U (m)  U (m+1) and V (m)  V (m+1) , the equality (3.22) implies

    Γ (0) U (∗m+1) − U (m+1)  Γ (0) U (∗m) − U (m) − T U (∗m) + T U (m) + V (∗m) − V (m) .

(3.23)

In view of (U (0) , V (0) ) = (U ∗ , V ∗ ) = (U˜ , V˜ ), we have from the above relations (3.21) and (3.23) for m = 0 that (0)



J + h2 B Γ



(0)



 Γ (0) U (∗1) − U (1)  0.

V (∗1) − V (1) = 0,

The nonnegative property of ( J + h2 B Γ )−1 and (Γ (0) )−1 yields that V (1) = V ∗ and (3.23) we get from (3.5) and (3.6) that

(1)



J + h2 B Γ





(2)





 Γ (0) U (∗2) − U (2)  0.

V (∗2) − V (2)  0,

This proves V (2)  V ∗



(1)

and U (1)  U ∗ . Using this result in (3.21)

(2)

and U (2)  U ∗ . An induction argument shows that



U (m) , V (m)  U (∗m) , V (∗m) ,

m = 0, 1 , 2, . . . .

Another inequality in (3.20) can be proved in the same manner.

2

Remark 3.3. The comparison result in Theorem 3.3 implies that with the same initial iteration, which is either an upper solution or a lower solution, the sequence given by (3.3) converges faster than the one by (3.19). 4. Convergence of the method In this section, we deal with the convergence of finite difference scheme (2.10) (or (2.13)), and show its fourth-order accuracy. Let (u (xi ), v (xi )) be the value of the solution of (2.2) at xi , and let (u i , v i ) stand for the solution of (2.10). Since, by (2.2), (2.5) and (2.10),









g xi , u ( xi ) = γ ( xi ) v ( xi ) + f xi , u ( xi ) ,

g i (u i ) = γi v i + f (xi , u i ),

we have from (2.9) and (2.10) that the solutions (u (xi ), v (xi )) and (u i , v i ) satisfy, respectively, the systems



  −δh2 v (xi ) = Ph f xi , u (xi ) + O h6 , v (x0 ) = T (α ),

and





1  i  N − 1,

v (x N ) = T (β),

−δh2 v i = Ph f (xi , u i ), v 0 = T (α ),



T u (xi ) = v (xi ),

T (u i ) = v i ,

1  i  N − 1,

v N = T (β).

Based on (4.1) and (4.2), we show the convergence of (u i , v i ) to (u (xi ), v (xi )) as h → 0.

(4.1)

(4.2)

46

Y.-M. Wang / Applied Numerical Mathematics 61 (2011) 38–52

Theorem 4.1. Assume that u i , u (xi ) ∈ Λ ⊂ I for each i and some interval Λ in R, and assume

f u (x, ξ )

< π 2,

k(η)

0  x  1, (ξ, η) ∈ Λ × Λ.

(4.3)

Then for sufficiently small h,







maxu (xi ) − u i   ch4 ,



max v (xi ) − v i   ch4 ,

i

(4.4)

i

where c is a positive constant independent of h. Proof. Let e i = u (xi ) − u i and e i = v (xi ) − v i . A subtraction of the corresponding equations in the systems (4.1) and (4.2) and using the mean-value theorem lead to

  −δh2 e i = Ph f u (xi , ξi )e i + O h6 ,

e

0

= e

N



k ξi e i = e i ,

1  i  N − 1,

= 0,

(4.5)

where ξi and ξi are two intermediate values between u (xi ) and u i . Define

T



E  = e 1 , e 2 , . . . , e N −1 ,

E = (e 1 , e 2 , . . . , e N −1 )T ,





F u (Ξ ) = diag f u (x1 , ξ1 ), f u (x2 , ξ2 ), . . . , f u (x N −1 , ξ N −1 ) ,





  





T u Ξ  = diag k ξ1 , k ξ2 , . . . , k ξ N −1 . Then the vector form of (4.5) is given by



JE = h2 BF u (Ξ ) E + O h6 ,





(4.6)

T u Ξ  E = E . This yields







J − h2 BF u (Ξ ) T u Ξ 

−1



E  = O h6 .

(4.7)

By (4.3), the matrix − F u (Ξ )( T u (Ξ  ))−1 satisfies the condition (2.15). An application of Lemma 2.2 to (4.7) shows that for sufficiently small h,  E  ∞  c 1 h4 where c 1 is a positive constant independent of h. Hence  E ∞  c 1 ( T u (Ξ  ))−1 ∞ h4  1 4 c 1 k− 0 h . This proves the estimate (4.4). 2 Remark 4.1. (a) Theorem 4.1 shows that as h → 0, the solution (u i , v i ) of (2.10) converges to the solution (u (xi ), v (xi )) of (2.2) with the convergence order O (h4 ) in the L ∞ -norm. (b) If k(u ) ≡ 1, the convergence result in Theorem 4.1 coincides with that in [45] for the semilinear problem. 5. Applications and numerical results In this section, we apply the results of the previous sections to several model problems arising from heat-conduction, population dynamics and chemical engineering. We present some numerical results to demonstrate the monotone property of the linear iterative algorithm (3.3) and the high accuracy of the numerical solution. All computations are carried out by using a MATLAB subroutine on a Pentium-4 computer with 2G memory. Since the solution of (1.1) is usually nonnegative for practical problems, we assume the boundary values α  0 and β  0 in this section. Also we assume that the functions f (·, u ) and k(u ) are C 1 -functions of u, and there exists a positive constant ρ  max{α , β} such that k(u )  k0 > 0 for all u ∈ [0, ρ ]. This implies that the functions f (x, u ) and k(u ) satisfy hypothesis ( H ) with the set I = [0, ρ ]. For this situation, we always take δ = 0 in the definition of transformation T . Under these assumptions, the zero vector (Uˆ , Vˆ ) = (0, 0) ∈ R N −1 × R N −1 is a lower solution of (2.13) if

f (0, α )  0,

f (1, β)  0,

f (x, 0)  0 (0 < x < 1).

(5.1)

To construct an upper solution of (2.13), we let (U˜ , V˜ ) = (ρ E , T (ρ )E ), where E = (1, 1, . . . , 1)T ∈ R N −1 . A simple calculation shows that this vector is an upper if

JT (ρ )E  h2 BF (ρ E ) + H , where F (u ) is defined by (2.19). The above inequality is satisfied by any

ρ satisfying

Y.-M. Wang / Applied Numerical Mathematics 61 (2011) 38–52

⎧ f (x, ρ )  0 (0  x  1), ⎪ ⎪ ⎨ ρ

α ⎪ ⎪ ⎩

k(s) ds  max 0



β 2

(5.2)

2

k(s) ds + h f (0, α )/12, 0

47

k(s) ds + h f (1, β)/12 . 0

We now consider three examples to illustrate the applications of the above construction and the results in the previous sections. Example 5.1. Consider the heat transfer in a rod with the given temperatures α  0 and β  0 at two ends. If the thermal conductivity is considered temperature dependent and the form of the source function is based on the so-called Boltzmann fourth-power law, the equation governing the steady-state temperature distribution u (x) in the rod is given by (1.1) with





f (x, u ) = σ a4 (x) − u 4 (x) ,

(5.3)

where σ is a positive constant and a(x)  0 denotes the surrounding temperature (see [20,30]). Assume that a(0)  α and a(1)  β , and let a  max0x1 a(x). Since f u (x, u ) = −4σ u 3 , the function γ (x) in (2.3) can be chosen as γ (x) ≡ γ = 4σρ 3 /k0 . In this case, the function g (x, u ) in (2.5) and the difference scheme (2.10) are, respectively, reduced to

g (x, u ) = γ T (u ) + f (x, u ), and



(5.4)

−δh2 v i + γ Ph v i = Ph g i (u i ), v 0 = T (α ),

T (u i ) = v i ,

1  i  N − 1,

(5.5)

v N = T (β).

The vector form of (5.5) is given by





J + γ h2 B V = h2 BG(U ) + H ,

(5.6)

T (U ) = V ,

where H = ( T (α ) + h2 ( g (0, α ) − γ T (α ))/12, 0, . . . , 0, T (β) + h2 ( g (1, β) − γ T (β))/12)T . Since a(0)  α , a(1)  β and a(x)  0, the condition (5.1) is satisfied by the function f (x, u ) in (5.3). This ensures that (Uˆ , Vˆ ) = (0, 0) is a lower solution of (5.6). On the other hand, the condition (5.2) for the present problem holds if

⎧ ρ  a, ⎪ ⎪ ⎨ ρ

α   β  4  4 2 4 2 4 ⎪ k(s) ds  max k(s) ds + h σ a (0) − α /12, k(s) ds + h σ a (1) − β /12 . ⎪ ⎩ 0

0

(5.7)

0

This requirement can be satisfied by a sufficiently large ρ or by taking mesh size h suitably small. Therefore, (U˜ , V˜ ) = (ρ E , T (ρ )E ) and (Uˆ , Vˆ ) = (0, 0) form a pair of ordered upper and lower solutions of (5.6). With respect to this pair, the constants σ and σ in (3.13) is given by

σ = 4σρ 3 max k−1 (η), 

σ = 0.

η∈0,ρ

(5.8)



Let h < min{ 12/γ , 12/σ }. Using (U (0) , V (0) ) = (ρ E , T (ρ )E ), (U (0) , V (0) ) = (0, 0) and the functions in (5.3) and (5.4) in the linear iterative scheme (3.3) with Γ = γ I (where I denotes the identity matrix), we can compute the sequences {(U (m) , V (m) )} and {(U (m) , V (m) )}. By Theorem 3.2, these sequences converge monotonically to the unique solution (U ∗ , V ∗ ) of (5.6) in S where

   S = (U , V ) ∈ R N −1 × R N −1 ; (0, 0)  (U , V )  ρ E , T (ρ )E .

(5.9)

Since f u (x, u )  0 for all u  0 and all 0  x  1, the condition (4.3) is trivially satisfied. By Theorem 4.1, the solution (u ∗i , v ∗i ) of (5.5) converges to the continuous solution (u ∗ (xi ), v ∗ (xi )) of (2.2), (5.3) at every mesh point xi with the convergence order O (h4 ), as h → 0. To give some numerical results, we consider a special case of the above model problem where α = β = 0, σ = 1 and

k (u ) =

1 1 + u3



,

1/4

π 2 (sin π x + sin4 π x + 3 sin2 π x cos2 π x) a(x) = + sin4 π x (1 + sin3 π x)2

.

(5.10)

48

Y.-M. Wang / Applied Numerical Mathematics 61 (2011) 38–52

Fig. 1. The monotone convergence of the sequences at xi = 0.5 for Example 5.1.

Table 1 Numerical solution (U ∗ , V ∗ ) of Example 5.1. xi

u ∗i

v ∗i

0.1 0.2 0.3 0.4 0.5

0.30901814073003 0.58778747153044 0.80901821366068 0.95105728817513 1.00000076672804

0.30677611425942 0.56098264809778 0.72574878708737 0.81026505270532 0.83564923162852

It can be checked that u = sin π x is a solution of this example, and the transformation T in (2.1) is given by

u v = T (u ) = 0

√ 2u − 1 (1 + u )2 1 3π ds = ln 2 + √ arctan √ + . 3 6 18 1+s u −u+1 3 3 1

1

(5.11)

We see that it is very difficult to obtain the inverse T −1 . Since 0  a(x)  1.8 for all 0  x  1, the constants appeared above may be taken as





k 0 = 1/ 1 + ρ 3 ,

a = ρ = 1.8,





γ = σ = 4ρ 3 1 + ρ 3 .

(5.12)

With Γ = γ I in the linear iterative scheme (3.3), we compute the sequences {(U (m) , V (m) )} and {(U (m) , V (m) )} with (U (0) , V (0) ) = (ρ E , T (ρ )E ) and (U (0) , V (0) ) = (0, 0) for various values of h. In all the numerical computations, the basic feature of monotone convergence of the sequences was observed. Let h = 1/40. In Fig. 1, we present some numerical results of these sequences at xi = 0.5, where the solid line denotes the sequence {U (m) } or { V (m) } and the dash-dot line stands for the sequence {U (m) } or { V (m) }. As expected from our theoretical analysis in Theorem 3.2, the monotone convergence of the sequences holds and their common limit (U ∗ , V ∗ ) is the unique solution of (5.6) in S given by (5.9). More numerical results of (U ∗ , V ∗ ) at various xi are explicitly given in Table 1. To demonstrate the accuracy of the numerical solution (U ∗ , V ∗ ), we calculate the order of maximum error which is defined by

 orderα (h) = log2

errorα (h)

errorα (h/2)   erroru (h) = maxu ∗i − u ∗ (xi ), i

 ,

α = u, v , 



error v (h) = max v ∗i − v ∗ (xi ), i

(5.13)

where (u ∗i , v ∗i ) is the i-th component of (U ∗ , V ∗ ) and (u ∗ (xi ), v ∗ (xi )) is the continuous solution of (2.2), (5.3) at xi . The termination criterion of iterations is given by

 (m)  U − U (m) 

∞

  +  V (m) − V (m) ∞ < 10−13 .

(5.14)

In Table 2, we list the maximum errors erroru (h) and error v (h) as well as the orders of them for various values of h. The data in this table demonstrate that the numerical solution (U ∗ , V ∗ ) has the fourth-order accuracy. This coincides with the analysis very well.

Y.-M. Wang / Applied Numerical Mathematics 61 (2011) 38–52

49

Table 2 The accuracy of the numerical solution (U ∗ , V ∗ ) of Example 5.1. h

erroru (h)

orderu (h)

error v (h)

order v (h)

1/4 1/8 1/16 1/32 1/64 1/128 1/256 1/512

2.81395566738641e–002 1.24652816805571e–003 9.06433894494185e–005 5.51294410744418e–006 3.42257746477337e–007 2.13555905181906e–008 1.33351585329677e–009 8.10035372111884e–011

4.49661221054842 3.78156986465225 4.03930692923534 4.00966787470430 4.00239749843280 4.00130691012906 4.04110625795553

2.04635438672907e–002 1.18007005963311e–003 7.73722447332537e–005 4.70594950063852e–006 2.97095818502235e–007 1.85340933711586e–008 1.15732573524596e–009 7.04389879757628e–011

4.11611159332585 3.93091257322739 4.03925838585124 3.98548563403334 4.00267485004741 4.00131466422679 4.03827698371076

Example 5.2. Our second example is a logistic equation in population dynamics problems where the diffusion is density dependent. This equation in the one-dimensional domain [0, 1] is governed by (1.1) with

f (x, u ) = σ u (1 − u ) + q(x),

(5.15)

where σ is a positive constant and q(x)  0 is a possible internal source (see [7,11]). For this example, the function γ (x) in (2.3) can be chosen as γ (x) ≡ γ = σ (2ρ − 1)/k0 . In this case, the function g (x, u ) and the difference scheme (2.13) are reduced to (5.4) and (5.6), respectively, where f (x, u ) is replaced by (5.15). On the other hand, the conditions (5.1) and (5.2) for this example are, respectively, reduced to

σ α (1 − α ) + q(0)  0, and

σ β(1 − β) + q(1)  0,

(5.16)

⎧ σρ (1 − ρ ) + q(x)  0 (0  x  1), ⎪ ⎪ ⎨

α  ρ β 2 2 ⎪ k(s) ds  max k(s) ds + h f (0, α )/12, k(s) ds + h f (1, β)/12 . ⎪ ⎩ 0

0

(5.17)

0

We now consider the above problem with

k(u ) = arctan(1 + u ),

α = β = 0, σ = 2 and

q(x) = 2 arctan(1 + z) −

(2x − 1)2 − 2z(1 − z), 2 + 2z + z2

(5.18)

where z = x(1 − x). In this case, u = z is a solution of the problem, and the transformation T in (2.1) is given by

u v = T (u ) =

arctan(1 + s) ds = (1 + u ) arctan(1 + u ) −

1 2





ln 2 + 2u + u 2 −

π 4

√ + ln 2.

(5.19)

0

Clearly, it is difficult to find the inverse T −1 explicitly. With respect to the above parameters, the condition (5.16) is obviously satisfied, and the condition (5.17) holds by ρ = 2. This implies that the pair (U˜ , V˜ ) = (2E , T (2)E ) and (Uˆ , Vˆ ) = (0, 0) are a pair of ordered upper and lower solutions. Since ρ = 2 and arctan(1 + u )  π /4 for all u  0, we have k0 = π /4 and γ = 24/π . Using (U (0) , V (0) ) = (2E , T (2)E ) and (U (0) , V (0) ) = (0, 0), we compute the corresponding sequences {(U (m) , V (m) )} and {(U (m) , V (m) )} from (3.3) with Γ = γ I for the present problem, where γ = 24/π . Let h = 1/40. Some numerical results of these sequences at xi = 0.5 are plotted in Fig. 2, where the solid line denotes the sequence {U (m) } or { V (m) } and the dashdot line stands for the sequence {U (m) } or { V (m) }. As expected from our theoretical analysis in Theorem 3.2, the monotone convergence of the sequences holds and their common limit (U ∗ , V ∗ ) is the unique solution of the present problem in S given by (5.9) where ρ = 2. The maximum errors erroru (h) and error v (h) of the numerical solution (U ∗ , V ∗ ) as well as the orders of them for various values of h are presented in Table 3. In our computations, the termination criterion of iterations is still determined by (5.14). We see from Table 3 that the numerical solution (U ∗ , V ∗ ) has the fourth-order accuracy. Example 5.3. Our last example is given by (1.1) with





f (x, u ) = σ (1 − u ) exp −ν /(1 + u ) + q(x),

(5.20)

where σ and ν are positive constants, and q(x) is a nonnegative function. This problem arises from a chemical reactor with first-order reaction (cf. [10,12,20,31,32,34]). To find γ (x) in (2.3) we observe from (5.20) that





f u (x, u ) = −σ exp −ν /(1 + u )

 1−

ν (1 − u ) (1 + u )2

 .

(5.21)

50

Y.-M. Wang / Applied Numerical Mathematics 61 (2011) 38–52

Fig. 2. The monotone convergence of the sequences at xi = 0.5 for Example 5.2. Table 3 The accuracy of the numerical solution (U ∗ , V ∗ ) of Example 5.2. h

erroru (h)

orderu (h)

error v (h)

order v (h)

1/4 1/8 1/16 1/32 1/64 1/128 1/256 1/512

1.58684599293107e–004 1.01956076546672e–005 6.41229860487425e–007 4.01376447634050e–008 2.50951628921747e–009 1.56848256605002e–010 9.50076128880539e–012 5.40345546085064e–013

3.96014245092926 3.99096234202004 3.99781374138991 3.99947474754768 3.99996791204800 4.04518256310889 4.13608892220934

1.42185276064216e–004 9.13580885505971e–006 5.74577388973152e–007 3.59655524528879e–008 2.24866553155678e–009 1.40544686999533e–010 8.51330117512816e–012 4.84223772190262e–013

3.96009579872077 3.99095934055636 3.99781355033980 3.99947477046760 3.99996826760600 4.04516643685877 4.13597285658087

This implies that − f u (x, u )  σ (1 + ν ) for all u  0, and therefore, γ (x) can be chosen as γ (x) ≡ γ = σ (1 + ν )/k0 . Then the function g (x, u ) and the difference scheme (2.13) for this example are reduced to (5.4) and (5.6), respectively, where f (x, u ) is replaced by (5.20). To accommodate the exact solution of u = 1 − x2 /2, we let

k (u ) =



1 + u2 ,

q(x) =



1 + z2 − √

x2 z 1 + z2

 − σ (1 − z) exp −ν /(1 + z) ,

The above choice implies that α = 1, β = 1/2, k0 = 1 and q(x)  0 for all x ∈ [0, 1] if transformation T in (2.1) is given by

v = T (u ) =

u 

1 + s2 ds =

   1  u 1 + u 2 + lnu + 1 + u 2  . 2

z=1−

x2 2

(5.22)

. √

σ  3 exp(ν /2)/ 5. Moreover, the

(5.23)

0

As in the previous examples, we can prove that for any



√

√

√ 

σ + 2 exp(ν /2) 13 2 21 5 , ρ  max , , σ 12 40 ˆ ˆ the pair (U˜ , V˜ ) = (ρ E , T (ρ )E ) and √ (U , V ) = (0, 0) are a pair of ordered upper and lower solutions of the present problem. Let σ = 1, ν = 0.5 and ρ = 1 + 2 exp(1/4). We compute the corresponding numerical solution (U ∗ , V ∗ ) by the linear iterative scheme (3.3) with the initial iterations (U (0) , V (0) ) = (ρ E , T (ρ )E ) and (U (0) , V (0) ) = (0, 0). The termination criterion is the same as in Example 5.1. Fig. 3 shows the monotone convergence of the sequences {(U (m) , V (m) )} and {(U (m) , V (m) )} at xi = 0.5, where the solid line denotes the sequence {U (m) } or { V (m) } and the dash-dot line stands for the sequence {U (m) } or { V (m) } as before. The data in Table 4 show the maximum errors erroru (h) and error v (h) of the numerical solution (U ∗ , V ∗ ) as well as the orders of them for various values of h. We see that the fourth-order accuracy of the numerical solution (U ∗ , V ∗ ) is demonstrated. 6. Concluding remarks In this paper, we gave a numerical treatment for a class of strongly nonlinear two-point boundary value problems by Numerov’s method. The main contribution is to develop precisely a technique for computing accurately the numerical

Y.-M. Wang / Applied Numerical Mathematics 61 (2011) 38–52

51

Fig. 3. The monotone convergence of the sequences at xi = 0.5 for Example 5.3. Table 4 The accuracy of the numerical solution (U ∗ , V ∗ ) of Example 5.3. h

erroru (h)

orderu (h)

error v (h)

order v (h)

1/4 1/8 1/16 1/32 1/64 1/128 1/256 1/512

2.54113143148826e–005 1.68847875603451e–006 1.05727220867102e–007 6.62888643976345e–009 4.14343670485096e–010 2.58856269752528e–011 1.61481938931729e–012 1.14130926931466e–013

3.91167506288797 3.99730525735689 3.99543651708451 3.99986674925372 4.00060471650773 4.00270654236108 3.82261112679379

3.23999380957840e–005 2.16726185686689e–006 1.35707057924428e–007 8.48539494224809e–009 5.30407606724737e–010 3.31368266159871e–011 2.06989980711114e–012 1.42330591756945e–013

3.90204567565112 3.99730581562862 3.99937013513628 3.99980843898597 4.00059411265700 4.00080260662483 3.86224325092020

solutions of the problems by Numerov’s method without finding any inverse function. We generalized the method of upper and lower solutions to a coupled system of a nonlinear algebraic equation and a nonlinear functional equation. A linear monotone iterative algorithm was presented for the solutions of the resulting discrete problems. The monotone property of the iterations gives improved upper and lower bounds of the solution in each iteration. The proposed method complements the known method in [44] where it is necessary to find explicitly an inverse function, and so the full potential of Numerov’s method for strongly nonlinear two-point boundary value problems is realized. The numerical results coincide with the analysis very well. We conclude by taking note that the present study serves as an initial exploration of the application of higher-order compact finite difference methods to solve strongly nonlinear differential equations. The analysis is limited to a ordinary differential equation. The extension to partial differential equations is an interesting issue, and more challenging problems concerning this extension will be studied in the future. Acknowledgements The author would like to thank the referees for their valuable comments and suggestions which improved the presentation of the paper. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]

E. Adams, W.F. Ames, On contracting interval iteration for nonlinear problems in Rn : Part II: Applications, Nonlinear Anal. 5 (1981) 525–542. R.P. Agarwal, On Numerov’s method for solving two point boundary value problems, Util. Math. 28 (1985) 159–174. R.P. Agarwal, Difference Equations and Inequalities, Marcel Dekker, New York, 1992. R.P. Agarwal, Y.-M. Wang, The recent developments of the Numerov’s method, Comput. Math. Appl. 42 (2001) 561–592. J.H. Ahlberg, E.N. Nilson, J.L. Walsh, The Theory of Splines and Their Applications, Academic Press, New York, 1967. W.F. Ames, Nonlinear Partial Differential Equations in Engineering, Academic Press, New York, 1972. W.F. Ames, Numerical Methods for Partial Differential Equations, Academic Press, New York, 1977. Z.A. Anastassi, T.E. Simos, A family of two-stage two-step methods for the numerical integration of the Schrödinger equation and related IVPs with oscillating solution, J. Math. Chem. 45 (2009) 1102–1129. Z.A. Anastassi, T.E. Simos, Numerical multistep methods for the efficient solution of quantum mechanics and related problems, Phys. Rep. 482–483 (2009) 1–240. R. Aris, The Mathematical Theory of Diffusion and Reaction in Permeable Catalysts, Clarendon Press, Oxford, 1975. D.G. Aronson, H.F. Weinberger, Nonlinear diffusion in population genetics, combustion, and nerve propagation, in: Partial Differential Equations and Related Topics, in: Lecture Notes in Math., vol. 446, Springer, New York, 1975, pp. 5–49. J.W. Bebernes, D. Eberly, Mathematical Problems from Combustion Theory, Springer-Verlag, New York, 1989.

52

Y.-M. Wang / Applied Numerical Mathematics 61 (2011) 38–52

[13] L.K. Bieniasz, Use of the Numerov method to improve the accuracy of the spatial discretization in finite-difference electrochemical kinetic simulations, Comput. Chem. 26 (2002) 633–644. [14] R.L. Burden, J.D. Faires, Numerical Analysis, Prindle, Weber and Schmidt, Boston, 1993. [15] M.M. Chawla, P.N. Shivakumar, Numerov’s method for non-linear two-point boundary value problems, Int. J. Comput. Math. 17 (1985) 167–176. [16] M.M. Chawla, P.N. Shivakumar, A new fourth order cubic spline method for non-linear two-point boundary value problems, Int. J. Comput. Math. 22 (1987) 321–341. [17] J.P. Coleman, L.Gr. Ixaru, P -stability and exponential-fitting methods for y  = f (x, y ), IMA J. Numer. Anal. 16 (1996) 179–199. [18] L. Collatz, The Numerical Treatment of Differential Equations, third ed., Springer, Berlin, 1960. [19] Y. Fang, X. Wu, A trigonometrically fitted explicit Numerov-type method for second-order initial value problems with oscillating solutions, Appl. Numer. Math. 58 (2008) 341–351. [20] D.A. Frank-Kamenetskii, Diffusion and Heat Transfer in Chemical Kinetics, Plenum Press, New York, 1969. [21] B.-Y. Guo, J.J.H. Miller, Iterative and Petrov–Galerkin methods for solving a system of one-dimensional nonlinear elliptic equations, Math. Comp. 58 (1992) 531–547. [22] R.C. Gupta, Ravi. P. Agarwal, A new shooting method for multi-point discrete boundary value problems, J. Math. Anal. Appl. 112 (1985) 210–220. [23] P. Henrici, Discrete Variable Methods in Ordinary Differential Equations, John Wiley, New York, 1962. [24] S.R.K. Iyengar, P. Jain, Spline finite difference methods for singular two point boundary value problems, Numer. Math. 50 (1987) 363–376. [25] H.B. Keller, Elliptic boundary value problems suggested by nonlinear diffusion processes, Arch. Ration. Mech. Anal. 35 (1969) 363–381. [26] J.B. Keller, W.E. Olmstead, Temperature of a nonlinear radiating semi-infinite solid, Quart. Appl. Math. 29 (1972) 559–566. [27] A. Konguetsof, T.E. Simos, An exponentially fitted and trigonometrically fitted method for the numerical solution of periodic initial value problems, Comput. Math. Appl. 45 (2003) 547–554. [28] B.V. Numerov, A method of extrapolation of perturbations, Monthly Notices Royal Astron. Soc. 84 (1924) 592–601. [29] W.E. Olmstead, Temperature distribution in a convex solid with nonlinear radiation boundary condition, J. Math. Mech. 15 (1966) 899–908. [30] M.N. Özisik, Boundary Value Problems of Heat Conduction, Dover, New York, 1989. [31] C.V. Pao, Monotone iterative methods for finite difference system of reaction diffusion equations, Numer. Math. 46 (1985) 571–586. [32] C.V. Pao, Nonlinear Parabolic and Elliptic Equations, Plenum Press, New York, 1992. [33] C.V. Pao, Finite difference reaction diffusion equations with nonlinear boundary conditions, Numer. Methods Partial Differential Equations 11 (1995) 355–374. [34] S.V. Parter, Solutions of a differential arising in chemical reactor process, SIAM J. Appl. Math. 26 (1974) 687–716. [35] G. Psihoyios, T.E. Simos, A fourth algebraic order trigonometrically fitted predictor-corrector scheme for IVPs with oscillating solutions, J. Comput. Appl. Math. 175 (2005) 137–147. [36] S.M. Roberts, J.S. Shipman, Two Point Boundary Value Problems: Shooting Methods, American Elsevier, New York, 1972. [37] T.E. Simos, A Numerov-type method for the numerical-solution of the radial Schrödinger-equation, Appl. Numer. Math. 7 (1991) 201–206. [38] T.E. Simos, A high-order predictor-corrector method for periodic IVPS, Appl. Math. Lett. 6 (1993) 9–12. [39] T.E. Simos, Exponentially-fitted and trigonometrically-fitted symmetric linear multistep methods for the numerical integration of orbital problems, Phys. Lett. A 315 (2003) 437–446. [40] S. Stavroyiannis, T.E. Simos, Optimization as a function of the phase-lag order of nonlinear explicit two-step P -stable method for linear periodic IVPs, Appl. Numer. Math. 59 (2009) 2467–2474. [41] E.H. Twizell, S.I.A. Tirmizi, Multiderivative methods for non-linear second-order boundary value problems, J. Comput. Appl. Math. 17 (1987) 299–307. [42] H. Van de Vyver, Phase-fitted and amplification-fitted two-step hybrid methods for y  = f (x, y ), J. Comput. Appl. Math. 209 (2007) 33–53. [43] Y.-M. Wang, Monotone methods for a boundary value problem of second-order discrete equation, Comput. Math. Appl. 36 (1998) 77–92. [44] Y.-M. Wang, Numerov’s method for strongly nonlinear two-point boundary value problems, Comput. Math. Appl. 45 (2003) 759–763. [45] Y.-M. Wang, The extrapolation of Numerov’s scheme for nonlinear two-point boundary value problems, Appl. Numer. Math. 57 (2007) 253–269. [46] Y.-M. Wang, B.-Y. Guo, On Numerov scheme for nonlinear two-points boundary value problem, J. Comput. Math. 16 (1998) 345–366.