1Every 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.
3The domain is the L,
5$$\Omega = (-1,1)^2 \setminus \{x > 0,\ y < 0\},$$
7with 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.
9### What the corner does
11Near 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
13$$s_k(r,\theta) = r^{2k/3}\sin\frac{2k\theta}{3}, \qquad k = 1, 2, 3, \dots$$
15measuring $\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.
17<!--fig:reentrant-corner-->
19This 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.
21Three ground truths are available in the labs, so that the effect can be isolated:
23- `sing`, the singular function $s_1$ alone;
24- `smooth`, $\log|z - z_0|$ with $z_0$ outside, which is harmonic and smooth right up to and including the corner;
25- `both`, the sum, which is the generic case.
27<!--lab:field-->
29Look 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.
31### The rate
33<!--lab:sweep-->
35Three numbers from the sweep, on the maximum error over interior points at a distance of at least $1/8$ from the corner:
37| ground truth | observed order of the difference scheme |
38|---|---|
39| `smooth` | 2.00 |
40| `both` | about 1.36 |
41| `sing` | about 1.40 |
43The 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.
45The 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.
47### Putting the singularity in the basis
49If 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
51$$u \;\approx\; \sum_{m=0}^{M} \big(a_m \operatorname{Re} z^m + b_m \operatorname{Im} z^m\big) \;+\; \sum_{k} c_k\, r^{2k/3}\sin\frac{2k\theta}{3},$$
53every 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.
55The 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.
57Three honest qualifications, because a gap that size always has conditions attached.
59- **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.
60- **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.
61- **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.