A Fast Iterative Algorithm for Elliptic Interface Problems.
A fast, second-order accurate iterative method is proposed for the elliptic equation \[ \grad\cdot(\beta(x,y) \grad u) =f(x,y) \] in a rectangular region $\Omega$ in two-space dimensions. We assume that there is an irregular interface across which the coefficient $\beta$, the solution u and its derivatives, and/or the source term f may have jumps. We are especially interested in the cases where the coefficients $\beta$ are piecewise constant and the  jump in $\beta$ is large. The interface may or may not align with an underlying Cartesian grid. The idea in our approach is to precondition the differential equation before applying the immersed interface method proposed by LeVeque and Li [ SIAM J. Numer. Anal., 4 (1994), pp. 1019--1044]. In order to take advantage of fast Poisson solvers on a rectangular region, an intermediate unknown function, the jump in the normal derivative across the interface,  is introduced. Our discretization is equivalent to using a second-order difference scheme for a corresponding Poisson equation in the region, and a second-order discretization for a Neumann-like interface condition. Thus second-order accuracy is guaranteed. A GMRES iteration is employed to solve the Schur complement system derived from the discretization.  A new weighted least squares method is also proposed to approximate interface quantities from a grid function.  Numerical experiments are provided and analyzed. The number of iterations in solving the Schur complement system appears to be independent of both the jump in the coefficient and the mesh size.
Introduction
. Consider the elliptic equation
r (fi(x;
2\Omega
Given BC on
in a rectangular
domain\Omega in two space dimensions. Within the region, suppose there
is an irregular interface \Gamma across which the coefficient fi is discontinuous. Referring to
Fig 1, we assume that fi(x; y) has a constant value in each sub-domain,
The interface \Gamma may or may not align with a underline Cartesian grid.
Depending on the properties of the source term f(x; y), we usually have jump
conditions across the interface \Gamma:
This work was supported by URI grant #N00014092-J-1890 from ARPA, NSF Grant DMS-
9303404, and DOE Grant DE-FG06-93ER25181.
y Department of Mathematics, University of California at Los Angeles, Los Angeles, CA 90095.
(zhilin@math.ucla.edu).
Z. LI
(a)
-0.4
Interface
(b)
Fig. 1. Two typical computational domains and interfaces with uniform Cartesian grids.
where (X(s); Y (s)) is the arc-length parameterization of \Gamma, the superscripts \Gamma or
denotes the limiting values of a function from one side or the other of the interface.
These two jump conditions can be either obtained by physical reasoning or derived
from the differential equation, see [2, 9, 14] etc. Note in potential theory, v(s) 6j 0
corresponds to a single layer source along the interface \Gamma, while w(s) 6j 0 corresponds
to a double layer source. The normal derivative u n usually has a kink across the
interface due to the discontinuity in the coefficient fi. If w(s) 6j 0; then the solution
would be discontinuous across the interface.
There are many applications in solving elliptic equations with discontinuous coef-
ficients, for example, steady state heat diffusion or electrostatic problems, multi-phase
and porous flow, solidification problems, and bubble computations etc. There are two
main concerns in solving (1.1)-(1.4) numerically:
ffl How to discretize (1.1)-(1.4) to certain accuracy. It is difficult to study the
consistency and the stability of a numerical scheme because of the discontinuities
across the interface.
ffl How to solve the resulting linear system efficiently. Usually if the jump in the
coefficient is large, then the resulting linear system is ill-conditioned, and the
number of iterations in solving such a linear system is large and proportional
to the jump in the coefficient.
There are a few numerical methods designed to solve elliptic equations with discontinuous
coefficients, for example, harmonic averaging, smoothing method, and finite
element approach etc., see [2] for a brief review of different methods. Most of these
methods can be second order accurate in the l-1 or the l-2 norm, but not in the l-1
norm, since they may smooth out the solution near the interface.
A. Mayo and A. Greenbaum [14, 15] have derived an integral equation for elliptic
interface problems with piecewise coefficients. By solving the integral equation, they
can solve such interface problems to second order accuracy in the l-1 norm using
the techniques developed by A. Mayo in [13, 14] for solving Poisson and biharmonic
equations on irregular regions. The total cost includes solving the integral equation
and a regular Poisson equation using a fast solver, so this gives a fast algorithm. The
possibility of extension to variable coefficients is mentioned in [14].
R.J. LeVeque and Z. Li have recently developed a different approach for discretizing
A FAST ALGORITHM FOR INTERFACE PROBLEMS 3
elliptic problems with irregular interfaces called the immersed interface method (IIM)
[2, 9], which can handle both discontinuous coefficients and singular sources. This
approach has also been applied to three dimensional elliptic equations [7], parabolic
wave equations with discontinuous coefficients [4, 5],
and the incompressible Stokes flow problems with moving interfaces [3, 6]. L. Adams [1]
has successfully implemented a multi-grid algorithm for the immersed interface method.
However, there are some numerical examples with large jumps in the coefficients in
which the immersed interface method may fail to give accurate answers or converge
very slowly.
In this paper, we propose a fast algorithm for elliptic equations with large jumps
in the coefficients. The idea is to precondition the elliptic equation before using the
immersed interface method. In order to take advantage of fast Poisson solvers on
rectangular regions, we introduce an intermediate unknown function [u n ](s) which is
defined only on the interface. Then we discretize a corresponding Poisson equation,
which has different sources from the original one, using the standard five-point stencil
with some modification in the right hand side. Our discretization is equivalent to using
a second order difference scheme to approximate the Poisson equation in the interior
a second order discretization for the Neumann-like interface
condition
Thus from the error analysis for elliptic equations with Neumann boundary conditions,
for example, see [17], we would have second order accurate solution at all grid points
including those near or on the interface. A GMRES method is employed to solve
the Schur complement system derived from the discretization. A new weighted least
squares method is proposed to approximate interface quantities such as u \Sigma
n from a
grid function defined on the entire domain. This new technique has been successfully
applied in the multi-grid method for interpolating the grid function between different
levels [1] with remarkable improvement in the computed solution. These ideas will be
discussed in detail in the following sections. The method described in this paper seems
to be very promising not only because it is second order accurate, but also because the
number of iterations in solving the Schur complement system is almost independent
of both the jumps in the coefficients and the mesh size. This has been observed from
our numerous numerical experiments, though we have not been able to prove this
theoretically. Our new method has been used successfully for the computation of some
inverse problems [20].
This paper is organized as follows. In Section 2, we precondition (1.1)-(1.4) to
get an equivalent problem. In Section 3, we use the IIM idea to discretize the equivalent
problem and derive the Schur complement system. The weighted least squares
approach to approximate u \Sigma
n from the grid function u ij is discussed in Section 4. Some
implementation details are addressed in Section 5. Brief convergence analysis is given
in Section 6. An efficient preconditioner for the Schur complement system is proposed
in Section 7. Numerical experiments and analysis can be found in Section 8. Some
new approaches in the error analysis involving interfaces are also introduced there.
2. Preconditioning the PDE to an equivalent problem. The problem we
intend to solve is the following:
4 Z. LI
Problem (I).
r (fi(x;
(2.6a)
Given BC on
with specified jump conditions along the interface \Gamma
(2.7a)
Consider the solution set u g (x; y) of the following problem as a functional of g(s).
Problem (II).
f
if x
f
if x
(2.8a)
Given BC on
with specified jump conditions 1
(2.9a)
Let the solution of Problem (I) be u   (x; y), and define
along the interface \Gamma. Then u   (x; y) satisfies the elliptic equation (2.8a)-(2.8b) and
jump conditions (2.9a)-(2.9b) with g(s) j g   . In other words, u g   (x; y) j u   (x; y),
and
@n
is satisfied. Therefore, solving Problem (I) is equivalent to finding the corresponding
and then u g   (x; y) in Problem (II). Notice that g   is only defined along the
interface, so it is one dimension lower than u(x; y): Problem (II) is an elliptic interface
problem which is much easier to solve because the jump condition [u n ] is given instead
of [fiu n ]. With the immersed interface method, it is very easier to construct a second
order scheme which also satisfies the conditions of the maximum principle. In this
paper, we suppose fi is piecewise constant as in (1.2), so Problem (II) is a Poisson
equation with a discontinuous source term and given jump conditions. We can then
use the standard five-point stencil to discretize the left hand side of (2.8a), but modify
the right hand side to get a second order scheme, see [2, 9] for the detail. Thus we
can take advantage of fast Poisson solvers for the discrete system. The cost in solving
1 The jump conditions (2.9a) and (2.9b) depend on the singularities of the source term f(x; y)
along the interface. However, in the expression of (2.8a), we do not need information of f(x; y) on
the interface \Gamma, so there is no need to write f(x; y) differently.
A FAST ALGORITHM FOR INTERFACE PROBLEMS 5
Problem (II) is just a little more than that in solving a regular Poisson equation on the
rectangle with a smooth solution. For more general variable coefficient, the discussions
in this paper are still valid except we can not use a fast Poisson solver because of the
convection term (rfi \Delta ru)=fi in (2.8a). However a multi-grid approach developed by
L. Adams [1] perhaps can be used to solve Problem (II).
We wish to find numerical methods with which we can compute u g   (x; y) to second
order accuracy. We also hope that the total cost in computing g   and u g   is less
than that in computing u g   through the original Problem (I). The key to success is
computing g   efficiently. Below we begin to describe our method to solve g   . Once
is found, we just need one more fast Poisson solver to get the solution u   (x; y).
3. Discretization. We use a uniform grid on the rectangle [a; b] \Theta [c; d] where
the Problem (I) is defined:
We assume that simplicity. We use a cubic spline ~
passing through a number of control points (X
express the immersed interface, where s is the arc-length of the interface and (X
the position of the k-th point on the interface \Gamma. Other representations of the interface
are possible. A level set formulation is currently under investigation.
Any other quantity q(s) defined on the interface such as w(s) and g(s) can also
be expressed in terms of a cubic spline with the same parameter s. Since cubic splines
are twice differentiable we can gain access to the value of q(s) and its first or second
derivatives at any point on the interface in a continuous manner.
We use upper case letters to indicates the solution of the discrete problem and
lower case letters for the continuous solutions.
Given W k and G k , the discrete form of jump conditions (2.9a) and (2.9b), with
the immersed interface method, the discrete form of (2.8a) can be written as
where
is the discrete Laplacian operator using the standard five-point stencil. Note that if
happens to be on the interface, then f ij =fi ij is defined as the limiting value
from a pre-chosen side of the interface. C ij is zero except at those irregular grid
points where the interface cuts through the five-point stencil. A fast Poisson solver,
for example, FFT, ADI, Cyclic reduction, or Multi-grid, can be applied to solve (3.12).
The solution U ij depends on G k , W k , continuously. In matrix and vector
form we have
is the discrete linear system for the Poisson equation when W k , G k
are all zero. The solution is smooth for such a system. B(W;G) is a mapping from
(3.12). From [2, 9] we
6 Z. LI
know that B(W;G) depends on the first and second derivatives of w(s), and the first
derivative of g(s), where the differentiation is carried out along the interface. At this
stage we do not know whether such a mapping is linear or not. However in the discrete
case, all the derivatives are obtained by differentiating the corresponding splines which
are linear combination of the values on those control points. Therefore B(W;G) is
indeed linear function of W and G and can be written as
are two matrices with real entries. Thus (3.13) becomes
The solution U of the equation above certainly depends on G and we are interested in
finding G   which satisfies the discrete form of (2.7b)
where the components of the vectors U
are discrete approximation of the
normal derivative at control points from each side of the interface. In the next section,
we will discuss how to use the known jump G, and sometimes also V , to interpolate
U ij to get U
n in detail. As we will see in the next section, U
depend
on U , G and V linearly
where E, D, and -
are some matrices and
P . Combine (3.14) and (3.16) to
obtain the linear system of equations for U and G:
G
F
The solution U and G are the discrete forms of u g   (x; y) and g   , the solution of
Problem (II) which satisfies (2.11).
The next question is how to solve (3.17) efficiently. The GMRES method applied
to (3.17) directly or the multi-grid approach [1] are two attractive choices. However,
in order to take advantage of fast Poisson solvers, we have decided to solve G in (3.17)
first, and then to find the solution U by using one more fast Poisson solver. Eliminating
U from (3.17) gives a linear system for G
(D
F
This is an n b \Theta n b system for G, a much smaller linear system compared to the one
for U . The coefficient matrix is the Schur complement of D in (3.17). In practice,
the matrices A, B, E, D, P , and the vectors -
F are never formed. The matrix and
vector form are merely for theoretical purposes. Thus an iterative method, such as
the GMRES iteration [18], is preferred. The way we compute (3.16) will dramatically
change the condition number of (3.18).
A FAST ALGORITHM FOR INTERFACE PROBLEMS 7
4. A weighted least squares approach for computing interface quantities
from a grid function. When we apply the GMRES method to solve the Schur
complement system of (3.18) for G   , we need to compute the matrix-vector multiplica-
tion, which is equivalent to computing U \Gamma
n with the knowledge of U ij and the
jump condition [U n ]. This turns out to be a crucial step in solving the linear system
(3.18) for G   . Our approach is based on a weighted least squares formulation. The
idea described here can also be, and has been, applied to the case, where we want to
approximate some quantities on the interface from a grid function. For example, interpolating
U ij to the interface to get U \Gamma (X; Y ) or U
on the interface. This new approach has also been successfully applied to the multi-grid
method for interpolating the grid function between different levels by L. Adams[1]
with remarkable improvement in the computed solution.
We start from the continuous situation, the discrete version can be obtained ac-
cordingly. Let u(x; y) be a piecewise smooth function, with discontinuities only along
the interface. We want to interpolate u(x get approximations to the normal
derivatives
are only defined on the interface, to second
order.
Our approach is inspired by Peskin's method in interpolating a velocity field u(x; y)
to get the velocity of the interface using a discrete delta function. The continuous and
discrete forms are the following:
ZZ\Omega
where ~
is a discrete delta function. A commonly used one is
Notice that ffi h (x) is a smooth function of x. Peskin's approach is very robust and only
a few neighboring grid points near ~
are involved. However this approach is only first
order accurate and may smear out the solution near the interface.
Our interpolation formula for
n , for example, can be written in the following form
(j ~
where d ff (r) is a function of the distance measured from the point ~
X,
d ff
Q is a correction term which can be determined once fl ij are known. Although we are
trying to approximate the normal derivatives here, the same principle also applies to
the function values U as well with different choices of fl ij and Q. Note no extra effort
is needed to decide which grid points should be involved. Therefore, expression (4.21)
is robust and depends on the the grid function u ij continuously, two very attractive
8 Z. LI
properties of Peskin's formula (4.20). In addition to the advantages of Peskin's ap-
proach, we also have flexibility in choosing the coefficients fl ij and the correction term
Q to achieve second order accuracy. The parameter ff in (4.21) can be fixed or chosen
according to problems, see Section 8.
Below we discuss how to use the immersed interface method to determine the
coefficients and the correction term Q. They are different from point to point on
the interface. So they should really be labeled as fl ij; ~
etc. But for simplicity of
notation we will concentrate on a single point ~
drop the subscript ~
X.
Since the jump condition is given in the normal direction, we introduce local
coordinates at (X; Y ),
where ' is the angle between the x-axis and the normal direction. Under such new
coordinates, the interface can be parameterized by -(j); j. Note that
and, provided the interface is smooth at (X; Y well. The solution of
the Poisson equation Problem (II) will satisfies the following interface relations, see
[2, 9] for the derivation,
be the -j coordinates of
where the sign depends on whether (- lies on the side of \Gamma.
Using Taylor expansion of (4.25) about (X; Y ) in the new coordinates, after collecting
terms we have
(j ~
a 9
where the a j are given by
A FAST ALGORITHM FOR INTERFACE PROBLEMS 9
(j ~
a
(j ~
a
(j ~
a
(j ~
a
(j ~
a
(j ~
(j ~
a
(j ~
a
(j ~
a
(j ~
(j ~
a
(j ~
From the interface relations (4.24) we know that all the jumps in the expression above
can be expressed in terms of the known information. Since
- , we obtain the
linear system of equations for the coefficients
a 9
Note that we would have the exact same equation if we want to interpolate a smooth
function to get an approximation u n at ~
X to second order accuracy. The discontinuities
across the interface only contribute to the correction term Q. This agrees
with our analysis in [2, 9] for Poisson equations with discontinuous and/or singular
sources, where we can still use the classical five-point scheme but add a correction
term to the right hand side at irregular grid points.
If the linear system (4.26) has a solution, then we can obtain a second order
approximation to the normal derivative
n by choosing an appropriate correction term
Q. Therefore we want to choose ff big enough, say ff - 1:6h, such that at least six
grid points are involved. Usually we have an under-determined linear system which
has infinitely many solutions. We should then choose the one fl
ij with the minimal
2-norm
subject to (4:26):
For such a solution, each fl
ij will have roughly the same magnitude O(1=h); so
ij d ff (j ~
is roughly a decreasing function of the distance measured from ~
X.
This is one of desired properties of our interpolation. In practice, only a hand full of
grid points, controlled by the parameter ff, are involved. Those grid points which are
closer to (X; Y ) have more influence than others which are further away.
Z. LI
When we know the coefficients we also know the a k 's. From the a k 's and the
interface relations (4.24), we can determine the correction term Q easily,
Thus we are able to compute
n to second order accuracy. We can derive a formula
for
n in exactly the same way. However, with the relation u
g, we can write
down a second order interpolation scheme for
immediately
(j ~
is the solution we computed for
n . In the next section, we will explain an
important modification of either (4.21) or (4.28) depends on the magnitude of fi \Gamma or
We should mention another intuitive approach, one-sided interpolation, in which
we only use grid points on the proper side of the interface in computing a limiting
value at the interface:
This approach does not make use of the interface relations (4.24), so we have to have
at least six points from each side in order to have a second order scheme. Note that
we can also use the least squares technique described in this section for one-sided
interpolation. This approach has been tested already. The weighted least squares
approach using the interface relations (4.24) appears superior in practice. It has the
following advantages:
ffl Fewer grid points are involved. When we make use of the interface relation,
compared to the one-sided interpolation, the number of grid points which are
involved is reduced roughly by half.
ffl Second order accuracy with smaller error constant. The grid points involved
in our approach are clustered around the point (X; Y ) on the interface, and
those which are closer to (X; Y ) have more influence than those which are
further away in our weighted least squares approach. We have smaller error
constant in the Taylor expansions compared to the one-sided interpolation.
The error constant can be as much as 8 - 27 times smaller as the one-sided
interpolation. In two dimensional computation, we can not take m and n to
be very large, to have a smaller error constant sometimes is as important as
to have a high order accurate method.
ffl Robust and smoother error distribution. We have a robust way in choosing
the grid points which are involved. The interpolation formulas (4.21) and
depend continuously on the location (X; Y ) and the grid points
and so does the truncation error for these two interpolation schemes. In other
words, we will have a smooth error distribution. This is very important in
moving interface problems where we do no want to introduce any non-physical
oscillations.
A FAST ALGORITHM FOR INTERFACE PROBLEMS 11
downs. In one-sided interpolation, sometimes we can not find enough
grid points in one particular side of the interface, then the one-sided interpolation
will break down. In our approach, every grid point on one side is
connected to the other by the interface relations (4.24). So no break down will
occur.
ffl Trade off or disadvantages. The only trade off of our weighted least squares
approach is that we have to solve a under-determined 6 by p linear system of
equation (where p - instead of solving one that is 6 by 6. The larger ff is,
the more computational cost in solving (4.26). Fortunately, the linear system
has full row rank and can be solved by the LR-RU method [8] or other efficient
least squares solvers.
5. Some details in implementation. The main process of our algorithm is to
solve the Schur complement system (3.18) using the GMRES method with an initial
guess
G (0)
We need to derive the right hand side, and explain how to compute the matrix-vector
multiplication of the system without explicitly forming the coefficient matrix. The
right hand side needs to be computed just once which is described below.
5.1. Computing the right hand side of the Schur complement system.
If we take apply one step of the immersed interface method to solve
Problem (II) to get U(0), then
With the knowledge of U(0) and G = 0; we can compute the normal derivatives on
each side of the interface to get U \Sigma
using the approach described in the previous
section. Thus the right hand side of the Schur complement system is
F
The last two equalities are obtained from (3.16) and (3.18) with Now we are
able to compute the right hand side of the Schur complement system.
5.2. Computing the matrix-vector multiplication of the Schur comple-
ment. Now consider the matrix-vector multiplication
(D
of the Schur complement, where Q is an arbitrary vector of dimension n b . This involves
essentially two steps.
1. A fast Poisson solver for computing
which is the solution of Problem (II) with
Z. LI
2. The weighted least squares interpolation to compute U
The residual vector in the flux jump condition is
which is the same residual vector of the second equation in (3.17) from our definition.
In other words, see also (3.16)
The matrix-vector multiplication (5.29) then can be computed from the last equality
of the following derivation:
F from (5.30)
V from (5.32):
Note that from the second line to the third line we have used the following
which is defined in (3.18).
It worth to point out that once our algorithm is successfully terminated, which
means that the residual vector is close to the zero vector, we not only have an approximation
Q to the solution G   , an approximation U(Q) to the solution U , bult also
approximations U \Sigma
n (Q) to the normal derivatives from each side of the interface. The
normal derivative information is very useful for some moving interface problems where
the velocity of the interface depends on the normal derivative of the pressure.
6. Convergence Analysis. As to this point, we have had a complete algorithm
for solving the original elliptic equations of the form Problem (I). We have transformed
the original elliptic equation to a corresponding Poisson equation with different source
term and jump conditions, or internal boundary conditions, (2.9b) and (2.11). The
jump condition (2.11) is Neumann-like boundary condition which involves the normal
derivatives from each side of the interface. In our algorithm, the classical five-point
difference scheme at regular grid points is used. This discetization is second order
accurate. As discussed in Section 4, the Neumann-like internal boundary condition
(2.11) is also discretized to second order. So from the analysis in Chapter 6 of [17] on
Neumann conditions, we should be able to conclude second order convergence globally
for our computed solution, provide that we can solve the Poisson Problem (II) to second
order accuracy. This is confirmed in our early work [2, 9]. Numerical experiments have
confirmed second order accuracy of the computed solution for numerous test problems,
see Section 8.
7. An efficient preconditioner for the Schur complement system. With
the algorithm described in previous sections, we are able to solve Problem I to second
order accuracy. In each iteration we need to solve a Poisson equation with a
modified right hand side. A fast Poisson solver such as a fast Fourier transformation
method etc. [19], can be used. The number of iterations of
A FAST ALGORITHM FOR INTERFACE PROBLEMS 13
the GMRES method depends on the condition number of the Schur complement. If
we make use of both (4.21) and (4.28) to compute U \Sigma
n , the condition number seems
to be proportional to 1=h. Therefore the number of iterations will grow linearly as we
increase the number of grid points.
Below we propose a modification in the way of computing U \Sigma
which seems to
improve the condition number of the Schur complement system dramatically.
If
n and
n are the exact solutions, that is
then we can solve
n or u
n in terms of v, It is easy to get
or
The idea is simple and intuitive. We use one of the formulas (4.21) or (4.28) obtained
from the weighted least squares interpolation to approximate
n or u
n , and then use
or (7.33) to approximate u
n or
n to force the solution to satisfy the flux jump
condition. This is actually an acceleration process, or a preconditioner for the Schur
complement system (3.18).
With this modification, the number of iterations for solving the Schur complement
system seems to be independent of the mesh size h, and almost independent of the
jump [fi] in the coefficient as well, see the next section for more details. Although
we have not been able to prove this claim, the algorithm seems to be extraordinary
successful.
Whether we use the pair (4.21) and (7.34) or the other (4.28) and (7.33) have only
a little affect on the accuracy of the computed solutions and the number of iterations.
The algorithm otherwise behaves the same and the analysis in the next section seems
to be true no mater what pair we choose.
We have been using the following criteria to choose the desired pair
n is determined by (4.28)
n is determined by (4.21)
which seems always better than the choice of the other way around.
8. Numerical Experiments. We have done many numerical experiments with
different interfaces and various test functions. Since our scheme can handle jumps in
the solution, we have great flexibility in choosing test problems. From the numerical
tests we intend to determine:
ffl The accuracy of computed solutions. Are they second order accurate?
14 Z. LI
ffl The numbers of iterations as we change the mesh size h and the ratio of the
discontinuous coefficient,
ffl The ability of the algorithm to deal with complicated interfaces and large
jumps in the coefficient.
All the experiments are computed with double precision. The computational parameters
include:
Computational rectangle [a; b] \Theta [c; d].
ffl The number of grid points m and n in the x- and y- directions respectively,
we assume that is the mesh size.
ffl The number of control points n b . The interface is expressed in terms of cubic
splines passing through the control points.
ffl The parameter ff in the weighted least squares interpolation. We take
specified differently.
The maximum norm is used to measure the errors in the computed solution U ij ,
and the normal derivatives U \Gamma
n p from each side of the interface at the p-th
control point. The relative errors are defined as follows
where ~
is one of control points on the interface. Each grid point is labeled
as either in the
or the
outside\Omega + of the interface and the exact solution is
determined accordingly. In other words, the exact solution is not determined from the
exact interface but the discrete one. In Table 1, r i , 3 is the ratio of successive
errors. A ratio of 4 corresponds to second-order accuracy. In Table 1, k is the number
of iterations required in solving the Schur complement system (3.18). The ratio of
coefficients is defined . In the figures, we use S to express the
slopes of least squares line of experimental data (log(h i ); log(E i )).
Example 1. Consider the following interface
where
within the computational domain Fig 1 shows some interfaces with
different parameters r 0 , . Dirichlet boundary conditions, as well as
the jump conditions [u] and [fiu n ] along the interface, are determined from the exact
A FAST ALGORITHM FOR INTERFACE PROBLEMS 15
solution
if (x; y)
r 4
if (x; y)
. The source term can be determined accordingly:
if (x; y)
if (x; y)
We provide numerical results for three typical cases below.
Case A. The interface parameters are chosen as r
the interface is a circle centered at the origin, see Fig 2(a). With C
the solution is continuous everywhere, but u n and fiu n are discontinuous across the
circle. It is easy to verify that [fiu n when we take C Fig 3(a) is the
plot of the solution \Gammau with
Case B. The interface parameters are chosen as r
20,
Fig 2(a). We shift the center of the interface a little bit to have a non-symmetric
solution relative to the grid. We want our test problems to be as general
as possible. The interface is irregular but the curvature has modest magnitude. So
with a reasonable number of points on the interface, we can express it well. Now it is
almost impossible to find an exact solution which is continuous but not smooth across
the interface, so we simply set C Fig 3(b) is the plot of the
solution \Gammau with
Case C. The interface parameters are chosen as r
20,
Fig 2(b). Now the magnitude of the curvature is very large at some
points on the interface and we have to have enough control points to resolve it. The
solution parameters are set to be the same as in Case B.
Fig 4-6 and Table 1 are some plots and data from the computed solutions which
we will analyze below.
8.1. Accuracy. Table 1 shows the results of grid refinement analysis for Case
A with two very different ratio 0:5, the ratio r i
are very close to 4 indicating second order convergence. With
the error in the solution drops much more rapidly. This is because the solution in
approaches a constant as fi becomes large, and it is quadratic
order accurate method would give high accurate solution in both regions. So it is not
surprising to see the ratio r 1 is much larger than 4. For the normal derivatives, we
expect second order accuracy again since fi
n is not quadratic and has magnitude of
O(1). This agrees with the results r 2 and r 3 in Table 1.
In Fig 4 we consider the opposite case when fi . In this case
the solution is not quadratic so we see the expected second order accuracy. Fig 4(a)
Z. LI
(a)
-0.4
A
(b)
-0.4
Fig. 2. Different interfaces of Example 1. (a) Case A and B. (b) Case C.
(a)
-0.4
-0.3
-0.2
-0.4
-0.3
-0.2
-0.1Fig. 3. The solutions \Gammau of Example 1 with 1. (a) Case A, a circular
interface where the solution is continuous but [fiun Case B, an irregular interface where
both the solution and the flux [fiun ] are discontinuous.

Table
Numerical results and convergence analysis for Case A with
A FAST ALGORITHM FOR INTERFACE PROBLEMS 17
(a)
(b)
-5
log(h)
log(E
Fig. 4. (a): Error distribution for Case A. (b): Errors E i vs the mesh size h in log-log scale for
Case A with
(a)
log(h)
log(E
(b)
-6.4 -6.2 -6 -5.8 -5.6 -5.4 -5.2 -5 -4.8
-5
log(h)
log(E)
Fig. 5. Errors E i vs the mesh size h in log-log scale for Case B with
The solid line: n
. The dotted line: n (b) The solid line and
dotted line are the same as in (a) but on a different scale. The dash-dotted line is the result obtained
with
log(h)
log(E
2:71 2:07
Fig. 6. Errors E i vs the mesh size h in log-log scale for Case C with 1. The
solid line: fixed n b (n 520). The dotted line: n
Z. LI
plots the error distribution over the region. The error seems to change continuously
even though the maximum error occurs on or near the interface. Usually if the curvature
is very big in some part of an interface, for example, near a corner, then we would
observe large errors over the neighborhood of that part of the interface.
For interface problems, the errors usually do not decrease monotonously as we
refine the grid unless the interface is aligned with one of the axes. We need to study the
asymptotic convergence rate which is usually defined as the slope of the least squares
line of the experimental data (log(h i ); log(E i )). Fig 4(b) plots the errors versus the
mesh size h in log-log scale for the case n. The asymptotic convergence rate is
about 2:62 compared to 2 for a second order method. As h gets smaller we can see
the curves for the errors become flatter indicating the asymptotic convergence rate will
approach 2.
The dotted curves in Fig 5 and Fig 6 are the results for case B and C, where
the interfaces are more complicated compared to case A. Again we take
The asymptotic convergence rates are far more than two. Such behavior can also be
observed from Example 4 in [16]. Does it mean that our method is better than second
order? Certainly this is not true from our discretization. Below we explain what is
happening.
For interface problems, the errors depend on the solution u(x; y), the mesh size h,
the jump in the coefficient [fi], the interface \Gamma and its relative position to the grid, and
the maximal distance between control points on the interface, h b . We can write the
error in the solution, for example, as follows
The first term in the right hand side of (8.4) is the error from the discretization of
the differential equation. The term C (u; h; h b ; [fi]; \Gamma) has magnitude O(1) but does
not necessarily approach to a constant. The second term in the right hand side of
is the error from the discretization of the interface \Gamma. If we use a cubic spline
interpolation, then q ? 2. For Case A, the interface is well expressed and the first term
in (8.4) is dominant, we have clearly second order convergence. For Case B and C,
the interfaces are more complicated and the second term in (8.4) is dominant. That is
why we have higher than second order convergence. Eventually, the error in the first
term will dominate and we will then observe second order convergence.
To further verify the arguments above, we did some tests with fixed number of
control points n b . For example, we take n the solid line in
Fig 6. Presumably the interface is expressed well enough and the second term in (8.4)
is negligible. We see the slopes of the least squares line of the errors E 1 and E 2 are
2:15 and 2:07 respectively indicating second order convergence. Usually
the error in the normal derivatives
n and u
behaves the same, so we only need to
study one of them. If we let n b change with the same speed as the number of grids m
and n, then the second term in (8.4) is dominant and the slopes of the least squares
line of the errors E 1 and E 2 are 2:71 and 2:69 respectively. Once n b is large enough,
the first term will dominate in (8.4) and the error will decrease quadratically. This
can also be seen roughly from Fig 6. Note that the errors oscillate as n gets large
whether we fix n b or not. But the fluctuation becomes smaller as we refine the grid.
The upper envelop of E 1 behaves the same as the least squares line of the experimental
data (log(h i ); log(E i )). So it is reasonable to use the asymptotic convergence rate to
discuss the accuracy when errors do not behave monotonously.
A FAST ALGORITHM FOR INTERFACE PROBLEMS 19
As another test, we let n b change slower than m and n. The solid lines in Fig 5(a)
are obtained with
Now we have roughly
and the errors decrease quadratically with the mesh size h, but not h b . The
slopes of the least squares line of the errors E 1 and E 2 are 2:23 and 2:22 respectively.
We now discuss the effect of the different choice of the parameter ff in the least
squares approximation described in Section 4 on the solution. Most of the computations
are done with Fig 5(b), the dash-dotted line where
As we can expect, the smaller ff is, the higher accuracy in the computed solution because
the points involved are clustered together and the error in the Taylor expansion
will be smaller. However, the smaller ff is, the more oscillatory in the error as we
refine the grid. For larger ff, the computation cost increases quickly, but the error
behaves much smoother with the mesh size h. Usually we can take small ff for smooth
interfaces, and larger ff if we want a smoother error distribution for more complicated
interfaces.
8.2. The number of iterations versus the mesh size h. Fig 7(a), also see
Fig 10(a) for Example 2, shows the number of iterations versus the number of grid
points m and n for case A, B, and C. It is not surprising to see that the number of
iterations depends on the shape of the interface. The number of iterations required for
Case C is larger than that for Case A and B. But it is wonderful to see that the number
of iterations is almost independent of the mesh size h. For Case A, where the interface
is a circle, we only need about iterations for all choices of the mesh size h for two
extreme cases We will see in the next paragraph that this is also
true for different choices of the ratio . Note that the number of iterations is
about two or three fewer than the numbers of calls of the fast Poisson solver. We need
two or three of them for initial set up of the Schur complement system. In Fig 7(a),
the lowest curve corresponds to case A with the lowest but
the second curve corresponds to 1. For case B, the number of
iterations required is about 17 - 21 for respectively. For case C,
the most complicated interface, the number of iterations is about 46 with reasonable
number of control points on the interface for
8.3. The number of iterations versus the jump ratio Fig 7(b),
also see Fig 10(b) for Example 2 , plots the number of iterations versus the jump ratio
ae in log-log scale with fixed number of grids goes away from
the unit we have larger jump relatively in the coefficient. The number of iterations
increases proportional to jlog(ae)j when ae is small but soon reaches a point after which
the number of iterations will remain as a constant. Such points depends on the shape
of the interface. For Case A, it requires only about 5 - 6 iterations at the most for
iterations for ae ? 1 in solving the Schur complement system
using the GMRES method. For Case B, the numbers are about 17 - 22. For Case C,
the most complicated interface in our examples, the numbers are about 47 - 69. As
we mentioned in the previous paragraph, also see Fig 7(a), for Case C, with only 160
control points we can not express the complicated interface Fig 2(b) very well. If we
take more control points on the interface, then the number of iterations will be about
Z. LI
(a)
A
(b)
A
6 91769
Fig. 7. The number of iterations for Example 1 with vs the number of grids
n. Case A: lower curve, Case B: lower curve,
the ratio of
jumps in log-log scale with
when
Example 2. This geometry of this example is adapted from Problem 3 of [1]. The
solution domain is the [\Gamma1; 1] \Theta [0; 3] rectangle and the interface is determined by
Fig 8(a) show the solution domain and the interface \Gamma with x
Again Dirichlet boundary condition, as well as the jump conditions [u] and [fiu n ] are
determined from the exact solution
The source term can be determined accordingly:
Fig 8(b) is a plot of the computed solution.
This example is different from Example 1 in several ways. The solution is independent
of the coefficient fi. But the magnitude of the jump [fiu n ] and the source
increase with the magnitude of the jump [fi]. However we have observed
similar behaviors in the numerical results as we discussed in Example 1. Example 1
and Example 2 are two extreme samples of elliptic interface problems. So we should
be able to get some insights about the method proposed in this paper.
Fig 9 shows errors E i versus mesh size h in log-log scale with different choice of
b . In Fig 9(a), . The solid lines correspond to a fixed discretization
of the interface, n As we expected, the asymptotic convergence rate for
are 2:1272. They are all close to 2 indicating
A FAST ALGORITHM FOR INTERFACE PROBLEMS 21
(a)
(b)0.51.52.5-11
y
x
Fig. 8. (a) The interface of Example 2. (b) The solution of Example 2.
second order accuracy. The dotted line in Fig 9(a) correspond to a variable n b which
changes in the same rate as the number of grid point m in x-direction. The asymptotic
convergence rate of E i for are S 3:3473. They are
Fig 9(b).
These numbers are all larger than 2 similar to the cases we saw in Fig 5, and Fig 6.
We have explained such phenomena already.
Fig 10(a) plots the number of iteration versus the number of grids n with
Again we consider two extreme cases,
with Once the interface is well expressed somewhere after n ? 180, the
number of iteration will slightly decrease to a constant which is about 28 for
and 34 for Fig 10(b) plots the number of iteration versus the ratio ae with
fixed grid
ae ? 1. We observe the same behavior as in Fig 7(b). Initially the number of iterations
increases proportional to jlog(ae)j as ae goes away from the unit, but it soon approaches
a constant which is about 28 for ae ! 1 and 34 for ae ? 1.
(a)
-6.2 -6 -5.8 -5.6 -5.4 -5.2 -5 -4.8 -4.6
log(h)
3:35
2:14
(b)
-6 -5.9 -5.8 -5.7 -5.6 -5.5 -5.4 -5.3 -5.2 -5.1
log(E
Fig. 9. Errors E i vs the mesh size h in log-log scale for Example 2 with
1: The solid line: fixed n b , dotted line: n
Summary of the numerical experiments. In our computations, the largest
error usually occurs at those points which are close to the part of the interface which
22 Z. LI
(a)
(b)
)28Fig. 10. The number of iterations for Example 2 with vs the number of
grids n. Lower curve: vs the ratio
of jumps in log-log scale with fixed grid
has large curvature. Depending on the shape of the interface, we should take enough
control points on the interface so the error in expressing the interface does not dominate
the global error. However, once such a critical number is decided, we do not need to
double it as we double the number of grid points, which saves some computational cost.
We should still be able to maintain second order accuracy. The number of iteration for
solving the Schur complement system using a GMRES method is almost independent
of both the mesh size h as well as the jump in the coefficient.
9. Conclusions. We have developed a second-order accurate fast algorithm for
a type of elliptic interface problems with large jumps in the coefficient across some
irregular interface. We precondition the original partial differential equation to obtain
an equivalent Poisson problem with different source terms and a Neumann-like interface
condition. The fast Poisson solver proposed in [2, 9] can be employed to solve
the Schur complement system for the intermediate unknown, the jump in the normal
derivative along the interface. Then we proposed a preconditioning technique for the
Schur complement system which seems to be very successful. Numerical tests revealed
that the number of iterations in solving the Schur complement system is independent
of both the mesh size h and the jump in the coefficient, though we have not proved this
strictly in theory. The idea introduced in this paper might be applicable to other related
problem, for example, to domain decomposition techniques. A new least squares
approach to approximate interface quantities from a grid function is also proposed. By
analyzing the numerical experiments, we have discussed some issues in error analysis
involving interfaces.
There is still a lot of room for improving the method described in this paper. For
example, we have used cubic spline interpolations for closed interfaces. There are some
advantages of this approach. But large errors can occur at the connection of the first
and the last control points when we try to make the curve closed. That might also
be one of reasons why the error does not decrease monotonously. As an alternative, a
level set formulation is under investigation.
The next project following this paper is to study the case with variable coefficients.
A FAST ALGORITHM FOR INTERFACE PROBLEMS 23
We can rewrite (1.1) either as (2.8a) or
r
if x
r
f
if x
are the averages of the coefficients fi from each side of the interface.
Whether (2.8a) or (9.5) is used, we shall still introduce an intermediate unknown, the
jump in the normal derivative across the interface if the jump condition is given in the
form of [fiu n ]. In this way, the coefficients of the difference scheme would be very close
to those obtained form the classical five-point stencil. We can not take advantage of
the fast Poisson solvers for variable coefficient anymore, but we can make use of the
multi-grid method developed by L. Adams in [1].
10.

Acknowledgments

. It is my pleasure to acknowledge the encouragements
and advice from various people including Prof. Randy LeVeque, Stanley Osher, Tony
Chan, Loyce Adams, Jun Zou and Barry Merriman. Thanks also to Prof. Yousef Saad
and Dr. Victor Eijkhout for helping me to implement and understand the GMRES
method.



--R

A multigrid algorithm for immersed interface problems.
The immersed interface method for elliptic equations with discontinuous coefficients and singular sources.
Simulation of bubbles in creeping flow using the immersed interface method.
Immersed interface methods for wave equations with discontinuous coefficients.
Finite difference methods for wave equations with discontinuous coefficients.
Immersed interface method for Stokes flow with elastic boundaries or surface tension.
A note on immersed interface methods for three dimensional elliptic equations.
Uniform treatment of linear systems - algorithm and numerical stability
The Immersed Interface Method - A Numerical Approach for Partial Differential Equations with Interfaces
Immersed interface method for moving interface problems.
ADI methods for heat equations with discontinuties along an arbitrary interface.
On the rapid evaluation of heat potentials on general regions.
The fast solution of Poisson's and the biharmonic equations on irregular regions.
The rapid evaluation of Volume
Fast parallel iterative solution of Poisson's and the biharmonic equations on irregular regions.
A fast poisson solver for complex geometries.
Numerical Solution of Partial Differential Equations.
GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems.
Fast Poisson solver.
Computing some inverse problems.
--TR

--CTR
Kazufumi Ito , Zhilin Li, Solving a Nonlinear Problem in Magneto-Rheological Fluids Using the Immersed Interface Method, Journal of Scientific Computing, v.19 n.1-3, p.253-266, December
Songming Hou , Xu-Dong Liu, A numerical method for solving variable coefficient elliptic equation with interfaces, Journal of Computational Physics, v.202 n.2, p.411-445, 20 January 2005
M. Oevermann , R. Klein, A Cartesian grid finite volume method for elliptic equations with variable coefficients and embedded interfaces, Journal of Computational Physics, v.219 n.2, p.749-769, December, 2006
Petter Andreas Berthelsen, A decomposed immersed interface method for variable coefficient elliptic equations with non-smooth and discontinuous solutions, Journal of Computational Physics, v.197 n.1, p.364-386, 10 June 2004
Shaozhong Deng , Kazufumi Ito , Zhilin Li, Three-dimensional elliptic solvers for interface problems and applications, Journal of Computational Physics, v.184 n.1, p.215-243,
Peter Schwartz , Michael Barad , Phillip Colella , Terry Ligocki, A Cartesian grid embedded boundary method for the heat equation and Poisson's equation in three dimensions, Journal of Computational Physics, v.211 n.2, p.531-550, 20 January 2006
Xu-Dong Liu , Thomas C. Sideris, Convergence of the ghost fluid method for elliptic equations with interfaces, Mathematics of Computation, v.72 n.244, p.1731-1746, October
Shi Jin , Xuelei Wang, Robust numerical simulation of porosity evolution in chemical vapor infiltration: II. Two-dimensional anisotropic fronts, Journal of Computational Physics, v.179 n.2, p.557-577, July 2002
Ming-Chih Lai , Zhilin Li , Xiaobiao Lin, Fast solvers for 3D Poisson equations involving interfaces in a finite or the infinite domain, Journal of Computational and Applied Mathematics, v.191 n.1, p.106-125, 15 June 2006
Do Wan Kim , Young-Cheol Yoon , Wing Kam Liu , Ted Belytschko, Extrinsic meshfree approximation using asymptotic expansion for interfacial discontinuity of derivative, Journal of Computational Physics, v.221 n.1, p.370-394, January, 2007
I. Klapper , T. Shaw, A large jump asymptotic framework for solving elliptic and parabolic equations with interfaces and strong coefficient discontinuities, Applied Numerical Mathematics, v.57 n.5-7, p.657-671, May, 2007
Carlos J. Garca-Cervera , Zydrunas Gimbutas , Weinan E., Accurate numerical methods for micromagnetics simulations with general geometries, Journal of Computational Physics, v.184 n.1, p.37-52,
B. P. Lamichhane , B. I. Wohlmuth, Mortar finite elements for interface problems, Computing, v.72 n.3-4, p.333-348, May 2004
Chohong Min , Frdric Gibou , Hector D. Ceniceros, A supra-convergent finite difference scheme for the variable coefficient Poisson equation on non-graded grids, Journal of Computational Physics, v.218 n.1, p.123-140, 10 October 2006
Frederic Gibou , Ronald P. Fedkiw , Li-Tien Cheng , Myungjoo Kang, A second-order-accurate symmetric discretization of the Poisson equation on irregular domains, Journal of Computational Physics, v.176 n.1, p.205-227, February 10, 2002
John K. Hunter , Zhilin Li , Hongkai Zhao, Reactive autophobic spreading of drops, Journal of Computational Physics, v.183 n.2, p.335-366, December 10
Y. C. Zhou , G. W. Wei, On the fictitious-domain and interpolation formulations of the matched interface and boundary (MIB) method, Journal of Computational Physics, v.219 n.1, p.228-246, 20 November 2006
Y. C. Zhou , Shan Zhao , Michael Feig , G. W. Wei, High order matched interface and boundary method for elliptic equations with discotinuous coefficients and singular sources, Journal of Computational Physics, v.213 n.1, p.1-30, 20 March 2006
Xiaolin Zhong, A new high-order immersed interface method for solving elliptic equations with imbedded interface of discontinuity, Journal of Computational Physics, v.225 n.1, p.1066-1099, July, 2007
George Biros , Lexing Ying , Denis Zorin, A fast solver for the Stokes equations with distributed forces in complex geometries, Journal of Computational Physics, v.193 n.1, p.317-348, January 2004
Sining Yu , Yongcheng Zhou , G. W. Wei, Matched interface and boundary (MIB) method for elliptic problems with sharp-edged interfaces, Journal of Computational Physics, v.224 n.2, p.729-756, June, 2007
GMRES method;immersed interface method;discontinuous coefficients;elliptic equation;schur complement;preconditioning;cartesian grid