1// Exact solution, geometry, and sampling rules for laplace-dirichlet-2d,
2// implemented independently of the MATLAB side (build_problem.m). The two
3// implementations of these simple formulas check each other: a solver can
4// only reach high accuracy if both agree.
6import type { Laplace2dInstance } from "./spec";
8export interface Source {
9 x: number;
10 y: number;
11 c: number;
12}
14/** Boundary radius r(t) = 1 + a cos(k t). */
15export function boundaryR(inst: Laplace2dInstance, t: number): number {
16 return 1 + inst.a * Math.cos(inst.k * t);
17}
19/** Boundary point at parameter t. */
20export function boundaryPoint(inst: Laplace2dInstance, t: number) {
21 const r = boundaryR(inst, t);
22 return { x: r * Math.cos(t), y: r * Math.sin(t) };
23}
25/** The three exact-solution sources: boundary points at phi_j pushed a
26 * distance d along the outward normal, strengths [1.0, -0.6, 0.8]. */
27export function sources(inst: Laplace2dInstance): Source[] {
28 const { a, k, d } = inst;
29 const strengths = [1.0, -0.6, 0.8];
30 return strengths.map((c, j) => {
31 const phi = (2 * Math.PI * j) / 3 + 0.4;
32 const r = 1 + a * Math.cos(k * phi);
33 const dr = -a * k * Math.sin(k * phi);
34 const bx = r * Math.cos(phi);
35 const by = r * Math.sin(phi);
36 const dx = dr * Math.cos(phi) - r * Math.sin(phi);
37 const dy = dr * Math.sin(phi) + r * Math.cos(phi);
38 const sp = Math.hypot(dx, dy);
39 return { x: bx + (d * dy) / sp, y: by - (d * dx) / sp, c };
40 });
41}
43/** Exact solution u(x, y) = sum_j c_j log|x - s_j|. */
44export function exactU(inst: Laplace2dInstance, x: number, y: number): number {
45 let u = 0;
46 for (const s of sources(inst)) {
47 u += s.c * 0.5 * Math.log((x - s.x) ** 2 + (y - s.y) ** 2);
48 }
49 return u;
50}
52/** The 289 evaluation points: 32 rays, radial fractions 0.1..0.9, plus
53 * the origin. Order matches build_problem.m: radius outer, angle inner,
54 * origin last. */
55export function evalPoints(inst: Laplace2dInstance): { x: number; y: number }[] {
56 const rho = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9];
57 const pts: { x: number; y: number }[] = [];
58 for (const r of rho) {
59 for (let j = 0; j < 32; j++) {
60 const th = (2 * Math.PI * j) / 32 + 0.13;
61 const rr = r * boundaryR(inst, th);
62 pts.push({ x: rr * Math.cos(th), y: rr * Math.sin(th) });
63 }
64 }
65 pts.push({ x: 0, y: 0 });
66 return pts;
67}
69/** Exact solution at the evaluation points. */
70export function exactAtEvalPoints(inst: Laplace2dInstance): Float64Array {
71 const pts = evalPoints(inst);
72 const u = new Float64Array(pts.length);
73 pts.forEach((p, i) => (u[i] = exactU(inst, p.x, p.y)));
74 return u;
75}
77/** The visualization grid: ngrid x ngrid points over [-R, R]^2 with
78 * R = 1.05 (1 + |a|). Flat index p = ix * ngrid + iy with x = xs[ix],
79 * y = xs[iy] (y varies fastest), matching build_problem.m's meshgrid
80 * column order. */
81export const VIZ_NGRID = 200;
83export function vizGrid(inst: Laplace2dInstance) {
84 const R = 1.05 * (1 + Math.abs(inst.a));
85 const xs = new Float64Array(VIZ_NGRID);
86 for (let i = 0; i < VIZ_NGRID; i++) {
87 xs[i] = -R + (2 * R * i) / (VIZ_NGRID - 1);
88 }
89 return { R, ngrid: VIZ_NGRID, xs };
90}
92/** Whether (x, y) is inside the domain (used only for display masking,
93 * never for scoring, so float tie-breaks at the boundary are harmless). */
94export function insideDomain(inst: Laplace2dInstance, x: number, y: number): boolean {
95 const rr = Math.hypot(x, y);
96 const th = Math.atan2(y, x);
97 return rr < boundaryR(inst, th);
98}
100/** Relative errors of numeric values against the exact solution at the
101 * evaluation points: max and L2, both relative to the exact values. */
102export function evalErrors(
103 inst: Laplace2dInstance,
104 uNum: ArrayLike<number>
105): { relMax: number; relL2: number } {
106 const uEx = exactAtEvalPoints(inst);
107 if (uNum.length !== uEx.length) {
108 throw new Error(`expected ${uEx.length} values, got ${uNum.length}`);
109 }
110 let maxDiff = 0;
111 let maxEx = 0;
112 let sumDiff2 = 0;
113 let sumEx2 = 0;
114 for (let i = 0; i < uEx.length; i++) {
115 const diff = Math.abs(uNum[i] - uEx[i]);
116 maxDiff = Math.max(maxDiff, diff);
117 maxEx = Math.max(maxEx, Math.abs(uEx[i]));
118 sumDiff2 += diff * diff;
119 sumEx2 += uEx[i] * uEx[i];
120 }
121 return {
122 relMax: maxDiff / maxEx,
123 relL2: Math.sqrt(sumDiff2 / sumEx2),
124 };
125}