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), entry by entry. A
62% double loop is about twice as fast here as the same expression written with
63% implicit expansion, and several times faster than forming the coordinate
64% differences with outer products against ones(1,n); it also writes each
65% entry once instead of passing over four n-by-n temporaries.
66A = zeros(n, n);
67c = h / (2*pi);
68for j = 1:n
69 xj = x(j);
70 yj = y(j);
71 nxj = nx(j);
72 nyj = ny(j);
73 sj = sp(j);
74 for i = 1:n
75 dx = x(i) - xj;
76 dy = y(i) - yj;
77 A(i, j) = c * (nxj*dx + nyj*dy) / (dx*dx + dy*dy) * sj;
78 end
79 % the diagonal limit, not a singularity, plus the jump relation's -1/2
80 A(j, j) = -0.5 - c * kappa(j) * sj / 2;
81end
83sigma = A \ lp_bexact(x, y);
85% evaluate the layer potential at the target points, by the same loop
86nt = numel(px);
87v = zeros(nt, 1);
88for k = 1:nt
89 xk = px(k);
90 yk = py(k);
91 s = 0;
92 for j = 1:n
93 dx = xk - x(j);
94 dy = yk - y(j);
95 s = s + (nx(j)*dx + ny(j)*dy) / (dx*dx + dy*dy) * sp(j) * sigma(j);
96 end
97 v(k) = c * s;
98end
99end
101function v = mfs(n, a, nlobe, par, px, py)
102% n sources on a circle of radius par, fitted at 3n boundary points.
103ts = 2*pi*(0:n-1)' / n;
104sx = par * cos(ts);
105sy = par * sin(ts);
107nb = 3 * n;
108tb = 2*pi*(0:nb-1)' / nb;
109C = lp_curve(a, nlobe, tb);
110bx = C(:,1);
111by = C(:,2);
113A = source_matrix(bx, by, sx, sy);
114c = A \ lp_bexact(bx, by);
115v = source_matrix(px(:), py(:), sx, sy) * c;
116end
118function A = source_matrix(x, y, sx, sy)
119% log|z - z_j| for each source, plus a constant column
120np = numel(x);
121n = numel(sx);
122A = zeros(np, n + 1);
123for j = 1:n
124 A(:, j) = 0.5 * log((x - sx(j)).^2 + (y - sy(j)).^2);
125end
126A(:, n+1) = 1;
127end