088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 1% Two ways to solve the Dirichlet problem on the L-shaped domain.
2%
3% 1. the five-point scheme on a uniform grid, solved by conjugate gradients
4% until the residual has fallen by 1e-9, which is far below the
5% discretisation error, so that what is measured is the discretisation
6% and not the solver;
7% 2. an expansion in harmonic functions, fitted to the boundary data in the
8% least-squares sense: the method of particular solutions.
9%
10% The second method's basis is the interesting part. Harmonic polynomials
11% Re z^m and Im z^m handle everything that is smooth up to the boundary, and
12% they alone would converge slowly, because the solution is not smooth at the
13% re-entrant corner. Adding the corner's own singular functions
14%
15% r^(2k/3) sin(2k theta / 3), k = 1, 2, 4, 5, 7, ...
16%
17% (the multiples of 3 are omitted because they are already harmonic
18% polynomials) puts the singularity into the basis instead of trying to
19% approximate it. This is the Fox-Henrici-Moler construction, and it is why a
20% few dozen basis functions can beat a grid of a hundred thousand points.
22function names = method_names()
23names = {'5-point + CG', 'harmonic basis, least squares'};
24end
26function v = solve_l(mid, n, g, px, py)
27if mid == 1
28 v = fd_l(n, g, px, py);
29else
30 v = mps_l(n, g, px, py);
31end
32end
34% --------------------------------------------------- finite differences -----
36function v = fd_l(N, g, px, py)
37% N intervals per unit length, so the grid is (-1 : 1/N : 1) in both
38% directions and every corner of the domain is a grid point.
39h = 1 / N;
40m = 2*N + 1;
41xs = (-N:N)' * h;
42[X, Y] = meshgrid(xs, xs);
44inL = ~(X > 1e-12 & Y < -1e-12); % the closed domain
45% a node is an unknown when it lies in the open domain: inside the square, in
46% the L, and not on either edge of the slit
47onslit = (abs(Y) < 1e-12 & X > -1e-12) | (abs(X) < 1e-12 & Y < 1e-12);
48unknown = inL & abs(X) < 1 - 1e-12 & abs(Y) < 1 - 1e-12 & ~onslit;
49mask = double(unknown);
51U = zeros(m, m);
52U(~unknown) = g(X(~unknown), Y(~unknown)); % boundary values
53U(~inL) = 0; % outside: never read
55% Note the fresh name: assigning the result of cg_solve back into U, which was
56% built here with zeros(), would give U two different inferred types in this
57% function, and numbl would then pass it to cg_solve as an untyped value and
58% run the whole solve about twenty times slower.
59S = cg_solve(U, mask, m, 1e-9, 4000);
61% the requested points are grid points, so sampling is a lookup
62v = zeros(numel(px), 1);
63for k = 1:numel(px)
64 i = round(py(k)/h) + N + 1;
65 j = round(px(k)/h) + N + 1;
66 v(k) = S(i, j);
67end
68end
70function U = cg_solve(U, mask, m, tol, maxit)
71R = neg_residual(U, mask, m);
72P = R;
73rs = dot2(R, R);
74stop = tol * tol * rs;
75for it = 1:maxit
76 AP = applyA(P, mask, m);
77 al = rs / dot2(P, AP);
78 U = U + al * P;
79 R = R - al * AP;
80 rsn = dot2(R, R);
81 if rsn < stop
82 break;
83 end
84 P = R + (rsn / rs) * P;
85 rs = rsn;
86end
87end
89function Y = applyA(V, mask, m)
90Y = zeros(m, m);
91for j = 2:m-1
92 for i = 2:m-1
93 Y(i,j) = mask(i,j) * (4*V(i,j) - V(i-1,j) - V(i+1,j) - V(i,j-1) - V(i,j+1));
94 end
95end
96end
98function Y = neg_residual(V, mask, m)
99Y = zeros(m, m);
100for j = 2:m-1
101 for i = 2:m-1
102 Y(i,j) = -mask(i,j) * (4*V(i,j) - V(i-1,j) - V(i+1,j) - V(i,j-1) - V(i,j+1));
103 end
104end
105end
107function s = dot2(A, B)
108s = 0;
109n = numel(A);
110for k = 1:n
111 s = s + A(k) * B(k);
112end
113end
115% ------------------------------------------- the method of particular sols --
117function v = mps_l(n, g, px, py)
118M = n; % harmonic polynomials up to degree M
119K = n; % corner functions r^(2k/3) sin(2k theta/3)
121nb = (2*M + 1) + K;
122[bx, by] = boundary_points(max(8, ceil(nb / 2))); % about three points per unknown
123A = basis_matrix(bx, by, M, K);
124c = A \ g(bx, by); % least squares: many more points than unknowns
126v = basis_matrix(px, py, M, K) * c;
127end
129function A = basis_matrix(x, y, M, K)
130np = numel(x);
131nb = (2*M + 1) + numel(corner_orders(K));
132A = zeros(np, nb);
133z = x(:) + 1i*y(:);
134col = 1;
135p = ones(np, 1);
136for mdeg = 0:M
137 A(:, col) = real(p);
138 col = col + 1;
139 if mdeg > 0
140 A(:, col) = imag(p);
141 col = col + 1;
142 end
143 p = p .* z;
144end
145th = atan2(y(:), x(:));
146th(th < -1e-13) = th(th < -1e-13) + 2*pi;
147r = abs(z);
148ks = corner_orders(K);
149for q = 1:numel(ks)
150 k = ks(q);
151 A(:, col) = r.^(2*k/3) .* sin(2*k*th/3);
152 col = col + 1;
153end
154end
156function ks = corner_orders(K)
157% k = 1, 2, 4, 5, 7, ...: the multiples of 3 are harmonic polynomials already
158% in the basis, and including them twice would make the least-squares matrix
159% exactly rank deficient.
160ks = [];
161k = 1;
162while numel(ks) < K
163 if mod(k, 3) ~= 0
164 ks = [ks, k];
165 end
166 k = k + 1;
167end
168end
170function [bx, by] = boundary_points(nper)
171% Points along the six edges of the boundary, walking round from the origin:
172%
173% 1: y = 0, x from 0 to 1 the upper edge of the slit
174% 2: x = 1, y from 0 to 1
175% 3: y = 1, x from 1 to -1
176% 4: x = -1, y from 1 to -1
177% 5: y = -1, x from -1 to 0
178% 6: x = 0, y from -1 to 0 the left edge of the slit
179%
180% The two edges that meet at the re-entrant corner carry points clustered
181% towards it, where the data varies fastest.
182t = ((1:nper)' - 0.5) / nper;
183c = t .^ 2;
184z = zeros(nper, 1);
185o = ones(nper, 1);
186bx = [c; o; 1 - 2*t; -o; -1 + t; z];
187by = [z; t; o; 1 - 2*t; -o; -c];
188end