088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 1% Two boundary methods for the interior Dirichlet problem on a smooth curve.
2%
3% 1. NYSTROM, double layer. Look for the solution as the potential of a dipole
4% density sigma spread over the curve,
5%
6% u(x) = (1/2pi) INT_Gamma (n_y . (x - y)) / |x - y|^2 sigma(y) ds(y),
7%
8% which is harmonic inside for any sigma. Letting x approach the curve from
9% inside picks up a jump, and the boundary condition becomes a Fredholm
10% equation of the second kind,
11%
12% -sigma(x)/2 + INT_Gamma K(x, y) sigma(y) ds(y) = g(x).
13%
14% Second kind is the important phrase: the operator is identity plus
15% compact, its eigenvalues cluster away from zero, and the discrete system
16% is well conditioned no matter how fine the discretisation. Compare the
17% finite difference matrices, whose condition number grows like N^2.
18%
19% The kernel K looks singular but is not: as y approaches x it tends to
20% -kappa(x)/(4pi), where kappa is the curvature. So the integrand is smooth
21% and periodic, and the plain trapezoid rule, which is spectrally accurate
22% for smooth periodic integrands, converges geometrically. That single fact
23% is what makes this method fast.
24%
25% 2. FUNDAMENTAL SOLUTIONS. Place point sources on a circle of radius par
26% outside the curve and fit their combined field to the boundary data by
27% least squares. Every basis function is harmonic, so the equation is
28% satisfied exactly and only the boundary condition is approximated. It is
29% three lines of code and it also converges geometrically, but the matrix is
30% exponentially ill conditioned in the number of sources, and the accuracy
31% it reaches depends on where the sources are put.
33function names = method_names()
34names = {'double-layer Nystrom', 'fundamental solutions'};
35end
37function v = solve_bie(mid, n, a, nlobe, par, px, py)
38if mid == 1
39 v = nystrom(n, a, nlobe, px, py);
40else
41 v = mfs(n, a, nlobe, par, px, py);
42end
43end
45function v = nystrom(n, a, nlobe, px, py)
46t = 2*pi*(0:n-1)' / n;
47h = 2*pi / n;
48C = lp_curve(a, nlobe, t);
49x = C(:,1);
50y = C(:,2);
51xp = C(:,3);
52yp = C(:,4);
53xpp = C(:,5);
54ypp = C(:,6);
56sp = sqrt(xp.^2 + yp.^2); % speed |gamma'(t)|
57nx = yp ./ sp; % outward normal, counterclockwise curve
58ny = -xp ./ sp;
59kappa = (xp .* ypp - yp .* xpp) ./ sp.^3;
61% the Nystrom matrix: A = -I/2 + h * K * diag(speed)
62DX = x * ones(1, n) - ones(n, 1) * x';
63DY = y * ones(1, n) - ones(n, 1) * y';
64R2 = DX.^2 + DY.^2;
65K = ((ones(n, 1) * nx') .* DX + (ones(n, 1) * ny') .* DY) ./ (2*pi*R2);
66K(1:n+1:end) = -kappa / (4*pi); % the diagonal limit, not a singularity
67A = -0.5*eye(n) + h * K .* (ones(n, 1) * sp');
69sigma = A \ lp_bexact(x, y);
71% evaluate the layer potential at the target points
72np = numel(px);
73EX = px(:) * ones(1, n) - ones(np, 1) * x';
74EY = py(:) * ones(1, n) - ones(np, 1) * y';
75KE = ((ones(np, 1) * nx') .* EX + (ones(np, 1) * ny') .* EY) ./ (2*pi*(EX.^2 + EY.^2));
76v = (h * KE .* (ones(np, 1) * sp')) * sigma;
77end
79function v = mfs(n, a, nlobe, par, px, py)
80% n sources on a circle of radius par, fitted at 3n boundary points.
81ts = 2*pi*(0:n-1)' / n;
82sx = par * cos(ts);
83sy = par * sin(ts);
85nb = 3 * n;
86tb = 2*pi*(0:nb-1)' / nb;
87C = lp_curve(a, nlobe, tb);
88bx = C(:,1);
89by = C(:,2);
91A = source_matrix(bx, by, sx, sy);
92c = A \ lp_bexact(bx, by);
93v = source_matrix(px(:), py(:), sx, sy) * c;
94end
96function A = source_matrix(x, y, sx, sy)
97% log|z - z_j| for each source, plus a constant column
98np = numel(x);
99n = numel(sx);
100A = zeros(np, n + 1);
101for j = 1:n
102 A(:, j) = 0.5 * log((x - sx(j)).^2 + (y - sy(j)).^2);
103end
104A(:, n+1) = 1;
105end