Immersed Boundary Smooth Extension: A highorder method for solving PDE on arbitrary smooth domains using Fourier spectral methods
Abstract
The Immersed Boundary method is a simple, efficient, and robust numerical scheme for solving PDE in general domains, yet it only achieves firstorder spatial accuracy near embedded boundaries. In this paper, we introduce a new highorder numerical method which we call the Immersed Boundary Smooth Extension (IBSE) method. The IBSE method achieves highorder accuracy by smoothly extending the unknown solution of the PDE from a given smooth domain to a larger computational domain, enabling the use of simple Cartesiangrid discretizations (e.g. Fourier spectral methods). The method preserves much of the flexibility and robustness of the original IB method. In particular, it requires minimal geometric information to describe the boundary and relies only on convolution with regularized deltafunctions to communicate information between the computational grid and the boundary. We present a fast algorithm for solving elliptic equations, which forms the basis for simple, highorder implicittime methods for parabolic PDE and implicitexplicit methods for related nonlinear PDE. We apply the IBSE method to solve the Poisson, heat, Burgers’, and FitzhughNagumo equations, and demonstrate fourthorder pointwise convergence for Dirichlet problems and thirdorder pointwise convergence for Neumann problems.
keywords:
Embedded boundary, Immersed Boundary, Fourier spectral method, Complex geometry, Partial Differential Equations, HighorderremarkRemark
1 Introduction
The Immersed Boundary (IB) method was originally developed for the study of moving, deformable structures immersed in a fluid, and it has been widely applied to such problems since its introduction [1, 2, 3]. Recently, the method has been adapted to more general fluidstructure problems, including the motion of rigid bodies immersed in a fluid [4] and fluid flow through a domain with either stationary boundaries or boundaries with prescribed motion [5, 6]. In this broadened context, we use the term Immersed Boundary method to refer only to methods in which (i) the boundary is treated as a Lagrangian structure embedded in a geometrically simple domain, (ii) the background PDE (e.g. the NavierStokes equations) are solved on a Cartesian grid everywhere in that domain, and (iii) all communication between the Lagrangian structure and the underlying PDE is mediated only by convolutions with regularized functions. These methods have many desirable properties: they make use of robust and efficient Cartesiangrid methods for solving the underlying PDE, are flexible to a wide range of problems, and are simple to implement, requiring minimal geometric information and processing to describe the boundary.
The IB method belongs to the broad category of methods known as embedded boundary (EB) methods, including the Immersed Interface [7], Ghost Fluid [8], and Volume Penalty methods [9]. These methods share a common feature: they enable solutions to PDE on nontrivial domains to be computed using efficient and robust structuredgrid discretizations; yet these methods differ largely in how boundary conditions are enforced and whether or not the solution is produced in the entirety of a simple domain. Methods which compute the solution everywhere in a dimensional rectangle admit the simplest discretizations and enable the use of highorder discretizations such as Fourier spectral methods. Unfortunately, this simplicity comes coupled with a fundamental difficulty: the analytic solution to these problems is rarely globally smooth on the entire domain. Consider the onedimensional Poisson problem on the periodic interval with Dirichlet boundary conditions for . Even if , the solution will typically display jumps in its derivative at the values and (see Figure 2). The lack of regularity in the analytic problem leads to loworder convergence in many numerical schemes, including the Immersed Boundary method.
The advantages of EB methods are substantial enough that significant effort has been expended on improving their accuracy [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. Two different approaches are generally taken. The first approach involves locally altering the discretization of the PDE in the vicinity of the boundary to accommodate the lack of smoothness in the solution. One example of this approach is the Immersed Interface method [7]. Such approaches are particularly useful for interface problems where the solution is required on both sides of the embedded boundary. When the solution is only needed on one side of the interface, a second approach may be taken in which variables are redefined outside of the domain of interest to obtain higher global regularity. Improved convergence rates are achieved as a natural consequence of the properties of the discretization scheme when applied to smooth problems. Variations of this basic idea have been used by the Fourier Continuation (FC) [24, 25, 23] and the Active Penalty (AP) [26] methods to provide high order solutions to PDE on general domains.
In this paper, we introduce a new method termed the Immersed Boundary Smooth Extension (IBSE) method. This method is a first step in remedying two deficiencies of Immersed Boundary methods:

For generic problems, the IB method produces discretizations that are only firstorder accurate in the vicinity of domain boundaries (or fluidstructure interfaces) [27].

The IB method is only able to specify Dirichlet boundary conditions (noslip, for fluid problems) without specialized interpolation schemes subordinate to the geometry [28]. Although this is not of interest when solving traditional IB fluidinterface problems, it allows the IB method to be generalized for solving other PDE (i.e. reactiondiffusion equations).
In this paper we restrict our attention to PDE set on stationary domains, and consider only PDE without a global divergence constraint. Directforcing IB methods produce solutions to PDE that are , with jumps in the normalderivative of the solution across the boundary, and converge at first order in the gridspacing [4, 5]. Drawing on ideas from the AP and FC methods, we use Fourier spectral methods to obtain a highly accurate discretization, while adding a volumetric forcing to nonphysical portions of the computational domain to force the solution to be globally . We will use the shorthand IBSE to refer to our method when we need to explicitly denote the global regularity of the solution that is enforced. Highorder accuracy is achieved naturally, as a simple consequence of the convergence properties of the Fourier transform.
The IBSE method retains the essential robustness and simplicity of the original IB method. All communication between the Lagrangian boundary and the underlying Cartesian grid is achieved by convolution with regularized functions or normal derivatives of those functions. This allows an absolute minimum of geometric information to be used. In the traditional IB method, only the position of the Lagrangian structure must be known; the IBSE method additionally requires normals to that structure and an indicator variable denoting whether Cartesian grid points lie inside or outside of the physical domain where the PDE is defined. Additionally, since normal derivatives can be accurately approximated by convolution with derivatives of regularized functions, Neumann and Robin boundary conditions can be imposed in the same way that Dirichlet boundary conditions are in directforcing IB approaches.
In contrast to other approaches based on the smooth extension of the forcing function or the solution [24, 25, 23, 26, 19], we do not extend the forcing function or solution from a previous timestep. Instead, we smoothly extend the unknown solution to the entire computational domain. This approach directly enforces smoothness of the solution, and it allows for highorder implicittime discretizations for parabolic equations and implicitexplicit discretizations for many nonlinear PDE. Remarkably, the coupled problem for the solution to an elliptic PDE and that solution’s smooth extension can be reduced to the solution of a relatively small dense system of equations (with size a small multiple of the number of points used to discretize the boundary), along with several FFTs. The dense system depends only on the boundary and the discretization, and so it can be formed and prefactored to allow for efficient timestepping of parabolic equations.
This paper is organized as follows. In Section 2, we introduce the method, ignoring the details of the numerical implementation. The fundamental contributions of this paper are contained in the modification to the IB discretization that enforces smoothness of the solution, and in the system of equations (given in Section 2.2) that allows for the simultaneous solution of the PDE along with the smooth extension of that unknown solution. In Section 3, we discuss the particular numerical implementation that we choose, detailing a fast algorithm for solving elliptic equations. In Section 4, we compute solutions to the Poisson problem in one and two dimensions, demonstrating highorder pointwise convergence: up to fourthorder for Dirichlet problems and thirdorder for Neumann problems. In Section 5, we discuss discretization of the heat equation, detailing a fast algorithm for implicit timestepping, and demonstrate fourthorder convergence in space and time. These numerical tests include direct comparisons to the Fourier Continuation [24] and Active Penalty [26] methods. Finally, in Section 6, we solve two nonlinear problems: a 2D Burgers’ equation with homogeneous Dirichlet boundary conditions, and the FitzhughNagumo equations with homogeneous Neumann boundary conditions. We demonstrate fourthorder convergence in space and time for Burgers’ equation and thirdorder convergence for the FitzhughNagumo equations.
2 Methods
Let , , and . We will assume that is smooth and does not selfintersect. Two typical domains are shown in Figure 1. We refer to as the physical domain, as the extension domain, and as the computational domain. We first consider the Poisson problem with Dirichlet boundary conditions: \cref@addtoresetequationparentequation
(1a)  
(1b) 
For now, we assume that and are smooth () functions defined in and respectively. Let be a smooth extension of to , that is, and for all . The direct forcing formulation of the Immersed Boundary method provides a way to solve Equation 1: \cref@addtoresetequationparentequation
(2a)  
(2b) 
Equation 1a is solved in the entire computational domain , while the boundary condition is enforced by the addition of a singular force supported on the boundary that acts as a Lagrange multiplier. This singular force term leads to jumps in the normal derivative of the solution across the boundary . This lack of smoothness restricts Immersed Boundary formulations such as Equation 2 to firstorder convergence in the vicinity of the boundary unless onesided discrete functions are used, which make the implementation of the method more difficult [10, 27].
We define an example onedimensional problem to illustrate the limited regularity of solutions produced by Equation 2. Let the computational domain be , where the onedimensional torus is identified with the periodic interval , and let the physical domain be given by the interval . The extension domain for this problem is . On this domain, we will solve the problem: \cref@addtoresetequationparentequation
(3a)  
(3b) 
As is chosen to be in , we may choose for all . The analytic solution to this problem is shown in Figure 2, along with the solution to the extended problem given by Equation 2. Note that the solution is only continuous despite the fact that . Choosing to be smooth leads to low regularity in the solution when the solution to the associated homogeneous problem is nontrivial, as is the case here. Without modification, direct discretizations of this problem will exhibit slow convergence due to the limited regularity of in the vicinity of the boundary.
Rather than choose a smooth extension to the forcing function , we choose an extension of the forcing function that gives a smooth solution on the entire computational domain. Let be any smooth extension of the unknown solution into . We can compute a forcing function associated with :
(4) 
The extended forcing function is then defined to be
(5) 
where denotes the characteristic function of the domain . Figure 2 shows the extended forcing function , along with the associated smooth solution to Equation 3. Since the solution is smooth, we expect that standard discretizations of Equation 2 will converge more rapidly when computed using the extension than when using the extension . We emphasize that the extended forcing function depends on the unknown solution .
This motivates us to define a reformulated version of Equation 2, with the extended forcing function given by as defined in Equation 5: \cref@addtoresetequationparentequation
(6a)  
(6b) 
Here is any smooth extension operator that satisfies \cref@addtoresetequationparentequation
(7a)  
(7b) 
When restricted to the physical domain , Equation 6a reduces to , the problem of physical interest. When restricted to the extension domain , Equation 6a reduces to . As long as , then Equation 7b ensures that , and hence in . Again by Equation 7b, in , and so in . Since by Equation 7a, then . This additional regularity of the solution allows standard discretizations of Equation 6 to converge at a faster rate than discretizations of Equation 2. Using and a simple Fourierspectral discretization, numerical solutions to Equation 6 should converge at [29].
In the remainder of this section, we lay out the remaining components needed to fully specify the IBSE method that do not depend on the particular choice of numerical implementation. In Section 2.1, we discuss how to smoothly extend a known function from to . In Section 2.2, we outline a system of equations that allows us to solve for and its smooth extension simultaneously. A numerical implementation of this method will be discussed in Section 3.
2.1 Smooth extension of a known function from to
One way to construct a smooth extension to a function is to solve a highorder PDE. Compared to the localized and effectively onedimensional extension strategies taken in [24, 26], the decision to extend a function by solving a fully dimensional PDE may appear to be needlessly complex. However, we will show in Sections 3 and 2.2 that this choice leads to a robust method that is straightforward to implement and requires relatively low computational cost.
In particular, to compute a extension to a function , we solve the following equation in the extension domain : \cref@addtoresetequationparentequation
(8a)  
(8b) 
Here denotes the normal derivative of on the boundary ; a computational formula for is given later in Equation 22. A simple choice of the operator is the polyharmonic operator . From an analytic perspective, this choice is sufficient, although we will show in Section 3.3 that other choices of will produce better results due to numerical issues. The smooth extension may then be constructed as
(9) 
Our choice of smoothed forcing only depends on in the extension region (see Equation 5), so there is no need to compute the actual extension since the values of in are irrelevant. For our purposes, we can think of an ‘extension’ as simply a function in the computational domain that shares its first normal derivatives with on the boundary , i.e. the function .
As with Equation 1, we may solve Equation 8 in by extending the problem to all of . This can be accomplished by adding singular forces supported along the boundary: \cref@addtoresetequationparentequation
(10a)  
(10b) 
The integral on the lefthand side of Equation 10b is notation for the action of the distribution on the smooth function . The boundary condition given in Equation 10b forces the solution to be , so there is no need to alter this formulation to provide additional regularity.
For notational convenience, we define the operators and by: \cref@addtoresetequationparentequation
(11a)  
(11b) 
In the language of the Immersed Boundary method, is an interpolation operator; is a spread operator. The notation is suggestive of the fact that these operators obey the adjoint property
(12) 
for all sufficiently smooth and , where the notation denotes the inner product on . Using the operators and , we may represent Equation 10 as
(13) 
2.2 A coupled system of equations for and its extension
Let the spread () and interpolation () operators be defined as \cref@addtoresetequationparentequation
(14a)  
(14b) 
Notice that and . We can now represent Equation 6 simply as \cref@addtoresetequationparentequation
(15a)  
(15b) 
where is defined by the extension equation
(16) 
along with the constraint that
(17) 
Equations 17, 16 and 15 can be combined into one system of equations:
(18) 
Equation 18 is the system of equations that allows the IBSE method to solve for and its extension simultaneously. The remainder of this paper will be concerned with the discretization and inversion of Equation 18 and numerical studies demonstrating the accuracy of solutions produced by those discretizations. We remark that Equation 18 is equally valid for operators of the Helmholtz type () that appear in the discretization of parabolic equations, which will be discussed in Section 5.
3 Numerical implementation
For concreteness, we will discuss the numerical implementation in the context of Fourier spectral methods, with the computational domain given by . We will discretize the domain using a regular Cartesian mesh with points discretizing each dimension: and . Differential operators are discretized in the usual way.
Few of the details depend upon the choice to use Fourier spectral methods, and they are chosen because of their simple implementation, computational speed, and high order of accuracy. Inversion of the elliptic operator , as well as the extension operator is also greatly simplified when using Fourier spectral methods. However, any discretization based on a regular Cartesian mesh could be used with minimal modifications to the method.
In order to fully describe a discretization and solution strategy to Equation 18, we must describe several key elements.

In Section 3.1, we discretize the spread (, ) and interpolation (, ) operators.

In Section 3.2 we define a new regularization of the function that is accurate and smooth.

In Section 3.3, we define our extension operator and show how this choice of controls the numerical conditioning of the extension problem.

In Section 3.4, we describe an efficient inversion strategy for the IBSE system given by Equation 18.

Finally, in Section 3.5, we provide implementational details and briefly discuss the computational complexity of our inversion scheme.
3.1 Discretization of , , , and
Let the boundary be parametrized by the function . In all examples in this manuscript, the boundary is onedimensional and closed; the single parameter is defined on the periodic interval . The spread operator is defined as
(19) 
The discrete version of this operator requires a regularized function and a discretization of the integral over . Construction of a regularized function with the properties necessary for the IBSE method is nontrivial; we will delay discussion of this choice until Section 3.2 and simply denote the regularized function as . Multivariate functions are computed as Cartesian products of the univariate and also denoted by . Discretization of the integral over is made by choosing quadrature nodes , equally spaced in the parameter interval so that and . Quadrature weights are computed at each quadrature node to be ; this is a spectrally accurate quadrature rule for the integral of smooth periodic functions on . The discrete spread operator maps functions sampled at points in to by
(20) 
We do not adopt explicit notation to distinguish between the analytic and discretized operators. The number of points in the quadrature is chosen so that . This choice of nodespacing is wider than that recommended for the traditional IB method [1] but has been observed empirically to be the optimal choice in other studies of directforcing IB methods [4]. The interpolation operator may be defined by the adjoint property , but we note that the discrete interpolation operator produces a discrete function
(21) 
Discrete integrals over are straightforward sums computed over the underlying uniform Cartesian mesh.
Normal derivatives of are computed by the formula [30]
(22) 
where the Einstein summation convention has been used to indicate sums over repeated indices and is computed as Cartesian products of the appropriate derivatives of the onedimensional . For example, in two dimensions is computed as
(23) 
Analogous to the definition of , we define the spread operator for the normal derivative, , as
(24) 
and again define the interpolation operator by the adjoint property . Analogous to the action of , the operator produces approximations of the normal derivative of a function defined on at the quadrature nodes of . Finally, the composite operator is defined by:
(25) 
while .
3.2 Construction of a smooth discretization of
The IBSE method is capable of producing solutions that converge at . In order to achieve this accuracy, our choice of regularized function must satisfy several conditions. For a general value of :

Its interpolation accuracy must be at least .

It must be at least , so that its normal derivative is continuous, allowing it to be used in the discretization of and .
In addition, we will require that the function has compact (and small) support so that the spread and interpolation operators , , , and may be rapidly applied. The commonly used ‘four point’ function is , has a support of four gridpoints, and produces approximations of , , and that are accurate to [31]. The use of this function regularization is sufficient for the implementation of IBSE, but it is not sufficient for higher order methods, due to both its limited accuracy and regularity.
In this manuscript, we discretize the IBSE method for , , and , corresponding to second, third, and fourth order accuracy in for Dirichlet problems. We therefore need a regularization of that is and has an interpolation accuracy of for smooth functions. We are not aware of any functions with these properties currently defined in the literature, and so we construct such a function here. For simplicity, we do not attempt to impose other conditions that are often imposed on functions, such as the evenodd condition or sum of squares condition [32]. The strategy that we follow is to choose a function with sufficient accuracy but limited regularity and convolve it against itself to increase its smoothness. A similar approach was used in [33] to generate a smoother version of the traditional Peskin four point function [1]. In order to provide an analytic formula, a base function with a simple functional form must be used. We start with the function
(26) 
which has interpolation accuracy, regularity, and a support of four gridpoints [34, 35]. Define:
(27) 
where denotes convolution. It is clear that , its support is gridpoints, and it is a simple exercise to show that convolution preserves interpolation accuracy. The formula for is given in Appendix A.
While our choice of function regularization limits the spatial discretization in this paper to fourthorder accuracy, this is not a fundamental limitation of the method. One way to generalize to higher orders would be to construct smoother and more accurate regularizations of the function similar to the one that we construct, which maintain finite support for computational efficiency. An alternative approach would be to use globally supported regularizations of the function that are (i.e. sinc or Gaussian functions) and make use of fast sinc transforms or fast Gaussian transforms in order to rapidly apply the spread and interpolation operators [36, 37].
Our construction of is likely not optimal, in that there probably exist function regularizations accurate to with a support of less than 16 gridpoints (e.g. a function with regularity and accuracy of with support of only six gridpoints is defined in [38]). The construction of such a function would be worthwhile, allowing for coarser discretizations of problems with closely spaced boundaries and faster application of the operators , , , and .
3.3 Choice of extension operator
Due to the nullspace of the periodic polyharmonic operator , we choose the operator used in the extension problem given by Equation 8 to be the invertible operator
(28) 
Here is a positive scalar function that depends on the smoothness of the extension () and the number of Fourier modes () used in the discretization of the problem. The function is chosen to mitigate the numerical condition number of the operator :
(29) 
Minimizing the condition number (by taking to be large) must be balanced against the need to resolve the intrinsic length scale introduced to the problem. In practice, we find that taking to be
(30) 
where is machine precision ( for double precision and for quadruple precision) and produces stable and accurate solutions. Numerical results are not particularly sensitive to the choice of (see Figure 3). This stability is demonstrated later (in Section 4.1) with highly accurate solutions shown for the IBSE, IBSE, and IBSE methods up to . All numerical results in this paper are produced using to construct .
For large values of and , the value of can substantially impact the overall accuracy of the IBSE algorithm. We demonstrate this effect in Figure 3 for , , and solutions to Equation 3. In Figure 3, we plot the error in the solution as a function of across 32 orders of magnitude. The minimum error observed is , at , which is within an order of magnitude of (denoted by a vertical grey line). Although there is some sensitivity to near the minimum, it is small: the error is bounded by over 6 orders of magnitude of . In Figure 3, we show a zoom of the extension near the domain boundary at for the solution produced using , which has an error of . The extension decays rapidly to , but is well resolved by the discretization: there are gridpoints between the boundary at and the peak at . In Figure 3, we show a zoom of the extension produced when , which gives a solution with an error of , four orders of magnitude worse than the error produced when . Although taking this large leads to a very wellconditioned operator , the extension decays very rapidly to , and the discretization is no longer able to fully resolve the length scales: there are only three gridpoints between the boundary at and the peak at . Finally, in Figure 3, we use . Although there are no fine lengthscales to be resolved, illconditioning in the operator leads to an error of .
3.4 Inversion of the IBSE system (18)
In this section, we consider solving Equation 18 for an invertible elliptic operator in place of . The case where is the periodic Laplacian is complicated by the nullspace of : even though Equation 18 is invertible, the method for inversion that we outline in this section is not directly applicable; details for the resolution of this problem are provided in Appendix B. The system of equations is
(31) 
This system is large: in two spatial dimensions, the number of equations in the system is . By block Gaussian elimination, we can find a Schurcomplement for this matrix, which we label :
(32) 
along with an associated system of equations for the Lagrange multipliers and :
(33) 
The size of is comparatively small, only square. This equation is the key to the efficiency of the algorithm, as it allows the Lagrange multipliers and to be computed rapidly without first solving for or . Once and are determined, we may solve for :
(34) 
and once is known, may be computed as
(35) 
The computation of and requires only a series of fast operations: FFTs, multiplies, and applications of the spread operators.
Rapid inversion of the system of equations for and given in Equation 33 is not a trivial task. Because we have restricted to problems set on stationary domains and the size of is small, it is feasible to proceed using dense linear algebra. We form by repeatedly applying it to basis vectors. This operation is expensive: for two dimensions it is in the total number of unknowns . The Schurcomplement is then factored by the LU algorithm provided by LAPACK [39]. Once this factorization is computed Equation 33 can be solved rapidly. This Schurcomplement depends only on the domain and the discretization, so the LUdecomposition can be reused to solve multiple problems on the same domain or in each timestep when solving timedependent problems. See Section 3.5 for a discussion of the computational complexity of the method and Table 2 for the actual numerical cost of the method applied to a test problem.
3.5 Summary of algorithm and numerical complexity
For a given geometry, choice of elliptic operator , and discretization size , the implementation of the IBSE method proceeds as follows.

Setup

The boundary is discretized as described in Section 3.1, and stencils for the evaluation of and its first normal derivatives are computed at all nodes of the discrete boundary in to allow for the rapid application of the , , , and operators.

The Schurcomplement matrix (Equation 32 for invertible , and Equation 70 for or other noninvertible ) is formed columnwise, by repeatedly applying it to basis vectors. For invertible , this task may be completed as follows (for noninvertible , the process is nearly identical):

One element of either or is set to , all other elements are set to .

These forces are spread to the Eulerian grid, by the application of and .

, and then , are computed by applying Equation 34 followed by Equation 35.

The quantities and are computed, and placed in the appropriate rows and column of the Schurcomplement.
Note that item (ii) corresponds to application of the rightmost matrix in Equation 32, item (iii) corresponds to inverting the middle matrix in Equation 32, while item (iv) corresponds to the application of the leftmost matrix in Equation 32. We remark that although this operation is expensive (see Table 2), it is trivially parallelizable, as the task of applying the Schurcomplement to each individual basis vector is independent.


The LU decomposition of the Schurcomplement is computed.

Optionally, these items may be saved, bypassing the expensive computation and factorization of the Schurcomplement in future uses of the IBSE method when applied to the same geometry, , and discretization size .


Solve

The righthand side of Equation 33 (or equivalently, Equation 70 when is not invertible) is computed.

Using the precomputed LU decomposition of the Schur complement matrix, the Lagrange multipliers and are found by solving Equation 33 (or equivalently, Equation 70 when is not invertible).

The smooth extension , and finally , are computed via Equations 35 and 34 (or equivalently, Equation 72 when is not invertible)

In items 1(b), 2(a), and 2(c), application of the FFT and IFFT is used to move variables between frequency and physical space, while differential operators are applied in frequency space and multiplications with characteristic functions are carried out in physical space.
The computational work for this method can be broken into two main tasks: (i) the setup items, listed above, which are dominated by the formation and factorization of the Schurcomplement defined in Equation 32 and (ii) the production of one solution given the Schurcomplement, its factorization, the body forcing , and the boundary conditions . We summarize the scaling of the algorithm in Table 1; see Table 2 in Section 5.1 for the actual numerical cost of the method applied to a test problem. In two dimensions, the cost of a solve scales like the cost of an FFT; in three dimensions the cost is slightly worse: in the total number of unknowns . The higher cost of the algorithm in threedimensions is due to the cost of factoring and inverting the Schurcomplement, however, these matrices are highly structured. It is quite possible that the cost of inversion could be reduced with an appropriately preconditioned iterative method, by directly exploiting the structure of the matrix, e.g. with an inverse Multipole Method [40], or by indirectly exploiting that structure using an algorithm like HODLR [41] or Algebraic Multigrid [42].
Dimension  Cost of FFT  Cost of Setup  Cost of Solve 

2  
3 
4 Results: Poisson equation
4.1 Onedimensional test problem
To demonstrate both the power and the limitations of the IBSE method, we compute the solution to the onedimensional example defined in Equation 3. Solutions are computed on grids ranging from to points, in both double and quadrupleprecision arithmetic, for the IB, IBSE, IBSE, and IBSE methods. Figure 4 shows a refinement study demonstrating the expected first, second, third, and fourthorder accuracy in for the IB, IBSE, IBSE, and IBSE methods, respectively. The IB method achieves an error of with grid points. Imposing additional smoothness on the solution allows this error to be matched by the IBSE method at , by the IBSE method at , and by the IBSE method at . For the IBSE method, this is a factor of nearly 10000 less gridpoints; equivalent to a factor of 100 million less gridpoints for twodimensional problems. In practice, obtaining solutions accurate to six digits with the traditional IB method is often impractical; with the IBSE method this kind of accuracy can be achieved on reasonably sized grids.
In double precision, all of the methods achieve the expected convergence rate up to some value of . The doubleprecision solutions begin to diverge from the quadrupleprecision solutions at for IBSE, at for IBSE, and at for IBSE. Using (see Section 3.3) to control the condition number of the extension operator provides sufficient numerical stability at all values of to enable the computation of solutions accurate to 1011 digits. In quadruple precision, all methods exhibit the expected order of convergence across a wide range of values of , although the fourth order method begins to fall short of the expected path of convergence for extremely fine grids ().
We note that this test problem, as well as the test problems in Sections 6.2, 6.1 and 5.1, are exterior problems set on periodic domains. Periodicity of solutions in the physical domain is enforced naturally in these problems through the use of the Fourier basis.
4.2 Twodimensional test problem: solution inside a circle
Let be the circle of radius centered at , and identify with the square . We will solve the Dirichlet problem: \cref@addtoresetequationparentequation
(36a)  
(36b) 
which has an analytic solution given by
(37) 
We solve this problem in doubleprecision arithmetic with the IBSE, IBSE, and IBSE methods for to . In Figure 5, we show error as a function of , indicating second, third, and fourthorder accuracy for IBSE, IBSE, and IBSE respectively.
Using the solutions to this problem, we demonstrate the property that the IBSE method produces globally smooth solutions. In Figure 6, we show the solution to Equation 36, along with its first, second, and third derivatives in the direction, produced using the IB, IBSE, IBSE, and IBSE methods with . For simplicity, all functions are shown only along the slice , and the intersection of this slice with the boundary is shown as gray vertical lines. Derivatives are computed spectrally using the FFT. We expect solutions produced by the IBSE method to exhibit discontinuities in the normal derivative across interfaces, and for all derivatives of lower order to be continuous. We see this expectation validated in the solutions presented in Figure 6: functions shown in plots with shading behind them in Figure 5(a) are at least continuous, all others are discontinuous. This property of global smoothness of the solutions enables the Fourierspectral discretization to obtain highorder accuracy.
4.3 Twodimensional test problem with Neumann boundary conditions
Traditional Immersed Boundary approaches are unable to provide solutions to Neumann problems since convolutionstyle estimation of the normal derivative is not convergent at the boundary due to the lack of regularity of the solution. The IBSE method does not have this limitation because of the additional smoothness of the solution. The only modification necessary to adapt Equation 18 to solve Neumann problems is to change the and operators that enforce and specify the boundary conditions. Consider the general Neumann problem: \cref@addtoresetequationparentequation
(38a)  
(38b) 
Define \cref@addtoresetequationparentequation
(39a)  
(39b) 
The analogue to Equation 18 is
(40) 
It is just as simple to modify the IBSE method to allow the solution of Robin problems. Rather than replacing the upperright and lowerleft operators in Equation 18 with and , they should instead be replaced with linear combinations of , , , and .
Let . We will solve the problem \cref@addtoresetequationparentequation
(41a)  
(41b) 
which has the analytic solution
(42) 
We solve this problem with the IBSE, IBSE, and IBSE methods, for which we expect first, second, and thirdorder convergence, respectively. The reason for the lower order of convergence than in Dirichlet problems is that our discretization of is one order less accurate than when acting on functions of the same regularity. Despite this, the solution produced by IBSE still has global regularity. Solutions are computed in double precision, for to . Figure 7 shows error as a function of . Second and thirdorder convergence is observed for the IBSE and IBSE methods. The IBSE method converges at a rate that is slightly less than firstorder.
5 Results: heat equation
We now consider solving the heat equation with Dirichlet boundary conditions: \cref@addtoresetequationparentequation
(43a)  
(43b) 
Since the IBSE method can directly invert elliptic equations, it can easily be adapted to provide highorder implicittime discretizations for the heat equation. Many implicit methods could be used; we choose BDF4 [43]:
(44) 
Applying this scheme to Equation 43b is equivalent to solving the problem: \cref@addtoresetequationparentequation
(45a)  
(45b) 
where we define: \cref@addtoresetequationparentequation
(46a)  
(46b) 
This may be solved using the algorithm described in Section 3.4. There are several comments worth making regarding this discretization:

The use of an implicit scheme allows for large timesteps to be taken. For all test problems presented in this section, we choose to be:
(47) where is the ceiling function, so that .

The fact that depends on implies that the Schurcomplement depends on as well. This means that the Schurcomplement must be reformed and refactored whenever is changed, complicating the use of a method that uses adaptive timestepping.

Modifying this formulation to solve the heat equation with Neumann and Robin boundary conditions may be done in the same way as shown for the Poisson equation in Section 4.3.
5.1 Solution in a periodic domain outside an obstacle
In this section we demonstrate highorder convergence to a heatequation set on a periodic domain with a simple obstacle in it. To provide direct numerical comparison with results from the Active Penalty (AP) method, the following problem is from Section 6.2 of [26]: \cref@addtoresetequationparentequation
(48a)  
(48b)  
(48c) 
where the physical domain is . Equation 48 has the analytic solution
(49) 
The initial condition is integrated from to by solving Equation 45 at each timestep. Startup values for BDF4 are computed by evaluating the analytic solution at , , and . In Figure 8, errors at for solutions generated with the IBSE, IBSE, and IBSE methods are shown, demonstrating second, third, and fourthorder convergence in space and time, respectively. On this plot, we also show errors from the second and thirdorder AP methods with , the finest resolution reported in [26]. Our method yields errors lower by several orders of magnitude.
In contrast to the AP method [26], we are able to use an implicit time discretization, enabling the use of large timesteps. The ability to take large timesteps efficiently comes at the cost of a significant precomputation for the IBSE method, a cost that the AP method does not incur. This precomputation prevents the IBSE method as described in this paper from being efficiently applied to moving boundary problems, a limitation not faced by the AP method.
In Table 2, we provide the computational time required by the IBSE, IBSE, and IBSE methods to setup the problem (dominated by the formation and factorization of Equation 32), the time required to take one timestep, and the number of timesteps taken to advance the equation from to . Computational times are given as a multiple of the time required to take one FFT. In our code, 8 FFTs are used per timestep regardless of ; the time required to take a timestep does not vary significantly across the methods. The precomputation time increases approximately linearly in , but even for the finest discretization presented, the precomputation time is not prohibitive: for and the IBSE method, the precomputation requires about 50 minutes using serial code that has not been carefully optimized.
IBSE  IBSE  IBSE  timesteps to  

prep time  timestep  prep time  timestep  prep time  timestep  
674  39  895  39  1451  43  2  
957  34  1776  39  2371  39  3  
1310  30  2207  31  3241  37  5  
1833  21  2833  24  3917  23  9  
3478  21  4493  22  5782  20  17  
5328  19  7521  19  10187  19  33  
9951  18  13707  18  18686  18  66 
5.2 Solution inside a parametrically defined region
In this section, we demonstrate highorder convergence to a heat equation set inside a parametrically defined domain. This problem is from Section 6.1 of [24], allowing direct comparison of the IBSE method with the Fourier Continuation method [24]. For convenience, we rescale the domain of the problem. We will solve \cref@addtoresetequationparentequation
(50a)  
(50b)  
(50c) 
where
(51) 
The domain is the region inside the boundary defined by the parametric equations: \cref@addtoresetequationparentequation
(52a)  
(52b) 
for . Equation 50 has the analytic solution . The initial condition is integrated from