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

Zhilin Li Department of Mathematics, North Carolina State University, Raleigh,
NC 27695-8205, USA, Email<EMAIL_ADDRESS>Kejia Pan School of
Mathematics and Statistics, HNP-LAMA, Central South University, Changsha,
Hunan, 410083, China, Email<EMAIL_ADDRESS>Juan Ruiz-Álvarez
Departamento de Matemática Aplicada y Estadística. Universidad Politécnica de
Cartagena, Spain, Email<EMAIL_ADDRESS>
###### Abstract

Second order accurate Cartesian grid methods have been well developed for
interface problems in the literature. However, it is challenging to develop
third or higher order accurate methods for problems with curved interfaces and
internal boundaries. In this paper, alternative approaches based on two
different grids are developed for some interface and internal layer problems,
which are different from adaptive mesh refinement (AMR) techniques. For one
dimensional, or two-dimensional problems with straight interfaces or boundary
layers that are parallel to one of the axes, the discussion is relatively
easy. One of challenges is how to construct a fourth order compact finite
difference scheme at boarder grid points that connect two meshes. A two-grid
method that employs a second order discretization near the interface in the
fine mesh and a fourth order discretization away from the interface in the
coarse and boarder grid points is proposed. For two dimensional problems with
a curved interface or an internal layer, a level set representation is
utilized for which we can build a fine mesh within a tube $|\varphi({\bf
x})|\leq\delta h$ of the interface. A new super-third seven-point
discretization that can guarantee the discrete maximum principle has been
developed at hanging nodes. The coefficient matrices of the finite difference
equations developed in this paper are M-matrices, which leads to the
convergence of the finite difference schemes. Non-trivial numerical examples
presented in this paper have confirmed the desired accuracy and convergence.

keywords: Two-grid method, interface problems, internal layer problems, finite
difference method, discrete maximum principle, hanging nodes.

AMS Subject Classification 2000 65N06, 65N50.

##  1 Introduction

High order compact schemes are useful or even required for many applications
such as high wave numbers, oscillatory solutions, and problems on large or
infinite domains. In general, high order compact schemes are suitable for
rectangular meshes and boundaries. For interface problems, the solutions often
depend on the curvature of the interface, which presents challenges to develp
high order discretizations near the interface. Various efforts have been made
recently. For singular source only interface problems, fourth order schemes
for Poisson equations have been developed, see for example [33, 25]; and a
sixth order method if the interface has an analytic expressions [7]. A third
order scheme for elliptic partial differential equations (PDEs) with piecewise
constant but discontinuous coefficients in the presence of fixed interfaces in
[25]. The implementation for high order schemes for interface problems with
curved interface is rather complicated. Note that the curvature of an
interface involves second order derivatives, and is non-linear. A related
topic is how to resolve the solution accurately near an internal layer, which
is a long outstanding challenging problem in computational fluid dynamics and
other applications especially if the boundary layer is not parallel to one of
the axes. It is well-known that often the accuracy of the pressure is one
order lower than that of the velocity due to the boundary layer effect in
solving incompressible Naiver-Stokes equations [4, 3, 5]. In some sense, to
resolve the solution near an internal (boundary) layer is like solving an
interface problem in a structured mesh.

Our idea in this paper is to use two-grid methods to solve these challenging
problems. The idea is different from adaptive mesh refinement (AMR) techniques
[23, 22, 2, 10, 9, 18, 19]. For one dimensional (1D) interface problems, we
use a fourth order method on a mesh with a step size $h$ away from the
interface; and a second order method with a finer mesh with step size $h/r$
near the interface. Thus, we can get an approximation with
$O(h^{4})+O((h/r)^{2})$ global accuracy. In the left diagram in Figure 1, we
show a diagram of a two-grid for an interface problem (top); and a boundary
layer problem (bottom). In the right plots in Figure 1, we show the solution
and error plots for a famous boundary layer problem. In Figure 5, we show the
solution and a two-grid for a two-dimensional example with a line interface;
and in Figure 6 a two-grid for a circular interface. Interested readers can
get an idea of the proposed two-grid methods from those plots.

Note also that there are different interpretations of ‘two-grid methods’ in
the literature including different meshes for different variables [34, 6, 1];
different meshes for interior and boundary domains [32, 1, 35, 31]. For
boundary layers that are parallel to one of the axes, a Shishkin mesh has been
employed to capture the boundary layer effect, see for example, [8, 11, 12,
30].

In this paper, we develop high order two-grid methods that are as compact as
possible with some new ideas. One of the keys in our methods is that we try to
enforce the sign properties for elliptic differential equations so that the
discrete maximum principle can be preserved, which leads to the convergence of
the finite difference methods. For one-dimensional problems, we have developed
a fourth-order compact scheme for Poison/Helmholtz equations with non-uniform
grids. For two-dimensional (2D) interface problems with a line interface that
is parallel to one of the axes, we have developed a new finite difference
method that is fourth order away from the interface, second order in one
coordinate direction with a fine mesh, and fourth order in the other with a
coarse mesh. For 2D general interface problems, the algorithm and discussion
are more complicated. One of difficulties is how to develop high order
discretizations at hanging nodes, which can preserve the discrete maximum
principle. Different from a traditional approach using ghost points and
interpolations, which cannot guarantee the stability of the algorithm, we use
a new undetermined coefficients method so that the sign property can be
satisfied with the highest possible accuracy. The convergence of the developed
two-grid methods is guaranteed with the sign property and the local truncation
error analysis for elliptic interface problems and internal layer problems.

