% Three methods on the same problem: the two finite difference schemes of the
% previous chapters, and Chebyshev collocation.
%
% Collocation drops the idea of a local stencil. Put the N+1 Chebyshev points
%
%     x_j = ( 1 - cos(pi j / N) ) / 2,     j = 0, ..., N,
%
% on [0,1], let p be the polynomial of degree N through the unknown values, and
% require -p'' = f at the interior points. Differentiation of a polynomial
% given by its values at fixed nodes is a matrix, so the discrete operator is
% the dense matrix D^2 rather than a three-term recurrence. Every unknown is
% coupled to every other, the solve costs O(N^3) instead of O(N), and in
% exchange the error decays faster than any power of N for analytic solutions.
%
% The clustering of the nodes near the ends is not decoration. It is what makes
% polynomial interpolation of high degree stable, and equally spaced points of
% the same degree would be useless.

function names = method_names()
names = {'3-point differences', 'compact 4th order', 'Chebyshev collocation'};
end

function [x, u] = solve1d(mid, N, f, ua, ub)
if mid == 3
    [x, u] = cheb_solve(N, f, ua, ub);
else
    [x, u] = fd_solve(mid, N, f, ua, ub);
end
end

function [x, u] = cheb_solve(N, f, ua, ub)
[D, s] = lp_cheb(N);
x = (1 - s) / 2;          % Chebyshev points, ascending, on [0, 1]
Dx = -2 * D;              % d/dx after the affine map
L = -(Dx * Dx);           % the operator -d^2/dx^2, as a dense matrix

ii = 2:N;                 % interior collocation points
A = L(ii, ii);
b = f(x(ii)) - L(ii, 1) * ua - L(ii, N+1) * ub;
u = [ua; A \ b; ub];
end

function [x, u] = fd_solve(mid, N, f, ua, ub)
h = 1 / N;
x = (0:N)' / N;
fx = f(x);
if mid == 1
    r = h^2 * fx(2:N);
else
    r = h^2 * (fx(1:N-1) + 10*fx(2:N) + fx(3:N+1)) / 12;
end
r(1) = r(1) + ua;
r(end) = r(end) + ub;
u = [ua; thomas(N-1, r); ub];
end

function w = thomas(n, d)
c = zeros(n, 1);
e = zeros(n, 1);
b = 2;
c(1) = -1 / b;
e(1) = d(1) / b;
for i = 2:n
    b = 2 + c(i-1);
    c(i) = -1 / b;
    e(i) = (d(i) + e(i-1)) / b;
end
w = zeros(n, 1);
w(n) = e(n);
for i = n-1:-1:1
    w(i) = e(i) - c(i) * w(i+1);
end
end