MONOTONE FINITE DIFFERENCE DOMAIN DECOMPOSITION ALGORITHMS AND

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 difference iterative algorithms for solving non- linear singularly perturbed reaction-diffusion 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-diffusion problems.

The first problem considered corresponds to the singularly perturbed reaction-diffu- sion problem of elliptic type

−μ²uxx+uyy

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

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

(1.1)

whereμis a small positive parameter,c∗>0 is a constant,∂ω is the boundary ofω. If f andgare sufficiently smooth, then under suitable continuity and compatibility condi- tions on the data, a unique solutionuof (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 widthO(μ|lnμ|) near∂ω(see [1] for details).

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

The second problem considered corresponds to the singularly perturbed reaction- diffusion problem of parabolic type

−μ²uxx+uyy

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

fu≥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⁰(x,y), (x,y)∈ω,

u(x,y,t)=g(x,y,t), (x,y,t)∈∂ω×(0,T]. (1.3) The functions f,g, andu⁰are sufficiently smooth. Under suitable continuity and com- patibility conditions on the data, a unique solutionuof (1.2) exists (see [5] for details).

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

In solving such nonlinear singularly perturbed problems by the finite difference method, the corresponding discrete problem is usually formulated as a system of non- linear algebraic equations. One then requires a reliable and efficient 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-diffusion 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 difference 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 efficient 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 difference 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 difference scheme. Each function in the sequence is generated as the so- lution of a linear difference 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 difference 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. Difference schemes for solving (1.1) and (1.2)

Onωand [0,T] introduce nonuniform meshesω^h=ω^hx×ω^hyandω^τ: ω^hx=

xi, 0≤i≤Nx;x0=0,xNx=1;hxi=xi+1−xi , ω^hy=

yj, 0≤j≤Ny; y0=0, yNy=1; hy j=yj+1−yj , ω^τ=

tk=kτ, 0≤k≤Nτ,Nττ=T.

(2.1)

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

ᏸ^hU+f(P,U)=0, P∈ω^h,U=gon∂ω^h, (2.2)

whereᏸ^hUis defined by

ᏸ^hU= −μ²Ᏸ²x+Ᏸ²y

U, (2.3)

andᏰ²xU(P),Ᏸ²yU(P) are the central difference approximations to the second derivatives Ᏸ²xUij=

xi₋1

Ui+1,j−Uij hxi₋1

−

Uij−Ui−1,j hx,i−1

₋1 , Ᏸ²yUij=

y j−1Ui,j+1−Uijhy j−1

−Uij−Ui,j−1hy,j−1−1 , xi=2⁻¹hx,i−1+hxi

, y j=2⁻¹hy,j−1+hy j ,

(2.4)

whereP=(xi,yj)∈ω^handUij=U(xi,yj).

To approximate the parabolic problem (1.2), we use the implicit difference scheme ᏸ^hτU(P,t) + f(P,t,U)=τ⁻¹U(P,t−τ), (P,t)∈ω^h×ω^τ,

ᏸ^hτU(P,t)≡ᏸ^hU(P,t) +τ⁻¹U(P,t),

U(P, 0)=u⁰(P), P∈ω^h, U(P,t)=g(P,t), (P,t)∈∂ω^h×ω^τ,

(2.5)

whereᏸ^his 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⁰(P), P∈∂ω^h, c(P)≥c0>0, P∈ω^h,c0=const,

(2.6)

whereᏸ=ᏸ^hfor (2.2) andᏸ=ᏸ^hτfor (2.5). Now we formulate the maximum principle for the difference operatorᏸ+cand 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∈∂ω^h, (2.7) thenW(P)≥0(≤0),P∈ω^h.

(ii) The following estimate of the solution to (2.6) holds true W _ω^h≤max W⁰ _∂ωh, F ω^h/c0+βτ⁻¹, W⁰ _∂ωh≡max

P_∈∂ω^hW⁰(P), F ω^h≡max

P_∈ω^hF(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=

xm−1,xm

× yl−1,yl

, x0=0, xM=1, y0=0, yL=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 < xm< x^e_m, γ⁰_m=∂ω∩∂θm,

(3.2)

θ_m−1 xm−1

ωm,l−1

xm

ϑ_l−1 y^b_l−1

y_l−1 y^e_l−1

ωml

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

y^b_l

yl

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 < yl< y_l^e, ρ_l⁰=∂ω∩∂ϑ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_m,xm,x^e_m^M_m=⁻₁¹∈ω^hx, y_l^b,yl,y_l^eL_l=⁻₁¹∈ω^hy,

(3.4)

withω^hx,ω^hyfrom (2.1).

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θ^hm,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 functionV⁽⁰⁾(P), P∈ω^hsatisfying the boundary conditionsV⁽⁰⁾(P)=g(P) on∂ω^h.

Step 2. On subdomainsω^h_ml,m=1,...,M,l=1,...,L, compute mesh functionsV_ml⁽ⁿ⁺¹⁾(P) (here the indexnstands for a number of iterative steps) satisfying the following difference problems

ᏸ^h+c^∗Z_ml⁽ⁿ⁺¹⁾= −G⁽ⁿ⁾(P), P∈ω^h_ml,

G⁽ⁿ⁾(P)≡ᏸ^hV⁽ⁿ⁾+ fP,V⁽ⁿ⁾, Z⁽ⁿ⁺¹⁾_ml (P)=0, P∈∂ω_ml^h , V_ml⁽ⁿ⁺¹⁾(P)=V⁽ⁿ⁾(P) +Z⁽ⁿ⁺¹⁾_ml (P), P∈ω^h_ml.

(3.5)

Step 3. On the vertical interfacial subdomainsθ^h_m,m=1,...,M−1, compute the differ- ence problems

ᏸ^h+c^∗Z_m⁽ⁿ⁺¹⁾= −G⁽ⁿ⁾(P), P∈θ_m^h,

Z⁽ⁿ⁺¹⁾_m (P)=

⎧⎪

⎪⎪

⎪⎨

⎪⎪

⎪⎪

⎩

0, P∈γ_m^h0;

Z_ml⁽ⁿ⁺¹⁾(P), P∈γ_m^hb∩ω^h_ml,l=1,...,L;

Z_m+1,l⁽ⁿ⁺¹⁾(P), P∈γ_m^he∩ω^h_m+1,l,l=1,...,L, V_m⁽ⁿ⁺¹⁾(P)=V⁽ⁿ⁾(P) +Z_m⁽ⁿ⁺¹⁾(P), P∈θ^h_m,

(3.6)

where we use the notation

γ^h0_m =γ⁰_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 difference problems

ᏸ^h+c^∗ Z_l⁽ⁿ⁺¹⁾= −G⁽ⁿ⁾(P), P∈ϑ^h_l,

Z_l⁽ⁿ⁺¹⁾(P)=

⎧⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎨

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎩

0, P∈ρ^h0_l ; Z_ml⁽ⁿ⁺¹⁾(P), P∈

ρ^hb_l \θ^h∩ω^h_ml,m=1,...,M;

Z_m,l+1⁽ⁿ⁺¹⁾(P), P∈

ρ^he_l \θ^h∩ω^h_m,l+1,m=1,...,M;

Zm⁽ⁿ⁺¹⁾(P), P∈∂ϑ^h_l∩θ^hm,m=1,...,M−1, V_l⁽ⁿ⁺¹⁾(P)=V⁽ⁿ⁾(P) +Z⁽ⁿ⁺¹⁾_l (P), P∈ϑ^hl,

(3.8)

where we use the notation

θ^h=

M−1 m=1

θ^hm, ϑ^h=

L−1 l=1

ϑ^hl,

ρ^h0_l =ρ⁰_l∩∂ω^h, ρ_l^hb=ρ^b_l∩ϑ^h_l, ρ_l^he=ρ^e_l∩ϑ^h_l.

(3.9)

Step 5. Compute the mesh functionV⁽ⁿ⁺¹⁾(P),P∈ω^hby piecing together the solutions on the subdomains

V⁽ⁿ⁺¹⁾(P)=

⎧⎪

⎪⎪

⎪⎨

⎪⎪

⎪⎪

⎩

V_ml⁽ⁿ⁺¹⁾(P), P∈ω^h_ml\

θ^h∪ϑ^h;

Vm⁽ⁿ⁺¹⁾(P), P∈θ^h_m\ϑ^h,m=1,...,M−1;

V_l⁽ⁿ⁺¹⁾(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 different 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_∗≤fu≤c^∗, c_∗,c^∗=const. (3.11) We say thatV(P) is an upper solution of (2.2) if it satisfies the inequalities

ᏸ^hV+f(P,V)≥0, P∈ω^h,V≥gon∂ω^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+fu(P)δV(P)≥0, P∈ω^h,

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

where fu(P)≡ fu[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⁽⁰⁾andV⁽⁰⁾be upper and lower solutions of (2.2), and let f(x,y,u) sat- isfy (3.11). Then the upper sequence{V⁽ⁿ⁾}generated by (3.5)–(3.10) converges monotoni- cally from above to the unique solutionUof (2.2), and the lower sequence{V⁽ⁿ⁾}generated by (3.5)–(3.10) converges monotonically from below toU:

V⁽⁰⁾≤V⁽ⁿ⁾≤V⁽ⁿ⁺¹⁾≤U≤V⁽ⁿ⁺¹⁾≤V⁽ⁿ⁾≤V⁽⁰⁾, inω^h. (3.15)

Proof. We consider only the case of the upper sequence. LetV⁽ⁿ⁾be an upper solution.

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

Z_ml⁽ⁿ⁺¹⁾(P)≤0, P∈ω^h_ml,m=1,...,M,l=1,...,L. (3.16)

Using the mean-value theorem and the equation forZ_ml⁽ⁿ⁺¹⁾(P), we obtain the difference equation forV_ml⁽ⁿ⁺¹⁾

ᏸ^hV_ml⁽ⁿ⁺¹⁾+fP,V_ml⁽ⁿ⁺¹⁾= −

c^∗−f_u,ml⁽ⁿ⁾(P)Z_ml⁽ⁿ⁺¹⁾(P)≥0, P∈ω^h_ml, f_u,ml⁽ⁿ⁾(P)≡fuP,V⁽ⁿ⁾(P) +Θ⁽ⁿ⁾_ml(P)Z_ml⁽ⁿ⁺¹⁾(P), 0<Θ⁽ⁿ⁾_ml(P)<1,

V_ml⁽ⁿ⁺¹⁾(P)=V⁽ⁿ⁾(P), P∈∂ω_ml^h ,

(3.17)

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

Taking into account (3.16) andV⁽ⁿ⁾is an upper solution, by the maximum principle inLemma 2.1, from (3.6) and (3.8) it follows that

Z_m⁽ⁿ⁺¹⁾(P)≤0, P∈θ^h_m,m=1,...,M−1,

Z_l⁽ⁿ⁺¹⁾(P)≤0, P∈ϑ^h_l,l=1,...,L−1. (3.18)

Similar to (3.17), we obtain the difference problems forVm⁽ⁿ⁺¹⁾

ᏸ^hV_m⁽ⁿ⁺¹⁾+fP,V_m⁽ⁿ⁺¹⁾= −c^∗−f_u,m⁽ⁿ⁾(P)Z_m⁽ⁿ⁺¹⁾(P)≥0, P∈θ_m^h,

V_m⁽ⁿ⁺¹⁾(P)=

⎧⎪

⎪⎪

⎪⎨

⎪⎪

⎪⎪

⎩

g(P), P∈γ^h0m;

V_ml⁽ⁿ⁺¹⁾(P), P∈γ^hb_m ∩ω^h_ml,l=1,...,L;

V_m+1,l⁽ⁿ⁺¹⁾(P), P∈γ^he_m∩ω^h_m+1,l,l=1,...,L,

(3.19)

and forV_l⁽ⁿ⁺¹⁾

ᏸ^hV_l⁽ⁿ⁺¹⁾+fP,V_l⁽ⁿ⁺¹⁾= −

c^∗−f_u,l⁽ⁿ⁾(P) Z_l⁽ⁿ⁺¹⁾(P)≥0, P∈ϑ^h_l,

V_l⁽ⁿ⁺¹⁾(P)=

⎧⎪

⎪⎪

⎪⎪

⎪⎪

⎨

⎪⎪

⎪⎪

⎪⎪

⎪⎩

g(P), P∈ρ_l^h0; V_ml⁽ⁿ⁺¹⁾(P), P∈

ρ^hb_l \θ^h∩ω^h_ml,m=1,...,M;

V_m,l+1⁽ⁿ⁺¹⁾(P), P∈

ρ^he_l \θ^h∩ω^h_m,l+1,m=1,...,M;

Vm⁽ⁿ⁺¹⁾(P), P∈∂ϑ^h_l ∩θ^h_m,m=1,...,M−1,

(3.20)

where nonnegativeness of the right-hand sides of the difference equations follows from (3.11) and (3.18). Now we verify that the mesh functionV⁽ⁿ⁺¹⁾defined by (3.10) is an up- per solution. From the boundary conditions forV_ml⁽ⁿ⁺¹⁾,Vm⁽ⁿ⁺¹⁾andV_l⁽ⁿ⁾, it follows that V⁽ⁿ⁺¹⁾satisfies the boundary condition in (2.2). Now from here, (3.17), (3.19), (3.20) and the definition ofV⁽ⁿ⁺¹⁾in (3.10), we conclude that

G⁽ⁿ⁺¹⁾(P)=ᏸ^hV⁽ⁿ⁺¹⁾+ fP,V⁽ⁿ⁺¹⁾≥0, P∈ω^h\

γ^h∪ρ^h, γ^hb,e_ml =

xi=x^b,em ,y^e_l−1< yj< y_l^b, γ^hb,em = L l=1

γ^hb,e_ml , y0^e=0, yL^b=1, γ^h=

M−1 m=1

γ^hb,e_m , ρ^h=

L−1 l=1

ρ^hb,e_l .

(3.21)

To prove thatV⁽ⁿ⁺¹⁾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,eandρ^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_mlis checked in a similar way. From (3.5), (3.6), and (3.18), we conclude that the mesh functionW_ml⁽ⁿ⁺¹⁾=V_ml⁽ⁿ⁺¹⁾−Vm⁽ⁿ⁺¹⁾satisfies the difference problem

ᏸ^h+c^∗W_ml⁽ⁿ⁺¹⁾=0, P∈θ_ml^h =ω^h_ml∩θ_m^h, W_ml⁽ⁿ⁺¹⁾(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⁽ⁿ⁺¹⁾(P)−V_m⁽ⁿ⁺¹⁾(P)≥0, P∈θ^h_ml. (3.23) By (3.6),Vm⁽ⁿ⁺¹⁾(P)=V_ml⁽ⁿ⁺¹⁾(P),P∈γ^hb_ml, and from (3.10) and (3.23), it follows that

−μ²Ᏸ²yV_ml⁽ⁿ⁺¹⁾(P)= −μ²Ᏸ²yV⁽ⁿ⁺¹⁾(P), P∈γ_ml^hb,

−μ²Ᏸ²xV_ml⁽ⁿ⁺¹⁾(P)≤ −μ²Ᏸ²xV⁽ⁿ⁺¹⁾(P), P∈γ^hb_ml. (3.24)

Thus, using (3.17), we conclude

G⁽ⁿ⁺¹⁾(P)≥ᏸ^hV_ml⁽ⁿ⁺¹⁾(P) +fP,V_ml⁽ⁿ⁺¹⁾≥0, P∈γ_ml^hb. (3.25) Now we verify the inequalityG⁽ⁿ⁺¹⁾(P)≥0 on the interfacial boundaryρ^hb_l , and the case withρ^he_l is checked in a similar way. From (3.5), (3.8), (3.18), and (3.23), we conclude that the mesh functionW_ml⁽ⁿ⁺¹⁾=V_ml⁽ⁿ⁺¹⁾−V_l⁽ⁿ⁺¹⁾satisfies the difference problem

ᏸ^h+c^∗ W_ml⁽ⁿ⁺¹⁾=0, P∈ϑ^h_ml=ω^h_ml∩ϑ^h_l,

W_ml⁽ⁿ⁺¹⁾(P)=

⎧⎨

⎩

0, P∈ρ^hb_ml=

x^e_m−1< xi< x^b_m,yj=y_l^b;

≥0, P∈∂ϑ^h_ml\ρ^hb_ml.

(3.26)

By the maximum principle inLemma 2.1,

V_ml⁽ⁿ⁺¹⁾(P)−V_l⁽ⁿ⁺¹⁾(P)≥0, P∈ϑ^h_ml. (3.27) By (3.8), V_l⁽ⁿ⁺¹⁾(P)=V_ml⁽ⁿ⁺¹⁾(P), P∈ρ^hb_ml∪ {(x^e_m−1,y_l^b), (x^b_m,y_l^b)}, and from (3.10) and (3.27), it follows that

−μ²Ᏸ²xV_ml⁽ⁿ⁺¹⁾(P)= −μ²Ᏸ²xV⁽ⁿ⁺¹⁾(P), P∈ρ^hb_ml,

−μ²Ᏸ²yV_ml⁽ⁿ⁺¹⁾(P)≤ −μ²Ᏸ²yV⁽ⁿ⁺¹⁾(P), P∈ρ^hb_ml. (3.28) Thus, using (3.17), we conclude

G⁽ⁿ⁺¹⁾(P)≥ᏸ^hV_ml⁽ⁿ⁺¹⁾+fP,V_ml⁽ⁿ⁺¹⁾≥0, P∈ρ^hb_ml. (3.29) From (3.6), (3.8), and (3.27), the mesh function ˆW_ml⁽ⁿ⁺¹⁾=Vm⁽ⁿ⁺¹⁾−V_l⁽ⁿ⁺¹⁾satisfies the difference problem

ᏸ^h+c^∗ W_ml⁽ⁿ⁺¹⁾=0, P∈τ_ml^h =θ^h_m∩ϑ^h_l, W_ml⁽ⁿ⁺¹⁾(P)=

⎧⎨

⎩

0, P∈ρ^hb,e_ml =x^b_m< xi< x^e_m,yj=y^b,e_l ;

≥0, P∈∂τ_ml^h \

ρ^hb_ml∪ρ^he_ml.

(3.30)

By the maximum principle inLemma 2.1,

Vm⁽ⁿ⁺¹⁾(P)−V_l⁽ⁿ⁺¹⁾(P)≥0, P∈τ^h_ml. (3.31) By (3.8),V_l⁽ⁿ⁺¹⁾(P)=Vm⁽ⁿ⁺¹⁾(P),P∈ρ^hb_ml∪ {(x^e_m,y_l^b), (x^b_m,y^b_l)}, and from (3.10) and (3.31), it follows that

−μ²Ᏸ²xVm⁽ⁿ⁺¹⁾(P)= −μ²Ᏸ²xV⁽ⁿ⁺¹⁾(P), P∈ρ^hb_ml,

−μ²Ᏸ²yV_m⁽ⁿ⁺¹⁾(P)≤ −μ²Ᏸ²yV⁽ⁿ⁺¹⁾(P), P∈ρ^hb_ml. (3.32)

Thus, using (3.19), we conclude

G⁽ⁿ⁺¹⁾(P)≥ᏸ^hVm⁽ⁿ⁺¹⁾+fP,Vm⁽ⁿ⁺¹⁾

≥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 = ∪^Mm=⁻1¹(x^b,em, y_l^b). AtP^b_ml=(x^b_m,y_l^b), we have

V_ml⁽ⁿ⁺¹⁾P_ml^b =V_m⁽ⁿ⁺¹⁾P_ml^b =V_l⁽ⁿ⁺¹⁾P_ml^b , (3.34) and from (3.10), it follows that

−μ²Ᏸ²xV⁽ⁿ⁺¹⁾P_ml^b =− μ^{2 b}xm

Vm⁽ⁿ⁺¹⁾

P_ml^bx+−V_ml⁽ⁿ⁺¹⁾P_ml^b

h^b+_xm ⁻

V_ml⁽ⁿ⁺¹⁾P_ml^b −V_ml⁽ⁿ⁺¹⁾P_ml^bx− h^b_xm⁻

,

−μ²Ᏸ²yV⁽ⁿ⁺¹⁾P_ml^b =−μ^{2 b}_yl

V_l⁽ⁿ⁺¹⁾P^by+_ml −V_ml⁽ⁿ⁺¹⁾P_ml^b

h^b+_yl ⁻

V_ml⁽ⁿ⁺¹⁾P_ml^b −V_ml⁽ⁿ⁺¹⁾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⁻¹h^b_xm⁻+h^b+_xm, ^b_yl=2⁻¹h^b_yl⁻+h^b+_yl,

(3.35) whereh^b+_xm,h^b_xm⁻are the mesh step sizes on the left and right fromP_ml^b , andh^b+_yl,h^b_yl⁻are the mesh step sizes on the top and bottom fromP^b_ml. From here, (3.17), (3.23) and (3.27), we conclude

G⁽ⁿ⁺¹⁾(P)≥ᏸ^hV_ml⁽ⁿ⁺¹⁾+fP,V_ml⁽ⁿ⁺¹⁾≥0, P=P^b_ml. (3.36) With a similar argument for mesh pointP_l^e=(x^em,y_l^b), we prove thatV⁽ⁿ⁺¹⁾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⁽ⁿ⁾(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⁽ⁿ⁾=0, limZ_l⁽ⁿ⁾=0 and limZ_l⁽ⁿ⁾=0 asn→ ∞. Now by linearity of the operatorᏸ^h and the continuity of f, we have also from (3.5)–(3.8) that the mesh functionUdefined by

U(P)=_nlim

→∞V⁽ⁿ⁾(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 solutionsU1and U2to (2.8), then by the mean-value theorem, the differenceδU=U1−U2satisfies the following difference problem

ᏸ^hδU+fuδU=0, P∈ω^h, δU=0, P∈∂ω^h. (3.38)