1 Introduction

A high-order Nyström discretization scheme for boundary integral equations defined on rotationally symmetric surfaces

P. Young, S.  Hao, P.G. Martinsson

Abstract: A scheme for rapidly and accurately computing solutions to boundary integral equations (BIEs) on rotationally symmetric surfaces in is presented. The scheme uses the Fourier transform to reduce the original BIE defined on a surface to a sequence of BIEs defined on a generating curve for the surface. It can handle loads that are not necessarily rotationally symmetric. Nyström discretization is used to discretize the BIEs on the generating curve. The quadrature is a high-order Gaussian rule that is modified near the diagonal to retain high-order accuracy for singular kernels. The reduction in dimensionality, along with the use of high-order accurate quadratures, leads to small linear systems that can be inverted directly via, e.g., Gaussian elimination. This makes the scheme particularly fast in environments involving multiple right hand sides. It is demonstrated that for BIEs associated with the Laplace and Helmholtz equations, the kernel in the reduced equations can be evaluated very rapidly by exploiting recursion relations for Legendre functions. Numerical examples illustrate the performance of the scheme; in particular, it is demonstrated that for a BIE associated with Laplace’s equation on a surface discretized using points, the set-up phase of the algorithm takes 1 minute on a standard laptop, and then solves can be executed in 0.5 seconds.

1. Introduction

The premise of the paper is that it is much easier to solve a boundary integral equation (BIE) defined on a curve in than one defined on a surface in . With the development of high order accurate Nyström discretization techniques [14, 2, 15, 5, 13], it has become possible to attain close to double precision accuracy in 2D using only a very moderate number of degrees of freedom. This opens up the possibility of solving a BIE on a rotationally symmetric surface with the same efficiency since such an equation can in principle be written as a sequence of BIEs defined on a generating curve. However, there is a practical obstacle: The kernels in the BIEs on the generating curve are given via Fourier integrals that cannot be evaluated analytically. The principal contribution of the present paper is to describe a set of fast methods for constructing approximations to these kernels.

1.1. Problem formulation

This paper presents a numerical technique for solving boundary integral equations (BIEs) defined on axisymmetric surfaces in . Specifically, we consider second kind Fredholm equations of the form


under two assumptions: First, that is a surface in obtained by rotating a curve about an axis. Second, that the kernel is invariant under rotation about the symmetry axis in the sense that


where and are cylindrical coordinates for and , respectively,


see Figure 1. Under these assumptions, the equation (1), which is defined on the two-dimensional surface , can via a Fourier transform in the azimuthal variable be recast as a sequence of equations defined on the one-dimensional curve . To be precise, letting , , and denote the Fourier coefficients of , , and , respectively (so that (9), (10), and (11) hold), the equation (1) is equivalent to the sequence of equations


Whenever can be represented with a moderate number of Fourier modes, the formula (5) provides an efficient technique for computing the corresponding modes of . The conversion of (1) to (5) appears in, e.g., [20], and is described in detail in Section 2. Note that the conversion procedure does not require the data function to be rotationally symmetric.

1.2. Applications and prior work

Equations like (1) arise in many areas of mathematical physics and engineering, commonly as reformulations of elliptic partial differential equations. Advantages of a BIE approach include a reduction in dimensionality, often a radical improvement in the conditioning of the mathematical equation to be solved, a natural way of handling problems defined on exterior domains, and a relative ease in implementing high-order discretization schemes, see, e.g., [3].

The observation that BIEs on rotationally symmetric surfaces can conveniently be solved by recasting them as a sequence of BIEs on a generating curve has previously been exploited in the context of stress analysis [4], scattering [9, 17, 22, 23, 24], and potential theory [12, 19, 20, 21]. Most of these approaches have relied on collocation or Galerkin discretizations and have generally used low-order accurate discretizations.

1.3. Kernel evaluations

A complication of the axisymmetric formulation is the need to determine the kernels in (5). Each kernel is defined as a Fourier integral of the original kernel function in the azimuthal variable in (2), cf. (8), that cannot be evaluated analytically, and would be too expensive to approximate via standard quadrature techniques. The FFT can be used to a accelerate the computation in certain regimes. For many points, however, the function which is to be transformed is sharply peaked, and the FFT would at these points yield inaccurate results.

1.4. Principal contributions of present work

This paper resolves the difficulty of computing the kernel functions described in Section 1.3. For the case of kernels associated with the Laplace equation, it provides analytic recursion relations that are valid precisely in the regions where the FFT loses accuracy. The kernels associated with the Helmholtz equation can then be obtained via a perturbation technique.