The rest of paper is organized as follows. In the next section, we discuss the
two-grid method for one-dimensional interface problems, and two-dimensional
interface problems with line interfaces that are parallel to one of the axes.
Numerical examples are also presented. In Section 3, we discuss the two-grid
method for general elliptic interface problems, some implementation details,
and a new finite difference discretization at hanging nodes, which can
preserve the discrete maximum principle. Several benchmark examples from the
literature are tested and analyzed. We also present a numerical experiment of
a non-trivial internal layer problem. We conclude in the last section.

##  2 Two-grid methods for 1D interface problems and 2D problems with line
interfaces

We first show a simple and well-known one-dimensional boundary layer example,

$\displaystyle\epsilon u^{\prime\prime}-u^{\prime}+u=f(x),\quad
0<x<1,\quad\quad u(0)=1,\quad u(1)=3.$ (1)

The solution exhibits a boundary layer near $x=1$ if $\epsilon$ is small.

Figure 1: Diagrams of the two-grids method in one-dimension. Left-top: for an
interface problem; Left-bottom: for a boundary layer problem. Right plots: A
computed boundary layer solution and the error ($\sim 10^{-6}$) with a coarse
mesh size $h=1/N=1/80$, and a fine mesh $h_{f}=h^{2}$ in $(1-2h,1)$.

For this example, we show the computed solution and error plot in Figure 1
using a standard centered finite difference discretization. Let $h$ be a
coarse grid size and $N$ be the number of grid points in the interval
$(0,\,1)$. In the interval $[1-2h,1]$, where $h=1/N$, we use a fine mesh
$h_{f}=h^{2}$ for the discretization so that the central discretization has
little effect on the solution even though that $|u^{\prime}|$ is relatively
large.

Next, we use one-dimensional interface problems,

$\displaystyle(\kappa\,u^{\prime})^{\prime}+Ku=f(x)+C\delta(x-\alpha)+\bar{C}\delta^{\prime}(x-\alpha),\quad\quad
a<x<b,\quad a<\alpha<b,$ (2)

with two linear boundary conditions (Dirichlet, Neumann, Robin) at $x=a$ and
$x=b$ to demonstrate the two-grid method. The problem is equivalent to

$\displaystyle(\kappa
u^{\prime})^{\prime}+Ku=f(x),\quad\quad[u]_{\alpha}=\frac{2\bar{C}}{\kappa^{-}+\kappa^{+}},\quad[\kappa
u^{\prime}]_{\alpha}=C,\quad a<x<b,$ (3)

see [20]. We assume that the coefficient $\beta$ is a piecewise constant,

$\displaystyle\kappa=\left\\{\begin{array}[]{ll}\kappa^{-}&\mbox{if
$a<x<\alpha$ }\\\ \vskip 5.0pt\cr\kappa^{+}&\mbox{ if $\alpha<x<b$.
}\end{array}\right.$ (6)

Given a coarse grid parameter $N$, a grid ratio $r$, and a width parameter
$\lambda>0$, we can easily generate a two-grid. First, we define $h=(b-a)/N$,
and let the fine mesh size be $h_{f}=h/r$. The fine mesh is defined as
$x_{j}=a+jh_{f}$ if $|x_{j}-\alpha|\leq\lambda h$. The coarse mesh is defined
as those $x_{i}=a+hi$ and $|x_{i}-\alpha|>\lambda h$, see the left plot in
Figure 3 for an illustration.

The standard fourth order compact scheme at a grid point $x_{i}$ for a uniform
mesh, i.e. $x_{i+1}-x_{i}=x_{i}-x_{i-1}=h$, is

$\displaystyle\kappa\frac{U_{i-1}-2U_{i}+U_{i+1}}{h^{2}}=\frac{f_{i-1}+10f_{i}+f_{i+1}}{12}.$
(7)

###  2.1 The IIM scheme near the interface

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.

Figure 2: A diagram of some grid points and the interface $\alpha$ that is
between $x_{j}$ and $x_{j+1}$.

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].

###  2.2 High-order compact FD discretization at boarder grid points

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.

###  2.3 A one-dimensional numerical example

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.$

A two-grid along the $x$-axis, the exact solution (solid line) and computed
solution (little o’s) with $N=10$, $\lambda=2$, and $r=8$. The maximum error
is $4.5343\times 10^{-7}$.

$N$ $\|E\|_{\infty}$ (uniform) $\|E\|_{\infty}$ (two-grid) 10 $4.4333\times
10^{-4}$ $7.2861\times 10^{-5}$ $1.1428\times 10^{-5}$ $4.5343\times 10^{-6}$
$6.1344\times 10^{-7}$ 20 $1.3470\times 10^{-4}$ $1.0335\times 10^{-5}$
$4.5343\times 10^{-6}$ $6.1344\times 10^{-7}$ $2.8222\times 10^{-7}$ 40
$4.1331\times 10^{-5}$ $4.4768\times 10^{-6}$ $6.1344\times 10^{-7}$
$2.8222\times 10^{-7}$ $3.9592\times 10^{-8}$

Figure 3: Computed results of the 1D example.

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.

###  2.4 2D problems with line interfaces that are parallel to one of the axes

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.

Figure 4: A two-grid with a line interface (red solid). There are three types
of grid points: coarse grid points (black dotted), fine grid points in the
$x$-direction (blue dotted), and boarder grid points (red dotted) that connec
the coarse and fine meshes in the $x$-direction. The refinement ratio is
$r=5$, and the refinement tube width is $\lambda=1.5$ in the figure.

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})$.

