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. Accuracy
13% of the evaluated potential degrades for targets very close to the
14% boundary (the close-evaluation problem); the evaluation points of
15% this problem stay a modest distance inside.
16%
17% n : number of boundary quadrature nodes.
19h = 2*pi/n;
20t = h*(0:n-1)';
21xb = prob.curve(t);
22dxb = prob.curveD(t);
23ddxb = prob.curveDD(t);
24sp = sqrt(dxb(:,1).^2 + dxb(:,2).^2);
25nqx = dxb(:,2)./sp;
26nqy = -dxb(:,1)./sp;
28% M(i,j) = (1/2pi) (x_i - x_j).n(x_j) / |x_i - x_j|^2 * |x'(t_j)|
29dx = repmat(xb(:,1), 1, n) - repmat(xb(:,1)', n, 1);
30dy = repmat(xb(:,2), 1, n) - repmat(xb(:,2)', n, 1);
31r2 = dx.^2 + dy.^2;
32num = dx.*repmat(nqx', n, 1) + dy.*repmat(nqy', n, 1);
33M = (num./r2).*repmat(sp', n, 1)/(2*pi);
35% Diagonal limit: -kappa/2 * |x'| / (2 pi), kappa the signed curvature.
36kap = (dxb(:,1).*ddxb(:,2) - dxb(:,2).*ddxb(:,1))./sp.^3;
37md = -(kap/2).*sp/(2*pi);
38for i = 1:n
39 M(i, i) = md(i);
40end
42sigma = (h*M - 0.5*eye(n)) \ prob.g(t);
44out = struct();
45out.uEval = dlp_eval(prob.evalXY, xb, nqx, nqy, h*sp.*sigma);
46if size(prob.vizXY, 1) > 0
47 out.uGrid = dlp_eval(prob.vizXY, xb, nqx, nqy, h*sp.*sigma);
48else
49 out.uGrid = zeros(0, 1);
50end
52end
54function u = dlp_eval(XY, xb, nqx, nqy, w)
55% Evaluate the double-layer potential with combined weights w at the
56% rows of XY, in blocks to bound memory.
57m = size(XY, 1);
58nb = size(xb, 1);
59u = zeros(m, 1);
60B = 4000;
61for i0 = 1:B:m
62 i1 = min(i0 + B - 1, m);
63 mm = i1 - i0 + 1;
64 dx = repmat(XY(i0:i1, 1), 1, nb) - repmat(xb(:,1)', mm, 1);
65 dy = repmat(XY(i0:i1, 2), 1, nb) - repmat(xb(:,2)', mm, 1);
66 r2 = dx.^2 + dy.^2;
67 K = (dx.*repmat(nqx', mm, 1) + dy.*repmat(nqy', mm, 1))./r2/(2*pi);
68 u(i0:i1) = K*w;
69end
70end