Solving integral equations on piecewise smooth boundaries using the RCIP method: a tutorial

Recursively compressed inverse preconditioning (RCIP) is a kernel-independent and purely numerical method for solving Fredholm second kind boundary integral equations in situations where the boundary shape induces a non-smooth behavior in the solution. The method originated in 2008 within a scheme for Laplace's equation in two-dimensional domains with corners. In a series of subsequent papers the method was then refined and extended as to apply to integral equation formulations of a broad range of boundary value problems in physics and engineering. The purpose of the present tutorial is threefold: First, to review the RCIP method in a simple setting. Second, to show how easily the method can be implemented in Matlab. Third, to present new applications.


Introduction
This tutorial is about an efficient numerical solver for elliptic boundary value problems in domains whose boundaries contain some sort of singular points. Such a solver is useful for applications in physics and engineering, where computational domains of interest often have corners, triple junctions, and close-to-touching boundary parts. Furthermore, these problems are difficult to solve irrespective of what numerical method is used. The reason being that the solution, or the quantity representing the solution, often exhibits a non-smooth behavior close to boundary singularities. That behavior is hard to resolve by polynomials, which underlie most approximation schemes. Mesh refinement is needed. This is costly and may lead to artificial illconditioning and the loss of accuracy.
The numerical solver we propose takes its starting point in an integral equation reformulation of the boundary value problem at hand. We assume that the problem can be modeled as a Fredholm second kind integral equation with compact integral operators away from singular boundary points and whose solution is a layer density representing the solution to the original problem. We seek a discrete approximation to the layer density using Nyström discretization [2,Chapter 4]. At the heart of the solver lies an integral transform whose inverse modifies the kernels of the integral operators in such a way that the layer density becomes piecewise smooth and simple to resolve by polynomials. The inverse is constructed recursively on small, locally refined, temporary meshes. Conceptually, this corresponds to applying a fast direct solver [45] locally to regions with troublesome geometry. A global iterative method is then applied. Finally, the original layer density is reconstructed by running the recursion backwards, should it be explicitly needed. This gives us many of the advantages of fast direct methods, for example the ability to deal with certain classes of operators whose spectra make them unsuitable for iterative methods. In addition, the approach is typically much faster than using only a fast direct solver.
Our method, or scheme, has been referred to as recursive compressed inverse preconditioning [20,21,33,34] and there is a good reason for that name: the scheme relies on applying a relieving right inverse to the integral equation; on compressing this inverse to a low-dimensional subspace; and on carrying out the compression in a recursive manner. Still, the name recursive(ly) compressed inverse preconditioning is a bit awkward and we will here simply use the acronym RCIP.
A strong motivation for writing the present tutorial is that the original references [20,21,33,34] are hard to read. Certain derivations in [20,21,33,34] use complicated intermediary constructions, application specific issues obscure the general picture, and the notation has evolved from paper to paper. Here we focus on the method itself, on how it works and how it can be implemented, and refer to the original research papers for details. Demo programs in Matlab, updated as of December 2021, are a part of the exposition and can be downloaded from the web page: http://www.maths.lth.se/na/staff/helsing/Tutor/ Section 2 provides a historical background. Section 3 is a summary of the main features of RCIP. The basics of the method are then explained by solving a simple model problem in Sections 4-7. Sections 8-17 review general algorithmic improvements. Sections 18-21 contain applications to scattering problems. Sections 22-27 deal with close-to-touching objects, mixed (Zaremba) boundary conditions, Steklov eigenvalue problems, limit polarizability, vertex singularity exponents, and planar crack problems. Some of this material is new and has not been published elsewhere.

Background
The line of research on fast solvers for elliptic boundary value problems in piecewise smooth domains, leading up to the RCIP method, grew out of work in computational fracture mechanics. Early efforts concerned finding efficient integral equation formulations. Corner singularities were either resolved by brute force or by using special basis functions [19,27,37]. Such strategies, in combination with fast multipole [14] accelerated iterative solvers, work well for simple small-scale problems.
Real world physics is more complicated and, for example, the study [12] on a high-order time-stepping scheme for crack propagation (a series of biharmonic problems for an evolving piecewise smooth surface) shows that radically better methods are needed. Special basis functions are too complicated to construct and brute force is not economical -merely storing the discretized solution becomes too costly in a large-scale simulation.
A breakthrough came in 2007, when a scheme was created that resolves virtually any problem for Laplace's equation in piecewise smooth twodimensional domains in a way that is fully automatic, fast, stable, memory efficient, and whose computational cost scales linearly with the number of corners in the computational domain. The resulting paper [33] constitutes the origin of the RCIP method. Unfortunately, however, there are some flaws in [33]. For example, the expressions in [33,Section 9] are not generally valid and the paper fails to apply RCIP in its entirety to the biharmonic problem of [33,Section 3], which was the ultimate goal.
The second paper on RCIP [34] deals with elastic grains. The part [34,Appendix B], on speedup and enhanced stability, is particularly useful.
The third paper on RCIP [20] contains improvement relative to the earlier papers, both in the notation and in the discretization of singular operators. The overall theme is mixed boundary conditions, which pose similar difficulties as do piecewise smooth boundaries.
The fourth paper on RCIP [21], finally, solves the problem of [33,Section 3] in a broad setting, involving dominant integral operators with nonzero Fredholm indices and compositions of integral operators. In this context, too, some subsequent improvements have been made. See [28] and Sections 16,19.2,and 27,below. Further work on developing RCIP deal with more general boundary conditions [51], with problem-specific stabilization techniques [25], with singular right-hand sides [26], with problems in three dimensions [29,35,36], and with large-scale applications to aggregates of millions of grains [22,23].
We end this retrospection by noting that several research groups in recent years have proposed numerical schemes for integral equations stemming from elliptic partial differential equations (PDEs) in domains with boundary singularities. See, for example, [1,3,4,5,6,7,8,40,41,52,54]. There is also a widespread notion that a slight rounding of corners is a good idea for numerics. While rounding may work in particular situations, we do not believe it is a generally viable method. For one thing, how does one round a triple junction?

