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%
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 8% where H has, in its three rows and one column per iteration, the relative
9% residual, the maximum error, and the elapsed milliseconds. Columns past the
10% iteration the method stopped at are NaN.
12% and the driver supplies the starting array U0, which holds the exact
13% boundary data on the boundary and zeros inside, and the exact solution UE for
14% measurement only. For reference it also solves the same system exactly with
15% the sine transform of Chapter 6, which gives both the discretisation error,
16% the accuracy no iterative method can beat, and the time a direct solve takes.
17p = lp_params();
18kind = p.kind;
19par = p.par;
20N = p.N;
21names = method_names();
22nm = numel(names);
24x = (0:N)' / N;
25[X, Y] = meshgrid(x, x);
26UE = lp_exact2d(kind, X, Y, par);
28U0 = zeros(N+1, N+1);
29U0(1, :) = UE(1, :);
30U0(N+1, :) = UE(N+1, :);
31U0(:, 1) = UE(:, 1);
32U0(:, N+1) = UE(:, N+1);
34% the exact solution of the discrete system, by the fast Poisson solver
35t0 = tic;
36Ud = drv_fast(U0, N);
37ms_direct = toc(t0) * 1000;
38disc_err = max(max(abs(Ud - UE)));
40nit = p.nit;
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 41tol = p.tol;
088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 42R = NaN(nm, nit+1);
43E = NaN(nm, nit+1);
44T = NaN(nm, nit+1);
45for k = 1:nm
46 H = run_iter(k, N, U0, UE, nit, tol);
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 47 R(k, :) = H(1, :);
48 E(k, :) = H(2, :);
49 T(k, :) = H(3, :);
52out = struct();
53out.names = names;
54out.N = N;
55out.its = 0:nit;
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 56out.R = lp_rows(R);
57out.E = lp_rows(E);
088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 58out.T = lp_rows(T);
59out.discErr = disc_err;
60out.directMs = ms_direct;
61out.unknowns = (N-1)^2;
62lp_write(out);
64function U = drv_fast(U, N)
65n = N - 1;
66B = zeros(n, n);
67B(1, :) = B(1, :) + U(1, 2:N);
68B(n, :) = B(n, :) + U(N+1, 2:N);
69B(:, 1) = B(:, 1) + U(2:N, 1);
70B(:, n) = B(:, n) + U(2:N, N+1);
71lam = 2 * (1 - cos(pi * (1:n)' / N));
72D = lam * ones(1, n) + ones(n, 1) * lam';
73Bh = lp_dst1(lp_dst1(B).').';
74U(2:N, 2:N) = lp_dst1(lp_dst1(Bh ./ D).').';
75end