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
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 55S = cg_solve(U, mask, m, 1e-9, 4000); % U holds the boundary data, S the solution
57% the requested points are grid points, so sampling is a lookup
58v = zeros(numel(px), 1);
59for k = 1:numel(px)
60 i = round(py(k)/h) + N + 1;
61 j = round(px(k)/h) + N + 1;
62 v(k) = S(i, j);
63end
64end
66function U = cg_solve(U, mask, m, tol, maxit)
67R = neg_residual(U, mask, m);
68P = R;
69rs = dot2(R, R);
70stop = tol * tol * rs;
71for it = 1:maxit
72 AP = applyA(P, mask, m);
73 al = rs / dot2(P, AP);
74 U = U + al * P;
75 R = R - al * AP;
76 rsn = dot2(R, R);
77 if rsn < stop
78 break;
79 end
80 P = R + (rsn / rs) * P;
81 rs = rsn;
82end
83end
85function Y = applyA(V, mask, m)
86Y = zeros(m, m);
87for j = 2:m-1
88 for i = 2:m-1
89 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));
90 end
91end
92end
94function Y = neg_residual(V, mask, m)
95Y = zeros(m, m);
96for j = 2:m-1
97 for i = 2:m-1
98 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));
99 end
100end
101end
103function s = dot2(A, B)
104s = 0;
105n = numel(A);
106for k = 1:n
107 s = s + A(k) * B(k);
108end
109end
111% ------------------------------------------- the method of particular sols --
113function v = mps_l(n, g, px, py)
114M = n; % harmonic polynomials up to degree M
115K = n; % corner functions r^(2k/3) sin(2k theta/3)
117nb = (2*M + 1) + K;
118[bx, by] = boundary_points(max(8, ceil(nb / 2))); % about three points per unknown
119A = basis_matrix(bx, by, M, K);
120c = A \ g(bx, by); % least squares: many more points than unknowns
122v = basis_matrix(px, py, M, K) * c;
123end
125function A = basis_matrix(x, y, M, K)
126np = numel(x);
127nb = (2*M + 1) + numel(corner_orders(K));
128A = zeros(np, nb);
129z = x(:) + 1i*y(:);
130col = 1;
131p = ones(np, 1);
132for mdeg = 0:M
133 A(:, col) = real(p);
134 col = col + 1;
135 if mdeg > 0
136 A(:, col) = imag(p);
137 col = col + 1;
138 end
139 p = p .* z;
140end
141th = atan2(y(:), x(:));
142th(th < -1e-13) = th(th < -1e-13) + 2*pi;
143r = abs(z);
144ks = corner_orders(K);
145for q = 1:numel(ks)
146 k = ks(q);
147 A(:, col) = r.^(2*k/3) .* sin(2*k*th/3);
148 col = col + 1;
149end
150end
152function ks = corner_orders(K)
153% k = 1, 2, 4, 5, 7, ...: the multiples of 3 are harmonic polynomials already
154% in the basis, and including them twice would make the least-squares matrix
155% exactly rank deficient.
156ks = [];
157k = 1;
158while numel(ks) < K
159 if mod(k, 3) ~= 0
160 ks = [ks, k];
161 end
162 k = k + 1;
163end
164end
166function [bx, by] = boundary_points(nper)
167% Points along the six edges of the boundary, walking round from the origin:
168%
169% 1: y = 0, x from 0 to 1 the upper edge of the slit
170% 2: x = 1, y from 0 to 1
171% 3: y = 1, x from 1 to -1
172% 4: x = -1, y from 1 to -1
173% 5: y = -1, x from -1 to 0
174% 6: x = 0, y from -1 to 0 the left edge of the slit
175%
176% The two edges that meet at the re-entrant corner carry points clustered
177% towards it, where the data varies fastest.
178t = ((1:nper)' - 0.5) / nper;
179c = t .^ 2;
180z = zeros(nper, 1);
181o = ones(nper, 1);
182bx = [c; o; 1 - 2*t; -o; -1 + t; z];
183by = [z; t; o; 1 - 2*t; -o; -c];
184end