Summary of RCIP
This section summarizes Sections 4-14, below, and reviews the most important features of the RCIP method.
The starting point is an integral equation on a boundary Γ containing a corner (I + λK) ρ(r) = h(r) , r ∈ Γ .
Here I is the identity, λ is a parameter, K is an integral operator that is compact away from the corner, h(r) is a piecewise smooth right-hand side, and ρ(r) is an unknown layer density to be solved for.
Let the operator K be split into two parts where K describes the kernel interaction close to the corner and K • is a compact operator. Now introduce the transformed densitỹ ρ(r) = (I + λK ) ρ(r) .
Then use (2) and (3) to rewrite (1) as Although (4) looks similar to (1), there are advantages with using (4) from a numerical point of view. The RCIP method discretizes (4) chiefly on a grid on a coarse mesh on Γ that is sufficient to resolve K • and h(r). Only (I + λK ) −1 needs a grid on a locally refined fine mesh. Nyström discretization is used. The discretization of (4) assumes the form where R is a sparse block matrix called the compressed inverse. Note that (5) is a discrete system on the coarse grid only. The power of RCIP lies in the construction of R. In theory, R corresponds to a discretization of (I + λK ) −1 on the fine grid, followed by a lossless compression to the coarse grid. In practice, R is constructed via a forward recursion (29) where refinement and compression occur in tandem. The recursion starts on the smallest panels in a hierarchy of nested meshes around the corner, gradually moves up the hierarchy, and finally reaches the coarse mesh. At each refinement level a small matrix K • ib is needed as input and a small matrix R i is generated as output. The computational cost grows, at most, linearly with the number of refinement levels. Now, with access to the coarse-grid quantities R andρ coa only, surprisingly much is known about the solution ρ(r) on the fine grid. For example, the weight-corrected densityρ can be used to compute numerical approximations of integrals of ρ(r) against smooth functions f (r) as if they were carried out on the fine mesh Here w Γfin j and w Γcoa j are quadrature weights suitable for integrating polynomials on the fine grid and the coarse grid, respectively. With access also to the matrices K • ib and R i , everything is known about ρ(r) on the fine grid: The discrete density ρ fin can be reconstructed via a backward recursion (40); The eigenvalues of a certain backward recursion submatrix C contains information about the asymptotics of ρ(r) close to the corner vertex; The backward recursion acting on smooth basis functions automatically generates a tailor-made (singular) basis for ρ(r). Such a basis is helpful when RCIP is used for integral equations on non-smooth domains with edges in three dimensions.
It is important to observe that the forward recursion (29) is fast. It can be executed on the fly, even when the layer density ρ(r) is strongly singular and does not lie in any usual L p space. Deep inte the corner, the sequence of matrices K • ib have often converged to a beforehand given precision and (29) assumes the form of a fixed-point iteration. This opens up for the use of Newton's method. The computational cost for obtaining R can then be said to grow sub-linearly with respect to the number of refinement levels.
Let G(r, r ) be the fundamental solution to Laplace's equation in the plane: We shall solve the integral equation numerically for the unknown layer density ρ(r). Here ν is the exterior unit normal at r ∈ Γ, d is an element of arc length, λ is a parameter, e is a unit vector, and x y Figure 1: The contour Γ of (8) with a corner at the origin. The solid curve corresponds to opening angle θ = π/3. The dashed curve has θ = 4π/3.
The equation (10) models an electrostatic transmission problem [33] where e is an applied electric field. Using complex notation, where vectors r, r , ν, and e in the real plane R 2 correspond to points z, τ , n, and e in the complex plane C, one can write (10) as where the overbar symbol denotes the complex conjugate. Equation (12) is a simplification over (10) from a programming point of view. In many contexts it is advantageous to abbreviate (10) as where I is the identity. If Γ is smooth, then (13) is a Fredholm second kind integral equation with a compact, non-self-adjoint, integral operator K whose spectrum is discrete, bounded by one in modulus, and accumulates at zero. We also need a way to monitor the convergence of solutions ρ(r) to (13). For this purpose we introduce a quantity q, which corresponds to dipole moment or (un-normalized) polarizability [35] q Remark: Existence issues are important. Loosely speaking, the boundary value problem modeled by (10) has a unique finite-energy solution for a Γ ⋆ Figure 2: Left: A coarse mesh with ten panels on the contour Γ of (8) with opening angle θ = π/2. A subset of Γ, called Γ , covers four coarse panels as indicated by the dashed curve. Right: A fine mesh created from the coarse mesh by subdividing the panels closest to the corner n sub = 3 times.
large class of non-smooth Γ when λ is either off the real axis or when λ is real and λ ∈ [−1, 1). See [35] for sharper statements. The precise meaning of a numerical solution to an integral equation such as (10) also deserves comment. In this paper, a numerical solution refers to approximate values of ρ(r) at a discrete set of points r i ∈ Γ. The values ρ(r i ) should, in a postprocessor, enable the extraction of quantities of interest including values of ρ(r) at arbitrary points r ∈ Γ, functionals of ρ(r) such as q of (14), and the solution to the underlying boundary value problem at points in the domain where that problem was set.

Discretization on two meshes
We discretize (13) using standard Nyström discretization based on composite 16-point Gauss-Legendre quadrature on two different meshes: a coarse mesh with n pan quadrature panels and a fine mesh which is constructed from the coarse mesh by n sub times dyadically subdividing the panels closest to the corner in a direction toward the corner. The discretization is in parameter. The four panels on the coarse mesh that are closest to the corner should be equi-sized in parameter. These innermost four panels form a subset of Γ called Γ . See Figure 2. The linear systems resulting from the discretization on the coarse mesh and on the fine mesh can be written formally as where I and K are square matrices and ρ and g are column vectors. The subscripts fin and coa indicate what type of mesh is used. Discretization points on a mesh are said to constitute a grid. The coarse grid has n p = 16n pan points. The fine grid has n p = 16(n pan + 2n sub ) points.   (14) using (16) and the program demo1b.m (a loop of demo1.m) with lambda=0.999, theta=pi/2, npan=10, and evec=1. The reference value is q = 1.1300163213105365. There are n p = 160 + 32n sub unknowns in the main linear system. Left: Convergence with n sub . Right: The number of iterations needed to meet an estimated relative residual of mach .
The discretization of (13) is carried out by first rewriting (12) as whereτ (t) = dτ (t)/dt. Then Nyström discretization with n p points z i and weights w i on Γ gives The program demo1.m sets up the system (16), solves it using the GM-RES iterative solver [53] incorporating a low-threshold stagnation avoiding technique [32,Section 8], and computes q of (14). The user has to specify the opening angle θ, the parameter λ, the number n pan of coarse panels on Γ, the unit vector e and the number of subdivisions n sub . The opening angle should be in the interval π/3 ≤ θ ≤ 5π/3. We choose λ = 0.999, θ = π/2, n pan = 10, and e = (1, 0). The quantity q converges initially as n sub is increased, but for n sub > 44 the results start to get worse. See Figure 3. This is related to the fact, pointed out by Bremer [4], that standard Nyström discretization captures the L ∞ behavior of the solution ρ, while our ρ is unbounded. See, further, Appendix E.

Compressed inverse preconditioning
Let us split the matrices K coa and K fin of (15) and (16) into two parts each Here the superscript indicates that only entries of a matrix K ij whose indices i and j correspond to points z i and z j that both belong to the boundary subset Γ are retained. The remaining entries are zero. Now we introduce two diagonal matrices W coa and W fin which have the quadrature weights w i on the diagonal. Furthermore, we need a prolongation matrix P which interpolates functions known at points on the coarse grid to points on the fine grid. The construction of P relies on panelwise 15-degree polynomial interpolation in parameter using Vandermonde matrices. We also construct a weighted prolongation matrix P W via The matrices P and P W share the same sparsity pattern. They are rectangular matrices, similar to the identity matrix, but with one full (4 + 2n sub )16 × 64 block. Let superscript T denote the transpose. Then holds exactly. See Appendix A and [21,Section 4.3]. Equipped with P and P W we are ready to compress (16) on the fine grid to an equation essentially on the coarse grid. This compression is done without the loss of accuracy -the discretization error in the solution is unaffected and no information is lost. The compression relies on the variable substitution (I fin + λK fin ) ρ fin = Pρ coa .
Hereρ coa is the discretization of a piecewise smooth transformed density. The compression also uses the low-rank decomposition which should hold to about machine precision. The compressed version of (16) reads where the compressed weighted inverse R is given by  See Appendix B for details on the derivation. The compressed weighted inverse R, for Γ of (8), is a block diagonal matrix with one full 64 × 64 block and the remaining entries coinciding with those of the identity matrix. After having solved (25) forρ coa , the density ρ fin can easily be reconstructed fromρ coa in a post-processor, see Section 10. It is important to observe, however, that ρ fin is not always needed. For example, the quantity q of (14) can be computed directly fromρ coa . Let ζ coa be a column vector which contains values of |ż i | multiplied with {ēz i }. Then See Appendix C for a proof.

