Poisson 方程在不規則區域的高階緊緻差分法

全文

(1)國立交通大學. 應用數學系碩士論文. Poisson 方程在不規則區域的高階緊緻差分法 Higher-order compact difference scheme for the Poisson equation on irregular domains. 研究生：陳莉君指導老師：賴明治教授. 中華民國九十四年六月.

(2) Poisson 方程在不規則區域的高階緊緻差分法 Higher-order compact difference scheme for the Poisson equation on irregular domains. 研究生：陳莉君. Student：Li-Jun Chen. 指導教授：賴明治. Advisor：Ming-Chih Lai. 國立交通大學應用數學系碩士論文. A Thesis Submitted to Department of Applied Mathematics College of Science National Chiao Tung University in Partial Fulfillment of the Requirements for the Degree of Master in Applied Mathematics June 2005 Hsinchu, Taiwan, Republic of China. 中華民國九十四年六月.

(3) Poisson 方程在不規則區域的高階緊緻差分法. 學生：陳莉君. 國立交通大學. 指導教授：賴明治. 應用數學. 摘. 學系﹙研究所﹚碩士班. 要. 在這篇論文中，我們考慮在不規則區域上有著 Dirichlet 邊界條件的二維 Poisson 方程，而此不規則區域將被限制在一個矩形區域中來做計算。在原不規則區域內的網格點，我們的高階精度方法是用標準的緊緻九點格式來離散此 Poisson 方程，但對於靠近邊界的點將要做些特別處理。在這些需要另外做處理的點，我們利用外插法來造出仿真的值。使用的外插法有常數、線性以及二次的外插法，並將分別可以得到一階、二階以及三階的精度。其中常數及線性的外插法，可以使要解的線性系統仍保持對稱性。. i.

(4) Higher-order compact difference scheme for the Poisson equation on irregular domains student：Li-Jun Chen. Advisors：Dr. Ming-Chih Lai. Department﹙Institute﹚of Applied Mathematics National Chiao Tung University. ABSTRACT. In this thesis, we consider the 2D Poisson equation subject to Dirichlet boundary conditions on an irregular domain. The region of interest is embedded in a rectangular domain. For our higher-order accurate scheme, at internal grid points, the Poisson equation is discretized with the standard compact nine point stencil with special treatment at the edges. At the irregular point, we define ghost value constructed by extrapolations. This yields first, second and third order accuracy in the case of the constant, linear and quadratic extrapolations, respectively. In the case of constant and linear extrapolations, the linear system is symmetric.. ii.

(5) 誌. 謝. 這篇論文的完成必須感謝許多協助與支持我的人。首先，感謝我的指導教授賴明治老師這兩年來的指導與勉勵，在學問與待人處事方面，都讓我受惠良多，謹此致上我最誠摯的敬意與謝意。口試期間，承蒙葉立明老師、黃聰明老師及吳金典老師費心審閱並提供許多寶貴之意見，使本論文之完稿得以更加齊備，永誌於心。研究所求學過程中，感謝昱豪學長總在我遇到問題的時候，給予我意見以及耐心的指導；接著感謝同窗好友雅靜與研究所同學們的關心與協助，這些日子裡的互相幫忙、互相砥礪，深厚的情誼永難忘懷。此外，謝謝球隊的舉卿學姊、雅文學姊、曲敏、乃心以及學妹們，因為有你們碩士班兩年生涯能夠如此充實愉快的度過，與你們一起的點點滴滴，將是我永遠的回憶。最後，要感謝我最敬愛的父母，感謝你們在我學習的路上給予我最大的支持以及自由的學習環境，讓我可以無後顧之憂的學習。再次感謝所有幫助過我、關心我的人，謝謝你們！. iii.

(6) 目中文提要英文提要誌謝目錄 1. 2. 2.1 2.2 3. 4. 5. Reference. 錄. ……………………………………………………………… ……………………………………………………………… ……………………………………………………………… ……………………………………………………………… Introduction …………………………………………… One-dimensional Poisson equation ………………… Dirichlet boundary condition ……………………… Neumann boundary condition ………………………… Two-dimensional Poisson equation ………………… Examples ………………………………………………… Conclusion ……………………………………………… ………………………………………………………………. iv. i ii iii iv 1 2 3 6 9 12 20 21.

