High order two-grid finite difference methods for interface and internal layer problems

While it is possible to have fourth order schemes for regular one dimensional problems, see for example, Chapter 7 in [16], it is difficult to obtain fourth order accurate schemes in two and three dimensions for discontinuous coefficients and curved interfaces. One of the purposes of the two-grid method is to provide an alternative approach that can obtain comparable high order accuracy near the interface. We use a second order method near the interface with a finer mesh but a fourth order scheme away from the interface.

Assume that the interface $\alpha$ is between two grid points $x_{j}$ and $x_{j+1}$, that is, $\leq\alpha<x_{j+1}$. The interface is one of the grid points when $x_{j}=\alpha$, see Figure 2 for an illustration. The finite difference equations at two irregular grid points $x_{j}$ and $x_{j+1}$ have the following forms, see [16],

$$\begin{array}[]{rcl}&\displaystyle\gamma_{j,1}U_{j-1}+\gamma_{j,2}U_{j}+\gamma_{j,3}U_{j+1}+KU_{j}=f_{j}+C_{j},\\ \vskip 5.0pt\cr&\displaystyle\gamma_{j+1,1}U_{j}+\gamma_{j+1,2}U_{j+1}+\gamma_{j+1,3}U_{j+2}+KU_{j+1}=f_{j+1}+C_{j+1}.\end{array}$$

(8)

The coefficients of the finite difference equations for the special case $\bar{C}=0$ have the following closed form:

