088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 1% Driver for the iterative-solver chapter.
2%
3% The editable script supplies
4%
5% names = method_names()
6% H = run_iter(mid, N, U0, UE, nit, tol)
7%
8% where H has, in its three rows and one column per iteration, the squared
9% relative residual, the squared maximum error, and the elapsed milliseconds.
10% The squares are not an oversight: taking the square roots in the MATLAB layer
11% means either a sqrt inside the iteration loop or a call after it, and either
12% one makes numbl compile that loop through its slow path. The page roots them.
13%
14% and the driver supplies the starting array U0, which holds the exact
15% boundary data on the boundary and zeros inside, and the exact solution UE for
16% measurement only. For reference it also solves the same system exactly with
17% the sine transform of Chapter 6, which gives both the discretisation error,
18% the accuracy no iterative method can beat, and the time a direct solve takes.
19p = lp_params();
20kind = p.kind;
21par = p.par;
22N = p.N;
23names = method_names();
24nm = numel(names);
26x = (0:N)' / N;
27[X, Y] = meshgrid(x, x);
28UE = lp_exact2d(kind, X, Y, par);
30U0 = zeros(N+1, N+1);
31U0(1, :) = UE(1, :);
32U0(N+1, :) = UE(N+1, :);
33U0(:, 1) = UE(:, 1);
34U0(:, N+1) = UE(:, N+1);
36% the exact solution of the discrete system, by the fast Poisson solver
37t0 = tic;
38Ud = drv_fast(U0, N);
39ms_direct = toc(t0) * 1000;
40disc_err = max(max(abs(Ud - UE)));
42nit = p.nit;
43tol = p.tol; % hoisted: a struct field read inside a loop puts numbl's
44 % compiler into its slow path, and the loop below is the timed one
45R = NaN(nm, nit+1);
46E = NaN(nm, nit+1);
47T = NaN(nm, nit+1);
48for k = 1:nm
49 % the copy goes through a function call: a slice assignment written here in
50 % the loop body would put the loop, and everything it calls, on numbl's
51 % slow path
52 H = run_iter(k, N, U0, UE, nit, tol);
53 R = drv_store(R, H, 1, k);
54 E = drv_store(E, H, 2, k);
55 T = drv_store(T, H, 3, k);
56end
58out = struct();
59out.names = names;
60out.N = N;
61out.its = 0:nit;
62out.R2 = lp_rows(R);
63out.E2 = lp_rows(E);
64out.T = lp_rows(T);
65out.discErr = disc_err;
66out.directMs = ms_direct;
67out.unknowns = (N-1)^2;
68lp_write(out);
70function A = drv_store(A, H, row, k)
71for c = 1:size(A, 2)
72 A(k, c) = H(row, c);
73end
74end
76function U = drv_fast(U, N)
77n = N - 1;
78B = zeros(n, n);
79B(1, :) = B(1, :) + U(1, 2:N);
80B(n, :) = B(n, :) + U(N+1, 2:N);
81B(:, 1) = B(:, 1) + U(2:N, 1);
82B(:, n) = B(:, n) + U(2:N, N+1);
83lam = 2 * (1 - cos(pi * (1:n)' / N));
84D = lam * ones(1, n) + ones(n, 1) * lam';
85Bh = lp_dst1(lp_dst1(B).').';
86U(2:N, 2:N) = lp_dst1(lp_dst1(Bh ./ D).').';
87end