1% Second-order finite differences for -u'' = f on (0, 1), u(0) = ua, u(1) = ub.
2%
3% Replace u''(x_i) by the three-point second difference
4%
5% u''(x_i) ~ ( u(i-1) - 2 u(i) + u(i+1) ) / h^2 ,
6%
7% which is exact for quadratics and wrong by (h^2/12) u''''(x_i) for anything
8% else. Collecting the N-1 interior equations gives a tridiagonal system whose
9% matrix is (1/h^2) * tridiag(-1, 2, -1); the two known boundary values move to
10% the right-hand side.
12function names = method_names()
13names = {'3-point differences'};
14end
16function [x, u] = solve1d(mid, N, f, ua, ub)
17h = 1 / N;
18x = (0:N)' / N;
20r = h^2 * f(x(2:N)); % scale by h^2 so the matrix is tridiag(-1, 2, -1)
21r(1) = r(1) + ua; % the neighbour at x = 0 is known
22r(end) = r(end) + ub; % and so is the one at x = 1
24u = [ua; thomas(N-1, r); ub];
25end
27function w = thomas(n, d)
28% Solve tridiag(-1, 2, -1) w = d by Gaussian elimination without pivoting,
29% which is safe here because the matrix is symmetric positive definite and
30% diagonally dominant. Two sweeps, O(n) work and O(n) storage: the whole
31% reason a one-dimensional problem is cheap at any resolution.
32c = zeros(n, 1);
33e = zeros(n, 1);
34b = 2;
35c(1) = -1 / b;
36e(1) = d(1) / b;
37for i = 2:n
38 b = 2 + c(i-1); % = 2 - (-1)*c(i-1)
39 c(i) = -1 / b;
40 e(i) = (d(i) + e(i-1)) / b;
41end
42w = zeros(n, 1);
43w(n) = e(n);
44for i = n-1:-1:1
45 w(i) = e(i) - c(i) * w(i+1);
46end
47end