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. numbl's
60 % jsondecode turns a JSON array of arrays into a cell of column vectors
61 % where MATLAB would give a matrix, so 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 = NaN(nm, nn);
70 T = NaN(nm, nn);
71 Q = NaN(nm, nn);
72 for k = 1:nm
73 vw = solve_l(k, NS(k, 1), g, px(1), py(1)); % warm-up
74 for j = 1:nn
75 n = NS(k, j);
76 if ~isfinite(n) || n <= 0
77 continue;
78 end
79 t0 = tic;
80 v = solve_l(k, n, g, px, py);
81 T(k, j) = toc(t0) * 1000;
82 E(k, j) = max(abs(v - ve));
83 end
84 Q(k, :) = lp_orders(E(k, :), NS(k, :));
85 end
86 out.ns = NS(1, :);
87 out.NS = lp_rows(NS);
88 out.E = lp_rows(E);
89 out.T = lp_rows(T);
90 out.orders = lp_rows(Q);
91end
93lp_write(out);