088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 1% The same five-point system, solved two ways.
2%
3% The matrix is A = kron(I,T) + kron(T,I) with T = tridiag(-1,2,-1), and the
4% eigenvectors of T are known exactly: the sine vectors
5%
6% v_k(j) = sin(pi j k / N), T v_k = lambda_k v_k,
7% lambda_k = 2 (1 - cos(pi k / N)), k = 1, ..., N-1.
8%
9% So the eigenvectors of A are the products v_k(i) v_l(j) and its eigenvalues
10% are the sums lambda_k + lambda_l. Expanding the right-hand side in that basis,
11% dividing by the eigenvalues and transforming back solves the system exactly,
12% and the expansion is a discrete sine transform, which one FFT computes in
13% O(N log N) per line. Total: O(N^2 log N) for N^2 unknowns, no matrix formed
14% and nothing stored beyond the grid itself.
15%
16% This is not an iterative method and not an approximation. It returns the same
17% vector the dense solve returns, to within roundoff, and it does so for grids
18% the dense solve cannot represent.
20function names = method_names()
21names = {'5-point, dense solve', '5-point, sine transform'};
22end
24function [x, y, U] = solve2d(mid, N, g)
25h = 1 / N;
26x = (0:N)' / N;
27y = x;
28n = N - 1;
30U = zeros(N+1, N+1);
31U(1, :) = g(x', y(1));
32U(N+1, :) = g(x', y(N+1));
33U(:, 1) = g(x(1), y);
34U(:, N+1) = g(x(N+1), y);
36B = zeros(n, n);
37B(1, :) = B(1, :) + U(1, 2:N);
38B(n, :) = B(n, :) + U(N+1, 2:N);
39B(:, 1) = B(:, 1) + U(2:N, 1);
40B(:, n) = B(:, n) + U(2:N, N+1);
42if mid == 1
43 T = 2*eye(n) - diag(ones(n-1, 1), 1) - diag(ones(n-1, 1), -1);
44 A = kron(eye(n), T) + kron(T, eye(n));
45 U(2:N, 2:N) = reshape(A \ B(:), n, n);
46else
47 lam = 2 * (1 - cos(pi * (1:n)' / N)); % eigenvalues of T
48 D = lam * ones(1, n) + ones(n, 1) * lam'; % eigenvalues of A
49 Bh = lp_dst1(lp_dst1(B).').'; % to the sine basis, both ways
50 U(2:N, 2:N) = lp_dst1(lp_dst1(Bh ./ D).').'; % divide, and back again
51end
52end