Slevinsky JCP 332 2017 290

Journal of Computational Physics 332 (2017) 290–315
Contents lists available at ScienceDirect
Journal of Computational Physics

www.elsevier.com/locate/jcp
A fast and well-conditioned spectral method for singular

integral equations
Richard Mikael Slevinsky a,∗ , Sheehan Olver b
a
Department of Mathematics, University of Manitoba, Winnipeg, Canada
b
School of Mathematics and Statistics, The University of Sydney, Sydney, Australia
a r t i c l e i n f o a b s t r a c t
Article history: We develop a spectral method for solving univariate singular integral equations over
Received 3 August 2015 unions of intervals by utilizing Chebyshev and ultraspherical polynomials to reformulate
Received in revised form 4 October 2016 the equations as almost-banded infinite-dimensional systems. This is accomplished by
Accepted 3 December 2016
utilizing low rank approximations for sparse representations of the bivariate kernels. The
Available online 9 December 2016
resulting system can be solved in O (m2 n) operations using an adaptive QR factorization,
Keywords: where m is the bandwidth and n is the optimal number of unknowns needed to resolve
Spectral method the true solution. The complexity is reduced to O (mn) operations by pre-caching the QR
Ultraspherical polynomials factorization when the same operator is used for multiple right-hand sides. Stability is
Singular integral equations proved by showing that the resulting linear operator can be diagonally preconditioned to
be a compact perturbation of the identity. Applications considered include the Faraday
cage, and acoustic scattering for the Helmholtz and gravity Helmholtz equations, including
spectrally accurate numerical evaluation of the far- and near-field solution. The Julia
software package SingularIntegralEquations.jl implements our method with a
convenient, user-friendly interface.
© 2016 Elsevier Inc. All rights reserved.
1. Introduction
Singular integral equations are prevalent in the study of fracture mechanics [1], acoustic scattering problems [2–6], Stokes
flow [7], Riemann–Hilbert problems [8], and beam physics [9,10]. We develop a fast and stable algorithm for the solution of
univariate singular integral equations of general form [11]

× K (x, y )u ( y ) d y = f (x), for x ∈ , B u = c, (1)

where K (x, y ) is singular along the line y = x, the × in the integral sign denotes either the Cauchy principal value or the
Hadamard finite-part, is a union of bounded smooth open arcs in R2 , and B is a list of functionals. To be precise, we
consider the prototypical singular integral equations on [−1, 1] given by:
* Corresponding author.
E-mail addresses: Richard.Slevinsky@umanitoba.ca (R.M. Slevinsky), Sheehan.Olver@sydney.edu.au (S. Olver).
http://dx.doi.org/10.1016/j.jcp.2016.12.009
0021-9991/© 2016 Elsevier Inc. All rights reserved.
R.M. Slevinsky, S. Olver / Journal of Computational Physics 332 (2017) 290–315 291
1
1 K 1 (x, y ) K 2 (x, y )
= + + log | y − x| K 3 ( x, y ) + K 4 ( x, y ) u ( y ) d y = f (x), for x ∈ [−1, 1],
π ( y − x)2 y−x
−1
where K 1 , . . . , K 4 are continuous bivariate kernels.

In this work, we use several remarkable properties of Chebyshev polynomials including their spectral convergence, ex-
plicit formulæ for their Hilbert and Cauchy transforms, and low rank bivariate approximations to construct a fast and
well-conditioned spectral method for solving univariate singular integral equations. Chebyshev and ultraspherical polynomi-
als are utilized to convert singular integral operators into numerically banded infinite-dimensional operators. To represent
bivariate kernels, we use the low rank approximations of [12], where expansions in Chebyshev polynomials are constructed
via sums of outer products of univariate Chebyshev expansions. The minimal solution to the recurrence relation is auto-
matically revealed by the adaptive QR factorization of [13]. Diagonal right preconditioners are derived for integral equations
encoding Dirichlet and Neumann boundary conditions such that the preconditioned operators are compact perturbations of
the identity.
The inspiration behind the proposed numerical method is the ultraspherical spectral method for solving ordinary dif-
ferential equations [13], where ordinary differential equations are converted to infinite-dimensional almost banded linear
systems (an almost banded operator is a banded operator apart from a finite number of dense rows). These systems
can be solved in infinite-dimensions, i.e., without truncating the operators [14], as implemented in ApproxFun.jl
[15] in the Julia programming language [16,17]. The Julia software package SingularIntegralEquations.jl [18]
implements our method with a convenient, user-friendly interface. As an extension of this framework for infinite-
dimensional linear algebra, mixed equations involving derivatives and singular integral operators can be solved in a unified
way.
Several classical numerical methods exist for singular Fredholm integral equations of the first kind. These include: the
Nyström method [19–21], whereby integral operators are approximated by quadrature rules; the collocation method [22,23],
where approximate solutions in a finite-dimensional subspace are required to satisfy the integral equation at a finite number
of collocation points; and the Galerkin method [24,25], where the approximate solution is sought from an orthogonal
subspace and is minimal in the energy norm. The use of hybrid Gauss-trapezoidal quadrature rules [26–29] can significantly
increase the convergence rates when treating weakly singular kernels.
Numerous methods have exploited the underlying structure of the linear systems arising from discretizing inte-
gral equations. The most celebrated of these is the Fast Multipole Method of Greengard and Rokhlin [30]. Other
characterizations in terms of semi-separability or other hierarchies have also gained prominence [31–33]. Exploiting
the matrix structure allows for fast matrix–vector products, which then allows for Krylov subspace methods [34]
to be extremely competitive. For scattering of the Helmholtz equation in very special geometries, hybrid numerical-
asymptotic methods have been derived for frequency-independent solutions to the Dirichlet and Neumann problems [4,
35,6,36].
Previous works on Chebyshev-based methods for singular integral equations include Frenkel [37], which derives recur-
rence relations for the Chebyshev expansion of a singular integral equation after expanding the bivariate kernel in a basis
of Chebyshev polynomials of the first kind in both variables, and Chan et al. [38,39] in fracture mechanics, among others.
A similar analysis in [40] is used for hypersingular integrodifferential equations by expanding the bivariate kernel in a basis
of Chebyshev polynomials of the second kind. This paper is an extension of these ideas with essential practical numerical
considerations.
Remarks.
1. Combined with fast multiplication of Chebyshev series, our method is suitable for use in iterative Krylov subspace
methods.
2. There is a great diversity of integral equation formulations. The choice of formulation depends on many properties,
including for example, whether the boundary is open or closed and whether there are resonances. Most equations
involve operators that contain manipulations of the fundamental solution, which would still satisfy the requirements of
our method. However, we focus on the direct integral equations to retain a simple exposition.
2. Boundary integral equations in two dimensions
In two dimensions, let x = (x1 , x2 ) and y = ( y 1 , y 2 ). Positive definite second-order linear elliptic partial differential oper-
ators (PDOs) with variable coefficients are always reducible to the following canonical form [41]:
∂u ∂u
L{u } = u + a +b + cu . (2)
∂ x1 ∂ x2
Let (x, y) denote the positive definite fundamental solution of (2) satisfying the formal partial differential equation (PDE)
292 R.M. Slevinsky, S. Olver / Journal of Computational Physics 332 (2017) 290–315
Lx {} = −δ(x − y), (3)
where δ is the two-dimensional Dirac delta distribution and the subscript indicates that L is acting in the x variable.
2.1. Exterior scattering problems
Let be a union of disjoint bounded smooth open arcs in R2 and let := R2 \ .

1 −ix·y u (y) dy be the standard Fourier transform in R2 . Then for s ∈ R, H s (R2 ) defines the Bessel
Let û (x) := 2 e
2π R
potential space as the set of normed tempered distributions equipped with:

u 2H s (R2 ) = (1 + |x|2 )s |û (x)|2 dx. (4)
R2

For bounded with non-empty interior, the space H s () := u | : u ∈ H s (R2 ) . Furthermore, let H loc
s
() denote the spaced

of locally integrable functions in H s () and lastly, let H̃ s () := u ∈ H s (R2 ) : supp u ⊂ .
Definition 1. For x ∈ , let S and D define the single- and double-layer potentials:

1
S u (x) = (x, y)u (y) d(y) : H̃ − 2 () → H loc
1
(), (5)

∂(x, y) 1
1
D u (x) = u (y) d(y) : H̃ 2 () → H loc (). (6)
∂ n(y)

Definition 2 (Radiation condition at infinity [42]). We say that u ∈ H loc

1
() satisfies the radiation condition at infinity if:

lim S|y|=ρ [∂ u /∂ n] (x) − D|y|=ρ [u] (x) = 0, for x ∈ . (7)
ρ →+∞
For solutions of the homogeneous equation L{u } = 0 satisfying the radiation condition at infinity, Green’s representation
theorem allows for the determination of the exterior solutions given data on the boundary :
u (x) = −S [∂ u /∂ n] (x) + D [u] (x), for x ∈ . (8)
Here, [u ] denotes the jump in u along and [∂ u /∂ n] the jump in its normal derivative. These are formally defined by
the Dirichlet trace and conormal derivative [25], or in the case of the Laplace equation, simply as the difference between
the limiting values on as we approach from the left and the right. This identity can be interpreted as representing u
in terms of the potential of a distribution of poles on through the single-layer and normal dipoles on through the
double layer. With either Dirichlet or Neumann boundary conditions, we restrict (8) to the boundary and solve for the
unknown boundary value. Once both quantities on the boundary are determined, the solution to the exterior problem is
readily available in integral form.
1
Dirichlet Problem. Given an incident wave u i (x) ∈ H 2 () satisfying L{u i } = 0, find u s (x) ∈ H loc
1
() satisfying L{u s } = 0,
the radiation condition at infinity, and
u i (x) + u s (x) = 0, for x ∈ . (9)
∂ u i (x) 1
Neumann Problem. Given ∈ H − 2 () satisfying L{u i } = 0, find u s (x) ∈ H loc
1
() satisfying L{u s } = 0, the radiation
∂ n(x)
condition at infinity, and
∂ i
u (x) + u s (x) = 0, for x ∈ . (10)
∂ n(x)
For the case of the Laplace and Helmholtz equations, the Dirichlet problem is originally formulated in [43, Eqs. (1.1) &
(1.2), (1.6) & (1.7)]. Similarly, the Neumann problem is originally formulated in [44, Eqs. (1.1) & (1.2)].
Dirichlet Solution. The Dirichlet problem is solved by (8) where [u s ] = 0, and the scattered solution is represented
1
everywhere by the single-layer potential. The density [∂ u s /∂ n] ∈ H̃ − 2 () in (8) satisfies:

∂ us
(x, y) d(y) = u i (x), for x ∈ . (11)
∂n

Neumann Solution. The Neumann problem is solved by (8) where [∂ u s /∂ n] = 0, and the scattered solution is represented
1
everywhere by the double-layer potential. The density [u s ] ∈ H̃ 2 () in (8) satisfies:

∂ ∂(x, y) s ∂ u i (x)
u d(y) = − , for x ∈ . (12)
∂ n(x) ∂ n(y) ∂ n(x)