(7) 1. Introduction. In this thesis, we consider the solution of the Poisson equation on an irregular domain, subject to Dirichlet boundary conditions. The Poisson equation subject to Dirichlet boundary conditions on an irregular domain can be treated by embedding the region in a rectangular domain and solving by finite differences over the domain. The crucial issue is the discretization of the boundaries of the irregular domain. There are many other approaches to this problem in the literature. In [4], the authors solved a variable coefficient Poisson equation in the presence of an irregular interface where Dirichlet boundary conditions were imposed. They used a finite volume method that results in a non-symmetric discretization matrix. Both multigrid methods and adaptive mesh refinement were used. In [5], this non-symmetric discretization was coupled to a volume of fluid front tracking method in order to solve Stefan problem. In [9], the basic idea of the ghost fluid method [7] was employed to develop a first-order-accurate symmetric finite difference scheme based on the Cartesian grid to solve a variable Poisson equation in the presence of an irregular interface. Subsequently in [1], the approach in [9] was modified to obtain a second order accurate symmetric finite difference scheme based on the Cartesian grid to solve a variable Poisson equation with a Dirichlet boundary condition. The modification used the signed distance level set function to obtain a linear interpolation from the boundary value and the solution values in coordinate-wise directions to determine the ghost fluid values. The intention of this paper is to extend the idea of [2,11]. In [2], the authors exploit the methodology of [1] to derive a fourth order accurate finite difference discretization for the Laplace equation on irregular domains. But in order to guarantee a fourth order accurate, the difference scheme that the authors used was a standard long stencil. The primary objective of this paper is to keep higher-order accurate and to make the scheme to be a compact one. The rest of the paper is organized as follows: In Section 2 we deal with the 1-D Poisson equation with Dirichlet boundary conditions and try to solve Neumann boundary problem as well. In Section 3 we extend the methodology discussed in Section 2 to two spatial dimensions. Numerical examples are presented in Section 4 before we conclude with a summary in Section 5.. 1.

(8) 2. One-dimensional Poisson equation. We consider a Cartesian computational domain, Ω = [a, b], with a lower dimensional interface, Γ, that divides the computational domain into disjoint pieces, Ω− and Ω+ . The 1-D Poisson equation is given by Txx = f, x ∈ Ω− = [a, xI ].. (2.1). A uniform grid is taken over [a, b] . Dirichlet boundary conditions or Neumann boundary conditions are assumed given at two boundary points x = a and xI , xI typically is not grid point. In the other subdomain we set T = 0 , so that we have ½ Txx = f x ∈ Ω− . (2.2) T =0 x ∈ Ω+ In general there is a discontinuity at xI . The solution to the Poisson equation is computed at the grid points and is written as Ti = T (xi ). We consider the fourth order discretization : Ti+1 − 2Ti + Ti−1 ∆x2 4 + O(∆x ) = (T ) + (Txxxx )i xx i ∆x2 12 Ti+1 − 2Ti + Ti−1 fi+1 + 10fi + fi−1 ≈ (2.3) 2 ∆x 12 For each unknown, Ti , Eq(2.3) is used to fill in one row of a matrix creating a linear system of equations. This discretization is valid if all the node values belong to the same domain, but needs to be modified otherwise. For example, suppose the interface location, xI is located in between the nodes xi and xi+1 (see Fig.1) and suppose that we seek to write the equation satisfied by Ti . Since the solution is not defined across the interface, we need valid values for Ti+1 and fi+1 that emulate the behavior of the solution defined to the left of the interface. We achieve this G G by defining ghost values Ti+1 and fi+1 constructed by extrapolating the value of T across the interface. The discretization for such points in the neighborhood of the interface is rewritten as G G + 10fi + fi−1 fi+1 − 2Ti + Ti−1 Ti+1 ≈ 2 ∆x 12. (2.4). In the remainder of this section, we describe how to construct the ghost values G more precisely. and fi+1. G Ti+1. 2.

