088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 1% Driver for the boundary integral chapter.
2%
3% The editable script supplies
4%
5% names = method_names()
6% v = solve_bie(mid, n, a, nlobe, par, px, py)
7%
8% returning the computed solution at the points (px, py) inside the curve
9% r(t) = 1 + a cos(nlobe t). n is the method's resolution: boundary nodes for
10% the Nystrom method, source points for the method of fundamental solutions.
11% par is the auxiliary parameter each method needs (the radius of the source
12% circle, for the second one).
13%
14% The Dirichlet data is the trace of lp_bexact, so the exact solution is known
15% everywhere inside and the error needs no reference computation.
16p = lp_params();
17a = p.a;
18nlobe = p.nlobe;
19names = method_names();
20nm = numel(names);
22out = struct();
23out.names = names;
25% the curve itself, for the pictures
26tt = linspace(0, 2*pi, 401)';
27C = lp_curve(a, nlobe, tt);
28out.cx = C(:,1)';
29out.cy = C(:,2)';
31if strcmp(p.mode, 'field')
32 % a Cartesian grid over the bounding box, masked to the interior
33 m = 96;
34 rmax = 1 + abs(a) + 0.05;
35 xs = linspace(-rmax, rmax, m)';
36 [X, Y] = meshgrid(xs, xs);
37 th = atan2(Y, X);
38 rr = sqrt(X.^2 + Y.^2);
39 rb = 1 + a * cos(nlobe * th);
40 inside = rr < rb - 1e-9;
41 px = X(inside);
42 py = Y(inside);
44 t0 = tic;
45 v = solve_bie(p.mid, p.n, a, nlobe, p.par, px, py);
46 out.ms = toc(t0) * 1000;
47 ve = lp_bexact(px, py);
48 out.maxerr = max(abs(v - ve));
50 U = NaN(m, m);
51 E = NaN(m, m);
52 U(inside) = v;
53 E(inside) = abs(v - ve);
54 out.x = xs';
55 out.y = xs';
56 out.U = lp_rows(U);
57 out.Err = lp_rows(E);
58 out.unknowns = p.n;
59else
60 % errors on a polar lattice: four rings, the last of them close enough to
61 % the boundary to show what close evaluation costs
62 nth = 60;
63 th = 2*pi*(0:nth-1)'/nth;
64 rb = 1 + a * cos(nlobe * th);
65 rho = [0.25; 0.5; 0.75];
66 px = [rho(1)*rb .* cos(th); rho(2)*rb .* cos(th); rho(3)*rb .* cos(th)];
67 py = [rho(1)*rb .* sin(th); rho(2)*rb .* sin(th); rho(3)*rb .* sin(th)];
68 ve = lp_bexact(px, py);
69 % a ring at 98 per cent of the radius, where a layer potential is hardest
70 % to evaluate accurately
71 qx = 0.98 * rb .* cos(th);
72 qy = 0.98 * rb .* sin(th);
73 qe = lp_bexact(qx, qy);
75 ns = p.ns(:)';
76 nn = numel(ns);
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 77 E = NaN(nm, nn);
78 T = NaN(nm, nn);
79 Q = NaN(nm, nn);
80 EN = NaN(nm, nn);
088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 81 for k = 1:nm
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 82 vw = solve_bie(k, 16, a, nlobe, p.par, px(1), py(1)); % warm-up
83 for j = 1:nn
84 n = ns(j);
1a8bb6aMake the cross-language benchmark measure the languagesJeremy Magland 85 % Keep the fastest of several runs: three at least whenever a
86 % solve is cheap enough to allow it, many more when it is very
87 % cheap, one when it is expensive. A single measurement of a
88 % twenty-millisecond solve varies by a factor of two with whatever
89 % else the machine is doing, and the numbers here are quoted in the
90 % prose. Chapter 11 times itself the same way.
91 el = Inf;
92 tot = 0;
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 93 rep = 0;
1a8bb6aMake the cross-language benchmark measure the languagesJeremy Magland 94 while rep == 0 || (rep < 50 && tot < 0.6 && (rep < 3 || tot < 0.1))
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 95 t0 = tic;
96 v = solve_bie(k, n, a, nlobe, p.par, px, py);
98 el = min(el, dt);
99 tot = tot + dt;
100 rep = rep + 1;
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 101 end
102 q = solve_bie(k, n, a, nlobe, p.par, qx, qy);
103 E(k, j) = max(abs(v - ve));
104 EN(k, j) = max(abs(q - qe));
105 T(k, j) = el * 1000;
106 end
088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 107 Q(k, :) = lp_orders(E(k, :), ns);
108 end
109 out.Enear = lp_rows(EN);
110 out.ns = ns;
111 out.E = lp_rows(E);
112 out.T = lp_rows(T);
113 out.orders = lp_rows(Q);
114end
116lp_write(out);