For the case of the Laplace and Helmholtz equations, the Dirichlet solution is originally proved in [43, Theorems 1.4
& 1.7]. Similarly, the Neumann solution is originally proved in [44, Theorem 1.3]. Furthermore, by appealing to the theory
of Mellin transforms, inverse square root singular behaviour is derived for the open ends of in the Dirichlet problem [43,
Theorem 2.3], and square root singular behaviour is derived for the open ends of in the Neumann problem [44, Theo-
rem 1.8].
For elliptic PDOs with variable coefficients, we use the solutions provided by the singular integral equations (11) and (12),
and while the scope of this paper is numerical, we conjecture that they hold more generally.
2.2. Riemann functions
In addition to the PDO in (2), consider its adjoint:
∂(av ) ∂(bv )
L∗ { v } = v − − + cv . (13)
∂ x1 ∂ x2
With the change to complex characteristic variables:
z = x1 + ix2 , ζ = x1 − ix2 , z0 = y 1 + i y 2 , ζ0 = y 1 − i y 2 , (14)
L and L∗ take the form:
∂2U ∂U ∂U
L̂{U } = +A +B + CU, (15)
∂ z∂ζ ∂z ∂ζ
∂2V ∂( A V ) ∂( B V )
L̂∗ { V } = − − + CV, (16)
∂ z∂ζ ∂z ∂ζ
where:

1 z+ζ z−ζ z+ζ z−ζ
A ( z, ζ ) = a , + ib , , (17)
4 2 2i 2 2i

1 z+ζ z−ζ z+ζ z−ζ
B ( z, ζ ) = a , − ib , , (18)
4 2 2i 2 2i

1 z+ζ z−ζ
C ( z, ζ ) = c , . (19)
4 2 2i
Theorem 3 (Vekua and Garabedian [41,45]). For analytic functions (17)–(19), there exist analytic functions R( z, ζ, z0 , ζ0 ) and
g 0 ( z, ζ, z0 , ζ0 ) such that:
1
( z, ζ, z0 , ζ0 ) = − R( z, ζ, z0 , ζ0 ) log[( z − z0 )(ζ − ζ0 )] + g 0 ( z, ζ, z0 , ζ0 ), (20)
4π
where L̂{} = 0 in ( z, ζ ) and L̂∗ {} = 0 in ( z0 , ζ0 ) so long as z = z0 and ζ = ζ0 . In (20), R is the Riemann function of the operator L
satisfying:
L̂∗ {R} = 0, (21)

⎧ ⎫
⎪
⎨ζ ⎪
⎬
R( z0 , ζ, z0 , ζ0 ) = exp A ( z 0 , τ ) dτ , and (22)
⎪
⎩ ⎪
⎭
ζ0
⎧ z ⎫
⎨ ⎬
R( z, ζ0 , z0 , ζ0 ) = exp B (t , ζ0 ) dt . (23)
⎩ ⎭
z0
Remarks. It is straightforward to reformulate (21)–(23) to the following integral equation:
z ζ
R( z, ζ, z0 , ζ0 ) − B (t , ζ )R(t , ζ , z0 , ζ0 ) dt − A ( z, τ )R( z, τ , z0 , ζ0 ) dτ
z0 ζ0
z ζ
+ C (t , τ )R(t , τ , z0 , ζ0 ) dτ dt = 1. (24)
z0 ζ0
Returning to the original coordinates x and y, fundamental solutions for elliptic PDOs with analytic coefficients can be
written as:
(x, y) = A (x, y) log |x − y| + B (x, y), (25)

where A and B are both analytic functions of x and y and where A (x, y) = R( z, ζ, z0 , ζ0 ) implying A (x, x) = −(2π )−1 .
− 21π
If, furthermore, the PDO is formally self-adjoint, then A and B are also symmetric functions of x and y.
3. Practical approximation theory
Chebyshev approximation theory is a very rich subject that has seen numerous exceptional contributions: see [46–48]
and the references therein. In this section, we describe some approximation spaces for one-dimensional intervals and two-
dimensional squares. For every approximation space, one may consider the interpolants, which are equal to the function at
a set of interpolation points, and the projections, which are truncations of the function’s expansion. Unless an extraordinary
amount of analytic information is known about a function, interpolants are generally easier to construct.
We consider an approximation space practical if there is a fast way to transform the interpolation condition into approx-
imate projections. While a few methods exist to create fast transforms, all the practical approximation spaces we consider
resort to some variation of the fast Fourier transform (FFT) [49,50] to reduce O (n2 ) complexity to O (n log n). Other prop-
erties which make an approximation space practical are: O (n) evaluation; a low Lebesgue constant; absolute, uniform, and
geometric convergence with analyticity; and, easy manipulation for the development of new properties. For approximation
on the canonical unit interval I := [−1, 1], we will make our statements precise in the following subsection.
3.1. One dimension
Let K be the field of R or C. A function f : I → K is of bounded total variation if:

Vf = | f ( z)| dz < +∞. (26)
I
Chebyshev polynomials of the first kind are defined by [47]:
T n (x) = cos(n cos−1 (x)), for n ∈ N0 , and x ∈ I. (27)

A Chebyshev interpolant to a continuous function f : I → K is the approximation

N −1
p N (x) = cn T n (x), x ∈ I, (28)
n =0
which interpolates f at the Chebyshev points of the first kind:

2n + 1
p N (xn ) = f (xn ) where xn = cos π , for n = 0, . . . , N − 1. (29)
2N
The Chebyshev basis has fast transforms between values at Chebyshev points and coefficients via fast implementations
of the discrete cosine transforms (DCTs). The (orthogonal) Chebyshev polynomials satisfy a three-term recurrence relation
that can be used in Clenshaw’s algorithm [51] for O (n) evaluation of interpolants. Compared with the best polynomial
approximants, Chebyshev interpolants are near-best in the sense that their Lebesgue constants exhibit similar logarithmic
growth.
Theorem 4 (Battles and Trefethen [52]). Let f be a continuous function on I, p N its N-point polynomial interpolant in the Chebyshev
points of the first kind and p N its best degree-N − 1 polynomial approximation. Then:

1. f − p N ∞ ≤ 2 + π2 log N − 1 f − p N ∞ ;
2. if f has a kth derivative in I of bounded variation for some k ≥ 1, f − p N ∞ = O ( N −k ) as N → ∞; and,

3. if f is analytic in a neighbourhood of I, f − p N ∞ = O (C N ) as N → ∞ for some C < 1; in particular we may take C =
1/( M + m) if f is analytic in the closed Bernstein ellipse with foci ±1 and semimajor and semiminor axis lengths M ≥ 1 and
m ≥ 0.
An interpolant can be constructed to any relative or absolute tolerance by successively doubling the number of inter-
polation conditions, transforming values to coefficients, and determining an acceptable degree.1
3.2. Two dimensions
Numerous methods have been devised to approximate functions in more than one dimension. The straightforward gen-
eralization of the one-dimensional approach is to sample the function on a tensor of one-dimensional interpolation points
and to adaptively truncate coefficients below a certain threshold.
Consider the function f : I2 → K , whose two-dimensional Chebyshev interpolant takes the form:

m −1 n
−1
pm,n (x, y ) = A i , j T i (x) T j ( y ). (30)
i =0 j =0
While the tensor approach in general suffers from the curse of dimensionality, it can still be competitive in two dimensions,
scaling with O (mn) function samples and O (min(mn log n, nm log m)) arithmetic via fast two-dimensional transforms.
The singular value decomposition of an m × n matrix A over K is the factorization [53]:
A = U V∗ , (31)
where U is an m × m unitary matrix over K , is an m × n diagonal matrix of non-negative singular values, and V∗ is an n × n
unitary matrix over K . The singular value decomposition reveals the rank of a matrix as the number of nonzero singular
values.
If we perform the singular value decomposition of the matrix of coefficients in (30), the approximation to f can be
re-expressed as:

k
p SVD (x, y ) = σi u i (x) v ∗i ( y ), (32)
i =1
where σi are the singular values, and u i (x) and v ∗i ( y ) are univariate Chebyshev approximants with coefficients from the
columns of U and the rows of V∗ , respectively, and where A is of rank k. It follows that p SVD is the best rank-k approximant
in L 2 (I2 ) to f that can be obtained for the original two-dimensional interpolant. For any given tolerance > 0, a function
f has numerical rank k if [54]

k = inf inf f − f k ∞ ≤ f ∞ , (33)
k∈N fk
where the inner infimum is taken over all rank-k functions.
Definition 5 (Townsend Definition 3.1 [54]). For some > 0, let k be the numerical rank of f : I2 → K , and m and n be the
maximal degrees of the univariate approximations in the x and y variables. If k (m + n ) < m n , we say the function f is
numerically of low rank, and if k ≈ min(m , n ), then the function f is numerically of full rank.
A particularly attractive scheme for calculating low rank approximation in two dimensions can be described as a contin-
uous analogue of Gaussian elimination [54] and is a direct extension of the greedy algorithm in one dimension [48, Chapter
5]. This algorithm is studied in depth in Townsend’s D.Phil. thesis and implementations are found in Chebfun [55] and
ApproxFun.jl [15]. In this algorithm, the function is initially sampled on a grid to locate its approximate absolute maxi-
mum. Two one-dimensional approximations are created in the x and y variables to interpolate the function along the row
and column that intersect at the approximate absolute maximum. After subtracting this rank-one approximation, the algo-
rithm continues its search for the next approximate absolute maximum. After k iterations, it is clear that the approximant

k
p GE (x, y ) = A i (x) B i ( y ), (34)
i =1
1
This heuristic determination is usually based on, among other things, the relative and absolute magnitudes of initial and final coefficients, the decay
rate of the coefficients, an estimate of the condition number of the function, and an estimate of the Lebesgue constant for a given degree.
coincides with f in the k rows and columns whose intersections coincide with an iteration’s approximate absolute maxi-
mum. As the size of the sampling grid increases, the approximate absolute maxima will converge to the true absolute max-
ima and in this sense we reproduce close approximations to p SVD . In terms of the degrees of the one-dimensional approx-
imations m, n and the rank k, the algorithm scales with a search over O (mn) function samples and O (k (m log m + n log n))
arithmetic via fast one-dimensional transforms.
Definition 6 (Townsend Definition 4.11 [54]). The function f : I2 → K is Hermitian if it satisfies the conjugate symmetry
f (x, y ) = f ∗ ( y , x) and it is non-negative definite, i.e.:

a∗ ( y ) f ( y , x)a(x) d y dx ≥ 0, (35)
I2
for all a(x) ∈ C (I).
When a bivariate function is Hermitian, even further savings can be obtained by drawing the analogy to the Cholesky
factorization of a Hermitian matrix [56]:

k
p Cholesky (x, y ) = A i (x) A ∗i ( y ). (36)
i =1
In this case, it is known that the function’s absolute maxima after every iteration are on the diagonal line y = x, leading to
a reduction in the dimension of the search space. In addition, as they are conjugates only either the row or column slices
may be computed and stored.
3.3. An algorithm to extract the splitting of a fundamental solution
Accurate numerical evaluation of a fundamental solution on or near the singular diagonal may not always be possible
or may be more expensive [57]. To avoid the numerical problems associated with the singular diagonal, we use Chebyshev
points of the first kind in one direction and Chebyshev points of the second kind [47] in the other direction. This ensures
that the diagonal is never sampled. In terms of the DCTs, taking 2n points of the first kind is optimal and taking 2n + 1
points of the second kind is nearly optimal.
When both A (x, y) and B (x, y) in (25) are not known a priori, but the fundamental solution itself can be evaluated, we
can use such skewed grids in combination with the Riemann function R to:
1. approximate A (x, y) ≡ − 21π R(x1 + ix2 , x1 − ix2 , y 1 + i y 2 , y 1 − i y 2 ); and subsequently,

2. approximate the difference B (x, y) ≡ (x, y) − A (x, y) log |x − y|.
4. The ultraspherical spectral method
The ultraspherical spectral method of Olver and Townsend [13] represents solutions of linear ordinary differential equa-
tions of the form
Au = f , Bu = c, (37)
where A is a linear operator of the form
dN d
A = a N (x) + · · · + a1 (x) + a0 (x), (38)
dx N dx
and B contains N linear functionals. Typically, B encodes boundary conditions such as Dirichlet or Neumann conditions. We
consider u (x) in its Chebyshev expansion
∞

