- Email: [email protected]
A new finite difference representation for Poisson’s equation on R3 from a contour integral
Applied Mathematics and Computation 217 (2010) 3624–3634
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A new finite difference representation for Poisson’s equation on R3 from a contour integral q Anour F.A. Dafa-Alla a, Kyung-Won Hwang b, Soyoung Ahn c, Philsu Kim c,⇑ a
Chungbuk National University, 410 Seongbongro, Heungduk-gu, Cheongju, Chungbuk 361-763, Republic of Korea Department of Mathematics, Dong-A University, Busan 604-714, Republic of Korea c Department of Mathematics, Kyungpook National University, 1370 Sankyuk-Dong, Buk-gu, Daegu 702-701, Republic of Korea b
a r t i c l e
i n f o
Keywords: Finite difference Contour integral Boundary integral equation Green’s representation formula Poisson’s equation
a b s t r a c t In this paper, we propose a new finite difference representation for solving a Dirichlet problem of Poisson’s equation on R3 . The key idea of the new approach is to represent the solution with a contour integral connecting the nodal values of each local domain centered at each isolated grid node, which is based on the boundary integral equation on the local domain, and calculate the contour integral using a piecewise linear interpolation of the solution. A superconvergence of the scheme is analyzed using a maximum principle and a priori estimate for the finite difference operator. The convergence behavior is comparable to that of standard finite difference methods on rectangle grids, and a superconver Þ [ C 3 ðX Þ, gence property is attained when the solution u is in the function class C 2;a ðX Þ, the standard O(h2) convergence is obtained. 0 < a < 1. Also, if u 2 C 3;1 ðX 2010 Elsevier Inc. All rights reserved.
1. Introduction The model problem to be considered is the Dirichlet problem for Poisson’s equation
Mu ¼ f
in X;
u¼g
on C;
ð1:1Þ
where X is a bounded domain in R3 with the boundary C and is explained in Theorem 4.6, f and g are given functions de Þ as the set of functions whose lth order partial derivscribed below. For a nonnegative integer l, we use the notation C l;a ðX . Then it is known [1,2] that if f 2 C0,a(X), g 2 C(C), then there exists a unique solution atives are Hölder continuous in X Þ of the problem (1.1). Furthermore, if f 2 C l;a ðX Þ, g 2 C lþ2;a ðX Þ, and X is a C lþ2;a ðX Þ domain, then u 2 C 2;a ðXÞ \ CðX lþ2;a u2C ðXÞ. The finite difference method (FDM) is one of the popular methods for solving the problem (1.1) because it has a merit that the discrete system computing can easily be set up with less amount of calculation. Therefore, it has been used as a fundamental technique in many fields of science and engineering. However, it has a restriction that the difference grids are restricted to be the coordinate lines, which are parallel to the coordinate axes, and so it is not useful on the grids of a domain with irregular boundaries which occur in the majority of the practical problems.
q
This work was supported by the Korea Research Foundation Grant No. KRF-2003-041-C00039.
⇑ Corresponding author.
E-mail addresses: [email protected] (A.F.A. Dafa-Alla), [email protected] (K.-W. Hwang), [email protected] (S. Ahn), kimps@ knu.ac.kr (P. Kim). 0096-3003/$ - see front matter 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.10.017
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
3625
Many researches has been devoted to find an effective reformulation the classical FDM [3,4]. The representative reformulation methods are as follows:
the the the the
box integration method, balance method, finite volume method, multi-element balancing method.
These difference methods constructed through the integration interpolation over irregular networks posses many advantages of the finite element method and are effective methods for the numerical computation of partial differential equations. One of the aims of this article is to introduce a new finite difference representation on rectangular grids, which is extendable into any types of grids flexibly. The Shortley–Weller (S–W) approximation [5] is known as the most popular finite difference scheme for the non-uniform grid. Recently, Yamamoto and his colleagues report some new discoveries for the S–W approximation of several second-order elliptic equations on R2 . In particular, they applied the S–W approximation on the rectangular difference grids to the Poisson-type Dirichlet problems [6,7], and showed a super-convergence property. That is, the approximation solution by the S–W method gives the convergence order O(h3) near the boundary and O(h2) far from the boundary when the exact solu Þ. In particular, under the conditions that X is a C2,a domain and the solution belongs to Cl+2,a, tion has the regularity C 3;1 ðX l = 0 or l = 1 and 0 < a < 1, Matsunaga and Yamamoto [6] showed the super-convergence order O(hl+a) far from the boundary and O(hl+1+a) near the boundary. There were several related and continued more works (see [8–12]). The boundary element method (BEM) is also one of popular approximate schemes for solving the elliptic boundary value problem (see [13–16]). In particular, there are many research manuscripts using a combination of FEM and FDM, or FEM and BEM for solving the elliptic boundary value problem. On the other hand, there have been few usage of the combination of FDM and BEM. Most recently, the authors [17,18] tried to find a relation between FDM and BEM, and showed the central finite difference formula for the Laplace’s operator on R2 can be derived from the BEM on the local uniform grid. In particular, the authors introduced a new finite difference representation on a nonuniform rectangle grids in R2 using a local boundary integral representation on each local domain and analyzed that it has a superconvergence property like S–W method [18]. The aim of this paper is to introduce a new finite difference representation for solving the three dimensional Poisson’s Eq. (1.1) by extending the idea of the article [18]. The key of the approach is to represent the solution with a contour integral using the Green’s representation formula on local domains centered at each isolated grid node. Then, the contour integral is calculated exactly by assuming that the solution is a piecewise linear interpolation on each local triangular pyramid. This procedure gives a new seven points representation of a finite difference type. Even though we only discussed the derivation of a seven points formula on rectangle grids, the usage of an arbitrary shape of grid can be allowed, and it is so flexible to mesh. We show that the new scheme satisfies a maximum principle and a priori estimate for the finite difference operator. Based on this facts, we also prove that the scheme has a comparable convergence behavior with that of standard finite dif Þ, then the approximate solution gives O(h2) accuracy at every ference methods. That is, if the exact solution belongs to C 3 ðX point near to the boundary, and O(h) accuracy at every point far away from the boundary. While if the exact solution belongs Þ, 0 < a 6 1, then the approximate solution gives O(h1+a) accuracy at every point near to the boundary, and O(ha) to C 2;a ðX Þ, then the standard convergence accuracy at every point far from the boundary. Also, if the solution belongs to C 3;1 ðX 2 O(h ) is proved. The rest of the manuscript is written as follows. In Section 2, we consider the derivation of a new finite difference representation for the problem (1.1). In Section 3, several properties of the new finite difference operator are examined, and the truncation errors for the difference operator are observed. In Section 4, we analyze the superconvergence with a maximum principle and a priori estimate for the difference operator. We conclude with a discussion of future work that needs to be done. 2. A new finite difference representation Throughout this paper, we assume that the domain X is partitioned with quasiuniform rectangular grids based on the equilateral grid template (for example, see [19,20]). We denote the difference nodes by p1, p2, . . . , pm such that p1 ; . . . ; pm0 2 X and pm0 þ1 ; . . . ; pm 2 C for fixed positive integers m0 and m.!For each isolated node pi (i = 1,!. . . , m0), we ! relate ðkÞ ðkÞ ðkÞ ðkþ2Þ it to six different neighbor nodes pi , k ! = 1, 2, . . . , 6 such that six vectors pi pi are mutually orthogonal; pi pi and pi pi are ! ð5Þ ð6Þ parallel to xk-axis, k = 1, 2; pi pi and pi pi are parallel to x3-axis;
! ! ! ðkÞ ðkþ1Þ ð5Þ pi pi pi pi ¼ lk pi pi ;
! ! ! ðkþ1Þ ðkÞ ð6Þ pi pi pi ¼ kk pi pi ; pi
k ¼ 1; 2; 3;
ð2:1Þ
ðkÞ where lk and kk are positive constants and pi pi ¼ hik ; k ¼ 1; 2; . . . ; 6, where we use the notation jabjfor two points a and b as the Euclidean distance between a and b, and simply let jaj = j0aj. Now, for each nodes pi, we define the local octahedron sub-domain D :¼ Di consisting of eight triangular pyramids T ik , k = 1, 2, . . . , 8, such that four vertices of T i2k are h i h i ð5Þ ðkÞ ðkþ1Þ ð6Þ ðkÞ ðkþ1Þ pi ; pi ; pi ; pi and while four vertices of T i2k1 are pi ; pi ; pi ; pi as shown in Fig. 1(a). Then, the node pi is the
3626
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
Fig. 1. Local domain related to the isolated node pi.
common point of all triangular pyramids in D. In this sense, hereafter, we simply call pi as the center (or center point) of the local domain D. Also, let h be the maximum mesh size among hik for any i and k. 2.1. Contour integrals In this subsection we will introduce a contour integral transformation for a surface integral which is required in the evaluation of surface integrals in the boundary integral formulation on each local octahedron subdomain for the solution of (1.1) ! (see (2.9)). As shown in Fig. 1(b), let T be a triangular pyramid with four vertices [pi, a(1), a(2), a(3)], where three vectors pi aðkÞ , ! ! ! k = 1, 2, 3, are mutually orthogonal satisfying pi að1Þ pi að2Þ ¼ l ! pi að3Þ , where l is a positive constant and let T be the triangular surface of T with three vertices [a(1), a(2), a(3)], where we assume that the boundary @T of T has a counterclockwise ! orientation. For simplicity of notations, we let pi aðkÞ ¼ hk , k = 1, 2, 3, as shown in Fig. 1(b). Let us consider the surface integral
I¼
Z @u @w dsq ; w u @nq @nq T
ð2:2Þ
1 where wðqÞ ¼ jqp and nq is the unit outward normal to the triangle T, and thus a constant vector on the triangle surface T. ij Now we assume that u is a piecewise linear polynomial in the pyramid T interpolating the nodal values ui = u(pi) and ðkÞ ui ¼ uðaðkÞ Þ, (k = 1, 2, 3), and thus ru is constant on T. A general concept for transforming a surface integral like (2.2) into a contour integral can be found in the manuscript [21]. Let b = (b1, b2, b3) be the constant vector ru on the surface T. Then, the linear polynomial u(q) can be written as
uðqÞ ¼ uðpi Þ þ b ðq pi Þ ¼ uðpi Þ þ
3 X
bm ðq pi Þm ;
q2T;
ð2:3Þ
m¼1
where (q)j denotes the jth component of the point q. Also, if we let
"
U
ðmÞ
ðmÞ
ðqÞ ¼ wðqÞ d
þ
ðq pi Þm jq pi j2
#
ðq pi Þ ;
m ¼ 1; 2; 3;
where d(m) = (dm1, dm2, dm3) and dij denotes the Kronecker delta function, then we can check that 3 X
"
bm U ðmÞ ðqÞ ¼ wðqÞ b þ
m¼1
b ðq pi Þ jq pi j2
#
ðq pi Þ :
ð2:4Þ
Also, U(m)(q) satisfies the following lemma. Lemma 2.1. If we let V(m)(q) = U(m) (q pi), then it holds
r V ðmÞ ðqÞ ¼ U ðmÞ ðqÞ;
m ¼ 1; 2; 3:
Proof. For the proof, we consider only the case m = 1. For the other cases m = 2, 3, we may prove the required identities by similar processes. From the definition of U(1), V(1)(q) directly becomes
V ð1Þ ðqÞ ¼ wðqÞdð1Þ ðq pi Þ: By the definitions of cross product and d(1), we get
V ð1Þ ¼ wðqÞð0; ðq pi Þ3 ; ðq pi Þ2 Þ:
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
3627
Hence a straightforward calculation gives
! 2jQ j2 ðQ Þ22 ðQ Þ23 ðQÞ2 ðQ Þ1 ðQ Þ3 ðQ Þ1 ; ; ; jQ j2 jQ j2 jQ j2
r V ð1Þ ¼ wðqÞ Q ¼ q pi :
Therefore, some simple manipulation gives the required identity for the case m = 1. h Due to the identity (2.4) and Lemma 2.1, the integrand of the surface integral I of (2.2) can be written as
w
3 3 @u @w X @w X @w u ¼ bm U ðmÞ nq ui ¼ bm r V ðmÞ nq ui : @nq @nq m¼1 @nq m¼1 @nq
Thus, we can get a contour integral representation by applying Stoke’s theorem to the surface integral I
I¼
3 X
bm IðmÞ ui
IðmÞ ¼
@w dsq ; @nq
T
m¼1
where
Z
I
V ðmÞ ðqÞ Tds;
ð2:5Þ
m ¼ 1; 2; 3:
ð2:6Þ
@T
H
Here, @T denotes the contour integral on the boundary @T of the triangle surface T and T denotes the unit tangent vector on the boundary @T. For the calculation for the coefficients bm of (2.3) (or (2.5)), we try to evaluate the Eq. (2.3) at three vertices points a(k), k = 1, 2, 3. Then the constant vector b must satisfy the linear system
0
ð1Þ
ui ui
1
0 ð1Þ a pi 1
B A ¼ @ að2Þ pi 1 ð3Þ
a pi 1
C B ð2Þ T C Ab ¼ B @ ui ui A; ð3Þ
ui ui
ð1Þ
a pi 2 ð2Þ
a pi 2 ð3Þ
a pi 2
1 ð2Þ
3 C a p i 3 A: ð3Þ
a pi 3 að1Þ pi
! Since three vectors pi aðkÞ , k = 1, 2, 3, are mutually orthogonal, we can see
2 2 2 AAT ¼ diag h1 ; h2 ; h3 ;
hk ¼ pi aðkÞ ;
k ¼ 1; 2; 3:
Thus we can have the following explicit formula for the coefficients bj.
0
0
ð1Þ
ui ui
11
B B h21 CC ðmÞ
ðmÞ B B CC 3 a pi j ui ui X B T B uð2Þ ui CC CC B i bj ¼ B ; BA B h22 CC ¼ 2 hm B B CC m¼1 @ @ uð3Þ u AA i
j ¼ 1; 2; 3:
ð2:7Þ
i
h23
j
Now, we look over the explicit formula for the contour integral I(m) for each m = 1, 2, 3. First we note that
dðmÞ
hj hjþ1 aðjÞ pi aðjþ1Þ pi ¼ aðjþ2Þ pi m ; hjþ2
! where h4 = h1, h5 = h2, and we used the orthogonality of vectors pi aðjÞ . Hence from the simple parametrization q(t) = (a(j+1) a(j))t + a(j), t 2 [0, 1] for each point q on the line aðjÞ aðjþ1Þ , we can see that
V ðmÞ ðqÞ Tds ¼
dðmÞ
ðjþ2Þ
a pi m hj hjþ1 dt aðjÞ pi aðjþ1Þ pi dt qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi: ffi ¼ 2 2 2 2 hj ð1 tÞ2 þ hjþ1 t 2 hjþ2 hj ð1 tÞ2 þ hjþ1 t2
Finally, from the integral identity
Z
1
0
dt 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi log 2 2 2 2 a þb b ð1 tÞ þ a2 t 2
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 a2 þ b þ a a2 þ b þ b ab
;
the contour integral I(m) can be calculated as
I
ðmÞ
ðjþ2Þ
ðjÞ ðjþ1Þ
3 a a þ hj ðaðjÞ aðjþ1Þ þ hjþ1 Þ X a pi m hj hjþ1 log ¼ ; hj hjþ1 hjþ2 jaðjÞ aðjþ1Þ j j¼1
ð2:8Þ
where a(4) = a(1), a(5) = a(2). Summarizing the Eqs. (2.5)–(2.8), we have the following explicit formula for the boundary integral I of (2.2).
3628
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
Lemma 2.2. Assume that u is a piecewise linear polynomial in the pyramid T interpolating the nodal values ui = u(pi) and ðkÞ ui ¼ uðaðkÞ Þ, (k = 1, 2, 3). Then the surface integral I of (2.2) has the explicit formula 3 X ðkÞ ui ui
I¼
k¼1
where a
(0)
ðkþ1Þ ðk1Þ
Z a þ hkþ1 aðkþ1Þ aðk1Þ þ hk1 a hk1 hkþ1 @w log u dsq ; i hk1 hkþ1 hk jaðk1Þ aðkþ1Þ j T @nq
= a(3), a(4) = a(1), h0 = h3 and h4 = h1.
Proof. By substituting the identities (2.7) and (2.8) into the Eq. (2.5), we have
ðjÞ ðjþ1Þ
ðlÞ Z 3 X 3 hj hjþ1 u ui a a þ hj aðjÞ aðjþ1Þ þ hjþ1 X ðlÞ
ðjþ2Þ
i @w a I¼ p p ui dsq : a log i i 2 ðjÞ ðjþ1Þ @n h h q j jþ1 T h h j a a j jþ2 l j¼1 l¼1 !
Since the vectors pi aðjÞ are mutually orthogonal, the above double summation becomes a single one, and hence we can get the required result. h 2.2. Finite difference representation In this subsection we will consider a new finite difference representation with the results of the previous subsection for solving the Dirichlet problem of Poisson’s Eq. (1.1). First we multiply the governing equation of (1.1) by the fundamental 1 solution wðqÞ ¼ jqp , and then integrate over each local domain D ¼ Di , i.e., j i
Z
ðwðqÞMuðqÞ þ f ðqÞwðqÞÞdq ¼ 0;
i ¼ 1; 2; . . . ; m0 :
D
Then we apply Green’s representation formula to the domain integral and get a boundary element formulation (see [13–16]): for each i = 1, 2, . . . , m0,
Z
wðqÞf ðqÞdq ¼ 4puðpi Þ þ D
Z
w
@D
8 Z X @u @w @u @w dsq ¼ 4puðpi Þ þ dsq ; u w u @nq @nq @nq @nq e ik @T k¼1
ð2:9Þ
h i f i2k1 f ik denote triangle surfaces such that three vertices of @T f i2k are pð5Þ ; pðkÞ ; pðkþ1Þ , and while three vertices of @T where @T i i i h i ð6Þ ðkÞ ðkþ1Þ . To approximate the right hand side of (2.9), we assume that u is a piecewise linear polynomial on each are pi ; pi ; pi ðkÞ ðkÞ pyramid T ik interpolating the nodal values ui = u(pi), ui ¼ u pi . Then applying the result of Lemma 2.2 to surface integrals on each pyramid T ik , we can get the following approximation formula. ðkÞ
Theorem 2.3. Let D ¼ Di be the local octahedral sub-domain centered at the node pi with the six vertices pi , k = 1, 2, . . . , 6. Assume that u is a piecewise linear polynomial on each pyramid T ik interpolating the nodal values. Then the right hand side of (2.9) can be calculated as
4puðpi Þ þ
Z @D
@u @w dsQ ¼ w u @nQ @nQ
! ð1Þ ð3Þ ui ui ui ui þ ðGðhi2 ; hi5 ; hi6 Þ þ Gðhi4 ; hi5 ; hi6 ÞÞ hi1 hi3 ! ð2Þ ð4Þ ui ui ui ui ðGðhi1 ; hi5 ; hi6 Þ þ Gðhi3 ; hi5 ; hi6 ÞÞ þ þ hi2 hi4 ! ð5Þ ð6Þ ui ui ui ui þ ðGðhi2 ; hi1 ; hi3 Þ þ Gðhi4 ; hi1 ; hi3 ÞÞ; þ hi5 hi6
where
0
y Gðx; y; zÞ ¼ x@pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi log x2 þ y 2
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2 þ y2 þ x x2 þ y 2 þ y xy
z þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi log x2 þ z2
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 x2 þ z2 þ x x2 þ z 2 þ z A: xz
ð2:10Þ
Proof. Recall that
Z @D
@w dsq ¼ 4p; @nq
which can be obtained by taking u = 1 in Green’s representation formula (2.9). Thus, from two facts (2.1) and Lemma 2.2 with some manipulations, we can get the required result. h
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
3629
For a simple approximation of the volume integral in the left hand side of (2.9), let us define
Wi ¼
Z D
1 dq: jq pi j
ð2:11Þ
Then the volume integral of (2.9) can be simply approximated as
Z
f ðqÞwðqÞdq f ðpi ÞW i :
D
¼ ðu1 ; . . . ; um Þ 2 Rm and each i = 1, 2, . . . , m0, we introduce a new finite difference operator Wi u defined Now, for any vector u by
! ! ð1Þ ð3Þ ð2Þ ð4Þ ui ui ui ui ui ui ui ui þ þ ðGðhi2 ; hi5 ; hi6 Þ þ Gðhi4 ; hi5 ; hi6 ÞÞ þ ðGðhi1 ; hi5 ; hi6 Þ þ Gðhi3 ; hi5 ; hi6 ÞÞ hi1 hi3 hi2 hi4 ! ð5Þ ð6Þ ui ui ui ui þ ðGðhi2 ; hi1 ; hi3 Þ þ Gðhi4 ; hi1 ; hi3 ÞÞ þ hi5 hi6
:¼ Wi u
¼
6 X ðkÞ ðkÞ ðui ui ÞGi ;
say;
ð2:12Þ
k¼1
where G is the function defined in (2.10). Then, Theorem 2.3 gives an approximated system for (2.9) as follows.
f ðpi ÞW i ; Wi u
i ¼ 1; 2; . . . ; m0 :
ð2:13Þ
An asymptotic behavior for the volume integral Wi of (2.11) can be observed by terms of Wi() as follows. Lemma 2.4. The volume integral Wi defined in (2.11) has an asymptotic behavior
Wi
Wi a^ 6
;
i ¼ 1; 2; . . . ; m0 ;
ð2:14Þ
^ = (jp1j, 2, . . . , jpmj2). where a Proof. If we let u(q) = jqj2, then Mu(q) = 6, q 2 X. Hence the approximated system (2.13) directly gives the asymptotic behavior (2.14). h From (2.13) and (2.14), we naturally introduce a new finite difference scheme for the problem (1.1) as follows. ¼ ðu1 ; . . . ; um Þ 2 Rm such that (Ph) Find u
(
^; ¼ 16 f ðpi ÞWi a Wi u ui ¼ gðpi Þ;
i ¼ 1; . . . ; m0 ;
i ¼ m0 þ 1; . . . ; m;
^ = (jp1j, 2, . . . , jpmj2). where a 3. Properties of Wi() and truncation error In this section, we will discuss some properties of the finite difference formula Wi() and then analyze a truncation error for the operator Wi. Let
j ¼ ðp1 Þj ; . . . ; ðpm Þj ; aj ¼ ððp1 Þ2j ; . . . ; ðpm Þ2j Þ; a ^ ¼ jp1 j2 ; . . . ; jpm j2 ; a jk ¼ ðp1 Þj ðp1 Þk ; . . . ; ðpm Þj ðpm Þk : a ðkÞ
Then we first note that from the definition of G in (2.10), the coefficients Gi can observe the following properties of Wi.
ðkÞ
of ui ui in (2.12) are always positive and we
^ and a jk , the difference operator Wi() of (2.12) satisfies the following identities: for each i = 1, 2, . . . , m0 Lemma 3.1. For aj ; aj ; a and j = 1, 2, 3,
Wi aj ¼
6 X ðkÞ ðkÞ pi pi Gi ¼ 0; k¼1
Wi aj ¼
j
6 2 X ðkÞ ðkÞ pi pi Gi ; k¼1
j
Wi a12 ¼ Wi a13 ¼ Wi a23 ¼ 0;
Wi a^ ¼
2 6 6 X X 2 ðkÞ ðkÞ ðkÞ hik Gi : pi pi Gi ¼ k¼1
k¼1
ð3:1Þ
3630
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634 ðkÞ
Proof. We first note that the definition of the nodes pi ð1Þ
ð2Þ
pi pi ¼ ðhi1 ; 0; 0Þ; ð3Þ
pi pi
ðkÞ
Also, the definition of Gi ð1Þ hi1 Gi
¼
ð3Þ hi3 Gi ;
ð5Þ
pi pi ¼ ð0; hi2 ; 0Þ; ð4Þ
¼ ðhi3 ; 0; 0Þ;
pi pi
gives
pi pi ¼ ð0; 0; hi5 Þ; ð6Þ
¼ ð0; hi4 ; 0Þ;
pi pi
¼ ð0; 0; hi6 Þ:
ð3:2Þ
shows that ð2Þ
hi2 Gi
ð4Þ
¼ hi4 Gi ;
ð5Þ
hi5 Gi
ð6Þ
¼ hi6 Gi :
Hence the Eq. (3.2) directly gives the first identity of (3.1). From the relations in (3.2), we see that and hence
ð3:3Þ
ðkÞ pi
ðkÞ pi pi pi ¼ 0 1
2
ðkÞ ðkÞ ðkÞ ðkÞ pi pi ðpi Þ1 ðpi Þ2 ¼ ðpi Þ2 pi pi þ ðpi Þ1 pi pi 1
2
1
2
for any k = 1, 2, . . . , 6. This equation and the first equation of (3.1) shows that
Wi a12 ¼ 0: 13 ¼ 0 ¼ Wi a 23 . The third and fourth identities of (3.1) can be proved from the first one of (3.1) Similarly, we can see that Wi a and the facts
2 2 ðkÞ ðkÞ ðkÞ pi pi ¼ pi ðpi Þ2j 2ðpi Þj pi pi ; j
j
j
Hence we complete the proof.
j ¼ 1; 2; 3:
h
Remark 3.2. From the asymptotic behavior of (2.14), if we define
¼ Mi u ^ ¼ where a
¼ Mi u
6
Wi a^
Wi u ;
i ¼ 1; 2; . . . ; m0 ;
ð3:4Þ
¼ ðu1 ; . . . ; um Þ, then the difference operator Mi() becomes jp21 ; . . . ; jpm j2 and ðu
6 1 X ðkÞ ui ui 2 h k¼1
under the uniform grid (that is, hi1 = hi2 = hi3 = hi4 = ni1 = ni2 = h for some fixed mesh size h). This formula is the same with the is a new genclassical central finite difference formula for Laplace’s operator. So, we may say that the difference formula Wiu eralized one on a non-uniform grid corresponding to Laplace’s operator. ^ and aj, the difference operator Wi() has the following asymptotic properties. Lemma 3.3. For the vectors a
Wi a^ 3Wi aj ¼ Oðh2 Þ;
j ¼ 1; 2; 3:
ð3:5Þ
In particular, for the uniform case hik = h, k = 1, 2, . . . , 6,
Wi a^ 3Wi aj ¼ 0;
j ¼ 1; 2; 3:
ð3:6Þ ðkÞ
Proof. From the definition of G defined in (2.10), we directly see that Gi ¼ Oð1Þ. Thus the equations of (3.1) pffiffiffithirdand pffiffiffi fourth ðkÞ give the asymptotic behavior (3.5). For the uniform case, Gi has the constant value 2 2 log 2 þ 1 . Thus the third and fourth equations of (3.1) show
^¼ 3Wi a1 Wi a
6 2 2 2 X ðkÞ ðkÞ ðkÞ ðkÞ Gi ¼ 0: 2 pi pi pi pi pi pi k¼1
1
2
3
Similarly, we can prove the Eq. (3.6) for the cases j = 2, 3. h of the problem (Ph), we now analyze the asymptotic For the solution u of the problem (1.1) and an approximate solution u ~u Þ, where u ~ ¼ ðuðp1 Þ . . . ; uðpm ÞÞ. For simplicity of the notation, we let hk ¼ pðkÞ behavior of the truncation error Wi ðu i pi . Then Taylor series of u becomes 3 1X @2u ðkÞ u pi ðhk Þi ðhk Þj ðp Þ þ Rk ; uðpi Þ ¼ ruðpi Þ hk þ 2 i;j¼1 @xi @xj i
ð3:7Þ
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
where
Rk ¼
8 3 R1 P > > > > < i;j;l¼1ðhk Þi ðhk Þj ðhk Þl 0
ð1tÞ2 @3 u @xi @xj @xl 2
ðpi þ thk Þdt;
3631
Þ; u 2 C 3 ðX
h 2 i 3 > R1 P 2 > u u > ðhk Þi ðhk Þj 0 ð1 tÞ @x@ @x ðpi þ thk Þ @x@ @x ðpi Þ dt; > : i j i j
Þ: u 2 C 2;a ðX
i;j¼1
ðkÞ
Multiplying both sides of (3.7) by Gi
and summing over the index k, then Lemma 3.1 shows
3 1X @2u 1 ^ ~¼ þ eðpi ; hÞ; a þ Wi u Wi u ðp Þ W a W i j i i 2 j¼1 @x2j 3 or equivalently
3 1X @2u 1 ^ a þ eðPi ; hÞ; ðp Þ W a W i j i 2 j¼1 @x2j i 3
Wi ðu~ uÞ ¼
where e(pi, h) is defined by eðpi ; hÞ ¼
(
eðpi ; hÞ ¼
3
Oðh
ðkÞ k¼1 Rk Gi .
ðkÞ
Since Gi
¼ Oð1Þ, we see that
3
Þ; u 2 C ðX
Oðh Þ; 2þa
P6
ð3:8Þ
Þ;
ð3:9Þ
Þ: u 2 C 2;a ðX
Þ, then from the Taylor expansion Similarly, if u 2 C 3;1 ðX 3 3 1X @2u 1 X @3u 4 ðkÞ uðpi Þ ¼ ruðpi Þ hk þ u pi ðhk Þi ðhk Þj ðpi Þ þ ðhk Þi ðhk Þj ðhk Þl ðp Þ þ Oðh Þ; 2 i;j¼1 @xi @xj 6 i;j;l¼1 @xi @xj @xl i
ð3:10Þ
we get
Wi ðu~ uÞ ¼
3 1X @2u 1 1X @3u 4 ^ a þ ðp Þ W a W C ðpi Þ þ Oðh Þ; i j i jkj i 2 j¼1 @x2j 3 6 jkj¼3 @xk11 @xk22 @xk33
ð3:11Þ
where k = (k1, k2, k3), jkj = k1 + k2 + k3, ki 2 {0, 1, 2, 3} and
C jkj ¼
6 X
ðkÞ
Gi
k1 k2 k3 ðkÞ ðkÞ ðkÞ pi pi pi pi pi pi : 1
k¼1
2
ð3:12Þ
3
ðkÞ
Here we note that if hik = h, k = 1, . . . , 6, then the Eq. (3.2) and the fact that Gi
C jkj ¼ 0;
is constant show
8jkj ¼ 3:
ð3:13Þ
Hence the Eqs. (3.8), (3.9), (3.11), (3.13) and Lemma 3.3 give the following lemma. ¼ ðu1 ; . . . ; um Þ be the solution of the Lemma 3.4 (Truncation error). Let u be the solution of the Dirichlet problem (1.1) and u ~ ¼ ðuðp1 Þ; . . . ; uðpm ÞÞ. Then for a nonuniform case hik – h for some k = 1, . . . , 6, the truncation error discrete system (Ph). Let u becomes 2
Wi ðu~ uÞ ¼ Oðh Þ: On the other hand, for the uniform case hik = h for all k = 1, . . . , 6, it becomes
8 2þa 2;a > > < Oðh Þ; if u 2 C ðXÞ; 3 3 Þ; ~u Þ ¼ Oðh Þ; if u 2 C ðX Wi ðu > > : 4 Þ: Oðh Þ; if u 2 C 3;1 ðX
4. Convergence analysis In this section,we will show the solvability of the discrete system (Ph) using a maximum principle of the difference operator Wi, and then analyze a superconvergence property for the approximate solution. To demonstrate that the discrete system (Ph) has a unique solution, we shall prove that the corresponding homogeneous system has only the trivial solution. ¼ ðv 1 ; . . . ; v m Þ 2 Rm . Lemma 4.1 (Maximum principle). Let us consider v satisfies (a) If v
Wi v P 0;
8i ¼ 1; . . . ; m0 ;
3632
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
then it holds
vi
0
¼ max v i ) i0 2 fm0 þ 1; . . . ; mg: 16i6m
(b) Alternatively, if v
Wi v 6 0;
8i ¼ 1; . . . ; m0 ;
then it holds
vi
0
¼ min
16i6m
v i ) i0 2 fm0 þ 1; . . . ; mg: ðkÞ
Proof. We prove part (a) by contradiction. Let us assume i0 2 {1, . . . , m0}. Then we review that the coefficients Gi in (2.12) < 0 since v iðkÞ v i0 < 0. This is a are always positive for sufficiently small mesh size. Hence the definition of i0 gives that Wi0 v 0 contradiction. We can prove the part (b) with a similar argument to the above. h Remark 4.2. Let us consider the homogeneous system corresponding to the discrete system (Ph), i.e., g(pi) = 0, i = m0 + 1, . . . , m. Then Lemma 4.1 shows that the maximum and minimum of the solution of this homogeneous system vanish. Hence the only solution is the trivial one. Thus it follows by the alternative principle for linear systems that the discrete system (Ph) has a unique solution for arbitrary boundary condition g. Now we observe a priori estimate of the difference operator Wi for the error analysis. For simplicity, we assume that there ðkÞ is a fixed positive integer n such that n < m0; for the indices i = 1, 2, . . . , n, pi 2 X and hik = h (" k = 1, 2 ,. . . , 6); for ðkÞ i = n + 1, . . . , m0, there is at least one index 1 6 k 6 6 with pi 2 @ X and hik – h. Consider the Dirichlet problem of
M/ ¼ 1 in X;
ð4:1Þ
/ ¼ 0 on C; ¼ ð/ ; . . . ; / Þ 2 Rm be the solution of the discrete problem (Ph) corresponding to it. Then Wi ð/Þ satisfies and let / 1 m
Wi ða^Þ
¼ Wi /
6
> 0;
i ¼ 1; 2; . . . ; m0 ;
ð4:2Þ
and hence Lemma 4.1(Maximum principle) shows
/j 6 0;
j ¼ 1; . . . ; m0 ;
and /j ¼ 0;
j ¼ m0 þ 1; . . . ; m:
ð4:3Þ
Using these facts, we can get the following a priori estimate for the difference operator Wi. ¼ ðv 1 ; . . . ; v m Þ 2 Rm , it holds that the difference operator Wi() satisfies the Theorem 4.3 (A priori estimate). For every vector v following inequality
jv j j 6
max jv i j þ K 1 j/j j max jWi v j þ K 2 max jWi v j;
m0 þ16i6m
16i6n
nþ16i6m0
j ¼ 1; . . . ; m;
¼ ð/1 ; . . . ; /m Þ 2 Rm is the solution of the discrete problem (Ph) corresponding to the problem (4.1) and where /
^j K 1 ¼ 6 min jWi a 16i6n
0
1 ;
B K 2 ¼ @ min
nþ16i6m0
X ðkÞ pi 2C
11 ðkÞ C Gi A :
ð4:4Þ
¼ ðb1 ; . . . ; bm Þ 2 Rm be a vector such that Proof. Let b
bj ¼
1; 0;
j ¼ 1; 2; . . . ; m0 ; j ¼ m0 þ 1; . . . ; m;
and þ N2 b ¼ v . Further, we let w ¼ v þ N1 / and v
N1 ¼ 6
max16i6n jWi v j ; ^ min16i6n Wi a
N2 ¼
maxnþ16i6m0 jWi v j : P ðkÞ minnþ16i6m0 pðkÞ 2C Gi i
From the facts of (4.3), bj = 0 = /j for any j = m0 + 1, . . . , m. Thus we see w j ¼ v j for any j = m0 + 1, . . . , m. Also Wi b becomes
8 > < 0; i ¼ 1; 2; . . . ; n; P ðkÞ Wi b ¼ Gi ; i ¼ n þ 1; . . . ; m0 : > : ðkÞ pi 2C
3633
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
Then the Eq. (4.2) shows that
P max jWj v j þ N1 Wi / ¼ max jWj v j 1 þ ¼ Wi v þ N1 Wi / Wi w 16j6n
16j6n
6Wi / P 0; ^ min16j6n Wj a
i ¼ 1; 2; . . . ; n
and
P max jW v j þ N1 W / þ N 2 Wi b ¼ Wi v þ N1 Wi / Wi w j i þ N 2 Wi b nþ16j6m0
0 B ¼ max jWj v j@ nþ16j6m0
Wi b minnþ16j6m0
P
ðkÞ
ðkÞ
pj 2C
Gj
1 C > 0; 1A þ N1 Wi /
i ¼ n þ 1; . . . ; m0 :
to obtain that for all j = 1, 2, . . . , m, Hence we can apply the maximum principle of Lemma 4.1 to w
v j 6 wj þ N1 j/j j N2 bj 6 max wi þ N1 j/j j þ N2 6 16i6m
¼
max wi þ N1 j/j j þ N2 ¼
m0 þ16i6m
max ðv i Þ þ j/j jK 1 max jWi v j þ K 2 max jWi v j ¼
m0 þ16i6m
16i6n
nþ16i6m0
max
m0 þ16i6m
v i þ N1 j/j j þ N2
max jv i j þ j/j jK 1 max jWi v j þ K 2 max jWi v j;
m0 þ16i6m
16i6n
nþ16i6m0
±
where denotes the jth component of w and Kj, j = 1, 2, are defined in (4.4). From the definition of v±, for all j = 1, 2, . . . , m, we have w j
v j 6 v j 6
max jv i j þ j/j jK 1 max jWi v j þ K 2 max jWi v j
m0 þ16i6m
16i6n
nþ16i6m0
which completes the proof. h ¼ ð/ ; . . . ; / Þ 2 Rm is the solution of the discrete problem (Ph) Corollary 4.4. Assume that X is a C2,a (0 < a 6 1) domain and let / 1 m satisfies corresponding to the problem (4.1). Then for sufficiently small mesh size h, /
/j ¼ OðhÞ;
j ¼ n þ 1; . . . ; m0 :
ð4:5Þ
Þ and / = 0 on C. Thus if we let / ~ ¼ ð/ðp Þ; . . . ; /ðp ÞÞ, Lemma 3.4 Proof. First note that the solution / of (4.1) is in C 2;a ðX 1 m shows
( ~ /Þ ¼ Wi ð/
2
Oðh Þ; 2þa
Oðh
hik – h; Þ;
hik ¼ h:
^ ¼ Oðh2 Þ, and also the definition of G defined in (2.10) shows that the coefOn the other hand, Lemma 3.1 shows that Wi a ðkÞ ðkÞ ðkÞ ficients Gi of ui ui defined in (2.12) have an asymptotic behavior Gi ¼ Oð1Þ. Thus the constants K1 and K2 defined in (4.4) have the asymptotic behaviors 2
K 1 ¼ Oðh Þ;
K 2 ¼ Oð1Þ:
Hence Theorem 4.3 shows that there are some constant Ck independent of h such that
2 ~ /Þ þ C 2 max Wi ð/ ~ /Þ þ j/ðp Þj j/j j 6 /j /ðpj Þ þ /ðpj Þ 6 C 1 h j/j j max Wi ð/ j 16i6n
a
2
¼ C 1 h j/j j þ C 2 h þ j/ðpj Þj;
nþ16i6m0
j ¼ 1; . . . ; m0 :
Þ and / = 0 on C, it shows that On the other hand, since the solution / of (4.1) is in CðX
/ðpj Þ ¼ OðhÞ;
j ¼ n þ 1; . . . ; m0 :
Hence for sufficiently enough small h, we can complete the proof. h From this Corollary, the above theorem can be restated as follows. 2,a Corollary 4.5. Assume that X is a C , 0 6 a 6 1, domain and let p1, . . . , pm be the node points of a quasi-uniform rectangle grid ^ ¼ jp1 j2 ; . . . ; jpm j2 . Then for every vector v ¼ ðv 1 ; . . . ; v m Þ 2 Rm , the difference operator Wi() satisfies the following and a
8 > 2 > > jv i j þ h max jWi v j þ max jWi v j ; < O m max 16i6n nþ16i6m0 0 þ16i6m jv j j ¼ > 1 > > : O max jv i j þ h max jWi v j þ max jWi v j ; m0 þ16i6m
16i6n
nþ16i6m0
j ¼ 1; . . . ; n; j ¼ n þ 1; . . . ; m0 :
Hence by Lemma 3.4 and the above Corollary 4.5, we can obtain the following superconvergence property for the discrete solution of (Ph).
3634
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
¼ ðu1 ; . . . ; um Þ be the solutions of (1.1) and (Ph), Theorem 4.6 (Superconvergence). Assume that X is a C2,a domain. Let u and u Þ, then respectively. Then the discrete problem (Ph) has the following convergence property. If u 2 C 3 ðX
(
juðpi Þ ui j ¼ while if u 2 C
2;a
Oðh Þ;
i ¼ n þ 1; . . . ; m0 ;
Þð0 < a 6 1Þ, ðX
juðpi Þ ui j ¼ Also, if u 2 C
i ¼ 1; 2; . . . ; n;
2
(
3;1
OðhÞ;
a
Oðh Þ; Oðh
i ¼ 1; 2; . . . ; n;
aþ1
Þ;
i ¼ n þ 1; . . . ; m0 ;
Þ, we get ðX 2
juðpi Þ ui j ¼ Oðh Þ;
i ¼ 1; 2; . . . ; m0 :
5. Conclusion We developed a new finite difference representation for solving the Dirichlet problem of Poisson’s equation in R3 based on a contour integration and the piecewise linear interpolation of the solution. We also analyze the convergence and show that it has a superconvergence property like two dimensional scheme. Even though we only discussed the derivation of a seven points difference formula on rectangle grid and its analysis, we think that it may be possible to derive a general n-points finite difference formula with n P 7 on a general shape of polyhedron, where each local subdomain consists of a triangular pyramid, because the scheme is based on the boundary element formulation which can be constructed on an arbitrary shape of polyhedron. The work presented here is primarily theoretical, and further numerical calculations need to be done to guarantee its efficiency. This will be the subject of further research. Acknowledgements The authors would like to express our sincere gratitude to reviewers for many helpful comments and valuable suggestions on the first draft of this paper. In particular crucial improvements of Section 4 and several other places are due to anonymous referee’s encouragements and advices. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]
R. Courant, D. Hilbert, Methods of Mathematical Physics, vol. 2, Interscience, New York, 1962. D. Gilbarg, N.S. Trudinger, Elliptic Partial Differential Equations of Second Order, Springer, Berlin, 1977. R. Li, Z. Chen, W. Wu, Generalized Difference Methods for Differential Equations, Marcel Dekker, Inc., 2000. M. Shashkov, Conservative Finite Difference Methods on General Grids, CRC Press, 1996. G.H. Shortley, R. Weller, Numerical solution of Laplace equation, J. Appl. Phys. 9 (1938) 334–348. N. Matsunaga, T. Yamamoto, Superconvergence of the Shortley–Weller approximation for Dirichlet problems, J. Comput. Appl. Math. 116 (2000) 263– 273. T. Yamamoto, On the accuracy of finite difference solutions for Dirichlet problems, RIMS Kokyuroku 1040, RIMS, Kyoto University, 1998, pp. 135–142. Q. Fang, Convergence of finite difference methods for convection diffusion problems with singular solution, J. Comp. Appl. Math. 152 (2003) 119–131. Z. Li, T. Yamamoto, Q. Fang, Superconvergence of solution derivatives for the Shortley–Weller difference approximation of Poisson’s equation. part I: smoothness problem, J. Comp. Appl. Math. 151 (2003) 221–307. Z. Li, H. Hu, Q. Fang, T. Yamamoto, Q. Fang, Superconvergence of solution derivatives for the Shortley–Weller difference approximation of Poisson’s equation. part II: singularity problems, Numer. Funct. Anal. Optim. 24 (2003) 195–221. Z. Li, C. Chien, H. Huang, Effective condition number for finite difference method, J. Comp. Appl. Math. 198 (2007) 208–235. T. Yamamoto, Q. Fang, X. Chen, Superconvergence and non-superconvergence of the Shortley–Weller approximations for Dirichlet problems, Numer Funct. Anal. Optim. 22 (2001) 455–470. K.E. Atkinson, The Numerical Solution of Integral Equations of the Second Kind, Cambridge Univ. Press, 1997. P.K. Banerjee, The Boundary Element Methods in Engineering, Mcgraw-Hill Book Company, 1994. C.A. Brebbia, J.C.F. Telles, L.C. Wrobel, Boundary Element Techniques–Theory and Applications in Engineering, Springer-Verlag, 1984. W. Mclean, Strongly Elliptic Systems and Boundary Integral Equations, Cambridge Univ. Press, 2000. P. Kim, H. Choi, S. Ahn, A new approach for the derivation of a discrete approximation formula on uniform grid for harmonic functions, Kyungpook Math. J. 47 (2007) 529–548. S. Ahn, S. Kim, P. Kim, Local boundary element based a new finite difference scheme for Poisson equations, Appl. Math. Comp. AMC-D-09-01762, Submitted for publication, 2009. W.C. Thacker, A. Conzalez, G.E. Putland, A method for automeshing the construction of irregular computational grids for Storm Suage forecase models, J. Comp. Phys. 37 (1980) 371–387. M. Ribeiro Filho, J.T. Pinho, Automatic and adaptive mesh generation for FEM with applications in microwave, SBMO/IEEE MTT-S IMOC Proceedings 17–20 (2001) 371–387. E.D. Lutz, Systematic derivation of contour integration formulae for Laplace and elastostatic gradient BIE’s, Comp. Mech. 14 (1994) 339–353.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A new finite difference representation for Poisson’s equation on R3 from a contour integral q Anour F.A. Dafa-Alla a, Kyung-Won Hwang b, Soyoung Ahn c, Philsu Kim c,⇑ a
Chungbuk National University, 410 Seongbongro, Heungduk-gu, Cheongju, Chungbuk 361-763, Republic of Korea Department of Mathematics, Dong-A University, Busan 604-714, Republic of Korea c Department of Mathematics, Kyungpook National University, 1370 Sankyuk-Dong, Buk-gu, Daegu 702-701, Republic of Korea b
a r t i c l e
i n f o
Keywords: Finite difference Contour integral Boundary integral equation Green’s representation formula Poisson’s equation
a b s t r a c t In this paper, we propose a new finite difference representation for solving a Dirichlet problem of Poisson’s equation on R3 . The key idea of the new approach is to represent the solution with a contour integral connecting the nodal values of each local domain centered at each isolated grid node, which is based on the boundary integral equation on the local domain, and calculate the contour integral using a piecewise linear interpolation of the solution. A superconvergence of the scheme is analyzed using a maximum principle and a priori estimate for the finite difference operator. The convergence behavior is comparable to that of standard finite difference methods on rectangle grids, and a superconver Þ [ C 3 ðX Þ, gence property is attained when the solution u is in the function class C 2;a ðX Þ, the standard O(h2) convergence is obtained. 0 < a < 1. Also, if u 2 C 3;1 ðX 2010 Elsevier Inc. All rights reserved.
1. Introduction The model problem to be considered is the Dirichlet problem for Poisson’s equation
Mu ¼ f
in X;
u¼g
on C;
ð1:1Þ
where X is a bounded domain in R3 with the boundary C and is explained in Theorem 4.6, f and g are given functions de Þ as the set of functions whose lth order partial derivscribed below. For a nonnegative integer l, we use the notation C l;a ðX . Then it is known [1,2] that if f 2 C0,a(X), g 2 C(C), then there exists a unique solution atives are Hölder continuous in X Þ of the problem (1.1). Furthermore, if f 2 C l;a ðX Þ, g 2 C lþ2;a ðX Þ, and X is a C lþ2;a ðX Þ domain, then u 2 C 2;a ðXÞ \ CðX lþ2;a u2C ðXÞ. The finite difference method (FDM) is one of the popular methods for solving the problem (1.1) because it has a merit that the discrete system computing can easily be set up with less amount of calculation. Therefore, it has been used as a fundamental technique in many fields of science and engineering. However, it has a restriction that the difference grids are restricted to be the coordinate lines, which are parallel to the coordinate axes, and so it is not useful on the grids of a domain with irregular boundaries which occur in the majority of the practical problems.
q
This work was supported by the Korea Research Foundation Grant No. KRF-2003-041-C00039.
⇑ Corresponding author.
E-mail addresses: [email protected] (A.F.A. Dafa-Alla), [email protected] (K.-W. Hwang), [email protected] (S. Ahn), kimps@ knu.ac.kr (P. Kim). 0096-3003/$ - see front matter 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.10.017
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
3625
Many researches has been devoted to find an effective reformulation the classical FDM [3,4]. The representative reformulation methods are as follows:
the the the the
box integration method, balance method, finite volume method, multi-element balancing method.
These difference methods constructed through the integration interpolation over irregular networks posses many advantages of the finite element method and are effective methods for the numerical computation of partial differential equations. One of the aims of this article is to introduce a new finite difference representation on rectangular grids, which is extendable into any types of grids flexibly. The Shortley–Weller (S–W) approximation [5] is known as the most popular finite difference scheme for the non-uniform grid. Recently, Yamamoto and his colleagues report some new discoveries for the S–W approximation of several second-order elliptic equations on R2 . In particular, they applied the S–W approximation on the rectangular difference grids to the Poisson-type Dirichlet problems [6,7], and showed a super-convergence property. That is, the approximation solution by the S–W method gives the convergence order O(h3) near the boundary and O(h2) far from the boundary when the exact solu Þ. In particular, under the conditions that X is a C2,a domain and the solution belongs to Cl+2,a, tion has the regularity C 3;1 ðX l = 0 or l = 1 and 0 < a < 1, Matsunaga and Yamamoto [6] showed the super-convergence order O(hl+a) far from the boundary and O(hl+1+a) near the boundary. There were several related and continued more works (see [8–12]). The boundary element method (BEM) is also one of popular approximate schemes for solving the elliptic boundary value problem (see [13–16]). In particular, there are many research manuscripts using a combination of FEM and FDM, or FEM and BEM for solving the elliptic boundary value problem. On the other hand, there have been few usage of the combination of FDM and BEM. Most recently, the authors [17,18] tried to find a relation between FDM and BEM, and showed the central finite difference formula for the Laplace’s operator on R2 can be derived from the BEM on the local uniform grid. In particular, the authors introduced a new finite difference representation on a nonuniform rectangle grids in R2 using a local boundary integral representation on each local domain and analyzed that it has a superconvergence property like S–W method [18]. The aim of this paper is to introduce a new finite difference representation for solving the three dimensional Poisson’s Eq. (1.1) by extending the idea of the article [18]. The key of the approach is to represent the solution with a contour integral using the Green’s representation formula on local domains centered at each isolated grid node. Then, the contour integral is calculated exactly by assuming that the solution is a piecewise linear interpolation on each local triangular pyramid. This procedure gives a new seven points representation of a finite difference type. Even though we only discussed the derivation of a seven points formula on rectangle grids, the usage of an arbitrary shape of grid can be allowed, and it is so flexible to mesh. We show that the new scheme satisfies a maximum principle and a priori estimate for the finite difference operator. Based on this facts, we also prove that the scheme has a comparable convergence behavior with that of standard finite dif Þ, then the approximate solution gives O(h2) accuracy at every ference methods. That is, if the exact solution belongs to C 3 ðX point near to the boundary, and O(h) accuracy at every point far away from the boundary. While if the exact solution belongs Þ, 0 < a 6 1, then the approximate solution gives O(h1+a) accuracy at every point near to the boundary, and O(ha) to C 2;a ðX Þ, then the standard convergence accuracy at every point far from the boundary. Also, if the solution belongs to C 3;1 ðX 2 O(h ) is proved. The rest of the manuscript is written as follows. In Section 2, we consider the derivation of a new finite difference representation for the problem (1.1). In Section 3, several properties of the new finite difference operator are examined, and the truncation errors for the difference operator are observed. In Section 4, we analyze the superconvergence with a maximum principle and a priori estimate for the difference operator. We conclude with a discussion of future work that needs to be done. 2. A new finite difference representation Throughout this paper, we assume that the domain X is partitioned with quasiuniform rectangular grids based on the equilateral grid template (for example, see [19,20]). We denote the difference nodes by p1, p2, . . . , pm such that p1 ; . . . ; pm0 2 X and pm0 þ1 ; . . . ; pm 2 C for fixed positive integers m0 and m.!For each isolated node pi (i = 1,!. . . , m0), we ! relate ðkÞ ðkÞ ðkÞ ðkþ2Þ it to six different neighbor nodes pi , k ! = 1, 2, . . . , 6 such that six vectors pi pi are mutually orthogonal; pi pi and pi pi are ! ð5Þ ð6Þ parallel to xk-axis, k = 1, 2; pi pi and pi pi are parallel to x3-axis;
! ! ! ðkÞ ðkþ1Þ ð5Þ pi pi pi pi ¼ lk pi pi ;
! ! ! ðkþ1Þ ðkÞ ð6Þ pi pi pi ¼ kk pi pi ; pi
k ¼ 1; 2; 3;
ð2:1Þ
ðkÞ where lk and kk are positive constants and pi pi ¼ hik ; k ¼ 1; 2; . . . ; 6, where we use the notation jabjfor two points a and b as the Euclidean distance between a and b, and simply let jaj = j0aj. Now, for each nodes pi, we define the local octahedron sub-domain D :¼ Di consisting of eight triangular pyramids T ik , k = 1, 2, . . . , 8, such that four vertices of T i2k are h i h i ð5Þ ðkÞ ðkþ1Þ ð6Þ ðkÞ ðkþ1Þ pi ; pi ; pi ; pi and while four vertices of T i2k1 are pi ; pi ; pi ; pi as shown in Fig. 1(a). Then, the node pi is the
3626
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
Fig. 1. Local domain related to the isolated node pi.
common point of all triangular pyramids in D. In this sense, hereafter, we simply call pi as the center (or center point) of the local domain D. Also, let h be the maximum mesh size among hik for any i and k. 2.1. Contour integrals In this subsection we will introduce a contour integral transformation for a surface integral which is required in the evaluation of surface integrals in the boundary integral formulation on each local octahedron subdomain for the solution of (1.1) ! (see (2.9)). As shown in Fig. 1(b), let T be a triangular pyramid with four vertices [pi, a(1), a(2), a(3)], where three vectors pi aðkÞ , ! ! ! k = 1, 2, 3, are mutually orthogonal satisfying pi að1Þ pi að2Þ ¼ l ! pi að3Þ , where l is a positive constant and let T be the triangular surface of T with three vertices [a(1), a(2), a(3)], where we assume that the boundary @T of T has a counterclockwise ! orientation. For simplicity of notations, we let pi aðkÞ ¼ hk , k = 1, 2, 3, as shown in Fig. 1(b). Let us consider the surface integral
I¼
Z @u @w dsq ; w u @nq @nq T
ð2:2Þ
1 where wðqÞ ¼ jqp and nq is the unit outward normal to the triangle T, and thus a constant vector on the triangle surface T. ij Now we assume that u is a piecewise linear polynomial in the pyramid T interpolating the nodal values ui = u(pi) and ðkÞ ui ¼ uðaðkÞ Þ, (k = 1, 2, 3), and thus ru is constant on T. A general concept for transforming a surface integral like (2.2) into a contour integral can be found in the manuscript [21]. Let b = (b1, b2, b3) be the constant vector ru on the surface T. Then, the linear polynomial u(q) can be written as
uðqÞ ¼ uðpi Þ þ b ðq pi Þ ¼ uðpi Þ þ
3 X
bm ðq pi Þm ;
q2T;
ð2:3Þ
m¼1
where (q)j denotes the jth component of the point q. Also, if we let
"
U
ðmÞ
ðmÞ
ðqÞ ¼ wðqÞ d
þ
ðq pi Þm jq pi j2
#
ðq pi Þ ;
m ¼ 1; 2; 3;
where d(m) = (dm1, dm2, dm3) and dij denotes the Kronecker delta function, then we can check that 3 X
"
bm U ðmÞ ðqÞ ¼ wðqÞ b þ
m¼1
b ðq pi Þ jq pi j2
#
ðq pi Þ :
ð2:4Þ
Also, U(m)(q) satisfies the following lemma. Lemma 2.1. If we let V(m)(q) = U(m) (q pi), then it holds
r V ðmÞ ðqÞ ¼ U ðmÞ ðqÞ;
m ¼ 1; 2; 3:
Proof. For the proof, we consider only the case m = 1. For the other cases m = 2, 3, we may prove the required identities by similar processes. From the definition of U(1), V(1)(q) directly becomes
V ð1Þ ðqÞ ¼ wðqÞdð1Þ ðq pi Þ: By the definitions of cross product and d(1), we get
V ð1Þ ¼ wðqÞð0; ðq pi Þ3 ; ðq pi Þ2 Þ:
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
3627
Hence a straightforward calculation gives
! 2jQ j2 ðQ Þ22 ðQ Þ23 ðQÞ2 ðQ Þ1 ðQ Þ3 ðQ Þ1 ; ; ; jQ j2 jQ j2 jQ j2
r V ð1Þ ¼ wðqÞ Q ¼ q pi :
Therefore, some simple manipulation gives the required identity for the case m = 1. h Due to the identity (2.4) and Lemma 2.1, the integrand of the surface integral I of (2.2) can be written as
w
3 3 @u @w X @w X @w u ¼ bm U ðmÞ nq ui ¼ bm r V ðmÞ nq ui : @nq @nq m¼1 @nq m¼1 @nq
Thus, we can get a contour integral representation by applying Stoke’s theorem to the surface integral I
I¼
3 X
bm IðmÞ ui
IðmÞ ¼
@w dsq ; @nq
T
m¼1
where
Z
I
V ðmÞ ðqÞ Tds;
ð2:5Þ
m ¼ 1; 2; 3:
ð2:6Þ
@T
H
Here, @T denotes the contour integral on the boundary @T of the triangle surface T and T denotes the unit tangent vector on the boundary @T. For the calculation for the coefficients bm of (2.3) (or (2.5)), we try to evaluate the Eq. (2.3) at three vertices points a(k), k = 1, 2, 3. Then the constant vector b must satisfy the linear system
0
ð1Þ
ui ui
1
0 ð1Þ a pi 1
B A ¼ @ að2Þ pi 1 ð3Þ
a pi 1
C B ð2Þ T C Ab ¼ B @ ui ui A; ð3Þ
ui ui
ð1Þ
a pi 2 ð2Þ
a pi 2 ð3Þ
a pi 2
1 ð2Þ
3 C a p i 3 A: ð3Þ
a pi 3 að1Þ pi
! Since three vectors pi aðkÞ , k = 1, 2, 3, are mutually orthogonal, we can see
2 2 2 AAT ¼ diag h1 ; h2 ; h3 ;
hk ¼ pi aðkÞ ;
k ¼ 1; 2; 3:
Thus we can have the following explicit formula for the coefficients bj.
0
0
ð1Þ
ui ui
11
B B h21 CC ðmÞ
ðmÞ B B CC 3 a pi j ui ui X B T B uð2Þ ui CC CC B i bj ¼ B ; BA B h22 CC ¼ 2 hm B B CC m¼1 @ @ uð3Þ u AA i
j ¼ 1; 2; 3:
ð2:7Þ
i
h23
j
Now, we look over the explicit formula for the contour integral I(m) for each m = 1, 2, 3. First we note that
dðmÞ
hj hjþ1 aðjÞ pi aðjþ1Þ pi ¼ aðjþ2Þ pi m ; hjþ2
! where h4 = h1, h5 = h2, and we used the orthogonality of vectors pi aðjÞ . Hence from the simple parametrization q(t) = (a(j+1) a(j))t + a(j), t 2 [0, 1] for each point q on the line aðjÞ aðjþ1Þ , we can see that
V ðmÞ ðqÞ Tds ¼
dðmÞ
ðjþ2Þ
a pi m hj hjþ1 dt aðjÞ pi aðjþ1Þ pi dt qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi: ffi ¼ 2 2 2 2 hj ð1 tÞ2 þ hjþ1 t 2 hjþ2 hj ð1 tÞ2 þ hjþ1 t2
Finally, from the integral identity
Z
1
0
dt 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi log 2 2 2 2 a þb b ð1 tÞ þ a2 t 2
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 a2 þ b þ a a2 þ b þ b ab
;
the contour integral I(m) can be calculated as
I
ðmÞ
ðjþ2Þ
ðjÞ ðjþ1Þ
3 a a þ hj ðaðjÞ aðjþ1Þ þ hjþ1 Þ X a pi m hj hjþ1 log ¼ ; hj hjþ1 hjþ2 jaðjÞ aðjþ1Þ j j¼1
ð2:8Þ
where a(4) = a(1), a(5) = a(2). Summarizing the Eqs. (2.5)–(2.8), we have the following explicit formula for the boundary integral I of (2.2).
3628
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
Lemma 2.2. Assume that u is a piecewise linear polynomial in the pyramid T interpolating the nodal values ui = u(pi) and ðkÞ ui ¼ uðaðkÞ Þ, (k = 1, 2, 3). Then the surface integral I of (2.2) has the explicit formula 3 X ðkÞ ui ui
I¼
k¼1
where a
(0)
ðkþ1Þ ðk1Þ
Z a þ hkþ1 aðkþ1Þ aðk1Þ þ hk1 a hk1 hkþ1 @w log u dsq ; i hk1 hkþ1 hk jaðk1Þ aðkþ1Þ j T @nq
= a(3), a(4) = a(1), h0 = h3 and h4 = h1.
Proof. By substituting the identities (2.7) and (2.8) into the Eq. (2.5), we have
ðjÞ ðjþ1Þ
ðlÞ Z 3 X 3 hj hjþ1 u ui a a þ hj aðjÞ aðjþ1Þ þ hjþ1 X ðlÞ
ðjþ2Þ
i @w a I¼ p p ui dsq : a log i i 2 ðjÞ ðjþ1Þ @n h h q j jþ1 T h h j a a j jþ2 l j¼1 l¼1 !
Since the vectors pi aðjÞ are mutually orthogonal, the above double summation becomes a single one, and hence we can get the required result. h 2.2. Finite difference representation In this subsection we will consider a new finite difference representation with the results of the previous subsection for solving the Dirichlet problem of Poisson’s Eq. (1.1). First we multiply the governing equation of (1.1) by the fundamental 1 solution wðqÞ ¼ jqp , and then integrate over each local domain D ¼ Di , i.e., j i
Z
ðwðqÞMuðqÞ þ f ðqÞwðqÞÞdq ¼ 0;
i ¼ 1; 2; . . . ; m0 :
D
Then we apply Green’s representation formula to the domain integral and get a boundary element formulation (see [13–16]): for each i = 1, 2, . . . , m0,
Z
wðqÞf ðqÞdq ¼ 4puðpi Þ þ D
Z
w
@D
8 Z X @u @w @u @w dsq ¼ 4puðpi Þ þ dsq ; u w u @nq @nq @nq @nq e ik @T k¼1
ð2:9Þ
h i f i2k1 f ik denote triangle surfaces such that three vertices of @T f i2k are pð5Þ ; pðkÞ ; pðkþ1Þ , and while three vertices of @T where @T i i i h i ð6Þ ðkÞ ðkþ1Þ . To approximate the right hand side of (2.9), we assume that u is a piecewise linear polynomial on each are pi ; pi ; pi ðkÞ ðkÞ pyramid T ik interpolating the nodal values ui = u(pi), ui ¼ u pi . Then applying the result of Lemma 2.2 to surface integrals on each pyramid T ik , we can get the following approximation formula. ðkÞ
Theorem 2.3. Let D ¼ Di be the local octahedral sub-domain centered at the node pi with the six vertices pi , k = 1, 2, . . . , 6. Assume that u is a piecewise linear polynomial on each pyramid T ik interpolating the nodal values. Then the right hand side of (2.9) can be calculated as
4puðpi Þ þ
Z @D
@u @w dsQ ¼ w u @nQ @nQ
! ð1Þ ð3Þ ui ui ui ui þ ðGðhi2 ; hi5 ; hi6 Þ þ Gðhi4 ; hi5 ; hi6 ÞÞ hi1 hi3 ! ð2Þ ð4Þ ui ui ui ui ðGðhi1 ; hi5 ; hi6 Þ þ Gðhi3 ; hi5 ; hi6 ÞÞ þ þ hi2 hi4 ! ð5Þ ð6Þ ui ui ui ui þ ðGðhi2 ; hi1 ; hi3 Þ þ Gðhi4 ; hi1 ; hi3 ÞÞ; þ hi5 hi6
where
0
y Gðx; y; zÞ ¼ x@pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi log x2 þ y 2
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2 þ y2 þ x x2 þ y 2 þ y xy
z þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi log x2 þ z2
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 x2 þ z2 þ x x2 þ z 2 þ z A: xz
ð2:10Þ
Proof. Recall that
Z @D
@w dsq ¼ 4p; @nq
which can be obtained by taking u = 1 in Green’s representation formula (2.9). Thus, from two facts (2.1) and Lemma 2.2 with some manipulations, we can get the required result. h
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
3629
For a simple approximation of the volume integral in the left hand side of (2.9), let us define
Wi ¼
Z D
1 dq: jq pi j
ð2:11Þ
Then the volume integral of (2.9) can be simply approximated as
Z
f ðqÞwðqÞdq f ðpi ÞW i :
D
¼ ðu1 ; . . . ; um Þ 2 Rm and each i = 1, 2, . . . , m0, we introduce a new finite difference operator Wi u defined Now, for any vector u by
! ! ð1Þ ð3Þ ð2Þ ð4Þ ui ui ui ui ui ui ui ui þ þ ðGðhi2 ; hi5 ; hi6 Þ þ Gðhi4 ; hi5 ; hi6 ÞÞ þ ðGðhi1 ; hi5 ; hi6 Þ þ Gðhi3 ; hi5 ; hi6 ÞÞ hi1 hi3 hi2 hi4 ! ð5Þ ð6Þ ui ui ui ui þ ðGðhi2 ; hi1 ; hi3 Þ þ Gðhi4 ; hi1 ; hi3 ÞÞ þ hi5 hi6
:¼ Wi u
¼
6 X ðkÞ ðkÞ ðui ui ÞGi ;
say;
ð2:12Þ
k¼1
where G is the function defined in (2.10). Then, Theorem 2.3 gives an approximated system for (2.9) as follows.
f ðpi ÞW i ; Wi u
i ¼ 1; 2; . . . ; m0 :
ð2:13Þ
An asymptotic behavior for the volume integral Wi of (2.11) can be observed by terms of Wi() as follows. Lemma 2.4. The volume integral Wi defined in (2.11) has an asymptotic behavior
Wi
Wi a^ 6
;
i ¼ 1; 2; . . . ; m0 ;
ð2:14Þ
^ = (jp1j, 2, . . . , jpmj2). where a Proof. If we let u(q) = jqj2, then Mu(q) = 6, q 2 X. Hence the approximated system (2.13) directly gives the asymptotic behavior (2.14). h From (2.13) and (2.14), we naturally introduce a new finite difference scheme for the problem (1.1) as follows. ¼ ðu1 ; . . . ; um Þ 2 Rm such that (Ph) Find u
(
^; ¼ 16 f ðpi ÞWi a Wi u ui ¼ gðpi Þ;
i ¼ 1; . . . ; m0 ;
i ¼ m0 þ 1; . . . ; m;
^ = (jp1j, 2, . . . , jpmj2). where a 3. Properties of Wi() and truncation error In this section, we will discuss some properties of the finite difference formula Wi() and then analyze a truncation error for the operator Wi. Let
j ¼ ðp1 Þj ; . . . ; ðpm Þj ; aj ¼ ððp1 Þ2j ; . . . ; ðpm Þ2j Þ; a ^ ¼ jp1 j2 ; . . . ; jpm j2 ; a jk ¼ ðp1 Þj ðp1 Þk ; . . . ; ðpm Þj ðpm Þk : a ðkÞ
Then we first note that from the definition of G in (2.10), the coefficients Gi can observe the following properties of Wi.
ðkÞ
of ui ui in (2.12) are always positive and we
^ and a jk , the difference operator Wi() of (2.12) satisfies the following identities: for each i = 1, 2, . . . , m0 Lemma 3.1. For aj ; aj ; a and j = 1, 2, 3,
Wi aj ¼
6 X ðkÞ ðkÞ pi pi Gi ¼ 0; k¼1
Wi aj ¼
j
6 2 X ðkÞ ðkÞ pi pi Gi ; k¼1
j
Wi a12 ¼ Wi a13 ¼ Wi a23 ¼ 0;
Wi a^ ¼
2 6 6 X X 2 ðkÞ ðkÞ ðkÞ hik Gi : pi pi Gi ¼ k¼1
k¼1
ð3:1Þ
3630
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634 ðkÞ
Proof. We first note that the definition of the nodes pi ð1Þ
ð2Þ
pi pi ¼ ðhi1 ; 0; 0Þ; ð3Þ
pi pi
ðkÞ
Also, the definition of Gi ð1Þ hi1 Gi
¼
ð3Þ hi3 Gi ;
ð5Þ
pi pi ¼ ð0; hi2 ; 0Þ; ð4Þ
¼ ðhi3 ; 0; 0Þ;
pi pi
gives
pi pi ¼ ð0; 0; hi5 Þ; ð6Þ
¼ ð0; hi4 ; 0Þ;
pi pi
¼ ð0; 0; hi6 Þ:
ð3:2Þ
shows that ð2Þ
hi2 Gi
ð4Þ
¼ hi4 Gi ;
ð5Þ
hi5 Gi
ð6Þ
¼ hi6 Gi :
Hence the Eq. (3.2) directly gives the first identity of (3.1). From the relations in (3.2), we see that and hence
ð3:3Þ
ðkÞ pi
ðkÞ pi pi pi ¼ 0 1
2
ðkÞ ðkÞ ðkÞ ðkÞ pi pi ðpi Þ1 ðpi Þ2 ¼ ðpi Þ2 pi pi þ ðpi Þ1 pi pi 1
2
1
2
for any k = 1, 2, . . . , 6. This equation and the first equation of (3.1) shows that
Wi a12 ¼ 0: 13 ¼ 0 ¼ Wi a 23 . The third and fourth identities of (3.1) can be proved from the first one of (3.1) Similarly, we can see that Wi a and the facts
2 2 ðkÞ ðkÞ ðkÞ pi pi ¼ pi ðpi Þ2j 2ðpi Þj pi pi ; j
j
j
Hence we complete the proof.
j ¼ 1; 2; 3:
h
Remark 3.2. From the asymptotic behavior of (2.14), if we define
¼ Mi u ^ ¼ where a
¼ Mi u
6
Wi a^
Wi u ;
i ¼ 1; 2; . . . ; m0 ;
ð3:4Þ
¼ ðu1 ; . . . ; um Þ, then the difference operator Mi() becomes jp21 ; . . . ; jpm j2 and ðu
6 1 X ðkÞ ui ui 2 h k¼1
under the uniform grid (that is, hi1 = hi2 = hi3 = hi4 = ni1 = ni2 = h for some fixed mesh size h). This formula is the same with the is a new genclassical central finite difference formula for Laplace’s operator. So, we may say that the difference formula Wiu eralized one on a non-uniform grid corresponding to Laplace’s operator. ^ and aj, the difference operator Wi() has the following asymptotic properties. Lemma 3.3. For the vectors a
Wi a^ 3Wi aj ¼ Oðh2 Þ;
j ¼ 1; 2; 3:
ð3:5Þ
In particular, for the uniform case hik = h, k = 1, 2, . . . , 6,
Wi a^ 3Wi aj ¼ 0;
j ¼ 1; 2; 3:
ð3:6Þ ðkÞ
Proof. From the definition of G defined in (2.10), we directly see that Gi ¼ Oð1Þ. Thus the equations of (3.1) pffiffiffithirdand pffiffiffi fourth ðkÞ give the asymptotic behavior (3.5). For the uniform case, Gi has the constant value 2 2 log 2 þ 1 . Thus the third and fourth equations of (3.1) show
^¼ 3Wi a1 Wi a
6 2 2 2 X ðkÞ ðkÞ ðkÞ ðkÞ Gi ¼ 0: 2 pi pi pi pi pi pi k¼1
1
2
3
Similarly, we can prove the Eq. (3.6) for the cases j = 2, 3. h of the problem (Ph), we now analyze the asymptotic For the solution u of the problem (1.1) and an approximate solution u ~u Þ, where u ~ ¼ ðuðp1 Þ . . . ; uðpm ÞÞ. For simplicity of the notation, we let hk ¼ pðkÞ behavior of the truncation error Wi ðu i pi . Then Taylor series of u becomes 3 1X @2u ðkÞ u pi ðhk Þi ðhk Þj ðp Þ þ Rk ; uðpi Þ ¼ ruðpi Þ hk þ 2 i;j¼1 @xi @xj i
ð3:7Þ
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
where
Rk ¼
8 3 R1 P > > > > < i;j;l¼1ðhk Þi ðhk Þj ðhk Þl 0
ð1tÞ2 @3 u @xi @xj @xl 2
ðpi þ thk Þdt;
3631
Þ; u 2 C 3 ðX
h 2 i 3 > R1 P 2 > u u > ðhk Þi ðhk Þj 0 ð1 tÞ @x@ @x ðpi þ thk Þ @x@ @x ðpi Þ dt; > : i j i j
Þ: u 2 C 2;a ðX
i;j¼1
ðkÞ
Multiplying both sides of (3.7) by Gi
and summing over the index k, then Lemma 3.1 shows
3 1X @2u 1 ^ ~¼ þ eðpi ; hÞ; a þ Wi u Wi u ðp Þ W a W i j i i 2 j¼1 @x2j 3 or equivalently
3 1X @2u 1 ^ a þ eðPi ; hÞ; ðp Þ W a W i j i 2 j¼1 @x2j i 3
Wi ðu~ uÞ ¼
where e(pi, h) is defined by eðpi ; hÞ ¼
(
eðpi ; hÞ ¼
3
Oðh
ðkÞ k¼1 Rk Gi .
ðkÞ
Since Gi
¼ Oð1Þ, we see that
3
Þ; u 2 C ðX
Oðh Þ; 2þa
P6
ð3:8Þ
Þ;
ð3:9Þ
Þ: u 2 C 2;a ðX
Þ, then from the Taylor expansion Similarly, if u 2 C 3;1 ðX 3 3 1X @2u 1 X @3u 4 ðkÞ uðpi Þ ¼ ruðpi Þ hk þ u pi ðhk Þi ðhk Þj ðpi Þ þ ðhk Þi ðhk Þj ðhk Þl ðp Þ þ Oðh Þ; 2 i;j¼1 @xi @xj 6 i;j;l¼1 @xi @xj @xl i
ð3:10Þ
we get
Wi ðu~ uÞ ¼
3 1X @2u 1 1X @3u 4 ^ a þ ðp Þ W a W C ðpi Þ þ Oðh Þ; i j i jkj i 2 j¼1 @x2j 3 6 jkj¼3 @xk11 @xk22 @xk33
ð3:11Þ
where k = (k1, k2, k3), jkj = k1 + k2 + k3, ki 2 {0, 1, 2, 3} and
C jkj ¼
6 X
ðkÞ
Gi
k1 k2 k3 ðkÞ ðkÞ ðkÞ pi pi pi pi pi pi : 1
k¼1
2
ð3:12Þ
3
ðkÞ
Here we note that if hik = h, k = 1, . . . , 6, then the Eq. (3.2) and the fact that Gi
C jkj ¼ 0;
is constant show
8jkj ¼ 3:
ð3:13Þ
Hence the Eqs. (3.8), (3.9), (3.11), (3.13) and Lemma 3.3 give the following lemma. ¼ ðu1 ; . . . ; um Þ be the solution of the Lemma 3.4 (Truncation error). Let u be the solution of the Dirichlet problem (1.1) and u ~ ¼ ðuðp1 Þ; . . . ; uðpm ÞÞ. Then for a nonuniform case hik – h for some k = 1, . . . , 6, the truncation error discrete system (Ph). Let u becomes 2
Wi ðu~ uÞ ¼ Oðh Þ: On the other hand, for the uniform case hik = h for all k = 1, . . . , 6, it becomes
8 2þa 2;a > > < Oðh Þ; if u 2 C ðXÞ; 3 3 Þ; ~u Þ ¼ Oðh Þ; if u 2 C ðX Wi ðu > > : 4 Þ: Oðh Þ; if u 2 C 3;1 ðX
4. Convergence analysis In this section,we will show the solvability of the discrete system (Ph) using a maximum principle of the difference operator Wi, and then analyze a superconvergence property for the approximate solution. To demonstrate that the discrete system (Ph) has a unique solution, we shall prove that the corresponding homogeneous system has only the trivial solution. ¼ ðv 1 ; . . . ; v m Þ 2 Rm . Lemma 4.1 (Maximum principle). Let us consider v satisfies (a) If v
Wi v P 0;
8i ¼ 1; . . . ; m0 ;
3632
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
then it holds
vi
0
¼ max v i ) i0 2 fm0 þ 1; . . . ; mg: 16i6m
(b) Alternatively, if v
Wi v 6 0;
8i ¼ 1; . . . ; m0 ;
then it holds
vi
0
¼ min
16i6m
v i ) i0 2 fm0 þ 1; . . . ; mg: ðkÞ
Proof. We prove part (a) by contradiction. Let us assume i0 2 {1, . . . , m0}. Then we review that the coefficients Gi in (2.12) < 0 since v iðkÞ v i0 < 0. This is a are always positive for sufficiently small mesh size. Hence the definition of i0 gives that Wi0 v 0 contradiction. We can prove the part (b) with a similar argument to the above. h Remark 4.2. Let us consider the homogeneous system corresponding to the discrete system (Ph), i.e., g(pi) = 0, i = m0 + 1, . . . , m. Then Lemma 4.1 shows that the maximum and minimum of the solution of this homogeneous system vanish. Hence the only solution is the trivial one. Thus it follows by the alternative principle for linear systems that the discrete system (Ph) has a unique solution for arbitrary boundary condition g. Now we observe a priori estimate of the difference operator Wi for the error analysis. For simplicity, we assume that there ðkÞ is a fixed positive integer n such that n < m0; for the indices i = 1, 2, . . . , n, pi 2 X and hik = h (" k = 1, 2 ,. . . , 6); for ðkÞ i = n + 1, . . . , m0, there is at least one index 1 6 k 6 6 with pi 2 @ X and hik – h. Consider the Dirichlet problem of
M/ ¼ 1 in X;
ð4:1Þ
/ ¼ 0 on C; ¼ ð/ ; . . . ; / Þ 2 Rm be the solution of the discrete problem (Ph) corresponding to it. Then Wi ð/Þ satisfies and let / 1 m
Wi ða^Þ
¼ Wi /
6
> 0;
i ¼ 1; 2; . . . ; m0 ;
ð4:2Þ
and hence Lemma 4.1(Maximum principle) shows
/j 6 0;
j ¼ 1; . . . ; m0 ;
and /j ¼ 0;
j ¼ m0 þ 1; . . . ; m:
ð4:3Þ
Using these facts, we can get the following a priori estimate for the difference operator Wi. ¼ ðv 1 ; . . . ; v m Þ 2 Rm , it holds that the difference operator Wi() satisfies the Theorem 4.3 (A priori estimate). For every vector v following inequality
jv j j 6
max jv i j þ K 1 j/j j max jWi v j þ K 2 max jWi v j;
m0 þ16i6m
16i6n
nþ16i6m0
j ¼ 1; . . . ; m;
¼ ð/1 ; . . . ; /m Þ 2 Rm is the solution of the discrete problem (Ph) corresponding to the problem (4.1) and where /
^j K 1 ¼ 6 min jWi a 16i6n
0
1 ;
B K 2 ¼ @ min
nþ16i6m0
X ðkÞ pi 2C
11 ðkÞ C Gi A :
ð4:4Þ
¼ ðb1 ; . . . ; bm Þ 2 Rm be a vector such that Proof. Let b
bj ¼
1; 0;
j ¼ 1; 2; . . . ; m0 ; j ¼ m0 þ 1; . . . ; m;
and þ N2 b ¼ v . Further, we let w ¼ v þ N1 / and v
N1 ¼ 6
max16i6n jWi v j ; ^ min16i6n Wi a
N2 ¼
maxnþ16i6m0 jWi v j : P ðkÞ minnþ16i6m0 pðkÞ 2C Gi i
From the facts of (4.3), bj = 0 = /j for any j = m0 + 1, . . . , m. Thus we see w j ¼ v j for any j = m0 + 1, . . . , m. Also Wi b becomes
8 > < 0; i ¼ 1; 2; . . . ; n; P ðkÞ Wi b ¼ Gi ; i ¼ n þ 1; . . . ; m0 : > : ðkÞ pi 2C
3633
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
Then the Eq. (4.2) shows that
P max jWj v j þ N1 Wi / ¼ max jWj v j 1 þ ¼ Wi v þ N1 Wi / Wi w 16j6n
16j6n
6Wi / P 0; ^ min16j6n Wj a
i ¼ 1; 2; . . . ; n
and
P max jW v j þ N1 W / þ N 2 Wi b ¼ Wi v þ N1 Wi / Wi w j i þ N 2 Wi b nþ16j6m0
0 B ¼ max jWj v j@ nþ16j6m0
Wi b minnþ16j6m0
P
ðkÞ
ðkÞ
pj 2C
Gj
1 C > 0; 1A þ N1 Wi /
i ¼ n þ 1; . . . ; m0 :
to obtain that for all j = 1, 2, . . . , m, Hence we can apply the maximum principle of Lemma 4.1 to w
v j 6 wj þ N1 j/j j N2 bj 6 max wi þ N1 j/j j þ N2 6 16i6m
¼
max wi þ N1 j/j j þ N2 ¼
m0 þ16i6m
max ðv i Þ þ j/j jK 1 max jWi v j þ K 2 max jWi v j ¼
m0 þ16i6m
16i6n
nþ16i6m0
max
m0 þ16i6m
v i þ N1 j/j j þ N2
max jv i j þ j/j jK 1 max jWi v j þ K 2 max jWi v j;
m0 þ16i6m
16i6n
nþ16i6m0
±
where denotes the jth component of w and Kj, j = 1, 2, are defined in (4.4). From the definition of v±, for all j = 1, 2, . . . , m, we have w j
v j 6 v j 6
max jv i j þ j/j jK 1 max jWi v j þ K 2 max jWi v j
m0 þ16i6m
16i6n
nþ16i6m0
which completes the proof. h ¼ ð/ ; . . . ; / Þ 2 Rm is the solution of the discrete problem (Ph) Corollary 4.4. Assume that X is a C2,a (0 < a 6 1) domain and let / 1 m satisfies corresponding to the problem (4.1). Then for sufficiently small mesh size h, /
/j ¼ OðhÞ;
j ¼ n þ 1; . . . ; m0 :
ð4:5Þ
Þ and / = 0 on C. Thus if we let / ~ ¼ ð/ðp Þ; . . . ; /ðp ÞÞ, Lemma 3.4 Proof. First note that the solution / of (4.1) is in C 2;a ðX 1 m shows
( ~ /Þ ¼ Wi ð/
2
Oðh Þ; 2þa
Oðh
hik – h; Þ;
hik ¼ h:
^ ¼ Oðh2 Þ, and also the definition of G defined in (2.10) shows that the coefOn the other hand, Lemma 3.1 shows that Wi a ðkÞ ðkÞ ðkÞ ficients Gi of ui ui defined in (2.12) have an asymptotic behavior Gi ¼ Oð1Þ. Thus the constants K1 and K2 defined in (4.4) have the asymptotic behaviors 2
K 1 ¼ Oðh Þ;
K 2 ¼ Oð1Þ:
Hence Theorem 4.3 shows that there are some constant Ck independent of h such that
2 ~ /Þ þ C 2 max Wi ð/ ~ /Þ þ j/ðp Þj j/j j 6 /j /ðpj Þ þ /ðpj Þ 6 C 1 h j/j j max Wi ð/ j 16i6n
a
2
¼ C 1 h j/j j þ C 2 h þ j/ðpj Þj;
nþ16i6m0
j ¼ 1; . . . ; m0 :
Þ and / = 0 on C, it shows that On the other hand, since the solution / of (4.1) is in CðX
/ðpj Þ ¼ OðhÞ;
j ¼ n þ 1; . . . ; m0 :
Hence for sufficiently enough small h, we can complete the proof. h From this Corollary, the above theorem can be restated as follows. 2,a Corollary 4.5. Assume that X is a C , 0 6 a 6 1, domain and let p1, . . . , pm be the node points of a quasi-uniform rectangle grid ^ ¼ jp1 j2 ; . . . ; jpm j2 . Then for every vector v ¼ ðv 1 ; . . . ; v m Þ 2 Rm , the difference operator Wi() satisfies the following and a
8 > 2 > > jv i j þ h max jWi v j þ max jWi v j ; < O m max 16i6n nþ16i6m0 0 þ16i6m jv j j ¼ > 1 > > : O max jv i j þ h max jWi v j þ max jWi v j ; m0 þ16i6m
16i6n
nþ16i6m0
j ¼ 1; . . . ; n; j ¼ n þ 1; . . . ; m0 :
Hence by Lemma 3.4 and the above Corollary 4.5, we can obtain the following superconvergence property for the discrete solution of (Ph).
3634
A.F.A. Dafa-Alla et al. / Applied Mathematics and Computation 217 (2010) 3624–3634
¼ ðu1 ; . . . ; um Þ be the solutions of (1.1) and (Ph), Theorem 4.6 (Superconvergence). Assume that X is a C2,a domain. Let u and u Þ, then respectively. Then the discrete problem (Ph) has the following convergence property. If u 2 C 3 ðX
(
juðpi Þ ui j ¼ while if u 2 C
2;a
Oðh Þ;
i ¼ n þ 1; . . . ; m0 ;
Þð0 < a 6 1Þ, ðX
juðpi Þ ui j ¼ Also, if u 2 C
i ¼ 1; 2; . . . ; n;
2
(
3;1
OðhÞ;
a
Oðh Þ; Oðh
i ¼ 1; 2; . . . ; n;
aþ1
Þ;
i ¼ n þ 1; . . . ; m0 ;
Þ, we get ðX 2
juðpi Þ ui j ¼ Oðh Þ;
i ¼ 1; 2; . . . ; m0 :
5. Conclusion We developed a new finite difference representation for solving the Dirichlet problem of Poisson’s equation in R3 based on a contour integration and the piecewise linear interpolation of the solution. We also analyze the convergence and show that it has a superconvergence property like two dimensional scheme. Even though we only discussed the derivation of a seven points difference formula on rectangle grid and its analysis, we think that it may be possible to derive a general n-points finite difference formula with n P 7 on a general shape of polyhedron, where each local subdomain consists of a triangular pyramid, because the scheme is based on the boundary element formulation which can be constructed on an arbitrary shape of polyhedron. The work presented here is primarily theoretical, and further numerical calculations need to be done to guarantee its efficiency. This will be the subject of further research. Acknowledgements The authors would like to express our sincere gratitude to reviewers for many helpful comments and valuable suggestions on the first draft of this paper. In particular crucial improvements of Section 4 and several other places are due to anonymous referee’s encouragements and advices. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]
R. Courant, D. Hilbert, Methods of Mathematical Physics, vol. 2, Interscience, New York, 1962. D. Gilbarg, N.S. Trudinger, Elliptic Partial Differential Equations of Second Order, Springer, Berlin, 1977. R. Li, Z. Chen, W. Wu, Generalized Difference Methods for Differential Equations, Marcel Dekker, Inc., 2000. M. Shashkov, Conservative Finite Difference Methods on General Grids, CRC Press, 1996. G.H. Shortley, R. Weller, Numerical solution of Laplace equation, J. Appl. Phys. 9 (1938) 334–348. N. Matsunaga, T. Yamamoto, Superconvergence of the Shortley–Weller approximation for Dirichlet problems, J. Comput. Appl. Math. 116 (2000) 263– 273. T. Yamamoto, On the accuracy of finite difference solutions for Dirichlet problems, RIMS Kokyuroku 1040, RIMS, Kyoto University, 1998, pp. 135–142. Q. Fang, Convergence of finite difference methods for convection diffusion problems with singular solution, J. Comp. Appl. Math. 152 (2003) 119–131. Z. Li, T. Yamamoto, Q. Fang, Superconvergence of solution derivatives for the Shortley–Weller difference approximation of Poisson’s equation. part I: smoothness problem, J. Comp. Appl. Math. 151 (2003) 221–307. Z. Li, H. Hu, Q. Fang, T. Yamamoto, Q. Fang, Superconvergence of solution derivatives for the Shortley–Weller difference approximation of Poisson’s equation. part II: singularity problems, Numer. Funct. Anal. Optim. 24 (2003) 195–221. Z. Li, C. Chien, H. Huang, Effective condition number for finite difference method, J. Comp. Appl. Math. 198 (2007) 208–235. T. Yamamoto, Q. Fang, X. Chen, Superconvergence and non-superconvergence of the Shortley–Weller approximations for Dirichlet problems, Numer Funct. Anal. Optim. 22 (2001) 455–470. K.E. Atkinson, The Numerical Solution of Integral Equations of the Second Kind, Cambridge Univ. Press, 1997. P.K. Banerjee, The Boundary Element Methods in Engineering, Mcgraw-Hill Book Company, 1994. C.A. Brebbia, J.C.F. Telles, L.C. Wrobel, Boundary Element Techniques–Theory and Applications in Engineering, Springer-Verlag, 1984. W. Mclean, Strongly Elliptic Systems and Boundary Integral Equations, Cambridge Univ. Press, 2000. P. Kim, H. Choi, S. Ahn, A new approach for the derivation of a discrete approximation formula on uniform grid for harmonic functions, Kyungpook Math. J. 47 (2007) 529–548. S. Ahn, S. Kim, P. Kim, Local boundary element based a new finite difference scheme for Poisson equations, Appl. Math. Comp. AMC-D-09-01762, Submitted for publication, 2009. W.C. Thacker, A. Conzalez, G.E. Putland, A method for automeshing the construction of irregular computational grids for Storm Suage forecase models, J. Comp. Phys. 37 (1980) 371–387. M. Ribeiro Filho, J.T. Pinho, Automatic and adaptive mesh generation for FEM with applications in microwave, SBMO/IEEE MTT-S IMOC Proceedings 17–20 (2001) 371–387. E.D. Lutz, Systematic derivation of contour integration formulae for Laplace and elastostatic gradient BIE’s, Comp. Mech. 14 (1994) 339–353.