1% Two methods on the same tridiagonal solver:
2%
3% 1. the three-point scheme of the previous chapter;
4% 2. the same scheme run at N and at 2N, combined by Richardson
5% extrapolation to cancel the leading h^2 term of the error.
6%
7% If the error has an asymptotic expansion u_h = u + c h^2 + O(h^4), then
8%
9% (4 u_{h/2} - u_h) / 3 = u + O(h^4)
10%
11% at the points the two meshes share. Nothing about the discretisation
12% changes; the fourth-order accuracy comes out of two second-order solves and
13% one subtraction. The expansion is what makes it work, and the expansion
14% needs u to be smooth: on the low-regularity solution there is no h^2 term to
15% cancel and the extrapolation gains nothing.
17function names = method_names()
18names = {'3-point differences', 'Richardson from N and 2N'};
19end
21function [x, u] = solve1d(mid, N, f, ua, ub)
22if mid == 1
23 [x, u] = fd2(N, f, ua, ub);
24else
25 [x, uc] = fd2(N, f, ua, ub);
26 [~, uf] = fd2(2*N, f, ua, ub);
27 u = (4 * uf(1:2:end) - uc) / 3;
28end
29end
31function [x, u] = fd2(N, f, ua, ub)
32h = 1 / N;
33x = (0:N)' / N;
34r = h^2 * f(x(2:N));
35r(1) = r(1) + ua;
36r(end) = r(end) + ub;
37u = [ua; thomas(N-1, r); ub];
38end
40function w = thomas(n, d)
41c = zeros(n, 1);
42e = zeros(n, 1);
43b = 2;
44c(1) = -1 / b;
45e(1) = d(1) / b;
46for i = 2:n
47 b = 2 + c(i-1);
48 c(i) = -1 / b;
49 e(i) = (d(i) + e(i-1)) / b;
50end
51w = zeros(n, 1);
52w(n) = e(n);
53for i = n-1:-1:1
54 w(i) = e(i) - c(i) * w(i+1);
55end
56end