The recursion for R
The compressed weighted inverse R is costly to compute from its definition (26). As we saw in Section 5, the inversion of large matrices (I + K) on highly refined grids could also be unstable. Fortunately, the computation of R can be greatly sped up and stabilized via a recursion. In [33,Section 7.2] this recursion is derived in a roundabout way and uses a refined grid that differs from that of the present tutorial. A better derivation can be found in [21,Section 5], but there the setting is more general so that text could be hard to follow. Here we focus on results.

Basic prolongation matrices
Let P bc be a prolongation matrix, performing panelwise 15-degree polynomial interpolation in parameter from a 64-point grid on a four-panel mesh Figure 5: The boundary subsets Γ 3 , Γ 2 , and Γ 1 along with their corresponding type b meshes for n sub = 3.
to a 96-point grid on a six-panel mesh as shown in Figure 4. Let P W bc be a weighted prolongation matrix in the style of (21   Figure 5 for n sub = 3. Compare [20, Figure 2] and [21, Figure 5.1]. Let K ib denote the discretization of K on a type b mesh on Γ i . In the spirit of (19,20) we write where the superscript indicates that only entries with both indices corresponding to points on the four inner panels are retained.

The recursion proper
Now, let R n sub denote the full 64 × 64 diagonal block of R. The recursion for R n sub is derived in Appendix D and it reads where the operator F{·} expands its matrix argument by zero-padding (adding a frame of zeros of width 16 around it). Note that the initializer R 0 of (30) makes the recursion (29) take the first step The program demo2.m sets up the linear system (25), runs the recursion (29,30), and solves the linear system using the same techniques as demo1.m, see Section 5. In fact, the results produced by the two programs are very similar, at least up to n sub = 40. This supports the claim of Section 6 that the discretization error in the solution is unaffected by compression. Figure 6 demonstrates the power of RCIP: fewer unknowns and faster execution, better conditioning (the number of GMRES iterations does not grow), and higher achievable accuracy. Compare Figure 3. We emphasize that the number n sub of recursion steps (levels) used in (29) corresponds to the number of subdivisions n sub used to construct the fine mesh.

Schur-Banachiewicz speedup of the recursion
The recursion (29) can be sped up using the Schur-Banachiewicz inverse formula for partitioned matrices [38], which in this context can be written [34,  where A plays the role of R i−1 , P and P W are submatrices of P bc and P W bc , and U, V, and D refer to blocks of Besides, the integral equation (10) is replaced with which has the same solution ρ(r) but is more stable for λ close to one. For the discretization of (32) to fit the form (25), the last term on the left-hand side of (32) is added to the matrix λK • coa of (25). The execution of demo3.m is faster than that of demo2.m. Figure 7 shows that a couple of extra digits are gained by using (32) rather than (10) and that full machine accuracy is achieved for n sub > 60.

Various useful quantities
Let us introduce a new discrete densityρ coa viâ Rewriting (25) in terms ofρ coa gives which resembles the original equation (15). We see that K • coa , which is discretized using Gauss-Legendre quadrature, acts onρ coa . Therefore one  can interpretρ coa as pointwise values of the original density ρ(r), multiplied with weight corrections suitable for integration against polynomials. We refer toρ coa as a weight-corrected density. See, further, Appendix C.
Assume now that there is a square matrix S which mapsρ coa to discrete values ρ coa of the original density on the coarse grid The matrix S allows us to rewrite (25) as a system for the original density We can interpret the composition RS −1 as a matrix of multiplicative weight corrections that compensate for the singular behavior of ρ(r) on Γ when Gauss-Legendre quadrature is used. Let Y denote the rectangular matrix and let Q be a restriction operator which performs panelwise 15-degree polynomial interpolation in parameter from a grid on the fine mesh to a grid on a the coarse mesh. We see from (23) that Y is the mapping from ρ coa to ρ fin . Therefore the columns of Y can be interpreted as discrete basis functions for ρ(r). It holds by definition that The quantities and interpretations of this section come in handy in various situations, for example in 3D extensions of the RCIP method [35]. An efficient scheme for constructing S will be presented in Section 11.

Reconstruction of ρ fin fromρ coa
The action of Y onρ coa , which gives ρ fin , can be obtained by, in a sense, running the recursion (29) backwards. The process is described in detail in [20,Section 7]. Here we focus on results.
The backward recursion on Γ reads (40) Hereρ coa,i is a column vector with 64 elements. In particular,ρ coa,n sub is the restriction ofρ coa to Γ , whileρ coa,i are taken as elements {17 : 80} of ρ coa,i+1 for i < n sub . The elements {1 : 16} and {81 : 96} of ρ coa,i are the reconstructed values of ρ fin on the outermost panels of a type b mesh on Γ i . Outside of Γ , ρ fin coincides withρ coa .
When the recursion is completed, the reconstructed values of ρ fin on the four innermost panels are obtained from Should one wish to interrupt the recursion (40) prematurely, at step i = j say, then gives values of a weight-corrected density on the four innermost panels of a type b mesh on Γ j . That is, we have a part-way reconstructed weightcorrected densityρ part on a mesh that is n sub − j + 1 times refined. This observation is useful in the context of evaluating layer potentials close to their sources. If the memory permits, one can store the matrices K • ib and R i in the forward recursion (29) and reuse them in the backward recursion (40). Otherwise they may be computed afresh.
The program demo4.m builds on the program demo3.m, using (25) for (32). After the main linear system is solved forρ coa , a postprocessor reconstructs ρ fin via (40). Then a comparison is made with a solution ρ fin obtained by solving the un-compressed system (16). Figure 8 shows that for n sub < 10 the results are virtually identical. This verifies the correctness of (40). For n sub > 10 the result start to deviate. That illustrates the instabilities associated with solving (16) on a highly refined mesh. Compare Figure 3.
The program demo5.m investigates the effects of premature interruption of (40). The number of recursion steps is set to n sub = 100 and the recursion is interrupted at different levels. The density ρ fin is reconstructed on outer panels down to the level of interruption. Then a weight-corrected density is produced at the innermost four panels according to (42). Finally q of (14) is computed from this part-way reconstructed solution. The right image of Figure 8 shows that the quality of q is unaffected by the level of interruption.  (16) and ρ fin reconstructed fromρ coa of (25) via (40). Right: Relative accuracy in q of (14) from part-way reconstructed solutionŝ ρ part .

The construction of S
This section discusses the construction of S and other auxiliary matrices. Note that in many applications, these matrices are not needed.
The entries of the matrices P, P W , Q, R, S, and Y can only differ from those of the identity matrix when both indices correspond to discretization points on Γ . For example, the entries of R only differ from the identity matrix for the 64 × 64 block denoted R n sub in (29). In accordance with this notation we introduce P n sub , P W n sub , Q n sub , S n sub and Y n sub for the restriction of P, P W , Q, S and Y to Γ . In the codes of this section we often use this restricted type of matrices, leaving the identity part out.
We observe that S n sub is a square 64 × 64 matrix; P n sub , P W n sub and Y n sub are rectangular 16(4 + 2n sub ) × 64 matrices; and Q n sub is a rectangular 64 × 16(4 + 2n sub ) matrix. Furthermore, Q n sub is very sparse for large n sub . All columns of Q n sub with column indices corresponding to points on panels that result from more than eight subdivisions are identically zero.
The program demo6.m sets up P n sub , P W n sub and Q n sub , shows their sparsity patterns, and verifies the identities (22) and (38). The implementations for P n sub and P W n sub rely on repeated interpolation from coarser to finer intermediate grids. The implementation of Q n sub relies on keeping track of the relation between points on the original coarse and fine grids. Output from demo6.m is depicted in the left image of Figure 9. Note that the matrices P n sub and P W n sub are never needed in applications.
We are now ready to construct S. Section 10 presented a scheme for evaluating the action of Y n sub on discrete functions on the coarse grid on Γ . The matrix Y n sub , itself, can be constructed by applying this scheme to a 64 × 64 identity matrix. The matrix Q n sub was set up in demo6.m. Composing these two matrices gives S n sub , see (39). This is done in the program demo7.m, where the identity part is added as to get the entire matrix S. In previous work on RCIP we have found use for S in complex situations where (36) is preferable over (25), see [35,Section 9]. If one merely needs ρ coa fromρ coa in a post-processor, setting up S and using (35) is not worthwhile. It is cheaper to let Y act onρ coa and then let Q act on the resulting vector. Anyhow, demo7.m builds on demo4.m and gives as output ρ coa computed via (35), see the right image of Figure 9. For comparison, ρ fin , computed from (16), is also shown.
12 Initiating R using fixed-point iteration It often happens that Γ i is wedge-like. A corner of a polygon, for example, has wedge-like Γ i at all levels. If Γ is merely piecewise smooth, then the Γ i are wedge-like to double precision accuracy for n sub − i 1. Wedge-like sequences of Γ i open up for simplifications and speedup in the recursion (29,30). Particularly so if the kernel of the integral operator K of (13) is scale invariant on wedges. Then the matrix K • ib becomes independent of i. It can be denoted by K • b and needs only to be constructed once. Furthermore, the recursion (29,30) assumes the form of a fixed-point iteration The iteration (43) can be run until R i reaches its converged value R * . One need not know in advance how many iterations this takes. Choosing the number n sub of levels needed, in order to meet a beforehand given tolerance in R n sub , is otherwise a problem in connection with (29,30) and non-wedgelike Γ . This number has no general upper bound. Assume now that the kernel of K is scale invariant on wedges. If all Γ i are wedge-like, then (43,44) replaces (29,30) entirely. If Γ is merely piecewise smooth, then (43,44) can be run on a wedge with the same opening angle as Γ , to produce an initializer to (29). That initializer could often be more appropriate than R 0 of (30), which is plagued with a very large discretization error whenever (18) is used. The fixed-point initializer R * is implemented in the program demo8b.m, which is an upgrading of demo3b.m, and produces Figure 10. A comparison of Figure 10 with Figure 7 shows that the number n sub of levels needed for full convergence with the initializer R * is halved compared to when using the initializer R 0 of (30).
There are, generally speaking, several advantages with using the initializer R * , rather than R 0 of (30), in (29) on a non-wedge-like Γ : First, the number of different matrices R i and K • ib needed in (29) and in (40) is reduced as the recursions are shortened. This means savings in storage. Second, the number n sub of levels needed for full convergence in (29) seems to always be bounded. The hard work is done in (43). Third, Newton's method can be used to accelerate (43). That is the topic of Section 13.

Newton acceleration
When solving integral equations stemming from particularly challenging elliptic boundary value problems with solutions ρ(r) that are barely absolutely integrable, the fixed-point iteration (43,44) on wedge-like Γ may need a very large number of steps to reach full convergence. See [23, Section 6.3] for an example where 2 · 10 5 steps are needed.
Fortunately, (43) can be cast as a non-linear matrix equation where R * , as in Section 12, is the fixed-point solution and The non-linear equation (45), in turn, can be solved for R * with a variant of Newton's method. Let X be a matrix-valued perturbation of R * and expand G(R * + X) = 0 to first order in X. This gives a Sylvester-type matrix equation for the Newton update X. One can use the Matlab built-in function dlyap for (47), but GMRES seems to be more efficient and we use that method. Compare [23, Section 6.2]. Figure 11 shows a comparison between the fixed-point iteration (43,44) and Newton's method for computing the fixed-point solution R * to (45) on a wedge-like Γ . The program demo9.m is used and it incorporates Schur-Banachiewicz speedup in the style of Section 8. The wedge opening angle is θ = π/2, The integral operator K is the same as in (13), and λ = 0.999. The relative difference between the two converged solutions R * is 5 · 10 −16 . Figure 11 clearly shows that (43,44) converges linearly (in 68 iterations), while Newton's method has quadratic convergence. Only four iterations are needed. The computational cost per iteration is, of course, higher for Newton's method than for the fixed-point iteration, but it is the same at each step. Recall that the size of the underlying matrix I fin + λK fin , that is inverted according to (26), grows linearly with the number of steps needed in the fixed-point iteration. This example therefore demonstrates that one can invert and compress a linear system of the type (26) in sub-linear time.
14 The asymptotics of ρ(r) in the corner The recursion (40) provides a powerful tool for computing the asymptotics of ρ(r) close to the corner vertex: Deep in the corner, for large n sub and small i, the matrices R i−1 and K • ib can be replaced with their asymptotic counterparts R * and K • b , see Section 12, so that (40) reads where C is the constant 96 × 64 matrix Each step in (48) reconstructs ρ(r) on the outermost panels of a mesh on a subset Γ i that is half the size of the subset Γ i+1 in the previous step. Furthermore, the evolution ofρ coa,i is determined by power iteration applied to a submatrix of C given by row indices {17 : 80} and all columns. We denoted this submatrix by C . The eigenvalues of C control the behavior of ρ(r) as the distance to the corner vertex is halved. In particular, if γ is the arc length distance to the vertex then the leading asymptotic behavior where d max is the largest eigenvalue in modulus of C . For the opening angle θ = π/2 in the model problem of Section 4 it is possible to derive the closed-form expression [31,Eq. (13)] The program demo8c.m compares β computed from (50) to β from (51) with λ = 0.999. The relative difference is a mere 2 · 10 −15 , which means that RCIP provides a competitive alternative to traditional techniques, such as separation of variables [19,Section 2], also for asymptotic studies. The left image of Figure 12 shows the asymptotic behavior of ρ(γ) in the corner.
15 On the accuracy of "the solution" The integral equation (10) comes from a boundary value problem for Laplace's equation where the potential field U (r) at a point r in the plane is related to ρ(r) via see [33, Section 2.1]. The right image of Figure 12 shows how U (r) converges with mesh refinement at a point r = (0.4, 0.1) inside the contour Γ. We see that the accuracy in U (r) is slightly better than the accuracy of the dipole moment q of (14). One can say that measuring the field error at a point some distance away from the corner is more forgiving than measuring the dipole moment error. It is possible to construct examples where the difference in accuracy between field solutions and moments of layer densities are more pronounced and this raises the question of how the accuracy of integral equation solvers best should be measured.

Composed integral operators
Assume that we have a modification of (13) which reads Here K and g are as in (13), M is a new, bounded, integral operator, and ρ 1 is an unknown layer density to be solved for. This section shows how to apply RCIP to (53) using a simplified version of the scheme in [21]. Let us, temporarily, expand (53) into a system of equations by introducing a new layer density ρ 2 (r) = Kρ 1 (r). Then −Kρ 1 (r)+ ρ 2 (r) = 0 , and after discretization on the fine mesh Standard RCIP gives where the compressed inverse R is partitioned into four equi-sized blocks. Now we replaceρ 1coa andρ 2coa with a single unknownρ coa viã The change of variables (58,59) is chosen so that the second block-row of (57) is automatically satisfied. The first block-row of (57) becomes When (60) has been solved forρ coa , the weight-corrected version of the original density ρ 1 can be recovered aŝ (61) Figure 13 shows results for (53) with M being the double layer potential  Figure 13: Convergence for q of (14) with ρ = ρ 1 from (53). The curve Γ is as in (8) and theta=pi/2, npan=11, and evec=1. The reference values is taken as q = 1.95243329423584. Left: Results from the inner product preserving scheme of Appendix E produced with demo10.m. Right: Results with RCIP according to (60,61) produced with demo10b.m .
The left image shows the convergence of q of (14) with n sub using the inner product preserving discretization scheme of Appendix E for (53) as implemented in demo10.m. The right image shows q produced with RCIP according to (60,61) as implemented in demo10b.m. The reference value for q is computed with the program demo10c.m, which uses inner product preserving discretization together with compensated summation [39,42] in order to enhance the achievable accuracy. One can see that, in addition to being faster, RCIP gives and extra digit of accuracy. Actually, it seems as if the scheme in demo10.m converges to a q that is slightly wrong.
In conclusion, in this example and in terms of stability, the RCIP method is better than standard inner product preserving discretization and on par with inner product preserving discretization enhanced with compensated summation. In terms of computational economy and speed, RCIP greatly outperforms the two other schemes.

Nyström discretization of singular kernels
The Nyström scheme of Section 5 discretizes (13) using composite 16-point Gauss-Legendre quadrature. This works well when the kernel of the integral operator K is smooth on smooth Γ. When the kernel is not smooth on smooth Γ, then the quadrature fails and something better is needed. See [15] for a comparison of the performance of various modified high-order accurate Nyström discretizations for weakly singular kernels and [44] for a high-order general approach to the evaluation of layer potentials.
We are not sure what modified discretization is optimal in every situation. When nearly singular, weakly singular, and singular operators need to be discretized in the following, we chiefly use a modification to composite Gauss-Legendre quadrature called local panelwise evaluation. See [20, Section 2] and [24,30] for a description of this technique.

The exterior Dirichlet Helmholtz problem
Let D be the domain enclosed by the curve Γ and let E be the exterior to the closure of D. The exterior Dirichlet problem for the Helmholtz equation has a unique solution U (r) under mild assumptions on Γ and g(r) [49] and can be modeled using a combined integral representation [11, Chapter 3] where Φ ω (r, r ) is the fundamental solution to the Helmholtz equation in two dimensions Φ ω (r, r ) = i 4 H Here H 0 (·) is the zeroth order Hankel function of the first kind. Insertion of (66) into (64) gives the combined field integral equation where Figure 14 shows the performance of RCIP applied to (68) for 1000 different values of ω ∈ [1, 10 3 ]. The program demo11.m is used. This program has a fixed-point initializer R * , see Section 12, whose construction only takes the leading asymptotic behavior of I + K ω − iωS ω /2 at the corner vertex

The exterior Neumann Helmholtz problem
The exterior Neumann problem for the Helmholtz equation has a unique solution U (r) under mild assumptions on Γ and g(r) [49] and can be modeled as an integral equation in several ways. We shall consider two options: an "analogy with the standard approach for Laplace's equation", which is not necessarily uniquely solvable for all ω, and a "regularized combined field integral equation" which is always uniquely solvable. See, further, [3,9].

An analogy with the standard Laplace approach
Let K A ω be the adjoint to the double-layer integral operator K ω of (69) Insertion of the integral representation into (72) gives the integral equation Figure 15 shows the performance of RCIP applied to (76). The program demo12.m is used and the setup is the same as that for the Dirichlet problem in Section 18. A comparison between Figure 15 and Figure 14 shows that the number of GMRES iterations needed for full convergence now grows much faster with ω. Furthermore, the relative error in the solution to the Neumann problem is larger, particularly so when ω happens to be close to values for which the operator I − K A ω in (76) has a nontrivial null space. Recall that (68) is always uniquely solvable while (76) is not.   Figure 15, but RCIP is now applied to (81). The program demo13b.m is used.

A regularized combined field integral equation
The literature on regularized combined field integral equations for the exterior Neumann problem is rich and several formulations have been suggested. We shall use the representation [9] which after insertion into (72) gives the integral equation where The hypersingular operator T ω of (79) can be expressed as a sum of a simple operator and an operator that requires differentiation with respect to arc length only [46] T ω ρ(r) = 2ω 2 This makes it possible to write (78) in the form where A = −K A ω , B 2 = S iω , and the action of the operators B 1 , C 1 , and C 2 is given by All integral operators in (81) are such that their discretizations admit the low-rank decomposition (24). We use the temporary expansion technique of Section 16 for (81), with two new layer densities that are later eliminated, to arrive at a single compressed equation analogous to (60). That equation involves nine equi-sized blocks of the compressed inverse R. Solving the problem in the example of Section 19.1 again, we now take the number of panels on the coarse mesh as npan=0.6*omega+48 rounded to the nearest integer. Figure 16 shows results from the program demo13b.m. The resonances, visible in Figure 15, are now gone. It is interesting to observe in Figure 16 that, despite the presence of several singular operators and compositions in (81), the results produced with RCIP are essentially fully accurate and the number of GMRES iterations needed for convergence grows very slowly with ω.
The program demo13c.m differs from demo13b.m in that it uses local regularization [20, Section 2.1] for the Cauchy-singular operators of (83) and (84) rather than local panelwise evaluation. The results produced by the two programs are virtually identical. We do not show yet another figure.

Field evaluations
Strictly speaking, a boundary value problem is not properly solved until its solution can be accurately evaluated in the entire computational domain. The program demo11b.m is a continuation of demo11.m which, after solving (68) forρ coa with RCIP and formingρ coa via (33), computes the solution U (r) via (66) using three slightly different discretizations: (i) When r is away from Γ, 16-point Gauss-Legendre quadrature is used in (66) on all quadrature panels.
(ii) When r is close to Γ, but not close to a panel neighboring a corner, 16-point Gauss-Legendre quadrature is used in (66) on panels away from r and local panelwise evaluation is used for panels close to r.
(iii) When r is close to a panel neighboring a corner, the densityρ coa is first used to reconstructρ part according to Section 10. Then 16-point   Figure 17, but the exterior Neumann Helmholtz problem is solved using demo13d.m. The accuracy is even higher than in Figure 17. Gauss-Legendre quadrature is used in (66) on panels away from r and local panelwise evaluation is used for panels close to r.
The first two discretizations only use the coarse grid on Γ. The third discretization needs a grid on a partially refined mesh on Γ.
The program demo13d.m is a continuation of demo13b.m which, after solving (78) with RCIP as described in Section 19.2, computes the solution U (r) via (77) using the three discretizations of the previous paragraph. Figure 17 and 18 show that RCIP in conjunction with the quadrature of [20, Section 2] can produce very accurate solutions to exterior Helmholtz problems in, essentially, the entire computational domain.
The main source of error in the computed field U (r) of Figure 17 is cancellation in the evaluation of the difference r − r for r ∈ Γ. The program demo11b.m needs such differences in the discretized kernels of (69) and (70) and the vectors r and r are individually evaluated in global coordinates via (8). The program demo11e.m is the same as demo11b.m, but with r − r computed in local coordinates whenever r is close to r ∈ Γ. A comparison of Figure 19 with Figure 17 shows that the use of local coordinates on Γ lead to an improved quality in ρ(r) which, in turn, affects U (r). The improvement is most pronounced for U (r) close to Γ.
Further examples of Helmholtz problems in non-smooth exterior domains and more details on the discretization of Hankel kernels are found in [24,28].

A Helmholtz transmission problem
This section reviews some results from [30]. A transmission problem for the Helmholtz equation is formulated as where ε is a material parameter and ω 2 = √ εω 1 . We separate U (r) into an incident field U in (r) and a scattered field, represented by two layer densities µ and ρ and a uniqueness parameter c, so that for r ∈ E and for r ∈ D By this, the scattered field satisfies the outgoing radiation condition (65). Insertion of (89) and (90) into (87) and (88) gives the system of integral equations [43,Eq. (4.2)] with compact (differences of) operators and with We first apply RCIP to (91) for the purpose of computing eigenfields. The boundary Γ is as in (8) with θ = π/2. We set ε = 2.25, U in (r) = 0, and c = 1 (which is a common choice in the literature) and look for ω 1 , µ, ρ that are non-trivial solutions to the homogeneous system (91). Unfortunately, the system (91) admits false eigenwavenumbers, that is non-trivial solutions with {ω 1 } < 0 whose corresponding µ and ρ generate fields U (r) that vanish when inserted in (89) and (90). Nevertheless, a true eigenfield is found at ω 1 = 9.701129417644246 − 2.000374579086419i and shown in Figure 20 along with estimated field error. The coarse grid on Γ has 320 discretization points. The reference solution is computed with 50 per cent more points. The program demo19.m is used. We then compute the field {U (r)} in the limit of ε approaching the point ε = −1.1838 from above in the complex ε-plane. We set ω 1 = 18, c = −i, U in (r) = e iω 1 (r·d) with d = (cos(π/4), sin(π/4)), and use 800 discretization points on the coarse grid on Γ. The program, demo19b.m, is an extension of demo19.m: the construction of the initializer R * is accelerated using Newton's method, as described in Section 13, and a homotopy method is used for the limit {ε} → 0 + , see [23,Section 6.3]. Figure 21 shows results.

Close-to-touching objects
The usefulness of RCIP is not restricted to corner problems. RCIP works well also in more general contexts where solutions to integral equations exhibit some sort of (near) singularities. This section is about two such problems. First we compute the polarizability q of a pair of close-to-touching and highly conducting unit disks embedded in a background unit medium. Then we proceed to doubly periodic boundary conditions and compute the effective conductivity σ eff of a square array of conducting disks.

The two-disk problem
The setup is shown in Figure 22. This problem can be modeled with (10) and (14) and To avoid stability problems for λ close to one, we instead use an alternative formulation which in complex notation reads [33,Eqs. (9,10)]  Figure 22: Two unit disks with conductivity σ 2 > 1 are separated by a distance d and embedded in an infinite plane with conductivity σ 1 = 1. A unit field e is applied in the x-direction. Equi-sized quadrature panels are placed on the circles in such a way that there are breakpoints at r = (−d/2, 0) and r = (d/2, 0). type a type b type c Figure 23: Meshes of type a, type b, and type c on the boundary subset Γ i for the two-disk problem and for i = n sub = 3. The type a mesh has 8 + 4i panels. The type b mesh has twelve panels. The type c mesh has eight panels. The type a mesh is the restriction of the fine mesh to Γ i . For i = n sub , the type c mesh is the restriction of the coarse mesh to Γ . The type a mesh and the type b mesh coincide for i = 1.
RCIP can now be applied by considering r = (−d/2, 0) and r = (d/2, 0) to be singular boundary points treated in tandem. The subset Γ then covers the eight panels (four on each disk boundary) that are closest to the origin. Families of twelve-panel type b meshes are constructed in analogy with the procedure in Section 7.2. The superscript in K ib indicates that only entries with both indices corresponding to points on the eight innermost panels of a type b mesh are retained. The derivation of the recursion in Appendix D uses meshes of type a and type c with twice the number of panels compared to the single corner case. See Figure 23 and compare Figure 34.
The two-disk problem of Figure 22 is, in a sense, harder to solve than the one-corner model problem of Section 4. The reason being that the fine mesh on Γ for the two-disk problem has many panels that lie close to each other and where special quadratures techniques, see Section 17, need to be activated in the discretization of (95). This, in turn, slows down convergence and may even endanger the validity of the basic assumptions (24) and (D.3) upon which the entire RCIP scheme rests. The prolongation in (24) and (D.3) only holds on panels where standard quadrature is sufficient. In the one-corner model problem, on the other hand, special quadrature is barely needed. The basic assumptions (24) and (D.3) hold with, say, standard 16-point Gauss-Legendre quadrature provided that the opening angle θ is not too small.
The families of meshes introduced in Figure 23 are constructed with the validity of (24) and (D.3) in mind. When constructing K • ib in the recursion (29) and in (D.3), special quadrature may be activated in the discretization on meshes of type b if needed. Note, however, that the need for special quadrature will only arise for source points on the eight panels farthest away from the origin. For source points on the innermost four panels of type b meshes, standard quadrature is enough. The same is true for K • ia  Figure 24 illustrates the performance of RCIP applied to (95) and (96) for the two-disk problem. The program demo14.m is used. Convergence is immediate and it appears as if rather extreme cases can be treated accurately. Reference values for q can be found in demo14.m.
Remark: The two-disk problem was addressed in [33, Section 10.3], but not solved with RCIP in its entirety due to the too simplistic mesh construction technique used in [33].

The square array of disks
The left image of Figure 25 shows the geometry of this classic problem. The following modification of (95) and (96) is used for modeling where Γ per refers to all disk interfaces in the plane, Γ unit refers to the disk interface in the unit cell, see [18,Eqs. (13) and (14)], and σ eff is the effective conductivity. The applied electric field is chosen as e = (1, 0). The disk separation distance d of Figure 25 may be expressed in terms of a parameter c as The higher the number min{c, σ 2 /σ 1 }, the more difficult it is to compute σ eff via traditional numerical methods [48]. We solve (97,98) for three setups: σ 2 = 10 8 and c = 3 · 10 3 , which corresponds to the most extreme parameter choices in [10]; σ 2 = 10 4 and c = 10 4 , which is the hardest test case of [16, Table 2]; and σ 2 = 10 3 and c = 10 3 , which is used both in [16, Table 2] and [17, Table 1]. The right image of Figure 25 shows that the RCIP-accelerated Nyström solver demo14b.m resolves σ eff to full achievable accuracy already at 352 discretization points on the coarse grid on Γ unit and that the results are stable under mesh refinement. The number of converged digits compares favorably to what is reported in [10,16,17]. Reference values and a uniformly valid asymptotic expression [48] for σ eff are contained in demo14b.m.  (100) and (106). A coarse mesh is constructed on Γ. Two parts of the boundary, Γ 1 and Γ 2 , cover the four coarse panels closest to the singular boundary points γ 1 and γ 2 where the boundary conditions change type. The three sources S k of (107), for the generation of boundary conditions, are marked by green stars. Thousand target points in D are marked by tiny red dots.

Mixed boundary conditions
Elliptic PDEs with mixed boundary conditions, that is, Dirichlet conditions on parts of the boundary and Neumann conditions on the remaining contiguous parts (also known as Zaremba boundary conditions) can often be modeled using Fredholm second kind integral equations with operators that are smooth away from the points where the boundary conditions change type. In this context, too, RCIP improves the stability and greatly reduces the computational cost of Nyström discretization schemes.
The paper [20] shows how to apply RCIP to mixed planar harmonicand biharmonic problems. In this section we simply repeat two of the experiments in [20] for the purpose of disseminating the underlying Matlab programs (demo15.m and demo15b.m).
The interior mixed problem for Laplace's equation is solved on the domain D bounded by the contour Γ with the parameterization r(s) = (1 + 0.3 cos(5s))(cos(s), sin(s)) , −π ≤ s ≤ π . (100) We seek a function U (r), harmonic in D, such that where g D (r) is Dirichlet data on the boundary part Γ D , g N (r) is Neumann data on the boundary part Γ N , and Γ D ∪ Γ N = Γ. See Figure 26.
The solution U (r), r ∈ D ∪ Γ N , is represented by a density ρ(r), r ∈ Γ, (103) Insertion of (103) into (101) and (102) gives the system The boundary parts Γ D and Γ N are taken as r(s) ∈ Γ D , −π < s < − π 2 , and r(s) ∈ Γ N , − π 2 < s < π , (106) and the boundary conditions g D (r) and g N (r) are constructed from a closed form reference solution where S 1 = 1.4 + 1.4i, S 2 = −0.25 + 1.4i, and S 3 = −0.5 − 1.4i are sources outside of D, see Figure 26. Figure 27 illustrates the performance of RCIP applied to (104) and (105). The program demo15.m is used. The solution U (r) is evaluated via (103) at the 1000 target points in D indicated by red dots in Figure 26. The rapid convergence and high achievable accuracy seen in Figure 27 means that RCIP resolves the mixed problem very well.
The program demo15b.m is about reconstruction. It is a simplified version of a program used in [20]. Once the solutionρ coa is obtained, the discrete density ρ fin is reconstructed on the fine grid on Γ using (40) and (41). The program demo15b.m also constructs U (r) on the fine grid on Γ N , using [20,Eqs. (39), (40), and (49)], and then restricts U (r) to the coarse grid. Figure 28 shows results. The convergence and the achievable accuracy for U (r), r ∈ Γ N , is similar to that of U (r) with r some distance away from Γ. Compare the right image of Figure 28 with the left image of Figure 27.
Here the condition (109) on the boundary part Γ S is called a Steklov boundary condition and Γ = Γ D ∪ Γ S . Finding nontrivial harmonic solutions U (r) in D satisfying (108) and (109), along with associated values ς, is called a Steklov eigenvalue problem. We solve the mixed Steklov eigenvalue problem on the smooth domain D given by (100) using the same representation for U (r) as in (103) Insertion of (110) into (108) and (109) gives the homogeneous system The boundary parts Γ D and Γ S are taken from (100) as r(s) ∈ Γ D , −π < s < 0 , and r(s) ∈ Γ S , 0 < s < π , see the left image of Figure 29. Discretization of (111) and (112) together with RCIP leads to a linear system where R depends on ς, the matrix K • coa1 contains entries coming from the discretization of the integral operators in (111) and (112) that are not multiplied with ς, and the entries of K • coa2 come from the discretization of the remaining operators. Values of ς that correspond to a zero eigenvalue of the system matrix in (114) are solutions to the Steklov eigenvalue problem.
The right image of Figure 29, produced by the program demo16c.m, shows the three smallest system matrix eigenvalues of (114) as a function of ς. The program demo16d.m uses an eigenvalue search algorithm [29, Section 9.1] to generate a table of the first 50 Steklov eigenvalues. The estimated relative accuracy is about 10 −14 .

Pure Steklov eigenvalue problem on a square
Setups where Γ = Γ S , that is Γ D = ∅, and where Γ is only piecewise smooth are of particular interest in spectral theory. Recently some fascinating open problems have emerged [13]. The program demo16b.m computes the 20 first pure Steklov eigenvalues on the square D = (−1, 1) × (−1, 1) using a RCIPaccelerated solver very similar to that of demo16d. The numerical results are compared with results from the semi-analytic expressions of [13, Section 3.1]. The estimated relative accuracy is on the order of mach .

Limit polarizability
Let us return to (13) and write it in the form  Figure 30: Elements α + 11 (u) and α + 22 (u) of the limit polarizability tensor for the object enclosed by Γ of (8) and with θ = π/2. The program demo17.m is used.
is a new complex variable. Values of w for which (115) has no solution are points in the spectrum of K. The precise nature of this spectrum depends both on Γ and on the function space considered [35,36]. On the "energy space" H −1/2 (Γ), the spectrum of K is real and may have both discrete and continuous parts. We solve (115) and compute the normalized polarizability where |V | is the area enclosed by Γ. We are particularly interested in α + (u), that is, the limit of α(w) as v → 0 + . The programs used are extensions of demo8b.m. The construction of the initializer R * is accelerated using Newton's method and homotopy, compare Section 21.
The program demo17.m computes α + (u) for Γ as in (8) and with θ = π/2. The applied electric field is either e = (1, 0), giving the element α + 11 (u) of the limit polarizability tensor, or e = (0, 1), giving α + 22 (u). Figure 30 shows results. By varying θ in demo17.m, one can see that a continuous non-zero {α + (u)} is only possible in the interval −|1 − θ/π| < u < |1 − θ/π|. The program demo17b.m computes α + (u) for Γ being the unit square, compare [35, Figure 5(a)] where a similar Matlab program is used. In addition, demo17b.m also computes the singularity exponent β in the leading asymptotic behavior of ρ + (r) in the square corners, where γ is the arc length distance to the nearest corner vertex, see Section 14. Figure 31 shows results. We emphasize that for v = 0 and −0.5 < u < 0.5, there is no solution ρ(r) to (115). There is, however, a solution ρ(r) ∈ H −1/2 (Γ) for v arbitrarily close to 0 and it is the polarizability α(w) corresponding to this limit solution, with v → 0 + , that is depicted in the left image of Figure 31. Furthermore, the general polarizability α(w) is simply related to the limit polarizability {α + (u)} via see [35,Section 3]. One can say that {α + (u)}/π is the derivative of a spectral measure associated with α(w).

Some computations on the cube
The applicability of RCIP acceleration to Nyström discretization of Fredholm second kind integral equations is not restricted to planar problems. It extends also to 3D. Rotationally symmetric surfaces that are smooth aside from isolated sharp edges or conical points are particularly simple to deal with [29,36]. Surfaces that contain a mix of contiguous edges and corners require that RCIP is applied in a two-step manner: First it is used to find multiplicative weight corrections that capture the singular behavior of ρ(r) in directions perpendicular to the edges using the techniques of Section 9.
With these corrections incorporated into the standard quadrature, ρ(r) can be resolved on the coarse grid except for in the corners. There intense refinement has to take place and RCIP is used a second time.
Reference [35] details the two-step procedure for the solution of (115) on the surface of the unit cube. The action of the integral operator K in 3D is The left image of Figure 32, taken from [35], shows the limit polarizability α + (u) of the cube computed in this way. The right image of Figure 32, produced by the program demo18.m, shows the vertex singularity exponent β(u) of a cube corner. The function β(u) is interesting since for some of its arguments there exist a number of benchmarks. For example, the quantity 1 + β(−1) is a so-called "Fichera-type eigenvalue" for which the values 0.45417371 [56] and 0.454173734 [47] have been reported. Our estimate, produced by an upgraded version of demo18.m, is 0.45417373430 (14). The two digits within parenthesis are extrapolated.

Branched cracks in an elastic plane
As mentioned in Section 2, the RCIP method grew out of work in computational fracture mechanics. A particular goal was to find an efficient way to compute the so called normalized stress intensity factors at the tips γ ep 1 and γ ep 2 of a V-shaped crack in an elastic plane. This biharmonic boundary value problem can be modeled as a Fredholm second kind integral equation with composed operators in the form (53). The stress intensity factors are simple functionals of the layer density. See [21] for details. One purpose of the present section is to disseminate a Matlab program, demo20.m, that reproduces [21, Figure 7.1(a)]. This figure shows how the stress intensity factors converge with the number of levels n sub in the recursion for the compressed inverse R. The treatment of composed integral operators in [21] has proven to be unnecessarily complicated and is not used in demo20.m. Instead, demo20.m relies on the simplified expansion technique for (53) that was described in Section 16.
Introducing the fundamental function [50, Section 107] the actions of the operators corresponding to K and M of (53) can, for the class of crack problem studied in [21, Section 6], be expressed as and The right-hand side in (53) is Figure 33, top row, shows the V-shaped crack and output from demo20.m. The bottom row of Figure 33 shows convergence of the stress intensity factor at the tip γ ep 2 of the symmetrically branched crack. For this crack, the fundamental function assumes the more complicated form but the simplified expansion technique for (53), described in Section 16, still applies. The results in the bottom row of Figure 33 are obtained with demo20c.m and agree with those in [21, Table 7.2]. The program demo20b.m produces the same results (not shown), but uses the more involved original treatment of composed integral operators from [21].

RCIP in a Method-of-Moments context
The RCIP method is developed to accelerate and stabilize Nyström solvers in the presence of boundary singularities. While Nyström schemes are efficient, they are not the most common in, for example, computational electromagnetics, where Method-of-Moments (MoM) solvers dominate. A fair question to ask is therefore: does RCIP apply also to the MoM? The answer is "Yes". Roughly speaking, as we now show, the MoM amounts to a similarity transformation of the linear system resulting from Nyström discretization. RCIP still applies, unaffected by this. Let us return to (15) and (16) which, after the change of variables can be written in MoM-form as Here c is a column vector with n p entries (interpreted as coefficients) and L is a block diagonal matrix with identical blocks L j , j = 1, . . . , n pan and j = 1, . . . , n pan + 2n sub for (127) and (128), respectively. The L j is of size 16 × 16 and has the 16 first standard Legendre polynomials P n (x), n = 0, . . . , 15, evaluated at the nodes T16, see Section 7.1, as columns (interpreted as basis functions). Note that the inverse L −1 is simple to set up accurately since where N is a diagonal matrix of normalization constants (2n + 1)/2 , n = 0, . . . , 15 , and W is a diagonal matrix containing the quadrature weights W16. Now, repeating the derivation steps of Section 6, we arrive at (25) which assumes the form where andc coa is related to the weight-corrected densityρ coa of (33) viâ The recursion (29,30) assumes the form F{R −1 Here where L bb is L on a type b mesh, the matrix L bc maps coefficients on a type c mesh to values of Legendre polynomials on a grid on a type b mesh, the diagonal matrix W b can be constructed as and N c is N on a type c mesh. The program demo21.m implements an RCIP-accelerated MoM solver for the first numerical example in Section 20. Simply put: demo21.m is a conversion of demo11.b to the MoM framework, as outlined above. In particular, the Nyström matrix K has been replaced with the MoM matrix K. The results produced by demo11.b and demo21.m are almost identical, although demo11.b of course has a faster setup phase than demo21.m. Less work is required to obtain the entries of K than those of K. *** Appendicies *** Appendix A. Proof that P T W P = I coa Let f coa and g coa be two column vectors, corresponding to the discretization of two panelwise polynomials with panelwise degree 15 on the coarse mesh of Γ. Then because composite 16-point Gauss-Legendre quadrature has panelwise polynomial degree 31. The diagonal matrix W coa has size 16n pan × 16n pan .
Since there are 16n pan linearly independent choices of f coa and of g coa it follows from (A.1) that which, using (21), can be rewritten

Appendix B. Derivation of the compressed equation
The compression of (16), leading up to (25), was originally described in [33,Section 6.4]. Here we give a summary. The starting point is (13) which, using the operator split analogous to (19,20) K and the variable substitution gives the right preconditioned equatioñ Now, let us take a close look at (B.3). We observe that K • (I + λK ) −1 is an operator whose action on any function gives a function that is smooth on the innermost two panels of the coarse mesh on Γ . This is so since K • is constructed so that its action on any function gives a function that is smooth on the innermost two panels of the coarse mesh on Γ . Furthermore, the right-hand side λg(r) of (B.3) is assumed to be panelwise smooth on the coarse mesh. Using an argument of contradiction we see thatρ(r) has to be panelwise smooth on the innermost two panels of the coarse mesh on Γ .

Appendix D. Derivation of the recursion
The recursion (29) for the rapid construction of the diagonal blocks of the compressed weighted inverse R was originally derived in [33, Section 7] using different notation and different meshes than in the present tutorial. The recursion was derived a second time in [34,Section 7] using new meshes. Better notation was introduced in [20, Section 6]. A third derivation, in a general setting, takes place in [21,Section 5] and it uses the same notation and meshes as in the present tutorial.
A problem when explaining the derivation of (29) is that one needs to introduce intermediate meshes and matrices whose appearance may cause enervation at a first glance. Particularly so since these meshes and matrices are not needed in the final expression (29). We emphasize that the underlying matrix property that permits the recursion is the low rank of certain off-diagonal blocks in discretizations of K • of (B.1) on nested meshes. type a type b type c Figure 34: Meshes of type a, type b, and type c on the boundary subset Γ i for i = n sub = 3. The type a mesh has 4 + 2i panels. The type b mesh has six panels. The type c mesh has four panels. The type a mesh is the restriction of the fine mesh to Γ i . For i = n sub , the type c mesh is the restriction of the coarse mesh to Γ . The type a mesh and the type b mesh coincide for i = 1.
The recursion (29) only uses uses one type of mesh explicitly -the type b mesh of Figure 5. On each Γ i there is a type b mesh and a corresponding discretization of K • denoted K • ib . Here we need two new types of meshes denoted type a and type c, along with corresponding discrete operators. For example, K ia is the discretization of K on a type a mesh on Γ i . The three types of meshes are depicted in Figure 34. Actually, a straight type c mesh was already introduced in Figure 4. Now we define R i as where P W iac and P iac are prolongation operators (in parameter) from a grid on a type c mesh on Γ i to a grid on a type a mesh on Γ i . Note that R i for i = n sub , according to the definition (D.1), is identical to the full diagonal 64 × 64 block of R of (26). Note also that R 1 comes cheaply. The rest of this appendix is about finding an expression for R i in terms of R i−1 that is cheap to compute. Let us split K ia into two parts where K ia = F{K (i−1)a } and K • ia is such that holds to about machine precision, compare (24). The prolongation operators P iab and P W iab act from a grid on a type b mesh to a grid on a type a mesh. It holds that P iac = P iab P bc , (D.4) Summing up, we can rewrite (D.1) as (D.6) The subsequent steps in the derivation of (29) are to expand the inverse of the sum of matrices within parentheses in (D.6) using a Taylor series where A corresponds to the first two terms and B corresponds to the last term; multiply the terms in this series with P T W iab from the left and with P iab from the right; and bring the series back in closed form. The result is

Appendix E. An inner product preserving scheme
In [4], Bremer describes a scheme that stabilizes the solution to the discretized system (16) on the fine mesh. The scheme can be interpreted as an inner product preserving discretization. In practice it corresponds to making a similarity transformation of the system matrix. While inner product preserving Nyström discretization elegantly solves problems related to stability (the condition number of the system matrix is improved) it does not reduce the number of discretization points (unknowns) needed to achieve a given precision in the solution. Neither does it affect the spectrum of the system matrix (similarity transformations preserve eigenvalues) and hence it does not in any substantial way improve the convergence rate of the GMRES iterative method [55, Lecture 35].
For completeness, we have implemented inner product preserving Nyström discretization in the program demo1d.m. The program is a continuation of demo1b.m where we also have replaced (10) with the more stable integral equation (32). This should facilitate comparison with the program demo3b.m and the results shown in Figure 7.   Figure 35 shows results produced by demo1d.m. Beyond n sub = 60 one now achieves essentially full machine precision in q of (14). Despite this success, inner product preserving Nyström discretization can perhaps not quite compete with the RCIP method in this example. The differences in performance relate to issues of memory and speed. The RCIP method uses a much smaller linear system (16n pan unknowns) than does inner product preserving Nyström discretization (16(n pan + 2n sub ) unknowns). Besides, the RCIP method converges in only eight GMRES iterations, irrespective of n sub . See Figure 7.