NUMERICAL METHODS FOR COUPLED SYSTEMS OF QUASILINEAR ELLIPTIC EQUATIONS WITH NONLINEAR BOUNDARY CONDITIONS

This paper is concerned with numerical solutions of a coupled system of arbitrary number of quasilinear elliptic equations under combined Dirichlet and nonlinear boundary conditions. A finite difference system for a transformed system of the quasilinear equations is formulated, and three monotone iterative schemes for the computation of numerical solutions are given using the method of upper and lower solutions. It is shown that each of the three monotone iterations converges to a minimal solution or a maximal solution depending on whether the initial iteration is a lower solution or an upper solution. A comparison result among the three iterative schemes is given. Also shown is the convergence of the minimal and maximal discrete solutions to the corresponding minimal and maximal solutions of the continuous system as the mesh size tends to zero. These results are applied to a heat transfer problem with temperature dependent thermal conductivity and a Lotka-Volterra cooperation system with degenerate diffusion. This degenerate property leads to some interesting distinct property of the system when compared with the non-degenerate semilinear systems. Numerical results are given to the above problems, and in each problem an explicit continuous solution is constructed and is used to compare with the computed solution.


Introduction
In the theory of heat transfer if the thermal conductivity is temperature dependent and the boundary is subjected to a Boltzmann fourth-power radiation law then for certain participating media such as glass, fibers, and powders, the steady-state temperature u(x) is governed by the quasilinear boundary-value problem D(u) and f (x, u) are given by for some α ≥ 1, where k c and k r are the thermal conductivity constants based on conduction and radiation, respectively, while p(x) and c(x) are positive functions representing the internal source and evaporation coefficients(cf [15, p249]). The function D(u) in (1.1a) is due to the effect of simultaneous conduction and radiation in the participating medium.
On the other hand, in the Lotka-Volterra cooperation system of two-cooperating species the steady-state density functions of the species u(x), v(x) are governed by the coupled system − ∆u m = u(a (1) where m > 1, n > 1 and (a (l) , b (l) , c (l) ), l = 1, 2, are positive constants(cf [17,19]). The terms ∆u m , ∆v n with m > 1, n > 1 implies that the diffusion rate of movement from high density region to low density region is slow, and it models a tendency to avoid crowding(cf. [12,13,23]). Problem (1.2) and its corresponding time-dependent system have been extensively investigated in the field of ecology, but most of the investigations are devoted to the semilinear case m = n = 1 (cf [6-8, 13, 27, 28]). This system has also been used as a mathematical model in economics where u and v represent the wealth of two cooperating nations or regions (cf. [1]).
In particular, the boundary condition in (1.3) is of Dirichlet type for all l if n 0 = N + 1, and it is nonlinear (or linear Neumann-Robin type) if n 0 = 1. We also allow D (l) (u (l) ) = d (l) to be independent of u (l) for some or all l. In the later case, problem (1.3) becomes a coupled system of semilinear elliptic boundary problems.
The purpose of this paper is the following: (a) To present three monotone iterative schemes for the computation of positive minimal and maximal solutions of a discrete system of (1.3), including the existence of these solutions. (b) To show the convergence of the discrete solutions to the corresponding solutions of the continuous problem as the mesh size decreases to zero. (c) To apply the results in (a) and (b) to the heat transfer problem (1.1) and the cooperation system (1.2), including some numerical results for these two problems. In the computation of numerical solutions we first construct a source function so that a continuous solution of the problem is explicitly known and is used to compare with the computed solution by the monotone iterative schemes. Our numerical results demonstrate excellent agreement between the computed solution and the true continuous solution.
Literature dealing with numerical solutions of quasilinear elliptic equations is extensive and various topics of the problem, such as method of computation, error estimate, and convergence of the discrete solution have been discussed (cf. [2,9,11,21,26]). Most of the discussions in the above works are for scalar equations with linear boundary condition. The work in [21] treated a scalar quasilinear equation with either Dirichlet or nonlinear boundary condition. Similar discussions for semilinear elliptic equations, including coupled system of two or more equation have been treated by either the finite difference method or the finite element method (cf. [10,16,18,19]). For numerical solutions of quasilinear elliptic equations, it is often assumed that the equation is non-degenerate and the problem has a unique solution. An example of this situation is the heat transfer problem (1.1) where the thermal conductivity is temperature dependent. This equation with D(u) = k c + k r u 3 and the physical meaning of the constants k c and k r have been derived in [3,15]. Other models and many relevant references in this area can be found in chapter 9 of [15]. On the other hand, in many reaction-diffusion problems with density dependent diffusion coefficients the governing equations are either degenerate or it is nondegenerate but the solution is not unique. This type of equations, including the system (1.3) have been treated recently in [22]. The emphasis in this paper is on the degenerate property of the system. This property and the non-uniqueness of the solution is especially true for coupled system of equations, including semilinear systems. For example, in the semilinear system (1.2) where m = n = 1, the problem has only the trivial solution (0, 0) if a (l) ≤ λ 0 , l = 1, 2, where λ 0 > 0 is the smallest eigenvalue of the eigenvalue problem ∆φ + λφ = 0 in Ω, φ = 0 on ∂Ω.
(1. 4) In the case a (l) > λ 0 the problem has two semitrivial solutions in the form (u s , 0), (0, v s ), where u s > 0 and v s > 0 in Ω. In addition, it has at least one positive solution if b (1) b (2) > c (1) c (2) , and it has no positive solution if b (1) b (2) < c (1) c (2) (cf. [17, p678]). However, the situation is quite different for the quasilinear system (1.2) and its corresponding finite difference system when m > 1, n > 1. In fact, we show that for any positive constants (a (l) , b (l) , c (l) ), l = 1, 2, the corresponding finite difference system of (1.2) has two semitrivial solutions (u(x i ), 0), (0, v(x i )), where x i denotes a mesh point inΩ. If, in addition, 1/m + 1/n < 1 then the problem has a positive minimal solution (u(x i ), v(x i )) and a positive maximal solution (u(x i ), v(x i )). If, in addition, (u(x i ), v(x i )) = (u(x i ), v(x i )) then their common value is the unique positive finite difference solution (see Theorem 5.3). The above conclusions hold true for every positive constants (a (l) , b (l) , c (l) ), including the case a (l) ≤ λ 0 or b (1) b (2) ≤ c (1) c (2) . These results are in sharp contrast to the semilinear case m = n = 1.
To compute numerical solutions of the above problems, including positive mini-mal and maximal solutions we use the method of upper and lower solutions and its associated monotone iterations for a finite difference system of the coupled system (1.3). This approach leads to some general results for the system which can be used not only for the problems in (1.1) and (1.2) but also for many other degenerate and non-generate reaction-convection-diffusion system with density-dependent diffusions. An important feature of this method is that it can be used to compute positive minimal and maximal solutions without the uniqueness assumption. In particular, if the solution is unique then the monotone convergence of the minimal and maximal sequences gives a simple and reliable error estimate of the solution without any explicit knowledge of the solution. The plan of the paper is as follows: In section 2, we transform the system (1.3) into a system of semilinear and algebraic equations which are discretised into a finite difference system. A monotone iterative scheme is given for computing positive solutions using the method of upper and lower solutions. Proofs of the monotone convergence of the sequence of iterations together with two other monotone iterative schemes are given in Section 3. In Section 4, we show the convergence of the finite difference solution to the continuous solution as the mesh size tends to zero. The above results are applied to problem (1.1) and (1.2) in Section 5 where explicit lower and upper solutions are constructed. Finally, in Section 6 we give some numerical results for these two problems, including the construction of known continuous solutions which are used to compare with the computed solutions.

Finite Difference Systems
To develop computational schemes for numerical solutions of (1.3) we make a transformation by letting (2.1) In view of dw (l) /du (l) = D (l) (u (l) ) > 0 for u (l) > 0 the inverse function of w (l) = I (l) [u (l) ], denoted by u (l) = q (l) (w (l) ), exists and dq (l) /dw (l) = 1/D (l) (u (l) ). Since for each l = 1, . . . , N , where η (l) = I (l) [ξ (l) ]. Although the above system may be written in a semilinear form in terms of w (l) by using u = q(w) = (q (1) (w (1) ), . . . , q (N ) (w (N ) )) we prefer to use its present form as a coupled system of w and u where w = (w (1) , . . . , w (N ) ). One reason for this is that in the degenerate case D (l) (0) = 0 the functions f (l) (x, q(w)), g (l) (x, q(w)) are not Lipschitz continuous in w (l) at w (l) = 0 because dw (l) /du (l) = 1/D (l) (u (l) ) which is not defined at u (l) = 0. It is obvious that problem (2.2) can be discretised into a discrete system in the same fashion as for semilinear elliptic boundary-value problems. In this paper, we use the finite difference method for the discrete system although it can also be formulated by the finite element method with a suitable choice of the basis functions (see Remark 3.1). Let i = (i 1 , . . . , i p ) be a multiple index with i ν = 1, . . . , M ν , and let x i = (x i1 , . . . , x ip ) be a mesh point inΩ ≡ Ω∪∂Ω, where ν = 1, . . . , p and M ν is the total number of intervals in the x ν -direction. Denote by Λ, ∂Λ andΛ the sets of mesh points in Ω, ∂Ω andΩ, respectively. When no confusion arises we write i ∈ Λ for x i ∈ Λ where Λ stands for Λ, ∂Λ orΛ. Let h ν be the spatial increment in the x νdirection, and let u where l = 1, . . . , N . Then by the central difference approximations and the boundary approximation where e ν is the unit vector in lR p with the ν-th component one and zero elsewhere andx i is a neighboring point of x (b) in Λ, we approximate the transformed system (2.2) by the finite difference system [5,24]).
To develop monotone iterative schemes for the computation of solutions of (2.5) we use the method of lower and upper solutions. The definition of these solutions, denoted by (û i ,ŵ i ) = ((û )) respectively, are defined in the following.
Similarly, (ũ i ,w i ) is called an upper solution of (2.5) if its components satisfy all the inequalities in (2.6) in reversed order. The For a given pair of ordered lower and upper solutions (û i ,ŵ i ), (ũ i ,w i ) we set In the following discussion we assume that a pair of ordered lower and upper solutions (û i ,ŵ i ), (ũ i ,w i ) exist and (û i ,ŵ i ) ≥ (0, 0). In addition, we make the following basic hypothesis for each l = 1, . . . , N : and f (l) (·, u), g (l) (·, u) are C 1 -functions of u for u ∈ S (1) . (iii) f (l) (·, u) and g (l) (·, u) are quasi-monotone nondecreasing in u for u ∈ S (1) and there exist nonnegative functions γ (l) (x) and γ (l) (x), not both identical zero, such that In the hypothesis (H 1 ) − (i) we allow D (l) (0) = 0 for some or all l. In this situation the system (2.5) is degenerate on the boundary if ξ (l) (x) = 0 for some or all It is easy to see from the hypothesis D (l) (u (l) ) > 0 for u (l) > 0 that condition (2.8) holds by any nonnegative functions γ (l) (x), γ (l) (x) if ∂f (l) (·, u)/∂u (l) ≥ 0 and ∂g (l) (·, u)/∂u (l) ≥ 0. It is also satisfied by some positive functions γ (l) (x), γ (l) (x) if D (l) (0) > 0. Hence this condition is needed only for the degenerate case D (l) (0) = 0. Define where w (l) i ] is given by (2.1). Then we may write problem (2.5) in the equivalent form It is obvious from (2.8), (2.9) and (d/du (l) )I[u By the quasi-monotone nondecreasing property of f (l) (·, u), g (l) (·, u) in the hypothesis (H 1 ) − (iii) we see that for every l = 1, . . . , N , i ) (k) ), l = 1, . . . , N , are the components of (u i } is well-defined and can be easily computed. Specifically, starting from any u i ) (1) of (2.12) for k = 1 by solving a linear finite difference system under Dirichlet boundary condition for l = 1, . . . , n 0 − 1 and under Robin boundary condition for l = n 0 , . . . , N . This is because for k = 1 the functions at the right-hand side of the first three equations are known. Knowing the value of (w (2.12). This gives the first iteration ((u (1) ) in (2.12) instead of u (0) the same computational procedure for k = 2 yields the second iteration ((u (l) i ) (2) , (w (l) i ) (2) ). A continuation of this process leads to the k-th iteration ((u i ) (k) ) for every k. It is clear that this sequence of iterations depends on the initial iteration u   i }, respectively. The following theorem gives the monotone convergence of these sequences. Theorem 2.1. Let (û i ,ŵ i ), (ũ i ,w i ) be a pair of ordered lower and upper solutions of (2.5), and let hypothesis (H 1 ) be satisfied. Then the following statements hold true: (a) The sequence {u i } converges to a maximal solution (u i , w i ), and they satisfy the relation (2) , if any, satisfies the relation is the unique solution of (2.5) in S (2) .
Although the above theorem can be proved directly from (2.12) we postpone its proof to the next section (along with two other monotone iterative schemes) using vector form of the finite difference system.

Monotone Sequences
To show the monotone convergence of the sequence given by (2.12) and the sequences governed by two other iterative schemes, called Gauss-Seidel and Jacobi iterations, we formulate the system (2.5) in vector form. Let and matrices  Since our main concern in the proof of monotone convergence of the sequences is the mathematical structure of the discrete system (2.10), detailed formulation of (3.3) is omitted (see [18] for some detailed discussion). However, it is to be noted from   It can be easily shown that if (û i ,ŵ i ) is a lower solution of (2.5), whereû i = (û i , · · · ,ŵ (N ) i ), then its vector representation (Û,Ŵ) witĥ U = (Û (1) , · · · ,Û (N ) ),Ŵ = (Ŵ (1) , · · · ,Ŵ (N ) ) is a lower solution of (3.3). The same is true for an upper solution. Since many of the techniques for the construction of lower and upper solutions of the differential system (2.2) can be used for the finite difference system (2.5) it is often convenient to use Definition 2.1 for the search of explicit lower and upper solutions of specific problems.
To ensure the existence of a positive solution of (3.3) and to obtain a computational algorithm for its numerical values we assume that a pair of ordered lower and upper solutions (Û,Ŵ), (Ũ,W) exist and impose the following hypothesis for each l = 1, . . . , N .
It is well-known that under the Hypothesis (H 2 ) − (i) the inverse matrix (A (l) ) −1 exists and is a positive matrix if strict inequality in (3.5) holds for at least one j (cf. [25,29]). This implies that for any nonnegative matrix Γ (l) , the inverse (A (l) ) −1 exists and is positive for is sometimes called a monotone matrix (cf. [14]). This type of matrices has been widely used in matrix theory and in linear and semilinear elliptic equations(cf. [14,16,18,19,21]). It is easy to verify from the approximation (2.3)-(2.4) and the assumption h ν < |b ij are satisfied. The same approximation and the connectedness assumption on Ω implies that A (l) is irreducible. Moreover, condition (2.11) ensures that the condition in (H 2 ) − (ii) on F(U) and G(U) are satisfied. Hence, under the hypothesis (H 1 ) all the requirements in (H 2 ) are fulfilled. We first consider some computational algorithms for numerical solutions.

Three monotone iterative schemes.
In this subsection we present three monotone iterations called Picard, Gauss-Seidel, and Jacobi iterations, using eitherÛ orŨ as the initial iteration.
(A) Picard iteration: For the Picard iteration the sequence is governed by where k = 1, 2, . . .. It is easily seen that the sequence given by (3.6) is simply a vector representation of the sequence governed by (2.12). Since starting from k = 1 the function U (0) is eitherÛ orŨ the right-hand side of the first two equations in (3.6) is known. This implies that the first iteration (W (l) ) (1) exists and can be computed by solving a linear algebraic equation for each l = 1, . . . , N . Knowing the values of (W (l) ) (1) we can compute the value of (U (l) ) (1) from the last equation in (3.6). Using the value of U (1) ≡ ((U (1) ) (1) , · · · , (U (N ) ) (1) ) instead of U (0) we can compute the solution (W (l) ) (2) and then (U (l) ) (2) for each l. Repeating this process leads to the sequence {U (k) , W (k) } for every k = 1, 2, . . .. Denote the sequence by , and refer to them as minimal and maximal sequence, respectively. The following theorem gives the monotone convergence of these sequences.
The positivity of (A (l) ) −1 and ( Moreover, by (3.6) and the condition in (H 2 ) − (ii), we have This leads to (W (l) ) (1) ≥ (W (l) ) (1) for every l = 1, . . . , N . This result and ( Assume, by induction, that for some k > 1. Then by (3.6) and (H 2 ) − (ii), The positivity of (A (l) ) −1 and (A (l) ) −1 implies that (W (l) ) (k+1) ≥ (W (l) ) (k) for every ). The monotone property of the minimal and maximal property in (3.7) follows from the principle of induction. This monotone property implies that for every l = 1, . . . , N , the limits exist and the vectors (U, W), satisfy the relation (3.7). Letting k → ∞ in (3.6) shows that both (U, W) and (U, W) are solutions of (3.3). (The minimal and maximal property of these solutions is a consequence of the result in (b) below). This proves the conclusion in (a).
(b). To show the relation (3.8) for any solution (U, W) ∈ S we observe from (H 2 ) − (ii) andÛ ≤ U ≤Ũ that the components (U (l) , W (l) ) of (U, W) satisfy the relation This proves the relation (3.8), and therefore the minimal and maximal property of (U, W) and (U, W).
The uniqueness of the solution (U * , W * ) in S follows from (3.8). This proves the theorem. To treat the system (3.3) by the Gauss-Seidel and Jacobi iterations we write Then Using the above notation we have the following additional iterative schemes: (B) Gauss-Seidel Iteration.
where (U (0) , W (0) ) is either (Û,Ŵ) or (Ũ,W), and k = 1, 2, . . .. Since by the Hypothesis (H 2 )−(i) the inverses (G (l) ) −1 , (J (l) ) −1 and (G (l) ) −1 , (J (l) ) −1 exist and are nonnegative, the sequences governed by (3.10) and (3.11) are well-defined and can be computed by solving a linear algebraic system with triangular and diagonal coefficient matrices, respectively. Denote the sequence again by , and refer to them as minimal and maximal sequences, respectively. The following theorem gives the monotone convergence of these sequences. Proof. Gauss-Seidel iteration We first show the monotone property of the minimal and maximal sequences governed by (3.10).
This gives (Z (l) ) (k+1) ≥ 0 which yields (W The first inequality in (3.13) for the minimal sequence follows from the principle of induction.
Using the relation (U This proves the second inequality of (3.13) for the minimal sequences.
Remark 3.1. (a). Theorem 3.2 implies that with the same initial iterationÛ or U the sequence governed by the Picard iteration converges faster than the Gauss-Seidel iteration which, in turn, converges faster than the Jacobi iteration. However, the Jacobi iteration is the simplest to use in practical computation while the Picard iteration may require additional iterations for each k if the spatial domain Ω is of two or higher dimension. On the other hand, since G (l) and G (l) are triangular matrices it is more suitable for practical computation when Ω is of higher dimension. (b). In Theorem 3.1 and Theorem 3.2 if U = U or W = W then the solution (U * , W * ) is unique in S and satisfies the relation (3.7). This implies that an error estimate for (U * , W * ) is given by and this error decreases to zero as k → ∞. This result and the monotone convergence of the minimal and maximal sequences to (U * , W * ) is very useful in practical computation since it does not require any explicit knowledge of the solution.

Convergence of finite difference solutions
Using the method of lower and upper solutions for both continuous and finite difference problems we can prove the convergence of the finite difference solution (u i , w i ) of (2.5) to the continuous solution (u(x i ), w(x i )) of (2.2) at every mesh point x i of a given partition Λ * . Recall that a smooth function (û(x),ŵ(x)) is called a lower solution of (2.2) if it satisfies (2.2) with all the equality sign "=" replaced by the inequality sign "≤". Similarly, (ũ(x),w(x)) is called an upper solution if it satisfies (2.2) with "=" replace by "≥" (cf [22]). This pair of lower and upper solutions are said to be ordered if (û(x),ŵ(x)) ≤ (ũ(x),w(x)) in Ω. Assume that problem (2.2) has a pair of ordered lower and upper solutions, and the functions c(x), η(x), f (x, u) and g(x, u) are smooth functions in their respective domains. Then by using either (û(x),ŵ(x)) or (ũ(x),w(x)) as the initial iteration (u (0) , w (0) ) we can construct a sequence {u (k) , w (k) } from the linear iteration process for k = 1, 2, . . ., where ((u (l) ) (k) , (w (l) ) (k) ), l = 1, . . . , N , are the components of (u (k) , w (k) ) and It is obvious that the sequence {u (k) , w (k) } governed by (4.1) is well-defined and can be obtained by solving a linear Dirichlet boundary problem for l = 1, . . . , n 0 − 1 and a Robin boundary problem for l = n 0 , . . . , N . We denote the sequence by , and refer to them as minimal and maximal sequence, respectively. In the following theorem we state the monotone convergence of these sequences from [22].
A proof of the above theorem can be found in a recent article [22]. By using the above monotone convergence theorem and the results in Theorem 2.1 we show the convergence of the minimal and maximal finite difference solutions to the corresponding continuous solutions at every mesh point x i ∈ Λ * , where Λ * is a fixed partition of Ω. It is assumed that every refinement of Λ * contains Λ * , and there exist a pair of ordered lower and upper solutions (û(x),ŵ(x)), (ũ(x),w(x)) of (2.2) and (û i ,ŵ i ), (ũ i ,w i ) of (2.5). It is also assumed that given any 0 > 0 there exists δ 0 > 0 such that for |h| < δ 0 and l = 1, . . . , N , where (û (l) (x i ),ŵ (l) (x i )) and (ũ (l) (x i ),w (l) (x i )) are the respective components of (û(x i ),ŵ(x i )) and (ũ(x i ),w(x i )) and |h| = h 1 + . . . + h p .
Since by the mean-value theorem, and similarly, and ξ, ξ andξ are some intermediate values between (u (j) (x i )) (k−1) and (u (j) i ) (k−1) and between (w (j) (x i )) (k) and (w (j) i ) (k) , respectively, we see that the above system is equivalent to Notice from the nondecreasing property of F Let (V (l) ) (k) , (Z (l) ) (k) be the vector representation of (v (l) i ) (k) and (z (l) i ) (k) respectively, and let Then in vector form we may write (4.7) as Then the above inequalities become It is well-known that given any 1 > 0 there exists a matrix norm and a vector norm such that for every Y ∈ lR M (l) and Y ∈ lR M (l) , where µ (l) ≥ 0 and µ (l) ≥ 0 are the respective principle eigenvalues of A (l) and A (l) (cf. [25,29]). We choose 1 sufficiently small (and choose γ (l) > 1 if µ (l) = 0) so that ρ (l) > 0 for each l = 1, . . . , N . A similar choice gives ρ (l) > 0. Let ρ 0 ≥ max{ρ (l) , ρ (l) } for all l. Then by (4.9) and (4.10), where o (k) (|h|) → 0 as |h| → 0 and Addition of the first two inequalities in (4.11) and also the last inequality in (4.11) over l yield where o(|h|) = maximal {o (k) (|h|), k ≤ k * }. Let Then (4.12) is satisfied if where M is the total number of components of U (0) , we see from (4.13) that for any 0 > 0 there exists δ 0 such that r (0) < 0 and o(|h|) < when |h| < δ 0 . This implies that Using the relation (4.13), an induction argument gives when |h| < δ 0 . Let k ≥ k * be fixed. Then by choosing 0 sufficiently small there exists δ > 0 such that s (k) + r (k) < /3 when |h| < δ. This is equivalent to which ensures that the relation in (4.5) holds. This proves the convergence of the minimal solution. The proof for the maximal solution is similar and is omitted.
If the minimal solution (u(x), w(x)) coincide with the maximal solution (u(x), w(x)), then Theorem 4.1 implies that both (u i , w i ) and (u i , w i ) converge to their common value (u * (x i ), w * (x i )). This observation leads to the following.

Applications
In this section, we given some applications of the existence theorem and the monotone iterative schemes to the model problems ( i (u (k−1) ), l = n 0 , . . . , N, For the heat-transfer problem (1.1) we impose the following hypothesis: and k c , k r , α and a 0 are positive constants with α ≥ 1 It is clear that problem (1.1) is a special case of (1.3) with N = n 0 = 1, u (1) = u, b (1) = 0 and This implies that w i = I[u i ] = k c u i + (k r /4)u 4 i , and the finite difference system for problem (1.1) becomes where u i is determined from the equation k c u i + (k r /4)u 4 i = w i . By (5.1) the iteration process for problem (5.3) is given by By Theorem 2.1 (or Theorem 3.1) the finite difference problem (5.3) has a minimal solution (u i , w i ) and a maximal solution (u i , w i ) such that The positive property of p i ensures that (u i , w i ) > (0, 0) in Λ (cf. [21]). Moreover, by the relation ∂f (1)   i } converge monotonically to (u * i , w * i ) and satisfy the relation   (B). The Lotka-Volterra cooperation system. We next consider an extended finite difference system of the cooperation model (1.2) which is given by i ≥ 0 are some nonnegative functions in Λ. The consideration of the above source functions is for the purpose of constructing a known continuous solution. It is clear that the above problem is a special case of (1.3) with N = 2, n 0 = 3, b (1) = 0, (u (1) , u (2) ) = (u, v), and i ] = v n i , the transformed system of (5.10) becomes Recall from Definition 2.1 that lower and upper solutions of (5.12), denoted by (û i ,ŵ i ) ≡ ((û i ,v i ), (ŵ i ,ẑ i )) and (ũ i ,ṽ i ) ≡ ((ũ i ,ṽ i ), (w i ,z i )), are required to satisfy (5.12) with all the equality sign "=" replaced by the respective inequality sign "≤" and "≥". Let such a pair be given. Then from the relation we see that all the conditions in Hypothesis (H 1 ) are satisfied by any γ (5.14) Using the functions γ (1) , γ (2) determined from (5.14) we compute the minimal and maximal sequences {û (k) ,ŵ (k) }, {ũ (k) ,w (k) } from the iteration process Hence to show the existence and the computation of a positive solution of (5.12) we need to find a pair of positive lower and upper solutions (û i ,ŵ i ), (ũ i ,w i ). In fact, it suffices to . Before doing this we show the existence of two semitrivial solutions of the system (5.10) for the case q Lemma 5.1. Let m > 1, n > 1 and q Then for any positive constants a (l) , b (l) and c (l) , l = 1, 2, problem (5.10) has two semitrivial solutions To show the existence of the semitrivial solution (u * i , 0) we apply Theorem 2.1 to the scalar problem It is obvious that for any constant ρ ≥ a (1) /b (1) , (ũ i ,w i ) = (ρ, ρ m ) is an upper solution. To find a positive lower solution we let λ * and φ i be the smallest eigenvalue and its corresponding (normalized) positive eigenfunction of the eigenvalue problem where λ * > 0. In terms of the matrix A (1) , λ * is the smallest eigenvalue of A (1) and Φ ≡ (φ 1 , · · · , φ M ) T is the positive eigenvector. It is clear that for any small In view of (5.17) the above inequality is equivalent to Since m > 1 and φ i > 0 there exists a small constant δ * > 0 such that the above inequality holds for every δ ≤ δ * . This shows that the pair (û i ,ŵ i ) = ((δφ i ) 1/m , δφ i ) and (ρ, ρ m ) are ordered lower and upper solutions of (5.16). By Theorem 2.1(or Theorem 3.1), this problem has a positive solution (u * i , w * i ) which ensures that (u * i , 0) is a semitrivial solution of (5.10) (for q Lemma 5.2. Given any m > 1, n > 1, there exists δ * > 0 such that for every δ ≤ δ * the function (û i ,ŵ i ), wherê is a positive lower solution of (5.12).

By (5.17) and q (l)
i ≥ 0 for l = 1, 2, the above inequalities are satisfied if In view of m > 1, n > 1, and 0 ≤ φ ≤ 1 there exists a small constant δ * > 0 such that the above relation holds for every δ ≤ δ * . This proves the lemma.
We next give some sufficient conditions for the construction of a positive upper solution (ũ i ,w i ) = ((ũ i ,ṽ i ), (w i ,z i )). This function is given in the form for a sufficiently large constant ρ, where ψ i is either the positive eigenfunction of (5.17) in a slightly large domainΛ containing Λ or ψ i = 1 in Λ. In fact, we choose , where ρ (1) and ρ (2) are given by and a (l) , l = 1, 2, are any constants satisfying a (l) ≥ a (l) + q (l) i . Specifically we have the following results.
i ≥ 0, problem (5.10) has a positive minimal solution (u i , v i ) and a positive maximal solution (u i , v i ) that satisfy the relation (5.27).  (A). The heat-transfer problem Consider problem (1.1) in a disk region of radius R = 1 for the case α = 1 and c 0 (x) = c 0 , where the solution depends only on the radial direction. The physical problem may be considered as a fuel rod of a fuel assembly in a nuclear reactor where the temperature in the rod due to fission depends on the radial direction only (cf. [3,17]). To demonstrate the accuracy and reliability of the monotone iterative schemes we construct a particular source function p(r) so that the solution u * (r) of the continuous problem is explicitly known. The values of this continuous solution will be used to compare with the computed finite difference solution u * i at every mesh point r i in a given partition Λ * . In the disk region the continuous problem (1.1) (for the case α = 1) becomes

Numerical results
where σ(0) = 0, σ 1 = σ(1) and c(r) ≡ c 0 > 0. By letting w = k c u + (k r /4)u 4 the transformed system of (6.1) is given by It is easy to see by the central difference approximation (2.3) that the finite difference system of the above boundary-value problem is given by In vector form, this system may be written as and q(w i ) = u i which is determined from the equation w i = k c u i + (k r /4)u 4 , for i = 1, . . . , N .
To construct a continuous solution of (6.1) for any physical parameters k c , k r , c 0 , σ 1 and a 0 we choose where a > 1 is a constant satisfying the relation 3) It is easy to verify that for any choice of p(r) and a > 1 in (6.2) and (6.3), the solution of (6.1) is given by In particular, if we choose a ≥ 4 then p(r) ≥ 0 and is bounded by 4(k c + k r a 3 ) + c 0 a for r ∈ [0, 1]. This implies that the constant P in (5.6) becomes P = a + (4/c 0 )(k c + k r a 3 ). (6.4) Using the above value of P and any ρ ≥ maximal {P , a 0 } in (5.8) for γ (1) and γ (2) we compute the minimal and maximal sequences {u i } from the iteration process (5.1) (or any one of the iteration process in (3.6), (3.10) and (3.11)). The initial iterations for these sequences are (u The physical constants in the above relation can be arbitrarily chosen. For example, by choosing  (B). The cooperation system.