An Augmented IB Method & Analysis for Elliptic BVP on Irregular Domains

: The immersed boundary method is well-known, popular, and has had vast areas of applications due to its simplicity and robustness even though it is only ﬁrst order accurate near the interface. In this paper, an immersed boundary-augmented method has been developed for linear elliptic boundary value problems on arbitrary domains (exterior or interior) with a Dirichlet boundary condition. The new method inherits the simplicity, robustness, and ﬁrst order convergence of the IB method but also provides asymptotic ﬁrst order convergence of partial derivatives. Numerical examples are provided to conﬁrm the analysis.


Introduction
In this paper, we develop an augmented immersed boundary method for linear elliptic boundary value problems (BVP) on arbitrary domains. Since its introduction in 1970's, the Immersed Boundary (IB) method [Peskin (1977)] has been flourished in many fields including mathematics, engineering, biology, computational fluid dynamics (CFD), see for example Peskin et al. [Peskin and McQueen (1995)] for reviews and references therein. The IB method is not only a mathematical modeling tool in which complicated boundary conditions can be treated as source distributions but also a numerical method in which a discrete delta function is utilized to spread the singular source to nearby grid points. The IB method is robust, simple, and has been applied to many problems. The problem of interest in this paper is to solve the following boundary value problem numerically x ∈ Ω, u| ∂Ω = u 0 (X).
It is challenging to solve the problem above both in discretization and solving the resulting linear system of equations after the discretization. Various efforts can be found in the literature including the capacitance matrix method, fictitious domain methods, boundary integral methods if the coefficients are constants, and the augmented immersed interface method [Hunter, Li and Zhao (2002); Li, Zhao and Gao (1999)] etc., and some immersed boundary (IB) based methods. Most of the methods mentioned are relatively sophisticated. The purpose of this paper is to develop the IB method for the irregular domain problem that has the same nature as the original IB method in terms of the accuracy and the simplicity. The idea is to extend the setting to a rectangular domain so that the Peskin's IB method can be applied. The source strength involving a two-dimensional Dirac delta function should be such a one that the boundary condition is satisfied, which once again fits the Peskin's interpolation formula that involves the Dirac delta function. We also propose a simple method to approximate the first order partial derivatives using the computed finite difference solution.
The rest of paper is organized as follows. In the next section, we explain the augmented idea that transforms the irregular domain problem to an interface problem with an unknown source strength v(∂Ω). The solution of the interface problem should satisfy the boundary condition using the distribution theory. The discrete version is also presented. In Section 3, we present the convergence analysis that confirms the first order convergence of solution.
In Section 4, we describe an algorithm to approximate the first order partial derivatives followed by the error analysis. In Section 5, we present numerical examples of both interior and exterior problems with a general boundary. We conclude and acknowledge in the last section.

An augmented IB method for linear elliptic BVP
We use an interior problem to explain the idea as illustrated in Fig. 1. We first embed the domain Ω into a rectangular domain R = [a, b] × [c, d]. We also assume that ∂Ω ∩ ∂R is empty, see Fig. 1 for an illustration. In discretization, the discrete boundaries ∂Ω h and ∂R h are at least several grids apart (≥ 5h). The idea is to extend the boundary value problem to a rectangular domain R; then re-write the problem using the Peskin's immersed boundary formulation, where x = (x, y) is a point in the domain R ⊃ Ω , X = (X(s), Y (s)) is a point on ∂Ω assuming that ∂Ω has a parametric representation ∂Ω = (X(s), Y (s)). The modified source term F (x) has the following form, The coefficientsβ,σ are non-negative extensions of β, σ from Ω to R. In this way, we have transformed the original irregular domain problem to an interface problem with an additional boundary variable v(s) which has been shown to be the flux jump across the boundary ∂Ω, The extended problem is still wellposed. The restriction of the solution of u(x, y) on Ω is the solution to the original problem. Note that, when β and σ are constants, we use the same β and σ in the entire rectangular domain R.

