{'fraction_non_alphanumeric': 0.0412509534706331, 'fraction_numerical': 0.020747520976353927, 'mean_word_length': 4.2593068035943515, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Direct Solver for the Rapid Solution of Boundary Integral Equations on Axisymmetric Surfaces in Three Dimensions

A scheme for rapidly and accurately computing solutions to boundary integral equations (BIEs) on rotationally symmetric surfaces in three dimensions 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. Nystrom discretization is used to discretize the BIEs on the generating curve. The quadrature used 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 Laplace's equation, 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 320 000 points, the set-up phase of the algorithm takes 2 minutes on a standard desktop, and then solves can be executed in 0.5 seconds.


Introduction
This paper presents a numerical technique for solving boundary integral equations (BIEs) defined on axisymmetric surfaces in R 3 . Specifically, we consider second kind Fredholm equations of the form under two assumptions: First, that Γ is a surface in R 3 obtained by rotating a curve γ about an axis. Second, that the kernel k is invariant under rotation about the symmetry axis in the sense that (1.2) k(x, x ) = k(θ − θ , r, z, r , z ), where (r, z, θ) and (r , z , θ ) are cylindrical coordinates for x and x , respectively, x = (r cos θ, r sin θ, z), x = (r cos θ , r sin θ , z ), see Figure 1.1. Under these assumptions, the equation (1.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 one-dimensional curve γ. To be precise, letting σ n , f n , and k n denote the Fourier coefficients of σ, f , and k, respectively (so that (2.6), (2.7), and (2.8) hold), the equation (1.1) is equivalent to the sequence of equations (1.3) σ n (r, z) + √ 2π γ k n (r, z, r , z ) σ n (r , z ) r dl(r , z ) = f n (r, z), (r, z) ∈ γ, n ∈ Z.
Whenever f can be represented with a moderate number of Fourier modes, the formula (1.3) provides an efficient technique for computing the corresponding modes of σ. The conversion of (1.1) to (1.3) appears in, e.g., [19], and is described in detail in Section 2.
Equations of the type (1.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., [2].
The numerical solution of BIEs such as (1.1) poses certain difficulties, the foremost being that the discretizations generally involve dense matrices. Until the 1980s, this issue often times made it prohibitively expensive to use BIE formulations as numerical tools. However, with the advent of "fast" algorithms (the Fast Multipole Method [9,10], panel clustering [13], etc.) for matrix-vector multiplication and the inversion of dense matrices arising from the discretization of BIE operators, these problems have largely been overcome for problems in two dimensions. This is not necessarily the case in three dimensions; issues such as surface representation and the construction of quadrature rules in a three dimensional environment still pose unresolved questions. The point of recasting the single BIE (1.1) defined on a surface as the sequence of BIEs (1.3) defined on a curve is in part to avoid these difficulties in discretizing surfaces, and in part to exploit the exceptionally high speed of the Fast Fourier Transform (FFT).
The reduction of (1.1) to (1.3) is only applicable when the geometry of the boundary is axisymmetric, but presents no such restriction in regard to the boundary load. Formulations of this kind have been known for a long time, and have been applied to problems in stress analysis [3], scattering [7,16,21,22,23], and potential theory [12,18,19,20]. Most of these approaches have relied on collocation or Galerkin discretizations and have generally relied on low-order accurate discretizations. A complication of the axisymmetric formulation is the need to determine the kernels k n for a large number of Fourier modes n, since direct integration of (1.2) through the azimuthal variable tends to be prohibitively expensive. When k is smooth, this calculation can rapidly be accomplished using the FFT, but when k is near-singular, other techniques are required (quadrature, local refinement, etc.) that can lead to significant slowdown in the construction of the linear systems.
The technique described in this paper improves upon previous work in terms of both accuracy and speed. The gain in accuracy is attained by constructing a high-order quadrature scheme for kernels with integrable singularities. This quadrature is obtained by locally modifying a Guassian quadrature scheme, in a manner similar to that of [4,5]. Numerical experiments indicate that for simple surfaces, a relative accuracy of 10 −10 is obtained using as few as a hundred points along the generating curve. 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 and when the linear system is ill-conditioned. To describe the asymptotic complexity of the method, we need to introduce some notation. We let N P denote the number of panels used to discretize the generating curve γ, we let N G denote the number of Gaussian points in each panel, and we let N F denote the number of Fourier modes included in the calculation. Splitting the computational cost into a "set-up" cost that needs to be incurred only once for a given geometry and given discretization parameters, and a "solve" cost representing the time required to process each right hand side, we have The technique described gets particularly efficient for problems of the form (1.1) in which the kernel k is either the single or the double layer potential associated with Laplace's equation. We demonstrate that for such problems, it is possible to exploit recursion relations for Legendre functions to very rapidly construct the Fourier coefficients k n in (1.3). This reduces the computational complexity of the setup (which requires the construction of a sequence of dense matrices) from (1.4) to Numerical experiments demonstrate that for a problem with N P = 80, N G = 10, and N F = 400 (for a total of 80×10×400 = 320 000 degrees of freedom), this accelerated scheme requires only 2.2 minutes for the setup, and 0.46 seconds for each solve when implemented on a standard desktop PC.
The technique described in this paper can be accelerated further by combining it with a fast solver applied to each of the equations in (1.3), such as those based on the Fast Multipole Method, or the fast direct solver of [17]. This would result in a highly accurate scheme with near optimal complexity. The paper is organized as follows: Section 2 describes the reduction of (1.1) to (1.3) and quantifies the error incurred by truncating the Fourier series. Section 3 presents the Nyström discretization of the reduced equations using high-order quadrature applicable to kernels with integrable singularities, and the construction of the resulting linear systems. Section 4 summarizes the algorithm for the numerical solution of (1.3) and describes its computational costs. Section 5 presents the application of the algorithm for BIE formulations of Laplace's equation and describes the rapid calculation of k n in this setting. Section 6 presents numerical examples applied to problems from potential theory, and Section 7 gives conclusions and possible extensions and generalizations.

Fourier representation of BIE
2.1. Problem formulation. Suppose that Γ is a surface in R 3 obtained by rotating a smooth contour γ about a fixed axis and consider the boundary integral equation In this section, we will demonstrate that if the kernel k is rotationally symmetric in a sense to be made precise, then by taking the Fourier transform in the azimuthal variable, (2.1) can be recast as a sequence of BIEs defined on the curve γ. To this end, we introduce a Cartesian coordinate system in R 3 with the third coordinate axis being the axis of symmetry. Then cylindrical coordinates (r, z, θ) are defined such that

Separation of variables.
We define for n ∈ Z the functions f n , σ n , and k n via 2π k(θ, r, z, r , z ) dθ.
k(x, x ) = k(θ − θ , r, z, r , z ) = n∈Z e in(θ−θ ) √ 2π k n (r, z, r , z ). (2.8) To determine the Fourier representation of (2.1), we multiply the equation by e −inθ / √ 2π and integrate θ over T (for our purposes, we can think of T as simply the interval [−π, π]). Equation (2.1) can then be said to be equivalent to the sequence of equations Invoking (2.8), we evaluate the bracketed factor in (2.9) as Inserting (2.10) into (2.9) and executing the integration of θ over T, we find that (2.1) is equivalent to the sequence of equations For future reference, we define for n ∈ Z the boundary integral operators K n via Then equation (2.11) can be written When each operator I + K n is continuously invertible, we can write the solution of (2.1) as

2.3.
Truncation of the Fourier series. When evaluating the solution operator (2.14) in practice, we will choose a truncation parameter N F , and evaluate only the lowest 2N F + 1 Fourier modes. If N F is chosen so that the given function f is well-represented by its lowest 2N F + 1 Fourier modes, then in typical environments the solution obtained by truncating the sum (2.14) will also be accurate. To substantiate this claim, suppose that ε is a given tolerance, and that N F 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 k(x, x ) has sufficient smoothness such that the Fourier modes k n (r, z, r , z ) decay as n → ∞. Then ||K n || → 0 as n → ∞ and ||(I + K n ) −1 || → 1. Thus, an accurate approximation of f leads to an approximation in σ that is of the same order of accuracy. Figure 6.3 illustrates that when k is the double layer kernel associated with the Laplace equation, and γ is a simple curve, then ||(I + K n ) −1 || → 1 with rapid convergence.

Discretization of BIEs in two dimensions
The technique in Section 2 reduces the BIE (1.1) defined on an axisymmetric surface Γ = γ × T contained in R 3 , to a sequence of BIEs defined on the curve γ contained in R 2 . These equations take the form where the kernel k n is defined as in (2.5). In this section, we describe some standard techniques for discretizing an equation such as (3.1). For simplicity, we limit attention to the case where γ is a smooth closed curve, but extensions to non-smooth curves can be handled by slight variations of the techniques described here, [4,5,14].
To keep our formulas uncluttered, we suppress the parameterization of the curve and the dependence on n and introduce a new kernel as well as the functions ϕ(t) = σ(τ (t)) and ψ(t) = f (τ (t)).
Then techniques for solving where ψ is given and ϕ is to be determined, will be equally applicable to (3.2).
3.2. Nyström method. We will discretize (3.4) via Nyström discretization on standard Gaussian quadrature nodes, see [2]. To this end, we divide the interval Ω = [0, T ] into a disjoint partition of N P intervals, where each Ω p is a subinterval called a panel. On each panel Ω p , we place the nodes of a standard N G -point Gaussian quadrature rule {t The idea is now to enforce (3.4) at each of the N P N G nodes: To obtain a numerical method, suppose that we can construct for p, q ∈ {1, 2, . . . , Then the Nyström method is given by solving the linear system where A is a matrix formed by N P × N P blocks, each of size N G × N G . We let A (p,q) denote the block of A representing the interactions between the panels Ω p and Ω q .
3.3. Quadrature and interpolation. We need to determine the numbers A (p,q) i,j such that (3.5) holds. The detailed construction is given in Section 3.4, and utilizes some well-known techniques of quadrature and interpolation, which we review in this section.

Standard Gaussian quadratures. Given an interval [0, h] and a positive integer
whenever g is a polynomial of degree at most 2N G − 1, and such that whenever g is a function with 2N G continuous derivatives, see [1].

3.3.2.
Quadrature rules for singular functions. Now suppose that given an interval [0, h] and a point t ∈ [−h, 2h], we seek to integrate over [0, h] functions g that take the form where φ 1 and φ 2 are polynomials of degree at most 2N G −1. Standard Gaussian quadrature would be highly inaccurate if applied to integrate (3.7). Rather, seek a N G -node quadrature that will evaluate exactly. Techniques for constructing such generalized quadratures are readily available in the literature, see for example [15]. These quadratures will be of degree 2N G − 1, just as with standard Gaussian quadratures and exhibit comparable accuracy, although in general N G > N G . The generalized quadratures used in this paper were determined using the techniques of [15], and can be found in the appendix. We observe that the quadrature nodes constructed by such methods are typically different from the nodes of the standard Gaussian quadrature. This complicates the the construction of the matrix A, as described in Section 3.4.
denote the nodes of a N G -point Gaussian quadrature rule on the interval [0, h]. If the values of a polynomial g of degree at most N G − 1 are specified at these nodes, the entire polynomial g can be recovered via the formula where the functions L j are the Lagrange interpolating polynomials If g is a smooth function with N G continuous derivatives that is not a polynomial, then the Lagrange interpolant provides an approximation to g satisfying 3.4. Constructing the matrix A. Using the tools reviewed in Section 3.3, we are now in position to construct numbers A When Ω p and Ω q are well separated, the integrand in (3.9) is smooth, and our task is easily solved using standard Gaussian quadrature (as described in Section 3.3.1): It directly follows that the ij entry of the block A (p,q) takes the form . Complications arise when we seek to form a diagonal block A (p,p) , or even a block that is adjacent to a diagonal block. The difficulty is that the kernel K(t, s) has a singularity as s → t. To be precise, for any fixed t, there exist smooth functions u t and v t such that We see that when t (p) i is a point in Ω q the integrand in (3.9) becomes singular. When t (p) i is a point in a panel neighboring Ω q , the problem is less severe, but Gaussian quadrature would still be inaccurate. To maintain full accuracy, we use the modified quadrature rules described in Section 3.3.2. For every node t In order to have a quadrature evaluated at the Gaussian nodes t (q) j ∈ Ω q , we next use Inserting (3.12) into (3.11), we find that We now find that the block A (p,q) of A has entries We observe that the formula (3.13) is quite expensive to evaluate; in addition to the summation, it requires the construction of a quadrature rule for each point t (p) i and evaluation of Lagrange interpolants. Fortunately, this process must be executed for at most three blocks in each row of blocks of A.

A general algorithm
4.1. Summary. At this point, we have shown how to convert a BIE defined on an axisymmetric surface in R 3 to a sequence of equations defined on a curve in R 2 (Section 2), and then how to discretize each of these reduced equations (Section 3). Putting everything together, we obtain the following algorithm for solving (1.1): (1) Given the right hand side f , and a computational tolerance ε, determine a truncation parameter N F such that (2.15) holds.
(2) Form for n = −N F , −N F + 1, −N F + 2, . . . , N F the matrix A n discretizing the equation (2.13) encapsulating the n'th Fourier mode. The matrix is formed via Nyström discretization as described in Section 3 with the discretization parameters N P and N G chosen to meet the computational tolerance ε. The construction of the matrices A n in Step 2 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, and given discretization parameters N F , N P , and N G . The method therefore becomes particularly efficient when (1.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 of each matrix I + A n .

4.2.
Techniques for forming the matrices. We need to construct for each Fourier mode n, a matrix A n consisting of N P × N P blocks A (p,q) n , each of size N G × N G . Constructing an off-diagonal block A (p,q) n when Ω p and Ω q are not directly adjacent is relatively straightforward. For any pair of nodes t (p) i ∈ Ω p and t (q) j ∈ Ω q , we need to construct the numbers, cf. (3.3) and (3.10), j )/ds|, for n = −N F , −N F + 1, . . . , N F , where τ is a parameterization of γ (see Section 3.1) and the kernel k n is defined by (2.5). Fortunately, we do not need to explicitly evaluate the integrals in (2.5) since all the 2N F + 1 numbers can be evaluated by a single application of the FFT to the function j ), the function in (4.2) is smooth, and the trapezoidal rule implicit in applying the FFT is highly accurate.
Evaluating the blocks on the diagonal, or directly adjacent to the diagonal is somewhat more involved. The matrix entries are now given by the formula, cf.
where τ and k n are as in (4.1). To further complicate things, the points τ (t ) are now in close proximity to each other, and so the functions have a sharp peak around the point θ = 0. They are typically still easy to integrate away from the origin, so the integrals in (2.5) can for a general kernel be evaluated relatively efficiently using quadratures that are adaptively refined near the origin. Even with the accelerations described in this section, the cost of forming the matrices A n tends to dominate the computation whenever the kernels k n must be evaluated via formula (2.5). In particular environments, it is possible to side-step this problem by evaluating the integral in (2.5) analytically. That this can be done for the single and double layer kernels associated with Laplace's equation is demonstrated in Section 5.

Computational costs.
The asymptotic cost of the algorithm described in Section 4.1 has three components: (a) the cost of forming the matrices {A n } N F n=−N F , (b) the cost of transforming functions from physical space to Fourier space and back, and (c) the cost of solving the linear systems (I + A n ) σ n = f n . In this section, we investigate the asymptotic cost of these steps. We consider a situation where N F Fourier modes need to be resolved, and where N P × N G nodes are used to discretize the curve γ.
(a) Cost of forming the linear systems: Suppose first that we have an analytic formula for each kernel k n . (As we do, e.g., when the original BIE (1.1) involves either the single or the double layer kernel associated with Laplace's equation, see Section 5.) Then the cost T mat of forming the matrices satisfies When the kernels have to be evaluated numerically via formula (2.5), the cost of forming the matrices is still moderate. In the rare situations where the kernel is smooth, standard Gaussian quadrature can be used everywhere and the FFT acceleration described in Section 4.2 can be used for all entries. In this situation, In the more typical situation where each kernel k n involves an integrable singularity at the diagonal, the FFT acceleration can still be used to rapidly evaluate all entries well-removed from the diagonal. However, entries close to the diagonal must be formed via the composite quadrature rule combined with numerical evaluation of k n via an adaptive quadrature. In this situation, (b) Cost of Fourier transforms: The boundary data defined on the surface must be converted into the Fourier domain. This is executed via the FFT at a cost T fft satisfying We observe that the constant of proportionality in (4.5) is very small, and the cost of this step is typically negligible compared to the costs of the other steps.
(c) Cost of linear solves: Using standard Gaussian elimination, the cost T solve of solving N F linear systems (I + A n ) σ n = f n , each of size N P N G × N P N G , satisfies In situations where the equations need to be solved for multiple right hand sides, it pays off to first compute the inverses (I + A n ) −1 , and then simply apply these to each right hand side (or, alternatively, to form the LU factorizations, and then perform triangular solves). The cost T inv of computing the inverses, and the cost T apply of applying them then satisfy We make some practical observations: • The cost of forming the matrices by far dominates the other costs unless the kernel is either smooth, or analytic formulas for k n are available.
• The scheme is highly efficient in situations where the same equation needs to be solved for a sequence of different right hand sides. Given an additional right hand side, the added cost T solve is given by with a very small constant of proportionality. We note that this cost remains small even if an analytic formula for k n is not available.
• The system matrices I + A n often have internal structure that allow them to be inverted using "fast methods" such as, e.g., those in [17]. The cost of inversion and application can then be accelerated to near optimal complexity.

5.
Simplifications for the double layer kernels associated with Laplace's equation The solutions to (5.1) and (5.2) can be written in the respective forms where σ is a boundary charge distribution that can be determined using the boundary conditions. The resulting equations are 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 paper is equally applicable to single-layer formulations that lead to first kind Fredholm BIEs.
Equivalently, (5.5) and (5.6) can be arrived at by considering Laplace's equation written in cylindrical coordinates, Taking the Fourier transform of u with respect to theta θ gives n for n = −N F , −N F +1, . . . , N F need to be computed efficiently and with high accuracy to construct the Nyström discretization of (5.5) and (5.6). Note that the integrands of d n can be written in a similar form. This integrand is oscillatory and increasingly peaked at the origin as both r → r and z → z. As long as r and r as well as z and z 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 d n 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 n required and becomes prohibitively expensive if N F is large.
Fortunately, an analytical solution to (5.7) exists. As noted in [6], the single-layer kernel can be expanded with respect to the azimuthal variable as where Q n−1/2 is the half-integer degree Legendre function of the second kind and In light of this expansion, single-layer kernel can similarly be written as To find an analytical form for (5.7), first note that in cylindrical coordinates the doublelayer 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 d (i) n , which can be written using the previous equation as 8π 3 rr n r ∂Q n−1/2 (χ) ∂χ ∂χ ∂r − Q n−1/2 (χ) 2r + n z ∂Q n−1/2 (χ) ∂χ ∂χ ∂z .
To utilize this form of d (i) n , set µ = 2 χ+1 and note that where K and E 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 d Remark 5.2. Note that the forward recursion relation for the Legendre functions Q n−1/2 (χ) is unstable when χ > 1. In practice, the instability is mild when χ is near 1 and can still be employed to accurately compute values in this regime. Additionally, if stability becomes and issue, Miller's algorithm [8] can be used to calculate the values of the Legendre functions using the backwards recursion relation, which is stable for χ > 1.

Numerical results
This section describes several numerical experiments performed to assess the efficiency and accuracy of the the numerical scheme outlined in Section 4.1. All experiments were executed for the double layer kernels associated with Laplace's equation, calculated using the recursion relation described in Section 5.3. Note that the kernels in this case give us the property that A −n = A n , and so we need only to invert N F + 1 matrices. Further, the FFT used here is complex-valued, and a real-valued FFT would yield a significant decrease in computation time. The geometries investigated are described in Figure 6.1. The generating curves were parameterized by arc length, and split into N P equisized panels. A 10-point Gaussian quadrature has been used along each panel, with the modified quadratures used to handle the integrable singularities in the kernel. These quadratures are listed in the appendix. The algorithm was implemented in MatLab and the experiments were run on a 2.66GHz Intel Quad Core with 6Gb of RAM. All timings were averaged over 10 runs. N P the number of panels used to discretize the contour N F the Fourier truncation parameter (we keep 2N F + 1 modes) T setup time to setup the discrete system, excluding construction of the linear systems T mat time to construct the linear systems (utilizing the recursion relation) T inv time to invert the linear systems T fft time to Fourier transform the right hand side and the solution T apply time to apply the inverse to the right hand side The most expensive component of the calculation is the construction of the linear systems. This is primarily a result of the cost of evaluating the kernel and applying the modified quadrature rules. Table 2 compares the use of the recursion relation in evaluating the kernel when it is near-singular to using an adaptive Gaussian quadrature. The efficiency of the recursion relation is clearly evident in this case. Figure 6.2 plots the time to construct the linear systems as the number of panels and as the number of Fourier modes increases. The time to construct the systems as the number of panels increases grows as O(N 2 P N 2 G ) as expected, and the timings go as O(N F ) as the number of Fourier modes increase. Note that it is rather inexpensive to increase the number of Fourier modes used, making the scheme particularly attractive in environments where a large number is required. As an artifact of the small lengths of the vectors being Fourier transformed, the asymptotic behavior of the timings as the number of Fourier modes increases is not evident until we reach a large number of modes.
We observe that the largest problem reported in Table 1 involves 320 000 degrees of freedom. The method requires 2.2 minutes of pre-computation for this example, and is then capable of computing a solution u from a given data function f in 0.46 seconds.  Table 1. Timing results in seconds performed for the domain given in Figure 6.1(a) for the interior Dirichlet problem.
6.2. Accuracy and conditioning of discretization. The accuracy of the discretization has been tested using the interior and exterior Dirichlet problems on the domains given in Figure 6.1. Exact solutions were generated by placing a few random point charges outside of the domain where the solution was calculated. The solution was evaluated at points defined on a sphere encompassing (or interior to) the boundary. The errors reported in Tables 3, 4, and 5 are relative errors measured in the l ∞ -norm, ||u − u|| ∞ /||u|| ∞ , where u is the exact potential and u is the potential obtained from the numerical solution.
In all cases, 10 digits of accuracy has been obtained from a discretization involving a relatively small number of panels, due to the rapid convergence of the Gaussian quadrature. This is especially advantageous, as the most expensive component of the algorithm is the construction of the linear systems, the majority of the cost being directly related to the  Table 2. Timing comparison in seconds for constructing the matrices (I + A n ) using composite Gaussian quadrature and the recursion relation described in Section 5.3 to evaluate k n for diagonal and near diagonal blocks. The FFT is used to evaluate k n at all other entries. 2N F + 1 is the total number of Fourier modes used. 5 panels were used to discretize the boundary.
(a) (b) number of panels used. Further, the number of Fourier modes required to obtain 10 digits of accuracy is on the order of 100 modes. Although not investigated here, the discretization technique naturally lends itself to nonuniform refinement of the surface, allowing one to resolve features of the surface that require finer resolution. The number of correct digits obtained as the number of panels and number of Fourier modes increases eventually stalls. This is a result of a loss of precision in determining the kernels, as well as cancelation errors incurred when evaluating interactions between nearby points. This is especially prominent with the use of Gaussian quadratures, as points cluster near the ends of the panels. If more digits are required, high precision arithmetic can be employed in the setup phase of the algorithm.  Figure 6.1(a). The integral equations of this paper are second kind Fredholm equations, and generally lead to well-conditioned systems. As seen in Figure 6.3, this hold true for the discretization presented in this paper. N P 2N F + 1 -25 50 100 200 400 5 1.93869e-04 4.10935e-07 5.37883e-08 5.37880e-08 5.37880e-08 10 1.93869e-04 4.10513e-07 3.27169e-12 6.72270e-13 6.72270e-13 20 1.93869e-04 4.10513e-07 3.30601e-12 1.66132e-13 1.66132e-13 40 1.93869e-04 4.10513e-07 3.23162e-12 8.28568e-14 8.28568e-14 80 1.93869e-04 4.10512e-07 2.92918e-12 2.92091e-13 2.92091e-13 Table 3. Error in internal Dirichlet problem solved on domain (a) in Figure 6.1.

Generalizations and Conclusions
This paper describes a numerical technique for computing solutions to boundary integral equations defined on axisymmetric surfaces in R 3 with no assumption on the loads being axisymmetric. The technique is introduced as a generic method with only very mild conditions imposed on the kernel; specifically, we assume that the kernel has an integrable singularity at the diagonal, and that it is rotationally symmetric (in the sense that (2.2) holds). The technique described improves upon previous work in two regards: (1) A highly accurate quadrature scheme for kernels with integrable singularities is introduced. Numerical experiments indicate that solutions with a relative accuracy of 10 −10 or better can easily be constructed. (2) A rapid technique for numerically constructing the kernel functions k n in (1.3) is introduced. It works when k is either the single or the double layer potential associated with Laplace's equation. The technique is a hybrid scheme that relies on the FFT when possible, and uses recursion relations for Legendre functions when not. The resulting scheme is fast enough that a problem involving 320 000 degrees of freedom can be solved in 2.2 minutes on a standard desktop PC. Once one problem has been solved, additional right hand sides can be processed in 0.46 seconds.