(9) Solution Profile. Ti. TI. xi. xi 1. xI T'x. :. . Interface. :. Position. . Ti G1. Figure 1: Definition of the ghost cells with linear extrapolation. 2.1. Dirichlet boundary condition. In this section, we consider the situation that Dirichlet boundary conditions are given at two boundary points, T (a) = Ta and T (xI ) = TI . We first construct an interpolant Te(x) of T (x) on the left of the interface, such that Te(0) = Ti , and then G we define Ti+1 = Te(∆x). Fig.1 illustrates the definition of the ghost cell in the case of the linear extrapolation. We consider constant, linear, quadratic and cubic extrapolations defined by: • Constant extrapolation Take Te(x) = d with: 1. d = TI • Linear extrapolation Take Te(x) = cx + d with: 1. Te(0) = Ti 2. Te(θ∆x) = TI. 3.

(10) • Quadratic extrapolation Take Te(x) = bx2 + cx + d with: 1. Te(−∆x) = Ti−1 2. Te(0) = Ti 3. Te(θ∆x) = TI • Cubic extrapolation Take Te(x) = ax3 + bx2 + cx + d with: 1. Te(−∆x) = Ti−1 2. Te(0) = Ti 3. Te(θ∆x) = TI 4. Te00 (θ∆x) = fI Then we get G Ti+1 = TI. G Ti+1. 1 1 G Ti+1 = (1 − )Ti + TI θ θ 1−θ 2(θ − 1) 2 = Ti−1 + Ti + TI 1+θ θ θ(1 + θ). (2.5) (2.6) (2.7). and G Ti+1 =. (1 − θ)(2θ − 1) 4(θ − 1) 6 1−θ 2 Ti−1 + Ti + TI + h fI (2.8) (2θ + 1)(θ + 1) 2θ + 1 (2θ + 1)(θ + 1) 2θ + 1. which are defined by constant, linear, quadratic and cubic extrapolations, respectively. Similarly, we construct an interpolant fe(x) of f (x). The definitions of constant, linear and quadratic extrapolations are the same as Te. So, we have G fi+1 = fI. 1 1 G = (1 − )fi + fI fi+1 θ θ and G = fi+1. 1−θ 2(θ − 1) 2 fi−1 + fi + fI 1+θ θ θ(1 + θ). (2.9) (2.10) (2.11). which are defined by constant, linear and quadratic extrapolations, respectively. But there is a little different to the cubic extrapolation :. 4.

(11) • Cubic extrapolation Take fe(x) = ax3 + bx2 + cx + d with: 1. fe(−2∆x) = fi−2 2. fe(−∆x) = fi−1 3. fe(0) = fi 4. fe(θ∆x) = fI then we have G fi+1 =. θ−1 (3 − 3θ) 3θ − 3 6 fi−2 + fi−1 + fi + 3 fI θ+2 θ+1 θ θ + 3θ2 + 2θ. (2.12). which is defined by cubic extrapolation. In these equations θ = (xI − xi )/∆x refers to the cell fraction occupied by the subdomain Ω− . This yields first, second, third and fourth order accuracy in the case of the constant, linear, quadratic and cubic extrapolations, respectively. If we were solving for the domain Ω+ , the equation satisfied by Ti+1 and fi+1 requires the definition of the ghost cells TiG and fiG . In this case, we write TiG = Te(∆x) and fiG = fe(∆x) with the definition for Te modified as follows: θ = (xi+1 − xI )/∆x, Ti is replaced by Ti+1 , fi is replaced by fi+1 , Ti−1 is replaced by Ti+2 and fi−1 is replaced by fi+2 . We note that the construction of Te and fe cannot be arbitrary. It is obviously limited by the number of points within the domain, but also by how close the interface from a grid node. The latter restriction comes from the fact that, as θ → 0, the behavior of the interpolant deteriorates. In the case of the constant extrapolation, the corresponding matrix   −2 1  1 −2 1      . . . . . . AConstant =   . . .    1 −2 0  → ith row → (i + 1)th row 0 1 is symmetric and diagonally dominant. In the case of the linear extrapolation, the corresponding matrix   −2 1   1 −2 1     .. .. .. ALinear =   . . .   1  1 −1 − θ 0  → ith row → (i + 1)th row 0 1. 5.

