1function out = solver(prob, n)
2% Method of fundamental solutions (MFS) for the interior Dirichlet
3% Laplace problem.
4%
5% The solution is represented as a sum of n logarithmic point charges
6% placed on a fictitious curve outside the domain: each charge sits at
7% the boundary point of parameter t_j pushed a fixed distance delta
8% along the outward normal. Charge strengths are found by collocation
9% of the Dirichlet data at the n boundary points t_j (a square dense
10% system). Convergence is geometric when the solution continues
11% harmonically past the charge curve; when the data has singularities
12% closer to the boundary than delta, the method stagnates. The system
13% grows exponentially ill-conditioned with n, which caps the attainable
14% accuracy near 1e-10.
15%
16% n : number of charges = number of collocation points.
18delta = 0.3;
20t = 2*pi*(0:n-1)'/n;
21xb = prob.curve(t);
22dxb = prob.curveD(t);
23sp = sqrt(dxb(:,1).^2 + dxb(:,2).^2);
24nx = dxb(:,2)./sp;
25ny = -dxb(:,1)./sp;
26qx = xb(:,1) + delta*nx;
27qy = xb(:,2) + delta*ny;
29A = kernel_matrix(xb(:,1), xb(:,2), qx, qy);
30coef = A \ prob.g(t);
32out = struct();
33out.uEval = apply_potential(prob.evalXY, qx, qy, coef);
34if size(prob.vizXY, 1) > 0
35 out.uGrid = apply_potential(prob.vizXY, qx, qy, coef);
36else
37 out.uGrid = zeros(0, 1);
38end
40end
42function A = kernel_matrix(px, py, qx, qy)
43% A(i,j) = log|p_i - q_j|
44m = numel(px);
45nq = numel(qx);
46dx = repmat(px, 1, nq) - repmat(qx', m, 1);
47dy = repmat(py, 1, nq) - repmat(qy', m, 1);
48A = 0.5*log(dx.^2 + dy.^2);
49end
51function u = apply_potential(XY, qx, qy, coef)
52% Evaluate the charge sum at the rows of XY, in blocks to bound memory.
53m = size(XY, 1);
54u = zeros(m, 1);
55B = 4000;
56for i0 = 1:B:m
57 i1 = min(i0 + B - 1, m);
58 u(i0:i1) = kernel_matrix(XY(i0:i1, 1), XY(i0:i1, 2), qx, qy)*coef;
59end
60end