A highorder 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 highorder Gaussian rule that is modified near the diagonal to retain highorder accuracy for singular kernels. The reduction in dimensionality, along with the use of highorder 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 setup 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
(1) 
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
(2) 
where and are cylindrical coordinates for and , respectively,
(3)  
(4) 
see Figure 1. Under these assumptions, the equation (1), which is defined on the twodimensional surface , can via a Fourier transform in the azimuthal variable be recast as a sequence of equations defined on the onedimensional 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
(5) 
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 highorder 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 loworder 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 highorder 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 righthand 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 highorder discretization scheme is used, even complicated geometries can be resolved to high accuracy with a moderate number of points.
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 onedimensional torus, usually parameterized by .
2.1. Separation of Variables
We define for the functions , , and via
(6)  
(7)  
(8) 
Formulas (6), (7), and (8) define , , and as the coefficients in the Fourier series of the functions , , and about the azimuthal variable,
(9)  
(10)  
(11) 
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
(12) 
Invoking (11), we evaluate the bracketed factor in (12) as
(13) 
Inserting (13) into (12) and executing the integration of over , we find that (1) is equivalent to the sequence of equations
(14) 
For future reference, we define for the boundary integral operators via
(15) 
Then equation (14) can be written
(16) 
When each operator is continuously invertible, we write the solution of (1) as
(17) 
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 wellrepresented 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
(18) 
We define an approximate solution via
(19) 
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 ,
(20) 
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 piecewise smooth.
3.2. A simplistic Nyström scheme
The Nyström discretization of (20) corresponding to a quadrature with nodes takes the form
(21) 
where are coefficients such that
(22) 
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
(23) 
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. Highorder accurate Nyström discretization
It is possible to construct a highorder 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 near
(24) 
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 :

Given , determine a truncation parameter such that .

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 .

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

Solve the equation for .

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 righthand sides. In this case, it may be worth the cost to precompute the inverse (or LUfactorization) 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
(25) 
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 wellseparated on the curve .
For two nodes and that are not wellseparated on , evaluating is dicier. The first complication is that we must now use the corrected formula, cf. (24),
(26) 
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:

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 .

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 .

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 setup 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
(29)  
(30) 
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
(31)  
(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 singlelayer 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
Further,
and
Then for a point , the kernel of the internal Dirichlet problem can be expanded as
where
Similarly, the kernel of the external Dirichlet problem can be written as
with
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
(33) 
and
(34) 
respectively for . Note that the kernels and contain a logsingularity 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
(35) 
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 singlelayer kernel can be expanded with respect to the azimuthal variable as
where
is the halfinteger degree Legendre function of the second kind, and
To find an analytical form for (35), first note that in cylindrical coordinates the doublelayer kernel can be written in terms of the singlelayer kernel,
The coefficients of the Fourier series expansion of the doublelayer 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
(36) 
where