(12) is symmetric and diagonally dominant. This allows for the use of fast iterative solvers such as preconditioned conjugate gradient. But in the case of the quadratic and cubic extrapolations, the corresponding matrices   −2 1  1 −2  1     . . . . . . AQuadratic =  . . .    2(θ+1) 1−θ  1 + 1+θ −2 + θ 0  → ith row → (i + 1)th row 0 1 and . ACubic. −2 1  1 −2 1   .. .. .. . . . =   (1−θ)(2θ−1) 4(θ−1)  1 + (2θ+1)(θ+1) −2 + 2θ+1 0 0 1.        → ith row → (i + 1)th row. are non-symmetric. So the non-symmetric linear system is solved with a BiCGSTAB (see e.g. [10] ) using an incomplete LU factorization for the preconditioner. For the linear solver that require an initial guess, setting all Ti identically zero is usually sufficient.. 2.2. Neumann boundary condition. We take little effort to replace the Dirichlet boundary conditions with Neumann boundary conditions, reformulating the Poisson equation as ½ Txx = f, x ∈ [a, xI ] . (2.13) Tx (a) = α, Tx (xI ) = β So, we needs to do some modifications of the extrapolations Te(x). We discuss linear, quadratic and cubic extrapolations defined by: • Linear extrapolation Take Te(x) = cx + d with: 1. Te(0) = Ti 2. Te0 (θ∆x) = β • Quadratic extrapolation Take Te(x) = bx2 + cx + d with:. 6.

(13) 1. Te(−∆x) = Ti−1 2. Te(0) = Ti 3. Te0 (θ∆x) = β • Cubic extrapolation Take Te(x) = ax3 + bx2 + cx + d with: 1. Te(−∆x) = Ti−1 2. Te(0) = Ti 3. Te0 (θ∆x) = β 4. Te00 (θ∆x) = fI The discussion above leads naturally to G Ti+1 = Ti + β∆x G Ti+1 =. (2.14). 1 − 2θ 4θ 2β∆x Ti−1 + Ti + 1 + 2θ 1 + 2θ 1 + 2θ. and. (2.15). −3θ2 + 3θ − 1 6θ2 + 2 6θ T + Ti + 2 β∆x i−1 2 2 3θ + 3θ + 1 3θ + 3θ + 1 3θ + 3θ + 1 −3θ2 + 1 + 2 fI ∆x2 (2.16) 3θ + 3θ + 1 which are defined by linear, quadratic and cubic extrapolations, respectively. And the coefficient matrices of the linear system are   −2 1  1 −2 1      . . . . . . ALinear =   . . .    1 −1 0  → ith row → (i + 1)th row 0 1 G Ti+1 =. . AQuadratic. −2 1  1 −2 1   . . .. .. .. =  .  4θ 1−2θ  0 1 + 1+2θ −2 + 1+2θ 0 1. 7.        → ith row → (i + 1)th row.

(14) and . ACubic. −2 1  1 −2 1   .. .. .. . . . =   2 +2) 2 +3θ−1 6θ −3θ  1 + 3θ2 +3θ+1 −2 + 3θ2 +3θ+1 0 0 1.        → ith row → (i + 1)th row. which are obtained by linear, quadratic and cubic extrapolations, respectively.The extrapolations fe of f don’t have any modifications. Because we have the value of f on the boundary exactly. This yields first, second and third order accuracy in the case of the linear, quadratic and cubic extrapolations, respectively. But this method has some shortcomings. It is unable to expand this method to solve two-dimensional Poisson equation with Neumann boundary conditions. Because we do not have enough information on the boundary to use this method. The overall accuracy for T and the nature of the resulting linear system is determined by the degree of the interpolation function Te, which is summarized in Tab. 1.. Dirichlet boundary condition Degree of extrapolation Constant Linear Quadratic Cubic. Order of accuracy Linear system First Symmetric Second Symmetric Third Non-symmetric Fourth Non-symmetric. Neumann boundary condition Degree of extrapolation Linear Quadratic Cubic. Order of accuracy Linear system First Symmetric Second Non-symmetric Third Non-symmetric. Table 1: Order of accuracy and nature of the linear system corresponding to the constant, linear, quadratic and cubic case. 8.

