1function out = solver(prob, n)
2% Nystrom discretization of the double-layer boundary integral equation
3% for the interior Dirichlet Laplace problem.
4%
5% The solution is represented as a double-layer potential
6% u(x) = (1/2pi) int_Gamma sigma(y) (x - y).n(y) / |x - y|^2 ds(y)
7% whose interior boundary limit gives the second-kind equation
8% (W - I/2) sigma = g.
9% W is discretized with the periodic trapezoid rule at n equispaced
10% parameter nodes; the kernel is smooth on a smooth curve, with the
11% diagonal limit -kappa(t) |x'(t)| / (4 pi). Convergence is geometric,
12% at a rate set by how far the data g continues analytically. No
13% correction is applied for targets close to the boundary, so the error
14% of the evaluated potential at a target a distance delta inside behaves
15% like exp(-2 pi delta / h) in the node spacing h (the close-evaluation
16% problem). On most instances the evaluation points stay a modest
17% distance inside and this costs nothing; on star-nearfield, whose
18% closest target is 0.005 in, it is what limits the solver.
19%
20% n : number of boundary quadrature nodes.
22h = 2*pi/n;
23t = h*(0:n-1)';
24xb = prob.curve(t);
25dxb = prob.curveD(t);
26ddxb = prob.curveDD(t);
27sp = sqrt(dxb(:,1).^2 + dxb(:,2).^2);
28nqx = dxb(:,2)./sp;
29nqy = -dxb(:,1)./sp;
31% M(i,j) = (1/2pi) (x_i - x_j).n(x_j) / |x_i - x_j|^2 * |x'(t_j)|
32dx = repmat(xb(:,1), 1, n) - repmat(xb(:,1)', n, 1);
33dy = repmat(xb(:,2), 1, n) - repmat(xb(:,2)', n, 1);
34r2 = dx.^2 + dy.^2;
35num = dx.*repmat(nqx', n, 1) + dy.*repmat(nqy', n, 1);
36M = (num./r2).*repmat(sp', n, 1)/(2*pi);
38% Diagonal limit: -kappa/2 * |x'| / (2 pi), kappa the signed curvature.
39kap = (dxb(:,1).*ddxb(:,2) - dxb(:,2).*ddxb(:,1))./sp.^3;
40md = -(kap/2).*sp/(2*pi);
41for i = 1:n
42 M(i, i) = md(i);
43end
45sigma = (h*M - 0.5*eye(n)) \ prob.g(t);
47out = struct();
48out.uEval = dlp_eval(prob.evalXY, xb, nqx, nqy, h*sp.*sigma);
49if size(prob.vizXY, 1) > 0
50 out.uGrid = dlp_eval(prob.vizXY, xb, nqx, nqy, h*sp.*sigma);
51else
52 out.uGrid = zeros(0, 1);
53end
55end
57function u = dlp_eval(XY, xb, nqx, nqy, w)
58% Evaluate the double-layer potential with combined weights w at the
59% rows of XY, in blocks to bound memory.
60m = size(XY, 1);
61nb = size(xb, 1);
62u = zeros(m, 1);
63B = 4000;
64for i0 = 1:B:m
65 i1 = min(i0 + B - 1, m);
66 mm = i1 - i0 + 1;
67 dx = repmat(XY(i0:i1, 1), 1, nb) - repmat(xb(:,1)', mm, 1);
68 dy = repmat(XY(i0:i1, 2), 1, nb) - repmat(xb(:,2)', mm, 1);
69 r2 = dx.^2 + dy.^2;
70 K = (dx.*repmat(nqx', mm, 1) + dy.*repmat(nqy', mm, 1))./r2/(2*pi);
71 u(i0:i1) = K*w;
72end
73end