The paper also describes a high-order Nyström discretization of the BIEs (5) that provides far higher accuracy and speed than previously published methods. The discretization scheme converges fast enough that for simple generating curves, a relative accuracy of is obtained using as few as a hundred points, cf. Section 8. The rapid convergence of the discretization leads to linear systems of small size that can be solved directly via, e.g., Gaussian elimination, making the algorithm particularly effective in environments involving multiple right hand sides or when the linear system is challenging for iterative solvers (as happens for many scattering problems).

Finally, the efficient techniques for evaluating the fundamental solutions to the Laplace and Helmholtz equations in an axisymmetric environment have applications beyond solving boundary integral equations, for details see Section 7.

1.5. Asymptotic costs

To describe the asymptotic complexity of the method, we let denote the total number of discretization points, and assume that as grows, the number of Fourier modes required to resolve the solution scales proportionate to the number of discretization points required along the generating curve . Then the asymptotic cost of solving (1) for a single right-hand side is . If additional right hand sides are given, any subsequent solve requires only operations. Numerical experiments presented in Section 8 indicate that the constants of proportionality in the two estimates are moderate. For instance, in a simulation with , the scheme requires 1 minute for the first solve, and seconds for each additional right hand side (on a standard laptop).

Observe that since a high-order discretization scheme is used, even complicated geometries can be resolved to high accuracy with a moderate number of points.

Figure 1. The axisymmetric domain generated by the curve .

1.6. Outline

Section 2 provides details on the conversion of a BIE on a rotationally symmetric surface to a sequence of BIEs on a generating curve. Section 3 describes a high order methods for discretizing a BIE on a curve. Section 4 summarizes the algorithm and estimates its computational costs. Section 5 describes how to rapidly evaluate the kernels associated with the Laplace equation, and then Section 6 deals with the Helmholtz case. Section 7 describes other applications of the kernel evaluation techniques. Section 8 illustrates the performance of the proposed method via numerical experiments.

2. Fourier Representation of BIE

Consider the BIE (1) under the assumptions on rotational symmetry stated in Section 1.1 (i.e.  is a rotationally symmetric surface generated by a curve and that is a rotationally symmetric kernel). Cylindrical coordinates are introduced as specified in (3). We write where is the one-dimensional torus, usually parameterized by .

2.1. Separation of Variables

We define for the functions , , and via


Formulas (6), (7), and (8) define , , and as the coefficients in the Fourier series of the functions , , and about the azimuthal variable,


To determine the Fourier representation of (1), we multiply the equation by and integrate over . Equation (1) can then be said to be equivalent to the sequence of equations


Invoking (11), we evaluate the bracketed factor in (12) as


Inserting (13) into (12) and executing the integration of over , we find that (1) is equivalent to the sequence of equations


For future reference, we define for the boundary integral operators via


Then equation (14) can be written


When each operator is continuously invertible, we write the solution of (1) as


2.2. Truncation of the Fourier series

When evaluating the solution operator (17) in practice, we will choose a truncation parameter , and evaluate only the lowest Fourier modes. If is chosen so that the given function is well-represented by its lowest Fourier modes, then in typical environments the solution obtained by truncating the sum (17) will also be accurate. To substantiate this claim, suppose that is a given tolerance, and that has been chosen so that


We define an approximate solution via


From Parseval’s identity, we then find that the error in the solution satisfies

It is typically the case that the kernel has enough smoothness that the Fourier modes decay as . Then as and . Thus, an accurate approximation of leads to an approximation in that is of the same order of accuracy. Figure 4 illustrates how fast this convergence is for Laplace’s equation (note that in the case illustrated, the original equation is , and it is shown that ).

3. Nyström discretization of BIEs on the generating curve

We discretize the BIEs (5) defined on the generating curve using a Nyström scheme. In describing the scheme, we keep the formulas uncluttered by discussing a generic integral equation

where is a simple smooth curve in the plane, and is a weakly singular kernel function.

3.1. Quadrature nodes

Consider a quadrature rule on with nodes and weights . In other words, for a sufficiently smooth function on ,


For the experiments in this paper, we use a composite Gaussian rule with points per panel. Such a rule admits for local refinement, and can easily be modified to accommodate contours with corners that are only piece-wise smooth.

3.2. A simplistic Nyström scheme

The Nyström discretization of (20) corresponding to a quadrature with nodes takes the form


where are coefficients such that


