---
title: Method of Fundamental Solutions
url: https://www.emergentmind.com/topics/method-of-fundamental-solutions
type: topic
---

# Method of Fundamental Solutions

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$ [2404.07914, 2006.12763, 2009.13892, 2508.00515, 2409.04215].

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$ [2404.07914].

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 [2504.18261, 2508.00515, 2302.01551, 2106.04807].

## 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 [2404.07914, 2109.01886, 2408.00212]. 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$ [2404.07914, 2109.01886, 1006.0329, 2408.00212].
- **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 [2006.12763, 2006.15992].

## 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 [2109.01886, 2408.00212, 1006.0329]. 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 [2109.01886].

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$ [2408.00212].

### 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 [2504.09458].

## 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 [2006.12763, 2006.15992]. The solution ansatz then involves a linear combination of periodic fundamental solutions [2006.12763]:
\[
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 [2508.00515]. 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 [2504.18261, 2302.01551, 2106.04807]. 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 [2203.00996].

## 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 [2409.04215]:

- **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 [2409.04215].

## 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 [2109.01886, 2601.11041].
- **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 [2601.11041].
- **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) [2504.09458].

### 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 [2601.11041].
- **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 [2105.14811, 2601.11041].
- **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) [1909.02402, 2508.00515].

## 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 [2404.07914, 2508.00515, 2504.18261, 2302.01551].
- Multibody, multiphase, and multiphysics problems in complex, periodic, and doubly-periodic domains [2006.12763, 2409.04215].
- Large-scale particle, suspension, and interface-tracking simulations (e.g., evolving meshes, particle-fluid interaction) [2601.11041, 2409.04215].

### 7.2. Limitations

- Ill-conditioning is fundamental for classical MFS for large $N$ unless regularized [2109.01886, 2408.00212].
- Knowledge of analyticity radius/singularity location is often unavailable in complex domains; source placement may rely on a posteriori diagnostics [2404.07914].
- Corners and nonanalytic boundaries necessitate specialized source distributions or bases (e.g., graded, wavelet, or Whitney frames) [2504.09458].
- For exterior or unbounded domains, radiation, decay, or net-flux constraints require modified basis, completion flows, or global neutrality conditions [1006.0329, 2409.04215].

### 7.3. Current and Emerging Themes

- Robust well-conditioned MFS via harmonic basis orthogonalization, SVD/Arnoldi, and adaptive graded distributions [2109.01886, 2408.00212, 2504.09458].
- Automated or a posteriori source placement strategies guided by cross-validation or error indicators [2601.11041].
- MFS for nonlinear, time-dependent, and rarefied kinetic models (R13, CCR), including direct extension to 3D matrix-operator kernels [2504.18261, 2106.04807, 2302.01551].
- Large-scale parallelization and FMM-acceleration for multibody and suspension flows [2409.04215].

## References

For comprehensive methodology, convergence analysis, and recent stabilization frameworks: [2404.07914], [2109.01886], [2408.00212], [2504.09458], [1006.0329]. For extensions to periodic domains, elasticity, and multi-physics models: [2006.12763], [2006.15992], [2508.00515], [2504.18261], [2302.01551]. For applications to moving domains and mesh motion: [2601.11041]. For large-scale many-body and FMM-accelerated MFS: [2409.04215].

Source: https://www.emergentmind.com/topics/method-of-fundamental-solutions