u (x) = un T n (x), (39)
n =0
so that u (x) can be identified by a vector of its Chebyshev coefficients u = (u 0 , u 1 , . . .) .

To solve such a problem efficiently, a change of basis occurs for each order of spectral differentiation, using the formula:

dλ T n (x) 0, 0 ≤ n ≤ λ − 1,
= (λ) (40)
dxλ 2λ−1 (λ − 1)! n C n−λ (x), n ≥ λ,
where C nλ represents the ultraspherical polynomial of integral order λ and of degree n. This sparse differentiation has the
operator representation:
⎛ ⎞
λ times

⎜0 ··· 0 λ ⎟
⎜ ⎟
⎜ λ+1 ⎟
Dλ = 2λ−1 (λ − 1)! ⎜ ⎟, λ ≥ 1, (41)
⎜ λ+2 ⎟
⎝ ⎠
..
.
and maps the Chebyshev coefficients to the λth order ultraspherical coefficients.
Since in (38), each derivative maps to a different ultraspherical basis, the sparse differentiation operators are accompa-
nied by sparse conversion operators such that A can be expressed completely in the basis of highest order N:
⎛ ⎞ ⎛ ⎞
1 0 − 12 1 0 − λ+
λ
2
⎜ 1
− 12 ⎟ ⎜ λ
0 − λ+
λ ⎟
⎜ 2
0 ⎟ ⎜ λ+1 3 ⎟
S0 = ⎜
⎜ .. ⎟ , Sλ = ⎜
⎜ .. ⎟ , λ ≥ 1. (42)
⎝
1
2
0 .⎟⎠ ⎝
λ
λ+2 0 .⎟⎠
.. .. .. ..
. . . .
Here, S0 maps the Chebyshev coefficients to the first order ultraspherical coefficients and Sλ maps the λth order ultras-
pherical coefficients to the (λ + 1)th order ultraspherical coefficients. Therefore, the conversion and differentiation operators
can be combined in A as follows:
(a N D N + a N −1 S N −1 D N −1 + · · · + a0 S N −1 · · · S0 ) u = S N −1 · · · S0 f, (43)
where u and f are vectors of Chebyshev expansion coefficients. Were the coefficients ai (x), i = 0, . . . , N, constant, then (43)
would represent a linear recurrence relation in the coefficients u of length at most 2N + 1. However, the coefficients are in
general not constants, so the multiplication operators in Chebyshev and ultraspherical bases are also investigated in [13].
Let
∞

a(x) = an T n (x). (44)
n =0
Then it is shown in [13] that multiplication can be represented as a Toeplitz-plus-Hankel-plus-rank-one operator:

⎡⎛
2a0 a1 a2 ···
⎞ ⎛ ⎞⎤
0 0 ···0
⎢⎜
1 ⎢⎜ a1 2a0 a1
.. ⎟ ⎜a
.⎟ ⎜ 1 a2 · · · ⎟⎥
a3 ⎥
M0 [a] = ⎢⎜ .. ⎟ +⎜ .⎟⎟⎥ . (45)
2⎢ ⎜
⎣⎝ a2 a1 2a0 .⎟⎠ ⎝ a2 a3 a4 . . ⎠⎥ ⎦
.. .. .. .. . . .
. . .
..
. . .. .. ..
For λ > 0, an explicit formula for the entries is given in [13] and a three-term recurrence relation is shown in [54, Chap.
6]. By the associative and distributive properties of multiplication, the recurrence relation for the multiplication operators is
derived from the recurrence relation for the ultraspherical polynomials:
(λ) 2(n + λ) (λ) n + 2λ − 1 (λ)

