Method of Fundamental Solutions

• Method of Fundamental Solutions is a meshless boundary collocation method that approximates linear PDE solutions using a superposition of fundamental solutions with sources placed outside the domain.
• It achieves spectral or exponential convergence by leveraging analytic continuation, though proper source placement and regularization are essential to counteract severe ill-conditioning.
• Modern applications include elliptic, periodic, elasticity, and multiphysics problems, enhanced by techniques like SVD-based stabilization and Fast Multipole Method acceleration.

The Method of Fundamental Solutions (MFS), also referred to as the Method of Auxiliary Sources (MAS), is a meshless boundary collocation approach for the numerical solution of linear partial differential equations (PDEs), especially elliptic equations such as Laplace, Helmholtz, modified Helmholtz, biharmonic, elasticity, and extended hydrodynamics (including the Grad and CCR models). MFS constructs the approximate solution by a linear combination of fundamental solutions with source points placed outside the physical domain—thus avoiding singular or weakly singular volume or boundary integrals. This technique achieves spectral or exponential convergence when the solution is analytic up to a suitable exterior region, but robust stabilization and conditioning are significant issues in high-accuracy or large-scale contexts.

1. Mathematical Formulation and Classical Framework

MFS seeks a solution to a linear PDE in a domain $\Omega \subset \mathbb{R}^d$ (e.g., $\Delta u = 0$ or $(-\Delta+\alpha^2)u=0$) by representing $u$ as a superposition of fundamental solutions

$u_N(x) = \sum_{j=1}^N a_j G(x, y_j)$

