1% Driver for the L-shaped domain.
2%
3% The editable script supplies
4%
5% names = method_names()
6% v = solve_l(mid, n, g, px, py)
7%
8% and must return the computed solution at the points (px, py), whatever grid
9% or basis it used internally. Comparing two methods at a common set of points,
10% rather than each on its own grid, is what makes a finite difference scheme
11% and a basis expansion comparable at all.
12%
13% g(x, y) is the Dirichlet data, which here is the exact solution restricted to
14% the boundary. n is the method's own resolution parameter: intervals per unit
15% length for the difference scheme, number of basis terms for the expansion.
16p = lp_params();
17kind = p.kind;
18names = method_names();
19nm = numel(names);
20g = @(xx, yy) lp_lexact(kind, xx, yy);
22out = struct();
23out.names = names;
25if strcmp(p.mode, 'field')
26 % a picture on a fixed grid of spacing 1/32, masked outside the domain
27 m = 32;
28 xs = (-m:m)' / m;
29 [X, Y] = meshgrid(xs, xs);
30 inside = ~(X > 1e-12 & Y < -1e-12);
31 px = X(inside);
32 py = Y(inside);
33 t0 = tic;
34 v = solve_l(p.mid, p.n, g, px, py);
35 out.ms = toc(t0) * 1000;
36 ve = lp_lexact(kind, px, py);
37 out.maxerr = max(abs(v - ve));
39 U = NaN(2*m+1, 2*m+1);
40 E = NaN(2*m+1, 2*m+1);
41 U(inside) = v;
42 E(inside) = abs(v - ve);
43 out.x = xs';
44 out.y = xs';
45 out.U = lp_rows(U);
46 out.Err = lp_rows(E);
47else
48 % errors at the nodes of a fixed grid of spacing 1/8, well away from the
49 % plotting grid and coarse enough that every difference grid contains it
50 m = 8;
51 xs = (-m:m)' / m;
52 [X, Y] = meshgrid(xs, xs);
53 inside = ~(X > 1e-12 & Y < -1e-12) & abs(X) < 1 - 1e-12 & abs(Y) < 1 - 1e-12;
54 px = X(inside);
55 py = Y(inside);
56 ve = lp_lexact(kind, px, py);
58 % each method gets its own list of resolutions: intervals per unit length
59 % for the difference scheme, basis terms for the expansion
60 % numbl's jsondecode turns a JSON array of arrays into a cell of column
61 % vectors, where MATLAB would give a matrix; unpack it by hand
62 NSc = p.nsm;
63 nn = numel(NSc{1});
64 NS = zeros(nm, nn);
65 for k = 1:nm
66 row = NSc{k};
67 NS(k, :) = row(:)';
68 end
69 E = zeros(nm, nn);
70 T = zeros(nm, nn);
71 Q = zeros(nm, nn);
72 ETall = drv_methodsl(1, NaN(2*nm, nn), NS, g, px, py, ve, nm);
73 for k = 1:nm
74 E(k, :) = ETall(2*k-1, :);
75 T(k, :) = ETall(2*k, :);
76 end
77 for k = 1:nm
78 Q(k, :) = lp_orders(E(k, :), NS(k, :));
79 end
80 out.ns = NS(1, :);
81 out.NS = lp_rows(NS);
82 out.E = lp_rows(E);
83 out.T = lp_rows(T);
84 out.orders = lp_rows(Q);
85end
87lp_write(out);
89function ET = drv_methodsl(k, ET, NS, g, px, py, ve, nm)
90% Recursion rather than a loop: see the note in the one-dimensional driver.
91if k > nm
92 return;
93end
94vw = solve_l(k, NS(k, 1), g, px(1), py(1)); % warm-up
95R = drv_resl(1, NaN(2, size(NS, 2)), k, NS(k, :), g, px, py, ve);
96ET(2*k-1, :) = R(1, :);
97ET(2*k, :) = R(2, :);
98ET = drv_methodsl(k+1, ET, NS, g, px, py, ve, nm);
99end
101function R = drv_resl(j, R, mid, ns, g, px, py, ve)
102if j > numel(ns)
103 return;
104end
105n = ns(j);
106if isfinite(n) && n > 0
107 t0 = tic;
108 v = solve_l(mid, n, g, px, py);
109 R(2, j) = toc(t0) * 1000;
110 R(1, j) = max(abs(v - ve));
111end
112R = drv_resl(j+1, R, mid, ns, g, px, py, ve);
113end
115function A = drv_store(A, ET, row, k)
116for c = 1:size(A, 2)
117 A(k, c) = ET(row, c);
118end
119end