Rational Approximation with blahpack

The Laplace solver follows [C20], decomposing the solution into a smooth polynomial part and a singular rational part:

where is the domain centroid, are poles outside the domain, and is a constant that zeros the imaginary part at . The physical solution is .

Step 1: Polynomial fit. Given boundary samples , fit complex polynomial coefficients via real least squares. Since and , this separates into a real overdetermined system of size . We solve this with Householder QR, which is the standard approach for dense least squares. MATLAB’s \ operator uses the same method internally for overdetermined systems.

Step 2: AAA rational approximation. The boundary residual is approximated by the AAA algorithm [NST18], which builds a rational function in barycentric form:

where are greedily selected support points, , and are barycentric weights computed via SVD. At each iteration, the algorithm selects the point of maximum error, forms a Loewner matrix, and takes its smallest right singular vector as the new weights.

Our implementation vs. MATLAB. The reference MATLAB code calls svd(A, 0) — a thin SVD via LAPACK’s dgesdd (divide-and-conquer). Our JavaScript implementation calls zgesvd (complex bidiagonal QR iteration), translated from the reference Fortran LAPACK into JavaScript as part of the blahpack project. The key difference is algorithmic: dgesdd is typically faster for large matrices but zgesvd is simpler and equally accurate. Both compute the same singular values and vectors to machine precision.

We use the complex SVD (zgesvd) even though the AAA algorithm as stated in [NST18] operates on real-valued functions. This is because the Loewner matrix is complex (the sample points lie in ), so the SVD must handle complex entries. MATLAB’s svd dispatches to the complex routine automatically; our code calls zgesvd explicitly.

Step 3: Pole extraction. Poles and zeros are extracted from the barycentric form by solving a generalized eigenvalue problem on a companion-like matrix pencil of size , using the QZ algorithm (zggev). This is identical to the MATLAB approach (eig(E, B)), again translated from reference LAPACK Fortran.

Step 4: Pole filtering. Only poles exterior to the domain are retained (tested via ray-casting point-in-polygon). Interior poles would create non-physical singularities in the harmonic solution.

Step 5: Singular coefficient fit. With filtered exterior poles, fit complex coefficients and constant via least squares on the Cauchy system . This is again a real overdetermined system solved by Householder QR.

GPU evaluation. The final solution is evaluated per-pixel in a WebGPU fragment shader using f32 arithmetic. The polynomial part is evaluated in monomial form (18 terms) and the singular part as a direct partial-fraction sum (40–90 terms depending on mmax). At mmax ≥ 200, the singular coefficients grow to with massive cancellation, exceeding f32’s ~7-digit precision. At mmax = 100, the evaluation is clean.