$$\begin{array}[]{rcl}&&\left\{\begin{array}[]{l}\displaystyle\gamma_{j,1}=(\kappa^{-}-[\kappa]{(x_{j}-\alpha)}/{h_{f}})/D_{j},\\ \vskip 5.0pt\cr\displaystyle\gamma_{j,2}=(-2\kappa^{-}+[\kappa]{(x_{j-1}-\alpha)}/{h_{f}})/D_{j},\\ \vskip 5.0pt\cr\displaystyle\gamma_{j,3}={\kappa^{+}}/D_{j},\end{array}\right.\\ \vskip 5.0pt\cr\vskip 5.0pt\cr&&\left\{\begin{array}[]{l}\displaystyle\gamma_{j+1,1}={\kappa^{-}}/D_{j+1},\\ \vskip 5.0pt\cr\displaystyle\gamma_{j+1,2}=(-2\kappa^{+}+[\kappa]{(x_{j+2}-\alpha)}/h_{f})/D_{j+1},\\ \vskip 5.0pt\cr\displaystyle\gamma_{j+1,3}=(\kappa^{+}-[\kappa]{(x_{j+1}-\alpha)}{}/h_{f})/D_{j+1},\end{array}\right.\end{array}$$

(9)

where

$$\begin{array}[]{rcl}D_{j}&=&\displaystyle h_{f}^{2}+[\kappa](x_{j-1}-\alpha)(x_{j}-\alpha)/2\kappa^{-},\\ \vskip 5.0pt\cr D_{j+1}&=&\displaystyle h_{f}^{2}-[\kappa](x_{j+2}-\alpha)(x_{j+1}-\alpha)/2\kappa^{+}.\end{array}$$

(10)

The correction terms are given by

$$C_{j}=\gamma_{j,3}\,(x_{j+1}-\alpha)\,{\frac{C}{\kappa^{+}}},\quad\quad C_{j+1}=\gamma_{j+1,1}\,(\alpha-x_{j})\,{\frac{C}{\kappa^{-}}}.$$

(11)

The stability and consistency of the finite difference scheme are discussed in [16].

One of challenges in two-grid methods is to design a compatible high order finite difference discretization at grid points that border the coarse and fine grids. In this subsection, we assume that $\kappa$ and $K$ are constants. Let $x_{i}$ be such a grid. Assume that $x_{i}-x_{i-1}=h_{1}$ and $x_{i+1}-x_{i}=h_{2}$. We design the following compact finite difference scheme,

$\displaystyle\sum_{k=-1}^{1}\alpha_{k}U_{i+k}=\sum_{k=-1}^{1}\beta_{k}f_{i+k},\quad\quad\sum_{k=-1}^{1}\beta_{k}=1.$

(12)

There are six unknowns, $\alpha_{k}$ and $\beta_{k}$, $k=-1,0,\,1$. Define the ‘local truncation error’ as

$\displaystyle T_{i}=\sum_{k=-1}^{1}\alpha_{k}u(x_{i+k})-\sum_{k=-1}^{1}\beta_{k}f(x_{i+k}),$

(13)

assuming that $u(x)$ is the true solution. We want the finite difference scheme to be exact for all fourth order polynomials if possible. To achieve this, we expand the right hand side above at $x_{i}$ to get

$\displaystyle T_{i}=\sum_{k=-1}^{1}\alpha_{k}\sum_{l=0}^{4}\frac{u^{(l)}(x_{i})}{l!}\hat{h}_{k}^{l}-\sum_{k=-1}^{1}\beta_{k}\sum_{l=0}^{2}\frac{f^{(l)}(x_{i})}{l!}\hat{h}_{k}^{l}+O\left(\|{\mbox{\boldmath$\alpha$}}\|_{\infty}h^{5}+\|{\mbox{\boldmath$\beta$}}\|_{\infty}h^{3}\right),$

where $\hat{h}_{-1}=-h_{1}$, $\hat{h}_{0}=0$, $\hat{h}_{1}=h_{2}$, $h=\max\{h_{1},h_{2}\}$, and $\|{\mbox{\boldmath$\alpha$}}\|_{\infty}=\max\left\{|\alpha_{i}|\right\}$ and so on. Note that the first and second order derivatives of $f$ can be expressed as the derivatives of the solution below,

$\displaystyle f^{(1)}(x_{i})=\kappa u^{(3)}(x_{i})+Ku^{\prime}(x_{i}),\quad\quad f^{(2)}(x_{i})=\kappa u^{(4)}(x_{i})+Ku^{\prime\prime}(x_{i}),$

(14)

from the differential equation. Thus, by collecting corresponding terms involving $u^{(l)}$, $0\leq l\leq 4$, we get a system of equations $Ax=b$ for the coefficients $\alpha_{k}$ and $\beta_{k}$, $k=-1,0,1$. Note that the constraint $\sum\beta_{k}=1$ is necessary to avoid the trivial solution. Thus, in the 1D case, we have the same number of unknowns and constraints and we have uniquely determined the coefficients. We list the coefficient matrix and the right hand side below for the special case $\kappa=1$, $K=0$.

$\displaystyle A=\left(\begin{array}[]{cccccc}1&1&1&0&0&0\\ \vskip 5.0pt\cr-h_{1}&0&h_{2}&0&0&0\\ \vskip 5.0pt\cr\frac{h_{1}^{2}}{2}&0&\frac{h_{2}^{2}}{2}&-1&-1&-1\\ \vskip 5.0pt\cr-\frac{h_{1}^{3}}{6}&0&\frac{h_{2}^{3}}{6}&h_{1}&0&-h_{2}\\ \vskip 5.0pt\cr\frac{h_{1}^{4}}{24}&0&\frac{h_{2}^{4}}{24}&-\frac{h_{1}^{2}}{2}&0&-\frac{h_{2}^{2}}{2}\\ \vskip 5.0pt\cr 0&0&0&1&1&1\end{array}\right),\quad b=\left(\begin{array}[]{c}0\\ 0\\ 0\\ 0\\ 0\\ 1\end{array}\right).$

(27)

The solution (the set of coefficients) has an analytic form that is obtained from the Maple symbolic package,

$$\begin{array}[]{rcl}&&\displaystyle\alpha_{-1}=\frac{2}{h_{1}(h_{1}+h_{2})},\quad\quad\alpha_{1}=\frac{2}{h_{2}(h_{1}+h_{2})},\quad\quad\alpha_{0}=-\left(\alpha_{l}+\alpha_{r}\right)=-\frac{2}{h_{1}h_{2}},\\ \vskip 5.0pt\cr\vskip 5.0pt\cr&&\displaystyle\beta_{-1}=\frac{h_{1}^{2}-h_{2}^{2}+h_{1}h_{2}}{6h_{1}(h_{1}+h_{2})},\quad\quad\beta_{1}=\frac{-h_{1}^{2}+h_{2}^{2}+h_{1}h_{2}}{6h_{2}(h_{1}+h_{2})},\quad\quad\beta_{0}=\frac{h_{2}^{2}+h_{1}^{2}+3h_{1}h_{2}}{6h_{1}h_{2}}.\end{array}$$

(28)

We can see the ‘symmetry’ between the left and right grid points as expected. Note that when $h_{1}=h_{2}$ we recover the standard fourth-order compact formula. The finite difference coefficients satisfy the sign property needed for an M-matrix conditions. When $K\neq 0$, we can treat $Ku$ as a source term like $f$, that is, we add the term

$\displaystyle K(\beta_{-1}U_{i-1}+\beta_{0}U_{i}+\beta_{1}U_{i+1})$

(29)

to the left hand side of the finite difference equation.

This example is from [16]:

	$\displaystyle(\kappa u_{x})_{x}=12x^{2},\quad 0<x<1,\quad\kappa=\left\{\begin{array}[]{ll}\kappa^{-}&\mbox{if $x<\alpha$}\\ \vskip 5.0pt\cr\kappa^{+}&\mbox{if $x>\alpha$},\end{array}\right.$
	$\displaystyle\displaystyle u(0)=0,\quad u(1)=\frac{1}{\kappa^{+}}+\left(\frac{1}{\kappa^{-}}-\frac{1}{\kappa^{+}}\right)\alpha^{4}.$

In this example, $f(x)=12x^{2}$ is continuous and the natural jump conditions $[u]_{\alpha}=0$ and $[\kappa u_{x}]_{\alpha}=0$ are satisfied across the interface $\alpha$. The exact solution is

$\displaystyle u(x)=\left\{\begin{array}[]{ll}\displaystyle\frac{x^{4}}{\kappa^{-}},&\mbox{if $x<\alpha$},\\ \vskip 5.0pt\cr\displaystyle\frac{x^{4}}{\kappa^{+}},&\mbox{if $x>\alpha$}.\end{array}\right.$

We use the two-grid method to solve the problem and show the numerical results in Figure 3. In the left plot, we can clearly visualize the two grids in which the parameters are $N=10,\alpha=17/30,\lambda=2$, that is, we use a fine mesh in the interval $|x-\alpha|\leq 2h$; with the refinement ratio $r=h/h_{f}=8$, and $\kappa^{-}=4,\kappa^{+}=50$. The total number of grids points is $16$, and the maximum error is $4.5343\times 10^{-7}$. In the right table, we carry out some grid refinement analysis with various coarse meshes and mesh ratio $r=2,4,8,16$. The average convergence order is $2.0875$ with respect to $h/r$ as expected. We observe better than fourth order convergence from $r=2$ to $r=8$. Note that the fine mesh is within a tube and thus, it is more efficient in terms of the accuracy and complexity compared with that of a fourth high method using a uniform mesh.

Now we consider two dimensional Poisson equations with a line interface that is parallel to the $y$-axis

$\displaystyle u_{xx}+u_{yy}=f(x,y)+C\delta(x-\alpha)+\bar{C}\delta^{\prime}(x-\alpha),\quad\quad a<x<b,\quad c<y<d,$

(32)

where $C$ and $\bar{C}$ correspond to the jumps of $u$ and $\frac{\partial u}{\partial x}$ across $\alpha$, $[u]_{\alpha}=\bar{C}$ and $[u_{x}]_{\alpha}={C}$. That is, all quantities in the $y$-direction are continuous. More general problems with piecewise constant coefficients and curved interfaces will be discussed in the next section.

In the literature, a two-grid in one space dimension is often referred to as the Shishkin mesh that has been employed to capture boundary layer effects, see for example, [8, 11, 12, 30]. To solve (32), we use a two-grid in the $x$-direction as discussed in the previous sub-section; and a uniform mesh in the $y$-direction. Then, there are three types of grid points: coarse grid points (black dotted), fine grid points (blue dotted), and boarder grid points (red dotted), as demonstrated in Figure 4.

For a coarse gird point $(x_{i},y_{j})$ that is outside of the strip $|x_{i}-\alpha|>\lambda h$, we use the standard compact nine-point finite difference scheme, see for example, [29, 17]. The fourth order compact scheme for $u_{xx}+u_{yy}+KU=f(x,y)$ assuming the same step size in the $x$ and $y$ directions can be written as

$\displaystyle\left(L_{h}+KM_{h}\right)U_{ij}=M_{h}f_{ij},$

(33)

where $L_{h}$ is the discrete nine-point Laplacian operator whose coefficients at the four corners are $\frac{1}{6h^{2}}$, at the east-north-south-west grid points are $\frac{4}{6h^{2}}$, and at the center is $-\frac{20}{\,6h^{2}}$; $M_{h}$ is an averaging operator,

$\displaystyle M_{h}f_{ij}=\frac{1}{12}\left(\frac{\hbox{}}{\hbox{}}f_{i-1,j}+f_{i+1,j}+f_{i,j-1}+f_{i,j+1}+8f_{i,j}\right),$

(34)

where $f_{ij}=f(x_{i},y_{j})$ and so on.

In the strip $|x_{i}-\alpha|\leq\lambda h$ as shown in Figure 4, that is, a fine grid in the $x$ direction but the same coarse grid size in the $y$ direction, with the immersed interface method (IIM) as discussed in Subsection 2.1, the finite difference equation that is second order in $x$ and fourth order in $y$ in the tube can be written as

$\displaystyle\frac{U_{i-1,j}-2U_{ij}+U_{i+1,j}}{h_{f}^{2}}+\left(1+{\frac{h^{2}}{12}}\delta_{yy}^{2}\right)^{-1}\delta_{yy}^{2}U_{ij}=f_{ij}+C_{ij},$

(35)

where $\delta_{yy}^{2}$ is the standard three-point central finite difference operator in the $y$-direction and $C_{ij}$ is the correction term which is zero almost everywhere except at the two grid points ($x_{j}$ and $x_{j+1}$) neighboring the interface $\alpha$. Note also that $C_{ij}$ is independent of $y$. Thus, the discretization can further be written as

$\displaystyle\left(1+{\frac{h^{2}}{12}}\delta_{yy}^{2}\right)\frac{U_{i-1,j}-2U_{ij}+U_{i+1,j}}{h_{f}^{2}}+\delta_{yy}^{2}U_{ij}=\left(1+{\frac{h^{2}}{12}}\delta_{yy}^{2}\right)\left({\frac{\hbox{}}{\hbox{}}}f_{ij}+C_{ij}\right),$

(36)

which leads to

	$\displaystyle\!\!\!\!\!\frac{U_{i-1,j}-2U_{ij}+U_{i+1,j}}{h_{f}^{2}}+C_{ij}+\frac{h_{y}^{2}}{12h_{f}^{2}}\left(\frac{\hbox{}}{\hbox{}}U_{i-1,j-1}-2U_{i-1,j}+U_{i-1,j+1}-2U_{i,j-1}\right.$
	$\displaystyle\quad\left.\frac{\hbox{}}{\hbox{}}+4U_{i,j}-2U_{i,j+1}+U_{i+1,j-1}-2U_{i+1,j}+U_{i+1,j+1}\right)=\frac{1}{12}\left(f_{i,j-1}+10f_{i,j}+f_{i,j+1}\frac{\hbox{}}{\hbox{}}\right).$

The finite difference equation involves the nine-point compact stencil, and is second order accurate in $x$ and fourth order accurate in $y$, respectively.

For the boarder grid points as shown in Figure 4, following the similar derivation process for 1D boarder grid points in subsection 2.2, we can obtain a compact finite difference equation at a boarder grid point $(x_{i},y_{j})$. The finite difference coefficients for $U_{i,j}$ and the linear combination of $f_{i,j}$ are listed below,

$$U_{i,j}:\quad\frac{1}{6\sigma h_{y}^{2}}\begin{bmatrix}h_{2}(h_{1}h_{2}+h_{y}^{2}+\mu)&(h_{1}+h_{2})(3h_{1}h_{2}-h_{y}^{2}+\nu)&h_{1}(h_{1}h_{2}+h_{y}^{2}-\mu)\\ -2h_{2}(h_{1}h_{2}-5h_{y}^{2}+\mu)&-2(h_{1}+h_{2})(3h_{1}h_{2}+5h_{y}^{2}+\nu)&-2h_{1}(h_{1}h_{2}-5h_{y}^{2}-\mu)\\ h_{2}(h_{1}h_{2}+h_{y}^{2}+\mu)&(h_{1}+h_{2})(3h_{1}h_{2}-h_{y}^{2}+\nu)&h_{1}(h_{1}h_{2}+h_{y}^{2}-\mu)\end{bmatrix},$$

$$f_{i,j}:\quad\frac{1}{12\sigma}\begin{bmatrix}0&\sigma&0\\ 2h_{2}(h_{1}h_{2}-\nu)&2(h_{1}+h_{2})^{3}&2h_{1}(h_{1}h_{2}+\nu)\\ 0&\sigma&0\end{bmatrix},$$

where $h_{1}=x_{i}-x_{i-1},h_{2}=x_{i+1}-x_{i},h_{y}=y_{j}-y_{j-1}=y_{j+1}-y_{j}$ and $\sigma=h_{1}h_{2}(h_{1}+h_{2}),\ \mu=h_{1}^{2}-h_{2}^{2},\ \nu=h_{1}^{2}+h_{2}^{2}.$

Now we have defined finite difference equations at all grid points with local truncation errors being $O(h^{4}+h_{f}^{2})$. The discretize system of equations satisfy the discrete maximum principle, and thus, the global accuracy is also of $O(h^{4}+h_{f}^{2})$.

In the left plot of Figure 5, we show an error plot with $N=42$ and $h_{f}=h^{2}$ in which two grids with different mesh sizes are also visible. In the right table, we show a grid refinement analysis which indicates fourth order convergence in terms of $h$ when $h_{f}=h^{2}$ as expected.

One of the challenges for two-grid or adaptive mesh methods is how to treat hanging nodes, which has been intensively discussed in the literature. There are two main concerns about the discretization: accuracy, and the stability that is often more difficult to study. One sufficient (stronger) condition is to check the sign property of finite difference coefficients. Our new idea is to use the values of the solution, the source term $f$, and the partial differential equations to derive finite difference equations that can be accurate for highest polynomials possible while satisfying the sign property. We also want the finite difference scheme to be robust so that the coefficients just need to be computed once for PDEs with constant coefficients.

There are various approaches in designing finite difference equations at hanging nodes in the literature. One commonly used strategy is to use ghost grid points, and then apply interpolation schemes. We refer the readers to [18] and references therein. Often, the developed finite difference equations are consistent. However, it is difficult to show the stability, and the convergence of the finite difference method.

In this paper, we employ a new strategy, see also [21], and some new ideas to develop finite difference equations at hanging nodes for elliptic PDEs

$\displaystyle\kappa\Delta u+Ku=f(\mathbf{x}),$

(45)

assuming that $\kappa$ and $K$ are constants. We use the left diagram in Figure 7 to illustrate the method in which $r=4$. We develop a finite different equation for the hanging node $(x_{i_{c}},y_{j_{c}})$ (marked as a red five-star). In the figure, solid black circles are fine grid points, small rectangles are coarse ones. A grid point can be both fine and coarse. In fact, for the PDE above with constant coefficients, the coefficients of the finite difference equation and the weights for the source term are invariant in terms of the relative locations and the PDE. So we just need to use a reference geometry in the derivation as illustrated in Figure 7, in which there are three such hanging nodes along the line $x=0$, and another three along the line $y=0$.

After some research, testing, and reasoning based on symmetry arguments, we come to a conclusion that a seven-point stencil as demonstrated in Figure 7 is an ideal choice as we can see later in this section. In order to derive a high order compact (HOC) finite difference equation (7-point stencil) at the hanging node (marked as the red star), we use its neighboring six grid points (marked in blue square). The finite difference stencil has one coarse grid point, and six coarse and fine ones. Denote the coordinates at hanging node $(x_{i_{c}},y_{j_{c}})$ as $(h_{2},0)$. Then, the grid points involved are $(0,\pm h)$, $(h,\pm h)$, $(0,0)$, $(h,0)$ and $(h_{2},0)$. Note that $h_{2}=jh_{f}$ for $j=1,2,\cdots,r-1$. It is obvious that the finite difference equations for hanging points with $j=k$ and $j=r-k$ are the same using the symmetry arguments, which has been verified by symbolic computations. We design the finite difference equation as:

$$\begin{array}[]{rcl}&&\!\!\!\!\!\!\!\!\displaystyle\sum_{i_{k}=0}^{1}\sum_{j_{k}=-1}^{1}\!\!\alpha_{i_{k},j_{k}}U_{i+i_{k},j+j_{k}}+\alpha_{i_{c},j_{c}}U_{i_{c},j_{c}}=\sum_{i_{k}=0}^{1}\sum_{j_{k}=-1}^{1}\!\beta_{i_{k},j_{k}}\,f(x_{i+i_{k}},y_{j+j_{k}})+\beta_{i_{c},j_{c}}\,f(x_{i_{c}},y_{j_{c}}),\\ \vskip 5.0pt\cr&&\displaystyle\sum_{i_{k}=-1}^{1}\sum_{j_{k}=-1}^{1}\!\!\!\beta_{i_{k},j_{k}}+\beta_{i_{c},j_{c}}=1.\end{array}$$

(46)

Note that the indexes  are not real ones for actual algorithms but just for the derivation of the finite difference equation. The coefficients are just needed to be computed once for all for constant $\kappa$ and $K$.

Define the ‘local truncation error’ as

$$\begin{array}[]{rcl}T_{{i_{c}},{j_{c}}}&=&\displaystyle\sum_{i_{k}=0}^{1}\sum_{j_{k}=-1}^{1}\!\!\alpha_{i_{k},j_{k}}u(x_{i+i_{k}},y_{j+j_{k}})+\alpha_{i_{c},j_{c}}u(x_{i_{c}},y_{j_{c}})\\ &&\hbox{}\quad\quad\quad\displaystyle-\sum_{i_{k}=0}^{1}\sum_{j_{k}=-1}^{1}\!\beta_{i_{k},j_{k}}\,f(x_{i+i_{k}},y_{j+j_{k}})-\beta_{i_{c},j_{c}}\,f(x_{i_{c}},y_{j_{c}}),\end{array}$$

(47)

where $u(x,y)$ is the true solution to the boundary value problem. To determine the finite difference coefficients $\alpha$’s and the weight coefficients $\beta$’s, we expand all $u$ and $f$ terms at $(x_{i_{c}},y_{j_{c}})$ up to all fourth order derivatives of $u(x,y)$ and second order derivatives of $f(x,y)$; then we represent $f(x,y)$ terms in terms of $u(x_{c},y_{c})$ and its partial derivatives at the $(x_{c},y_{c})$ using the following relations,

$$\frac{\partial^{k_{1}+k_{2}}f}{\partial x^{k_{1}}\partial y^{k_{2}}}=\kappa\left(\frac{\partial^{k_{1}+k_{2}+2}u}{\partial x^{k_{1}+2}\partial y^{k_{2}}}+\frac{\partial^{k_{1}+k_{2}+2}u}{\partial x^{k_{1}}\partial y^{k_{2}+2}}\right)+K\frac{\partial^{k_{1}+k_{2}}u}{\partial x^{k_{1}}\partial y^{k_{2}}},\quad 0\leq k_{1}+k_{2}\leq 2.$$

(48)

With the above relations, we can rewrite the local truncation error as

$\displaystyle T_{{i_{c}},{j_{c}}}=\sum_{0\leq k_{1}+k_{2}\leq 4}\!\!\!\!\!L_{k_{1},k_{2}}\,\left.\frac{\partial^{k_{1}+k_{2}}u}{\partial x^{k_{1}}\partial y^{k_{2}}}\right|_{(x_{c},y_{c})}+O\left(\|{\mbox{\boldmath$\alpha$}}\|_{\infty}h^{5}+\|{\mbox{\boldmath$\beta$}}\|_{\infty}h^{3}\right),$

(49)

where $\|{\mbox{\boldmath$\alpha$}}\|_{\infty}$ is the maximum norm of all $\alpha$’s coefficients and so on. We want $T_{{i_{c}},{j_{c}}}$ to be as small as possible in magnitude in terms of $h$. So we set the coefficients of $\frac{\partial^{k_{1}+k_{2}}u}{\partial x^{k_{1}}\partial y^{k_{2}}}$ to be zero for all the $h^{k}$ coefficients, which leads to a linear system of equations for the coefficients $\alpha$’s and the weight coefficients $\beta$’s. The equality constraint in (46) prevents the trivial solution (all zero coefficients). If we make all the coefficients of $\frac{\partial^{k_{1}+k_{2}}u}{\partial x^{k_{1}}\partial y^{k_{2}}}$ to be zero for all $0\leq k_{1}+k_{2}\leq 4$, we would obtain 16 equations with 14 unknowns (7 $\alpha$’s and 7 $\beta$’s). It turns out that the system of equations is not consistent. Thus, we give up two equations corresponding to $x^{4}$ and $y^{4}$ terms to have $14$ equations and $14$ unknowns. The system of equations then is consistent and have infinity number of solutions. We can impose some additional conditions such as the symmetry, sign property to get desired solutions.

With the help of the Matlab symbolic package, we have analytically found one ‘good’ set of finite difference coefficients and weights of the source term when $\kappa=1$ and $K=0$ as listed in Table 1 and Table 2 for $r=2,4,8$. The last column in Table 1 lists the coefficients $\alpha_{i_{c},j_{c}}$ at the diagonals. The coefficients for $r=16$ are listed in Table 3. From these talbles, we have verified following properties of the finite difference scheme, which we referred as a ‘good’ set of the coefficients.

• •

  The finite difference coefficients satisfy the sign property, that is, all diagonals are negative while off-diagonals are positive. So the coefficient matrix for the finite difference equations is an M-matrix.

• •

  The finite difference coefficients and weights at hanging point $(jh_{f},0)$ are the same for any $j$ and $r$ with $j/r=\textrm{constant}$.

• •

  The finite difference coefficients and weights at hanging point $(jh_{f},0)$ for $j=k$ are the same as those for $j=r-k$ except they are in the reverse order.

• •

  Only two weight coefficients are nonzero: $\beta_{0,1}$ and $\beta_{1,1}$, which are positive and their sum is one.

• •

  The discretization at a hanging node is super-third order accurate. In our numerical tests, the errors from the discretization at hanging nodes are less dominant compared with those at irregular grid points near the interfaces in the finer mesh.

Now we discuss the convergence of the two-grid method. The proof is based on the discrete maximum principle characterized by the M-matrix property and the local truncation errors.