Special Issue Paper Received 20 February 2013
Published online in Wiley Online Library
(wileyonlinelibrary.com) DOI: 10.1002/mma.2865 MOS subject classification: 65M12; 65M20
Some new analysis results for a class of interface problems Zhilin Lia,b * † , Li Wangc , Eric Aspinwalld , Racheal Coopere , Paul Kuberryf , Ashley Sandersg and Ke Zengh Communicated by Q. Wang Interface problems modeled by differential equations have many applications in mathematical biology, fluid mechanics, material sciences, and many other areas. Typically, interface problems are characterized by discontinuities in the coefficients and/or the Dirac delta function singularities in the source term. Because of these irregularities, solutions to the differential equations are not smooth or discontinuous. In this paper, some new results on the jump conditions of the solution across the interface are derived using the distribution theory and the theory of weak solutions. Some theoretical results on the boundary singularity in which the singular delta function is at the boundary are obtained. Finally, the proof of the convergency of the immersed boundary (IB) method is presented. The IB method is shown to be first-order convergent in L1 norm. Copyright © 2013 John Wiley & Sons, Ltd. Keywords: immersed boundary (IB) method; immersed interface method (IIM); jump conditions; discontinuous coefficient; Dirac delta function; weak solution; boundary singularity; convergence of IB method; equivalent boundary conditions
1. Introduction Interface problems arise in many areas of science and mathematics. Whenever there are different materials such as oil and water, gas and water, or bubbles in a fluid, we have interface problems. Another kind of interface problems is the Peskin’s immersed boundary (IB) method (see, e.g., [1] for a review and references therein). The IB method is not only a mathematical modeling tool in which a complicated boundary condition can be treated as a source distribution but also a numerical method in which a discrete delta function is used. The IB method is robust, simple, and has been applied to many problems in mathematical biology, fluid mechanics, material science, and other areas. However, the IB method is only first-order accurate in general. The immersed interface method (IIM) [2–4] is a second-order accurate method particularly designed for interface problems. In the IIM, jump conditions are enforced near the interface to obtain high-order accuracy. Knowing the jump conditions is crucial in the IIM. Mathematically, there are two characteristics in an interface problem. One is that a discontinuity exists in the coefficient of the differential equation along a co-dimensional space. The other one is singular source terms such as the Dirac delta function or its derivatives along a co-dimensional space. The following equation is a model interface problem in one dimension (1D), .ˇ.x/u0 /0 D f .x/ C c ı.x ˛/ C cN ı 0 .x ˛/,
0 < x < 1,
0 < ˛ < 1,
(1)
where ˇ.x/ D
ˇ .x/ ˇ C .x/
if x < ˛, if x > ˛.
(2)
a Department of Mathematics, North Carolina State University, Raleigh, NC 27695, USA b School of Mathematical Sciences, Nanjing Normal University, Nanjing, China c School of Mathematical Sciences and Jiangsu Key Laboratory for NSLSCS, Nanjing Normal University, Nanjing, China d Department of Mathematics, Florida State University, Tallahassee, FL 32306, USA e Department of Mathematical Sciences, Virginia Commonwealth University, Richmond, VA 23284, USA f Department of Mathematical Sciences, Clemson University, Clemson, SC 29634, USA g Department of Mathematics, Jackson State University, Jackson, MS 39217, USA h Department of Mathematics, Zhejiang University, Hangzhou, Zhejiang Province, 310058, China *Correspondence to: Zhilin Li, Department of Mathematics, North Carolina State University, Raleigh, NC 27695, USA. † E-mail: [email protected]
Copyright © 2013 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2013
Z. LI ET AL. along with two boundary conditions at x D 0 and x D 1. In the IB method developed by Peskin [1, 5], the Dirac delta function is intensively used. Often, a delta function singularity corresponds to a source distribution, whereas the derivative of the delta function singularity corresponds to a dipole distribution. For Stokes and Navier–Stokes equations, a delta function distribution along a surface indicates that the derivative of the velocity has jump discontinuities across the interface, whereas the derivative of the delta function indicates that the pressure has a jump discontinuity. There are three goals in this paper. The first one is to derive the jump conditions for (1). The second one is to consider the limiting case as ˛ approaches one of boundary points. The third one is to provide a convergence proof for the IB method. Because there are no other regularities except at the interface ˛, we know that the solution is piecewise smooth excluding the interface ˛, that is, u.x/ 2 C 2 .0, 1/ n ˛. The discontinuity is at the interface ˛. We can derive the jump condition of the solution and its derivative across the interface. We define ˇ D ˇ.˛/ D
lim
x!˛,x˛
ˇ.x/,
and Œˇ D ˇ C ˇ .
(3)
Jump conditions are crucial in the immersed interface method [2–4]. If cN D 0 or ˇ is continuous at ˛, that is, ˇ D ˇ C , then the jump conditions have been derived (e.g., [4]). We will derive the more challenging case when cN ¤ 0 and ˇ ¤ ˇ C . To the best of our knowledge, we have not seen the jump conditions and the derivation in the literature for the difficult case. We also report theoretical results and the derivation for the boundary singularity problem for which ˛ is at one of the boundary points. We will see that the effect of the delta function singularity can be absorbed by suitable boundary conditions. In this paper, we will use the notations u0 .x/ D ux .x/ D du dx .x/ whenever it is convenient. The rest of the paper is organized as follows. In the next section, we define the weak solution and derive the jump conditions. In Section 3, we discuss the case when the delta function singularity appears at one of the boundaries. We show that it can be translated as a contribution from the boundary conditions. In Section 4, we show the convergence of the IB method. Our result confirms that the IB method is first-order accurate in L1 in general. We conclude in the last section.
2. Deriving the jump conditions To use the IIM, we must know or derive jump conditions. The jump conditions describe the discontinuity of the solution and its derivative. The information is necessary to obtain second-order accuracy of the finite difference scheme as seen in [3]. We first state one of the main results of this paper in the following theorem. Theorem 2.1 Assume ˇ.x/ 2 C 1 Œ0, ˛/ [ C 1 .˛, 1, f .x/ 2 C.0, 1/. Then, the following jump conditions hold: Œu D
2Nc , ˇ C ˇC
Œˇux D c.
(4)
We prove the theorem by considering the following cases although the last case can cover all the situations. The purpose is to show that for simpler cases, the proofs are also simpler. 2.1. The singularity from ı.x ˛/ only In this simplest case, we assume that Œˇ D 0, cN D 0. Thus, the ODE is simply .ˇux /x D f .x/ C cı.x ˛/. From the differential equation theory (e.g., [6]), we know that the solution is continuous, that is, Œu D 0. Integrating the differential equation from ˛ to ˛ C , we have Z
Z
˛C ˛
.ˇux /x dx D
Z
˛C
˛
f .x/dx C
˛C
˛
cı.x ˛/dx.
(5)
Notice that f .x/ is continuous in the neighborhood of ˛, and consequently, the first integration in the right-hand side vanishes. Also, from the definition of the Dirac delta function, we have the following property of the delta function: Z
C
g.x/ı.x ˛/dx D g.˛/
(6)
for any continuous function g.x/. Thus, from (5), we have ˇ C uC x .˛/ ˇ ux .˛/ D c,
or Œˇux D c.
Thus, the theorem is true for this case. 2.2. Singularities from ı.x ˛/ and ı 0 .x ˛/ In this case, we assume that Œˇ D 0 and both c and cN can be non-zero. Thus, the ODE is .ˇux /x D f .x/Cc ı.x ˛/C cN ı 0 .x ˛/. The presence of the derivative of the delta function corresponds to a double layer in the potential theory. Thus, theRsolution typically has a finite jump across the interface ˛. We cannot simply integrate the differential equation from ˛ to ˛ C because cN ı 0 .x ˛/dx does not exist. Copyright © 2013 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2013
Z. LI ET AL. Similar to the process in deriving the weak form in the finite element method, we multiply a testing function .x/ 2 C01 .0, 1/ to both sides of the differential equation and then integrate to obtain Z
Z
˛C ˛
.ˇux /x .x/dx D
Z
˛C
˛
f .x/.x/dx C
Z
˛C
˛
cı.x ˛/.x/dx C
˛C
˛
cN ı 0 .x ˛/.x/dx
(7)
D c.˛/ cN 0 .˛/. The right-hand side is obtained from the definition of the weak derivatives of the delta function. We manipulate the left-hand side of (7) as follows: ˇ˛ C Z ˇ LHS D .ˇux / ˇ ˛
˛C
˛
D Œˇux .˛/ ˇ.˛/
ˇux 0 .x/dx ˇ˛ C Z ˇ .u /ˇ 0
˛
˛C
˛
!
(8)
00
u .x/dx
D Œˇux .˛/ ˇ.˛/ Œu 0 .˛/. Because .x/ is arbitrary and Equation (7) is the same as (8), we have to have ˇ.˛/ Œu D cN ,
Œˇux D c,
that is,
Œu D
cN 2Nc . D ˇ.˛/ ˇ C ˇC
Thus, the theorem is also true for this case. 2.3. Singularities from ı.x ˛/, ı 0 .x ˛/, and a discontinuity in ˇ Now the coefficient ˇ.x/ has a discontinuity at ˛. For simplicity of proof, we assume that ˇ is piecewise constant if x < ˛, ˇ ˇ.x/ D ˇ C if x > ˛.
(9)
Before we derive the jump conditions, we define a weak solution using the standard regularization technique. We will see that the solution defined from the jump conditions is the same as the weak solution. First, we define a regularized Heaviside function H .x/ as follows: 8 < 0 if x < , 1 (10) H .x/ D 1 C x C 1 sin x if jxj , 2 : 1 if x > . Note that H .0/ D 12 . We define the regularized coefficient as ˇ .x/ D ˇ C Œˇ H .x ˛/. Note that ˇ .˛/ D
1 2
ˇ C ˇ
C
(11)
. We define the regularized solution u .x/ as the solution of the following boundary value problem:
.ˇ .x/u0 /0 D f .x/ C c ı.x ˛/ C cN ı 0 .x ˛/,
0 < x < 1,
0 < ˛ < 1.
(12)
Because ˇ .x/ is continuous now, we have Œu D
2Nc , ˇ .˛C/ C ˇ .˛/
ˇ u0 D c.
(13)
Note that we have lim Œu D lim
!0
!0
2Nc 2Nc D C . ˇ .˛C/ C ˇ .˛/ ˇ C ˇ
(14)
Definition 2.2 The weak solution to (1) is defined as uw .x/ D lim u .x/. !0
(15)
It is obvious that the weak solution satisfies the jump conditions (4). We now prove that from the original differential equations, we can derive the same jump conditions as that of the weak solution. Copyright © 2013 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2013
Z. LI ET AL. 2.3.1. The derivation of the general jump conditions. As before, we multiply a testing function .x/ 2 C01 .0, 1/ to both sides of the differential equation and then integrate to obtain Z
Z
˛C ˛
.ˇux /x .x/dx D
Z
˛C
f .x/.x/dx C
˛
Z
˛C
cı.x ˛/.x/dx C
˛
˛C
˛
cN ı 0 .x ˛/.x/dx
(16)
0
D c.˛/ cN .˛/. The right-hand side is obtained from the definition of the weak derivatives. We manipulate the left-hand side of (16) as follows: ˇ˛ C Z ˇ LHS D .ˇux / ˇ ˛
Z D Œˇux .˛/
˛C ˛
˛C
˛
ˇ.x/ux 0 .x/dx
ˇ.x/ux 0 .x/dx.
We cannot take ˇ.x/ out of the integration anymore because it is discontinuous. We represent ˇ.x/ and u.x/ as follows (see also [7]): ˇ.x/ D ˇ C Œˇ H.x ˛/, 1 u.x/ D u C Œu H.x ˛/ C Œux H.x ˛/.x ˛/ C Œuxx H.x ˛/.x ˛/2 C . 2 Differentiating the aforementioned representation, we obtain ux .x/ D Œu ı.x ˛/ C Œux ı.x ˛/.x ˛/ C Œux H.x ˛/ C . Plugging the aforementioned expressions into the integral of LHS, we have Z
˛C ˛
ˇ.x/ux 0 .x/dx D
D
Z
˛C
Z
˛C
ˇ C Œˇ H.x ˛/ Œu ı.x ˛/ C Œux ı.x ˛/.x ˛/ ˛ C Œux H.x ˛/ C . 0 .x/dx
˛
C
ˇ Œuı.x ˛/ 0 .x/dx C
Z
˛C ˛
0
Z
˛C
˛
ŒˇŒuH.x ˛/ı.x ˛/ 0 .x/dx
ˇ Œux H.x ˛/ .x/dx C Z
Z
˛C
˛
ŒˇŒux H2 .x ˛/ 0 .x/dx C
d ŒˇŒuH2 .x ˛/ 0 .x/dx C 0 C 0 C ˛ dx Z ˛C 1 D ˇ Œu 0 .˛/ C ŒˇŒuH2 .x ˛/ 0 .˛/ ŒˇŒuH2 .x ˛/ 00 .x/dx 2 ˛ 1 D ˇ Œu 0 .˛/ C ŒˇŒu 0 .˛/ C 0 2 ˇ C ˇC Œu 0 .˛/. D 2 1 D ˇ Œu .˛/ C 2
0
˛C
Because .x/ 2 C01 is arbitrary, by matching the coefficients in front of .˛/ and 0 .˛/ in (16) and using the aforementioned relation, we have the jump conditions (4).
3. The boundary singularity analysis In this section, we discuss what happens when a delta function appears at one of the boundary points. We call such an irregularity a boundary singularity for short. This problem is worth discussing for several reasons. First of all, in some real applications, the interface is quite close to the boundary. Second, for many applications, sources may indeed appear on parts of boundaries. Third, it is theoretically significant to change the irregular problem with a ı function on the right-hand side to a regular partial differential equation in order to make it more convenient to solve. The discussions here have some similarities of the limiting value theory (e.g., [8]). Our discussions in this section indicate that we can translate the delta function singularity on a boundary into a corresponding boundary condition. Without loss of generality, we consider that the boundary singularity is at the boundary x D 0. In this case, the discontinuity in the coefficient is irrelevant. Note that it is rarely seen a dipole (derivative of the delta function) singularities at boundaries. From another point of view, the two pieces of the solutions are almost arbitrarily independent because of arbitrary parameters c and cN for the dipole case. Therefore, we only consider the case when a delta function singularity appears at one of boundary points. Thus, the problem can be written as .ˇ.x/u0 /0 D f .x/ C c ı.x/, Copyright © 2013 John Wiley & Sons, Ltd.
0 < x < 1.
(17) Math. Meth. Appl. Sci. 2013
Z. LI ET AL. For simplicity, we assume a Dirichlet boundary condition at x D 1 with u.1/ D and a linear boundary condition at x D 0, that is, a Dirichlet, or Neumann, or a mixed boundary condition. Note that the solution can be decomposed as the sum of two problems: 0 (18) ˇ.x/u01 D f .x/ 0 < x < 1, with homogenous but the same type of boundary conditions, and 0 ˇ.x/u02 D c ı.x/,
0 < x < 1,
(19)
with the original boundary conditions. Because u1 .x/ is the solution to a standard problem, we can treat it easily. Thus, we will simply discuss the second equation. To further simplify the discussion, we assume that ˇ D 1. The problem is simplified to u00 D c ı.x/,
0 < x < 1.
(20)
We first define the weak solution of (20). Definition 3.1 We consider u00˛ D c ı.x ˛/,
0 < x < 1,
0 < ˛ < 1,
(21)
with a boundary condition at x D 0 and x D 1. The weak solution of (20) is defined as u.x/ D lim˛!0 u˛ .x/. To consider the effect of the delta function singularity at one boundary, we use the following procedure: use the jump condition and the boundary condition to obtain the analytical solution for u˛ .x/; take the limit of the solution on the right interval .˛, 1 as ˛ approaches zero; and use the limiting solution to derive the equivalent equations and equivalent boundary condition.
3.1. Dirichlet boundary conditions In this case, the interface problem has the following form: uxx .x/ D cı.x ˛/, u.0/ D
0 < ˛ < 1,
u.1/ D
(22)
0 < ˛ < 1.
Theorem 3.2 The aforementioned limiting solution lim u˛ .x/ is the solution of the following problem: ˛!0
uxx D 0,
u.0/ D ,
u.1/ D .
(23)
That is, for a Dirichlet boundary condition, the boundary singularity has no effect on the solution. Proof We know that the solution has the following form: u˛ .x/ D
v1 C v2 x, v3 C v4 x,
0 < x < ˛, ˛ < x < 1.
(24)
Using the known jump conditions Œu D 0 and Œux D c and the boundary condition u.0/ D , we have the following: v1 D v1 C v2 ˛ D v3 C v4 ˛ v2 C c D v4 v3 C v4 D . Solving the linear system, we have v1 D v2 D c˛ C c v3 D c˛ v4 D c˛ C . Consequently, the analytical solution is as follows: C .c˛ C c/x, u˛ .x/ D c˛ C .c˛ C /x, Copyright © 2013 John Wiley & Sons, Ltd.
0 < x < ˛, ˛ < x < 1.
(25)
Math. Meth. Appl. Sci. 2013
Z. LI ET AL. When ˛ approaches zero, considering the analytical result on the right interval .˛, 1, the solution is as follows: uN .x/ D C . /x.
(26)
It is the solution of the equivalent regular problem: uxx .x/ D 0 u.0/ D
(27)
u.1/ D . This completes the proof.
3.2. Neumann boundary conditions We consider the problem with a Neumann boundary condition at x D 0. uxx .x/ D cı.x ˛/, ux .0/ D , u.1/ D
0 < ˛ < 1,
(28)
0 < ˛ < 1.
Theorem 3.3 The aforementioned limiting solution lim u˛ .x/ is the solution of the following problem: ˛!0
uxx D 0,
ux .0/ D C c,
u.1/ D .
(29)
In other words, the boundary source strength was added to the derivative boundary condition. Proof We know that the solution has the following form: u˛ .x/ D
v1 C v2 x, v3 C v4 x,
0 < x < ˛, ˛ < x < 1.
Using the known jump conditions Œu D 0 and Œux D c and the boundary condition ux .0/ D , we have the following: v2 D v1 C v2 ˛ D v3 C v4 ˛ v2 C c D v4 v3 C v4 D . Solving the linear system, we have v1 D C c˛ c v2 D v3 D c v4 D C c. Consequently, the analytical solution is as follows: C c˛ c C x, u˛ .x/ D c C . C c/x,
0 < x < ˛, ˛ < x < 1.
When ˛ approaches zero, considering the right interval .˛, 1, the approximate solution is as follows: uN .x/ D c C . C c/x. This is the weak solution of the equivalent equation for the Neumann boundary condition without the ı function: uxx .x/ D 0, ux .0/ D C c,
(30)
u.1/ D . That completes the proof. Copyright © 2013 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2013
Z. LI ET AL. 3.3. Mixed boundary conditions We now consider a mixed (Robin) boundary condition at x D 0 of the following problem: uxx .x/ D cı.x ˛/, u.0/ C ˇux .0/ D ,
(31)
u.1/ D 0 < ˛ < 1. Theorem 3.4 The aforementioned limiting solution lim u˛ .x/ is the solution of the following problem: ˛!0
uxx D 0,
u.0/ C ˇux .0/ D cˇ C ,
u.1/ D .
(32)
Again, we see that the source of the boundary singularity was added to the derivative part of the mixed boundary condition. Proof We know that the solution has the following form: u˛ .x/ D
v1 C v2 x, v3 C v4 x,
0 < x < ˛, ˛ < x < 1.
Using the known jump conditions Œu D 0 and Œux D c and the boundary condition u.0/ C ˇux .0/ D , we have the following: v1 C ˇv2 D v1 C v2 ˛ D v3 C v4 ˛ v2 C c D v4 v3 C v4 D . Solving the linear system, we have v1 D
ˇ.Cc˛c/ ˇ
v2 D
C.cc˛/ ˇ
v3 D
ˇ.c/Cc˛ ˇ
v4 D
Ccˇ .Cc˛/ . ˇ
Therefore, the analytical solution is as follows: 8 < ˇ.Cc˛c/ C C.cc˛/ x ˇ ˇ u˛ .x/ D : ˇ.c/Cc˛ C Ccˇ .Cc˛/ x ˇ
0 < x < ˛, ˛ < x < 1.
ˇ
When ˛ approaches zero, using the same technique, the limiting solution is as follows: uN .x/ D
ˇ. c/ C cˇ C x. ˇ ˇ
(33)
Because uN .0/Cˇ uNx .0/ D cˇ C , the equivalent regular differential equation for the mixed boundary condition without the ı function is the following: uxx .x/ D 0 u.0/ C ˇux .0/ D cˇ C ,
u.1/ D .
This completes the proof.
4. A proof of the convergence of the immersed boundary method It is well known that Peskin’s IB method is only first-order accurate in L1 norm. In [9], the author has proved the first-order accuracy of the IB method for the Stokes equations (see also [10]). However, there are few theoretical proofs on the IB method for elliptic interface problems. This is the main motivation of this section. We will give a proof for the 1D model, u00 D c ı.x ˛/ 0 < x < 1, Copyright © 2013 John Wiley & Sons, Ltd.
0 < ˛ < 1,
u.0/ D u.1/ D 0,
(34) Math. Meth. Appl. Sci. 2013
Z. LI ET AL. in this section. The proof for the 2D problem is more challenging and is under investigation. Note that the analytic solution to the aforementioned equation cx .1 ˛/ if x xk , (35) u.x/ D c˛ .1 x/ otherwise. Given a uniform Cartesian grid xi D ih, i D 0, 1, : : : , n, h D 1=n, the IB method leads to the following system of linear equations: Ui1 2Ui C UiC1 D cıh .xi ˛/ , h2
i D 1, 2, : : : , n 1,
(36)
where Ui is the finite difference approximation of the solution u.xi /, and ıh .xi ˛/ is a discrete delta function applied to the grid point xi . In the matrix-vector form, the aforementioned finite difference equations can be written as AU D F, where A is the tridiagonal matrix formed by the discrete Laplacian. It is well known that A is a symmetric positive definite matrix and diagonally dominant. Note that a discrete delta function has a compact support in the neighborhood of the interface, that is, ıh .x/ ¤ 0, only if jxj Wh, where W is a constant. Commonly used discrete delta functions include the hat discrete delta function (ı hat .x/ with W D 1) and the cosine discrete delta function (ı cosine .x/ with W D 2) (e.g., [4]). Note that when we use the hat delta function, the result is the same as that of the IIM for the simple model. The solution to the finite difference equations is the same as a true solution if there are no round-off errors, that is, u.xi1 / 2u.xi / C u.xiC1 / D cıhhat .xi ˛/ , h2
i D 1, 2, : : : , n 1
(37)
(see, e.g., [4, 11]). But this is not the case for other discrete delta functions. We define the error vector as E D fEi g, where Ei D u.xi / Ui . The local truncation error is defined as T D fTi g, Ti D
u.xi1 / 2u.xi / C u.xiC1 / cıh .xi ˛/ . h2
(38)
With the definition, we have AU D F and Au D F C T, and therefore, AE D T. For the hat discrete delta function, we have jTi j D 0 for all i’s. For other discrete delta functions, we have jTj j D O.1=h/ for a few grid points neighboring the interface ˛. So, the interesting question is as follows: Why is the IB method still first-order accurate, that is, kEk1 D O.h/? To answer this question, we first introduce the following lemma. Lemma 4.1 Let Ay D ek , where ek is the kth unit base vector, then yi D
hxi .1 xk / hxk .1 xi /
if i k, otherwise.
(39)
The significance of this lemma is that the solution is order h smaller than the concentrated source. Proof We note the following identity: Ay D h
ek D h ıhhat ek . h
(40)
For this simple case, the IB and IIM are identical. Thus, from the IIM [4, 11], we know that y is the exact discrete solution at the grid points of the following boundary value problem u00 D h ı.x xk /,
0 < x < 1,
u.0/ D u.1/ D 0,
(41)
whose solution is y.x/ D
hx .1 xk / hxk .1 x/
if x xk , otherwise.
This completes the proof.
(42)
Note that jy.x/j h. From this lemma, we have the following corollary. Corollary 4.2 Let Ay D r, then jy.x/j h W, where W is the number of non-zero components of r. Notice that for a discrete delta function, it should satisfy at least the zeroth moment equation [7], that is, X ıh .xi ˛/ D 1,
(43)
i
corresponding to the continuous case Copyright © 2013 John Wiley & Sons, Ltd.
R
ıh .x ˛/dx D 1. Now, we are ready to prove the main theorem. Math. Meth. Appl. Sci. 2013
Z. LI ET AL. Theorem 4.3 The IB method for (34) is first-order accurate, that is, N kEk1 Ch,
(44)
where CN is a constant. Proof We can decompose the local truncation error into two groups T D Treg C Tirreg ,
(45)
where kTreg k1 D 0 corresponds to the local truncation errors at regular grid points where ıh .xi ˛/ D 0 and the true solution is piecewise linear. Note that we have n1 X
Ti D
X
Treg C
X
Tirreg D O C
X
Tirreg .
(46)
iD1
On the other hand, we also have n1 X
Ti D
iD1
n1 X
cıhhat .xi ˛/
n1 X
iD1
cıh .xi ˛/ D 0,
(47)
iD1
because the finite difference method using the discrete delta function gives the exact solution at the grid points. Thus, we have P irreg irreg irreg,C D 0. We divide Ti into two groups, one with all positive Ti ’s denoted as Ti ; the other one is all the negatives i D Ti P P irreg, irreg,C irreg, irreg,C irreg, . Because Ti C Ti D 0, Ti and Ti must have the same order of the magnitude ( O.1=h/) denoted as Ti although those indexed i are different except that jxi ˛j Wh is true for all irregular grid points. Because the solution is linear with c, we have (48) E D A1 T D A1 Tirreg,C C Tirreg, . From the solution expression, we know that, assuming that xk ˛, 0 Ek D @
X
irreg xi Ti .1 xk / C
X
1 irreg xk Ti .1 xi /A
xi >xk
xi xk
0 1 X irreg,C X irreg, A C O.Wh/ T C T D ˛.1 ˛/ @ i
i
j
j
D O.Wh/, after we expand all xi ’s and xk ’s at ˛ and because all related xi and xj are within Wh distance from the interface ˛. The proof when xk > ˛ is similar except that we need to use the solution for x > ˛. This completes the proof.
5. Conclusions In this paper, we have derived the new jump condition when both single-layer and double-layer singularities are present along with a discontinuity in the coefficient of a general 1D interface problem. We discussed the effect of the boundary delta function singularities on the solution and derived the equivalent regular differential equations and equivalent boundary conditions. We also gave a convergence proof of the IB method in the L1 norm. We showed that the IB method is first-order accurate if a general discrete delta function is used.
Acknowledgements The first author was partially supported by the US ARO grant 550694-MA, the AFSOR grant FA9550-09-1-0520, the US NSF grant DMS0911434, and the NIH grant 096195-01, and CNSF 11071123. The second author was partially supported by the National Natural Science Foundation of China under grant 10971102 and the Natural Science Foundation of Jiangsu Province of China under grant BK2009398. The third to seventh authors were participants in a summer Research Experiences for Undergraduate Students (REU) at North Carolina State University supported by National Science Foundation (NSF), USA. Copyright © 2013 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2013
Z. LI ET AL.
References 1. Peskin CS, McQueen DM. A general method for the computer simulation of biological systems interacting with fluids. Symposia of the Society for Experimental Biology 1995; 49:265–276. 2. Li Z. The immersed interface method – a numerical approach for partial differential equations with interfaces. PhD thesis, University of Washington, 1994. 3. LeVeque RJ, Li Z. The immersed interface method for elliptic equations with discontinuous coefficients and singular sources. SIAM Journal of Numerical Analysis 1994; 31:1019–1044. 4. Li Z, Ito K. The Immersed Interface Method: Numerical Solutions of PDEs Involving Interfaces and Irregular Domains, SIAM Frontier Series in Applied Mathematics. FR33: Philadelphia, 2006. 5. Peskin CS. Lectures on mathematical aspects of physiology. Lectures in Applied Mathematics 1981; 19. 6. Evans LC. Partial Differential Equations. AMS: Rhode Island, 1998. 7. Beyer RP, LeVeque RJ. Analysis of a one-dimensional model for the immersed boundary method. SIAM Journal on Numerical Analysis 1992; 29:332–364. 8. Kevorkian J. Partial Differential Equations, Second Edition. Springer-Verlag: New York, 2010. 9. Mori Y. Convergence proof of the velocity field for a Stokes flow immersed boundary method. Communication Pure Applied Mathematics 2008; 61:1213–1263. 10. Chen K-Y, Feng K-A, Kim Y, Lai M-C. A note on pressure accuracy in immersed boundary method for stokes flow. Journal of Computational Physics 2011; 230:4377–4383. 11. Beale JT, Layton AT. On the accuracy of finite difference methods for elliptic problems with interfaces. Mathematics and Computer Science 2006; 1:91–119.
Copyright © 2013 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2013
Published online in Wiley Online Library
(wileyonlinelibrary.com) DOI: 10.1002/mma.2865 MOS subject classification: 65M12; 65M20
Some new analysis results for a class of interface problems Zhilin Lia,b * † , Li Wangc , Eric Aspinwalld , Racheal Coopere , Paul Kuberryf , Ashley Sandersg and Ke Zengh Communicated by Q. Wang Interface problems modeled by differential equations have many applications in mathematical biology, fluid mechanics, material sciences, and many other areas. Typically, interface problems are characterized by discontinuities in the coefficients and/or the Dirac delta function singularities in the source term. Because of these irregularities, solutions to the differential equations are not smooth or discontinuous. In this paper, some new results on the jump conditions of the solution across the interface are derived using the distribution theory and the theory of weak solutions. Some theoretical results on the boundary singularity in which the singular delta function is at the boundary are obtained. Finally, the proof of the convergency of the immersed boundary (IB) method is presented. The IB method is shown to be first-order convergent in L1 norm. Copyright © 2013 John Wiley & Sons, Ltd. Keywords: immersed boundary (IB) method; immersed interface method (IIM); jump conditions; discontinuous coefficient; Dirac delta function; weak solution; boundary singularity; convergence of IB method; equivalent boundary conditions
1. Introduction Interface problems arise in many areas of science and mathematics. Whenever there are different materials such as oil and water, gas and water, or bubbles in a fluid, we have interface problems. Another kind of interface problems is the Peskin’s immersed boundary (IB) method (see, e.g., [1] for a review and references therein). The IB method is not only a mathematical modeling tool in which a complicated boundary condition can be treated as a source distribution but also a numerical method in which a discrete delta function is used. The IB method is robust, simple, and has been applied to many problems in mathematical biology, fluid mechanics, material science, and other areas. However, the IB method is only first-order accurate in general. The immersed interface method (IIM) [2–4] is a second-order accurate method particularly designed for interface problems. In the IIM, jump conditions are enforced near the interface to obtain high-order accuracy. Knowing the jump conditions is crucial in the IIM. Mathematically, there are two characteristics in an interface problem. One is that a discontinuity exists in the coefficient of the differential equation along a co-dimensional space. The other one is singular source terms such as the Dirac delta function or its derivatives along a co-dimensional space. The following equation is a model interface problem in one dimension (1D), .ˇ.x/u0 /0 D f .x/ C c ı.x ˛/ C cN ı 0 .x ˛/,
0 < x < 1,
0 < ˛ < 1,
(1)
where ˇ.x/ D
ˇ .x/ ˇ C .x/
if x < ˛, if x > ˛.
(2)
a Department of Mathematics, North Carolina State University, Raleigh, NC 27695, USA b School of Mathematical Sciences, Nanjing Normal University, Nanjing, China c School of Mathematical Sciences and Jiangsu Key Laboratory for NSLSCS, Nanjing Normal University, Nanjing, China d Department of Mathematics, Florida State University, Tallahassee, FL 32306, USA e Department of Mathematical Sciences, Virginia Commonwealth University, Richmond, VA 23284, USA f Department of Mathematical Sciences, Clemson University, Clemson, SC 29634, USA g Department of Mathematics, Jackson State University, Jackson, MS 39217, USA h Department of Mathematics, Zhejiang University, Hangzhou, Zhejiang Province, 310058, China *Correspondence to: Zhilin Li, Department of Mathematics, North Carolina State University, Raleigh, NC 27695, USA. † E-mail: [email protected]
Copyright © 2013 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2013
Z. LI ET AL. along with two boundary conditions at x D 0 and x D 1. In the IB method developed by Peskin [1, 5], the Dirac delta function is intensively used. Often, a delta function singularity corresponds to a source distribution, whereas the derivative of the delta function singularity corresponds to a dipole distribution. For Stokes and Navier–Stokes equations, a delta function distribution along a surface indicates that the derivative of the velocity has jump discontinuities across the interface, whereas the derivative of the delta function indicates that the pressure has a jump discontinuity. There are three goals in this paper. The first one is to derive the jump conditions for (1). The second one is to consider the limiting case as ˛ approaches one of boundary points. The third one is to provide a convergence proof for the IB method. Because there are no other regularities except at the interface ˛, we know that the solution is piecewise smooth excluding the interface ˛, that is, u.x/ 2 C 2 .0, 1/ n ˛. The discontinuity is at the interface ˛. We can derive the jump condition of the solution and its derivative across the interface. We define ˇ D ˇ.˛/ D
lim
x!˛,x˛
ˇ.x/,
and Œˇ D ˇ C ˇ .
(3)
Jump conditions are crucial in the immersed interface method [2–4]. If cN D 0 or ˇ is continuous at ˛, that is, ˇ D ˇ C , then the jump conditions have been derived (e.g., [4]). We will derive the more challenging case when cN ¤ 0 and ˇ ¤ ˇ C . To the best of our knowledge, we have not seen the jump conditions and the derivation in the literature for the difficult case. We also report theoretical results and the derivation for the boundary singularity problem for which ˛ is at one of the boundary points. We will see that the effect of the delta function singularity can be absorbed by suitable boundary conditions. In this paper, we will use the notations u0 .x/ D ux .x/ D du dx .x/ whenever it is convenient. The rest of the paper is organized as follows. In the next section, we define the weak solution and derive the jump conditions. In Section 3, we discuss the case when the delta function singularity appears at one of the boundaries. We show that it can be translated as a contribution from the boundary conditions. In Section 4, we show the convergence of the IB method. Our result confirms that the IB method is first-order accurate in L1 in general. We conclude in the last section.
2. Deriving the jump conditions To use the IIM, we must know or derive jump conditions. The jump conditions describe the discontinuity of the solution and its derivative. The information is necessary to obtain second-order accuracy of the finite difference scheme as seen in [3]. We first state one of the main results of this paper in the following theorem. Theorem 2.1 Assume ˇ.x/ 2 C 1 Œ0, ˛/ [ C 1 .˛, 1, f .x/ 2 C.0, 1/. Then, the following jump conditions hold: Œu D
2Nc , ˇ C ˇC
Œˇux D c.
(4)
We prove the theorem by considering the following cases although the last case can cover all the situations. The purpose is to show that for simpler cases, the proofs are also simpler. 2.1. The singularity from ı.x ˛/ only In this simplest case, we assume that Œˇ D 0, cN D 0. Thus, the ODE is simply .ˇux /x D f .x/ C cı.x ˛/. From the differential equation theory (e.g., [6]), we know that the solution is continuous, that is, Œu D 0. Integrating the differential equation from ˛ to ˛ C , we have Z
Z
˛C ˛
.ˇux /x dx D
Z
˛C
˛
f .x/dx C
˛C
˛
cı.x ˛/dx.
(5)
Notice that f .x/ is continuous in the neighborhood of ˛, and consequently, the first integration in the right-hand side vanishes. Also, from the definition of the Dirac delta function, we have the following property of the delta function: Z
C
g.x/ı.x ˛/dx D g.˛/
(6)
for any continuous function g.x/. Thus, from (5), we have ˇ C uC x .˛/ ˇ ux .˛/ D c,
or Œˇux D c.
Thus, the theorem is true for this case. 2.2. Singularities from ı.x ˛/ and ı 0 .x ˛/ In this case, we assume that Œˇ D 0 and both c and cN can be non-zero. Thus, the ODE is .ˇux /x D f .x/Cc ı.x ˛/C cN ı 0 .x ˛/. The presence of the derivative of the delta function corresponds to a double layer in the potential theory. Thus, theRsolution typically has a finite jump across the interface ˛. We cannot simply integrate the differential equation from ˛ to ˛ C because cN ı 0 .x ˛/dx does not exist. Copyright © 2013 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2013
Z. LI ET AL. Similar to the process in deriving the weak form in the finite element method, we multiply a testing function .x/ 2 C01 .0, 1/ to both sides of the differential equation and then integrate to obtain Z
Z
˛C ˛
.ˇux /x .x/dx D
Z
˛C
˛
f .x/.x/dx C
Z
˛C
˛
cı.x ˛/.x/dx C
˛C
˛
cN ı 0 .x ˛/.x/dx
(7)
D c.˛/ cN 0 .˛/. The right-hand side is obtained from the definition of the weak derivatives of the delta function. We manipulate the left-hand side of (7) as follows: ˇ˛ C Z ˇ LHS D .ˇux / ˇ ˛
˛C
˛
D Œˇux .˛/ ˇ.˛/
ˇux 0 .x/dx ˇ˛ C Z ˇ .u /ˇ 0
˛
˛C
˛
!
(8)
00
u .x/dx
D Œˇux .˛/ ˇ.˛/ Œu 0 .˛/. Because .x/ is arbitrary and Equation (7) is the same as (8), we have to have ˇ.˛/ Œu D cN ,
Œˇux D c,
that is,
Œu D
cN 2Nc . D ˇ.˛/ ˇ C ˇC
Thus, the theorem is also true for this case. 2.3. Singularities from ı.x ˛/, ı 0 .x ˛/, and a discontinuity in ˇ Now the coefficient ˇ.x/ has a discontinuity at ˛. For simplicity of proof, we assume that ˇ is piecewise constant if x < ˛, ˇ ˇ.x/ D ˇ C if x > ˛.
(9)
Before we derive the jump conditions, we define a weak solution using the standard regularization technique. We will see that the solution defined from the jump conditions is the same as the weak solution. First, we define a regularized Heaviside function H .x/ as follows: 8 < 0 if x < , 1 (10) H .x/ D 1 C x C 1 sin x if jxj , 2 : 1 if x > . Note that H .0/ D 12 . We define the regularized coefficient as ˇ .x/ D ˇ C Œˇ H .x ˛/. Note that ˇ .˛/ D
1 2
ˇ C ˇ
C
(11)
. We define the regularized solution u .x/ as the solution of the following boundary value problem:
.ˇ .x/u0 /0 D f .x/ C c ı.x ˛/ C cN ı 0 .x ˛/,
0 < x < 1,
0 < ˛ < 1.
(12)
Because ˇ .x/ is continuous now, we have Œu D
2Nc , ˇ .˛C/ C ˇ .˛/
ˇ u0 D c.
(13)
Note that we have lim Œu D lim
!0
!0
2Nc 2Nc D C . ˇ .˛C/ C ˇ .˛/ ˇ C ˇ
(14)
Definition 2.2 The weak solution to (1) is defined as uw .x/ D lim u .x/. !0
(15)
It is obvious that the weak solution satisfies the jump conditions (4). We now prove that from the original differential equations, we can derive the same jump conditions as that of the weak solution. Copyright © 2013 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2013
Z. LI ET AL. 2.3.1. The derivation of the general jump conditions. As before, we multiply a testing function .x/ 2 C01 .0, 1/ to both sides of the differential equation and then integrate to obtain Z
Z
˛C ˛
.ˇux /x .x/dx D
Z
˛C
f .x/.x/dx C
˛
Z
˛C
cı.x ˛/.x/dx C
˛
˛C
˛
cN ı 0 .x ˛/.x/dx
(16)
0
D c.˛/ cN .˛/. The right-hand side is obtained from the definition of the weak derivatives. We manipulate the left-hand side of (16) as follows: ˇ˛ C Z ˇ LHS D .ˇux / ˇ ˛
Z D Œˇux .˛/
˛C ˛
˛C
˛
ˇ.x/ux 0 .x/dx
ˇ.x/ux 0 .x/dx.
We cannot take ˇ.x/ out of the integration anymore because it is discontinuous. We represent ˇ.x/ and u.x/ as follows (see also [7]): ˇ.x/ D ˇ C Œˇ H.x ˛/, 1 u.x/ D u C Œu H.x ˛/ C Œux H.x ˛/.x ˛/ C Œuxx H.x ˛/.x ˛/2 C . 2 Differentiating the aforementioned representation, we obtain ux .x/ D Œu ı.x ˛/ C Œux ı.x ˛/.x ˛/ C Œux H.x ˛/ C . Plugging the aforementioned expressions into the integral of LHS, we have Z
˛C ˛
ˇ.x/ux 0 .x/dx D
D
Z
˛C
Z
˛C
ˇ C Œˇ H.x ˛/ Œu ı.x ˛/ C Œux ı.x ˛/.x ˛/ ˛ C Œux H.x ˛/ C . 0 .x/dx
˛
C
ˇ Œuı.x ˛/ 0 .x/dx C
Z
˛C ˛
0
Z
˛C
˛
ŒˇŒuH.x ˛/ı.x ˛/ 0 .x/dx
ˇ Œux H.x ˛/ .x/dx C Z
Z
˛C
˛
ŒˇŒux H2 .x ˛/ 0 .x/dx C
d ŒˇŒuH2 .x ˛/ 0 .x/dx C 0 C 0 C ˛ dx Z ˛C 1 D ˇ Œu 0 .˛/ C ŒˇŒuH2 .x ˛/ 0 .˛/ ŒˇŒuH2 .x ˛/ 00 .x/dx 2 ˛ 1 D ˇ Œu 0 .˛/ C ŒˇŒu 0 .˛/ C 0 2 ˇ C ˇC Œu 0 .˛/. D 2 1 D ˇ Œu .˛/ C 2
0
˛C
Because .x/ 2 C01 is arbitrary, by matching the coefficients in front of .˛/ and 0 .˛/ in (16) and using the aforementioned relation, we have the jump conditions (4).
3. The boundary singularity analysis In this section, we discuss what happens when a delta function appears at one of the boundary points. We call such an irregularity a boundary singularity for short. This problem is worth discussing for several reasons. First of all, in some real applications, the interface is quite close to the boundary. Second, for many applications, sources may indeed appear on parts of boundaries. Third, it is theoretically significant to change the irregular problem with a ı function on the right-hand side to a regular partial differential equation in order to make it more convenient to solve. The discussions here have some similarities of the limiting value theory (e.g., [8]). Our discussions in this section indicate that we can translate the delta function singularity on a boundary into a corresponding boundary condition. Without loss of generality, we consider that the boundary singularity is at the boundary x D 0. In this case, the discontinuity in the coefficient is irrelevant. Note that it is rarely seen a dipole (derivative of the delta function) singularities at boundaries. From another point of view, the two pieces of the solutions are almost arbitrarily independent because of arbitrary parameters c and cN for the dipole case. Therefore, we only consider the case when a delta function singularity appears at one of boundary points. Thus, the problem can be written as .ˇ.x/u0 /0 D f .x/ C c ı.x/, Copyright © 2013 John Wiley & Sons, Ltd.
0 < x < 1.
(17) Math. Meth. Appl. Sci. 2013
Z. LI ET AL. For simplicity, we assume a Dirichlet boundary condition at x D 1 with u.1/ D and a linear boundary condition at x D 0, that is, a Dirichlet, or Neumann, or a mixed boundary condition. Note that the solution can be decomposed as the sum of two problems: 0 (18) ˇ.x/u01 D f .x/ 0 < x < 1, with homogenous but the same type of boundary conditions, and 0 ˇ.x/u02 D c ı.x/,
0 < x < 1,
(19)
with the original boundary conditions. Because u1 .x/ is the solution to a standard problem, we can treat it easily. Thus, we will simply discuss the second equation. To further simplify the discussion, we assume that ˇ D 1. The problem is simplified to u00 D c ı.x/,
0 < x < 1.
(20)
We first define the weak solution of (20). Definition 3.1 We consider u00˛ D c ı.x ˛/,
0 < x < 1,
0 < ˛ < 1,
(21)
with a boundary condition at x D 0 and x D 1. The weak solution of (20) is defined as u.x/ D lim˛!0 u˛ .x/. To consider the effect of the delta function singularity at one boundary, we use the following procedure: use the jump condition and the boundary condition to obtain the analytical solution for u˛ .x/; take the limit of the solution on the right interval .˛, 1 as ˛ approaches zero; and use the limiting solution to derive the equivalent equations and equivalent boundary condition.
3.1. Dirichlet boundary conditions In this case, the interface problem has the following form: uxx .x/ D cı.x ˛/, u.0/ D
0 < ˛ < 1,
u.1/ D
(22)
0 < ˛ < 1.
Theorem 3.2 The aforementioned limiting solution lim u˛ .x/ is the solution of the following problem: ˛!0
uxx D 0,
u.0/ D ,
u.1/ D .
(23)
That is, for a Dirichlet boundary condition, the boundary singularity has no effect on the solution. Proof We know that the solution has the following form: u˛ .x/ D
v1 C v2 x, v3 C v4 x,
0 < x < ˛, ˛ < x < 1.
(24)
Using the known jump conditions Œu D 0 and Œux D c and the boundary condition u.0/ D , we have the following: v1 D v1 C v2 ˛ D v3 C v4 ˛ v2 C c D v4 v3 C v4 D . Solving the linear system, we have v1 D v2 D c˛ C c v3 D c˛ v4 D c˛ C . Consequently, the analytical solution is as follows: C .c˛ C c/x, u˛ .x/ D c˛ C .c˛ C /x, Copyright © 2013 John Wiley & Sons, Ltd.
0 < x < ˛, ˛ < x < 1.
(25)
Math. Meth. Appl. Sci. 2013
Z. LI ET AL. When ˛ approaches zero, considering the analytical result on the right interval .˛, 1, the solution is as follows: uN .x/ D C . /x.
(26)
It is the solution of the equivalent regular problem: uxx .x/ D 0 u.0/ D
(27)
u.1/ D . This completes the proof.
3.2. Neumann boundary conditions We consider the problem with a Neumann boundary condition at x D 0. uxx .x/ D cı.x ˛/, ux .0/ D , u.1/ D
0 < ˛ < 1,
(28)
0 < ˛ < 1.
Theorem 3.3 The aforementioned limiting solution lim u˛ .x/ is the solution of the following problem: ˛!0
uxx D 0,
ux .0/ D C c,
u.1/ D .
(29)
In other words, the boundary source strength was added to the derivative boundary condition. Proof We know that the solution has the following form: u˛ .x/ D
v1 C v2 x, v3 C v4 x,
0 < x < ˛, ˛ < x < 1.
Using the known jump conditions Œu D 0 and Œux D c and the boundary condition ux .0/ D , we have the following: v2 D v1 C v2 ˛ D v3 C v4 ˛ v2 C c D v4 v3 C v4 D . Solving the linear system, we have v1 D C c˛ c v2 D v3 D c v4 D C c. Consequently, the analytical solution is as follows: C c˛ c C x, u˛ .x/ D c C . C c/x,
0 < x < ˛, ˛ < x < 1.
When ˛ approaches zero, considering the right interval .˛, 1, the approximate solution is as follows: uN .x/ D c C . C c/x. This is the weak solution of the equivalent equation for the Neumann boundary condition without the ı function: uxx .x/ D 0, ux .0/ D C c,
(30)
u.1/ D . That completes the proof. Copyright © 2013 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2013
Z. LI ET AL. 3.3. Mixed boundary conditions We now consider a mixed (Robin) boundary condition at x D 0 of the following problem: uxx .x/ D cı.x ˛/, u.0/ C ˇux .0/ D ,
(31)
u.1/ D 0 < ˛ < 1. Theorem 3.4 The aforementioned limiting solution lim u˛ .x/ is the solution of the following problem: ˛!0
uxx D 0,
u.0/ C ˇux .0/ D cˇ C ,
u.1/ D .
(32)
Again, we see that the source of the boundary singularity was added to the derivative part of the mixed boundary condition. Proof We know that the solution has the following form: u˛ .x/ D
v1 C v2 x, v3 C v4 x,
0 < x < ˛, ˛ < x < 1.
Using the known jump conditions Œu D 0 and Œux D c and the boundary condition u.0/ C ˇux .0/ D , we have the following: v1 C ˇv2 D v1 C v2 ˛ D v3 C v4 ˛ v2 C c D v4 v3 C v4 D . Solving the linear system, we have v1 D
ˇ.Cc˛c/ ˇ
v2 D
C.cc˛/ ˇ
v3 D
ˇ.c/Cc˛ ˇ
v4 D
Ccˇ .Cc˛/ . ˇ
Therefore, the analytical solution is as follows: 8 < ˇ.Cc˛c/ C C.cc˛/ x ˇ ˇ u˛ .x/ D : ˇ.c/Cc˛ C Ccˇ .Cc˛/ x ˇ
0 < x < ˛, ˛ < x < 1.
ˇ
When ˛ approaches zero, using the same technique, the limiting solution is as follows: uN .x/ D
ˇ. c/ C cˇ C x. ˇ ˇ
(33)
Because uN .0/Cˇ uNx .0/ D cˇ C , the equivalent regular differential equation for the mixed boundary condition without the ı function is the following: uxx .x/ D 0 u.0/ C ˇux .0/ D cˇ C ,
u.1/ D .
This completes the proof.
4. A proof of the convergence of the immersed boundary method It is well known that Peskin’s IB method is only first-order accurate in L1 norm. In [9], the author has proved the first-order accuracy of the IB method for the Stokes equations (see also [10]). However, there are few theoretical proofs on the IB method for elliptic interface problems. This is the main motivation of this section. We will give a proof for the 1D model, u00 D c ı.x ˛/ 0 < x < 1, Copyright © 2013 John Wiley & Sons, Ltd.
0 < ˛ < 1,
u.0/ D u.1/ D 0,
(34) Math. Meth. Appl. Sci. 2013
Z. LI ET AL. in this section. The proof for the 2D problem is more challenging and is under investigation. Note that the analytic solution to the aforementioned equation cx .1 ˛/ if x xk , (35) u.x/ D c˛ .1 x/ otherwise. Given a uniform Cartesian grid xi D ih, i D 0, 1, : : : , n, h D 1=n, the IB method leads to the following system of linear equations: Ui1 2Ui C UiC1 D cıh .xi ˛/ , h2
i D 1, 2, : : : , n 1,
(36)
where Ui is the finite difference approximation of the solution u.xi /, and ıh .xi ˛/ is a discrete delta function applied to the grid point xi . In the matrix-vector form, the aforementioned finite difference equations can be written as AU D F, where A is the tridiagonal matrix formed by the discrete Laplacian. It is well known that A is a symmetric positive definite matrix and diagonally dominant. Note that a discrete delta function has a compact support in the neighborhood of the interface, that is, ıh .x/ ¤ 0, only if jxj Wh, where W is a constant. Commonly used discrete delta functions include the hat discrete delta function (ı hat .x/ with W D 1) and the cosine discrete delta function (ı cosine .x/ with W D 2) (e.g., [4]). Note that when we use the hat delta function, the result is the same as that of the IIM for the simple model. The solution to the finite difference equations is the same as a true solution if there are no round-off errors, that is, u.xi1 / 2u.xi / C u.xiC1 / D cıhhat .xi ˛/ , h2
i D 1, 2, : : : , n 1
(37)
(see, e.g., [4, 11]). But this is not the case for other discrete delta functions. We define the error vector as E D fEi g, where Ei D u.xi / Ui . The local truncation error is defined as T D fTi g, Ti D
u.xi1 / 2u.xi / C u.xiC1 / cıh .xi ˛/ . h2
(38)
With the definition, we have AU D F and Au D F C T, and therefore, AE D T. For the hat discrete delta function, we have jTi j D 0 for all i’s. For other discrete delta functions, we have jTj j D O.1=h/ for a few grid points neighboring the interface ˛. So, the interesting question is as follows: Why is the IB method still first-order accurate, that is, kEk1 D O.h/? To answer this question, we first introduce the following lemma. Lemma 4.1 Let Ay D ek , where ek is the kth unit base vector, then yi D
hxi .1 xk / hxk .1 xi /
if i k, otherwise.
(39)
The significance of this lemma is that the solution is order h smaller than the concentrated source. Proof We note the following identity: Ay D h
ek D h ıhhat ek . h
(40)
For this simple case, the IB and IIM are identical. Thus, from the IIM [4, 11], we know that y is the exact discrete solution at the grid points of the following boundary value problem u00 D h ı.x xk /,
0 < x < 1,
u.0/ D u.1/ D 0,
(41)
whose solution is y.x/ D
hx .1 xk / hxk .1 x/
if x xk , otherwise.
This completes the proof.
(42)
Note that jy.x/j h. From this lemma, we have the following corollary. Corollary 4.2 Let Ay D r, then jy.x/j h W, where W is the number of non-zero components of r. Notice that for a discrete delta function, it should satisfy at least the zeroth moment equation [7], that is, X ıh .xi ˛/ D 1,
(43)
i
corresponding to the continuous case Copyright © 2013 John Wiley & Sons, Ltd.
R
ıh .x ˛/dx D 1. Now, we are ready to prove the main theorem. Math. Meth. Appl. Sci. 2013
Z. LI ET AL. Theorem 4.3 The IB method for (34) is first-order accurate, that is, N kEk1 Ch,
(44)
where CN is a constant. Proof We can decompose the local truncation error into two groups T D Treg C Tirreg ,
(45)
where kTreg k1 D 0 corresponds to the local truncation errors at regular grid points where ıh .xi ˛/ D 0 and the true solution is piecewise linear. Note that we have n1 X
Ti D
X
Treg C
X
Tirreg D O C
X
Tirreg .
(46)
iD1
On the other hand, we also have n1 X
Ti D
iD1
n1 X
cıhhat .xi ˛/
n1 X
iD1
cıh .xi ˛/ D 0,
(47)
iD1
because the finite difference method using the discrete delta function gives the exact solution at the grid points. Thus, we have P irreg irreg irreg,C D 0. We divide Ti into two groups, one with all positive Ti ’s denoted as Ti ; the other one is all the negatives i D Ti P P irreg, irreg,C irreg, irreg,C irreg, . Because Ti C Ti D 0, Ti and Ti must have the same order of the magnitude ( O.1=h/) denoted as Ti although those indexed i are different except that jxi ˛j Wh is true for all irregular grid points. Because the solution is linear with c, we have (48) E D A1 T D A1 Tirreg,C C Tirreg, . From the solution expression, we know that, assuming that xk ˛, 0 Ek D @
X
irreg xi Ti .1 xk / C
X
1 irreg xk Ti .1 xi /A
xi >xk
xi xk
0 1 X irreg,C X irreg, A C O.Wh/ T C T D ˛.1 ˛/ @ i
i
j
j
D O.Wh/, after we expand all xi ’s and xk ’s at ˛ and because all related xi and xj are within Wh distance from the interface ˛. The proof when xk > ˛ is similar except that we need to use the solution for x > ˛. This completes the proof.
5. Conclusions In this paper, we have derived the new jump condition when both single-layer and double-layer singularities are present along with a discontinuity in the coefficient of a general 1D interface problem. We discussed the effect of the boundary delta function singularities on the solution and derived the equivalent regular differential equations and equivalent boundary conditions. We also gave a convergence proof of the IB method in the L1 norm. We showed that the IB method is first-order accurate if a general discrete delta function is used.
Acknowledgements The first author was partially supported by the US ARO grant 550694-MA, the AFSOR grant FA9550-09-1-0520, the US NSF grant DMS0911434, and the NIH grant 096195-01, and CNSF 11071123. The second author was partially supported by the National Natural Science Foundation of China under grant 10971102 and the Natural Science Foundation of Jiangsu Province of China under grant BK2009398. The third to seventh authors were participants in a summer Research Experiences for Undergraduate Students (REU) at North Carolina State University supported by National Science Foundation (NSF), USA. Copyright © 2013 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2013
Z. LI ET AL.
References 1. Peskin CS, McQueen DM. A general method for the computer simulation of biological systems interacting with fluids. Symposia of the Society for Experimental Biology 1995; 49:265–276. 2. Li Z. The immersed interface method – a numerical approach for partial differential equations with interfaces. PhD thesis, University of Washington, 1994. 3. LeVeque RJ, Li Z. The immersed interface method for elliptic equations with discontinuous coefficients and singular sources. SIAM Journal of Numerical Analysis 1994; 31:1019–1044. 4. Li Z, Ito K. The Immersed Interface Method: Numerical Solutions of PDEs Involving Interfaces and Irregular Domains, SIAM Frontier Series in Applied Mathematics. FR33: Philadelphia, 2006. 5. Peskin CS. Lectures on mathematical aspects of physiology. Lectures in Applied Mathematics 1981; 19. 6. Evans LC. Partial Differential Equations. AMS: Rhode Island, 1998. 7. Beyer RP, LeVeque RJ. Analysis of a one-dimensional model for the immersed boundary method. SIAM Journal on Numerical Analysis 1992; 29:332–364. 8. Kevorkian J. Partial Differential Equations, Second Edition. Springer-Verlag: New York, 2010. 9. Mori Y. Convergence proof of the velocity field for a Stokes flow immersed boundary method. Communication Pure Applied Mathematics 2008; 61:1213–1263. 10. Chen K-Y, Feng K-A, Kim Y, Lai M-C. A note on pressure accuracy in immersed boundary method for stokes flow. Journal of Computational Physics 2011; 230:4377–4383. 11. Beale JT, Layton AT. On the accuracy of finite difference methods for elliptic problems with interfaces. Mathematics and Computer Science 2006; 1:91–119.
Copyright © 2013 John Wiley & Sons, Ltd.
Math. Meth. Appl. Sci. 2013