The solution of (21) is a vector such that each is an approximation to .

A simplistic way to construct coefficients so that (22) holds is to simply apply the rule (20) to each function whence


This generally results in low accuracy since the kernel has a singularity at the diagonal. However, the formula (23) has the great advantage that constructing each costs no more than a kernel evaluation; we seek to preserve this property for as many elements as possible.

3.3. High-order accurate Nyström discretization

It is possible to construct a high-order discretization that preserves the simple formula (23) for the vast majority of coefficients [13]. The only coefficients that need to be modified are those for which the target point is near1 the panel holding . In this case, is conceptually constructed as follows: First map the pointwise values to their unique polynomial interpolant on , then integrate this polynomial against the singular (or sharply peaked) functions using the quadratures of [15] Operationally, the end result is that is given by


In (24), is a small integer (roughly equal to the order of the Gaussian quadrature), the numbers are coefficients that depend on but not on , and are points on located in the same panel as (in fact, when and belong to the same panel, so the number of kernel evaluations required is less than it seems). For details, see [13].

4. The full algorithm

4.1. Overview

At this point, we have shown how to convert a BIE defined on an axisymmetric surface in to a sequence of equations defined on a curve in (Section 2), and then how to discretize each of these reduced equations (Section 3). Putting the components together, we obtain the following algorithm for solving (1) to within some preset tolerance :

  1. Given , determine a truncation parameter such that .

  2. Fix a quadrature rule for with nodes and form for each Fourier mode the corresponding Nyström discretization as described in Section 3. The number of nodes must be picked to meet the computational tolerance . Denote the resulting coefficient matrices .

  3. Evaluate via the FFT the terms (as defined by (6)) for .

  4. Solve the equation for .

  5. Construct using formula (19) evaluated via the FFT.

The construction of the matrices in Step ii can be accelerated using the FFT (as described in Section 4.2), but even with such acceleration, it is typically by a wide margin the most expensive part of the algorithm. However, this step needs to be performed only once for any given geometry. The method therefore becomes particularly efficient when (1) needs to be solved for a sequence of right-hand sides. In this case, it may be worth the cost to pre-compute the inverse (or LU-factorization) of each matrix .

4.2. Cost of computing the coefficient matrices

For each of the Fourier modes, we need to construct an matrix with entries . These entries are derived from the kernel functions defined by (8). Note that whenever , the function is , but that as it develops a progressively sharper peak around .

For two nodes and that are “not near” (in the sense defined in Section 3.3) the matrix entries are given by the formula


where is given by (8). Using the FFT, all entries can be evaluated at once in operations. The FFT implicitly evaluates the integrals (8) via a trapezoidal rule which is highly accurate since the points and are well-separated on the curve .

For two nodes and that are not well-separated on , evaluating is dicier. The first complication is that we must now use the corrected formula, cf. (24),


The second complication is that the FFT acceleration for computing the kernels jointly no longer works since the integrand in (8) is too peaked for the simplistic trapezoidal rule implicit in the FFT. Fortunately, it turns out that for the kernels we are most interested in (the single and double layer kernels associated with the Laplace and Helmholtz equations), the sequence can be evaluated very efficiently via certain recurrence relations as described in Sections 5 and 6 (for the Laplace and Helmholtz equations, respectively). Happily, the analytic formulas are stable precisely in the region where the FFT becomes inaccurate.

4.3. Computational Costs

The asymptotic cost of the algorithm described in Section 4.1 will be expressed in terms of the number of Fourier modes required, and the number of discretization points required along . The total cost can be split into three components:

  1. Cost of forming the matrices : We need to form matrices, each of size . For entries in each matrix, the formula (25) applies and using the FFT, all entries can be computed at once at cost . For entries close to the diagonal, the formula (26) applies, and all the entries can be computed at once at cost using the recursion relations in Sections 5 and 6. The total cost of this step is therefore .

  2. Cost of transforming functions from physical space to Fourier space and back: The boundary data must be decomposed into Fourier modes , and after the linear systems have been solved, the Fourier modes must be transformed back to physical space. The asymptotic cost is .

  3. Cost of solving the linear systems : Using a direct solver such as Gaussian elimination, the asymptotic cost is .

