088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 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%
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 14% Timings come from tic/toc around the solve and nothing else. Each method gets
15% a warm-up call before its first timed one, because the first call also pays
16% for compiling the method, and anything under 5 ms is repeated with the best
17% kept: a sub-millisecond solve otherwise measures the clock.
088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 18p = lp_params();
19kind = p.kind;
20names = method_names();
21nm = numel(names);
22ub = lp_exact1d(kind, [0; 1]);
23f = @(xx) drv_rhs(kind, xx);
25out = struct();
26out.names = names;
27out.kind = kind;
29if strcmp(p.mode, 'single')
30 t0 = tic;
31 [x, u] = solve1d(p.mid, p.N, f, ub(1), ub(2));
32 out.ms = toc(t0) * 1000;
33 uex = lp_exact1d(kind, x(:));
34 err = u(:) - uex(:);
35 out.x = x(:)';
36 out.u = u(:)';
37 out.uex = uex(:)';
38 out.err = err';
39 out.maxerr = max(abs(err));
40 out.n = numel(x) - 1;
41 xf = linspace(0, 1, p.ngrid)';
42 out.xf = xf';
43 out.uf = lp_exact1d(kind, xf)';
44else
45 ns = p.ns(:)';
46 nn = numel(ns);
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 47 E = NaN(nm, nn);
48 T = NaN(nm, nn);
49 Q = NaN(nm, nn);
088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 50 nmax = p.nmaxs(:)';
51 timed = p.timed;
52 for k = 1:nm
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 53 [xw, uw] = solve1d(k, 8, f, ub(1), ub(2)); % warm-up
54 for j = 1:nn
55 N = ns(j);
56 if N > nmax(k)
57 continue;
58 end
59 t0 = tic;
60 [x, u] = solve1d(k, N, f, ub(1), ub(2));
61 el = toc(t0);
62 rep = 0;
63 while timed > 0 && el < 0.005 && rep < 4
64 rep = rep + 1;
65 t0 = tic;
66 [x, u] = solve1d(k, N, f, ub(1), ub(2));
67 el = min(el, toc(t0));
68 end
69 uex = lp_exact1d(kind, x(:));
70 E(k, j) = max(abs(u(:) - uex(:)));
71 T(k, j) = el * 1000;
72 end
088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 73 Q(k, :) = lp_orders(E(k, :), ns);
74 end
75 out.ns = ns;
76 out.E = lp_rows(E);
77 out.T = lp_rows(T);
78 out.orders = lp_rows(Q);
79end
81lp_write(out);
83function y = drv_rhs(kind, x)
84[~, y] = lp_exact1d(kind, x);
85end