(15) :. T. 0. 'T. f. :. Figure 2: An irregular interface Γ dividing the domain Ω into two subdomain Ω− and Ω+ .. 3. Two-dimensional Poisson equation Consider the two spatial dimension Poisson equation ∆T = f (x, y). and let Ω− be any irregular 2-D shape inscribed within a rectangle with boundary Γ at which Dirichlet conditions T (x, y) = g(x, y) are specified.We can regard the boundary Γ as a interface that divies the domain Ω into two disjoint pieces, Ω− and Ω+ (see Fig.2). As T = 0 outside the physical domain, there may be jumps on Γ. We use the standard compact nine point stencil scheme, the 9-point Laplacian, denote by ∆9 Ti,j , ∆9 Ti,j =. 1 [4Ti−1,j +4Ti+1,j +4Ti,j−1 +4Ti,j+1 +Ti−1,j−1 +Ti−1,j+1 +Ti+1,j−1 +Ti+1,j+1 −20Ti,j ] 6h2. If we apply this to the true solution and expand in Taylor series we find that ∆9 Ti,j = ∆Ti,j +. h2 [(Txxxx )i,j + 2(Txxyy )i,j + (Tyyyy )i,j ] + O(h4 ). 12. The additional terms lead to a very nice form for the dominant error term, since Txxxx + 2Txxyy + Tyyyy = ∆(∆T ).. 9.

(16) P1. ( xi 1 , y j 1 ). LI. P0. :. ( xi , y j ). L. :. . Figure 3: A diagram of an irregular grid point P0 .. This is the Laplacian of Laplacian of T which is known as the biharmonic. Because we are solving ∆T = f , then we have Txxxx + 2Txxyy + Tyyyy = ∆f. Hence we can compute the dominant term in the truncation error easily from the known function f without knowing the true solution T to the problem. So,we can obtain a fourth-order accurate method of the form ∆9 Ti,j ≈. fi−1,j + fi+1,j + fi,j−1 + fi,j+1 + 8fi,j . 12. (3.17). This discretization is valid if all the node values belong to the same domain, but needs to be modified otherwise. The methodology discussed in Section 2.1 extends naturally to two spatial dimensions. Before we modify the Eq.(3.17), we define a grid node (xi , yj ) as regular if all neighbouring nodes are on the same side of the interface.On contrary, a grid node (xi , yj ) as irregular if at least one adjacent node is on the other side of the interface. In order to illustrate our methodology, we suppose. 10.

(17) Ti G1, j. TI , j Tx. Ti 1, j. Ti , j xi xI 'x. Ti G1, j 1 TL. TL. LI P0 P1 P0. Ti , j Ti 1, j 1. G G Figure 4: The ghost values Ti−1,j and Ti−1,j+1 which were constructed along xdirection and the line segment L, respectively.. P0 = (xi , yj ) is a irregular grid node and the distribution of its neighbouring nodes were displayed in Fig.3. The discretization for the irregular point P0 = (xi , yj ) is then written as 1 G G [4T G +4Ti+1,j +4Ti,j−1 +4Ti,j+1 +Ti−1,j−1 +Ti−1,j+1 +Ti+1,j−1 +Ti+1,j+1 −20Ti,j ] 6h2 i−1,j G G fi−1,j + fi+1,j + fi,j−1 + fi,j+1 + 8fi,j . 12 G G About the ghost values Ti−1,j and fi−1,j ,we consider the left arm of the stencil, i.e. the line segment connecting (xi−1 , yj ) and (xi , yj ). We first find the interface location (xI , yj ) that is the intersection point of the interface and the line segment. In order to find (xI , yj ), we solve a nonlinear equation. In our example section we −xI . So use Newton’s Method to solve nonlinear equation. Then we define θx = xi∆x x x we can construct a constant, linear, or quadratic extrapolation Te and fe of T and G G are similar. and fi,j+1 f in the x-direction. The procedure to find Ti,j+1 G The remainder ghost value Ti−1,j+1 , we consider the line segment, L, connecting (xi−1 , yj+1 ) and (xi , yj ). Because the line segment and the interface are crossing in I −P0 | a point, so we have the interface locate LI (see Fig.3). Then we define θL = |L . |P −P1 | Finally, we can construct a constant, linear, or quadratic extrapolation TeL of T. ≈. 11.

