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);
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 45 E = NaN(nm, nn);
46 T = NaN(nm, nn);
47 Q = NaN(nm, nn);
088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 48 nmax = p.nmaxs(:)';
49 timed = p.timed;
50 for k = 1:nm
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 51 [xw, yw, Uw] = solve2d(k, 8, g); % warm-up
52 for j = 1:nn
53 N = ns(j);
54 if N > nmax(k)
55 continue;
56 end
57 t0 = tic;
58 [x, y, U] = solve2d(k, N, g);
59 el = toc(t0);
60 rep = 0;
61 while timed > 0 && el < 0.005 && rep < 4
62 rep = rep + 1;
63 t0 = tic;
64 [x, y, U] = solve2d(k, N, g);
65 el = min(el, toc(t0));
66 end
67 [X, Y] = meshgrid(x(:), y(:));
68 E(k, j) = max(max(abs(U - lp_exact2d(kind, X, Y, par))));
69 T(k, j) = el * 1000;
70 end
088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 71 Q(k, :) = lp_orders(E(k, :), ns);
72 end
73 out.ns = ns;
74 out.E = lp_rows(E);
75 out.T = lp_rows(T);
76 out.orders = lp_rows(Q);
77end
79lp_write(out);