where $G(x, y)$ is the fundamental solution for the operator (Green's function with $\delta(x-y)$ source), and $\{y_j\}$ are source points placed on a “fictitious” boundary $\Gamma_{\mathrm{aux}}$ outside (or inside, for exterior or unbounded problems) of $\Omega$ (Kolezas et al., 2024, Ogata, 2020, Ei et al., 2020, Barceló et al., 1 Aug 2025, Broms et al., 2024).

Boundary or interface conditions (Dirichlet, Neumann, Robin, or transmission) are enforced at collocation points $\{x_i\}$ on $\partial\Omega$ or interfaces, leading to a linear system for the unknown amplitudes $\{a_j\}$. For, e.g., Laplace-Neumann problems in 2D, one enforces $\partial u_N/\partial n(x_i)=g(x_i)$ at $x_i\in\partial\Omega$, where

$$\frac{\partial}{\partial n_{x}} G(x, y) = \frac{n_{x} \cdot (x - y)}{2\pi|x-y|^2}$$

leading to an $N\times N$ linear system $A a = g$ (Kolezas et al., 2024).

This framework generalizes naturally to systems (e.g., elasticity, Stokes, R13, CCR), fully vector-valued amplitudes, and more general PDEs by using matrix (or tensor) Green’s kernels and enforcing appropriate coupling conditions (Himanshi et al., 25 Apr 2025, Barceló et al., 1 Aug 2025, Himanshi et al., 2023, Rana et al., 2021).

2. Source Placement, Ill-Conditioning, and Analyticity

The convergence and conditioning of MFS heavily depend on source placement:

• Source distance and analyticity: The optimal convergence rate is dictated by the analytic continuation region of the true solution. For bounded domains with analytic boundary/solution, the radius of source placement $\rho_{\mathrm{aux}}$ must avoid crossing any singularities of the analytic continuation of the solution (Kolezas et al., 2024, Antunes, 2021, Sakakibara, 2024). Crossing such singularities triggers divergence and oscillation of the amplitudes, but the field $u_N$ may still converge by low-pass filtering.
• Exponential ill-conditioning: As source points approach the boundary, the collocation matrix’s smallest singular value decays exponentially, and the condition number scales like $\sim (r_{\mathrm{aux}}/R_{\mathrm{crit}})^N$ for Laplace-type problems. This exponential ill-conditioning is intrinsic and cannot be eliminated by standard collocation, limiting the practically achievable $N$ (Kolezas et al., 2024, Antunes, 2021, Shigeta et al., 2010, Sakakibara, 2024).
• Optimal positioning: Best practices recommend placing sources just outside the analytic region (e.g., $10$–$30\%$ outside the largest inscribed domain for bounded problems), avoiding crossing singularities. For periodic and multiply-connected problems, periodic analogs (e.g., Weierstrass, theta, sigma functions) are used (Ogata, 2020, Ogata, 2020).

3. Recent Theoretical Advances: Stabilized and Well-Conditioned MFS

To address the ill-conditioning, several stabilized MFS formulations have been developed.

3.1. Harmonic Basis Expansion and SVD/Arnoldi Regularization

Expanding shifted fundamental solutions in a suitable harmonic function basis (e.g., Fourier, Spherical Harmonics, Trefftz) followed by orthogonalization (QR, SVD, Arnoldi) yields a numerically stable and well-conditioned basis (Antunes, 2021, Sakakibara, 2024, Shigeta et al., 2010). In 2D, for Laplace’s equation,

$$G(x, y_j) = -\frac{1}{2\pi} \log|x-y_j| \approx \sum_{m=0}^p c_m(y_j) h_m(x)$$

for $h_m$ harmonic polynomials. The coefficient matrix is then regularized (MFS-SVD, MFS-Arnoldi) to produce a final collocation system with uniformly bounded condition number $O(1)$ for all $N$, provided the sources enclose the domain (Antunes, 2021).

For dipole-type MFS (DSM) based on singularity-removed double-layer potentials, the expansion contains additional cancellation structure, and QR regularization yields a uniformly well-conditioned system, even for high $N$ (Sakakibara, 2024).

3.2. Frames and Lusin Wavelet Bases for Lipschitz/Non-analytic Domains

For boundaries lacking analyticity or for singular solutions, MFS with source points distributed in a graded “Whitney” fashion and Lusin wavelet basis functions forms a frame for the appropriate Hardy space. This ensures numerical stability and spectral or high-order convergence, even in the presence of corners and singularities (Jonsson et al., 13 Apr 2025).

4. MFS for Complex, Periodic, and Multiphysics Problems

MFS extends flexibly to a broad array of model classes:

4.1. Periodic and Quasiperiodic Domains

For doubly-periodic potential problems, standard free-space Green’s functions are replaced by periodic analogs built from Jacobi theta-functions or Weierstrass sigma functions, ensuring the required periodicity of the solution and its derivatives (Ogata, 2020, Ogata, 2020). The solution ansatz then involves a linear combination of periodic fundamental solutions (Ogata, 2020): $$f_N(z) = U z - \frac{i}{2\pi} \sum_{j=1}^N Q_j \left\{ \log \vartheta_1\left(\frac{z-\zeta_j}{\omega_1} | \tau \right) - u_j z \right\}$$ with charge-neutrality and period conditions enforced in the linear system.

4.2. Vector PDEs, Elasticity, Hydrodynamics, and Extended Models

For the Navier–Lamé system (linear elasticity), the fundamental solution is a matrix-valued kernel with both P- and S-wave components. A recent addition theorem yields an explicit expansion solely in Bessel functions and scalar spherical harmonics, facilitating efficient MFS implementation in 3D with optimal complexity (Barceló et al., 1 Aug 2025). The MFS ansatz becomes

$u(x) \approx \sum_{j=1}^N \Phi(x, y_j) \alpha_j$

where $\Phi$ is the 3D elastic Green’s tensor.

For extended hydrodynamics (Grad’s 13-moment system, CCR theory) and rarefied gas dynamics, matrix-valued fundamental solutions are derived by Fourier analysis and then used as MFS kernels for boundary-driven microflows and phase-change problems (Himanshi et al., 25 Apr 2025, Himanshi et al., 2023, Rana et al., 2021). The approach supports rarefaction, slip, evaporation, temperature jump, and other kinetic interface effects.

4.3. Time-Dependent and Nonstandard Constraints

For timedomain wave and scattering problems, the MFS is coupled with convolution quadrature (CQ), discretizing time integrals via multistep methods. Modified CQ schemes further reduce dispersion/dissipation artifacts in wave propagation, using kernel shifts aligned with physical propagation speed (Aldahham et al., 2022).

5. Large-Scale, High-Dimensional, and Multibody Extensions

Recent scalable MFS variants enable solution of large-scale 3D elastance, mobility, and boundary value problems for systems of thousands of particles or conductors (Broms et al., 2024):

• Proxy source strategy: For each body or object, sources are placed on an “inner proxy” surface (e.g., slightly shrunk copy), and the degrees of freedom include both physical and auxiliary (completion) flows.
• Blockwise SVD preconditioning: Each object's "self" matrix is inverted via SVD once, and block-diagonal preconditioning yields global systems with identity diagonal blocks.
• Fast Multipole Method (FMM) acceleration: All matrix-vector products in GMRES iterations are performed with FMM, yielding $O(PN)$ cost for $P$ particles and $N$ proxies per particle.
• Constraint incorporation: Charges, forces, or torques are imposed via constraint orthogonal projections, guaranteeing physical net quantities.

This formulation achieves $O(N)$ complexity and enables solution of up to $10^7$ degrees of freedom on workstations with accuracy competitive with boundary integral methods (Broms et al., 2024).

6. Implementation, Conditioning, and Best Practices

6.1. Regularization and Error Indicators

• Conditioning management: Well-conditioned bases (via SVD or harmonic expansion), adaptive least-squares collocation (zero-padded), or regularization (Tikhonov, hat-matrix LOOCV) are used to mitigate roundoff errors (Antunes, 2021, Ammad et al., 16 Jan 2026).
• A posteriori error estimation: LOOCV using the hat matrix or direct use of the maximum boundary residual under the maximum principle provides robust error control (Ammad et al., 16 Jan 2026).
• Analytic or singular data: For solutions with analytic boundary data, MFS achieves exponential convergence. For nonanalytic data or corners, algebraic rates are observed but can be improved by local clustering or graded distributions (“Whitney” MFS, corner cones) (Jonsson et al., 13 Apr 2025).

6.2. Collocation, Overdetermination, and Geometric Adaptivity

• Collocation: Overdetermined least-squares collocation increases robustness without changing the source set; zero-padded schemes further regularize highly nonconvex or singularly perturbed scenarios (Ammad et al., 16 Jan 2026).
• Source/direct point pairing: For robustness, sources are often paired to collocation points by geometric proximity, normal direction, or conformal mapping, adapting to complex, moving, or time-dependent geometries (Sakakibara et al., 2021, Ammad et al., 16 Jan 2026).
• Boundary data and interface enforcement: For fluid–solid, multi-physics, or transmission problems, vector-valued collocation and multi-constraint interface systems are directly supported (e.g., two-phase or elastic transmission) (Kleefeld et al., 2019, Barceló et al., 1 Aug 2025).

7. Applications, Limitations, and Current Research Directions

7.1. Applicability

MFS provides high-accuracy meshless solutions for:

• Interior/exterior Laplace, Helmholtz, modified Helmholtz, biharmonic, Navier–Lamé, Stokes, and extended hydrodynamics PDEs (Kolezas et al., 2024, Barceló et al., 1 Aug 2025, Himanshi et al., 25 Apr 2025, Himanshi et al., 2023).
• Multibody, multiphase, and multiphysics problems in complex, periodic, and doubly-periodic domains (Ogata, 2020, Broms et al., 2024).
• Large-scale particle, suspension, and interface-tracking simulations (e.g., evolving meshes, particle-fluid interaction) (Ammad et al., 16 Jan 2026, Broms et al., 2024).

7.2. Limitations

• Ill-conditioning is fundamental for classical MFS for large $N$ unless regularized (Antunes, 2021, Sakakibara, 2024).
• Knowledge of analyticity radius/singularity location is often unavailable in complex domains; source placement may rely on a posteriori diagnostics (Kolezas et al., 2024).
• Corners and nonanalytic boundaries necessitate specialized source distributions or bases (e.g., graded, wavelet, or Whitney frames) (Jonsson et al., 13 Apr 2025).
• For exterior or unbounded domains, radiation, decay, or net-flux constraints require modified basis, completion flows, or global neutrality conditions (Shigeta et al., 2010, Broms et al., 2024).

7.3. Current and Emerging Themes

• Robust well-conditioned MFS via harmonic basis orthogonalization, SVD/Arnoldi, and adaptive graded distributions (Antunes, 2021, Sakakibara, 2024, Jonsson et al., 13 Apr 2025).
• Automated or a posteriori source placement strategies guided by cross-validation or error indicators (Ammad et al., 16 Jan 2026).
• MFS for nonlinear, time-dependent, and rarefied kinetic models (R13, CCR), including direct extension to 3D matrix-operator kernels (Himanshi et al., 25 Apr 2025, Rana et al., 2021, Himanshi et al., 2023).
• Large-scale parallelization and FMM-acceleration for multibody and suspension flows (Broms et al., 2024).

For comprehensive methodology, convergence analysis, and recent stabilization frameworks: (Kolezas et al., 2024, Antunes, 2021, Sakakibara, 2024, Jonsson et al., 13 Apr 2025, Shigeta et al., 2010). For extensions to periodic domains, elasticity, and multi-physics models: (Ogata, 2020, Ogata, 2020, Barceló et al., 1 Aug 2025, Himanshi et al., 25 Apr 2025, Himanshi et al., 2023). For applications to moving domains and mesh motion: (Ammad et al., 16 Jan 2026). For large-scale many-body and FMM-accelerated MFS: (Broms et al., 2024).