(18) 4. 3.5. 3. 2.5. 2. 1.5. 1. 0.5. 0. 0. 0.1. 0.2. 0.3. 0.4. 0.5. 0.6. 0.7. Figure 5: 1D Poisson equation, Txx = f ,on [0,0.5] with Dirichlet boundary conditions. The exact solution is T = x7 − x3 + 12x2 − 2.5x + 2. The grid size is 64 and the ghost cells are defined by quadratic extrapolation.. along the line segment L(see Fig.4). Note that on irregular domains, the number of available grid nodes within the domain might limit the extrapolation to a lower degree for some grid resolution.. 4. Examples. We test our methodology on the following examples. In the case where the linear system is symmetric, we use a preconditioned conjugate gradient method with an incomplete Cholesky preconditioner. In the case where the linear system is non-symmetric, we use the BICGSTAB method. The order of the scheme is given as ¯ ¯ ¯ log(kErrorn k∞ /kError2n k∞ ) ¯ ¯ ¯. order = ¯ ¯ log(2). 4.1. Example 1. Consider Txx = f on Ω = [0, 0.5] with an exact solution of T = x7 − x3 + 12x2 − 2.5x + 2. The computational domain is discretized into cell of size ∆x where the. 12.

(19) 0.4. 0.35. 0.3. 0.25. 0.2. 0.15. 0.1. 0.05. 0. −0.05. 0. 0.1. 0.2. 0.3. 0.4. 0.5. 0.6. 0.7. Figure 6: 1D Poisson equation, Txx = f ,on [0,0.5] with Neumann boundary conditions. The exact solution is T = 4x2 sin(2πx). The grid size is 64 and the ghost cells are defined by cubic extrapolation.. cell centers are referred to as grid nodes. The Dirichlet boundary conditions are specified. Tab.2 shows the results of the numerical accuracy test and the ghost cells are defined by constant, linear, quadratic and cubic extrapolations. Fig.5 shows the numerical solution with 64 grid points and the ghost cells are defined by quadratic extrapolation.. 4.2. Example 2. Consider Txx = f on Ω = [0, 0.5] with an exact solution of T = 4x2 sin(2πx). The computational domain is discretized into cell of size ∆x where the cell centers are referred to as grid nodes. The boundary conditions are given as T (0) = α and Tx (0.5) = β. Tab.3 shows the results of the numerical accuracy test and the ghost cells are defined by linear, quadratic and cubic extrapolations. Fig.6 shows the numerical solution with 64 grid points and the ghost cells are defined by cubic extrapolation.. 13.

(20) Constant extrapolation Number of points L∞ - error 16 0.13531 32 0.06839 64 3.4395E-02 128 1.7200E-02 256 8.6382E-03. Order 0.984 0.992 1.000 0.994. Linear extrapolation Number of points L∞ - error 16 2.9232E-03 32 7.3162E-04 64 1.8300E-04 128 4.5764E-05 256 1.1443E-05. Order 1.998 1.999 2.000 2.000. Quadratic extrapolation Number of points L∞ - error 16 1.1525E-05 32 1.5525E-06 64 2.0306E-07 128 2.5958E-08 256 3.2811E-09. Order 2.892 2.935 2.968 2.984. Cubic extrapolation Number of points L∞ - error 16 3.2479E-07 32 2.0668E-08 64 1.2974E-09 128 8.1219E-11 256 5.2820E-12. Order 3.974 3.994 3.998 3.943. Table 2: The results of the numerical accuracy test and the ghost cells are defined by constant, linear, quadratic and cuic extrapolations.. 14.

