088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 1% Three methods on the same problem: the two finite difference schemes of the
2% previous chapters, and Chebyshev collocation.
3%
4% Collocation drops the idea of a local stencil. Put the N+1 Chebyshev points
5%
6% x_j = ( 1 - cos(pi j / N) ) / 2, j = 0, ..., N,
7%
8% on [0,1], let p be the polynomial of degree N through the unknown values, and
9% require -p'' = f at the interior points. Differentiation of a polynomial
10% given by its values at fixed nodes is a matrix, so the discrete operator is
11% the dense matrix D^2 rather than a three-term recurrence. Every unknown is
12% coupled to every other, the solve costs O(N^3) instead of O(N), and in
13% exchange the error decays faster than any power of N for analytic solutions.
14%
15% The clustering of the nodes near the ends is not decoration. It is what makes
16% polynomial interpolation of high degree stable, and equally spaced points of
17% the same degree would be useless.
19function names = method_names()
20names = {'3-point differences', 'compact 4th order', 'Chebyshev collocation'};
21end
23function [x, u] = solve1d(mid, N, f, ua, ub)
24if mid == 3
25 [x, u] = cheb_solve(N, f, ua, ub);
26else
27 [x, u] = fd_solve(mid, N, f, ua, ub);
28end
29end
31function [x, u] = cheb_solve(N, f, ua, ub)
32[D, s] = lp_cheb(N);
33x = (1 - s) / 2; % Chebyshev points, ascending, on [0, 1]
34Dx = -2 * D; % d/dx after the affine map
35L = -(Dx * Dx); % the operator -d^2/dx^2, as a dense matrix
37ii = 2:N; % interior collocation points
38A = L(ii, ii);
39b = f(x(ii)) - L(ii, 1) * ua - L(ii, N+1) * ub;
40u = [ua; A \ b; ub];
41end
43function [x, u] = fd_solve(mid, N, f, ua, ub)
44h = 1 / N;
45x = (0:N)' / N;
46fx = f(x);
47if mid == 1
48 r = h^2 * fx(2:N);
49else
50 r = h^2 * (fx(1:N-1) + 10*fx(2:N) + fx(3:N+1)) / 12;
51end
52r(1) = r(1) + ua;
53r(end) = r(end) + ub;
54u = [ua; thomas(N-1, r); ub];
55end
57function w = thomas(n, d)
58c = zeros(n, 1);
59e = zeros(n, 1);
60b = 2;
61c(1) = -1 / b;
62e(1) = d(1) / b;
63for i = 2:n
64 b = 2 + c(i-1);
65 c(i) = -1 / b;
66 e(i) = (d(i) + e(i-1)) / b;
67end
68w = zeros(n, 1);
69w(n) = e(n);
70for i = n-1:-1:1
71 w(i) = e(i) - c(i) * w(i+1);
72end
73end