1% Driver for the one-dimensional labs.
2%
3% The editable script supplies two functions:
4%
5% names = method_names() a cell array of method labels
6% [x, u] = solve1d(mid, N, f, ua, ub) method mid, resolution N
7%
8% solve1d picks its own grid x (uniform, graded, Chebyshev, whatever the method
9% wants), returns the solution values on that grid, and the driver measures the
10% error against the manufactured solution at exactly those points. f is a
11% function handle giving the source term, so a method is free to evaluate it
12% anywhere.
13%
14% The refinement sweep is written as a recursion rather than a for loop. That
15% is not a stylistic choice: numbl compiles a call made from inside a loop
16% through a generic path whenever the callee itself contains loops that call
17% further functions, and the solves being timed here are exactly of that shape.
18% The penalty is a factor of five to twenty, it falls unevenly on different
19% methods, and it would therefore corrupt every timing in this tutorial.
20% Recursion is compiled normally. See CLAUDE.md.
21p = lp_params();
22kind = p.kind;
23names = method_names();
24nm = numel(names);
25ub = lp_exact1d(kind, [0; 1]);
26f = @(xx) drv_rhs(kind, xx);
28out = struct();
29out.names = names;
30out.kind = kind;
32if strcmp(p.mode, 'single')
33 t0 = tic;
34 [x, u] = solve1d(p.mid, p.N, f, ub(1), ub(2));
35 out.ms = toc(t0) * 1000;
36 uex = lp_exact1d(kind, x(:));
37 err = u(:) - uex(:);
38 out.x = x(:)';
39 out.u = u(:)';
40 out.uex = uex(:)';
41 out.err = err';
42 out.maxerr = max(abs(err));
43 out.n = numel(x) - 1;
44 xf = linspace(0, 1, p.ngrid)';
45 out.xf = xf';
46 out.uf = lp_exact1d(kind, xf)';
47else
48 ns = p.ns(:)';
49 nn = numel(ns);
50 E = zeros(nm, nn);
51 T = zeros(nm, nn);
52 Q = zeros(nm, nn);
53 nmax = p.nmaxs(:)';
54 timed = p.timed;
55 ET = drv_methods(1, NaN(2*nm, nn), ns, f, ub(1), ub(2), kind, timed, nmax, nm);
56 for k = 1:nm
57 E(k, :) = ET(2*k-1, :);
58 T(k, :) = ET(2*k, :);
59 Q(k, :) = lp_orders(E(k, :), ns);
60 end
61 out.ns = ns;
62 out.E = lp_rows(E);
63 out.T = lp_rows(T);
64 out.orders = lp_rows(Q);
65end
67lp_write(out);
69function ET = drv_methods(k, ET, ns, f, ua, ub, kind, timed, nmax, nm)
70% One method per level of the recursion; row 2k-1 of ET collects its errors and
71% row 2k its milliseconds.
72if k > nm
73 return;
74end
75% a warm-up call, so that the first timed solve is not also paying for the
76% compilation of the method
77[xw, uw] = solve1d(k, 8, f, ua, ub);
78R = drv_res(1, NaN(2, numel(ns)), k, ns, f, ua, ub, kind, timed, nmax(k));
79ET(2*k-1, :) = R(1, :);
80ET(2*k, :) = R(2, :);
81ET = drv_methods(k+1, ET, ns, f, ua, ub, kind, timed, nmax, nm);
82end
84function R = drv_res(j, R, mid, ns, f, ua, ub, kind, timed, nmax)
85% One resolution per level of the recursion.
86if j > numel(ns)
87 return;
88end
89N = ns(j);
90if N <= nmax
91 t0 = tic;
92 [x, u] = solve1d(mid, N, f, ua, ub);
93 el = toc(t0);
94 if timed > 0 && el < 0.005
95 % a sub-millisecond solve says more about the clock than about the
96 % method, so take the best of a few repeats
97 el = drv_best(4, el, mid, N, f, ua, ub);
98 end
99 uex = lp_exact1d(kind, x(:));
100 R(1, j) = max(abs(u(:) - uex(:)));
101 R(2, j) = el * 1000;
102end
103R = drv_res(j+1, R, mid, ns, f, ua, ub, kind, timed, nmax);
104end
106function el = drv_best(k, el, mid, N, f, ua, ub)
107if k == 0
108 return;
109end
110t0 = tic;
111[x, u] = solve1d(mid, N, f, ua, ub);
112e2 = toc(t0);
113if e2 < el
114 el = e2;
115end
116el = drv_best(k-1, el, mid, N, f, ua, ub);
117end
119function y = drv_rhs(kind, x)
120[~, y] = lp_exact1d(kind, x);
121end