(21) Linear extrapolation Number of points L∞ - error 16 8.2213E-04 32 3.4943E-04 64 1.0548E-04 128 2.8641E-05 256 7.4445E-06. Order 1.234 1.728 1.881 1.944. Quadratic extrapolation Number of points L∞ - error 16 1.9985E-03 32 5.0238E-04 64 1.2491E-04 128 3.1088E-05 256 7.7515E-06. Order 1.992 2.008 2.006 2.004. Cubic extrapolation Number of points L∞ - error 16 1.1580E-03 32 1.4326E-04 64 1.7887E-05 128 2.2373E-06 256 2.7984E-07. Order 3.015 3.002 2.999 2.999. Table 3: The results of the numerical accuracy test and the ghost cells are defined by linear, quadratic and cubic extrapolations.. 15.

(22) 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1.5 1 0.5 0 −0.5 −1 −1.5. −1.5. −0.5. −1. 0. 0.5. 1. 1.5. Figure 7: Solution of the Poisson equation on the unit circle. The exact solution is T = cos(x + y). The grid size is 64 × 64 and the ghost cells are defined by quadratic extrapolation.. 4.3. Example 3. Consider the Poisson equation ∆T = −2 cos(x + y) on the unit circle with Dirichlet boundary conditions. The exact solution is T = cos(x + y). The domain is embedded in a square. Outside the unit circle we set T = 0. Tab.4 shows the results of the numerical accuracy test and the ghost cells are defined by constant, linear, and quadratic extrapolations. Fig.7 depicts the solution on a 64 × 64 grid and the ghost cells are defined by quadratic extrapolation.. 4.4. Example 4. Consider the Poisson equation ∆T = −π 2 (sin(πx)+sin(πy)+cos(πx)+cos(πy)+ 30x4 + 30y 4 ) on an irregular domain Ω− with Dirichlet boundary conditions. The exact solution is T = sin(πx) + sin(πy) + cos(πx) + cos(πy) + x6 + y 6 . The domain is embedded in a square Ω. So we can regard the boundary of Ω− as a interface that divies Ω into two disjoint pieces, Ω− and Ω+ . Outside the interface we set T = 0.. 16.

(23) Constant extrapolation Number of points L∞ - error 16 6.6055E-02 32 4.4915E-02 64 2.3703E-02 128 1.2745E-02. Order 0.556 0.922 0.895. Linear extrapolation Number of points L∞ - error 16 1.5867E-03 32 3.3722E-04 64 9.6290E-05 128 2.2966E-05. Order 2.234 1.808 2.068. Quadratic extrapolation Number of points L∞ - error Order 16 3.2722E-04 32 2.9732E-05 3.46017 64 4.2617E-06 2.80252 128 4.8286E-07 3.14175 Table 4: The results of the numerical accuracy test and the ghost cells are defined by constant, linear, and quadratic extrapolations.. 17.

(24) 3. 2. 1. 0. −1. −2 0.5 0 −0.5 −0.8. −0.6. −0.4. −0.2. 0. 0.2. 0.4. 0.6. Figure 8: Solution of the Poisson equation on an irregular domain in two spatial dimensions. The exact solution is T = sin(πx)+sin(πy)+cos(πx)+cos(πy)+x6 +y 6 . The grid size is 64 × 64 and the ghost cells are defined by quadratic extrapolation.. The interface is parameterized by (x(α), y(α)),where √ ½ x(α) = 0.02√ 5 + (0.5 + 0.2 sin(5α)) cos(α), y(α) = 0.02 5 + (0.5 + 0.2 sin(5α)) sin(α),. (4.18). with α ∈ [0, 2π]. Tab.5 shows the results of the numerical accuracy test and the ghost cells are defined by constant, linear, and quadratic extrapolations. Fig.8 depicts the solution on a 64 × 64 grid and the ghost cells are defined by quadratic extrapolation. Note that on irregular domains, the number of available grid nodes within the domain might limit the extrapolation to a lower degree for some grid resolution. It is difficult to solve the 2D Poisson equation with Neumann boundary condition on irregular domain. Because we have no idea of the values of Txx and Tyy on the boundary. So we can’t extend the methodology discussed in Section 2.2 to two spatial dimensions.. 18.