Mλ [C n+1 ] = Mλ [x]Mλ [C n ] − Mλ [C n−1 ], n ≥ 1. (46)
n+1 n+1
Since we assume the coefficients ai (x) to be continuous functions with bounded variation on I, let m denote the highest
degree Chebyshev expansion such that for some > 0:
( (
(
m −1 (
( (
(ai (x) − ain T n (x)( ≤ ai (x)∞ , for i = 0, . . . , N . (47)
( (
n =0 ∞
Then in this way, the system
Bu = c,
(M N [a N ]D N + M N [a N −1 ]S N −1 D N −1 + · · · + M N [a0 ]S N −1 · · · S0 ) u = S N −1 · · · S0 f, (48)
is almost banded with bandwidth O (m). The proposed O (m n) solution process for such systems is the adaptive QR factor-
2
ization, generalizing (F. W. J.) Olver’s algorithm for second-order difference equations [58]. In this factorization, the forward
error is estimated at every step in the infinite-dimensional upper-triangularization to adaptively determine the minimal
order n required to resolve the solution below a pre-determined accuracy. Since the unitary transformations implied by Q
preserve the rank structure, the back substitution is also performed with O (m2 n) complexity.
Fig. 1 shows the typical structure of the system and an example of the type of singularly perturbed boundary value
problem that it can solve efficiently.
Fig. 1. Solution of ( + x2 )u (x) = x u (x), u (−1) = 1, u (1) = 0 via the ultraspherical spectral method. Left: the structure of the system. Right: a plot of the
solution for = 10−4 . In this case, a Chebyshev expansion of degree 3,276 is required to approximate the solution to double precision.
4.1. Almost-banded spectral methods in other bases
The key elements of the ultraspherical spectral method are a graded set of bases that permit banded differentiation and
conversion within the set of bases, and multiplication operators for variable coefficients. Other examples where a graded
basis can be exploited are the Jacobi polynomials (which include Legendre and ultraspherical polynomials as special cases),
and the generalized Laguerre polynomials. Hermite polynomials, which form an Appell sequence, satisfy H n (x) = 2nH n−1 (x),
and therefore do not require other bases for conversion.
From the three-term recurrence relation satisfied by orthogonal polynomials [59]:
xπn (x) = αn πn+1 (x) + βn πn (x) + γn πn−1 (x), (49)

it follows that multiplication by x is tridiagonal:
⎛ ⎞
β0 α0
⎜ γ1 β1 α1 ⎟
M[x] = ⎜
⎝ γ2 β2 α2
⎟.
⎠ (50)
.. .. ..
. . .
Therefore, banded multiplication operators in orthogonal bases can be derived from the recurrence relation:

M[x] − βn γn
M[πn+1 ] = M[πn ] − M[πn−1 ], n ≥ 1, (51)
αn αn
and consequently variable coefficients represented as interpolants have a finite-bandwidth operator form. To numerically
determine such variable coefficients practically requires fast transforms. Among the many possibilities, see [60] for a new
approach for a fast FFT-based discrete Legendre transform.
5. Ultraspherical spectral method for singular integral equations
In the following definitions, we identify C with R2 and let be bounded in C.
Definition 7 (Kress [5]). A real- or complex-, scalar- or vector-valued function f defined on is called uniformly Hölder
continuous with Hölder exponent 0 < α ≤ 1 if there exists a constant C such that
| f (x) − f (y)| ≤ C |x − y|α , for x, y ∈ . (52)
By C 0,α () we denote the space of all bounded and uniformly Hölder continuous functions with exponent α . For vectors,
we take | · | to be the Euclidean distance. With the norm
| f (x) − f (y)|
f 0,α := sup | f (x)| + sup , (53)
x∈ x,y∈ |x − y|α
x=y
the Hölder space is a Banach space, and we can further introduce C 1,α () as the space of all differentiable functions whose
gradient belongs to C 0,α ().
Definition 8. Let f ∈ C 0,α (). The Cauchy transform over is defined as:

1 f (ζ )
C f ( z) := dζ, for z ∈ C \ . (54)
2π i ζ −z

The Cauchy transform can be extended to z ∈ with integration understood as the Cauchy principal value.
Definition 9. Let f ∈ C 0,α (). The Hilbert transform over is defined as:

1
f (ζ )
H f ( z) := − dζ, for z ∈ , (55)
π ζ −z

where the integral is understood as the Cauchy principal value:

1 f (ζ ) 1 f (ζ )
− dζ = lim dζ, (56)
π ζ −z π ρ →0 ζ −z
\( z;ρ )
where ( z; ρ ) := {ζ ∈ : |ζ − z| ≤ ρ }.
Lemma 10 (Sokhotski–Plemelj [61,62]). If f ∈ C 0,α (), then:
H f ( z) = i[C + + C − ] f ( z), (57)

where C ± denotes the limit from the left/right of .
With the Hilbert and Cauchy transforms, further integrals with singularities can be defined.
Definition 11. For f ∈ C 0 () the log transform over is defined as:

1
L f ( z) := log |ζ − z| f (ζ ) dζ, for z ∈ C. (58)
π

For f ∈ C 1, α
() the derivative of the Hilbert transform is defined as:

1 f (ζ )
H f ( z) := = dζ, for z ∈ , (59)
π (ζ − z)2

where the integral is understood as the Hadamard finite-part [63,64]:

⎧ ⎫
⎪
⎨ ⎪
1 f (ζ ) 1 f (ζ ) 2 f ( z) ⎬
= dζ = lim dζ − , (60)
π (ζ − z)2 π ρ →0 ⎪
⎩ (ζ − z)2 ρ ⎪
⎭
\( z;ρ )
where ( z; ρ ) is defined as in Definition 9.
Remarks.
1. The Sokhotski–Plemelj lemma offers a convenient way to compute the Hilbert transform via the limit of two Cauchy
transforms.
2. The use of the Cauchy principal value and the Hadamard finite-part allows for the regularization of singular and hyper-
singular integral operators, respectively.
On a contour , we expand the kernel of the singular integral equation (1) in the following way:
Au = f , Bu = c, (61)
for
1
1 K 1 (x, y ) K 2 (x, y )
Au = = + + log | y − x| K 3 ( x, y ) + K 4 ( x, y ) u( y) d y,
π ( y − x)2 y−x
−1
where K 1 , K 2 , K 3 and K 4 are known continuous bivariate kernels, f is continuous, B contains N linear functionals, and u
is the unknown solution. If in (61), we replace the bivariate kernels with low rank approximations,

kλ
K λ (x, y ) ≈ A λ,i (x) B λ,i ( y ), for λ = 1, 2, 3, 4, (62)
i =1
we achieve at once two remarkable things: firstly, the approximations are compressed representations of the kernels; and
secondly, the separation of variables in the low rank approximation allows for the singular integral operators to be con-
structed via the Definitions 9 and 11.
In the following two subsections, we consider the case where is the unit interval, and emulate the construction of the
ultraspherical spectral method for ODEs to arrive at an almost-banded system to represent (61). In this setting, we must
use weighted Chebyshev bases to accomplish this task. Note that alternative spectral methods for open arcs are discussed
in [21].
5.1. Inverse square root endpoint singularities
Indeed, the Hilbert transform of weighted Chebyshev polynomials is known [65]:

T n (x) 0, n = 0,
H(−1,1) √ = (1 ) (63)
1−x 2 C n−1 (x), n ≥ 1.
This operation can then be expressed as the banded operator from the weighted Chebyshev coefficients to the ultraspherical
coefficients of order 1:
⎛ ⎞
0 1
⎜ 1 ⎟
H(−1,1) = ⎜
⎝ 1
⎟.
⎠ (64)
..
.
Upon integration with respect to x, we obtain an expression for the log transform:

)
T n (x) − log 2, n = 0,
L(−1,1) √ = T n (x) (65)
1 − x2 − , n ≥ 1,
n
or as an operator from the weighted Chebyshev coefficients to the Chebyshev coefficients:
⎛ ⎞
− log 2
⎜ −1 ⎟
⎜ ⎟
L(−1,1) = ⎜ − 12 ⎟. (66)
⎝ ⎠
..
.
In addition, upon differentiation with respect to x, we also obtain an expression for the derivative of the Hilbert transform:

T n (x) 0, n = 0, 1 ,
H(− √ = (2 ) (67)
1,1)
1 − x2 C n−2 (x), n ≥ 2.
This operation can then be expressed as the banded operator from the weighted Chebyshev coefficients to the ultraspherical
coefficients of order 2:
⎛ ⎞
0 0 1
⎜ 1 ⎟
⎜ ⎟.
H(− 1,1) = ⎝ 1 ⎠ (68)
..
.
Lastly, the orthogonality of the Chebyshev polynomials immediately yields for the functional

1
f := f (ζ ) dζ (69)
π

the following:

T n (x) 1 , n = 0,
(−1,1) √ = (70)
1 − x2 0, n ≥ 1 ,
or as a compact functional on the weighted Chebyshev coefficients:

* +
(−1,1) = 1 0 0 · · · . (71)
Combining the integral operators together with the bivariate approximations, we define:

k1

H(− 1,1) [ K 1 ] := M2 [ A 1,i (x)]H(− 1,1) M0 [ B 1,i ( y )], (72)
i =1

k2
H(−1,1) [ K 2 ] := M1 [ A 2,i (x)]H(−1,1) M0 [ B 2,i ( y )], (73)
i =1

k3
L(−1,1) [ K 3 ] := M0 [ A 3,i (x)]L(−1,1) M0 [ B 3,i ( y )], (74)
i =1

k4
(−1,1) [ K 4 ] := M0 [ A 4,i (x)](−1,1) M0 [ B 4,i ( y )]. (75)
i =1
Then, we can reduce singular integral equations of the form (61) into an infinite-dimensional almost-banded system:
Bu = c,

H(− 1,1) [ K 1 ] + S1 H(−1,1) [ K 2 ] + S1 S0 (L(−1,1) [ K 3 ] + (−1,1) [ K 4 ]) u = S1 S0 f. (76)
This system can be solved directly using the framework of infinite-dimensional linear algebra [14], built out of the adaptive
QR factorization introduced in [13].
5.2. Square root endpoint singularities
The Hilbert transform of weighted Chebyshev polynomials of the second kind is also known [65]:
, - .
HI U n (x) 1 − x2 = − T n+1 (x), n ≥ 0. (77)
This operation can then be expressed as the banded operator from the weighted ultraspherical coefficients of order 1 to the
Chebyshev coefficients:
⎛ ⎞
0
⎜ −1 ⎟
HI = ⎜
⎝ −1
⎟.
⎠ (78)
..
.
Upon integration with respect to x, we obtain an expression for the log transform:
⎧
⎪ 1 1
, - . ⎪
⎨ − log 2 + T 2 (x), n = 0,
2 4
LI U n (x) 1 − x2 = (79)
⎪
⎪ 1 T n+2 (x) T n (x)
⎩ − , n ≥ 1,
2 n+2 n
or as an operator from the weighted ultraspherical coefficients of order 1 to the Chebyshev coefficients:
⎛ ⎞
− 12 log 2
⎜ 0 − 12 ⎟
⎜ ⎟
LI = ⎜ 1
0 − 14 ⎟. (80)
⎝ 4 ⎠
.. .. ..
. . .
In addition, upon differentiation with respect to x, we also obtain an expression for the derivative of the Hilbert transform:
, - .
HI U n (x) 1 − x2 = −(n + 1)C n (x),
(1 )
n ≥ 0. (81)
This operation can then be expressed as the banded operator from the weighted ultraspherical coefficients of order 1 to the
ultraspherical coefficients of order 1:
⎛ ⎞
−1
⎜ −2 ⎟
HI = ⎜
⎝ −3
⎟.
⎠ (82)
..
.
Lastly, the orthogonality of the Chebyshev polynomials of the second kind immediately yields for I :
)
, - . 1
, n = 0,
I U n (x) 1 − x2 = 2 (83)
0, n ≥ 1,
or as a compact functional on the weighted Chebyshev basis:

*1 +
I = 2
0 0 ··· . (84)
Combining the integral operators together with the bivariate approximations, we define:

k1
HI [ K 1 ] := M1 [ A 1,i (x)]HI M1 [ B 1,i ( y )], (85)
i =1

k2
HI [ K 2 ] := M0 [ A 2,i (x)]HI M1 [ B 2,i ( y )], (86)
i =1

k3
LI [ K 3 ] := M0 [ A 3,i (x)]LI M1 [ B 3,i ( y )], (87)
i =1

k4
I [ K 4 ] := M0 [ A 4,i (x)]I M1 [ B 4,i ( y )], (88)
i =1
and in the framework of infinite-dimensional linear algebra [14], we may solve singular integral equations of the form (61)
via the almost-banded system:
Bu = c,
* +
HI [ K 1 ] + S0 (HI [ K 2 ] + LI [ K 3 ] + I [ K 4 ]) u = S0 f. (89)
Let mx + m y denote the largest sum of degrees of the bivariate Chebyshev expansions of the integral kernels such that
for some > 0:
( (
( kλ (
( (
( K λ (x, y ) − A ( x) B ( y )( ≤ K λ (x, y )∞ , for λ = 1, 2, 3, 4. (90)
( λ,i λ,i (
( i =1 (
∞
Then, the complexity of the adaptive QR factorization is O ((mx + m y )2 n) operations, where n is degree of the resulting
weighted Chebyshev expansion of the solution. This is reduced to O ((mx + m y )n) operations by pre-caching the QR factor-
ization.
Remarks.
1. The observation that |dζ | = dζ on I allows us to relate line integral formulations with the operators of Definitions 9
and 11.2
2. Mixed equations involving derivatives and singular integral operators are also covered in this framework.
3. It is straightforward to obtain the singular integral operators on arbitrary (complex) intervals (a, b) using an affine map.
2
Both variants of the singular integral operators are implemented in SingularIntegralEquations.jl.
5.3. Multiple disjoint contours
Singular integral equations on a union of disjoint intervals = 1 ∪ 2 ∪ · · · ∪ d are covered in this framework. We can
decompose (61) as
⎛ ⎞ ⎛ ⎞
B1 B2 · · · Bd ⎛ ⎞ c
u1
⎜ A1,1 A1,2 · · · A1,d ⎟ ⎜ f1 ⎟
⎜ ⎟ ⎜u ⎟ ⎜ ⎟
⎜ A2,1 A2,2 · · · A2,d ⎟ ⎜ .2 ⎟ = ⎜ f2 ⎟ , (91)
⎜ . .. .. ⎟ ⎝ .. ⎠ ⎜ . ⎟
⎝ .. .
..
. . ⎠ ⎝ .. ⎠
ud
Ad,1 Ad,2 · · · Ad,d fd
where each Bi is a set of linear functionals and Ai , j = Ai | j . The diagonal blocks are equivalent to the previous case
considered, hence result in banded representations. The off-diagonal blocks can be constructed directly by expanding the
entire non-singular kernel in low rank form and using the compact functionals (−1,1) or I . The resulting representation
is, in fact, finite-dimensional and hence every block is banded.
Here, we show how a block-almost-banded infinite-dimensional system can be interlaced to be re-written as a single
infinite-dimensional and almost-banded system. Re-ordering both vectors (u1 , u2 , . . . , ud ) and (f1 , f2 , . . . , fd ) to:
* +
U = u 1,0 u 2,0 · · · ud,0 u 1,1 u 2,1 · · · ud,1 · · · (92)
* +
F = f 1,0 f 2,0 · · · f d,0 f 1,1 f 2,1 · · · f d,1 · · · (93)
amounts to a permutation of almost every row and column in (91). Define each entry of B and A by:
Bi , j = B{(i −1) mod d}+1, i+d−1 , j , (94)

d
Ai , j = A{(i −1) mod d}+1,{( j −1) mod d}+1, i+d−1 , j+d−1 , (95)
d d
where the last two indices in each term on the right-hand sides denote the entries of the functional or operator. This perfect
shuffle allows for the system (91) to be re-written as the almost-banded system

B c
U= . (96)
A F
5.4. Diagonal preconditioners for compactness
We now show that our formulations leads to equations whose operators are compact perturbations of the identity. For
well-posed (integral) equations, this ensures convergence [13]. We show this for the singular operators in equations (11)
and (12) defined on the canonical unit interval and in suitably chosen spaces. Note that a similar analysis is performed
in [21]. Since we are working in coefficient space, we consider the problem as defined in 2λ spaces. In the case of Chebyshev
expansions, this corresponds to Sobolev spaces of the transformed function u (cos θ).
Definition 12 (Olver and Townsend [13]). The space 2λ ⊂ C∞ is defined as the Banach space with norm:
/
0∞
0
u2 = 1 |uk |2 (k + 1)2λ < ∞. (97)
λ
k =0
Let Pn = ( I n , 0) be the projection operator.
Lemma 13. For the Dirichlet problem singular integral operator in (11), if takes the form (25) with A and B analytic in both x and
y and if we take R to be
⎛ 1
⎞
log 2
⎜ ⎟
⎜ 1 ⎟
⎜ ⎟ 2
R = 2⎜ 2 ⎟ : λ → 2λ−1 , (98)
⎜ 3 ⎟
⎝ ⎠
..
.
then
* +
L(−1,1) [π A ] + (−1,1) [π B ] R = I + K, (99)
where K : λ → λ is compact for λ ∈ R.

2 2
Proof. Since A (x, x) = −(2π )−1 , we let Ã (x, y ) ≡ A (x, y ) − A (x, x) and separate the operator (11) as:
L(−1,1) [π A (x, x)] + L(−1,1) [π Ã (x, y )] + (−1,1) [π B ]. (100)
It is straightforward to show
L(−1,1) [π A (x, x)]R = I : 2λ → 2λ . (101)
Then, we need to show that the remainder is compact. Since:
Pn L(−1,1) Pn − L(−1,1) → 0 as n → ∞, (102)
L(−1,1) : 2λ → 2λ is compact. Compactness of (−1,1) is implied by its finite-rank. Expanding Ã and B in low rank Chebyshev
approximants, we have:
⎛ ⎞
k

Ã
kB
π⎝ M0 [ Ã 1,i (x)]L(−1,1) M0 [ Ã 2,i ( y )] + M0 [ B 1,i (x)](−1,1) M0 [ B 2,i ( y )]⎠ R. (103)
i =1 i =1
Since A and B are analytic with respect to y, then for every i and for every λ ∈ R:
M0 [ Ã 2,i ( y )] : 2λ−1 → 2λ =⇒ M0 [ Ã 2,i ( y )]R : 2λ → 2λ , (104a)

2 2 2 2
M0 [ B 2,i ( y )] : λ−1 → λ =⇒ M0 [ B 2,i ( y )]R : λ → λ , (104b)
are bounded. Compactness follows from the linear combination of a product of bounded and compact operators being
compact. 2
Lemma 14. For the Neumann problem singular integral operator in (12), if takes the form (25) with A and B analytic in both x and
y and if we take R to be
⎛ ⎞
1
⎜ 1 ⎟
⎜ 2 ⎟
⎜ 1 ⎟ 2
R = −2 ⎜ 3 ⎟ : λ → 2λ+1 , (105)
⎜ 1 ⎟
⎝ 4 ⎠
..
.
then:
* +
HI [−π A ] + LI [π A ] + I [π B ] R = I + K, (106)
where the two primes indicate:

2
∂ 2 A (x, y) 22
A (x, y ) = , (107)
∂ x2 ∂ y 2 2x,y=(x,0),( y ,0)
and where K : 2λ → 2λ is compact for λ ∈ R.
Proof. Since A (x, x) = −(2π )−1 , we let Ã (x, y ) ≡ A (x, y ) − A (x, x) and separate the operator (12) as:
HI [−π A (x, x)] + HI [−π Ã (x, y )] + LI [π A ] + I [π B ]. (108)
It is straightforward to show:
HI [−π A (x, x)]R = I : 2λ → 2λ . (109)
Then, we need to show that the remainder is compact. Since:
Pn RPn − R → 0 as n → ∞, (110)
R : 2λ → 2λ is compact. Furthermore, showing boundedness of S0 , LI , I : 2λ → 2λ and HI : 2λ+1 → 2λ is straightforward.
Expanding Ã, A and B in low rank Chebyshev and ultraspherical approximants, we have:
⎛
k

Ã
π ⎝− M1 [ Ã 1,i (x)]HI M1 [ Ã 2,i ( y )]

i =1
⎛ ⎞⎞
k A k B

+ S0 ⎝ M0 [ A 1,i (x)]LI M1 [ A 2,i ( y )] + M0 [ B 1,i (x)]I M1 [ B 2,i ( y )]⎠⎠ R. (111)
i =1 i =1
Since A and B are analytic with respect to y, then for every i and for every λ ∈ R:
M1 [ Ã 2,i ( y )] : 2λ → 2λ+1 =⇒ HI M1 [ Ã 2,i ( y )] : 2λ → 2λ , (112)
are bounded. Compactness follows from the linear combination of a product of bounded and compact operators being
compact. 2
Remarks.
1. For complicated fundamental solutions whose bivariate low rank Chebyshev approximants have large degrees, precon-
ditioners such as those in Lemmas 13 and 14 allow for continuous Krylov subspace methods or conjugate gradients
on the normal equations to converge in a relatively fewer number of iterations compared with the un-preconditioned
operators. Furthermore, the low rank Chebyshev approximants allow for the operator-function product to be carried
out in O ((m + n) log(m + n)), where m is the largest degree of a multiplication operator and n is the degree of the
Chebyshev approximant of the solution. Iterative solvers are outside the scope of this article, however.
2. Operator preconditioners [66] can also be derived which would yield similar I + K results. However, working in coeffi-
cient space allows for a simpler exposition.
5.5. Numerical evaluation of Cauchy and log transforms on intervals
Fast and spectrally accurate numerical evaluation of the scattered far-field can be derived from Clenshaw–Curtis inte-
gration of the fundamental solution multiplied by the density. For each evaluation point, the fundamental solution can be
sampled at the 2N roots of the 2Nth degree Chebyshev polynomial, where N is the length of the polynomial representation
of the density. Since the resulting density may be as complicated3 as the fundamental solution itself, doubling the length is
sufficient to resolve the coefficients of the fundamental solution multiplied by the density.
It is well known that such an evaluation technique is inaccurate near the boundary [67]. In the context of Riemann–
Hilbert problems, spectrally accurate evaluation near and far from can be obtained by exact integration of a modified
Chebyshev series that encodes vanishing conditions at the endpoints.
Consider the modified Chebyshev series:
T̂ 0 (x) := 1, T̂ 1 (x) := x, T̂ n (x) := T n (x) − T n−2 (x), n ≥ 2. (113)
If we expand u in a Chebyshev series and this modified Chebyshev series:
∞
∞

u (x) = un T n (x) = ûn T̂ n (x), (114)
n =0 n =0
then we have the relation:

⎛ ⎞ ⎛ ⎞⎛ ⎞
u0 1 0 −1 û 0
⎜ u1 ⎟ ⎜ 1 0 −1 ⎟ ⎜ û 1 ⎟
⎜ ⎟=⎜ ⎟⎜ ⎟. (115)
⎝ u2 ⎠ ⎝ 1 0 −1 ⎠ ⎝ û 2 ⎠
.. .. .. .. ..
. . . . .
Therefore, any finite sequences {un }nN=0 and {ûn }nN=0 can be transformed to the other in O ( N ) operations, either via forward
application of the banded operator, or via an in-place back substitution.
Lemmas 16 and 17 contain formulæ for Cauchy transforms of weighted Chebyshev polynomials evaluated in the com-
plex plane. These were originally derived in this form in [68], based on results in [69–72]. These formulæ are adapted in
Lemma 18 for the log transform as well.
3
In the Helmholtz equation, for example, both the density and the fundamental solution are oscillatory with the same wavenumber.
Definition 15. Define the Joukowsky transform:
z + z −1
J ( z) := , (116)
2
and one of its inverses:
−1
√ √
J+ ( z) := z − z − 1 z + 1, (117)
which maps the slit plane C \ I to the unit disk.
The Joukowsky transform is useful for proving and summarizing the following results.
Lemma 16 (Lemma 5.6 [68]). For k ≥ 0:

- i −1 k +1
CI [ 1 − 2 U k ]( z) = J + ( z) . (118)
2
Proof. We verify that the Sokhotski–Plemelj lemma is satisfied. Note that for x = cos θ we have:
-
−1
lim J + (x ± i ) = x ∓ i 1 − x2 = cos θ ∓ i sin θ = e∓iθ . (119)
0
It follows that:
−1 −1
J+ (cos θ + i )k+1 − J + (cos θ − i )k+1 e−i(k+1)θ − ei(k+1)θ
lim = ,
0 2i 2i
-
= − sin(k + 1)θ = −U k (cos θ) 1 − cos2 θ. 2 (120)
Lemma 17 (Lemma 5.11 [68]). For k ≥ 2:

1 i
C(−1,1) √ ( z) = √ √ , (121)
1 − 2 2 z−1 z+1

iz i
C(−1,1) √ ( z) = √ √ − , and (122)
1 − 2 2 z−1 z+1 2
3 4
T̂ k −1
C(−1,1) √ ( z) = −i J + ( z)k−1 . (123)
1 − 2
Proof. The first two parts follow immediately from the Sokhotski–Plemelj lemma. The last part follows since:
cos kθ − cos(k − 2)θ
sin(k − 1)θ = . 2 (124)
2 sin θ
We extend these results here to the log transform.
Lemma 18.
2 2
2 −1 2
,- . −1
J+ ( z)2 log 2 J + ( z)2 + log 2
LI 1 − 2 ( z ) = − , (125)
34 2
4
, - . 1
−1
J+ ( z)k+2 −1
J+ ( z)k
L I U k 1 − 2 ( z ) = − , (126)
2 k+2 k

2 2
1 2 −1 2
L(−1,1) √ ( z) = − log 2 J + ( z)2 − log 2, (127)
1 − 2

−1
L(−1,1) √ ( z) = − J + ( z), (128)
1− 2
3 4
2 2 J −1 ( z)2
T̂ 2 2 −1 2
L(−1,1) √ ( z) = log 2 J + ( z)2 + log 2 − + , and (129)
1 − 2 2
3 4 3 4
−1 −1
T̂ k J+ ( z)k−2 J+ ( z)k
L(−1,1) √ ( z) = − . (130)
1 − 2 k−2 k
Proof. These formulæ follow from integrating the formulæ for Cauchy transforms and taking the real part. We can compute
the indefinite integrals directly [68, §5.4.4]:
z
1 −1
√ √ dz = − log J + ( z), (131)
z−1 z+1
z
1−z −1
√ √ dz = J + ( z), (132)
z−1 z+1
z −1
−1 J+ ( z)2 −1
2 J+ ( z) dz = − log J + ( z), and (133)
2
z −1 −1
−1 J+ ( z)k+1 J+ ( z)k−1
2 J+ ( z)k dz = − , for k ≥ 2. (134)
k+1 k−1
We also have the normalization for z → +∞:
1 1
f (x) log( z − x) dx = log z f (x) dx + O ( z−1 ). (135)
−1 −1
Note that:
−1 1
J+ ( z) ∼ + O( z−3 ) as z → ∞, (136)
2z
hence:
−1
log J + ( z) = − log z − log 2 + O( z−1 ) as z → ∞. 2 (137)
These formulæ can be generalized to other intervals, including in the complex plane, by using a straightforward change
of variables:

|b − a| b+a b−a b + a − 2z
L(a,b) f ( z) = L(−1,1) f +
2 2 2 b−a

1
|b − a| |b − a| b+a b−a
+ log f + x dx. (138)
2π 2 2 2
−1
6. Applications
6.1. The Faraday cage
The Faraday cage effect describes how a wire mesh can reduce the strength of the electric field within its confinement.
This phenomenon was described as early as 1755 by Franklin [73, §2-18] and in 1836 by Faraday [74]. While the descrip-
tion of the phenomenon is quite prevalent in undergraduate material on electrostatics, a standard mathematical analysis has
been missing until only recently by Martin [75] and Chapman, Hewett and Trefethen [76]. In [76], three different approaches
are considered for numerical simulations: a collocated least squares direct numerical calculation, a homogenized approxi-
mation via coupling of the solutions at multiple scales, and an approximation by point charges determined by minimizing a
quadratic energy functional.
In [76], it is shown that the shielding of a Faraday cage of circular wires centred at the roots of unity is a linear
phenomenon instead of providing exponential shielding as the number of wires tends to infinity for geometrically feasible
radii, i.e. radii that prevent overlapping. In their synopsis, it is claimed that a Faraday cage with any arbitrarily shaped
objects will not provide considerably different shielding in the asymptotic limit. Here, we confirm this observation with
infinitesimally thin plates of the same electrostatic capacity as wires4 angled normal to the vector from the origin to their
centres. Our numerical results are in excellent asymptotic agreement with those presented in [76]. Departing from the
practical case of normal plates, we also consider infinitesimally thin plates angled tangential to the vector from the origin
to their centres. In this case, we escape the practical material limit on the number of shields as an infinite number of plates
can be modelled independent of radial parameter.
4
This corresponds to plates of width 4r where r is the wire radius.
Fig. 2. Left: a plot of the solution u (x) with 10 normal plates with radial parameter r = 10−1 . Right: a plot of the solution u (x) with 40 tangential plates
with the same radial parameter, surpassing the material limit in the original numerical experiments [76]. In both contour plots, 31 contours are linearly
spaced between −2 and +1.
We seek to find the solution to the Laplace equation such that, in addition:
u (x) = 0, for x ∈ , (139a)

u (x) = u 0 , for x ∈ , (139b)
u (x) = log |x − y| + O (1), as |x − y| → 0, (139c)
u (x) = log |x| + o(1), as |x| → ∞. (139d)
Since this is a Dirichlet problem, we begin by splitting the solution u = u i + u s , where:
u i (x) = log |x − y| = 2π (x, y), (140)

is the source term with strength 2π located at y = (2, 0), as in [76]. We represent u in terms of a density with the s
single-layer potential equal to the effect of the logarithmic source. Alone, this represents a solution to the Laplace equation
with Dirichlet boundary conditions on . To satisfy condition (139b), we augment our system to ensure there is a constant
charge of zero on the wires and plates:

∂ us
d(y) = 0, (141)
∂n

though each wire may individually carry a different charge, and the unknown constant u 0 to accommodate this condition.
Fig. 2 shows the numerical results for shielding by normal and tangential plates. Fig. 3 shows a plot of the convergence of
the density coefficients and the field strength at the origin for various parameter values.
6.2. Helmholtz equation with Neumann boundary conditions
The mathematical treatment of the scattering of time-harmonic acoustic waves by infinitely long sound-hard obstacles in
three dimensions with simply-connected bounded cross-sections leads to the exterior problem for the Helmholtz equation:
( + k2 )u (x) = 0, for k ∈ R, x ∈ , (142a)

∂ u (x)
= 0, for x ∈ , (142b)
∂ n(x)
s
√ ∂u
lim r − iku s = 0, for r := |x|. (142c)
r →+∞ ∂r
Equation (142b) enforces sound-hard obstacles, while equation (142c) is the Sommerfeld radiation condition [77], an explicit
radiation condition at infinity. Consider an incident wave with wavenumber k and unit direction d:
u i (x) = e ikd·x . (143)

We wish to find the scattered field u s such that the sum u = u i + u s satisfies the Helmholtz equation in the exterior.
The fundamental solution of the Helmholtz equation is proportional to the cylindrical Hankel function of the first kind
of order zero [78, §8.405]:
Fig. 3. Left: a plot of the Cauchy error of successive approximants for the solution of Laplace’s equation with 10 normal plates with r = 10−1 , corresponding
to the left plot in Fig. 2, where n is the total number of degrees
( of freedom. The + indicates
( where the adaptive
( QR factorization ( terminates in double
precision, and with this approximation the forward error is (u 0 + S [∂ u s /∂ n] d(y) − u i (2 = 2.90 × 10−15 and ( [∂ u s /∂ n] d(y)(2 = 1.47 × 10−15 . Right:
a plot of the field strength in the centre of the cage versus the number of plates. The dashed lines represent results for normal plates, while the solid lines
represent results for tangential plates of the same electrostatic capacity. The normal and tangential plates exhibit different asymptotic scalings.
√ √
Fig. 4. Acoustic scattering with Neumann boundary conditions from an incident wave with k = 100 and d = (1/ 2, −1/ 2). Left: a plot of the numerical
ranks of J 0 (k|x − y|) connecting domain i to domain j, where it can be seen that interaction between domains is relatively weaker than self-interaction.
Right: a plot of the total solution. 1,392 degrees of freedom are required to represent the piecewise density in double precision.
i (1 )
(x, y) = H 0 (k|x − y|), (144)
4
and the Riemann function is also well known [41] for the Helmholtz equation:
-
R( z, ζ, z0 , ζ0 ) = J 0 (k ( z − z0 )(ζ − ζ0 )). (145)
Fig. 4 shows the rank structure of the bivariate kernels and the total solution with a set of randomly generated screens
between [−3, 3]. N.B. it is known that [79] collinear screens have reduced off-diagonal numerical ranks compared with
randomly oriented screens.
6.3. Gravity Helmholtz equation with Dirichlet boundary conditions
The Helmholtz equation in a linearly stratified medium:
( + E + x2 )u (x) = 0, for E ∈ R, x ∈ , (146a)

u (x) = 0, for x ∈ , (146b)
2 -
22
1 2 ∂u 2
lim √ 2 − i E + x2 u 22 dx1 = 0, (146c)
x2 →+∞ E + x2 2 ∂ x2
R
2 2
2 ∂ u 22
lim |u |2 + 22 2 dx1 = 0, (146d)
x2 →−∞ ∂ x2 2
R
L 2 2
2 ∂ u 22
lim lim |u |2 + 22 2 dx2 = 0, (146e)
L →+∞ x1 →±∞ ∂ x1 2
−L
models quantum particles of fixed energy in a uniform gravitational field [57]. Equation (146b) enforces sound-soft obstacles,
while equations (146c)–(146e) form an explicit radiation condition at infinity derived in [57].
The fundamental solution of the Helmholtz equation in a linearly stratified medium is derived in [80]:
∞

1 |x − y|2 x2 + y 2 1 3 dt
(x, y) = exp i + E+ t− t . (147)
4π 4t 2 12 t
0
Numerical evaluation via the trapezoidal rule [81] along a contour of approximate steepest descent on the order of 105
evaluations per second is reported in [57]. This equation is also known as the gravity Helmholtz equation.
Consider an incident fundamental solution with energy E and source y:
u i (x) = (x, y). (148)
We wish to find the scattered field u such that the sum u = u + u satisfies the gravity Helmholtz equation in the exterior.
s i s
In addition to the fundamental solution, we require the Riemann function of the PDO. With the prospect of deriving a fast
numerical evaluation in future work, we prove the following theorem in Appendix A.
E z−ζ
Theorem 19. The Riemann function of the gravity Helmholtz equation, where c (x1 , x2 ) = E + x2 and therefore C ( z, ζ ) = 4
+ 8i
has
the power series:
∞
∞
R( z, ζ, z0 , ζ0 ) = 1 + A i , j ( z − z0 )i (ζ − ζ0 ) j , (149)
i =1 j =1
where the coefficients A i , j satisfy (A.3)–(A.5), and the integral representation:
γ+i∞
1 1 5
V (u , v ) = - exp 8i Ẽ (s2 − u /4i)1/2 − s
2π i s2 − u /4i
γ −i∞

8i
+ s3 − (s2 − u /4i)3/2 + ( v − u )s ds, (150)
3
E z0 − ζ0
where R( z, ζ, z0 , ζ0 ) = V ( z − z0 , ζ − ζ0 ) and where Ẽ = + .
4 8i
Fig. 5 shows the total solution to the gravity Helmholtz equation with Dirichlet boundary conditions and the 2-norm
condition number of the truncated and preconditioned system.
6.4. Helmholtz equation with nearly singular Dirichlet boundary data
In this application, we consider the Helmholtz equation with nearly singular Dirichlet boundary data. Consider the scat-
tering of a collection of point sources arbitrarily close to a sound-soft obstacle. If we parameterize the locations of the point
sources by the family of Bernstein ellipses E ρ , then we know that the Chebyshev series representation of the incident wave
will have degree which scales as n = O ((log ρ )−1 ) as ρ → 1. Additionally, as the point sources approach the boundary, the
integral operator has bandwidth O (k), independent of the Bernstein ellipse parameter. This is a challenging scenario for
conventional integral equation solvers since a piecewise polynomial approximation to the nearly singular boundary data
may not be much more efficient than a global representation. Furthermore, if the point sources are allowed to move freely
on the Bernstein ellipse, then no adaptivity may be used to uniformly accelerate the solvers. The demonstrations in this
section are also applicable to the important problem of many micro swimmers in Stokes flow approaching an obstacle, as
the swimmers can be modelled as point sources, see [82].
This set of problems completely demonstrates the scaling O ((mx + m y )2 n) of our algorithm: the bandwidth scales with
the wavenumber, and the degree scales with the reciprocal of the log of the Bernstein ellipse parameter. Additionally, a
Fig. 5. Acoustic scattering with Dirichlet boundary conditions from an incident fundamental solution (x, y) with y = (0, −5) and E = 20 against the
sound-soft intervals ((−10, −3), (−5, 0)) ∪ ((−2, 5), (2, 5)) ∪ ((5, 0), (10, −3)). Left: a plot of the 2-norm condition number of the truncated and precon-
ditioned system with n degrees of freedom. Right: a plot of the total solution. 332 degrees of freedom are required to represent the piecewise density in
double precision.
Fig. 6. Timings to solve the Helmholtz equation (left) with nearly singular boundary data, and a sample of the solution (right). Left: timings are illustrated
for three wavenumbers and multiple Bernstein ellipse parameters resulting in Chebyshev expansions reaching degrees on the order of half a million. The
high data set show the partial QR factorization and the back substitution, while the low data set show the reduced time with a cached QR factorization.
Right: a sample solution at k = 50, with 100 point sources uniformly distributed in angle on the top half of the Bernstein ellipse E 1.05 with charges +1 on
the right half and −1 on the left half, resulting in the symmetric output.
partial QR factorization of the singular integral operator may be cached or precomputed,5 resulting in the reduced O ((mx +
m y )n) complexity for additional solves. Fig. 6 shows the scalings of the computation for three wavenumbers and varying
Bernstein ellipse parameters. The figure also shows the solution of the Helmholtz equation with 100 nearby source terms.
7. Numerical discussion & outlook
The software package SingularIntegralEquations.jl [18] written in the Julia programming language [16,17]
implements the banded singular integral operators, methods relating to bivariate function approximation and construction
with diagonal singularities, fast & spectrally accurate numerical evaluation of scattered fields and several examples including
those described in this work. Built on top of ApproxFun.jl, SingularIntegralEquations.jl uses the adaptive
QR factorization described in [13] and acts as an extension to the framework for infinite-dimensional linear algebra. All
numerical simulations are performed on a MacBook Pro with a 2.8 GHz Intel Core i7-4980HQ processor and 16 GB of RAM.
While timings are continuously being improved, Table 1 shows the current timings to solve the problems in section 6.
All the numerical problems relating to our applications have been abstracted so that to explore a new elliptic PDE in
SingularIntegralEquations.jl, the user only needs a fast evaluation of the fundamental solution and its Riemann
function.
5
The cached QR factorization can be adaptively grown without re-computing from scratch by exploiting the fact that the operator is banded below, thus
the number of degrees of freedom (n) needed to resolve the solution within a prescribed tolerance need not be known a priori. This automatic caching of
the QR factorization is implemented in ApproxFun.jl.
Table 1
Calculation times in seconds to solve the problems in section 6. Evaluation of the scattered field is
reported per target. Timings for the Laplace equation are for 10 normal plates.
Kernel assembly Adaptive QR Evaluation of scattered field

Laplace 0.888 0.518 0.0000135
Helmholtz (k = 100) 1.73 67.6 0.00652
Gravity Helmholtz (E = 20) 3.11 1.20 0.0139
Fig. 7. Left: An illustration of a scenario that benefits from the abstraction of hierarchical matrix factorizations to our hierarchical operators. The Dirichlet
solution of the Helmholtz equation with k = 100 and one incident source in the North Sea. Right: Idealized fluid flow around three obstacles.
For problems involving a union of a considerably large number of domains, the current method of interlacing all opera-
tors can be improved. In future work on fractal screens motivated by [83], alternative algorithms based on hierarchical block
diagonalization via a symmetrized Schur complement [33] may be explored specifically exploiting the low rank off-diagonal
structure arising from coercive singular integral operators of elliptic PDOs. This is close in spirit to the Fast Multipole
Method [30], but applied to the banded representation of the singular integral operators, instead of discretizations arising
from quadrature rules. A preliminary result in this direction is shown in the left side of Fig. 7.
As illustrated in subsection 6.2 on the acoustic scattering of the Helmholtz equation with Neumann boundary conditions,
SingularIntegralEquations.jl supports higher order diagonal singularities. Future work may explore the feasibility
of combining automatic differentiation and differentiation of Chebyshev interpolants to automate the construction of the
operators with higher order singularities such that the user need only enter the fundamental solution with its logarithmic
splitting described by (25).
The approach developed in this article is also adaptable to other domains such as disjoint unions of circles and polyno-
mial maps of intervals and circles, see the right side of Fig. 7 for an example calculated using SingularIntegralEqua-
tions.jl of idealized fluid flow over three domains: an interval, a circle and a polynomial map of an interval. To take into
account circles, a similar analysis is straightforward with Laurent polynomials in place of weighted Chebyshev polynomials.
However, a combined field formulation is beneficial to ensure well-conditioning when the solution of the exterior problem
is near an eigenmode of the interior problem. Equations over maps of the unit interval and circle can also be reduced to
numerically banded singular integral operators via approximating the map by a polynomial and using the spectral mapping
theorem. The key formula in the Hilbert case is derived in [68, Theorem 5.32], which implies that the Hilbert transform
over a polynomial map of the unit interval can be reduced to a compact perturbation of the Hilbert transform over the unit
interval. Expanding on this result, as well as adapting the procedure to log transforms, will be the topic of a subsequent
publication. Future work may consider the use of these modified Chebyshev series for banded operators when two disjoint
contours are in close proximity. When two or more contours coalesce, banded singular integral operators will depend on
the ability to produce the orthogonal polynomials associated with that domain. Densities of the single- and double-layer
potentials will have singularities on domains with cusps. Such an analysis is undetermined.
As discussed in [57], the fundamental solution of the gravity Helmholtz equation has an analogy to the Schrödinger
equation with a linear potential. The Helmholtz equation with a parabolic refractive index shares the same analogy and
the fundamental solution is also known [9,10]. Parabolic refractive indices occur when considering the shielding of optical
fibres, leading to Gaussian beams. Scattering problems in this context may shed light on the effects when optical fibres are
occluded. Fast and accurate numerical evaluation of the fundamental solution as well as the Riemann function may also be
possible via the trapezoidal rule.
An important area of future research is extending the method to higher dimensional singular integral equations. The
ultraspherical spectral method was extended to automatically solve general linear partial differential equations on rectan-
gles [84] and the ideas used to do this successfully may well translate to singular integral equations.
Acknowledgements
We wish to thank Jared Aurentz, Folkmar Bornemann, Dave Hewett, Alex Townsend and Nick Trefethen for stimulating
discussions related to this work. We acknowledge the generous support of the Natural Sciences and Engineering Research
Council of Canada grant 6790 - 454127 - 2014 (RMS) and the Australian Research Council grant DE130100333 (SO).
Appendix A. Proof of Theorem 19
To immediately satisfy the boundary conditions (23), we start with the ansatz:
∞
∞
R( z, ζ, z0 , ζ0 ) = 1 + A i , j ( z − z0 )i (ζ − ζ0 ) j , (A.1)
i =1 j =1
and we insert it into the integral equation (24):

∞
∞
A i , j ( z − z0 )i (ζ − ζ0 ) j
i =1 j =1
⎛ ⎞
z ζ ∞
∞
E t−τ
+ + ⎝1 + A i , j (t − z0 )i (τ − ζ0 ) j ⎠ dτ dt = 0. (A.2)
4 8i
z0 ζ0 i =1 j =1
With the initial values:
E ( z 0 − ζ0 ) 1 1
A 1,1 = − − , A 2,1 = − , A 1,2 = , A 2,2 = A 21,1 /4, (A.3)
4 8i 16i 16i
and the additional values:
A i ,1 = A 1,i = 0, for i > 2, A i , j = 0, for i ≤ 0, j ≤ 0, (A.4)

the coefficients are found to satisfy in general:

E z 0 − ζ0 1 1
i j A i, j + + A i −1 , j −1 − A i −1 , j −2 + A i −2 , j −1 = 0. (A.5)
4 8i 8i 8i
The growth in the constant in front of A i , j ensures that coefficients decay at least exponentially fast, hence the power series
converges for all z and ζ .
To get an integral representation for the Riemann function, we start from the differential equation it satisfies after the
change of variables u = z − z0 and v = ζ − ζ0 :

∂2V u−v E z 0 − ζ0
+ Ẽ + V = 0, Ẽ = + , (A.6)
∂ u∂ v 8i 4 8i
together with V (0, v ) = V (u , 0) = 1.
Taking the Laplace transform:
∞
L{ f }(s) = f ( v )e −sv dv , (A.7)
0
of the differential equation, we obtain:
∂ V̂ 1 ∂ V̂ u 1
s + + Ẽ + V̂ = 0, V̂ (0, s) = . (A.8)
∂u 8i ∂ s 8i s
Using the method of characteristics for this first-order PDE, we obtain the general solution as:

8is3
V̂ (u , s) = exp −8i Ẽ s + − us + f (4is2 − u ) . (A.9)
3
The particular solution satisfying the initial condition is:

1 8i
V̂ (u , s) = - exp 8i Ẽ (s2 − u /4i)1/2 − s + s3 − (s2 − u /4i)3/2 − us . (A.10)
s2 − u /4i 3
Inverting the Laplace transform using the Bromwich integral, we find the solution (150).
References
[1] F. Erdogan, Fracture mechanics, Int. J. Solids Struct. 37 (2000) 171–183.

[2] D. Colton, R. Kress, Integral Equation Methods in Scattering Theory, Wiley, 1983.
[3] R. Kress, On the numerical solution of a hypersingular integral equation in scattering theory, J. Comput. Appl. Math. 61 (1995) 345–360.
[4] D. Huybrechs, S. Vandewalle, A sparse discretization for integral equation formulations of high frequency scattering problems, SIAM J. Sci. Comput. 29
(2007) 2305–2328.
[5] R. Kress, Linear Integral Equations, Applied Mathematical Sciences, vol. 82, Springer, 2010.
[6] D.P. Hewett, S. Langdon, J.M. Melenk, A high frequency hp boundary element method for scattering by convex polygons, SIAM J. Numer. Anal. 51 (2013)
629–653.
[7] D.L. Young, S.J. Jane, C.M. Fan, K. Murugesan, C.C. Tsai, The method of fundamental solutions for 2D and 3D Stokes problems, J. Comput. Phys. 211
(2006) 1–8.
[8] S. Olver, A general framework for solving Riemann–Hilbert problems numerically, Numer. Math. 122 (2012) 305–340.
[9] C.C. Constantinou, Path-Integral Analysis of Passive, Graded-Index Waveguides Applicable to Integrated Optics, Ph.D. thesis, University of Birmingham,
1991.
[10] E.J. Heller, Wavepacket dynamics and quantum chaology, in: Chaos et physique quantique (Les Houches, 1989), North Holland, Amsterdam, 1991,
pp. 547–664.
[11] N.I. Muskhelishvili, Singular Integral Equations, 2nd edition, Dover Publications Inc., P. Noordhoff, Groningen, Holland, 1953.
[12] A. Townsend, L.N. Trefethen, An extension of Chebfun to two dimensions, SIAM J. Sci. Comput. 35 (2013) C495–C518.
[13] S. Olver, A. Townsend, A fast and well-conditioned spectral method, SIAM Rev. 55 (2013) 462–489.
[14] S. Olver, A. Townsend, A practical framework for infinite-dimensional linear algebra, in: Proceedings of the First Workshop for High Performance
Technical Computing in Dynamic Languages, 2014, pp. 57–62.
[15] S. Olver, G. Goretkin, R.M. Slevinsky, A. Townsend, https://github.com/ApproxFun/ApproxFun.jl, GitHub.
[16] J. Bezanson, S. Karpinski, V.B. Shah, A. Edelman, Julia: a fast dynamic language for technical computing, arXiv:1209.5145, 2012.
[17] J. Bezanson, A. Edelman, S. Karpinski, V.B. Shah, Julia: a fresh approach to numerical computing, arXiv:1411.1607, 2014.
[18] S. Olver, R.M. Slevinsky, https://github.com/ApproxFun/SingularIntegralEquations.jl, GitHub.
[19] D. Berthold, P. Junghanns, New error bounds for the quadrature method for the solution of Cauchy singular integral equations, SIAM J. Numer. Anal. 30
(1993) 1351–1372.
[20] O.P. Bruno, S.K. Lintner, Second-kind integral solvers for TE and TM problems of diffraction by open arcs, arXiv:1204.3701, 2012.
[21] S.K. Lintner, O.P. Bruno, A generalized Calderón formula for open-arc diffraction problems: theoretical considerations, arXiv:1204.3699, 2012.
[22] D. Elliott, Orthogonal polynomials associated with singular integral equations having a Cauchy kernel, SIAM J. Math. Anal. 13 (1982) 1041–1052.
[23] D. Elliott, The classical collocation method for singular integral equations, SIAM J. Numer. Anal. 19 (1982) 816–832.
[24] B.D. Galerkin, Expansions in stability problems for elastic rods and plates (in Russian), Vestn. Inzh. 19 (1915) 897–908.
[25] W. Śmigaj, T. Betcke, S. Arridge, J. Phillips, M. Schweiger, Solving boundary integral problems with BEM++, ACM Trans. Math. Softw. 41 (2015) 6:1–6:40.
[26] R. Kress, Boundary integral equations in time-harmonic acoustic scattering, Math. Comput. Model. 15 (1991) 229–243.
[27] B.K. Alpert, Hybrid Gauss-trapezoidal quadrature rules, SIAM J. Sci. Comput. 20 (1999) 1551–1584.
[28] J. Helsing, R. Ojala, On the evaluation of layer potentials close to their sources, J. Comput. Phys. 227 (2008) 2899–2921.
[29] S. Hao, A.H. Barnett, P.G. Martinsson, P. Young, High-order accurate methods for Nyström discretization of integral equations on smooth curves in the
plane, Adv. Comput. Math. 40 (2014) 245–272.
[30] L. Greengard, V. Rokhlin, A fast algorithm for particle simulations, J. Comput. Phys. 73 (1987) 325–348.
[31] W. Hackbusch, Z.P. Nowak, On the fast matrix multiplication in the boundary element method by panel clustering, Numer. Math. 54 (1989) 463–491.
[32] S. Ambikasaran, E. Darve, An O ( N log N ) fast direct solver for partial hierarchically semi-separable matrices with application to radial basis function
interpolation, J. Sci. Comput. 57 (2013) 477–501.
[33] A. Aminfar, S. Ambikasaran, E. Darve, A fast block low-rank dense solver with applications to finite-element matrices, arXiv:1403.5337, 2014.
[34] A.N. Krylov, On the numerical solution of equations by which are determined in technical problems the frequencies of small vibrations of material
systems (in Russian), Izv. Akad. Nauk SSSR 7 (1931) 491–539.
[35] P.G. Martinsson, V. Rokhlin, A fast direct solver for scattering problems involving elongated structures, J. Comput. Phys. 221 (2007) 288–302.
[36] D.P. Hewett, S. Langdon, S.N. Chandler-Wilde, A frequency-independent boundary element method for scattering by two-dimensional screens and
apertures, arXiv:1401.2786, 2014.
[37] A. Frenkel, A Chebyshev expansion of singular integral equations with a logarithmic kernel, J. Comput. Phys. 51 (1983) 326–334.
[38] Y.-S. Chan, Hypersingular Integrodifferential Equations and Applications to Fracture Mechanics of Homogenous and Functionally Graded Materials with
Strain-Gradient Effects, Ph.D. thesis, University of California, 2001.
[39] Y.-S. Chan, A.C. Fannjiang, G.H. Paulino, Integral equations with hypersingular kernels – theory and applications to fracture mechanics, Int. J. Eng. Sci.
41 (2003) 683–720.
[40] A. Frenkel, A Chebyshev expansion of singular integrodifferential equations with a ∂ 2 ln |s − t |/∂ s∂ t kernel, J. Comput. Phys. 51 (1983) 335–342.
[41] I.N. Vekua, New Methods for Solving Elliptic Equations, North Holland, 1967.
[42] M. Costabel, M. Dauge, On representation formulas and radiation conditions, Math. Methods Appl. Sci. 20 (1997) 133–150.
[43] E.P. Stephan, W.L. Wendland, An augmented Galerkin procedure for the boundary integral method applied to two-dimensional screen and crack prob-
lems, Appl. Anal. 18 (1984) 183–219.
[44] W.L. Wendland, E.P. Stephan, A hypersingular boundary integral method for two-dimensional screen and crack problems, Arch. Ration. Mech. Anal. 112
(1990) 363–390.
[45] P.R. Garabedian, Partial Differential Equations, John Wiley & Sons, Inc., New York, 1964.
[46] J.P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd edition, Dover Publications Inc., 2000.
[47] J.C. Mason, D.C. Handscomb, Chebyshev Polynomials, CRC Press, 2002.
[48] L.N. Trefethen, Approximation Theory and Approximation Practice, SIAM, 2012.
[49] J.W. Cooley, J.W. Tukey, An algorithm for the machine calculation of complex Fourier series, Math. Comput. 19 (1965) 297–301.
[50] M. Frigo, S.G. Johnson, The design and implementation of FFTW3, Proc. IEEE 93 (2005) 216–231.
[51] C.W. Clenshaw, A note on the summation of Chebyshev series, Math. Comput. 9 (1955) 118–120.
[52] Z. Battles, L.N. Trefethen, An extension of Matlab to continuous functions and operators, SIAM J. Sci. Comput. 25 (2004) 1743–1770.
[53] D.S. Watkins, Fundamentals of Matrix Computations, third edition, Wiley, 2010.
[54] A. Townsend, Computing with Functions in Two Dimensions, Ph.D. thesis, University of Oxford, 2014.
[55] T.A. Driscoll, N. Hale, L.N. Trefethen (Eds.), Chebfun Guide, Pafnuty Publications, 2014.
[56] A. Townsend, L.N. Trefethen, Continuous analogues of matrix factorizations, Proc. R. Soc. A 471 (2015) 20140585.
[57] A.H. Barnett, B.J. Nelson, J.M. Mahoney, High-order boundary integral equation solution of high frequency wave scattering from obstacles in an un-
bounded linearly stratified medium, J. Comput. Phys. 297 (2015) 407–426.
[58] F.W.J. Olver, Numerical solution of second-order linear difference equations, J. Res. Natl. Bur. Stand. 71B (1967) 111–129.
[59] W. Gautschi, Orthogonal Polynomials: Computation and Approximation, Clarendon Press, Oxford, UK, 2004.
[60] N. Hale, A. Townsend, A fast FFT-based discrete Legendre transform, arXiv:1505.00354, 2015.
[61] Y.V. Sokhotski, On Definite Integrals and Functions Utilized for Expansions into Series (in Russian), Ph.D. thesis, University of St. Petersburg, 1873.
[62] J. Plemelj, Ein Ergänzungssatz zur Cauchyschen Integraldarstellung analytischer Funktionen, Randwerte betreffend, Monatshefte Math. Phys. 19 (1908)
205–210.
[63] P.A. Martin, Exact solution of a simple hypersingular integral equation, J. Integral Equ. Appl. 4 (1992) 197–204.
[64] G. Monegato, Numerical evaluation of hypersingular integrals, J. Comput. Appl. Math. 50 (1994) 9–31.
[65] F.W. King, Hilbert Transforms, vol. 1, Cambridge University Press, 2009.
[66] R. Hiptmair, Operator preconditioning, Comput. Math. Appl. 52 (2006) 699–706.
[67] A.H. Barnett, Evaluation of layer potentials close to the boundary for Laplace and Helmholtz problems on analytic planar domains, SIAM J. Sci. Comput.
36 (2014) A427–A451.
[68] T. Trogdon, S. Olver, Riemann–Hilbert Problems, Their Numerical Solution and the Computation of Nonlinear Special Functions, SIAM, 2015.
[69] S. Olver, Numerical solution of Riemann–Hilbert problems: Painlevé II, Found. Comput. Math. 11 (2011) 153–179.
[70] S. Olver, Computing the Hilbert transform and its inverse, Math. Comput. 80 (2011) 1745–1767.
[71] S. Olver, Computation of equilibrium measures, J. Approx. Theory 163 (2011) 1185–1207.
[72] S. Olver, T. Trogdon, Numerical solution of Riemann–Hilbert problems: random matrix theory and orthogonal polynomials, Constr. Approx. 39 (2013)
101–149.
[73] J.D. Kraus, Electromagnetics, 4th edition, McGraw–Hill, 1992.
[74] M. Faraday, Experimental Researches in Electricity, vol. 1, reprinted from Philosophical Transactions of 1831–1838
, Richard and John Edward Taylor, London, 1839.
[75] P.A. Martin, On acoustic and electric Faraday cages, Proc. R. Soc. A 470 (2014) 20140344.
[76] S.J. Chapman, D.P. Hewett, L.N. Trefethen, Mathematics of the Faraday cage, SIAM Rev. 57 (2015) 398–417.
[77] A. Sommerfeld, Partial Differential Equations in Physics, Pure and Applied Mathematics: A Series of Monographs and Textbooks, vol. 1, Academic Press,
New York, NY, 1949.
[78] I.S. Gradshteyn, I.M. Ryzhik, Table of Integrals, Series, and Products, 7th edition, Elsevier Academic Press, Burlington, MA, 2007.
[79] E. Michielssen, A. Boag, W.C. Chew, Scattering from elongated objects: direct solution in O ( N log2 N ) operations, IEE Proc., H Microw. Antennas Propag.
143 (1996) 277–283.
[80] C. Bracher, W. Becker, S.A. Gurvitz, M. Kleber, M.S. Marinov, Three-dimensional tunneling in quantum ballistic motion, Am. J. Phys. 66 (1998) 38–48.
[81] L.N. Trefethen, J.A.C. Weideman, The exponentially convergent trapezoidal rule, SIAM Rev. 56 (2014) 385–458.
[82] A.M.J. Davis, D.G. Crowdy, Matched asymptotics for a treadmilling low-Reynolds-number swimmer near a wall, Q. J. Mech. Appl. Math. 66 (2012)
53–73.
[83] S.N. Chandler-Wilde, D.P. Hewett, Acoustic scattering by fractal screens: mathematical formulations and wavenumber-explicit continuity and coercivity
estimates, arXiv:1401.2805, 2014.
[84] A. Townsend, S. Olver, The automatic solution of partial differential equations using a global spectral method, J. Comput. Phys. 299 (2015) 106–123.

Slevinsky JCP 332 2017 290

Uploaded by

Copyright:

Available Formats

Slevinsky JCP 332 2017 290

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Slevinsky JCP 332 2017 290

Uploaded by

Copyright:

Available Formats

Journal of Computational Physics 332 (2017) 290–315

Contents lists available at ScienceDirect

Journal of Computational Physics

A fast and well-conditioned spectral method for singular

where K 1 , . . . , K 4 are continuous bivariate kernels.

2. Boundary integral equations in two dimensions

Lx {} = −δ(x − y), (3)

2.1. Exterior scattering problems

Let be a union of disjoint bounded smooth open arcs in R2 and let  := R2 \ .

Deﬁnition 2 (Radiation condition at inﬁnity [42]). We say that u ∈ H loc

u (x) = −S [∂ u /∂ n] (x) + D [u] (x), for x ∈ . (8)

u i (x) + u s (x) = 0, for x ∈ . (9)

2.2. Riemann functions

In addition to the PDO in (2), consider its adjoint:

z = x1 + ix2 , ζ = x1 − ix2 , z0 = y 1 + i y 2 , ζ0 = y 1 − i y 2 , (14)

L and L∗ take the form:

L̂∗ {R} = 0, (21)

Remarks. It is straightforward to reformulate (21)–(23) to the following integral equation:

(x, y) = A (x, y) log |x − y| + B (x, y), (25)

3. Practical approximation theory

3.1. One dimension

Let K be the ﬁeld of R or C. A function f : I → K is of bounded total variation if:

T n (x) = cos(n cos−1 (x)), for n ∈ N0 , and x ∈ I. (27)

which interpolates f at the Chebyshev points of the ﬁrst kind:

2. if f has a kth derivative in I of bounded variation for some k ≥ 1, f − p N ∞ = O ( N −k ) as N → ∞; and,

3.2. Two dimensions

where the inner inﬁmum is taken over all rank-k functions.

3.3. An algorithm to extract the splitting of a fundamental solution

1. approximate A (x, y) ≡ − 21π R(x1 + ix2 , x1 − ix2 , y 1 + i y 2 , y 1 − i y 2 ); and subsequently,

4. The ultraspherical spectral method

so that u (x) can be identiﬁed by a vector of its Chebyshev coeﬃcients u = (u 0 , u 1 , . . .) .

Then it is shown in [13] that multiplication can be represented as a Toeplitz-plus-Hankel-plus-rank-one operator:

(λ) 2(n + λ) (λ) n + 2λ − 1 (λ)

4.1. Almost-banded spectral methods in other bases

xπn (x) = αn πn+1 (x) + βn πn (x) + γn πn−1 (x), (49)

5. Ultraspherical spectral method for singular integral equations

In the following deﬁnitions, we identify C with R2 and let be bounded in C.

| f (x) − f (y)| ≤ C |x − y|α , for x, y ∈ . (52)

where the integral is understood as the Cauchy principal value:

Lemma 10 (Sokhotski–Plemelj [61,62]). If f ∈ C 0,α (), then:

H f ( z) = i[C + + C − ] f ( z), (57)

where the integral is understood as the Hadamard ﬁnite-part [63,64]:

where ( z; ρ ) is deﬁned as in Deﬁnition 9.

5.1. Inverse square root endpoint singularities

Indeed, the Hilbert transform of weighted Chebyshev polynomials is known [65]:

or as a compact functional on the weighted Chebyshev coeﬃcients:

5.2. Square root endpoint singularities

or as a compact functional on the weighted Chebyshev basis:

5.3. Multiple disjoint contours

Bi , j = B{(i −1) mod d}+1, i+d−1 , j , (94)

5.4. Diagonal preconditioners for compactness

Let Pn = ( I n , 0) be the projection operator.

where K : λ → λ is compact for λ ∈ R.

L(−1,1) [π A (x, x)] + L(−1,1) [π Ã (x, y )] + (−1,1) [π B ]. (100)

L(−1,1) [π A (x, x)]R = I : 2λ → 2λ . (101)

Then, we need to show that the remainder is compact. Since:

Pn L(−1,1) Pn − L(−1,1) → 0 as n → ∞, (102)

M0 [ Ã 2,i ( y )] : 2λ−1 → 2λ =⇒ M0 [ Ã 2,i ( y )]R : 2λ → 2λ , (104a)

where the two primes indicate:

Lx {} = −δ(x − y), (3)

Let be a union of disjoint bounded smooth open arcs in R2 and let := R2 \ .

u (x) = −S [∂ u /∂ n] (x) + D [u] (x), for x ∈ . (8)

(x, y) = A (x, y) log |x − y| + B (x, y), (25)

so that u (x) can be identiﬁed by a vector of its Chebyshev coeﬃcients u = (u 0 , u 1 , . . .) .

Bi , j = B{(i −1) mod d}+1, i+d−1 , j , (94)

where K : λ → λ is compact for λ ∈ R.

L(−1,1) [π A (x, x)] + L(−1,1) [π Ã (x, y )] + (−1,1) [π B ]. (100)

L(−1,1) [π A (x, x)]R = I : 2λ → 2λ . (101)

Pn L(−1,1) Pn − L(−1,1) → 0 as n → ∞, (102)

M0 [ Ã 2,i ( y )] : 2λ−1 → 2λ =⇒ M0 [ Ã 2,i ( y )]R : 2λ → 2λ , (104a)

and where K : 2λ → 2λ is compact for λ ∈ R.

HI [−π A (x, x)] + HI [−π Ã (x, y )] + LI [π A ] + I [π B ]. (108)

HI [−π A (x, x)]R = I : 2λ → 2λ . (109)

Pn RPn − R → 0 as n → ∞, (110)

M1 [ Ã 2,i ( y )] : 2λ → 2λ+1 =⇒ HI M1 [ Ã 2,i ( y )] : 2λ → 2λ , (112)

u (x) = 0, for x ∈ , (139a)

u i (x) = log |x − y| = 2π (x, y), (140)

( + k2 )u (x) = 0, for k ∈ R, x ∈ , (142a)

( + E + x2 )u (x) = 0, for E ∈ R, x ∈ , (146a)

u i (x) = (x, y). (148)