Augmented IB discretization
There are variety of methods that can solve the above interface problem accurately, for example, the augmented immersed interface method. However, in this paper we wish to investigate how the IB method can be applied to solve the problem since it is straightforward, simple, and robust.
For simplification of discussion, we use a uniform mesh assuming R = (a, b) × (c, d), where M and N , M ∼ N , are the number of the grid lines in x and y directions, respectively. The boundary ∂Ω is represented by a set of ordered control points (X k , Y k ), k = 1, 2, · · · N b , which is then interpolated using a periodic cubic spline. We use U ij to represent the finite difference approximation to u( and so on. The immersed boundary discretization then is simply the followinĝ whereβ i+1/2,j =β(x i + h/2, y j ),σ i,j =σ(x i , y j ) and so on, and Note that, there are sophisticated discrete Dirac delta functions in the literature such as the discrete radial delta function, for example, for various purposes. Here we simply apply the original Peskin's IB method for a first order method. The matrix-vector form of the discretization (8), (9) can be written as

Solving the linear system of equations
If we can solve the system of Eq. (12), then we get an approximate solution to the original boundary value problems at the grid points, inside for interior problems; and outside for exterior problems. Another approach is to solve V first through the Schur complement of (12) if we apply the block Gaussian elimination method to get where A −1 h F 1 is the result of the elliptic solver on the rectangular domain given the source term F 1 . It has been shown in Li et al. [Li and Ito (2006)] and other related papers, that the matrix-vector multiplication SV given V is simply where R(V) = SV −F is the residual of the boundary condition given V. More precisely, the k-th component of SV is where U ij 's depend on V through (10). The right hand side of SV =F can be computed from −R(0) which corresponds to the residual of the boundary condition with zero component values of the augmented variable, which is the result from the boundary value problem on the rectangular domain. Since we know how to get the matrix-vector multiplication, we can use a direct or GMRES iterative method to solve the system (13).

Convergence analysis
In this section we show that like the original IB method, the augmented IB method is also first order accurate except for a factor of log h. We use the exterior problem to illustrate the proof. The proof follows the path of the proof in Li et al. [Li, Zhao and Gao (1999)] for different problems. We use U and u to represent the vectors of approximate and exact solution at grid points on entire R; T u and E u = U − u are the vectors of the local truncation errors and the global error. We have E u | ∂Rh = 0 for the values at grid points on the boundary ∂R. Similarly, we define T q and E q = V − v as the vectors of the local truncation error and the global error for the augmented variable in (9). According to the definitions, we have where the local truncation error vector T u is defined as T u = F 1 − A h u − Bv and so on.
Theorem 3.1. Assume that u(x) ∈ C 2 (Ω), then the computed solution of the finite difference equations (8)- (9) is also first order accurate except for a factor of log h, that is, E ∞ ≤ Ch log h.
Proof: From (16)-(17), we have From the consistency of the interpolation formula (9), we know that T q ∞ ≤ Ch. From the results in Li et al. [Li (2015)], we also know that A −1 h T u ∞ ≤ Ch log h. From the consistency of the interpolation scheme (9), we also conclude that This is because the IB method has been shown to be first order accurate except for a factor of log h. Thus, U(E v ) corresponds to the boundary value problem with u| ∂R = 0 and with u| ∂Ω = O(h log h). From the maximum principle of elliptic partial differential equations, we know that the solution is also bounded by O(h log h). Finally, the Schur complement matrix S is an interpolation scheme from U ij to U (X k , Y k ) on the boundary. From the first order convergence and the consistency condition we have This completed the proof.

Approximation of first order partial derivatives
We also propose a simple method to find approximations of the first order derivatives at the control points on the boundary ∂Ω using a one-sided interpolation. We think that with suitable choices of the parameters, the accuracy of the computed first order derivative can also be first order accurate except for a factor of log h. Note that the IB method is a smoothing method in which the discontinuity in the first order derivatives will be smoothed out. However, away from the boundary, the solution is smooth and the derivatives should be close to the true ones. We use the interior boundary value problem to explain how to compute the first order derivatives of the solution at a boundary point (X k , Y k ) to demonstrate the procedure.
Let X k = (X k , Y k ) be a selected point on the boundary ∂Ω. We select a region Ω k : δ 1,h ≤ |X k − x ik,jk | ≤ δ 2,h , where x ik,jk ∈ Ω are regular grid points used (outside of the support of the discrete delta function) in the interpolation, for example, δ 1,h = 3h, δ 1,h = 4h or larger so that the region contains at least one grid point in the domain. The selected region should be outside of the domain where the solution will be smoothed out by the discrete delta function. We then use the standard central finite difference scheme to compute the first order derivatives In our approach, we use the average of computed first order derivatives in the selected region, that is, where N δ is the number of grid points in Ω k , see Fig. 2 for an illustration. Figure 2: An illustration to compute an approximation of first order derivatives at (X k , Y k ) ∈ ∂Ω using a regular grid point We conjecture that the computed first order derivatives are asymptotically first order convergent to the true ones at X k with suitable choices of δ i . Below is the sketch of the proof. We know that the computed solution is continuous, and it is smooth away from the boundary, which means where ∂Ei k ,j k ∂x ∼ O(1) and ∂Ei k+1 ,j k ∂x − ∂Ei k−1 ,j k ∂x ∼ O(h) since the error is also smooth away from the boundary. Thus we have where |T Dx ik,jk | ≤ Ch 2 is the truncation error from the finite difference scheme. Thus we conclude that the error of the first order derivative is of O(h log h). The behavior is asymptotic since we need relatively fine meshes to make sure that the interpolation points are not far from boundary but are away from the smoothed region.
We present a couple of numerical experiments to show the performance of the new augmented IB method for solutions and their normal derivative at the boundary. All the experiments are computed with the double precision and are performed on a MacBook Pro with Pentium(R) Dual-Core CPU, 2.6 GHz, 8 GB memory. We present errors in the infinity norm (pointwise). In all tables listed in this section, we use M = N = N b . We use a constant β and σ so that a fast Poisson solver based on Fast Fourier Transform (FFT) can be utilized. We use the GMRES method to solve the Schur complement system for the augmented variable with the tolerance 10 −6 and 0 as the initial guess in all computations. The interface is a general curve r = 0.5 + sin(Lθ), 0 < < 0.5, L is an integer, We start with an interior problem with r = x 2 + y 2 ≤ 0.5 + 0.2 sin(5θ). The analytic solution is designed as u(x, y) = sin x cos y.
The source term and the Dirichlet boundary condition are determined from the PDE and the analytic solution. In Tab. 1, we present a grid refinement analysis when β = 1, σ = 1.5. The second column in the table shows the maximum errors of the finite difference solution while the third column shows the ratios of E N (u)/E 2N (u) which are close to number two for a first order accurate method as we double the N . The first order convergence of solution is clearly verified. The fourth column shows the errors of the normal derivative of the solution at the boundary ∂Ω pointing outward, assuming the normal derivative at (X k , Y k ) is n = (cos θ, sin θ). The fifth column shows the ratios of consecutive errors which approach number two as the N gets larger. Next we show the numerical results for an exterior problem with the same boundary. The analytic solution is designed as u(x, y) = r 4 + C 0 log(2r).
The source term and the boundary condition are determined from the PDE and the analytic solution. In Tab. 2, we present a grid refinement analysis when β = 2.5, σ = 1.5. We observe similar behavior as the previous example.  As last example, we show the case β = 1, σ = 0, the interior boundary is simply a circle r = 0.5, that is, we have a Poisson equation outside the circle. The solution is u(x, y) = sin x cos y. The grid refinement analysis results are presented in Tab. 3. We can see that while the errors for the solution are smaller than that of the normal derivatives, the convergence order is roughly the same. The last column shows the number of GMRES iterations without any preconditioning techniques. In this example, the curvature is a constant and the errors behave more uniformly. In this paper, we proposed a new augmented immersed boundary method for general elliptic boundary value problems with a Dirichlet boundary condition on irregular domains. The new method inherits many properties of the original IB method for interface problems such as simplicity, robustness, and first order convergence. We also proposed a simple method to estimate first order partial derivatives which is asymptotically first order accurate. The convergence of method for the solution and the first order derivatives are also proved under appropriate regularity assumptions. How to generalize the method to other boundary conditions and how to develop efficient preconditioning techniques for the Schur complement using only the matrix and vector multiplications are two open challenges.