6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 1% Four methods for the Dirichlet problem on the Pac-Man domain: the unit disk
2% with the quadrant x > 0, y < 0 removed. The boundary is a curved arc and two
3% straight edges; the edges meet the arc in two right-angle corners, and each
4% other in a re-entrant corner of interior angle 3*pi/2 at the origin. Every
5% difficulty this tutorial has met is present at once: a boundary that cuts
6% through grid cells (Chapter 10's curve), and a corner singularity r^(2/3)
7% (Chapter 9's corner).
8%
9% 1. SHORTLEY-WELLER + SOR. The five-point scheme, with the legs of the
10% stencil shortened wherever a neighbour falls outside the arc, so the
11% boundary condition is imposed on the arc itself and not on a staircase
12% approximation of it. The shortened stencil is no longer symmetric, so
13% conjugate gradients is not available; the solver is SOR with the
14% overrelaxation parameter near its optimum.
15%
16% 2. THE HARMONIC BASIS. The method of particular solutions from Chapter 9,
17% unchanged in spirit: harmonic polynomials Re z^m, Im z^m for the smooth
18% part, plus the corner's own functions r^(2k/3) sin(2k theta/3), fitted
19% to the boundary data by least squares.
20%
21% 3. DOUBLE LAYER, GRADED PANELS. The Nystrom method of Chapter 10 cannot
22% use the trapezoid rule here, because the boundary is not smooth: the
23% density has its own singularity at each corner. The quadrature is
24% composite Gauss-Legendre on panels refined dyadically towards the three
25% corners, which resolves the density algebraically panel by panel and
26% geometrically in the refinement depth.
27%
28% 4. CHEBYSHEV IN POLAR COORDINATES. In polar coordinates the Pac-Man is a
29% rectangle, (r, theta) in (0,1) x (0, 3*pi/2), and Laplace's equation
30% becomes r^2 u_rr + r u_r + u_tt = 0, which is separable: the matrix
31% decomposition of Chapter 7 applies verbatim. The whole side r = 0 of
32% the rectangle maps to the corner, where u is known, so the coordinate
33% degeneracy costs nothing to impose; what it costs is accuracy, because
34% the r^(2/3) singularity now sits on the edge of the expansion domain.
36function names = method_names()
37names = {'Shortley-Weller + SOR', 'harmonic basis, least squares', ...
38 'double layer, graded panels', 'Chebyshev, polar rectangle'};
39end
41function v = solve_pac(mid, n, g, px, py)
42if mid == 1
43 v = fd_pac(n, g, px, py);
44elseif mid == 2
45 v = mps_pac(n, g, px, py);
46elseif mid == 3
47 v = bie_pac(n, g, px, py);
48else
49 v = spec_pac(n, g, px, py);
50end
51end
53function d = dof_pac(mid, n)
54if mid == 1
55 d = fd_count(n);
56elseif mid == 2
57 d = 3*n + 1;
58elseif mid == 3
59 [qx, qy, qnx, qny, qw, qkap] = pac_panels(n);
60 d = numel(qx);
61else
62 d = (n - 1)^2;
63end
64end
66% --------------------------------------------------- finite differences -----
68function v = fd_pac(N, g, px, py)
69% N intervals per unit length on (-1:1/N:1)^2. Nodes on the two straight
70% edges are grid points and carry the data directly; nodes adjacent to the
71% arc get Shortley-Weller legs, shortened to where the arc actually is.
72h = 1 / N;
73m = 2*N + 1;
74xs = (-N:N)' * h;
75[X, Y] = meshgrid(xs, xs);
76R2 = X.^2 + Y.^2;
78% start from the data everywhere: boundary nodes keep it, and for the
79% unknowns it is a continuous extension, which is a decent initial guess
80U = g(X, Y);
81inquad = X > 1e-12 & Y < -1e-12;
82U(inquad) = 0;
83onedge = (abs(Y) < 1e-12 & X > -1e-12) | (abs(X) < 1e-12 & Y < 1e-12);
84unknown = R2 < 1 - 1e-12 & ~inquad & ~onedge;
85mask = double(unknown);
87% Shortley-Weller coefficients. Every leg is h except where the neighbour
88% lies beyond the arc; there the leg is shortened to the crossing and the
89% boundary value at the crossing moves to the right-hand side.
90CE = zeros(m, m);
91CW = zeros(m, m);
92CN = zeros(m, m);
93CS = zeros(m, m);
94CP = ones(m, m);
95B = zeros(m, m);
96for j = 2:m-1
97 for i = 2:m-1
98 if mask(i, j) == 0
99 continue;
100 end
101 x = xs(j);
102 y = xs(i);
103 hE = h; hW = h; hN = h; hS = h;
104 gE = 0; gW = 0; gN = 0; gS = 0;
105 if (x + h)^2 + y^2 >= 1 - 1e-12
106 hE = max(sqrt(1 - y^2) - x, 1e-8 * h);
107 gE = 1;
108 end
109 if (x - h)^2 + y^2 >= 1 - 1e-12
110 hW = max(x + sqrt(1 - y^2), 1e-8 * h);
111 gW = 1;
112 end
113 if x^2 + (y + h)^2 >= 1 - 1e-12
114 hN = max(sqrt(1 - x^2) - y, 1e-8 * h);
115 gN = 1;
116 end
117 if x^2 + (y - h)^2 >= 1 - 1e-12
118 hS = max(y + sqrt(1 - x^2), 1e-8 * h);
119 gS = 1;
120 end
121 cE = 2 / (hE * (hE + hW));
122 cW = 2 / (hW * (hE + hW));
123 cN = 2 / (hN * (hN + hS));
124 cS = 2 / (hS * (hN + hS));
125 CP(i, j) = cE + cW + cN + cS;
126 % a shortened leg reaches a point of the arc, not a grid node: its
127 % (known) value multiplies the coefficient and moves to the rhs
128 if gE == 1
129 B(i, j) = B(i, j) + cE * g(x + hE, y);
130 else
131 CE(i, j) = cE;
132 end
133 if gW == 1
134 B(i, j) = B(i, j) + cW * g(x - hW, y);
135 else
136 CW(i, j) = cW;
137 end
138 if gN == 1
139 B(i, j) = B(i, j) + cN * g(x, y + hN);
140 else
141 CN(i, j) = cN;
142 end
143 if gS == 1
144 B(i, j) = B(i, j) + cS * g(x, y - hS);
145 else
146 CS(i, j) = cS;
147 end
148 end
149end
151% SOR, swept until the largest Gauss-Seidel correction is far below the
152% discretisation error. The shortened stencils make the matrix
153% nonsymmetric, which is what rules out the CG solver of Chapter 9.
154om = 2 / (1 + sin(pi * h / 2));
155maxit = 60 * N;
156for it = 1:maxit
157 du = 0;
158 for j = 2:m-1
159 for i = 2:m-1
160 if mask(i, j) > 0
161 unew = (CE(i,j) * U(i,j+1) + CW(i,j) * U(i,j-1) ...
162 + CN(i,j) * U(i+1,j) + CS(i,j) * U(i-1,j) + B(i,j)) / CP(i,j);
163 d = unew - U(i, j);
164 U(i, j) = U(i, j) + om * d;
165 if abs(d) > du
166 du = abs(d);
167 end
168 end
169 end
170 end
171 if du < 1e-9
172 break;
173 end
174end
176% the requested points are grid points, so sampling is a lookup
177v = zeros(numel(px), 1);
178for k = 1:numel(px)
179 i = round(py(k) / h) + N + 1;
180 j = round(px(k) / h) + N + 1;
181 v(k) = U(i, j);
182end
183end
185function d = fd_count(N)
186h = 1 / N;
187xs = (-N:N)' * h;
188[X, Y] = meshgrid(xs, xs);
189R2 = X.^2 + Y.^2;
190inquad = X > 1e-12 & Y < -1e-12;
191onedge = (abs(Y) < 1e-12 & X > -1e-12) | (abs(X) < 1e-12 & Y < 1e-12);
192d = sum(sum(double(R2 < 1 - 1e-12 & ~inquad & ~onedge)));
193end
195% ------------------------------------------- the method of particular sols --
197function v = mps_pac(n, g, px, py)
198% Harmonic polynomials up to degree n plus n corner functions, fitted by
199% least squares at about three boundary points per unknown.
200nb = 3*n + 1;
201[bx, by] = pac_bpoints(3 * nb);
202A = pac_basis(bx, by, n, n);
203c = A \ g(bx, by);
204v = pac_basis(px, py, n, n) * c;
205end
207function A = pac_basis(x, y, M, K)
208% Re z^m, Im z^m for m = 0..M, then r^(2k/3) sin(2k theta/3) for K values of
209% k that are not multiples of 3 (those already are harmonic polynomials, and
210% listing them twice would make the matrix exactly rank deficient).
211np = numel(x);
212A = zeros(np, 2*M + 1 + K);
213z = x(:) + 1i * y(:);
214col = 1;
215p = ones(np, 1);
216for mdeg = 0:M
217 A(:, col) = real(p);
218 col = col + 1;
219 if mdeg > 0
220 A(:, col) = imag(p);
221 col = col + 1;
222 end
223 p = p .* z;
224end
225th = atan2(y(:), x(:));
226th(th < -1e-13) = th(th < -1e-13) + 2*pi;
227r = abs(z);
228k = 0;
229q = 0;
230while q < K
231 k = k + 1;
232 if mod(k, 3) ~= 0
233 q = q + 1;
234 A(:, col) = r.^(2*k/3) .* sin(2*k*th/3);
235 col = col + 1;
236 end
237end
238end
240function [bx, by] = pac_bpoints(ntot)
241% Boundary points: on each straight edge, clustered quadratically towards
242% the re-entrant corner, where the data varies fastest; on the arc, spaced
243% uniformly in angle. Points per piece in proportion to its length.
244ne = max(6, round(ntot / 6.7));
245na = max(12, ntot - 2*ne);
246t = ((1:ne)' - 0.5) / ne;
247c = t.^2;
248s = 3*pi/2 * ((1:na)' - 0.5) / na;
249bx = [c; cos(s); zeros(ne, 1)];
250by = [zeros(ne, 1); sin(s); -c];
251end
253% --------------------------------------------- boundary integral, panels ----
255function v = bie_pac(d, g, px, py)
256% The double-layer representation and jump relation of Chapter 10, with the
257% trapezoid rule replaced by 12-point Gauss-Legendre panels, refined
258% dyadically to depth d towards each of the three corners. On the straight
259% edges the kernel between two points of the same edge is identically zero,
260% and on the arc it is the constant -1/(4*pi); everything nontrivial in this
261% matrix happens where two pieces meet at a corner.
262[qx, qy, qnx, qny, qw, qkap] = pac_panels(d);
263n = numel(qx);
1a8bb6aMake the cross-language benchmark measure the languagesJeremy Magland 265% The matrix is built by a double loop rather than by broadcasting the
266% coordinate differences into n-by-n arrays. The loop is the faster shape in
267% every environment this benchmark runs in: about twice as fast as implicit
268% expansion under numbl and in MATLAB, and several times faster than forming
269% the differences with outer products against ones(1,n), which allocates four
270% n-by-n temporaries and multiplies two of them into being. It is also the
271% only shape that writes each entry exactly once.
272A = zeros(n, n);
273c = 1 / (2*pi);
274for j = 1:n
275 xj = qx(j);
276 yj = qy(j);
277 nxj = qnx(j);
278 nyj = qny(j);
279 wj = qw(j);
280 for i = 1:n
281 dx = qx(i) - xj;
282 dy = qy(i) - yj;
283 A(i, j) = c * (nxj*dx + nyj*dy) / (dx*dx + dy*dy) * wj;
284 end
285 % the diagonal is the smooth limit of the kernel, -kappa/(4 pi), not a
286 % singularity, and the jump relation contributes the -1/2
287 A(j, j) = -0.5 - qkap(j) * wj / (4*pi);
288end
290sigma = A \ g(qx, qy);
1a8bb6aMake the cross-language benchmark measure the languagesJeremy Magland 292% the layer potential at the target points, by the same loop
293nt = numel(px);
294v = zeros(nt, 1);
295for k = 1:nt
296 xk = px(k);
297 yk = py(k);
298 s = 0;
299 for j = 1:n
300 dx = xk - qx(j);
301 dy = yk - qy(j);
302 s = s + (qnx(j)*dx + qny(j)*dy) / (dx*dx + dy*dy) * qw(j) * sigma(j);
303 end
304 v(k) = c * s;
305end
308function [qx, qy, qnx, qny, qw, qkap] = pac_panels(d)
309% Panel endpoints on each piece: dyadic towards both ends, with no panel
310% longer than 0.35, then 12-point Gauss-Legendre on every panel. All three
311% pieces are parametrised by arclength, so the speed is 1 throughout.
312[gx, gw] = lp_gauss(12);
313be = graded_breaks(1, d, 0.35);
314ba = graded_breaks(3*pi/2, d, 0.35);
316[t1, w1] = panel_nodes(be, gx, gw);
317[t2, w2] = panel_nodes(ba, gx, gw);
318[t3, w3] = panel_nodes(be, gx, gw);
320% edge theta = 0 from the corner to (1,0); the arc; the edge theta = 3*pi/2
321% from (0,-1) back towards the corner. Outward normals for a domain kept on
322% the left of this traversal.
323qx = [t1; cos(t2); zeros(numel(t3), 1)];
324qy = [zeros(numel(t1), 1); sin(t2); t3 - 1];
325qnx = [zeros(numel(t1), 1); cos(t2); ones(numel(t3), 1)];
326qny = [-ones(numel(t1), 1); sin(t2); zeros(numel(t3), 1)];
327qw = [w1; w2; w3];
328qkap = [zeros(numel(t1), 1); ones(numel(t2), 1); zeros(numel(t3), 1)];
329end
331function bp = graded_breaks(T, d, base)
332% Breakpoints on [0, T]: halve towards each end d times, then split anything
333% still longer than base into equal pieces.
334half = T / 2;
335left = [0, half * 2.^(-d:0)];
336bp0 = [left, T - left(end-1:-1:1)];
337bp = bp0(1);
338for k = 2:numel(bp0)
339 len = bp0(k) - bp0(k-1);
340 ns = max(1, ceil(len / base));
341 for q = 1:ns
342 bp = [bp, bp0(k-1) + len * q / ns];
343 end
344end
345end
347function [t, w] = panel_nodes(bp, gx, gw)
348np = numel(bp) - 1;
349p = numel(gx);
350t = zeros(np * p, 1);
351w = zeros(np * p, 1);
352for k = 1:np
353 a = bp(k);
354 b = bp(k+1);
355 t((k-1)*p+1 : k*p) = (a + b)/2 + (b - a)/2 * gx;
356 w((k-1)*p+1 : k*p) = (b - a)/2 * gw;
357end
358end
360% --------------------------------------- Chebyshev on the polar rectangle ---
362function v = spec_pac(N, g, px, py)
363% Collocation on the tensor Chebyshev grid of the rectangle
364% (r, theta) in [0,1] x [0, 3*pi/2], where Laplace's equation reads
365% r^2 u_rr + r u_r + u_tt = 0. The side r = 1 carries the arc data, the
366% sides theta = 0 and theta = 3*pi/2 the edge data, and the side r = 0 is
367% the corner itself, where u equals the corner value of the data.
368[D, s] = lp_cheb(N);
369r = (s + 1) / 2; % from 1 down to 0
370th = (s + 1) * (3*pi/4); % from 3*pi/2 down to 0
371Dr = 2 * D;
372Dt = (4 / (3*pi)) * D;
373Ar = diag(r.^2) * (Dr * Dr) + diag(r) * Dr;
374Bt = Dt * Dt;
376G = zeros(N+1, N+1);
377G(1, :) = g(cos(th'), sin(th')); % r = 1: the arc
378G(N+1, :) = g(0, 0) * ones(1, N+1); % r = 0: the corner
379G(:, 1) = g(zeros(N+1, 1), -r); % theta = 3*pi/2
380G(:, N+1) = g(r, zeros(N+1, 1)); % theta = 0
382ii = 2:N;
383nn = N - 1;
384A = Ar(ii, ii);
385Bi = Bt(ii, ii);
386F = -(Ar(ii, [1, N+1]) * G([1, N+1], ii)) - (G(ii, [1, N+1]) * Bt(ii, [1, N+1])');
388% the matrix decomposition of Chapter 7, with two different one-dimensional
389% operators instead of one
390[V, LA] = eig(A);
391[W, LB] = eig(Bi);
392la = diag(LA);
393lb = diag(LB);
394Ft = V \ (F / W.');
395Ut = Ft ./ (la * ones(1, nn) + ones(nn, 1) * lb.');
396U = G;
397U(ii, ii) = real(V * Ut * W.');
399% evaluate off the grid by barycentric interpolation in each direction
400rr = sqrt(px.^2 + py.^2);
401tt = atan2(py, px);
402tt(tt < -1e-13) = tt(tt < -1e-13) + 2*pi;
403v = bary2(U, s, 2*rr - 1, tt / (3*pi/4) - 1);
404end
406function v = bary2(U, s, sr, st)
407% Tensor barycentric interpolation on the Chebyshev grid s (both
408% directions), rows of U indexed by the first coordinate.
409n1 = numel(s);
410wb = [0.5; ones(n1-2, 1); 0.5] .* (-1).^((0:n1-1)');
411np = numel(sr);
412v = zeros(np, 1);
413for k = 1:np
414 dr = sr(k) - s;
415 dt = st(k) - s;
416 hr = abs(dr) < 1e-14;
417 ht = abs(dt) < 1e-14;
418 if any(hr)
419 cr = double(hr);
420 else
421 cr = wb ./ dr;
422 cr = cr / sum(cr);
423 end
424 if any(ht)
425 ct = double(ht);
426 else
427 ct = wb ./ dt;
428 ct = ct / sum(ct);
429 end
430 v(k) = cr' * U * ct;
431end
432end