1% Three discretisations of the same Dirichlet problem.
2%
3% 1. the five-point stencil, solved by the sine transform (Chapter 6);
4% 2. the compact nine-point stencil, solved by the same transform;
5% 3. Chebyshev collocation in both directions, solved by diagonalising the
6% one-dimensional operator.
7%
8% The nine-point stencil weights the four edge neighbours by 4 and the four
9% corner neighbours by 1,
10%
11% ( 4(u_E + u_W + u_N + u_S) + (u_NE + u_NW + u_SE + u_SW) - 20 u_P ) / (6h^2),
12%
13% which is the Mehrstellen operator. Its truncation error is
14% (h^2/12) Delta^2 u + O(h^4), and Delta^2 u = 0 for a harmonic u, so on
15% Laplace's equation the leading two error terms both vanish and the scheme is
16% sixth-order accurate. On Poisson's equation the same stencil is fourth-order
17% unless the source is also averaged, which is the two-dimensional version of
18% the trick in Chapter 3.
19%
20% It is still diagonalised by the sine basis, because it is still a sum of
21% tensor products of tridiagonal operators: with a = pi k / N and b = pi l / N,
22% the eigenvalue is 20 - 8 cos a - 8 cos b - 4 cos a cos b. So the sixth-order
23% scheme costs exactly what the second-order one costs.
25function names = method_names()
26names = {'5-point (sine transform)', '9-point compact (sine transform)', 'Chebyshev collocation'};
27end
29function [x, y, U] = solve2d(mid, N, g)
30if mid == 3
31 [x, y, U] = cheb2d(N, g);
32else
33 [x, y, U] = fd2d(mid, N, g);
34end
35end
37function [x, y, U] = fd2d(mid, N, g)
38x = (0:N)' / N;
39y = x;
40n = N - 1;
41ii = 2:N;
43G = zeros(N+1, N+1); % boundary values, zero inside
44G(1, :) = g(x', y(1));
45G(N+1, :) = g(x', y(N+1));
46G(:, 1) = g(x(1), y);
47G(:, N+1) = g(x(N+1), y);
49a = pi * (1:n)' / N;
50ca = cos(a);
52if mid == 1
53 R = G(ii-1, ii) + G(ii+1, ii) + G(ii, ii-1) + G(ii, ii+1);
54 M = 4 - 2*ca*ones(1, n) - 2*ones(n, 1)*ca';
55else
56 R = 4 * (G(ii-1, ii) + G(ii+1, ii) + G(ii, ii-1) + G(ii, ii+1)) ...
57 + (G(ii-1, ii-1) + G(ii-1, ii+1) + G(ii+1, ii-1) + G(ii+1, ii+1));
58 M = 20 - 8*ca*ones(1, n) - 8*ones(n, 1)*ca' - 4*(ca*ca');
59end
61Rh = lp_dst1(lp_dst1(R).').';
62U = G;
63U(ii, ii) = lp_dst1(lp_dst1(Rh ./ M).').';
64end
66function [x, y, U] = cheb2d(N, g)
67% Collocation on the tensor product of Chebyshev grids. The discrete problem is
68% the Sylvester equation A*U + U*A' = F, where A is the one-dimensional
69% operator on interior points. Diagonalising A once, A = V diag(lam) inv(V),
70% turns it into a division: the two-dimensional analogue of Chapter 6, with the
71% eigenvectors computed rather than known.
72[D, s] = lp_cheb(N);
73x = (1 - s) / 2;
74y = x;
75Dx = -2 * D;
76L = -(Dx * Dx);
77n = N - 1;
78ii = 2:N;
80G = zeros(N+1, N+1);
81G(1, :) = g(x', y(1));
82G(N+1, :) = g(x', y(N+1));
83G(:, 1) = g(x(1), y);
84G(:, N+1) = g(x(N+1), y);
86A = L(ii, ii);
87F = -(L(ii, 1) * G(1, ii) + L(ii, N+1) * G(N+1, ii)) ...
88 - (G(ii, 1) * L(ii, 1)' + G(ii, N+1) * L(ii, N+1)');
90[V, Lam] = eig(A);
91lam = diag(Lam);
92Vi = inv(V);
93Ft = Vi * F * Vi.';
94Ut = Ft ./ (lam * ones(1, n) + ones(n, 1) * lam.');
96U = G;
97U(ii, ii) = real(V * Ut * V.');
98end