**MONOTONE FINITE DIFFERENCE DOMAIN** **DECOMPOSITION ALGORITHMS AND**

**APPLICATIONS TO NONLINEAR SINGULARLY** **PERTURBED REACTION-DIFFUSION PROBLEMS**

IGOR BOGLAEV AND MATTHEW HARDY

*Received 16 September 2004; Revised 21 December 2004; Accepted 11 January 2005*

This paper deals with monotone finite diﬀerence iterative algorithms for solving non- linear singularly perturbed reaction-diﬀusion problems of elliptic and parabolic types.

Monotone domain decomposition algorithms based on a Schwarz alternating method and on box-domain decomposition are constructed. These monotone algorithms solve only linear discrete systems at each iterative step and converge monotonically to the ex- act solution of the nonlinear discrete problems. The rate of convergence of the mono- tone domain decomposition algorithms are estimated. Numerical experiments are pre- sented.

Copyright © 2006 Hindawi Publishing Corporation. All rights reserved.

**1. Introduction**

We are interested in monotone discrete Schwarz alternating algorithms for solving non- linear singularly perturbed reaction-diﬀusion problems.

The first problem considered corresponds to the singularly perturbed reaction-diﬀu- sion problem of elliptic type

*−**μ*^{2}^{}*u**xx*+*u**yy*

+*f*(x,*y,u)**=*0, (x,*y)**∈**ω,*
*u**=**g* on*∂ω,* *ω**=**ω*^{x}*×**ω*^{y}*= {*0*< x <*1*} × {*0*< y <*1*}*,

*f**u**≥**c**∗*, (x,*y,u)**∈**ω**×*(*−∞*,*∞*), *f**u**≡**∂ f /∂u,*

(1.1)

where*μ*is a small positive parameter,*c**∗**>*0 is a constant,*∂ω* is the boundary of*ω. If*
*f* and*g*are suﬃciently smooth, then under suitable continuity and compatibility condi-
tions on the data, a unique solution*u*of (1.1) exists (see [6] for details). Furthermore,
for*μ*1, problem (1.1) is singularly perturbed and characterized by boundary layers
(i.e., regions with rapid change of the solution) of width*O(μ**|*ln*μ**|*) near*∂ω*(see [1] for
details).

Hindawi Publishing Corporation Advances in Diﬀerence Equations Volume 2006, Article ID 70325, Pages1–38 DOI10.1155/ADE/2006/70325

The second problem considered corresponds to the singularly perturbed reaction- diﬀusion problem of parabolic type

*−**μ*^{2}^{}*u**xx*+*u**yy*

+*f*(x,*y,t,u) +u**t**=*0, (x,*y)**∈**ω,t**∈*(0,T],

*f**u**≥*0, (x,*y,t,u)**∈**ω**×*[0,T]*×*(*−∞*,*∞*), (1.2)
where *ω**= {*0*< x <*1*} × {*0*< y <*1*}* and *μ*is a small positive parameter. The initial-
boundary conditions are defined by

*u(x,y, 0)**=**u*^{0}(x,*y),* (x,*y)**∈**ω,*

*u(x,y,t)**=**g(x,y,t),* (x,*y,t)**∈**∂ω**×*(0,T]. (1.3)
The functions *f*,*g*, and*u*^{0}are suﬃciently smooth. Under suitable continuity and com-
patibility conditions on the data, a unique solution*u*of (1.2) exists (see [5] for details).

For*μ*1, problem (1.2) is singularly perturbed and characterized by the boundary lay-
ers of width*O(μ**|*ln*μ**|*) at the boundary*∂ω* (see [2] for details). We mention that the
assumption *f**u**≥*0 in (1.2) can always be obtained via a change of variables.

In solving such nonlinear singularly perturbed problems by the finite diﬀerence method, the corresponding discrete problem is usually formulated as a system of non- linear algebraic equations. One then requires a reliable and eﬃcient computational algo- rithm for computing the solution. A fruitful method for the treatment of these nonlinear systems is the method of upper and lower solutions and its associated monotone itera- tions (in the case of unperturbed problems with reaction-diﬀusion equations see [8,9]

and the references therein). Since the initial iteration in the monotone iterative method is either an upper or lower solution constructed directly from the diﬀerence equations without any knowledge of the exact solution (see [3,4] for details), this method elimi- nates the search for the initial iteration as is often needed in Newton’s method. This gives a practical advantage in the computation of numerical solutions.

Iterative domain decomposition algorithms based on Schwarz-type alternating proce- dures have received much attention for their potential as eﬃcient algorithms for parallel computing. In [3,4], for solving the nonlinear problems (1.1) and (1.2), respectively, we proposed discrete iterative algorithms which combine the monotone approach and an iterative domain decomposition method based on the Schwarz alternating procedure.

The spatial computational domain is partitioned into many nonoverlapping subdomains
(vertical strips) with interface *γ. Small interfacial subdomains are introduced near the*
interface*γ, and approximate boundary values computed onγ*are used for solving prob-
lems on nonoverlapping subdomains. Thus, this approach may be considered as a vari-
ant of a block Gauss-Seidel iteration (or in the parallel context as a multicoloured al-
gorithm) for the subdomains with a Dirichlet-Dirichlet coupling through the interface
variables. In this paper, we generalize the monotone domain decomposition algorithms
from [3,4] and employ a box-domain decomposition of the spatial computational do-
main. This leads to vertical and horizontal interfaces*γ*and*ρ, and corresponding vertical*
and horizontal interfacial subdomain problems provide Dirichlet data on*γ*and*ρ*for the
problems on the nonoverlapping box-subdomains.

