1% The five-point Laplacian on the unit square, assembled as a matrix and
2% solved by Gaussian elimination.
3%
4% On a uniform grid of N intervals in each direction, replace
5%
6% u_xx + u_yy by ( u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1) - 4 u(i,j) ) / h^2
7%
8% at the (N-1)^2 interior points. Scaling by h^2 leaves the matrix
9%
10% A = kron(I, T) + kron(T, I), T = tridiag(-1, 2, -1),
11%
12% which is what the vec of the interior array is multiplied by. Boundary
13% neighbours are known, so they move to the right-hand side.
14%
15% This is the honest, naive version: the matrix is formed in full and handed to
16% backslash. It is also the version that stops working first. The unknown count
17% is (N-1)^2, so the dense matrix has (N-1)^4 entries and the elimination costs
18% (N-1)^6 / 3 operations: N = 32 needs 7 MB and a fraction of a second, N = 64
19% needs 450 MB, and N = 128 is out of reach on any machine.
21function names = method_names()
22names = {'5-point, dense solve'};
23end
25function [x, y, U] = solve2d(mid, N, g)
26h = 1 / N;
27x = (0:N)' / N;
28y = x;
29n = N - 1;
31% boundary values, and zeros in the interior
32U = zeros(N+1, N+1);
33U(1, :) = g(x', y(1));
34U(N+1, :) = g(x', y(N+1));
35U(:, 1) = g(x(1), y);
36U(:, N+1) = g(x(N+1), y);
38% right-hand side: h^2 f (zero here, since we are solving Laplace's equation)
39% plus the known neighbours of the points next to the boundary
40B = zeros(n, n);
41B(1, :) = B(1, :) + U(1, 2:N);
42B(n, :) = B(n, :) + U(N+1, 2:N);
43B(:, 1) = B(:, 1) + U(2:N, 1);
44B(:, n) = B(:, n) + U(2:N, N+1);
46T = 2*eye(n) - diag(ones(n-1, 1), 1) - diag(ones(n-1, 1), -1);
47I = eye(n);
48A = kron(I, T) + kron(T, I);
50W = A \ B(:);
51U(2:N, 2:N) = reshape(W, n, n);
52end