088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 1% Driver for the two-dimensional labs on the unit square.
2%
3% The editable script supplies
4%
5% names = method_names()
6% [x, y, U] = solve2d(mid, N, g)
7%
8% where g(x, y) is the Dirichlet data. Because every exact solution here is
9% harmonic on the whole plane, g happens to be defined everywhere, but a method
10% is only entitled to evaluate it on the boundary of the square; anything else
11% is looking up the answer.
12%
13% U is the solution on the tensor grid x by y, with U(i, j) at (x(j), y(i)), and
14% the driver measures the error against the exact solution on exactly that grid.
15p = lp_params();
16kind = p.kind;
17par = p.par;
18names = method_names();
19nm = numel(names);
20g = @(xx, yy) lp_exact2d(kind, xx, yy, par);
22out = struct();
23out.names = names;
25if strcmp(p.mode, 'single')
26 t0 = tic;
27 [x, y, U] = solve2d(p.mid, p.N, g);
28 out.ms = toc(t0) * 1000;
29 x = x(:);
30 y = y(:);
31 [X, Y] = meshgrid(x, y);
32 UE = lp_exact2d(kind, X, Y, par);
33 Err = U - UE;
34 out.maxerr = max(max(abs(Err)));
35 out.unknowns = (numel(x) - 2) * (numel(y) - 2);
36 [Us, xs, ys] = lp_thin(U, x, y, p.maxpix);
37 Es = lp_thin(abs(Err), x, y, p.maxpix);
38 out.x = xs(:)';
39 out.y = ys(:)';
40 out.U = lp_rows(Us);
41 out.Err = lp_rows(Es);
42else
43 ns = p.ns(:)';
44 nn = numel(ns);
45 E = zeros(nm, nn);
46 T = zeros(nm, nn);
47 Q = zeros(nm, nn);
48 nmax = p.nmaxs(:)';
49 timed = p.timed;
50 ET = drv_methods2(1, NaN(2*nm, nn), ns, g, kind, par, nmax, timed, nm);
51 for k = 1:nm
52 E(k, :) = ET(2*k-1, :);
53 T(k, :) = ET(2*k, :);
54 Q(k, :) = lp_orders(E(k, :), ns);
55 end
56 out.ns = ns;
57 out.E = lp_rows(E);
58 out.T = lp_rows(T);
59 out.orders = lp_rows(Q);
60end
62lp_write(out);
64function ET = drv_methods2(k, ET, ns, g, kind, par, nmax, timed, nm)
65% One method per level of the recursion. See the note in the one-dimensional
66% driver: a for loop here would put every timed solve on numbl's slow path.
67if k > nm
68 return;
69end
70[xw, yw, Uw] = solve2d(k, 8, g); % warm-up, so the first timing is honest
71R = drv_res2(1, NaN(2, numel(ns)), k, ns, g, kind, par, nmax(k), timed);
72ET(2*k-1, :) = R(1, :);
73ET(2*k, :) = R(2, :);
74ET = drv_methods2(k+1, ET, ns, g, kind, par, nmax, timed, nm);
75end
77function R = drv_res2(j, R, mid, ns, g, kind, par, nmax, timed)
78if j > numel(ns)
79 return;
80end
81N = ns(j);
82if N <= nmax
83 t0 = tic;
84 [x, y, U] = solve2d(mid, N, g);
85 el = toc(t0);
86 if timed > 0 && el < 0.005
87 el = drv_best2(4, el, mid, N, g);
88 end
89 [X, Y] = meshgrid(x(:), y(:));
90 R(1, j) = max(max(abs(U - lp_exact2d(kind, X, Y, par))));
91 R(2, j) = el * 1000;
92end
93R = drv_res2(j+1, R, mid, ns, g, kind, par, nmax, timed);
94end
96function el = drv_best2(k, el, mid, N, g)
97if k == 0
98 return;
99end
100t0 = tic;
101[x, y, U] = solve2d(mid, N, g);
102e2 = toc(t0);
103if e2 < el
104 el = e2;
105end
106el = drv_best2(k-1, el, mid, N, g);
107end