InSection 2, we introduce the classical nonlinear finite diﬀerence schemes for the nu- merical solution of (1.1) and (1.2). Iterative methods by which each of these schemes may be solved are presented in [3,4]. From an arbitrary initial mesh function, one may construct a sequence of functions which converges monotonically to the exact solution of the nonlinear diﬀerence scheme. Each function in the sequence is generated as the so- lution of a linear diﬀerence problem. InSection 3, we consider the elliptic problem and extend the monotone method to a box-decomposition of the computational domain.

We show that monotonic convergence is maintained under the proposed decomposition and associated algorithm. Further, we develop estimates of the rate of convergence. The box-decomposition of the spatial domain is applied to the parabolic nonlinear diﬀerence scheme inSection 4. Numerical experiments are presented inSection 5. These confirm the theoretical estimates of the earlier sections. Suggestions are made regarding future parallel implementation.

**2. Diﬀerence schemes for solving (1.1) and (1.2)**

On*ω*and [0,*T] introduce nonuniform meshesω*^{h}*=**ω*^{hx}*×**ω** ^{hy}*and

*ω*

*:*

^{τ}*ω*

^{hx}*=*

*x**i*, 0*≤**i**≤**N**x*;*x*0*=*0,*x**N**x**=*1;*h**xi**=**x**i+1**−**x**i*
,
*ω*^{hy}*=*

*y**j*, 0*≤**j**≤**N**y*; *y*0*=*0, *y**N**y**=*1; *h**y j**=**y**j+1**−**y**j*
,
*ω*^{τ}*=*

*t**k**=**kτ*, 0*≤**k**≤**N**τ*,*N**τ**τ**=**T*^{}*.*

(2.1)

For approximation of the elliptic problem (1.1), we use the classical diﬀerence scheme on nonuniform meshes

ᏸ^{h}*U*+*f*(P,U)*=*0, *P**∈**ω** ^{h}*,

*U*

*=*

*g*on

*∂ω*

*, (2.2)*

^{h}whereᏸ^{h}*U*is defined by

ᏸ^{h}*U**= −**μ*^{2}^{}Ᏸ^{2}*x*+Ᏸ^{2}*y*

*U,* (2.3)

andᏰ^{2}*x**U(P),*Ᏸ^{2}*y**U(P) are the central diﬀerence approximations to the second derivatives*
Ᏸ^{2}*x**U**ij**=*

*xi** _{−}*1

*U**i+1,j**−**U**ij*
*h**xi** _{−}*1

*−*

*U**ij**−**U**i**−*1,j
*h**x,i**−*1

* _{−}*1
,
Ᏸ

^{2}

*y*

*U*

*ij*

*=*

*y j**−*1*U**i,j+1**−**U**ij**h**y j**−*1

*−**U**ij**−**U**i,j**−*1*h**y,j**−*1*−*1
,
*xi**=*2^{−}^{1}^{}*h**x,i**−*1+*h**xi*

, *y j**=*2^{−}^{1}^{}*h**y,j**−*1+*h**y j*
,

(2.4)

where*P**=*(x*i*,*y**j*)*∈**ω** ^{h}*and

*U*

*ij*

*=*

*U(x*

*i*,

*y*

*j*).

To approximate the parabolic problem (1.2), we use the implicit diﬀerence scheme
ᏸ^{hτ}*U(P,t) +* *f*(P,t,U)*=**τ*^{−}^{1}*U(P,t**−**τ),* (P,t)*∈**ω*^{h}*×**ω** ^{τ}*,

ᏸ^{hτ}*U(P,t)**≡*ᏸ^{h}*U(P,t) +τ*^{−}^{1}*U(P,t),*

*U(P, 0)**=**u*^{0}(P), *P**∈**ω** ^{h}*,

*U(P,t)*

*=*

*g(P,t),*(P,

*t)*

*∈*

*∂ω*

^{h}*×*

*ω*

*,*

^{τ}(2.5)

whereᏸ* ^{h}*is defined in (2.3).

Consider the linear versions of problems (2.2) and (2.5)
ᏸW+*c(P*)W(P)*=**F(P*), *P**∈**ω** ^{h}*,

*W(P)**=**W*^{0}(P), *P**∈**∂ω** ^{h}*,

*c(P)*

*≥*

*c*0

*>*0,

*P*

*∈*

*ω*

*,*

^{h}*c*0

*=*const,

(2.6)

whereᏸ*=*ᏸ* ^{h}*for (2.2) andᏸ

*=*ᏸ

*for (2.5). Now we formulate the maximum principle for the diﬀerence operatorᏸ+*

^{hτ}*c*and give an estimate of the solution to (2.6).

*Lemma 2.1. (i) IfW(P) satisfies the conditions*

ᏸW+*c(P)W(P)**≥*0(*≤*0), *P**∈**ω** ^{h}*,

*W(P)*

*≥*0(

*≤*0),

*P*

*∈*

*∂ω*

*, (2.7)*

^{h}*thenW(P*)

*≥*0(

*≤*

*0),P*

*∈*

*ω*

^{h}*.*