We make some practical observations:

  • The cost of executing FFTs is minimal and is dwarfed by the remaining costs.

  • The scheme is highly efficient in situations where the same equation needs to be solved for a sequence of different right hand sides. In this situation, one factors the matrices once at cost , and then the cost of processing an additional right hand side is only with a very small constant of proportionality.

  • To elucidate the computational costs, let us express them in terms of the total number of discretization points under the simplifying assumption that . Since , we find and . Then:

    Cost of setting up linear systems:
    Cost of the first solve:
    Cost of subsequent solves:

    We observe that even though the asymptotic cost of forming the linear systems is less than the cost of factoring the matrices, the set-up phase tends to dominate unless is large. Moreover, the cost of the subsequent solves has a very small constant of proportionality.

  • The high order discretization employed achieves high accuracy with a small number of points. In practical terms, this means that despite the scaling, the scheme is very fast even for moderately complicated geometries.

  • The system matrices often have internal structure that allow them to be inverted using “fast methods” such as, e.g., those in [18]. The cost of inversion and application can in fact be accelerated to near optimal complexity.

5. Accelerations for the Single and Double Layer Kernels Associated with Laplace’s Equation

This section describes an efficient technique based on recursion relations for evaluating the kernel , cf. (8), when is either the single or double layer kernel associated with Laplace’s equation.

5.1. The Double Layer Kernels of Laplace’s Equation

Let be a bounded domain whose boundary is given by a smooth surface , let denote the domain exterior to , and let be the outward unit normal to . Consider the interior and exterior Dirichlet problems of potential theory [11],

(27) (interior Dirichlet problem)
(28) (exterior Dirichlet problem)

The solutions to (27) and (28) can be written in the respective forms


where is a boundary charge distribution that can be determined using the boundary conditions. The point can be placed at any suitable location in . The resulting equations are


where in (31) and (32).

Remark 1.

There are other integral formulations for the solution to Laplace’s equation. The double layer formulation presented here is a good choice in that it provides an integral operator that leads to well conditioned linear systems. However, the methodology of this chapter is equally applicable to single-layer formulations that lead to first kind Fredholm BIEs.

5.2. Separation of Variables

Using the procedure given in Section 2, if , then (27) and (28) can be recast as a series of BIEs defined along . We express in cylindrical coordinates as



Then for a point , the kernel of the internal Dirichlet problem can be expanded as


Similarly, the kernel of the external Dirichlet problem can be written as


where has been written in cylindrical coordinates as . With the expansions of the kernels available, the procedure described in Section 4 can be used to solve (31) and (32) by solving




respectively for . Note that the kernels and contain a log-singularity as .

5.3. Evaluation of Kernels

The values of and for need to be computed efficiently and with high accuracy to construct the Nyström discretization of (33) and (34). Note that the integrands of and are real valued and even functions on the interval . Therefore, can be written as


Note that can be written in a similar form.

This integrand is oscillatory and increasingly peaked at the origin as approaches . As long as and as well as and are well separated, the integrand does not experience peaks near the origin, and as mentioned before, the FFT provides a fast and accurate way for calculating and .

In regimes where the integrand is peaked, the FFT no longer provides a means of evaluating and with the desired accuracy. One possible solution to this issue is applying adaptive quadrature to fully resolve the peak. However, this must be done for each value of required and becomes prohibitively expensive if is large.

Fortunately, an analytical solution to (35) exists. As noted in [8], the single-layer kernel can be expanded with respect to the azimuthal variable as


is the half-integer degree Legendre function of the second kind, and

To find an analytical form for (35), first note that in cylindrical coordinates the double-layer kernel can be written in terms of the single-layer kernel,

The coefficients of the Fourier series expansion of the double-layer kernel are then given by , which can be written using the previous equation as

To utilize this form of , set and note that

where and are the complete elliptic integrals of the first and second kinds, respectively. The first two relations follow immediately from the definition of and the relations for the Legendre functions of the second kind can be found in [1]. With these relations in hand, the calculation of for can be done accurately and efficiently when and as well as and are in close proximity. The calculation of can be done analogously.

Remark 2.

Note that the forward recursion relation for the Legendre functions is unstable when . In practice, the instability is mild when is near and the recursion relation can still be employed to accurately compute values in this regime. Additionally, if stability becomes an issue, Miller’s algorithm [10] can be used to calculate the values of the Legendre functions using the backwards recursion relation, which is stable for .

6. Fast Kernel Evaluation for the Helmholtz Equation

Section 5 describes how to efficiently evaluate the kernels as defined by (8) for kernels associated with Laplace’s equation. This section generalizes these methods to a broad class of kernels that includes the single and double layer kernels associated with the Helmholtz equation.

6.1. Rapid Kernel Calculation via Convolution

Consider a kernel of the form