1function [u, f] = lp_exact1d(kind, x)
2% LP_EXACT1D A manufactured solution of -u'' = f on [0, 1], and its source.
3%
4% Pick u first, differentiate it by hand, and call the result f. Then u is by
5% construction the exact solution of the boundary value problem
6%
7% -u''(x) = f(x) on (0, 1), u(0), u(1) given by u itself,
8%
9% and any difference between a computed solution and u is error, known to full
10% precision at every point and at no cost. The five choices below differ only
11% in how hard they are to resolve.
12switch kind
13 case 'smooth'
14 % entire, moderate amplitude: the easy case
15 s = sin(3*pi*x);
16 u = exp(s);
17 f = (3*pi)^2 * (s - cos(3*pi*x).^2) .* u;
18 case 'oscil'
19 % eight wavelengths across the interval: still entire, but a mesh must
20 % resolve the wavelength before any convergence rate appears
21 u = sin(8*pi*x) + x;
22 f = (8*pi)^2 * sin(8*pi*x);
23 case 'peak'
24 % a Gaussian of width 1/20 at x = 0.4: entire, but with all of its
25 % structure in a tenth of the interval
26 a = 200;
27 u = exp(-a*(x - 0.4).^2);
28 f = (2*a - 4*a^2*(x - 0.4).^2) .* u;
29 case 'layer'
30 % a boundary layer of width eps at x = 1, as in a singular perturbation
31 ep = 0.02;
32 u = exp(-(1 - x)/ep);
33 f = -u/ep^2;
34 case 'lowreg'
35 % u is C^2 but not C^3: |x - 1/2|^(5/2) has a square root in u''.
36 % Every polynomial-based method feels this.
37 c = 0.5;
38 u = abs(x - c).^(5/2);
39 f = -(15/4) * sqrt(abs(x - c));
40 otherwise
41 error('lp_exact1d: unknown solution "%s"', kind);
42end
43end