*(ii) The following estimate of the solution to (2.6) holds true*
*W* _{ω}^{h}*≤*max^{ }*W*^{0}^{ }_{∂ω}*h*, *F* *ω*^{h}*/*^{}*c*0+*βτ*^{−}^{1}^{},
*W*^{0}^{ }_{∂ω}*h**≡*max

*P*_{∈}*∂ω*^{h}*W*^{0}(P)^{}, *F* *ω*^{h}*≡*max

*P*_{∈}*ω*^{h}*F(P)*^{}, (2.8)
*whereβ**=**0 for (2.2) andβ**=**1 for (2.5).*

The proof of the lemma can be found in [11].

**3. Monotone domain decomposition algorithm for the elliptic problem (1.1)**

We consider a rectangular decomposition of the spatial domain ¯*ω*into (M*×**L) nonover-*
lapping subdomains*ω**ml*,*m**=*1,*...,M*,*l**=*1,*...,L:*

*ω**ml**=*

*x**m**−*1,x*m*

*×*
*y**l**−*1,*y**l*

, *x*0*=*0, *x**M**=*1, *y*0*=*0, *y**L**=*1. (3.1)
Additionally, we introduce (M*−*1) interfacial subdomains*θ**m*,*m**=*1,...,M*−*1 (ver-
tical strips):

*θ**m**=**θ*^{x}_{m}*×**ω*^{y}*=*

*x*_{m}^{b}*< x < x*^{e}_{m}^{}*× {*0*< y <*1*}*, *θ**m**−*1*∩**θ**m**= ∅*,
*γ*_{m}^{b}*=*

*x**=**x*^{b}* _{m}*, 0

*≤*

*y*

*≤*1

^{},

*γ*

_{m}

^{e}*=*

*x**=**x*_{m}* ^{e}*, 0

*≤*

*y*

*≤*1

^{},

*x*

_{m}

^{b}*< x*

*m*

*< x*

^{e}*,*

_{m}*γ*

^{0}

_{m}*=*

*∂ω*

*∩*

*∂θ*

*m*,

(3.2)

*θ*_{m−1}*x**m−1*

*ω**m,l−1*

*x**m*

*ϑ*_{l−1}*y*^{b}_{l−1}

*y*_{l−1}*y*^{e}_{l−1}

*ω**ml*

*ω**m−1,l* *ω**m+1,l*

*y*^{b}_{l}

*y**l*

*y*^{e}_{l}*ϑ*_{l}*θ*_{m}

*x*^{b}_{m−1}*x*^{e}_{m−1}*x*^{b}_{m}*x*^{e}_{m}*ω*_{m,l+l}

Figure 3.1. Fragment of the domain decomposition.

and (L*−*1) interfacial subdomains*ϑ**l*,*l**=*1,*...,L**−*1 (horizontal strips):

*ϑ**l**=**ω*^{x}*×**ϑ*_{l}^{y}*= {*0*< x <*1*} ×*

*y*_{l}^{b}*< y < y*^{e}_{l}^{}, *ϑ**l**−*1*∩**ϑ**l**= ∅*,
*ρ*^{b}_{l}*=*

0*≤**x**≤*1, *y**=**y*_{l}^{b}^{}, *ρ*^{e}_{l}*=*

0*≤**x**≤*1, *y**=**y*_{l}^{e}^{},
*y*^{b}_{l}*< y**l**< y*_{l}* ^{e}*,

*ρ*

_{l}^{0}

*=*

*∂ω*

*∩*

*∂ϑ*

*l*

*.*

(3.3)

Figure 3.1illustrates a fragment of the domain decomposition.

On*ω**ml*,*m**=*1,...,M,*l**=*1,...,*L;θ**m*,*m**=*1,*...,M**−*1 and*ϑ**l*,*l**=*1,...,*L**−*1, intro-
duce meshes:

*ω*^{h}_{ml}*=**ω**ml**∩**ω** ^{h}*,

*θ*

^{h}

_{m}*=*

*θ*

*m*

*∩*

*ω*

*,*

^{h}*ϑ*

^{h}

_{l}*=*

*ϑ*

*l*

*∩*

*ω*

*,*

^{h}*x*

^{b}*,x*

_{m}*m*,x

^{e}

_{m}^{}

^{M}

_{m}

_{=}

^{−}_{1}

^{1}

*∈*

*ω*

*,*

^{hx}^{}

*y*

_{l}*,*

^{b}*y*

*l*,

*y*

_{l}

^{e}^{}

^{L}

_{l}

_{=}

^{−}_{1}

^{1}

*∈*

*ω*

*,*

^{hy}(3.4)

with*ω** ^{hx}*,

*ω*

*from (2.1).*

^{hy}**3.1. Statement of domain decomposition algorithm. We consider the following domain**
decomposition approach for solving (2.2). On each iterative step, we first solve problems
on the nonoverlapping subdomains*ω*^{h}* _{ml}*,

*m*

*=*1,

*...,M,l*

*=*1,...,Lwith Dirichlet bound- ary conditions passed from the previous iterate. Then Dirichlet data are passed from these subdomains to the vertical and horizontal interfacial subdomains

*θ*

^{h}*m*,

*m*

*=*1,...,M

*−*1 and

*ϑ*

^{h}*,*

_{l}*l*

*=*1,...,L

*−*1, respectively. Problems on the vertical interfacial subdomains are computed. Then Dirichlet data from these subdomains are passed to the horizontal inter- facial subdomains before the corresponding linear problems are solved. Finally, we piece together the solutions on the subdomains.

*Step 1. Initialization: On the whole meshω** ^{h}*, choose an initial mesh function

*V*

^{(0)}(P),

*P*

*∈*

*ω*

*satisfying the boundary conditions*

^{h}*V*

^{(0)}(P)

*=*

*g*(P) on

*∂ω*

*.*

^{h}*Step 2. On subdomainsω*^{h}* _{ml}*,

*m*

*=*1,

*...,M*,

*l*

*=*1,

*...,L, compute mesh functionsV*

_{ml}^{(n+1)}(P) (here the index

*n*stands for a number of iterative steps) satisfying the following diﬀerence problems

ᏸ* ^{h}*+

*c*

^{∗}^{}

*Z*

_{ml}^{(n+1)}

*= −*

*G*

^{(n)}(P),

*P*

*∈*

*ω*

^{h}*,*

_{ml}*G*^{(n)}(P)*≡*ᏸ^{h}*V*^{(n)}+ *f*^{}*P*,V^{(n)}^{}, *Z*^{(n+1)}* _{ml}* (P)

*=*0,

*P*

*∈*

*∂ω*

_{ml}*,*

^{h}*V*

_{ml}^{(n+1)}(P)

*=*

*V*

^{(n)}(P) +

*Z*

^{(n+1)}

*(P),*

_{ml}*P*

*∈*

*ω*

^{h}

_{ml}*.*

(3.5)

*Step 3. On the vertical interfacial subdomainsθ*^{h}* _{m}*,

*m*

*=*1,

*...,M*

*−*1, compute the diﬀer- ence problems

ᏸ* ^{h}*+

*c*

^{∗}^{}

*Z*

_{m}^{(n+1)}

*= −*

*G*

^{(n)}(P),

*P*

*∈*

*θ*

_{m}*,*

^{h}*Z*^{(n+1)}* _{m}* (P)

*=*

⎧⎪

⎪⎪

⎪⎨

⎪⎪

⎪⎪

⎩

0, *P**∈**γ*_{m}* ^{h0}*;

*Z*_{ml}^{(n+1)}(P), *P**∈**γ*_{m}^{hb}*∩**ω*^{h}* _{ml}*,

*l*

*=*1,...,L;

*Z*_{m+1,l}^{(n+1)}(P), *P**∈**γ*_{m}^{he}*∩**ω*^{h}* _{m+1,l}*,

*l*

*=*1,

*...,L,*

*V*

_{m}^{(n+1)}(P)

*=*

*V*

^{(n)}(P) +

*Z*

_{m}^{(n+1)}(P),

*P*

*∈*

*θ*

^{h}*,*

_{m}(3.6)

where we use the notation

*γ*^{h0}_{m}*=**γ*^{0}_{m}*∩**∂ω** ^{h}*,

*γ*

^{hb}

_{m}*=*

*γ*

^{b}

_{m}*∩*

*θ*

^{h}*,*

_{m}*γ*

^{he}

_{m}*=*

*γ*

^{e}

_{m}*∩*

*θ*

^{h}

_{m}*.*(3.7)

*Step 4. On the horizontal interfacial subdomainsϑ*

^{h}*,*

_{l}*l*

*=*1,...,L

*−*1, compute the follow- ing diﬀerence problems

ᏸ* ^{h}*+

*c*

^{∗}*Z*

_{l}^{(n+1)}

*= −*

*G*

^{(n)}(P),

*P*

*∈*

*ϑ*

^{h}*,*

_{l}*Z*_{l}^{(n+1)}(P)*=*

⎧⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎨

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎩

0, *P**∈**ρ*^{h0}* _{l}* ;

*Z*

_{ml}^{(n+1)}(P),

*P*

*∈*

*ρ*^{hb}_{l}*\**θ*^{h}^{}*∩**ω*^{h}* _{ml}*,

*m*

*=*1,

*...,M;*

*Z*_{m,l+1}^{(n+1)}(P), *P**∈*

*ρ*^{he}_{l}*\**θ*^{h}^{}*∩**ω*^{h}* _{m,l+1}*,

*m*

*=*1,...,M;

*Z**m*^{(n+1)}(P), *P**∈**∂ϑ*^{h}_{l}*∩**θ*^{h}*m*,*m**=*1,...,M*−*1,
*V*_{l}^{(n+1)}(P)*=**V*^{(n)}(P) +*Z*^{}^{(n+1)}* _{l}* (P),

*P*

*∈*

*ϑ*

^{h}*l*,

(3.8)

where we use the notation

*θ*^{h}*=*

*M**−*1
*m**=*1

*θ*^{h}*m*, *ϑ*^{h}*=*

*L**−*1
*l**=*1

*ϑ*^{h}*l*,

*ρ*^{h0}_{l}*=**ρ*^{0}_{l}*∩**∂ω** ^{h}*,

*ρ*

_{l}

^{hb}*=*

*ρ*

^{b}

_{l}*∩*

*ϑ*

^{h}*,*

_{l}*ρ*

_{l}

^{he}*=*

*ρ*

^{e}

_{l}*∩*

*ϑ*

^{h}

_{l}*.*

(3.9)

*Step 5. Compute the mesh functionV*^{(n+1)}(P),*P**∈**ω** ^{h}*by piecing together the solutions
on the subdomains

*V*^{(n+1)}(P)*=*

⎧⎪

⎪⎪

⎪⎨

⎪⎪

⎪⎪

⎩

*V*_{ml}^{(n+1)}(P), *P**∈**ω*^{h}_{ml}*\*

*θ*^{h}*∪**ϑ*^{h}^{};

*V**m*^{(n+1)}(P), *P**∈**θ*^{h}_{m}*\**ϑ** ^{h}*,

*m*

*=*1,

*...,M*

*−*1;

*V*_{l}^{(n+1)}(P), *P**∈**ϑ*^{h}* _{l}*,

*l*

*=*1,

*...,L*

*−*1.

(3.10)

*Step 6. Stopping criterion: If a prescribed accuracy is reached, then stop; otherwise go to*
Step 2.

Algorithm (3.5)–(3.10) can be carried out by parallel processing. Steps2, 3, and 4 must be performed sequentially, but on each step, the independent subproblems may be assigned to diﬀerent computational nodes.

*Remark 3.1. We note that the original Schwarz alternating algorithm with overlapping*
subdomains is a purely sequential algorithm. To obtain parallelism, one needs a subdo-
main colouring strategy, so that a set of independent subproblems can be introduced.

The modification of the Schwarz algorithm (3.5)–(3.10) can be considered as an additive Schwarz algorithm.

**3.2. Monotone convergence of algorithm (3.5)–(3.10). Additionally, we assume that** *f*
from (1.1) satisfies the two-sided constraints

0*< c*_{∗}*≤**f**u**≤**c** ^{∗}*,

*c*

*,*

_{∗}*c*

^{∗}*=*const. (3.11) We say that

*V*(P) is an upper solution of (2.2) if it satisfies the inequalities

ᏸ^{h}*V*+*f*(P,*V)**≥*0, *P**∈**ω** ^{h}*,

*V*

*≥*

*g*on

*∂ω*

^{h}*.*(3.12) Similarly,

*V*(P) is called a lower solution if it satisfies the reversed inequalities. Upper and lower solutions satisfy the following inequality

*V*(P)*≤**V*(P), *P**∈**ω** ^{h}*, (3.13)

since by the definitions of lower and upper solutions and the mean-value theorem, for
*δV**=**V**−**V* we have

ᏸ^{h}*δV*+*f**u*(P)δV(P)*≥*0, *P**∈**ω** ^{h}*,

*δV*(P)*≥*0, *P**∈**∂ω** ^{h}*, (3.14)

where *f**u*(P)*≡* *f**u*[P,*V*(P) +Θ(P)δV(P)], 0*<*Θ(P)*<*1. In view of the maximum princi-
ple inLemma 2.1, we conclude (3.13).

The following convergence property of algorithm (3.5)–(3.10) holds true.

*Theorem 3.2. LetV*^{(0)}*andV*^{(0)}*be upper and lower solutions of (2.2), and let* *f*(x,*y,u) sat-*
*isfy (3.11). Then the upper sequence**{**V*^{(n)}*}**generated by (3.5)–(3.10) converges monotoni-*
*cally from above to the unique solutionUof (2.2), and the lower sequence**{**V*^{(n)}*}**generated*
*by (3.5)–(3.10) converges monotonically from below toU:*

*V*^{(0)}*≤**V*^{(n)}*≤**V*^{(n+1)}*≤**U**≤**V*^{(n+1)}*≤**V*^{(n)}*≤**V*^{(0)}, *inω*^{h}*.* (3.15)

*Proof. We consider only the case of the upper sequence. LetV*^{(n)}be an upper solution.

Then by the maximum principle inLemma 2.1, from (3.5) we conclude that

*Z*_{ml}^{(n+1)}(P)*≤*0, *P**∈**ω*^{h}* _{ml}*,

*m*

*=*1,...,M,

*l*

*=*1,

*...,L.*(3.16)

Using the mean-value theorem and the equation for*Z*_{ml}^{(n+1)}(P), we obtain the diﬀerence
equation for*V*_{ml}^{(n+1)}

ᏸ^{h}*V*_{ml}^{(n+1)}+*f*^{}*P,V*_{ml}^{(n+1)}^{}*= −*

*c*^{∗}*−**f*_{u,ml}^{(n)}(P)^{}*Z*_{ml}^{(n+1)}(P)*≥*0, *P**∈**ω*^{h}* _{ml}*,

*f*

_{u,ml}^{(n)}(P)

*≡*

*f*

*u*

*P*,V

^{(n)}(P) +Θ

^{(n)}

*(P)Z*

_{ml}

_{ml}^{(n+1)}(P)

^{}, 0

*<*Θ

^{(n)}

*(P)*

_{ml}*<*1,

*V*_{ml}^{(n+1)}(P)*=**V*^{(n)}(P), *P**∈**∂ω*_{ml}* ^{h}* ,

(3.17)

where nonnegativeness of the right-hand side of the diﬀerence equation follows from (3.11) and (3.16).

Taking into account (3.16) and*V*^{(n)}is an upper solution, by the maximum principle
inLemma 2.1, from (3.6) and (3.8) it follows that

*Z*_{m}^{(n+1)}(P)*≤*0, *P**∈**θ*^{h}* _{m}*,

*m*

*=*1,

*...,M*

*−*1,

*Z*_{l}^{(n+1)}(P)*≤*0, *P**∈**ϑ*^{h}* _{l}*,

*l*

*=*1,

*...,L*

*−*1. (3.18)

Similar to (3.17), we obtain the diﬀerence problems for*V**m*^{(n+1)}

ᏸ^{h}*V*_{m}^{(n+1)}+*f*^{}*P,V*_{m}^{(n+1)}^{}*= −**c*^{∗}*−**f*_{u,m}^{(n)}(P)^{}*Z*_{m}^{(n+1)}(P)*≥*0, *P**∈**θ*_{m}* ^{h}*,

*V*_{m}^{(n+1)}(P)*=*

⎧⎪

⎪⎪

⎪⎨

⎪⎪

⎪⎪

⎩

*g(P),* *P**∈**γ*^{h0}*m*;

*V*_{ml}^{(n+1)}(P), *P**∈**γ*^{hb}_{m}*∩**ω*^{h}* _{ml}*,

*l*

*=*1,

*...,L;*

*V*_{m+1,l}^{(n+1)}(P), *P**∈**γ*^{he}_{m}*∩**ω*^{h}* _{m+1,l}*,

*l*

*=*1,...,L,

(3.19)

and for*V*_{l}^{(n+1)}

ᏸ^{h}*V*^{}_{l}^{(n+1)}+*f*^{}*P,V*^{}_{l}^{(n+1)}^{}*= −*

*c*^{∗}*−**f*_{u,l}^{(n)}(P) *Z*_{l}^{(n+1)}(P)*≥*0, *P**∈**ϑ*^{h}* _{l}*,

*V*_{l}^{(n+1)}(P)*=*

⎧⎪

⎪⎪

⎪⎪

⎪⎪

⎨

⎪⎪

⎪⎪

⎪⎪

⎪⎩

*g*(P), *P**∈**ρ*_{l}* ^{h0}*;

*V*

_{ml}^{(n+1)}(P),

*P*

*∈*

*ρ*^{hb}_{l}*\**θ*^{h}^{}*∩**ω*^{h}* _{ml}*,

*m*

*=*1,...,M;

*V*_{m,l+1}^{(n+1)}(P), *P**∈*

*ρ*^{he}_{l}*\**θ*^{h}^{}*∩**ω*^{h}* _{m,l+1}*,

*m*

*=*1,

*...,M;*

*V**m*^{(n+1)}(P), *P**∈**∂ϑ*^{h}_{l}*∩**θ*^{h}* _{m}*,

*m*

*=*1,...,M

*−*1,

(3.20)

where nonnegativeness of the right-hand sides of the diﬀerence equations follows from
(3.11) and (3.18). Now we verify that the mesh function*V*^{(n+1)}defined by (3.10) is an up-
per solution. From the boundary conditions for*V*_{ml}^{(n+1)},*V**m*^{(n+1)}and*V*^{}_{l}^{(n)}, it follows that
*V*^{(n+1)}satisfies the boundary condition in (2.2). Now from here, (3.17), (3.19), (3.20)
and the definition of*V*^{(n+1)}in (3.10), we conclude that

*G*^{(n+1)}(P)*=*ᏸ^{h}*V*^{(n+1)}+ *f*^{}*P*,V^{(n+1)}^{}*≥*0, *P**∈**ω*^{h}*\*

*γ*^{h}*∪**ρ*^{h}^{},
*γ*^{hb,e}_{ml}*=*

*x**i**=**x*^{b,e}*m* ,*y*^{e}_{l}* _{−}*1

*< y*

*j*

*< y*

_{l}

^{b}^{},

*γ*

^{hb,e}*m*

*=*

*L*

*l*

*=*1

*γ*^{hb,e}* _{ml}* ,

*y*0

^{e}*=*0,

*y*

*L*

^{b}*=*1,

*γ*

^{h}*=*

*M**−*1
*m**=*1

*γ*^{hb,e}* _{m}* ,

*ρ*

^{h}*=*

*L**−*1
*l**=*1

*ρ*^{hb,e}_{l}*.*

(3.21)

To prove that*V*^{(n+1)}is an upper solution of problem (2.2), we have to verify only that
the last inequality holds true on the interfacial boundaries*γ*_{ml}* ^{hb,e}*and

*ρ*

^{hb,e}*,*

_{l}*m*

*=*1,...,

*M*

*−*1,

*l*

*=*1,

*...,L*

*−*1.

We check this inequality in the case of the left interfacial boundary*γ*^{hb}* _{ml}*, since the case
with

*γ*

^{he}*is checked in a similar way. From (3.5), (3.6), and (3.18), we conclude that the mesh function*

_{ml}*W*

_{ml}^{(n+1)}

*=*

*V*

_{ml}^{(n+1)}

*−*

*V*

*m*

^{(n+1)}satisfies the diﬀerence problem

ᏸ* ^{h}*+

*c*

^{∗}^{}

*W*

_{ml}^{(n+1)}

*=*0,

*P*

*∈*

*θ*

_{ml}

^{h}*=*

*ω*

^{h}

_{ml}*∩*

*θ*

_{m}*,*

^{h}*W*

_{ml}^{(n+1)}(P)

*=*

⎧⎨

⎩

0, *P**∈**γ*^{hb}_{ml}*=**γ**m*^{hb}*∩**ω*^{h}* _{ml}*;

*≥*0, *P**∈**∂θ*_{ml}^{h}*\**γ*^{hb}_{ml}*.*

(3.22)

In view of the maximum principle inLemma 2.1,

*V*_{ml}^{(n+1)}(P)*−**V*_{m}^{(n+1)}(P)*≥*0, *P**∈**θ*^{h}_{ml}*.* (3.23)
By (3.6),*V**m*^{(n+1)}(P)*=**V*_{ml}^{(n+1)}(P),*P**∈**γ*^{hb}* _{ml}*, and from (3.10) and (3.23), it follows that

*−**μ*^{2}Ᏸ^{2}*y**V*_{ml}^{(n+1)}(P)*= −**μ*^{2}Ᏸ^{2}*y**V*^{(n+1)}(P), *P**∈**γ*_{ml}* ^{hb}*,

*−**μ*^{2}Ᏸ^{2}*x**V*_{ml}^{(n+1)}(P)*≤ −**μ*^{2}Ᏸ^{2}*x**V*^{(n+1)}(P), *P**∈**γ*^{hb}_{ml}*.* (3.24)

Thus, using (3.17), we conclude

*G*^{(n+1)}(P)*≥*ᏸ^{h}*V*_{ml}^{(n+1)}(P) +*f*^{}*P,V*_{ml}^{(n+1)}^{}*≥*0, *P**∈**γ*_{ml}^{hb}*.* (3.25)
Now we verify the inequality*G*^{(n+1)}(P)*≥*0 on the interfacial boundary*ρ*^{hb}* _{l}* , and the
case with

*ρ*

^{he}*is checked in a similar way. From (3.5), (3.8), (3.18), and (3.23), we conclude that the mesh function*

_{l}*W*

^{}

_{ml}^{(n+1)}

*=*

*V*

_{ml}^{(n+1)}

*−*

*V*

_{l}^{(n+1)}satisfies the diﬀerence problem

ᏸ* ^{h}*+

*c*

^{∗}*W*

_{ml}^{(n+1)}

*=*0,

*P*

*∈*

*ϑ*

^{h}

_{ml}*=*

*ω*

^{h}

_{ml}*∩*

*ϑ*

^{h}*,*

_{l}*W*_{ml}^{(n+1)}(P)*=*

⎧⎨

⎩

0, *P**∈**ρ*^{hb}_{ml}*=*

*x*^{e}_{m}* _{−}*1

*< x*

*i*

*< x*

^{b}*,*

_{m}*y*

*j*

*=*

*y*

_{l}

^{b}^{};

*≥*0, *P**∈**∂ϑ*^{h}_{ml}*\**ρ*^{hb}_{ml}*.*

(3.26)

By the maximum principle inLemma 2.1,

*V*_{ml}^{(n+1)}(P)*−**V*_{l}^{(n+1)}(P)*≥*0, *P**∈**ϑ*^{h}_{ml}*.* (3.27)
By (3.8), *V*^{}_{l}^{(n+1)}(P)*=**V*_{ml}^{(n+1)}(P), *P**∈**ρ*^{hb}_{ml}*∪ {*(x^{e}_{m}* _{−}*1,

*y*

_{l}*), (x*

^{b}

^{b}*,*

_{m}*y*

_{l}*)*

^{b}*}*, and from (3.10) and (3.27), it follows that

*−**μ*^{2}Ᏸ^{2}*x**V*_{ml}^{(n+1)}(P)*= −**μ*^{2}Ᏸ^{2}*x**V*^{(n+1)}(P), *P**∈**ρ*^{hb}* _{ml}*,

*−**μ*^{2}Ᏸ^{2}*y**V*_{ml}^{(n+1)}(P)*≤ −**μ*^{2}Ᏸ^{2}*y**V*^{(n+1)}(P), *P**∈**ρ*^{hb}_{ml}*.* (3.28)
Thus, using (3.17), we conclude

*G*^{(n+1)}(P)*≥*ᏸ^{h}*V*_{ml}^{(n+1)}+*f*^{}*P,V*_{ml}^{(n+1)}^{}*≥*0, *P**∈**ρ*^{hb}_{ml}*.* (3.29)
From (3.6), (3.8), and (3.27), the mesh function ˆ*W*_{ml}^{(n+1)}*=**V**m*^{(n+1)}*−**V*_{l}^{(n+1)}satisfies the
diﬀerence problem

ᏸ* ^{h}*+

*c*

^{∗}*W*

_{ml}^{(n+1)}

*=*0,

*P*

*∈*

*τ*

_{ml}

^{h}*=*

*θ*

^{h}

_{m}*∩*

*ϑ*

^{h}*,*

_{l}*W*

_{ml}^{(n+1)}(P)

*=*

⎧⎨

⎩

0, *P**∈**ρ*^{hb,e}_{ml}*=**x*^{b}_{m}*< x**i**< x*^{e}* _{m}*,

*y*

*j*

*=*

*y*

^{b,e}

_{l}^{};

*≥*0, *P**∈**∂τ*_{ml}^{h}*\*

*ρ*^{hb}_{ml}*∪**ρ*^{he}_{ml}^{}*.*

(3.30)

By the maximum principle inLemma 2.1,

*V**m*^{(n+1)}(P)*−**V*_{l}^{(n+1)}(P)*≥*0, *P**∈**τ*^{h}_{ml}*.* (3.31)
By (3.8),*V*^{}_{l}^{(n+1)}(P)*=**V**m*^{(n+1)}(P),*P**∈**ρ*^{hb}_{ml}*∪ {*(x^{e}* _{m}*,

*y*

_{l}*), (x*

^{b}

^{b}*,*

_{m}*y*

^{b}*)*

_{l}*}*, and from (3.10) and (3.31), it follows that

*−**μ*^{2}Ᏸ^{2}*x**V**m*^{(n+1)}(P)*= −**μ*^{2}Ᏸ^{2}*x**V*^{(n+1)}(P), *P**∈**ρ*^{hb}* _{ml}*,

*−**μ*^{2}Ᏸ^{2}*y**V*_{m}^{(n+1)}(P)*≤ −**μ*^{2}Ᏸ^{2}*y**V*^{(n+1)}(P), *P**∈**ρ*^{hb}_{ml}*.* (3.32)

Thus, using (3.19), we conclude

*G*^{(n+1)}(P)*≥*ᏸ^{h}*V**m*^{(n+1)}+*f*^{}*P,V**m*^{(n+1)}

*≥*0, *P**∈**ρ*^{hb}_{ml}*.* (3.33)
From here and (3.29), we conclude the required inequality on*ρ*^{hb}_{l}*\**P*^{b,e}* _{l}* ,

*P*

^{b,e}

_{l}*= ∪*

^{M}*m*

*=*

*1*

^{−}^{1}(x

^{b,e}*m*,

*y*

_{l}*). At*

^{b}*P*

^{b}

_{ml}*=*(x

^{b}*,*

_{m}*y*

_{l}*), we have*

^{b}*V*_{ml}^{(n+1)}^{}*P*_{ml}^{b}^{}*=**V*_{m}^{(n+1)}^{}*P*_{ml}^{b}^{}*=**V*_{l}^{(n+1)}^{}*P*_{ml}^{b}^{}, (3.34)
and from (3.10), it follows that

*−**μ*^{2}Ᏸ^{2}*x**V*^{(n+1)}^{}*P*_{ml}^{b}^{}*=−* *μ*^{2}
^{b}*xm*

*V**m*^{(n+1)}

*P*_{ml}^{bx+}^{}*−**V*_{ml}^{(n+1)}^{}*P*_{ml}^{b}^{}

*h*^{b+}_{xm}^{−}

*V*_{ml}^{(n+1)}^{}*P*_{ml}^{b}^{}*−**V*_{ml}^{(n+1)}^{}*P*_{ml}^{bx}^{−}^{}
*h*^{b}_{xm}^{−}

,

*−**μ*^{2}Ᏸ^{2}*y**V*^{(n+1)}^{}*P*_{ml}^{b}^{}*=−**μ*^{2}
^{b}_{yl}

*V*_{l}^{(n+1)}^{}*P*^{by+}_{ml}^{}*−**V*_{ml}^{(n+1)}^{}*P*_{ml}^{b}^{}

*h*^{b+}_{yl}^{−}

*V*_{ml}^{(n+1)}^{}*P*_{ml}^{b}^{}*−**V*_{ml}^{(n+1)}^{}*P*^{by}_{ml}^{−}^{}
*h*^{b}_{yl}^{−}

,
*P*_{ml}^{bx}^{±}*=*

*x*^{b}_{m}*±**h*^{b}_{xm}* ^{±}*,

*y*

^{b}

_{l}^{},

*P*

_{ml}

^{by}

^{±}*=*

*x*^{b}* _{m}*,

*y*

_{l}

^{b}*±*

*h*

^{b}

_{yl}

^{±}^{},

^{b}

_{xm}*=*2

^{−}^{1}

^{}

*h*

^{b}

_{xm}*+*

^{−}*h*

^{b+}

_{xm}^{},

^{b}

_{yl}*=*2

^{−}^{1}

^{}

*h*

^{b}

_{yl}*+*

^{−}*h*

^{b+}

_{yl}^{},

(3.35)
where*h*^{b+}* _{xm}*,

*h*

^{b}

_{xm}*are the mesh step sizes on the left and right from*

^{−}*P*

_{ml}*, and*

^{b}*h*

^{b+}*,*

_{yl}*h*

^{b}

_{yl}*are the mesh step sizes on the top and bottom from*

^{−}*P*

^{b}*. From here, (3.17), (3.23) and (3.27), we conclude*

_{ml}*G*^{(n+1)}(P)*≥*ᏸ^{h}*V*_{ml}^{(n+1)}+*f*^{}*P,V*_{ml}^{(n+1)}^{}*≥*0, *P**=**P*^{b}_{ml}*.* (3.36)
With a similar argument for mesh point*P*_{l}^{e}*=*(x^{e}*m*,*y*_{l}* ^{b}*), we prove that

*V*

^{(n+1)}is an upper solution of problem (2.2) on the whole computational domain

*ω*

*.*

^{h}For arbitrary *P**∈**ω** ^{h}*, it follows from (3.16), (3.18), and (3.13) that the sequence

*{*

*V*

^{(n)}(P)

*}*is monotonically decreasing and bounded below by

*V(P), where*

*V*is any lower solution. Therefore, the sequence is convergent and it follows from (3.5)–(3.8) that limZ

_{ml}^{(n)}

*=*0, lim

*Z*

_{l}^{(n)}

*=*0 and lim

*Z*

^{}

_{l}^{(n)}

*=*0 as

*n*

*→ ∞*. Now by linearity of the operatorᏸ

*and the continuity of*

^{h}*f*, we have also from (3.5)–(3.8) that the mesh function

*U*defined by

*U(P)**=** _{n}*lim

*→∞**V*^{(n)}(P), *P**∈**ω** ^{h}*, (3.37)

is an exact solution to (2.2). The uniqueness of the solution to (2.2) follows from esti-
mate (2.8). Indeed, if by contradiction, we assume that there exist two solutions*U*1and
*U*2to (2.8), then by the mean-value theorem, the diﬀerence*δU**=**U*1*−**U*2satisfies the
following diﬀerence problem

ᏸ^{h}*δU*+*f**u**δU**=*0, *P**∈**ω** ^{h}*,

*δU*

*=*0,

*P*

*∈*

*∂ω*

^{h}*.*(3.38)