(25) Constant extrapolation Number of points L∞ - error 16 0.2176 32 0.1306 64 6.5600E-02 128 3.4000E-02. Order 0.737 0.993 0.948. Linear extrapolation Number of points L∞ - error 16 7.2809E-02 32 3.4237E-03 64 8.2815E-04 128 2.0285E-04. Order 4.410 2.048 2.030. Quadratic extrapolation Number of points L∞ - error 16 6.4135E-02 32 3.4870E-04 64 5.1107E-05 128 6.1969E-06. Order 7.523 2.770 3.044. Table 5: The results of the numerical accuracy test and the ghost cells are defined by constant, linear, and quadratic extrapolations.. 19.

(26) 5. Conclusions. We have proposed a simple finite difference algorithm for obtaining higher-order accurate solutions for the Poisson equation subject to Dirichlet boundary conditions on irregular domains. The crucial issue is the discretization of the boundaries of the irregular domain. At the irregular point, we define ghost value constructed by extrapolation. In 1D Poisson equation with Dirichlet boundary condition problem, we get first, second, third and fourth order accuracy in the case of the constant, linear, quadratic and cubic extrapolations, respectively. In 1D Poisson equation with Neumann boundary condition problem, we get first, second and third order accuracy in the case of the linear and quadratic, and cubic extrapolations, respectively. In 2D Poisson equation with Dirichlet boundary condition problem, we get first, second and third order accuracy in the case of the constant, linear and quadratic extrapolations, respectively. And except the quadratic and cubic extrapolations, the linear system is symmetric.. 20.

(27) References [1] F. Gibou, R.P. Fedkiw, L.T Cheng, M Kang, A second-order accurate symmetric discretization of the Poisson equation on irregular domains. J. Comput. Phys. 176 (2002) 205–227. [2] Fr´ ed´ eric Gibou, Ranald Fedkiw. A fourth order accurate discretization for the Laplace and Heat equations on arbitrary domains, with applications to the Stenfan problem. J. Comput. Phys. 202 (2005) 577–601. [3] G. Golub and C. Van Loan, Matrix Computations. Johns Hopkins Univ. Press, Baltimore, 1989. [4] H. Johansen, P. Colella, A Cartesian grid embedded boundary method for Poisson’s equation on irregular domains. J. Comput. Phys. 147 (1998) 60–85. [5] H. Johansen, P. Colella, Cartesian grid embedded boundary finite difference methods for elliptic and parabolic differential equations on irregular domains. Ph. D. Thesis, University of California, Berkeley, CA, 1997. [6] Petter Andreas Berthelsen, A decomposed immersed interface method for variable coefficient elliptic equations with non-smooth and discontinuous solutions. J. Comput. Phys. 197 (2004) 364–386. [7] R. Fedkiw, T. Aslam, B Merriman, S. Osher ,A Non-oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid Method). J. Comput. Phys. 152 (1999) 457–492. [8] R. LeVeque, Z. Li, The immersed interface method for elliptic equations with discontinuous coefficients and singular sources. SIAM J. Numer. Anal 31 (1994) 1019–1044. [9] Xu-Dong Liu, Ronald P. Fedkiw, and Myungjoo Kang, A boundary condition capturing mtheod for Poisson’s equation on irregular domains. J. Comput. Phys. 160 (2000) 151–178. [10] Y. Saad, Iterative Methods for Spares Linear System, PWS. Publishing Company, Boston, 1996. [11] Z. Jomaa, C. Macaskill, The embedded finite difference method for the Poisson equation in a domain with an irregular boundary and Dirichlet boundary conditions. J. Comput. Phys. 202 (2005) 488–506.. 21.

(28)