######  Example 2.1.

This example is adapted from the 1D example in Subsection 2.3. A Poisson
equation with a line interface $x=\alpha=\frac{33}{70}$ on a domain $0<x<1$,
$0<y<1$ with the true solution,

$\displaystyle u(x,y)=\left\\{\begin{array}[]{ll}x(\alpha-1)+\sin(\pi
y)&\textrm{if $x\leq\alpha$}\\\ \vskip 5.0pt\cr\alpha(x-1)+\sin(\pi
y)&\textrm{if $x>\alpha$.}\end{array}\right.$ (39)

The source term is $f(x,y)=-\pi^{2}\sin(\pi y)$. Across the interface, the
jump conditions are $[u]_{x=\alpha}=0$ and $\left[\frac{\partial u}{\partial
x}\right]_{x=\alpha}=1$.

A two grid along the $x$-axis and the error plot with $N=42$, $\lambda=2$, and
$h_{f}=h^{2}$. The maximum error is $7.8476\times 10^{-8}$.

$N$ $\|E\|_{\infty}$(two-grid) order $6$ $1.6355\times 10^{-4}$ $12$
$1.1807\times 10^{-5}$ 3.9585 $24$ $7.3631\times 10^{-7}$ 4.0032 $42$
$7.8476\times 10^{-8}$ 4.1156 A grid refinement analysis for the 2D example
with a line interface. Fourth order convergence is confirmed when
$h_{f}=h^{2}$.

Figure 5: An error plot of the computed solution of the 2D example with a line
interface in which two grids in the $x$ direction are visible.

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.

##  3 A two-grid method for general 2D interface problems

For general two-dimensional elliptic interface problems, it is much more
challenging to develop two-grid methods with general curved interfaces. We
consider a general two-dimensional elliptic interface problem of the following
form,

$\displaystyle\displaystyle\nabla\cdot\left(\kappa\nabla
u\right)+Ku=f(\mathbf{x}),\quad\mathbf{x}\in R,$ (40)

together with jump conditions across interface $\Gamma$

$\displaystyle\displaystyle\left[u\right]_{\Gamma}=w(s),\quad\quad\left[\kappa\frac{\partial
u}{\partial n}\right]_{\Gamma}=v(s),$ (41)

and a boundary condition along $\partial R$, where $R$ is a rectangular
domain, $\Gamma\in C^{2}$ is a closed smooth interface that can be
parameterized by a one-dimensional variable $s$, say the arc-length. Within
the domain $R$; $w(s)\in C^{2}$ and $v(s)\in C^{1}$ are two functions defined
along the interface $\Gamma$. Note that, when $w(s)=0$, the problem can be
written as

$\displaystyle\displaystyle\nabla\cdot\left(\kappa\nabla
u\right)+Ku=f(\mathbf{x})+\int_{\Gamma}v(\mathbf{X}(s))\delta\left(\mathbf{x}-\mathbf{X}(s)\right)ds,\quad\mathbf{x}\in
R.$ (42)

Second order methods on Cartesian grids have been developed in the literature.
The challenge here is to develop higher (third or above) order compact
methods. It is more useful but challenging for the case when $\kappa$ is
discontinuous. In this paper, we propose an alternative approach to a uniform
mesh, that is, a two-grid method to obtain high order accurate methods, which
is not only simpler but also avoids higher order derivatives required for jump
conditions along the interface as well as the source terms and the solution.

We assume that the interface problem is defined in a rectangular domain
$\Omega=[a,\,b]\times[c,\,d]$. We start with a coarse Cartesian grid,
$x_{i}=a+ih$, $y_{j}=c+jh$, $i=0,1,\cdots,M$, $j=0,1,\cdots,N$. The interface
$\Gamma$ is implicitly represented by the zero level set of a Lipschitz
continuous function $\varphi(x,y)$:

$\displaystyle\Gamma=\left\\{\,(x,y),{\frac{\hbox{}}{\hbox{}}}\varphi(x,y)=0\right\\}.$
(43)

In the discrete case, $\varphi(x,y)$ is defined at grid points as
$\varphi(x_{i},y_{j})$. Often $\varphi(x,y)$ is the signed distance function
from $\Gamma$, or its approximations.

Figure 6: A two-grid with a circular interface (red solid). Parent grid points
(starred) are selected within the tube
$|\varphi_{ij}|=|\varphi(x_{i},y_{j})|\leq h$ (red dashed). The refinement
ratio is $r=2$, and the refinement tube width is $\lambda=1$.

To generate a finer mesh around the interface $\Gamma$, we first select parent
grid points within a tube of the interface according to

$\displaystyle|\,\varphi(x,y)\,|\leq\lambda h,$ (44)

where $\lambda$ is a control coefficient to adjust the width of the refinement
tube. When a grid point $\textbf{x}_{ij}=(x_{i},y_{j})$ within the tube is
selected, we build a refined mesh with a finer mesh $h/r$ ($r$ is refinement
ratio, $r=2$ or $4$, $\cdots$) within the square: $|x-x_{i}|\leq h$ and
$|y-y_{j}|\leq h$. Generating a refined square for every such grid points
yields a refined region around the interface. Figure 6 shows an example of the
refinement mesh around a circular interface. The diagram at the right in
Figure 7 is an illustration of some grid points in the coarse and finer grid
with $r=4$.

Once we have the two grids, we use the standard fourth-order compact finite
difference scheme at interior coarse grid points; and the standard second
order method at interior regular fine grid points; and second order maximum
principle immersed interface method (IIM) [15, 16] at the fine irregular grid
points within the tube. In all the three cases, the finite difference
coefficients satisfy the sign property that guarantee the discrete maximum
principle.

###  3.1 A new super-third order FD discretization for hanging nodes

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.

Figure 7: A schematic diagram of 7-point finite difference stencil for a
typical hanging node (marked as a red five-star) when $r=4$ ($h_{f}=h/4$).

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.

Table 1: Finite difference coefficients $h^{2}\alpha_{i_{k},j_{k}}$ at hanging point $(jh_{f},0)(j=1,2,\cdots,r-1)$ for $r=2$, $4$, $8$. The finite difference coefficients for $j=k$ are the same as those for $j=r-k$ except that they are in the reverse order, and the coefficients are the same for any $j$ and $r$’s when $j/r=\textrm{constant}$. $r$ | $j$ | $h^{2}\alpha_{i_{k},j_{k}}$ | $h^{2}\alpha_{i_{c},j_{c}}$
---|---|---|---
2 | 1 |  1/2 |  1/2 |  3 |  3 |  1/2 |  1/2 |  -8
4 | 1 | 7/12 | 5/12 | 41/6 | 11/6 | 7/12 | 5/12 |  -32/3
2 |  1/2 |  1/2 |  3 |  3 |  1/2 |  1/2 |  -8
3 | 5/12 | 7/12 | 11/6 | 41/6 | 5/12 | 7/12 |  -32/3
8 | 1 | 5/8 | 3/8 | 59/4 | 43/28 | 5/8 | 3/8 | -128/7
2 | 7/12 | 5/12 | 41/6 | 11/6 | 7/12 | 5/12 |  -32/3
3 | 13/24 | 11/24 | 17/4 | 137/60 | 13/24 | 11/24 | -128/15
4 |  1/2 |  1/2 |  3 |  3 |  1/2 |  1/2 |  -8
5 | 11/24 | 13/24 | 137/60 | 17/4 | 11/24 | 13/24 | -128/15
6 | 5/12 | 7/12 | 11/6 | 41/6 | 5/12 | 7/12 |  -32/3
7 | 3/8 | 5/8 | 43/28 | 59/4 | 3/8 | 5/8 | -128/7
Table 2: Weight coefficients $\beta_{i_{k},j_{k}}$ at hanging point $(jh_{f},0)$ for different $j/r(r=2,4,8)$.  $j/r$ | $\beta_{i_{k},j_{k}}$ | $\beta_{i_{c},j_{c}}$
---|---|---
1/2 | 0 | 0 | 1/2 | 1/2 | 0 | 0 | 0
1/4 | 0 | 0 | 7/12 | 5/12 | 0 | 0 | 0
1/8 | 0 | 0 | 5/8 | 3/8 | 0 | 0 | 0
3/8 | 0 | 0 | 13/24 | 11/24 | 0 | 0 | 0

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.

Table 3: Finite difference coefficients $h^{2}\alpha_{i_{k},j_{k}}$ and weight coefficients $\beta_{i_{k},j_{k}}$ at the hanging grid point $(jh_{f},0)(j=1,2,\cdots,r/2)$ when $r=16$. The coefficients for $j=9,10,\cdots,15$, which are not listed here, are the same as those for $j=7,6,\cdots,1$, respectively. Similar to the cases for $r=2,4,8$ (see Table 2), only two weight coefficients $\beta_{0,0}$ and $\beta_{1,0}$ are nonzero.  $j$ | $h^{2}\alpha_{i_{k},j_{k}}$ | $h^{2}\alpha_{i_{c},j_{c}}$ | $\beta_{0,0}$ | $\beta_{1,0}$
---|---|---|---|---
1 | 31/48 | 17/48 | 737/24 | 57/40 | 31/48 | 17/48 | -512/15 | 31/48 | 17/48
2 | 5/8 | 3/8 | 59/4 | 43/28 | 5/8 | 3/8 | -128/7 | 5/8 | 3/8
3 | 29/48 | 19/48 | 227/24 | 521/312 | 29/48 | 19/48 | -512/39 | 29/48 | 19/48
4 | 7/12 | 5/12 | 41/6 | 11/6 | 7/12 | 5/12 | -32/3 | 7/12 | 5/12
5 | 9/16 | 7/16 | 211/40 | 179/88 | 9/16 | 7/16 | -512/55 | 9/16 | 7/16
6 | 13/24 | 11/24 | 17/4 | 137/60 | 13/24 | 11/24 | -128/15 | 13/24 | 11/24
7 | 25/48 | 23/48 | 593/168 | 187/72 | 25/48 | 23/48 | -512/63 | 25/48 | 23/48
8 | 1/2 | 1/2 | 3 | 3 | 1/2 | 1/2 | -8 | 1/2 | 1/2

###  3.2 Convergence analysis

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.

######  Theorem 1.

Let $U_{ij}$ be the finite difference solution obtained from the two-grid
method for $\kappa\Delta u=f$ with a linear boundary condition (Dirichlet,
Neumann, or Robin), and at least one point Dirichlet boundary condition.
Assume that the finite difference coefficients satisfy the sign property at
all the grid points with $\|{\mbox{\boldmath$\alpha$}}\|_{\infty}\leq C/h^{2}$
and $\|{\mbox{\boldmath$\beta$}}\|\leq C$, and the solution to the elliptic
interface problem is piecewise smooth with the needed regularity, then,

$\displaystyle\left\|u(x_{i},y_{j})-U_{ij}\right\|_{\infty}\leq
C\left(O(h^{3})+O\left({h/r}\right)^{2}\right).$ (50)

Here, $C$ is a generic constant, $h$ is the step size of the coarse mesh and
$h/r$ is that of the fine mesh.

###### Proof.

Note that, in the interior domain, the regularity of the solution depends on
the source term $f$. In the neighborhood of coarse grid points, the needed
regularity is $u(x,y)\in C^{6}$ and we have $|T_{ij}|\leq Ch^{4}$. In the
neighborhood of hanging nodes, the needed regularity is $u(x,y)\in C^{5}$ and
we know that the local truncation errors are bounded by $|T_{ij}|\leq Ch^{3}$.
In the neighborhood of fine grid points including irregular grid points near
the interface, the needed regularity is $u(x,y)$ is piecewise $C^{3}$ and we
have $|T_{ij}|\leq Ch/r$. Note also that the coefficient matrix of the finite
difference equations is an M-matrix. Thus, from Theorem 6.1 and Theorem 6.2 of
Morton & Mayer’s book, [24] we obtain the conclusion.

######  Remark 1.

  * •

We have obtained explicit expressions for the coefficients and the weights for
$r=2,4,\cdots,16$ that satisfy the conditions of the theorem. Thus, the
convergence is guaranteed and confirmed for those $r$’s.

  * •

In our numerous numerical tests, the largest errors seem to occurr near the
interface rather than at the hanging nodes. We conjecture that the error
should be bounded by $\left\|u(x_{i},y_{j})-U_{ij}\right\|_{\infty}\leq
C\left(O(h^{4})+O\left({h/r}\right)^{2}\right)$ instead of $O(h^{3})$.

###  3.3 Numerical experiments for general 2D interface problems

We have tested the two-grid method for various benchmark problems in the
literature. We list some of the results in this subsection.

######  Example 3.1.

Peskin’s IB model in two dimensions.

This example was first presented in [13], and later was employed by the AMR-IB
method in [27, 28]. The partial differential equation is

$\displaystyle
u_{xx}+u_{yy}=\int_{\Gamma}{2\,\delta(x-X(s))\delta(y-Y(s))\,ds},$ (51)

where $\Gamma$ is a circle interface with the radius $r=0.5$ within the square
domain $R=(-1,1)\times(-1,1)$. The problem can be written equivalently as,

$\displaystyle
u_{xx}+u_{yy}=0,\quad\mbox{on}\quad\Omega\backslash\Gamma,\quad\quad[u]_{\Gamma}=0,\quad\left[\frac{\partial
u}{\partial n}\right]_{\Gamma}=2,$ (52)

in the domain $R$. The true solution to the interface problem is

$\displaystyle u(x,y)=\left\\{\begin{array}[]{ll}1&\textrm{if $r\leq 1/2$}\\\
\vskip 5.0pt\cr 1+\log(2r)&\textrm{if $r>1/2,$}\end{array}\right.$ (55)

where $r=\sqrt{x^{2}+y^{2}}$. The Dirichlet boundary condition is applied
according to the true solution along the boundary of the square.

Table 4: Comparison of numerical results using IB, IIM, and the two-grid method for Example 3.1, with $r=2,4,8$. $N$ | IB | IIM |  $\|E\|_{\infty}$ (coarse ) |  $\|E\|_{\infty}$ (fine)
---|---|---|---|---
20 | $3.614\times 10^{-1}$ | $2.3908\times 10^{-3}$ | $1.6088\times 10^{-4}$ | $3.1425\times 10^{-4}$
40 | $2.6467\times 10^{-2}$ | $8.3461\times 10^{-3}$ | $3.8203\times 10^{-5}$ | $7.8461\times 10^{-5}$
$5.9363\times 10^{-6}$ | $1.6361\times 10^{-5}$
$1.9566\times 10^{-6}$ | $4.6673\times 10^{-6}$
80 | $1.3204\times 10^{-2}$ | $2.4451\times 10^{-4}$ | $6.7665\times 10^{-6}$ | $1.5000\times 10^{-6}$
$2.0704\times 10^{-6}$ | $3.9415\times 10^{-6}$
$3.6112\times 10^{-7}$ | $8.1338\times 10^{-7}$
160 | $6.6847\times 10^{-3}$ | $6.6573\times 10^{-4}$ | $2.2361\times 10^{-6}$ | $3.7708\times 10^{-6}$
$3.8562\times 10^{-7}$ | $7.5454\times 10^{-7}$
$2.0934\times 10^{-7}$ | $2.5683\times 10^{-7}$
320 | $3.3393\times 10^{-3}$ | $1.5672\times 10^{-5}$ | $7.9409\times 10^{-7}$ | $9.3940\times 10^{-7}$
$1.8957\times 10^{-7}$ | $2.2149\times 10^{-7}$
5$.0883\times 10^{-8}$ | $5.9343\times 10^{-8}$

In Table 4, we show some numerical results. We can see that the the numerical
results of the two-grid method performs much better compared with the immersed
boundary (IB) method [26] and the IIM even with small $r=2,4,8$. The errors in
the coarse grid using the fourth-order method are comparable with that at the
fine grid using the second order IIM. From the result using a uniform mesh, we
observe fourth order convergence using the two-grid method when $r=4$.

######  Example 3.2.

An example of discontinuous coefficient and solution, and a complicated
interface.

This example is from [14, 16] where the augmented IIM was developed for
elliptic interface problems with piecewise constant coefficients. We consider
the elliptic interface problem $(\kappa u_{x})_{x}+(\kappa u_{y})_{y}=f(x,y)$
in the domain $(-1,1)\times(-1,1)$ with a flower shaped interface given by the
following equation in the polar coordinates

$\displaystyle r(\theta)=0.5+0.1\sin(8\theta),\quad 0\leq\theta<2\pi.$ (56)

The coefficient $\kappa$ is a piecewise constant: $\kappa^{-}$ inside the
interface and $\kappa^{+}$ outside. The right hand side $f$, the Dirichlet
boundary conditions, and the jump conditions are all determined from the exact
solution, see also Figure 8,

$\displaystyle
u(x,y)=\left\\{\begin{array}[]{ll}\displaystyle\frac{r^{2}}{\kappa^{-}}&\textrm{if
$(x,y)\in\Omega^{-},$}\\\ \vskip 5.0pt\cr\displaystyle\frac{r^{4}-0.1\log
2r}{\kappa^{+}}&\textrm{if $(x,y)\in\Omega^{+}$}.\end{array}\right.$ (59)

Table 5: Convergence results of two-grid method for Example 3.2. Left:
$\kappa^{-}=1$, $\kappa^{+}=10$; Right: $\kappa^{-}=50$, $\kappa^{+}=1$.

$N$ |  $\|E\|_{\infty}$ (coarse ) |  $\|E\|_{\infty}$ (fine)
---|---|---
40 | $2.6514\times 10^{-3}$ | $3.5781\times 10^{-4}$
$4.6060\times 10^{-4}$ | $6.3633\times 10^{-4}$
$4.1335\times 10^{-5}$ | 6.4903 $\times 10^{-5}$
80 | $4.9900\times 10^{-6}$ | $2.1077\times 10^{-5}$
$6.8967\times 10^{-6}$ | $2.3883\times 10^{-6}$
$2.3110\times 10^{-7}$ | $6.7858\times 10^{-7}$
160 | $9.8969\times 10^{-7}$ | $2.3783\times 10^{-6}$
$3.3377\times 10^{-7}$ | $6.7788\times 10^{-7}$
$7.6309\times 10^{-8}$ | $1.7808\times 10^{-7}$
320 | $2.8677\times 10^{-7}$ | $6.8180\times 10^{-7}$
$8.4336\times 10^{-8}$ | $1.7984\times 10^{-7}$
$2.1948\times 10^{-8}$ | $4.3751\times 10^{-8}$

$N$ |  $\|E\|_{\infty}$ (coarse ) |  $\|E\|_{\infty}$ (fine)
---|---|---
40 | $2.9005\times 10^{-3}$ | $5.1821\times 10^{-3}$
$2.4027\times 10^{-4}$ | $4.7331\times 10^{-4}$
$4.1345\times 10^{-5}$ | $8.0766\times 10^{-5}$
80 | $2.5074\times 10^{-4}$ | $4.7357\times 10^{-4}$
$4.3436\times 10^{-5}$ | $8.1008\times 10^{-5}$
$2.3557\times 10^{-6}$ | $4.9060\times 10^{-6}$
160 | $4.5306\times 10^{-5}$ | $8.1051\times 10^{-5}$
$2.4412\times 10^{-7}$ | $4.9550\times 10^{-6}$
$3.8420\times 10^{-7}$ | $1.1772\times 10^{-6}$
320 | $2.483\times 10^{-6}$ | $4.9601\times 10^{-6}$
$4.1551\times 10^{-7}$ | $1.1829\times 10^{-6}$
$1.5622\times 10^{-7}$ | $2.7301\times 10^{-7}$

Figure 8: The solution plot of the flower example with $N=40$, $r=2$, and
$\lambda=2$. Left plot: $u(x,y)$; Right plot: error plot with
$\|E\|_{\infty}=4.7331\times 10^{-4}$. We use blue dots to represent the
values at the coarse grid, and red dots for those at the fine grid.

This is a general example in which both the coefficients and the solution are
discontinuous. The interface is complicated with larger curvature at some part
of the interface, see Figure 8. In the figure, we use blue dots to represent
the values at the coarse grid, and red dots for those at the fine grid. In
Table 5, we show the numerical results with $r=2,4,8$. We have the same
convergence properties as before, that is, compared with the results using a
uniform mesh, we obtain roughly fourth order convergence when $r=4$. Note also
that the result of the two-grid method with $r=8$ has a fourth order
convergence compared with that obtained with $r=2$.

###  3.4 An internal layer example

Now we apply the developed two-grid method to an internal layer problem in
which the solution to the Poisson equation is

$\displaystyle
u(x,y)=\arctan\left(\frac{\sqrt{x^{2}+y^{2}+3/4}-1}{\epsilon}\right),\quad\quad-1<x,\;y<1.$
(60)

When $\epsilon$ is small, there is an internal layer centered at the circle
$x^{2}+y^{2}=1/4$. The source term below is obtained from Matlab symbolic
ToolBox,

$\begin{array}[]{rcl}&&\displaystyle
f(x,y)=\frac{1}{\epsilon}\left(\frac{2}{(\sigma_{3}^{2}+1)\sqrt{\sigma_{4}}}-\frac{2(x^{2}+y^{2})}{\sigma_{1}}-\frac{200(x^{2}+y^{2})\sigma_{3}}{\sigma_{2}}\right),\;\;\mbox{with}\\\
\vskip
5.0pt\cr&&\\!\\!\\!\displaystyle\sigma_{1}=(\sigma_{3}^{2}+1)\sigma_{4}^{3/2},\;\sigma_{2}=(\sigma_{3}^{2}+1)^{2}\sigma_{4},\;\sigma_{3}=\frac{1}{\epsilon}\left(\sqrt{\sigma_{4}}-1\right),\;\sigma_{4}=x^{2}+y^{2}+\frac{3}{4}.\end{array}$
(61)

The solution and the source term are displayed in Figure 9 at the coarse grid
(blue dots) and the fine grid (red dots) with $\epsilon=0.01$, $N=40$, $r=4$,
and $\lambda=2$. An internal boundary layer centered along the circle
$x^{2}+y^{2}=1/4$ can be seen clearly.

Figure 9: The solution plot of the internal layer example with
$\epsilon=0.01$, $N=40$, $r=4$, and $\lambda=2$. Left plot: $-u(x,y)$; Right
plot: $f(x,y)$ which has large magnitude and variations near the internal
layer around $x^{2}+y^{2}=1/4$. We use blue dots to represent the values at
the coarse grid, and red dots for those at the fine grid.

For this problem, we can apply the standard fourth order compact scheme both
in the fine and coarse meshes since there is no interface which does not
capture the internal layer well unless we use an unfordable fine mesh. The
developed two-grid method worked well for the internal layer problem. In Table
6, we show some numerical methods with different $N$ and $r$. In this example,
the largest errors have mostly appeared at hanging nodes since they have the
largest local truncation errors. We have the same convergence properties as
before, that is, from the uniform grid $(r=1)$ to that obtained from $r=4$, or
from $r=2$ to $r=8$, we obtain roughly fourth order convergence. The increased
computational complexity only happened in the refined tube.

Table 6: Convergence results of two-grid method for the internal layer example with $\epsilon=0.01$, and $r=2,4,8$.  $N$ |  $\|E\|_{\infty}$ (coarse ) |  $\|E\|_{\infty}$ (fine)
---|---|---
20 | $2.7557$ | $3.8289$
$3.0755\times 10^{-1}$ | $4.9217\times 10^{-1}$
$2.8506\times 10^{-3}$ | $8.0663\times 10^{-3}$
40 | $3.3413\times 10^{-1}$ | $4.9067\times 10^{-1}$
$5.3367\times 10^{-3}$ | $1.0304\times 10^{-2}$
$3.4447\times 10^{-5}$ | $2.0639\times 10^{-4}$
80 | $7.4351\times 10^{-3}$ | $7.4351\times 10^{-2}$
$1.3716\times 10^{-4}$ | $2.5445\times 10^{-4}$
$4.4450\times 10^{-6}$ | $1.4286\times 10^{-5}$
160 | $3.9407\times 10^{-3}$ | $4.5279\times 10^{-3}$
$1.4769\times 10^{-4}$ | $1.8625\times 10^{-4}$
$1.2278\times 10^{-6}$ | $1.6815\times 10^{-6}$

In Figure 9, we plot the computed solution (left) and the source term (right)
using blue dots to represent the values at the coarse grid, and red dots for
those at the fine grid. The internal layer is caused from large variations in
the source term near the center of the internal layer $x^{2}+y^{2}=1/4$.

##  4 Conclusions and Acknowledgements

In this paper, we have developed two-grid methods for some elliptic interface
and internal layer problems in one and two dimensions. The purpose is focused
on the accuracy of the solutions near the interface or internal layers so that
the accuracy is comparable with the solution in the coarse grid obtained from
fourth order compact finite difference schemes. New high order finite
difference schemes have been developed at boarder grid points neighboring two
grids for one-dimensional and two-dimensional problems, and at hanging nodes
for two dimensional problems. The developed two-grid methods preserve the
discrete maximum principle from the sign property, which leads to the
convergence of the finite difference methods.

We would like to thank Dr. Danfu Han of Hangzhou Normal University for the
internal layer example. Zhilin Li was partially supported by a Simons grant
633724. K. Pan is supported by Science Challenge Project (No. TZ2016002), the
National Natural Science Foundation of China (No. 41874086), the Excellent
Youth Foundation of Hunan Province of China (No. 2018JJ1042).

## References

  * [1] D. M. Bedivan, _A two-grid method for solving elliptic problems with inhomogeneous boundary conditions_ , Comput. Math. Appl. 29 (1995), no. 6, 59–66.
  * [2] M. Berger and P. Colella, _Local adaptive mesh refinement for shock hydrodynamics_ , J. Comput. Phys. 82 (1989), 64–84.
  * [3] D. L. Brown, R. Cortez, and M. L. Minion, _Accurate projection methods for the incompressible Navier-Stokes equations_ , J. Comput. Phys. 168 (2001), 464.
  * [4] A. J. Chorin, _On the convergence of discrete approximations of the Navier-Stokes equations_ , Math. Comp. 23 (1969), 341–353.
  * [5] R. Cortez and D. Varela, _Accurate projection methods for the incompressible Navier-Stokes equations_ , J. Comput. Phys. 138 (1997), 224–247.
  * [6] Binbin Du, Jianguo Huang, and Haibiao Zheng, _Two-grid Arrow-Hurwicz methods for the steady incompressible Navier-Stokes equations_ , J. Sci. Comput. 89 (2021), no. 1, Paper No. 24.
  * [7] Qiwei Feng, Bin Han, and Peter Minev, _Sixth order compact finite difference scheme for poisson interface problem with singular sources_ , Computers & Mathematics with Applications 99 (2021), 2–25.
  * [8] Bosco García-Archilla, _Shishkin mesh simulation: a new stabilization technique for convection-diffusion problems_ , Comput. Methods Appl. Mech. Engrg. 256 (2013), 1–16.
  * [9] B. E. Griffith, _A comparison of two adaptive versions of the immersed boundary method_ , Technical report, New York University (2009).
  * [10] B. E. Griffith, R. D. Hornung, D. M. Mcqueen, and C. S. Peskin, _An adaptive, formally second order accurate version of the immersed boundary method_ , J. Comput. Phys. 223 (2007), 10–49.
  * [11] Alan F. Hegarty, John J. H. Miller, Eugene O’Riordan, and G. I. Shishkin, _Special meshes for finite difference approximations to an advection-diffusion equation with parabolic layers_ , J. Comput. Phys. 117 (1995), no. 1, 47–54.
  * [12] Natalia Kopteva and Eugene O’Riordan, _Shishkin meshes in the numerical solution of singularly perturbed differential equations_ , Int. J. Numer. Anal. Model. 7 (2010), no. 3, 393–415.
  * [13] R. J. LeVeque and Z. Li, _The immersed interface method for elliptic equations with discontinuous coefficients and singular sources_ , SIAM J. Numer. Anal. 31 (1994), 1019–1044.
  * [14] Z. Li, _A fast iterative algorithm for elliptic interface problems_ , SIAM J. Numer. Anal. 35 (1998), 230–254.
  * [15] Z. Li and K. Ito, _Maximum principle preserving schemes for interface problems with discontinuous coefficients_ , SIAM J. Sci. Comput. 23 (2001), 1225–1242.
  * [16] , _The immersed interface method – numerical solutions of pdes involving interfaces and irregular domains_ , SIAM Frontier Series in Applied mathematics, FR33, 2006.
  * [17] Z. Li, Z. Qiao, and T. Tang, _An introduction to finite difference and finite element methods for ODE/PDEs of boundary value problems_ , Cambridge University Press, 2017.
  * [18] Z. Li and P. Song, _An adaptive mesh refinement strategy for immersed boundary/interface methods_ , Commun. Comput. Phys. 12 (2012), 515–527.
  * [19] , _Adaptive mesh refinement techniques for the immersed interface method applied to flow problems_ , Computers and Structures 122 (2013), 249–258.
  * [20] Z. Li, L. Wang, E. Aspinwall, R Cooper, P Kuberry, A Sanders, and K Zeng, _Some new analysis results for a class of interface problems_ , Mathematical Methods in the Applied Sciences 38(18) (2015), 4530–4539.
  * [21] Zhilin Li and Kejia Pan, _Can 4th-order compact schemes exist for flux type bcs?_ , submitted, 2021.
  * [22] William F. Mitchell, _A collection of 2D elliptic problems for testing adaptive grid refinement algorithms_ , Appl. Math. Comput. 220 (2013), 350–364.
  * [23] William F. Mitchell and Marjorie A. McClain, _A survey of $hp$-adaptive strategies for elliptic partial differential equations_, Recent advances in computational and applied mathematics, Springer, Dordrecht, 2011, pp. 227–258.
  * [24] K. W. Morton and D. F. Mayers, _Numerical solution of partial differential equations_ , Cambridge press, 1995.
  * [25] Kejia Pan, Dongdong He, and Zhilin Li, _A high order compact FD framework for elliptic BVPs involving singular sources, interfaces, and irregular domains_ , Journal of Scientific Computing 88 (2021).
  * [26] C. S. Peskin, _The immersed boundary method_ , Acta Numerica 11 (2002), 479–517.
  * [27] A. Roma, _A multi-level self adaptive version of the immersed boundary method_ , Ph.D. thesis, New York University, 1996.
  * [28] A. Roma, C. S. Peskin, and M. Berger, _An adaptive version of the immersed boundary method_ , J. Comput. Phys. 153 (1999), 509–534.
  * [29] J. C. Strikwerda, _Finite difference scheme and partial differential equations_ , Wadsworth & Brooks, 1989.
  * [30] S. V. Tikhovskaya, _A two-grid method for an elliptic equation with boundary layers on a Shishkin mesh_ , Lobachevskii J. Math. 35 (2014), no. 4, 409–415.
  * [31] Xin Tong, Lipo Wang, Jinqiang Chen, and Hua Ouyang, _A Cartesian grid method with improvement of resolving the boundary layer structure for two-dimensional incompressible flows_ , Internat. J. Numer. Methods Fluids 93 (2021), no. 8, 2637–2659.
  * [32] Yang Wang, Yanping Chen, and Yunqing Huang, _A two-grid method for semi-linear elliptic interface problems by partially penalized immersed finite element methods_ , Math. Comput. Simulation 169 (2020), 1–15.
  * [33] Y. Xie and W. Ying, _A fourth-order kernel-free boundary integral method for implicitly defined surfaces in three space dimensions_ , J. Comput. Phys. 415 (2020).
  * [34] Jinchao Xu, _Two-grid discretization techniques for linear and nonlinear PDEs_ , SIAM J. Numer. Anal. 33 (1996), no. 5, 1759–1777.
  * [35] A. I. Zadorin, S. V. Tikhovskaya, and N. A. Zadorin, _A two-grid method for elliptic problem with boundary layers_ , Appl. Numer. Math. 93 (2015), 270–278.