% Three discretisations of the same Dirichlet problem.
%
%   1. the five-point stencil, solved by the sine transform (Chapter 6);
%   2. the compact nine-point stencil, solved by the same transform;
%   3. Chebyshev collocation in both directions, solved by diagonalising the
%      one-dimensional operator.
%
% The nine-point stencil weights the four edge neighbours by 4 and the four
% corner neighbours by 1,
%
%     ( 4(u_E + u_W + u_N + u_S) + (u_NE + u_NW + u_SE + u_SW) - 20 u_P ) / (6h^2),
%
% which is the Mehrstellen operator. Its truncation error is
% (h^2/12) Delta^2 u + O(h^4), and Delta^2 u = 0 for a harmonic u, so on
% Laplace's equation the leading two error terms both vanish and the scheme is
% sixth-order accurate. On Poisson's equation the same stencil is fourth-order
% unless the source is also averaged, which is the two-dimensional version of
% the trick in Chapter 3.
%
% It is still diagonalised by the sine basis, because it is still a sum of
% tensor products of tridiagonal operators: with a = pi k / N and b = pi l / N,
% the eigenvalue is 20 - 8 cos a - 8 cos b - 4 cos a cos b. So the sixth-order
% scheme costs exactly what the second-order one costs.

function names = method_names()
names = {'5-point (sine transform)', '9-point compact (sine transform)', 'Chebyshev collocation'};
end

function [x, y, U] = solve2d(mid, N, g)
if mid == 3
    [x, y, U] = cheb2d(N, g);
else
    [x, y, U] = fd2d(mid, N, g);
end
end

function [x, y, U] = fd2d(mid, N, g)
x = (0:N)' / N;
y = x;
n = N - 1;
ii = 2:N;

G = zeros(N+1, N+1);            % boundary values, zero inside
G(1, :) = g(x', y(1));
G(N+1, :) = g(x', y(N+1));
G(:, 1) = g(x(1), y);
G(:, N+1) = g(x(N+1), y);

a = pi * (1:n)' / N;
ca = cos(a);

if mid == 1
    R = G(ii-1, ii) + G(ii+1, ii) + G(ii, ii-1) + G(ii, ii+1);
    M = 4 - 2*ca*ones(1, n) - 2*ones(n, 1)*ca';
else
    R = 4 * (G(ii-1, ii) + G(ii+1, ii) + G(ii, ii-1) + G(ii, ii+1)) ...
        + (G(ii-1, ii-1) + G(ii-1, ii+1) + G(ii+1, ii-1) + G(ii+1, ii+1));
    M = 20 - 8*ca*ones(1, n) - 8*ones(n, 1)*ca' - 4*(ca*ca');
end

Rh = lp_dst1(lp_dst1(R).').';
U = G;
U(ii, ii) = lp_dst1(lp_dst1(Rh ./ M).').';
end

function [x, y, U] = cheb2d(N, g)
% Collocation on the tensor product of Chebyshev grids. The discrete problem is
% the Sylvester equation A*U + U*A' = F, where A is the one-dimensional
% operator on interior points. Diagonalising A once, A = V diag(lam) inv(V),
% turns it into a division: the two-dimensional analogue of Chapter 6, with the
% eigenvectors computed rather than known.
[D, s] = lp_cheb(N);
x = (1 - s) / 2;
y = x;
Dx = -2 * D;
L = -(Dx * Dx);
n = N - 1;
ii = 2:N;

G = zeros(N+1, N+1);
G(1, :) = g(x', y(1));
G(N+1, :) = g(x', y(N+1));
G(:, 1) = g(x(1), y);
G(:, N+1) = g(x(N+1), y);

A = L(ii, ii);
F = -(L(ii, 1) * G(1, ii) + L(ii, N+1) * G(N+1, ii)) ...
    - (G(ii, 1) * L(ii, 1)' + G(ii, N+1) * L(ii, N+1)');

[V, Lam] = eig(A);
lam = diag(Lam);
Vi = inv(V);
Ft = Vi * F * Vi.';
Ut = Ft ./ (lam * ones(1, n) + ones(n, 1) * lam.');

U = G;
U(ii, ii) = real(V * Ut * V.');
end