Every method so far has been a tensor product on a square, and every convergence rate has come from a Taylor expansion of a smooth solution. This chapter removes one quarter of the square and watches both assumptions fail.

The domain is the L,

with a re-entrant corner at the origin of interior angle $3\pi/2$. It is the standard test problem for corner singularities, and its behaviour is representative of any domain with a non-convex corner.

What the corner does#

Near a corner of interior angle $\alpha$, harmonic functions with zero data on both edges behave like $r^{\pi/\alpha}$. For $\alpha = 3\pi/2$ the exponent is $2/3$, so the local behaviour is spanned by

measuring $\theta$ from the edge $y = 0$, $x > 0$. Each $s_k$ is harmonic and vanishes on both edges of the corner. The first one, $s_1 = r^{2/3}\sin(2\theta/3)$, is continuous but has an infinite gradient at the origin and no second derivative there.

This matters because a generic Dirichlet datum on $\partial\Omega$ produces a solution containing $s_1$. Not approximately: the solution genuinely has that term, with some coefficient determined non-locally by the data on the whole boundary. The finite difference scheme is then being asked to approximate a function that has no fourth derivative, no third, and not even a bounded first derivative, and the argument of Chapter 1 collapses.

Three ground truths are available in the labs, so that the effect can be isolated:

• sing, the singular function $s_1$ alone;
• smooth, $\log|z - z_0|$ with $z_0$ outside, which is harmonic and smooth right up to and including the corner;
• both, the sum, which is the generic case.

Look at the error picture on both with the difference method. The error is not spread over the domain; it is concentrated in a small neighbourhood of the origin, and refining the grid shrinks that neighbourhood without reducing the error in it very much. The smooth case, on the same domain and the same grid, has a smooth error of the usual shape. The corner is not the problem; the singular function is.

The rate#

Three numbers from the sweep, on the maximum error over interior points at a distance of at least $1/8$ from the corner:

ground truth	observed order of the difference scheme
smooth	2.00
both	about 1.36
sing	about 1.40

The classical result for the five-point scheme on this domain is $O(h^{4/3})$ away from the corner, which is what the last two rows show. Right at the corner the error is worse, $O(h^{2/3})$. Either way, second order is gone, and no solver recovers it: the discretisation, not the linear algebra, is the limit.

The standard remedies are worth naming even though the labs here take a different route. Local refinement with a grid graded like $h_j \sim j^{-3/2}$ towards the corner restores $O(h^2)$ at the cost of an unstructured or hierarchical grid, and is what an adaptive finite element code does automatically. Singular enrichment adds $s_1$ to the finite element space directly. Both are more work than the method below, on this problem.

Putting the singularity in the basis#

If the solution is a smooth harmonic function plus a few multiples of $s_k$, the obvious thing to do is to build an approximation out of exactly that. Take

every term of which is harmonic, so any combination satisfies the equation exactly in the interior. Nothing is left but the boundary condition, and that is imposed by least squares at a few hundred boundary points. This is the method of particular solutions, in the form Fox, Henrici and Moler used for the L-shaped eigenvalue problem in 1967, and it is the natural method for a problem like this one.

The results in the sweep are of a different kind from anything else in this tutorial. On sing, the exact solution is the first basis function, so the method is exact to roundoff with eight terms. On both it converges geometrically, reaching $1.6\times10^{-14}$ with 121 basis functions in about ten milliseconds, while the difference scheme with 49000 unknowns and four seconds of conjugate gradients is at $4\times10^{-4}$. That is ten orders of magnitude and a factor of four hundred in time, in the same direction.

Three honest qualifications, because a gap that size always has conditions attached.

• It needs the equation to be homogeneous and linear with known solutions. Every basis function must be harmonic. With a source term, a variable coefficient, or a nonlinearity, this method does not exist in this form.
• It needs to know where the singularities are. The exponents $2k/3$ come from the corner angle. On a domain with several re-entrant corners each one contributes its own family, which is manageable; on a problem where the singular behaviour is not known analytically it is not.
• The least-squares matrix is ill-conditioned, with a condition number that grows quickly in the number of terms. It is usable to about $10^{-14}$ here because the boundary data is smooth and the system is solved by QR; pushed further it stalls, and the stall is not visible in the residual of the fit. Adding terms past the point where the error stops improving makes things slowly worse.