A SEMI-DISCRETE TAILORED FINITE POINT METHOD FOR A ...
A SEMI-DISCRETE TAILORED FINITE POINT METHOD FOR A CLASS OF ANISOTROPIC DIFFUSION PROBLEMS HOUDE HAN, ZHONGYI HUANG† AND WENJUN YING Abstract. This work proposes a tailored finite point method (TFPM) for the numerical solution of an anisotropic diffusion problem, which has much smaller diffusion coefficient along one direction than the other on a rectangle domain. The paper includes analysis on the differentiability of the solution to the given problem under some compatibility conditions. It has detailed derivation for a semi-discrete TFPM for the given problem. This work also proves a uniform error estimate on the approximate solution. Numerical results show the TFPM is accurate as well as efficient for the strongly anisotropic diffusion problem. Examples include those that do not satisfy compatibility and regularity conditions. For the incompatible problems, numerical experiments indicate that the method proposed can still offer good numerical approximations.
1. Introduction In this paper, we consider the following anisotropic diffusion problem in two-dimensional space, −ε2
∂ 2 uε ∂ 2 uε − + a(x, y)uε = f (x, y), ∂x2 ∂y 2 uε ∂Ω = 0,
∀(x, y) ∈ Ω,
(1.1) (1.2)
with (cf. Fig. 1) { } Ω = (x, y) 0 < x < 1, 0 < y < 1 and its boundary ∂Ω = Γe ∪ Γw ∪ Γs ∪ Γn . Suppose that the functions a(x, y), f (x, y) satisfy ¯ a, f ∈ C (4,α) (Ω),
(1.3)
a(x, y) ≥ 0,
(1.4)
(α > 0), ¯ ∀(x, y) ∈ Ω,
Date: February 23, 2013. 1991 Mathematics Subject Classification. Primary: 65N35, 65N38; Secondary: 35L10. Key words and phrases. Tailored finite point method, anisotropic diffusion problem, boundary layer, compatibility conditions. † Corresponding author. H. Han was supported by the NSFC Project No. 10971116. Z. Huang was supported by the NSFC Project No. 11071139, the National Basic Research Program of China under the grant 2011CB309705. W. Ying was supported in part by the National Natural Science Foundation of China under Grant DMS– 11101278 and the Young Thousand Talents Program of China. 1
2
H.D. HAN, Z.Y. HUANG AND W.J. YING
p
4
n
p
3
w
e
p1
s
p2
Figure 1. The rectangle domain Ω = (0, 1) × (0, 1). ¯ with k be and 0 < ε ≤ 1. Here, C (k,α) (Ω) stands for the H¨older Space of index (k, α) on Ω a non-negative integer and α be a positive real number in the interval (0, 1]. Let {pl }4l=1 be ¯ (see Fig. 1). the four corners of the rectangle domain Ω When 0 < ε ≪ 1, problem (1.1)–(1.2) is an anisotropic diffusion problem and its solution may have boundary layers on a portion of ∂Ω, i.e. at x = 0 and x = 1. These layers are characterized by rapid transitions in the solution and difficult to capture in a numerical approximation without using a large number of unknowns. Also, such layers tend to cause spurious oscillations in a numerical solution to the problem. The anisotropic diffusion problem appears in a large variety of applications such as image processing [17, 22], chemical reactions [20] and anisotropic materials [15]. It has been numerically solved with asymptotic preserving method [2], least-square finite element method [16], mimetic finite difference methods [13], etc. A similar anisotropic but first-order ordinary differential system is solved with a finite difference method on Shishkin meshes [14]. This work proposes a tailored finite point method (TFPM) [6, 7, 8, 9, 10, 12] to numerically solve problem (1.1)–(1.2). The tailored finite point method was first proposed by Han, Huang and Kellogg [10] for the numerical solutions of singular perturbation problems of second-order elliptic equations with constant coefficients. Later, Han and Huang [6, 7, 8, 9] and Shih, Kellogg et al. [18, 19] systematically develop this method for the nonhomogeneous reactiondiffusion, convection-diffusion and convection-diffusion-reaction problems. The basic idea of the tailored finite point method is that the numerical scheme is tailor-made at each point/cell based on the local properties of the solution of the given problem. The method is expected to capture boundary layers associated with the problem (1.1)–(1.2) even with coarse grids. The remainder of this paper is organized as follows. Section 2 makes studies on the differentiability of the solution to problem (1.1)–(1.2) under some compatibility conditions. Section 3 derives a semi-discrete approximation of problem (1.1)–(1.2). The details of the TFPM for the semi-discrete and fully-discrete approximations of the problem are described in Section 4. Section 4 gives a proof for the uniform convergence property of the approximate solution. Finally, numerical examples supporting the theory and demonstrating the efficiency and reliability of the method are presented in Section 5.
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
3
2. The differentiability of the solution to problem (1.1)–(1.2) The differentiability of uε (x, y), the solution to problem (1.1)–(1.2), in the open domain Ω only depends on the differentiability of the given functions a(x, y) and f (x, y). By condition (1.3) and an interior a prior estimate for the elliptic equation (1.1), it is known that [11] uε ∈ C (6,α) (Ω).
(2.1)
¯ also But the differentiability of the solution to problem (1.1)–(1.2) on the closed domain Ω depends on some compatibility conditions at the corner points (refer to Grisvard [5] and Han & Kellogg [11]). For problem (1.1)–(1.2), suppose that the source term satisfies the following compatibility conditions } f (pl ) = 0, l = 1, · · · , 4, (2.2) fxx (pl ) = 0, fyy (pl ) = 0, l = 1, · · · , 4. ¯ Theorem 3.2 in Han & Kellogg [11] implies that uε (x, y) ∈ C (4,α) (Ω). Remark 2.1. If the compatibility conditions (2.2) (or part of them) do not hold, then the partial derivatives of the solution uε (x, y) may not be continuous at four corners [5, 11]. ¯ denote by ∥v∥∞ = max |v(x, y)| its For any function v(x, y) that is continuous on Ω, ¯ (x,y)∈Ω
maximum norm. Let Mk1 ,k2 Fk1 ,k2 Ak1 ,k2
k +k
∂ 1 2 uε
=
∂xk1 ∂y k2 ,
k +k ∞
∂ 1 2f
=
∂xk1 ∂y k2 ,
k +k ∞
∂ 1 2a
=
∂xk1 ∂y k2 , ∞
k1 ≥ 0, k2 ≥ 0, k1 + k2 ≤ 4, k1 ≥ 0, k2 ≥ 0, k1 + k2 ≤ 4, k1 ≥ 0, k2 ≥ 0, k1 + k2 ≤ 4,
and M0 = M0,0 , A0 = A0,0 , F0 = F0,0 . Assumption 2.1. In the rest of this section, assume that the constants Fk1 ,k2 and Ak1 ,k2 are bounded uniformly with respect to ε. Lemma 2.1. There exists a constant C independent of ε such that F0 . (2.3) M0 ≤ C ≡ 8 Proof. This result follows directly from the maximum principle of elliptic problem (1.1)–(1.2). Define the differential operator ∂2 ∂2 Lε ≡ −ε − +a ∂x2 ∂y 2 and introduce an auxiliary function F0 φ(y) = y(1 − y), 2 2
(2.4)
4
H.D. HAN, Z.Y. HUANG AND W.J. YING
which depends on variable y only. It is easy to check that 0 ≤ φ(y) ≤ C ≡
F0 , 8
∀y ∈ [0, 1],
and Lε φ(y) = F0 + a(x, y)φ(y) ≥ 0,
¯ ∀(x, y) ∈ Ω.
Let v± = uε ± φ. We have v+ |∂Ω ≥ 0, and
v− |∂Ω ≤ 0
( ) Lε v± (x, y) = f (x, y) ± Lε φ(y) = f (x, y) ± F0 + a(x, y)φ(y) .
By the definition of F0 and the assumption on the coefficient a(x, y), we get Lε v+ (x, y) ≥ 0,
Lε v− (x, y) ≤ 0,
∀(x, y) ∈ Ω.
The maximum principle further implies v+ (x, y) ≥ 0,
v− (x, y) ≤ 0,
¯ ∀(x, y) ∈ Ω.
As a result −φ(y) ≤ uε (x, y) ≤ φ(y),
¯ ∀(x, y) ∈ Ω. □
Then inequality (2.3) follows immediately. Lemma 2.2. The following inequalities hold 2 δ 1 M0,j−1 + M0,j+1 , ∀δ ∈ (0, ], j = 1, 2, 3; δ 2 2 2 j−1 δ j+1 1 ≤ ε Mj−1,0 + ε Mj+1,0 , ∀δ ∈ (0, ], j = 1, 2, 3. δ 2 2
M0,j ≤ εj Mj,0
(2.5) (2.6)
Proof. At first, we prove estimate (2.5). By the definition of M0,1 , there exists a point ¯ such that (x0 , y0 ) ∈ Ω, ∂uε (x0 , y0 ) . M0,1 = ∂y If 0 ≤ y0 ≤
1 2
(when
1 2
< y0 ≤ 1, the proof is analogous), for any δ ∈ (0, 12 ], we have
δ 2 ∂ 2 uε ∂uε (x0 , y0 ) + (x0 , y0 + ξ), uε (x0 , y0 + δ) = uε (x0 , y0 ) + δ ∂y 2 ∂y 2
with ξ ∈ [0, δ].
(2.7)
From (2.7), estimate (2.5) for j = 1 is obtained. For the case j = 2, 3, the proof is similar. Furthermore, following the same lines as those for (2.5), we can prove estimate (2.6). □ Let ∂ j uε (x, y) , ∂y j ∂ j uε (x, y) wj (x, y) = , ∂xj vj (x, y) =
j = 1, 2, 3, 4; j = 1, 2, 3, 4.
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
5
It can be verified that function v2 (x, y) satisfies the following problem ∂ 2 v2 ∂ 2 v2 −ε2 2 − + av2 = fyy (x, y) − ayy (x, y)uε − 2ay (x, y)v1 , ∂x ∂y 2 v2 Γe = v2 Γw = 0, v2 Γn = −f (x, 1), v2 Γs = −f (x, 0).
∀(x, y) ∈ Ω,
(2.8) (2.9) (2.10)
The boundary conditions (2.10) are from equation (1.1). Similar to the proof for Lemma 2.1, introducing an auxiliary function ψ(y) ≡ F0 + (F0,2 + A0,2 M0 + 2A0,1 M0,1 ) we have 0 ≤ ψ(y) ≤ F0 +
y(1 − y) , 2
F0,2 + A0,2 M0 + 2A0,1 M0,1 , 8
∀y ∈ [0, 1],
∀y ∈ [0, 1]
(2.11)
and Lε ψ(y) = F0,2 + A0,2 M0 + 2A0,1 M0,1 + a(x, y)ψ(y). Then it is easy to check that Lε (v2 + ψ) ≥ 0, and Lε (v2 − ψ) ≤ 0. The maximum principle implies M0,2 = max |v2 (x, y)| ≤ max |ψ(y)|. (x,y)∈Ω
y∈[0,1]
By (2.11) and (2.5) in Lemma 2.2, we obtain F0,2 + A0,2 M0 + 2A0,1 M0,1 8 F0,2 + A0,2 M0 A0,1 δA0,1 ≤ F0 + + M0 + M0,2 , 8 2δ 8
M0,2 ≤ F0 +
{
Let δ = min we arrive at
As we have (let δ =
1 2
F0,2 M0,2 ≤ 2F0 + + 4 in equation (2.5))
(
1 4 , 2 A0,1
(2.12)
} ,
A0,2 A0,1 + 4 δ
M0,1 ≤ 4M0 +
1 ∀δ ∈ (0, ]. 2
(2.13) ) M0 ≡ C2 .
(2.14)
M0,2 , 4
we get the following Lemma. Lemma 2.3. There exists a constant C independent of ε such that M0,1 ≤ C,
M0,2 ≤ C.
Furthermore, from equation (1.1) and inequality (2.3), we have
(2.15)
6
H.D. HAN, Z.Y. HUANG AND W.J. YING
Lemma 2.4. There exists a constant C independent of ε such that εj Mj,0 ≤ C,
j = 1, 2.
(2.16)
Now we try to get a prior estimate of v4 . First, we discuss the boundary condition of v4 at Γs and Γn . From equation (2.8), we have v4 =
4 ∂ 2 v2 2 ∂ uε = −ε + a(x, y)v2 − fyy (x, y) + ayy (x, y)uε + 2ay (x, y)v1 . ∂y 2 ∂y 2 ∂x2
(2.17)
Differentiating equation (1.1) with respect to x twice, we obtain ∂ 4 uε = −ε2 w4 + a(x, y)w2 + axx (x, y)uε + 2ax (x, y)w1 − fxx (x, y). ∂y 2 ∂x2
(2.18)
Substituting (2.18) into (2.17) and using boundary conditions wj |Γn = wj |Γs = 0 (j = 1, 2, 3, 4) yield ( ) v4 Γn = ε2 fxx + av2 − fyy + 2ay v1 Γn ≡ gn4 (x), (2.19) ( 2 ) v4 Γs = ε fxx + av2 − fyy + 2ay v1 Γn ≡ gs4 (x). (2.20) This means, by Lemma 2.3, we can find a constant C independent of ε such that ∥gn4 ∥Γn ,∞ ≤ C,
∥gs4 ∥Γs ,∞ ≤ C.
(2.21)
The function v4 satisfies ∂ 2 v4 ∂ 2 v4 −ε2 2 − + av4 = fyyyy − ayyyy uε − 4ayyy v1 − 6ayy v2 − 4ay v3 , ∀(x, y) ∈ Ω, (2.22) ∂x ∂y 2 v2 Γe = v2 Γw = 0, (2.23) v2 Γn = gn4 (x), v2 Γs = gs4 (x). (2.24) As before, applying the maximum principle to problem (2.22)–(2.24), we get Lemma 2.5. There exists a constant C independent of ε such that
4
∂ uε
∂y 4 ≤ C.
(2.25)
∞
Combining Lemma 2.1 – Lemma 2.5, we obtain Theorem 2.1. There exists a constant C independent of ε such that
j
∂ uε
∂y j ≤ C, j = 0, 1, · · · , 4.
(2.26)
∞
Similarly, we have Theorem 2.2. There exists a constant C independent of ε such that
j
j ∂ uε
ε
∂xj ≤ C, j = 0, 1, · · · , 4. ∞
Finally, we get the theorem
(2.27)
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
7
Theorem 2.3. There exists a constant C independent of ε such that
j ∂ j1 +j2 uε
ε 1
∂xj1 ∂y j2 ≤ C, 0 ≤ j1 + j2 ≤ 4.
(2.28)
∞
3. A semi-discrete approximation of problem (1.1)–(1.2) The solution of problem (1.1)–(1.2) is smooth in y-direction. It is expected the gradient of the solution in y-direction will in general have small variations provided that the coefficients aj (x) and the source term fj (x) are nice functions. We can discretize the problem in ydirection by the standard finite difference method, which allows us to transform the secondorder PDE into an ODE system. It is much like the method of lines. Taking a positive integer M , let ∆y = M1 and yj = j∆y,
j = 0, 1, · · · , M.
We discretize the second-order partial derivative with respect to the y variable in problem (1.1)-(1.2) with the centered three-point finite difference. This gives us the following semidiscrete system −ε2
d2 uεj (x) uεj−1 (x) − 2uεj (x) + uεj+1 (x) − + aj (x)uεj (x) = fj (x), 1 ≤ j ≤ M − 1, (3.1) dx2 (∆y)2 uεj (0) = uεj (1) = 0, 1 ≤ j ≤ M − 1, (3.2) uε0 (x) = uεM (x) = 0,
∀x ∈ [0, 1],
(3.3)
with aj (x) = a(x, yj ),
fj (x) = f (x, yj ),
j = 1, · · · , M − 1.
Here, uεj (x) is the finite difference approximation of uε (x, yj ) for each j ∈ {0, 1, · · · , M }. Let Ujε (x) = uε (x, yj ),
j = 0, 1, · · · , M
(3.4)
be the exact solution at the discrete points {yj }M The differentiability of functions j=0 . { ε } Uj (x), j = 0, 1, · · · , M that we derived in Section 2 allows us to make local Taylor expansions for them. As a result, the functions satisfy the modified semi-discrete system −ε2
ε ε (x) − 2Ujε (x) + Uj+1 d2 Ujε (x) Uj−1 (x) − + aj (x)Ujε (x) = fj (x) + Rjε (x), (3.5) dx2 (∆y)2 1 ≤ j ≤ M − 1,
Ujε (0) = Ujε (1) = 0, ε (x) = 0, U0ε (x) = UM
1 ≤ j ≤ M − 1, (3.6) ∀x ∈ [0, 1],
(3.7)
with the remainder term Rjε (x) on the order of ∆y 2 in that there is a constant C independent of ε such that (3.8) |Rjε (x)| ≤ C(∆y)2 , 1 ≤ j ≤ M − 1. The estimate (3.8) for the remainder term results from the application of Theorems 2.1–2.3.
8
H.D. HAN, Z.Y. HUANG AND W.J. YING
Lemma 3.1. Let ( )T ε vε (x) = v1ε (x), v2ε (x), · · · , vM −1 (x) be a vector-valued function. If for any 0 < x < 1 and 1 ≤ j ≤ M − 1, ε ε d2 vjε (x) vj−1 (x) − 2vjε (x) + vj+1 (x) −ε − + aj (x)vjε (x) ≤ 0 2 2 dx (∆y) 2
then we have { vjε (x)
≤ max
0,
max
{
1≤j≤M −1
vjε (0), vjε (1)
}
{ } ε , max v0ε (x), vM (x)
}
0≤x≤1
.
(3.9)
Proof. First, we prove that if ε ε d2 vjε (x) vj−1 (x) − 2vjε (x) + vj+1 (x) −ε − + aj (x)vjε (x) < 0, 2 2 dx (∆y) 2
(3.10)
then vjε (x) can not attain its maximum non-negative value in the interior of interval 0 < x < 1 and for 1 ≤ j ≤ M − 1. Assume that vjε (x) attains its maximum non-negative value at point x = x∗ for some 1 ≤ j = j ∗ ≤ M − 1. This indicates ε2
d2 vjε∗ (x∗ ) ≤ 0 and dx2
vjε∗ −1 (x∗ ) − 2vjε∗ (x∗ ) + vjε∗ +1 (x∗ ) ≤ 0. (∆y)2
The condition (3.10) further implies 0 ≤ aj ∗ (x∗ )vjε∗ (x∗ ) < ε2
d2 vjε∗ (x∗ ) vjε∗ −1 (x∗ ) − 2vjε∗ (x∗ ) + vjε∗ +1 (x∗ ) + ≤ 0. dx2 (∆y)2
This is a contradiction. As a result, the inequality (3.9) holds. In the next, we prove the lemma for the case ε ε d2 vjε (x) vj−1 (x) − 2vjε (x) + vj+1 (x) −ε − + aj (x)vjε (x) ≤ 0. dx2 (∆y)2 2
Let 1 1 ξj = − 8 2
( )2 1 yj − 2
(3.11)
and wjη (x) = vjε (x) − η ξj
with η > 0 be a small positive parameter. It obvious that ξj ≥ 0 and
ξj−1 − 2ξj + ξj+1 = −1. (∆y)2
(3.12)
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
9
We see that wjη (x) satisfies the condition (3.10) since η η (x) d2 wjη (x) wj−1 (x) − 2wjη (x) + wj+1 −ε − + aj (x)wjη (x) 2 2 dx (∆y) 2 ε ε ε d vj (x) vj−1 (x) − 2vjε (x) + vj+1 (x) ξj−1 − 2ξj + ξj+1 2 = −ε − + aj (x)vjε (x) + η − η aj ξj 2 2 dx (∆y) (∆y)2 ε ε d2 vjε (x) vj−1 (x) − 2vjε (x) + vj+1 (x) ξj−1 − 2ξj + ξj+1 ε ≤ −ε2 − + a (x)v (x) + η j j dx2 (∆y)2 (∆y)2 ε ε d2 vjε (x) vj−1 (x) − 2vjε (x) + vj+1 (x) = −ε2 − + aj (x)vjε (x) − η ≤ −η < 0. 2 2 dx (∆y) 2
Here, we used the condition (3.11). The previous conclusion implies { } { η } η η η η wj (x) ≤ max 0, max wj (0), wj (1) , max {w0 (x), wM (x)} . 1≤j≤M −1
0≤x≤1
(3.13)
Since the parameter η can be arbitrarily small, letting η → 0, we get the inequality (3.9) for vjε (x). □ Now we are ready to prove the following theorem. Theorem 3.1. Let )T uε1 (x), · · · , uεM −1 (x) , ( )T ε Uε (x) = U1ε (x), · · · , UM . −1 (x)
U ε (x) =
(
The following error estimate holds ∥U ε − Uε ∥∞ ≤ C(∆y)2 ,
(3.14)
where C is a constant independent of ε. Proof. Subtracting the semi-discrete system (3.1)-(3.3) by the semi-discrete system (3.5)(3.7) yields a new semi-discrete system −ε2
d2 eεj (x) eεj−1 (x) − 2eεj (x) + eεj+1 (x) − + aj (x)eεj (x) = Rjε (x), 1 ≤ j ≤ M − 1, dx2 (∆y)2 eεj (0) = eεj (1) = 0, 1 ≤ j ≤ M − 1, eε0 (x) = eεM (x) = 0,
∀x ∈ [0, 1],
for the vector-valued error function )T ( eε (x) = eε1 (x), eε2 (x), · · · , eεM −1 (x) = Uε − U ε with eεj (x) = Ujε (x) − uεj (x). Let
[
1 1 φj = − 8 2
( )2 ]
ε 1
Rj (x) yj − ∞ 2
10
H.D. HAN, Z.Y. HUANG AND W.J. YING
with Rjε (x) ∞ = −ε2
max
1≤j≤M −1, 0≤x≤1
|Rjε (x)|. Let vjε (x) = eεj (x) − φj . We have
ε ε
ε (x) − 2vjε (x) + vj+1 (x) d2 vjε (x) vj−1 ε ε
Rj (x) ≤ 0. − + a (x)v (x) ≤ R (x) − j j j ∞ dx2 (∆y)2
Lemma 3.1 implies vjε (x) ≤ 0 for all j ′ s and x ∈ [0, 1]. Thus we obtain eεj (x) ≤ φj ≤ max φj = 0≤j≤M
1
Rjε (x) . ∞ 8
(3.15)
Let wjε (x) = −eεj (x) − φj . We get ε ε
ε (x) − 2wjε (x) + wj+1 (x) d2 wjε (x) wj−1 ε ε
Rj (x) ≤ 0. − + a (x)w (x) ≤ −R (x) − −ε j j j ∞ dx2 (∆y)2 2
Similarly, Lemma 3.1 implies wjε (x) ≤ 0 for all j ′ s and x ∈ [0, 1]. We arrive at −eεj (x) ≤ φj ≤ max φj = 0≤j≤M
1
Rjε (x) . ∞ 8
Combining (3.15) and (3.16), we get
ε
ej (x) ≤ 1 Rjε (x) . ∞ ∞ 8 Together with (3.8), we finally obtain the error estimate (3.14).
(3.16)
□
4. A tailored finite point method for problem (3.1)–(3.3) The problem (3.1)–(3.3) can be rewritten as −ε2
d2 U ε + A(x) U ε = F (x), dx2 U ε (0) = U ε (1) = 0,
∀x ∈ (0, 1),
(4.1) (4.2)
where A(x) is an (M − 1) × (M − 1) matrix function, F (x) is an (M − 1) dimensional vector function, 2 1 − 0 (∆y)2 + a1 (x) (∆y)2 1 2 .. . − + a1 (x) 2 2 (∆y) (∆y) A(x) = 1 .. .. . . − 2 (∆y) 1 2 0 − + a (x) M −1 (∆y)2 (∆y)2 f1 (x) f2 (x) . F (x) = .. . fM −1 (x)
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
11
Taking a positive integer N , let ∆x =
1 N
and
xi = i∆x,
i = 0, 1, · · · , N.
We will introduce piecewise constant matrix functions Ah (x) and Fh (x) for 0 ≤ x ≤ 1, which are approximations of A(x) and F (x) respectively. On each interval [xi , xi+1 ], taking the midpoint x∗i ∈ [xi , xi+1 ] (i.e. , x∗i = xi +x2 i+1 ), define the matrix functions Ah (x) and Fh (x) by Ah (x) = A(x∗i ),
for x ∈ (xi , xi+1 ),
Fh (x) = F (x∗i ),
for x ∈ (xi , xi+1 ).
This means the matrix functions Ah (x) and Fh (x) are constants on each subinterval (xi , xi+1 ) (i = 0, 1, · · · , N − 1) and have jump discontinuities at points x = xi (i = 1, 2, · · · , N − 1). Now we consider the following approximation of problem (4.1)–(4.2), d2 Uhε + Ah (x) Uhε = Fh (x), 2 dx Uhε (0) = Uhε (1) = 0, [ ε] dUh ε [Uh ] x=xi = = 0, dx x=xi −ε2
∀x ∈ (xi , xi+1 ),
i = 0, · · · , N − 1
(4.3) (4.4) (4.5)
with the jumps defined by [Uhε ] x=xi = Uhε (xi + 0) − Uhε (xi − 0), [ ε] dUh dUhε dU ε = (xi + 0) − h (xi − 0). dx x=xi dx dx
(4.6) (4.7)
Here, the semi-discrete solution ( )T ε,h Uhε (x) = uε,h 1 (x), · · · , uM −1 (x) is an approximation of U ε (x). Correspondingly, each component uε,h j (x) of the vector-valued ε function is also an approximation of uj (x). Lemma 4.1. The problem (4.3)–(4.5) has a unique solution Uhε (x) and the solution satisfies the following estimate ∥Uhε ∥∞ ≤ C,
(4.8)
where C is a constant independent of ε. Lemma 4.2. The following uniform error estimate holds ∥U ε − Uhε ∥∞ ≤ C∆x, where C is a constant independent of ε and ∆x.
(4.9)
12
H.D. HAN, Z.Y. HUANG AND W.J. YING
For x ∈ [0, 1], we define yj+1 − y ε,h y − yj ε,h uε,h (x, y) = uj (x) + u (x), for y ∈ [yj , yj+1 ], ∆y ∆y j+1
(j = 0, · · · , M − 1)
where ε,h uε,h 0 (x) = uM (x) ≡ 0.
Combining with Theorem 3.1, we get the following error estimate: Theorem 4.1. Suppose that uε is the solution of problem (1.1)–(1.2), uε,h is defined above. The following error estimate holds { } ∥uε,h − uε ∥∞ ≤ C ∆x + (∆y)2 , (4.10) where C is a constant independent of ε, ∆x and ∆y. Now we discuss how to get the solution Uhε (x) exactly. We will propose a tailored finite point method to solve it [6, 7, 8, 9, 10, 12]. On each subinterval [xi , xi+1 ], Uhε (x) satisfies −ε2
d2 Uhε (x) + Ai Uhε (x) = Fi , dx2
for x ∈ (xi , xi+1 ),
(4.11)
with Ai = A(x∗i ),
Fi = F (x∗i ).
Let Uhf,i (x) = A−1 i Fi ,
(4.12)
which is a particular solution to the nonhomogeneous equation (4.11). Let Uhi (x) = ξi e
λx ε
∈ RM −1 ,
(4.13)
be a general solution of the homogeneous equation 2 ε 2 d Uh (x) −ε dx2
+ Ai Uhε (x) = 0,
for x ∈ (xi , xi+1 ).
(4.14)
Substituting (4.13) into (4.14), we arrive at A i ξ i = λ2 ξ i .
(4.15)
Since Ai is a symmetric definite real matrix, it has M − 1 positive real eigenvalues 0 < µi1 ≤ µi2 ≤ · · · ≤ µiM −1 , and corresponding eigenvectors Let λik =
√
ξ i1 , ξi2 , · · · , ξ iM −1 .
µik . The general solution of equation (4.11) on the interval [xi , xi+1 ] is given by ( i ( i ) )) M −1 ( ∑ λk (x − xi+1 ) −λk (x − xi ) f,i + − ε Uh (x) = Uh (x) + αk exp + αk exp ξ ik , (4.16) ε ε k=1
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
13
with some constants αk± (k = 1, · · · , M − 1). Let ( ) Ki = ξ i1 , · · · , ξ iM −1 , { ( i ) ( i )} λM −1 (x − xi+1 ) λ1 (x − xi+1 ) i Λ+ (x) = diag exp , · · · , exp , ε ε { ( i ) ( i )} −λM −1 (x − xi ) −λ1 (x − xi ) i Λ− (x) = diag exp , · · · , exp , ε ε ( )T + αi+ = α1+ , α2+ , · · · , αM , −1 ( ) T − αi− = α1− , α2− , · · · , αM . −1 Then equality (4.16) can be rewritten by ∀x ∈ [xi , xi+1 ].
Uhε (x) = Uhf,i (x) + Ki Λi+ (x)αi+ + Ki Λi− (x)αi− ,
(4.17)
At points x = xi , xi+1 , we obtain Uhε (xi ) = Uhf,i (xi ) + Ki Λi+ (xi )αi+ + Ki Λi− (xi )αi− ,
(4.18)
Uhε (xi+1 ) = Uhf,i (xi+1 ) + Ki Λi+ (xi+1 )αi+ + Ki Λi− (xi+1 )αi− . Let
( Ai =
Ki Λi+ (xi ) Ki Λi− (xi ) Ki Λi+ (xi+1 ) Ki Λi− (xi+1 )
(4.19)
) .
Then the coefficients αi± are given by ( ) ( ) f,i i ε ( ) α+ Uh (xi ) − Uh (xi ) −1 = Ai . i ε α− Uh (xi+1 ) − Uhf,i (xi+1 ) Assume that ( We have
(
αi+ αi−
)
( =
)−1 Ai ≡
(
Ai11 Ai21
Ai12 Ai22
) .
Ai11 Uhε (xi ) + Ai12 Uhε (xi+1 ) − (Ai11 + Ai12 ) Uhf,i Ai21 Uhε (xi ) + Ai22 Uhε (xi+1 ) − (Ai21 + Ai22 ) Uhf,i
) .
(4.20)
From (4.16), we get dUhε (xi + 0) d d = Ki Λi+ (xi )αi+ + Ki Λi− (xi )αi− , dx dx dx dUhε (xi − 0) d i−1 d i−1 i−1 . = Ki−1 Λi−1 Λi−1 (xi )α− + (xi )α+ + K dx dx dx − Continuity condition (4.7) on the first-order derivative implies Ki
(4.21) (4.22)
d d d i i−1 d i−1 i−1 . (4.23) Λ+ (xi )αi+ + Ki Λi− (xi )αi− = Ki−1 Λi−1 Λi−1 (xi )α− + (xi )α+ + K dx dx dx dx −
14
H.D. HAN, Z.Y. HUANG AND W.J. YING
Substituting (4.20) into (4.23), we finally get a numerical scheme Ai Uhε (xi−1 ) + Bi Uhε (xi ) + Ci Uhε (xi+1 ) = Fi ,
i = 1, 2, · · · , N − 1,
Uhε (0) = Uhε (1) = 0,
(4.24) (4.25)
where d i−1 i−1 d Λ+ (xi )Ai−1 Λi−1 (xi )Ai−1 11 − K 21 , dx dx − d d i−1 d i−1 d i−1 i−1 − Ki−1 Λ− , (xi )A12 (xi )A22 = Ki Λi+ (xi )Ai11 + Ki Λi− (xi )Ai21 − Ki−1 Λ+ dx dx dx dx d d = Ki Λi+ (xi )Ai12 + Ki Λi− (xi )Ai22 , dx dx ( ) ( ) ( i ) d d i i i i i i = K Λ (xi ) A11 + A12 + Λ− (xi ) A21 + A22 Uhf,i dx + dx ) ( ) ) ( i−1 ( i−1 d i−1 d i−1 i−1 i−1 i−1 Uhf,i−1 . −K Λ+ (xi ) A11 + A12 + Λ− (xi ) A21 + A22 dx dx
Ai = −Ki−1 Bi Ci Fi
5. Numerical examples This section presents numerical examples to demonstrate the efficiency and reliability of the tailored finite point method. In the examples, the error eh (x, y) of a numerical solution is computed in the discrete ℓ2 -norm defined by ∥eh ∥2l2
( ) 2 N ∑ M ∑ 1 1 eh xi + ∆x, yj + ∆y . = ∆x∆y 2 2 i=1 j=1
Example 5.1. First, let a(x, y) = 0, f (x, y) = 100 sin 2πy and ) ( 2π(x−1) 2πx 25 e− ε + e ε sin 2πy. uε (x, y) = 2 1 − π ) 1 − exp(− 2π ε
(5.1)
The function uε given in (5.1) is the exact solution to problem (1.1)–(1.2) with a and f given above. In this example, functions a and f satisfy the conditions given in Assumption 2.1. The numerical results are shown in Figures 2, 3 and Table 1. Second-order convergence rates for the numerical solutions are observed in the discrete l2 norm. Table 1. l2 errors of the numerical solutions for Example 5.1. ∆x and ∆y ε = 0.1 ∥uε − uε,h ∥l2 ε = 0.001 ∥uε − uε,h ∥l2
1 8
1 16
1 32
1 64
8.88E-2 2.23E-2 5.56E-3 1.39E-3 8.88E-2 2.25E-2 5.67E-3 1.43E-3
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
15 −3
x 10 8
2.5 2 3
1 0
6
1
0.008
0.5
0.006
5
0
0.004
4
0.002
3
−1 −0.5 −2 −1
−3 1
7
0.01
1.5
2
0 1
−1.5
2
1
1 −2
0.5
0.5
0.5 0
−2.5
0
0
1
0.5 0
0
(a)
(b) −3
x 10 2
2.5 −3
x 10
2 3 2
1
1
0.5
0
0
1.6
2
1.4
1.5
1.2 1
1
−1 −0.5
0.8
0.5
−2 −1
−3 1
1.8
2.5
1.5
0.6
0 1
−1.5 1
1 −2
0.5
0.5
0.5 0
−2.5
0
0.5 0
(c)
0.4 0.2 0
0
(d)
1 Figure 2. Example 5.1 with ε = 0.1: (a) uε,h with ∆x = ∆y = 32 ; (b) 1 1 |uε − uε,h | with ∆x = ∆y = 32 ; (c) uε,h with ∆x = ∆y = 64 ; (d) |uε − uε,h | 1 with ∆x = ∆y = 64 .
−1
10
−2
10
2 =0.1 =0.001
1 −3
10
−2
10
−1
10
0
10
Figure 3. Convergence rates of ∥uε − uε,h ∥l2 for Example 5.1.
Example 5.2. Now let a(x, y) = 0, f (x, y) = 1 and ) ( ) ( kπx kπ(x−1) ∞ y(1 − y) ∑ (−1)k − 1 e− ε + e ε uε (x, y) = + 2 sin kπy. 2 (kπ)3 1 + exp(− kπ ) ε k=1 The function uε given in (5.2) is the exact solution to problem (1.1)–(1.2).
(5.2)
16
H.D. HAN, Z.Y. HUANG AND W.J. YING
In this example, function f does not satisfy all the compatibility conditions (2.2). The numerical results are shown in Figure 4 (ε = 0.1) and Figure 5 (ε = 0.001) with mesh size ∆x = ∆y = 0.05. The method still yields good approximations on coarse grids (h ≫ ε). −5
x 10 9
0.12
−4
x 10 0.1
8
1 7 0.8
0.1
0.08
6 0.6 5
0.06
0.05
0.4 4 0.2 3
0.04 0 1
0 1 1 0.5
2
0.02
1
0
1
0.5
0.5 0
0
0.5 0
0
0
Figure 4. uε,h (left) and |uε − uε,h |(right), ε = 0.1. −5
x 10 9
0.12
−4
x 10 0.1
8
1 7 0.8
0.1
0.08
6 0.6 5
0.06
0.05
0.4 4 0.2 3
0.04 0 1
0 1 1 0.5 0
1 1
0.5
0.5 0
2
0.02
0
0.5 0
0
0
Figure 5. uε,h (left) and |uε − uε,h |(right) for example 5.2, ε = 0.001. ( )0.1 Example 5.3. Let a(x, y) = 1 and f (x, y) = 64x(1 − x)y(1 − y) . In this example, function f also does not satisfy all the compatibility conditions (2.2). The numerical results are shown in Figure 6 (ε = 0.1) and Figure 7 (ε = 0.001). The method also yields good numerical solutions when h ≫ ε. Example 5.4. Let a(x, y) = xy,
( ( )) ( x) x−1 + exp uε (x, y) = 100x(1 − x)y(1 − y) exp − , ε ε
and f (x, y) = Lε uε . In this example, the function f does not satisfy the conditions given in Assumption 2.1. The numerical results are shown in Figures 8, 9 and Table 2. The second-order convergence rates for the numerical solutions are also observed in the discrete l2 norm.
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
0.1
0.12
0.12
0.1
0.1 0.1
0.08
0.06
0.05
17
0.08
0.06
0.05
0.04
0.04
0 1
0 1 1 0.5
0.02
1 0.5
0.5 0
0
0
0
0
0
Figure 6. Graphs of uε,h for example 5.3, ε = 0.1, ∆y = 1 and ∆x = 200 (right).
1 : 20
∆x =
0.12
(left)
0.1 0.1
0.08
0.06
0.05
1 20
0.12
0.1 0.1
0.02
0.5
0.08
0.06
0.05
0.04
0.04
0 1
0 1 1 0.5
0.02
1 0.5
0.5 0
0
0
0.02
0.5 0
0
Figure 7. Graphs of uε,h for example 5.3, ε = 0.001, ∆y = 1 and ∆x = 200 (right).
1 : 20
∆x =
1 20
(left)
6. Conclusion This work proposes a tailored finite point method (TFPM) for the strongly anisotropic diffusion problem (1.1)–(1.2), which has much slower diffusion rate along one direction than the other in a rectangle domain. It presents analysis on the differentiability of the solution to problem (1.1)–(1.2) under some compatibility conditions and proves a uniform error estimate for the approximate solution to the discrete problem. Numerical experiments show that, for the problem that satisfies the compatibility conditions, the TFPM yields second-order convergence rates in the discrete l2 norm. For the problem that does not satisfy the compatibility and regularity conditions, the method proposed may still offer good approximate solutions.
Table 2. l2 errors of the numerical solutions for Example 5.4 with ∆y = ∆x
1 16
1 32
1 64
1 128
ε = 0.1 ∥uε − uε,h ∥l2 5.98E-3 1.45E-3 3.48E-4 8.58E-5 ε = 0.01 ∥uε − uε,h ∥l2 8.08E-3 2.12E-3 5.30E-4 1.27E-4
1 . 16
18
H.D. HAN, Z.Y. HUANG AND W.J. YING
In this work, for simplicity, the TFPM is proposed only for the problem defined in a rectangle domain. However, the method is by no means limited to the simple domain. Instead, the method can be generalized for problems on complex domains together with a structured grid method such as the kernel-free boundary integral method[21]. This extension will make the TFPM applicable to a larger class of problems. Another potential extension of the TFPM is its application for those problems where an advection or convection term is added into the PDE.
−3
x 10 0.02 4.5
0.018
−3
x 10
0.016
0.025 0.02
4
5 3.5
0.014
4
3
0.012 0.015
3
0.01
0.01 0.005 1
0 0
1
0.006
1
0.004
0 0
2 1.5 1
0.002
0.5
0.5
2.5
2
0.008
0.5
0
0
1
(a)
1 0.5
0.5
0
0
(b) −4
−3
x 10
x 10 1.2
3
−3
−4
x 10
x 10
1
1.2
3
0.8
2
1 0.8
2
0.6
0.6 0.4
1
0.4
1
0.2 1
0 0
1
0 0
0.2 0.5
0.5 1
0.5
0.5 0
0
1
(c)
0
0
(d)
Figure 8. Graphs of |uε − uε,h | for example 5.4, ∆y = 1 1 1 1 ∆x = 16 ; (b) ∆x = 32 ; (c) ∆x = 64 ; (d) ∆x = 128 .
1 , 16
ε = 0.1: (a)
−2
10
=0.1 =0.01
−3
10
2 −4
1
10
−3
10
−2
10
−1
10
Figure 9. Convergence rates of ∥uε − uε,h ∥l2 for Example 5.4.
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
19
References 1. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, Dover Publications, New York, 1964. 2. P. Degond, F. Deluzet, and C. Negulescu, An Asymptotic Preserving Scheme For Strongly Anisotropic Elliptic Problems, Multiscale Model. Simul. 8 (2010): 645–666. 3. L.C. Evans, Partial Differential Equations, American Mathematical Society, Providence, RI, 1998. 4. I.S. Gradshteyn and I.M. Ryzhik, Table of integrals, series and products, 6th Ed., Academic Press, San Diego, 2000. 5. P. Grisvard, Elliptic problems in nonsmooth domains, Pitman, Boston, 1985. 6. H. Han and Z. Huang, A tailored finite point method for the Helmholtz equation with high wave numbers in heterogeneous medium, J. Comp. Math, 26 (2008): 728-739. 7. H. Han and Z. Huang, Tailored finite point method for a singular perturbation problem with variable coefficients in two dimensions, J. Sci. Comp., 41 (2009): 200-220. 8. H. Han and Z. Huang, Tailored finite point method for steady-state reaction-diffusion equation, Comm. Math. Sci., 8 (2010): 887-899. 9. H. Han and Z. Huang, Tailored finite point method based on exponential bases for convection-diffusionreaction equation, Math. Comput., (to appear). 10. H. Han, Z. Huang and B. Kellogg, A Tailored finite point method for a singular perturbation problem on an unbounded domain, J. Sci. Comp. 36 (2008): 243-261. 11. H. Han and R. B. Kellogg, Differentiability properties of solutions of the equation −ε∆u + ru = f (x, y) in a square, SIAM J. Math. Anal. 21 (1990): 394-408. 12. Z. Huang, Tailored finite point method for the interface problem, Netw. Heterog. Media (2009) 4: 91-106. 13. J. Hyman, J. Morel, M. Shashkov and S. Steinberg, Mimetic finite difference methods for diffusion equations, Computational Geosciences 6 (2002): 333–352. 14. P. Maragatha Meenakshi, S. Valarmathi and J.J.H. Miller, Solving a partially singularly perturbed initial value problem on Shishkin meshes, Appl. Math. Comput., 215 (2010): 3170–3180. 15. J.G. Mullen, Effect of Bardeen-Herring Correlation on Vacancy Diffusion in Anisotropic Crystals, Phys. Rev., 124 (1961 ): 1723–1730. 16. J. Pasdunkorale and I.W. Turner, A second order control-volume finite-element least-squares strategy for simulating diffusion in strongly anisotropic media J comput. math. 23 (2005): 1–16. 17. P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Trans. Pattern Anal. Mach. Intel. 12 (1990): 629–639. 18. Y. Shih, R. B. Kellogg and P. Tsai, A Tailored Finite Point Method for Convection-Diffusion-Reaction Problems, Journal of Scientific Computing, 43 (2010): 239–260. 19. Y. Shih, R. B. Kellogg and Y. Chang, Characteristic Tailored Finite Point Method for ConvectionDominated Convection-Diffusion-Reaction Problems, Journal of Scientific Computing, 47 (2011): 198– 215. 20. G. Turk, Generating Textures on Arbitrary Surfaces Using Reaction-Diffusion, Computer Graphics, 25 (1991): 289–298. 21. W. Ying, C. S. Henriquez, A kernel-free boundary integral method for elliptic boundary value problems, J. Comp. Phy., 227 (2007): 1046-1074. 22. J. Weickert, Anisotropic Diffusion in Image Processing, Teubner, Stuttgart, 1998.
20
H.D. HAN, Z.Y. HUANG AND W.J. YING
(H.D. Han, Z.Y. Huang) Dept. of Mathematical Sciences, Tsinghua University, Beijing 100084, China E-mail address: [email protected] E-mail address: [email protected] (W.J. Ying) Dept. of Mathematics & Institute of Natural Sciences, Shanghai Jiao Tong University, Shanghai 200240, China E-mail address: [email protected]
1. Introduction In this paper, we consider the following anisotropic diffusion problem in two-dimensional space, −ε2
∂ 2 uε ∂ 2 uε − + a(x, y)uε = f (x, y), ∂x2 ∂y 2 uε ∂Ω = 0,
∀(x, y) ∈ Ω,
(1.1) (1.2)
with (cf. Fig. 1) { } Ω = (x, y) 0 < x < 1, 0 < y < 1 and its boundary ∂Ω = Γe ∪ Γw ∪ Γs ∪ Γn . Suppose that the functions a(x, y), f (x, y) satisfy ¯ a, f ∈ C (4,α) (Ω),
(1.3)
a(x, y) ≥ 0,
(1.4)
(α > 0), ¯ ∀(x, y) ∈ Ω,
Date: February 23, 2013. 1991 Mathematics Subject Classification. Primary: 65N35, 65N38; Secondary: 35L10. Key words and phrases. Tailored finite point method, anisotropic diffusion problem, boundary layer, compatibility conditions. † Corresponding author. H. Han was supported by the NSFC Project No. 10971116. Z. Huang was supported by the NSFC Project No. 11071139, the National Basic Research Program of China under the grant 2011CB309705. W. Ying was supported in part by the National Natural Science Foundation of China under Grant DMS– 11101278 and the Young Thousand Talents Program of China. 1
2
H.D. HAN, Z.Y. HUANG AND W.J. YING
p
4
n
p
3
w
e
p1
s
p2
Figure 1. The rectangle domain Ω = (0, 1) × (0, 1). ¯ with k be and 0 < ε ≤ 1. Here, C (k,α) (Ω) stands for the H¨older Space of index (k, α) on Ω a non-negative integer and α be a positive real number in the interval (0, 1]. Let {pl }4l=1 be ¯ (see Fig. 1). the four corners of the rectangle domain Ω When 0 < ε ≪ 1, problem (1.1)–(1.2) is an anisotropic diffusion problem and its solution may have boundary layers on a portion of ∂Ω, i.e. at x = 0 and x = 1. These layers are characterized by rapid transitions in the solution and difficult to capture in a numerical approximation without using a large number of unknowns. Also, such layers tend to cause spurious oscillations in a numerical solution to the problem. The anisotropic diffusion problem appears in a large variety of applications such as image processing [17, 22], chemical reactions [20] and anisotropic materials [15]. It has been numerically solved with asymptotic preserving method [2], least-square finite element method [16], mimetic finite difference methods [13], etc. A similar anisotropic but first-order ordinary differential system is solved with a finite difference method on Shishkin meshes [14]. This work proposes a tailored finite point method (TFPM) [6, 7, 8, 9, 10, 12] to numerically solve problem (1.1)–(1.2). The tailored finite point method was first proposed by Han, Huang and Kellogg [10] for the numerical solutions of singular perturbation problems of second-order elliptic equations with constant coefficients. Later, Han and Huang [6, 7, 8, 9] and Shih, Kellogg et al. [18, 19] systematically develop this method for the nonhomogeneous reactiondiffusion, convection-diffusion and convection-diffusion-reaction problems. The basic idea of the tailored finite point method is that the numerical scheme is tailor-made at each point/cell based on the local properties of the solution of the given problem. The method is expected to capture boundary layers associated with the problem (1.1)–(1.2) even with coarse grids. The remainder of this paper is organized as follows. Section 2 makes studies on the differentiability of the solution to problem (1.1)–(1.2) under some compatibility conditions. Section 3 derives a semi-discrete approximation of problem (1.1)–(1.2). The details of the TFPM for the semi-discrete and fully-discrete approximations of the problem are described in Section 4. Section 4 gives a proof for the uniform convergence property of the approximate solution. Finally, numerical examples supporting the theory and demonstrating the efficiency and reliability of the method are presented in Section 5.
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
3
2. The differentiability of the solution to problem (1.1)–(1.2) The differentiability of uε (x, y), the solution to problem (1.1)–(1.2), in the open domain Ω only depends on the differentiability of the given functions a(x, y) and f (x, y). By condition (1.3) and an interior a prior estimate for the elliptic equation (1.1), it is known that [11] uε ∈ C (6,α) (Ω).
(2.1)
¯ also But the differentiability of the solution to problem (1.1)–(1.2) on the closed domain Ω depends on some compatibility conditions at the corner points (refer to Grisvard [5] and Han & Kellogg [11]). For problem (1.1)–(1.2), suppose that the source term satisfies the following compatibility conditions } f (pl ) = 0, l = 1, · · · , 4, (2.2) fxx (pl ) = 0, fyy (pl ) = 0, l = 1, · · · , 4. ¯ Theorem 3.2 in Han & Kellogg [11] implies that uε (x, y) ∈ C (4,α) (Ω). Remark 2.1. If the compatibility conditions (2.2) (or part of them) do not hold, then the partial derivatives of the solution uε (x, y) may not be continuous at four corners [5, 11]. ¯ denote by ∥v∥∞ = max |v(x, y)| its For any function v(x, y) that is continuous on Ω, ¯ (x,y)∈Ω
maximum norm. Let Mk1 ,k2 Fk1 ,k2 Ak1 ,k2
k +k
∂ 1 2 uε
=
∂xk1 ∂y k2 ,
k +k ∞
∂ 1 2f
=
∂xk1 ∂y k2 ,
k +k ∞
∂ 1 2a
=
∂xk1 ∂y k2 , ∞
k1 ≥ 0, k2 ≥ 0, k1 + k2 ≤ 4, k1 ≥ 0, k2 ≥ 0, k1 + k2 ≤ 4, k1 ≥ 0, k2 ≥ 0, k1 + k2 ≤ 4,
and M0 = M0,0 , A0 = A0,0 , F0 = F0,0 . Assumption 2.1. In the rest of this section, assume that the constants Fk1 ,k2 and Ak1 ,k2 are bounded uniformly with respect to ε. Lemma 2.1. There exists a constant C independent of ε such that F0 . (2.3) M0 ≤ C ≡ 8 Proof. This result follows directly from the maximum principle of elliptic problem (1.1)–(1.2). Define the differential operator ∂2 ∂2 Lε ≡ −ε − +a ∂x2 ∂y 2 and introduce an auxiliary function F0 φ(y) = y(1 − y), 2 2
(2.4)
4
H.D. HAN, Z.Y. HUANG AND W.J. YING
which depends on variable y only. It is easy to check that 0 ≤ φ(y) ≤ C ≡
F0 , 8
∀y ∈ [0, 1],
and Lε φ(y) = F0 + a(x, y)φ(y) ≥ 0,
¯ ∀(x, y) ∈ Ω.
Let v± = uε ± φ. We have v+ |∂Ω ≥ 0, and
v− |∂Ω ≤ 0
( ) Lε v± (x, y) = f (x, y) ± Lε φ(y) = f (x, y) ± F0 + a(x, y)φ(y) .
By the definition of F0 and the assumption on the coefficient a(x, y), we get Lε v+ (x, y) ≥ 0,
Lε v− (x, y) ≤ 0,
∀(x, y) ∈ Ω.
The maximum principle further implies v+ (x, y) ≥ 0,
v− (x, y) ≤ 0,
¯ ∀(x, y) ∈ Ω.
As a result −φ(y) ≤ uε (x, y) ≤ φ(y),
¯ ∀(x, y) ∈ Ω. □
Then inequality (2.3) follows immediately. Lemma 2.2. The following inequalities hold 2 δ 1 M0,j−1 + M0,j+1 , ∀δ ∈ (0, ], j = 1, 2, 3; δ 2 2 2 j−1 δ j+1 1 ≤ ε Mj−1,0 + ε Mj+1,0 , ∀δ ∈ (0, ], j = 1, 2, 3. δ 2 2
M0,j ≤ εj Mj,0
(2.5) (2.6)
Proof. At first, we prove estimate (2.5). By the definition of M0,1 , there exists a point ¯ such that (x0 , y0 ) ∈ Ω, ∂uε (x0 , y0 ) . M0,1 = ∂y If 0 ≤ y0 ≤
1 2
(when
1 2
< y0 ≤ 1, the proof is analogous), for any δ ∈ (0, 12 ], we have
δ 2 ∂ 2 uε ∂uε (x0 , y0 ) + (x0 , y0 + ξ), uε (x0 , y0 + δ) = uε (x0 , y0 ) + δ ∂y 2 ∂y 2
with ξ ∈ [0, δ].
(2.7)
From (2.7), estimate (2.5) for j = 1 is obtained. For the case j = 2, 3, the proof is similar. Furthermore, following the same lines as those for (2.5), we can prove estimate (2.6). □ Let ∂ j uε (x, y) , ∂y j ∂ j uε (x, y) wj (x, y) = , ∂xj vj (x, y) =
j = 1, 2, 3, 4; j = 1, 2, 3, 4.
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
5
It can be verified that function v2 (x, y) satisfies the following problem ∂ 2 v2 ∂ 2 v2 −ε2 2 − + av2 = fyy (x, y) − ayy (x, y)uε − 2ay (x, y)v1 , ∂x ∂y 2 v2 Γe = v2 Γw = 0, v2 Γn = −f (x, 1), v2 Γs = −f (x, 0).
∀(x, y) ∈ Ω,
(2.8) (2.9) (2.10)
The boundary conditions (2.10) are from equation (1.1). Similar to the proof for Lemma 2.1, introducing an auxiliary function ψ(y) ≡ F0 + (F0,2 + A0,2 M0 + 2A0,1 M0,1 ) we have 0 ≤ ψ(y) ≤ F0 +
y(1 − y) , 2
F0,2 + A0,2 M0 + 2A0,1 M0,1 , 8
∀y ∈ [0, 1],
∀y ∈ [0, 1]
(2.11)
and Lε ψ(y) = F0,2 + A0,2 M0 + 2A0,1 M0,1 + a(x, y)ψ(y). Then it is easy to check that Lε (v2 + ψ) ≥ 0, and Lε (v2 − ψ) ≤ 0. The maximum principle implies M0,2 = max |v2 (x, y)| ≤ max |ψ(y)|. (x,y)∈Ω
y∈[0,1]
By (2.11) and (2.5) in Lemma 2.2, we obtain F0,2 + A0,2 M0 + 2A0,1 M0,1 8 F0,2 + A0,2 M0 A0,1 δA0,1 ≤ F0 + + M0 + M0,2 , 8 2δ 8
M0,2 ≤ F0 +
{
Let δ = min we arrive at
As we have (let δ =
1 2
F0,2 M0,2 ≤ 2F0 + + 4 in equation (2.5))
(
1 4 , 2 A0,1
(2.12)
} ,
A0,2 A0,1 + 4 δ
M0,1 ≤ 4M0 +
1 ∀δ ∈ (0, ]. 2
(2.13) ) M0 ≡ C2 .
(2.14)
M0,2 , 4
we get the following Lemma. Lemma 2.3. There exists a constant C independent of ε such that M0,1 ≤ C,
M0,2 ≤ C.
Furthermore, from equation (1.1) and inequality (2.3), we have
(2.15)
6
H.D. HAN, Z.Y. HUANG AND W.J. YING
Lemma 2.4. There exists a constant C independent of ε such that εj Mj,0 ≤ C,
j = 1, 2.
(2.16)
Now we try to get a prior estimate of v4 . First, we discuss the boundary condition of v4 at Γs and Γn . From equation (2.8), we have v4 =
4 ∂ 2 v2 2 ∂ uε = −ε + a(x, y)v2 − fyy (x, y) + ayy (x, y)uε + 2ay (x, y)v1 . ∂y 2 ∂y 2 ∂x2
(2.17)
Differentiating equation (1.1) with respect to x twice, we obtain ∂ 4 uε = −ε2 w4 + a(x, y)w2 + axx (x, y)uε + 2ax (x, y)w1 − fxx (x, y). ∂y 2 ∂x2
(2.18)
Substituting (2.18) into (2.17) and using boundary conditions wj |Γn = wj |Γs = 0 (j = 1, 2, 3, 4) yield ( ) v4 Γn = ε2 fxx + av2 − fyy + 2ay v1 Γn ≡ gn4 (x), (2.19) ( 2 ) v4 Γs = ε fxx + av2 − fyy + 2ay v1 Γn ≡ gs4 (x). (2.20) This means, by Lemma 2.3, we can find a constant C independent of ε such that ∥gn4 ∥Γn ,∞ ≤ C,
∥gs4 ∥Γs ,∞ ≤ C.
(2.21)
The function v4 satisfies ∂ 2 v4 ∂ 2 v4 −ε2 2 − + av4 = fyyyy − ayyyy uε − 4ayyy v1 − 6ayy v2 − 4ay v3 , ∀(x, y) ∈ Ω, (2.22) ∂x ∂y 2 v2 Γe = v2 Γw = 0, (2.23) v2 Γn = gn4 (x), v2 Γs = gs4 (x). (2.24) As before, applying the maximum principle to problem (2.22)–(2.24), we get Lemma 2.5. There exists a constant C independent of ε such that
4
∂ uε
∂y 4 ≤ C.
(2.25)
∞
Combining Lemma 2.1 – Lemma 2.5, we obtain Theorem 2.1. There exists a constant C independent of ε such that
j
∂ uε
∂y j ≤ C, j = 0, 1, · · · , 4.
(2.26)
∞
Similarly, we have Theorem 2.2. There exists a constant C independent of ε such that
j
j ∂ uε
ε
∂xj ≤ C, j = 0, 1, · · · , 4. ∞
Finally, we get the theorem
(2.27)
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
7
Theorem 2.3. There exists a constant C independent of ε such that
j ∂ j1 +j2 uε
ε 1
∂xj1 ∂y j2 ≤ C, 0 ≤ j1 + j2 ≤ 4.
(2.28)
∞
3. A semi-discrete approximation of problem (1.1)–(1.2) The solution of problem (1.1)–(1.2) is smooth in y-direction. It is expected the gradient of the solution in y-direction will in general have small variations provided that the coefficients aj (x) and the source term fj (x) are nice functions. We can discretize the problem in ydirection by the standard finite difference method, which allows us to transform the secondorder PDE into an ODE system. It is much like the method of lines. Taking a positive integer M , let ∆y = M1 and yj = j∆y,
j = 0, 1, · · · , M.
We discretize the second-order partial derivative with respect to the y variable in problem (1.1)-(1.2) with the centered three-point finite difference. This gives us the following semidiscrete system −ε2
d2 uεj (x) uεj−1 (x) − 2uεj (x) + uεj+1 (x) − + aj (x)uεj (x) = fj (x), 1 ≤ j ≤ M − 1, (3.1) dx2 (∆y)2 uεj (0) = uεj (1) = 0, 1 ≤ j ≤ M − 1, (3.2) uε0 (x) = uεM (x) = 0,
∀x ∈ [0, 1],
(3.3)
with aj (x) = a(x, yj ),
fj (x) = f (x, yj ),
j = 1, · · · , M − 1.
Here, uεj (x) is the finite difference approximation of uε (x, yj ) for each j ∈ {0, 1, · · · , M }. Let Ujε (x) = uε (x, yj ),
j = 0, 1, · · · , M
(3.4)
be the exact solution at the discrete points {yj }M The differentiability of functions j=0 . { ε } Uj (x), j = 0, 1, · · · , M that we derived in Section 2 allows us to make local Taylor expansions for them. As a result, the functions satisfy the modified semi-discrete system −ε2
ε ε (x) − 2Ujε (x) + Uj+1 d2 Ujε (x) Uj−1 (x) − + aj (x)Ujε (x) = fj (x) + Rjε (x), (3.5) dx2 (∆y)2 1 ≤ j ≤ M − 1,
Ujε (0) = Ujε (1) = 0, ε (x) = 0, U0ε (x) = UM
1 ≤ j ≤ M − 1, (3.6) ∀x ∈ [0, 1],
(3.7)
with the remainder term Rjε (x) on the order of ∆y 2 in that there is a constant C independent of ε such that (3.8) |Rjε (x)| ≤ C(∆y)2 , 1 ≤ j ≤ M − 1. The estimate (3.8) for the remainder term results from the application of Theorems 2.1–2.3.
8
H.D. HAN, Z.Y. HUANG AND W.J. YING
Lemma 3.1. Let ( )T ε vε (x) = v1ε (x), v2ε (x), · · · , vM −1 (x) be a vector-valued function. If for any 0 < x < 1 and 1 ≤ j ≤ M − 1, ε ε d2 vjε (x) vj−1 (x) − 2vjε (x) + vj+1 (x) −ε − + aj (x)vjε (x) ≤ 0 2 2 dx (∆y) 2
then we have { vjε (x)
≤ max
0,
max
{
1≤j≤M −1
vjε (0), vjε (1)
}
{ } ε , max v0ε (x), vM (x)
}
0≤x≤1
.
(3.9)
Proof. First, we prove that if ε ε d2 vjε (x) vj−1 (x) − 2vjε (x) + vj+1 (x) −ε − + aj (x)vjε (x) < 0, 2 2 dx (∆y) 2
(3.10)
then vjε (x) can not attain its maximum non-negative value in the interior of interval 0 < x < 1 and for 1 ≤ j ≤ M − 1. Assume that vjε (x) attains its maximum non-negative value at point x = x∗ for some 1 ≤ j = j ∗ ≤ M − 1. This indicates ε2
d2 vjε∗ (x∗ ) ≤ 0 and dx2
vjε∗ −1 (x∗ ) − 2vjε∗ (x∗ ) + vjε∗ +1 (x∗ ) ≤ 0. (∆y)2
The condition (3.10) further implies 0 ≤ aj ∗ (x∗ )vjε∗ (x∗ ) < ε2
d2 vjε∗ (x∗ ) vjε∗ −1 (x∗ ) − 2vjε∗ (x∗ ) + vjε∗ +1 (x∗ ) + ≤ 0. dx2 (∆y)2
This is a contradiction. As a result, the inequality (3.9) holds. In the next, we prove the lemma for the case ε ε d2 vjε (x) vj−1 (x) − 2vjε (x) + vj+1 (x) −ε − + aj (x)vjε (x) ≤ 0. dx2 (∆y)2 2
Let 1 1 ξj = − 8 2
( )2 1 yj − 2
(3.11)
and wjη (x) = vjε (x) − η ξj
with η > 0 be a small positive parameter. It obvious that ξj ≥ 0 and
ξj−1 − 2ξj + ξj+1 = −1. (∆y)2
(3.12)
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
9
We see that wjη (x) satisfies the condition (3.10) since η η (x) d2 wjη (x) wj−1 (x) − 2wjη (x) + wj+1 −ε − + aj (x)wjη (x) 2 2 dx (∆y) 2 ε ε ε d vj (x) vj−1 (x) − 2vjε (x) + vj+1 (x) ξj−1 − 2ξj + ξj+1 2 = −ε − + aj (x)vjε (x) + η − η aj ξj 2 2 dx (∆y) (∆y)2 ε ε d2 vjε (x) vj−1 (x) − 2vjε (x) + vj+1 (x) ξj−1 − 2ξj + ξj+1 ε ≤ −ε2 − + a (x)v (x) + η j j dx2 (∆y)2 (∆y)2 ε ε d2 vjε (x) vj−1 (x) − 2vjε (x) + vj+1 (x) = −ε2 − + aj (x)vjε (x) − η ≤ −η < 0. 2 2 dx (∆y) 2
Here, we used the condition (3.11). The previous conclusion implies { } { η } η η η η wj (x) ≤ max 0, max wj (0), wj (1) , max {w0 (x), wM (x)} . 1≤j≤M −1
0≤x≤1
(3.13)
Since the parameter η can be arbitrarily small, letting η → 0, we get the inequality (3.9) for vjε (x). □ Now we are ready to prove the following theorem. Theorem 3.1. Let )T uε1 (x), · · · , uεM −1 (x) , ( )T ε Uε (x) = U1ε (x), · · · , UM . −1 (x)
U ε (x) =
(
The following error estimate holds ∥U ε − Uε ∥∞ ≤ C(∆y)2 ,
(3.14)
where C is a constant independent of ε. Proof. Subtracting the semi-discrete system (3.1)-(3.3) by the semi-discrete system (3.5)(3.7) yields a new semi-discrete system −ε2
d2 eεj (x) eεj−1 (x) − 2eεj (x) + eεj+1 (x) − + aj (x)eεj (x) = Rjε (x), 1 ≤ j ≤ M − 1, dx2 (∆y)2 eεj (0) = eεj (1) = 0, 1 ≤ j ≤ M − 1, eε0 (x) = eεM (x) = 0,
∀x ∈ [0, 1],
for the vector-valued error function )T ( eε (x) = eε1 (x), eε2 (x), · · · , eεM −1 (x) = Uε − U ε with eεj (x) = Ujε (x) − uεj (x). Let
[
1 1 φj = − 8 2
( )2 ]
ε 1
Rj (x) yj − ∞ 2
10
H.D. HAN, Z.Y. HUANG AND W.J. YING
with Rjε (x) ∞ = −ε2
max
1≤j≤M −1, 0≤x≤1
|Rjε (x)|. Let vjε (x) = eεj (x) − φj . We have
ε ε
ε (x) − 2vjε (x) + vj+1 (x) d2 vjε (x) vj−1 ε ε
Rj (x) ≤ 0. − + a (x)v (x) ≤ R (x) − j j j ∞ dx2 (∆y)2
Lemma 3.1 implies vjε (x) ≤ 0 for all j ′ s and x ∈ [0, 1]. Thus we obtain eεj (x) ≤ φj ≤ max φj = 0≤j≤M
1
Rjε (x) . ∞ 8
(3.15)
Let wjε (x) = −eεj (x) − φj . We get ε ε
ε (x) − 2wjε (x) + wj+1 (x) d2 wjε (x) wj−1 ε ε
Rj (x) ≤ 0. − + a (x)w (x) ≤ −R (x) − −ε j j j ∞ dx2 (∆y)2 2
Similarly, Lemma 3.1 implies wjε (x) ≤ 0 for all j ′ s and x ∈ [0, 1]. We arrive at −eεj (x) ≤ φj ≤ max φj = 0≤j≤M
1
Rjε (x) . ∞ 8
Combining (3.15) and (3.16), we get
ε
ej (x) ≤ 1 Rjε (x) . ∞ ∞ 8 Together with (3.8), we finally obtain the error estimate (3.14).
(3.16)
□
4. A tailored finite point method for problem (3.1)–(3.3) The problem (3.1)–(3.3) can be rewritten as −ε2
d2 U ε + A(x) U ε = F (x), dx2 U ε (0) = U ε (1) = 0,
∀x ∈ (0, 1),
(4.1) (4.2)
where A(x) is an (M − 1) × (M − 1) matrix function, F (x) is an (M − 1) dimensional vector function, 2 1 − 0 (∆y)2 + a1 (x) (∆y)2 1 2 .. . − + a1 (x) 2 2 (∆y) (∆y) A(x) = 1 .. .. . . − 2 (∆y) 1 2 0 − + a (x) M −1 (∆y)2 (∆y)2 f1 (x) f2 (x) . F (x) = .. . fM −1 (x)
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
11
Taking a positive integer N , let ∆x =
1 N
and
xi = i∆x,
i = 0, 1, · · · , N.
We will introduce piecewise constant matrix functions Ah (x) and Fh (x) for 0 ≤ x ≤ 1, which are approximations of A(x) and F (x) respectively. On each interval [xi , xi+1 ], taking the midpoint x∗i ∈ [xi , xi+1 ] (i.e. , x∗i = xi +x2 i+1 ), define the matrix functions Ah (x) and Fh (x) by Ah (x) = A(x∗i ),
for x ∈ (xi , xi+1 ),
Fh (x) = F (x∗i ),
for x ∈ (xi , xi+1 ).
This means the matrix functions Ah (x) and Fh (x) are constants on each subinterval (xi , xi+1 ) (i = 0, 1, · · · , N − 1) and have jump discontinuities at points x = xi (i = 1, 2, · · · , N − 1). Now we consider the following approximation of problem (4.1)–(4.2), d2 Uhε + Ah (x) Uhε = Fh (x), 2 dx Uhε (0) = Uhε (1) = 0, [ ε] dUh ε [Uh ] x=xi = = 0, dx x=xi −ε2
∀x ∈ (xi , xi+1 ),
i = 0, · · · , N − 1
(4.3) (4.4) (4.5)
with the jumps defined by [Uhε ] x=xi = Uhε (xi + 0) − Uhε (xi − 0), [ ε] dUh dUhε dU ε = (xi + 0) − h (xi − 0). dx x=xi dx dx
(4.6) (4.7)
Here, the semi-discrete solution ( )T ε,h Uhε (x) = uε,h 1 (x), · · · , uM −1 (x) is an approximation of U ε (x). Correspondingly, each component uε,h j (x) of the vector-valued ε function is also an approximation of uj (x). Lemma 4.1. The problem (4.3)–(4.5) has a unique solution Uhε (x) and the solution satisfies the following estimate ∥Uhε ∥∞ ≤ C,
(4.8)
where C is a constant independent of ε. Lemma 4.2. The following uniform error estimate holds ∥U ε − Uhε ∥∞ ≤ C∆x, where C is a constant independent of ε and ∆x.
(4.9)
12
H.D. HAN, Z.Y. HUANG AND W.J. YING
For x ∈ [0, 1], we define yj+1 − y ε,h y − yj ε,h uε,h (x, y) = uj (x) + u (x), for y ∈ [yj , yj+1 ], ∆y ∆y j+1
(j = 0, · · · , M − 1)
where ε,h uε,h 0 (x) = uM (x) ≡ 0.
Combining with Theorem 3.1, we get the following error estimate: Theorem 4.1. Suppose that uε is the solution of problem (1.1)–(1.2), uε,h is defined above. The following error estimate holds { } ∥uε,h − uε ∥∞ ≤ C ∆x + (∆y)2 , (4.10) where C is a constant independent of ε, ∆x and ∆y. Now we discuss how to get the solution Uhε (x) exactly. We will propose a tailored finite point method to solve it [6, 7, 8, 9, 10, 12]. On each subinterval [xi , xi+1 ], Uhε (x) satisfies −ε2
d2 Uhε (x) + Ai Uhε (x) = Fi , dx2
for x ∈ (xi , xi+1 ),
(4.11)
with Ai = A(x∗i ),
Fi = F (x∗i ).
Let Uhf,i (x) = A−1 i Fi ,
(4.12)
which is a particular solution to the nonhomogeneous equation (4.11). Let Uhi (x) = ξi e
λx ε
∈ RM −1 ,
(4.13)
be a general solution of the homogeneous equation 2 ε 2 d Uh (x) −ε dx2
+ Ai Uhε (x) = 0,
for x ∈ (xi , xi+1 ).
(4.14)
Substituting (4.13) into (4.14), we arrive at A i ξ i = λ2 ξ i .
(4.15)
Since Ai is a symmetric definite real matrix, it has M − 1 positive real eigenvalues 0 < µi1 ≤ µi2 ≤ · · · ≤ µiM −1 , and corresponding eigenvectors Let λik =
√
ξ i1 , ξi2 , · · · , ξ iM −1 .
µik . The general solution of equation (4.11) on the interval [xi , xi+1 ] is given by ( i ( i ) )) M −1 ( ∑ λk (x − xi+1 ) −λk (x − xi ) f,i + − ε Uh (x) = Uh (x) + αk exp + αk exp ξ ik , (4.16) ε ε k=1
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
13
with some constants αk± (k = 1, · · · , M − 1). Let ( ) Ki = ξ i1 , · · · , ξ iM −1 , { ( i ) ( i )} λM −1 (x − xi+1 ) λ1 (x − xi+1 ) i Λ+ (x) = diag exp , · · · , exp , ε ε { ( i ) ( i )} −λM −1 (x − xi ) −λ1 (x − xi ) i Λ− (x) = diag exp , · · · , exp , ε ε ( )T + αi+ = α1+ , α2+ , · · · , αM , −1 ( ) T − αi− = α1− , α2− , · · · , αM . −1 Then equality (4.16) can be rewritten by ∀x ∈ [xi , xi+1 ].
Uhε (x) = Uhf,i (x) + Ki Λi+ (x)αi+ + Ki Λi− (x)αi− ,
(4.17)
At points x = xi , xi+1 , we obtain Uhε (xi ) = Uhf,i (xi ) + Ki Λi+ (xi )αi+ + Ki Λi− (xi )αi− ,
(4.18)
Uhε (xi+1 ) = Uhf,i (xi+1 ) + Ki Λi+ (xi+1 )αi+ + Ki Λi− (xi+1 )αi− . Let
( Ai =
Ki Λi+ (xi ) Ki Λi− (xi ) Ki Λi+ (xi+1 ) Ki Λi− (xi+1 )
(4.19)
) .
Then the coefficients αi± are given by ( ) ( ) f,i i ε ( ) α+ Uh (xi ) − Uh (xi ) −1 = Ai . i ε α− Uh (xi+1 ) − Uhf,i (xi+1 ) Assume that ( We have
(
αi+ αi−
)
( =
)−1 Ai ≡
(
Ai11 Ai21
Ai12 Ai22
) .
Ai11 Uhε (xi ) + Ai12 Uhε (xi+1 ) − (Ai11 + Ai12 ) Uhf,i Ai21 Uhε (xi ) + Ai22 Uhε (xi+1 ) − (Ai21 + Ai22 ) Uhf,i
) .
(4.20)
From (4.16), we get dUhε (xi + 0) d d = Ki Λi+ (xi )αi+ + Ki Λi− (xi )αi− , dx dx dx dUhε (xi − 0) d i−1 d i−1 i−1 . = Ki−1 Λi−1 Λi−1 (xi )α− + (xi )α+ + K dx dx dx − Continuity condition (4.7) on the first-order derivative implies Ki
(4.21) (4.22)
d d d i i−1 d i−1 i−1 . (4.23) Λ+ (xi )αi+ + Ki Λi− (xi )αi− = Ki−1 Λi−1 Λi−1 (xi )α− + (xi )α+ + K dx dx dx dx −
14
H.D. HAN, Z.Y. HUANG AND W.J. YING
Substituting (4.20) into (4.23), we finally get a numerical scheme Ai Uhε (xi−1 ) + Bi Uhε (xi ) + Ci Uhε (xi+1 ) = Fi ,
i = 1, 2, · · · , N − 1,
Uhε (0) = Uhε (1) = 0,
(4.24) (4.25)
where d i−1 i−1 d Λ+ (xi )Ai−1 Λi−1 (xi )Ai−1 11 − K 21 , dx dx − d d i−1 d i−1 d i−1 i−1 − Ki−1 Λ− , (xi )A12 (xi )A22 = Ki Λi+ (xi )Ai11 + Ki Λi− (xi )Ai21 − Ki−1 Λ+ dx dx dx dx d d = Ki Λi+ (xi )Ai12 + Ki Λi− (xi )Ai22 , dx dx ( ) ( ) ( i ) d d i i i i i i = K Λ (xi ) A11 + A12 + Λ− (xi ) A21 + A22 Uhf,i dx + dx ) ( ) ) ( i−1 ( i−1 d i−1 d i−1 i−1 i−1 i−1 Uhf,i−1 . −K Λ+ (xi ) A11 + A12 + Λ− (xi ) A21 + A22 dx dx
Ai = −Ki−1 Bi Ci Fi
5. Numerical examples This section presents numerical examples to demonstrate the efficiency and reliability of the tailored finite point method. In the examples, the error eh (x, y) of a numerical solution is computed in the discrete ℓ2 -norm defined by ∥eh ∥2l2
( ) 2 N ∑ M ∑ 1 1 eh xi + ∆x, yj + ∆y . = ∆x∆y 2 2 i=1 j=1
Example 5.1. First, let a(x, y) = 0, f (x, y) = 100 sin 2πy and ) ( 2π(x−1) 2πx 25 e− ε + e ε sin 2πy. uε (x, y) = 2 1 − π ) 1 − exp(− 2π ε
(5.1)
The function uε given in (5.1) is the exact solution to problem (1.1)–(1.2) with a and f given above. In this example, functions a and f satisfy the conditions given in Assumption 2.1. The numerical results are shown in Figures 2, 3 and Table 1. Second-order convergence rates for the numerical solutions are observed in the discrete l2 norm. Table 1. l2 errors of the numerical solutions for Example 5.1. ∆x and ∆y ε = 0.1 ∥uε − uε,h ∥l2 ε = 0.001 ∥uε − uε,h ∥l2
1 8
1 16
1 32
1 64
8.88E-2 2.23E-2 5.56E-3 1.39E-3 8.88E-2 2.25E-2 5.67E-3 1.43E-3
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
15 −3
x 10 8
2.5 2 3
1 0
6
1
0.008
0.5
0.006
5
0
0.004
4
0.002
3
−1 −0.5 −2 −1
−3 1
7
0.01
1.5
2
0 1
−1.5
2
1
1 −2
0.5
0.5
0.5 0
−2.5
0
0
1
0.5 0
0
(a)
(b) −3
x 10 2
2.5 −3
x 10
2 3 2
1
1
0.5
0
0
1.6
2
1.4
1.5
1.2 1
1
−1 −0.5
0.8
0.5
−2 −1
−3 1
1.8
2.5
1.5
0.6
0 1
−1.5 1
1 −2
0.5
0.5
0.5 0
−2.5
0
0.5 0
(c)
0.4 0.2 0
0
(d)
1 Figure 2. Example 5.1 with ε = 0.1: (a) uε,h with ∆x = ∆y = 32 ; (b) 1 1 |uε − uε,h | with ∆x = ∆y = 32 ; (c) uε,h with ∆x = ∆y = 64 ; (d) |uε − uε,h | 1 with ∆x = ∆y = 64 .
−1
10
−2
10
2 =0.1 =0.001
1 −3
10
−2
10
−1
10
0
10
Figure 3. Convergence rates of ∥uε − uε,h ∥l2 for Example 5.1.
Example 5.2. Now let a(x, y) = 0, f (x, y) = 1 and ) ( ) ( kπx kπ(x−1) ∞ y(1 − y) ∑ (−1)k − 1 e− ε + e ε uε (x, y) = + 2 sin kπy. 2 (kπ)3 1 + exp(− kπ ) ε k=1 The function uε given in (5.2) is the exact solution to problem (1.1)–(1.2).
(5.2)
16
H.D. HAN, Z.Y. HUANG AND W.J. YING
In this example, function f does not satisfy all the compatibility conditions (2.2). The numerical results are shown in Figure 4 (ε = 0.1) and Figure 5 (ε = 0.001) with mesh size ∆x = ∆y = 0.05. The method still yields good approximations on coarse grids (h ≫ ε). −5
x 10 9
0.12
−4
x 10 0.1
8
1 7 0.8
0.1
0.08
6 0.6 5
0.06
0.05
0.4 4 0.2 3
0.04 0 1
0 1 1 0.5
2
0.02
1
0
1
0.5
0.5 0
0
0.5 0
0
0
Figure 4. uε,h (left) and |uε − uε,h |(right), ε = 0.1. −5
x 10 9
0.12
−4
x 10 0.1
8
1 7 0.8
0.1
0.08
6 0.6 5
0.06
0.05
0.4 4 0.2 3
0.04 0 1
0 1 1 0.5 0
1 1
0.5
0.5 0
2
0.02
0
0.5 0
0
0
Figure 5. uε,h (left) and |uε − uε,h |(right) for example 5.2, ε = 0.001. ( )0.1 Example 5.3. Let a(x, y) = 1 and f (x, y) = 64x(1 − x)y(1 − y) . In this example, function f also does not satisfy all the compatibility conditions (2.2). The numerical results are shown in Figure 6 (ε = 0.1) and Figure 7 (ε = 0.001). The method also yields good numerical solutions when h ≫ ε. Example 5.4. Let a(x, y) = xy,
( ( )) ( x) x−1 + exp uε (x, y) = 100x(1 − x)y(1 − y) exp − , ε ε
and f (x, y) = Lε uε . In this example, the function f does not satisfy the conditions given in Assumption 2.1. The numerical results are shown in Figures 8, 9 and Table 2. The second-order convergence rates for the numerical solutions are also observed in the discrete l2 norm.
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
0.1
0.12
0.12
0.1
0.1 0.1
0.08
0.06
0.05
17
0.08
0.06
0.05
0.04
0.04
0 1
0 1 1 0.5
0.02
1 0.5
0.5 0
0
0
0
0
0
Figure 6. Graphs of uε,h for example 5.3, ε = 0.1, ∆y = 1 and ∆x = 200 (right).
1 : 20
∆x =
0.12
(left)
0.1 0.1
0.08
0.06
0.05
1 20
0.12
0.1 0.1
0.02
0.5
0.08
0.06
0.05
0.04
0.04
0 1
0 1 1 0.5
0.02
1 0.5
0.5 0
0
0
0.02
0.5 0
0
Figure 7. Graphs of uε,h for example 5.3, ε = 0.001, ∆y = 1 and ∆x = 200 (right).
1 : 20
∆x =
1 20
(left)
6. Conclusion This work proposes a tailored finite point method (TFPM) for the strongly anisotropic diffusion problem (1.1)–(1.2), which has much slower diffusion rate along one direction than the other in a rectangle domain. It presents analysis on the differentiability of the solution to problem (1.1)–(1.2) under some compatibility conditions and proves a uniform error estimate for the approximate solution to the discrete problem. Numerical experiments show that, for the problem that satisfies the compatibility conditions, the TFPM yields second-order convergence rates in the discrete l2 norm. For the problem that does not satisfy the compatibility and regularity conditions, the method proposed may still offer good approximate solutions.
Table 2. l2 errors of the numerical solutions for Example 5.4 with ∆y = ∆x
1 16
1 32
1 64
1 128
ε = 0.1 ∥uε − uε,h ∥l2 5.98E-3 1.45E-3 3.48E-4 8.58E-5 ε = 0.01 ∥uε − uε,h ∥l2 8.08E-3 2.12E-3 5.30E-4 1.27E-4
1 . 16
18
H.D. HAN, Z.Y. HUANG AND W.J. YING
In this work, for simplicity, the TFPM is proposed only for the problem defined in a rectangle domain. However, the method is by no means limited to the simple domain. Instead, the method can be generalized for problems on complex domains together with a structured grid method such as the kernel-free boundary integral method[21]. This extension will make the TFPM applicable to a larger class of problems. Another potential extension of the TFPM is its application for those problems where an advection or convection term is added into the PDE.
−3
x 10 0.02 4.5
0.018
−3
x 10
0.016
0.025 0.02
4
5 3.5
0.014
4
3
0.012 0.015
3
0.01
0.01 0.005 1
0 0
1
0.006
1
0.004
0 0
2 1.5 1
0.002
0.5
0.5
2.5
2
0.008
0.5
0
0
1
(a)
1 0.5
0.5
0
0
(b) −4
−3
x 10
x 10 1.2
3
−3
−4
x 10
x 10
1
1.2
3
0.8
2
1 0.8
2
0.6
0.6 0.4
1
0.4
1
0.2 1
0 0
1
0 0
0.2 0.5
0.5 1
0.5
0.5 0
0
1
(c)
0
0
(d)
Figure 8. Graphs of |uε − uε,h | for example 5.4, ∆y = 1 1 1 1 ∆x = 16 ; (b) ∆x = 32 ; (c) ∆x = 64 ; (d) ∆x = 128 .
1 , 16
ε = 0.1: (a)
−2
10
=0.1 =0.01
−3
10
2 −4
1
10
−3
10
−2
10
−1
10
Figure 9. Convergence rates of ∥uε − uε,h ∥l2 for Example 5.4.
TFPM FOR ANISOTROPIC DIFFUSION PROBLEM
19
References 1. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, Dover Publications, New York, 1964. 2. P. Degond, F. Deluzet, and C. Negulescu, An Asymptotic Preserving Scheme For Strongly Anisotropic Elliptic Problems, Multiscale Model. Simul. 8 (2010): 645–666. 3. L.C. Evans, Partial Differential Equations, American Mathematical Society, Providence, RI, 1998. 4. I.S. Gradshteyn and I.M. Ryzhik, Table of integrals, series and products, 6th Ed., Academic Press, San Diego, 2000. 5. P. Grisvard, Elliptic problems in nonsmooth domains, Pitman, Boston, 1985. 6. H. Han and Z. Huang, A tailored finite point method for the Helmholtz equation with high wave numbers in heterogeneous medium, J. Comp. Math, 26 (2008): 728-739. 7. H. Han and Z. Huang, Tailored finite point method for a singular perturbation problem with variable coefficients in two dimensions, J. Sci. Comp., 41 (2009): 200-220. 8. H. Han and Z. Huang, Tailored finite point method for steady-state reaction-diffusion equation, Comm. Math. Sci., 8 (2010): 887-899. 9. H. Han and Z. Huang, Tailored finite point method based on exponential bases for convection-diffusionreaction equation, Math. Comput., (to appear). 10. H. Han, Z. Huang and B. Kellogg, A Tailored finite point method for a singular perturbation problem on an unbounded domain, J. Sci. Comp. 36 (2008): 243-261. 11. H. Han and R. B. Kellogg, Differentiability properties of solutions of the equation −ε∆u + ru = f (x, y) in a square, SIAM J. Math. Anal. 21 (1990): 394-408. 12. Z. Huang, Tailored finite point method for the interface problem, Netw. Heterog. Media (2009) 4: 91-106. 13. J. Hyman, J. Morel, M. Shashkov and S. Steinberg, Mimetic finite difference methods for diffusion equations, Computational Geosciences 6 (2002): 333–352. 14. P. Maragatha Meenakshi, S. Valarmathi and J.J.H. Miller, Solving a partially singularly perturbed initial value problem on Shishkin meshes, Appl. Math. Comput., 215 (2010): 3170–3180. 15. J.G. Mullen, Effect of Bardeen-Herring Correlation on Vacancy Diffusion in Anisotropic Crystals, Phys. Rev., 124 (1961 ): 1723–1730. 16. J. Pasdunkorale and I.W. Turner, A second order control-volume finite-element least-squares strategy for simulating diffusion in strongly anisotropic media J comput. math. 23 (2005): 1–16. 17. P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Trans. Pattern Anal. Mach. Intel. 12 (1990): 629–639. 18. Y. Shih, R. B. Kellogg and P. Tsai, A Tailored Finite Point Method for Convection-Diffusion-Reaction Problems, Journal of Scientific Computing, 43 (2010): 239–260. 19. Y. Shih, R. B. Kellogg and Y. Chang, Characteristic Tailored Finite Point Method for ConvectionDominated Convection-Diffusion-Reaction Problems, Journal of Scientific Computing, 47 (2011): 198– 215. 20. G. Turk, Generating Textures on Arbitrary Surfaces Using Reaction-Diffusion, Computer Graphics, 25 (1991): 289–298. 21. W. Ying, C. S. Henriquez, A kernel-free boundary integral method for elliptic boundary value problems, J. Comp. Phy., 227 (2007): 1046-1074. 22. J. Weickert, Anisotropic Diffusion in Image Processing, Teubner, Stuttgart, 1998.
20
H.D. HAN, Z.Y. HUANG AND W.J. YING
(H.D. Han, Z.Y. Huang) Dept. of Mathematical Sciences, Tsinghua University, Beijing 100084, China E-mail address: [email protected] E-mail address: [email protected] (W.J. Ying) Dept. of Mathematics & Institute of Natural Sciences, Shanghai Jiao Tong University, Shanghai 200240, China E-mail address: [email protected]
