1% Two discretisations that share a matrix.
2%
3% 1. the three-point scheme: (u_{i-1} - 2u_i + u_{i+1}) / h^2 = -f_i
4% 2. the compact fourth-order scheme, due to Numerov:
5%
6% (u_{i-1} - 2u_i + u_{i+1}) / h^2 = -(f_{i-1} + 10 f_i + f_{i+1}) / 12
7%
8% The left-hand sides are identical. All that changes is that the source is
9% averaged over the same three points with weights 1/12, 10/12, 1/12, which is
10% exactly the combination that cancels the h^2 term of the truncation error.
11% The matrix is still tridiag(-1, 2, -1), the solve is still O(N), and the
12% method is fourth-order accurate: the extra two orders cost one weighted
13% average of numbers we already had.
14%
15% The catch is that the scheme wants f at the nodes and nothing else, so it
16% needs f to be smooth as well: the cancellation is a statement about the
17% Taylor expansion of f, not just of u.
19function names = method_names()
20names = {'3-point differences', 'compact 4th order (Numerov)'};
21end
23function [x, u] = solve1d(mid, N, f, ua, ub)
24h = 1 / N;
25x = (0:N)' / N;
26fx = f(x);
28if mid == 1
29 r = h^2 * fx(2:N);
30else
31 r = h^2 * (fx(1:N-1) + 10*fx(2:N) + fx(3:N+1)) / 12;
32end
33r(1) = r(1) + ua;
34r(end) = r(end) + ub;
36u = [ua; thomas(N-1, r); ub];
37end
39function w = thomas(n, d)
40c = zeros(n, 1);
41e = zeros(n, 1);
42b = 2;
43c(1) = -1 / b;
44e(1) = d(1) / b;
45for i = 2:n
46 b = 2 + c(i-1);
47 c(i) = -1 / b;
48 e(i) = (d(i) + e(i-1)) / b;
49end
50w = zeros(n, 1);
51w(n) = e(n);
52for i = n-1:-1:1
53 w(i) = e(i) - c(i) * w(i+1);
54end
55end