Kernel-Free Boundary Integral Method

• Kernel-Free Boundary Integral Method is a numerical technique that reformulates boundary value and interface problems as auxiliary PDE solves on structured grids.
• The method replaces explicit singular quadrature of layer and volume potentials with fast, grid-based solvers, preserving second-kind conditioning.
• KFBI has been successfully applied to bidomain, Stokes, Brinkman, and surface PDEs, achieving second-order accuracy and reduced computational costs.

Kernel-free boundary integral method (KFBI) denotes a class of potential-theory-based methods in which a boundary value or interface problem is first reformulated as a boundary integral equation, but the boundary and volume integrals are not evaluated by explicit singular-kernel quadrature. Instead, the action of layer and volume potentials is obtained from structured-grid solutions of equivalent interface problems, followed by interpolation or limiting to boundary or interface nodes. In this sense, the explicit analytical expression of the kernel function is not required when solving the boundary integral equations, even though Green’s functions remain part of the underlying theory (Gao et al., 2021, Yin et al., 22 Aug 2025). Recent formulations span bidomain equations, Stokes and Navier systems, two-phase Stokes and Brinkman interface problems, elliptic PDEs with implicitly defined interfaces, and elliptic equations posed on surfaces (Zhao et al., 2023, Dong et al., 2023).

1. Historical setting and conceptual scope

A useful prehistory appears in work on the relationship between boundary integral equations and radial basis functions. One 2002 survey established the RBF on numerical integration analysis based on an intrinsic relationship between the Green’s boundary integral representation and RBF, stated that the kernel function of integral equation is important to create efficient RBF, introduced the fundamental solution RBF as a strategy for constructing operator-dependent RBF, and presented the boundary knot method as a boundary-only, meshless, spectral convergent, integration-free RBF collocation technique [0207016]. Although this was not yet the later Cartesian-grid KFBI framework, it suggests an early shift from treating the kernel only as an object of singular quadrature to treating it as a constructive device for discretization design.

In the later KFBI literature, the central move is more specific: classical layer and volume potentials are retained at the analytical level, but their numerical evaluation is replaced by auxiliary PDE solves on a simple embedding box or on a parameter domain. The resulting methods preserve the dimension-reduction logic of boundary integral formulations while exploiting fast structured-grid solvers such as FFT-based Poisson solvers and geometric multigrid (Gao et al., 2021, Zhao et al., 2023).

The phrase “kernel-free” is not entirely uniform across adjacent literatures. In electromagnetic scattering from a perfect electric conductor, a field-only surface integral formulation was derived in which the surface integral equations do not contain divergent kernels, and all integrals can therefore be evaluated by standard quadrature after an analytical subtraction using auxiliary Helmholtz solutions (Sun et al., 2019). This is kernel-free in a nonsingular-integrand sense, whereas the dominant KFBI usage in elliptic and interface PDEs is kernel-free in the sense of replacing explicit kernel evaluation by equivalent interface solves.

2. Boundary-integral reformulations

KFBI inherits the standard potential-theoretic reformulation of elliptic problems as second-kind boundary integral equations. For scalar elliptic Dirichlet and Neumann problems with

$\nabla\cdot(\sigma\nabla u)-\kappa u=f,$

one formulation introduces the double-layer, single-layer, adjoint-double-layer, hypersingular, and volume operators and yields

$$\Bigl(\tfrac12+\mathcal K\Bigr)\varphi=g_D-\mathcal Gf, \qquad \Bigl(\tfrac12-\mathcal K'\Bigr)\psi=g_N-\sigma\,\partial_n\mathcal Gf,$$

with $\varphi=u|_\Gamma$ and $\psi=\sigma\partial_nu|_\Gamma$ (Zhou et al., 2023). For interface problems, two unknown boundary densities can be introduced to obtain a well-conditioned system of the second kind (Zhou et al., 2023).

The same pattern extends to systems. For the bidomain equations, after temporal discretization of the diffusion part, the homogeneous Neumann boundary value problem is embedded into a larger box, and the solution is represented as

$$u(p)=\int_\Omega G(q;p)\,f(q)\,dq-\int_{\partial\Omega}G(q;p)\,\psi(q)\,ds_q.$$

Enforcing the Neumann condition gives a second-kind Fredholm boundary integral equation for the unknown density $\psi$ (Gao et al., 2021).

For three-dimensional Stokes and Navier boundary value problems on irregular domains, the unknown boundary densities $\phi_s$ and $\phi_e$ satisfy

$$\tfrac12\,\phi_s-\mathcal M_v[\phi_s]=g_s-\mathcal G_v[f_s], \qquad \tfrac12\,\phi_e-\mathcal M[\phi_e]=g_e-\mathcal G[f_e],$$

and the corresponding discrete forms are described as well-conditioned and solved by GMRES (Zhao et al., 2023).

Interface formulations with discontinuous coefficients are especially prominent. For Brinkman- and Stokes-type interface problems with piecewise-constant coefficients, the original interface problem is recast as boundary integral equations in which the integral operators are interpreted as boundary data for potential functions satisfying simpler interface problems without coefficient discontinuities (Zhou et al., 16 Apr 2026). On surfaces, elliptic boundary value and interface problems admit analogous layer-potential representations in terms of the surface Green’s function $G(p,q)$, including

$$\Bigl(\tfrac12+\mathcal K\Bigr)\varphi=g_D-\mathcal Gf, \qquad \Bigl(\tfrac12-\mathcal K'\Bigr)\psi=g_N-\sigma\,\partial_n\mathcal Gf,$$0

for the Dirichlet problem and, in the special case $$\Bigl(\tfrac12+\mathcal K\Bigr)\varphi=g_D-\mathcal Gf, \qquad \Bigl(\tfrac12-\mathcal K'\Bigr)\psi=g_N-\sigma\,\partial_n\mathcal Gf,$$1,

$$\Bigl(\tfrac12+\mathcal K\Bigr)\varphi=g_D-\mathcal Gf, \qquad \Bigl(\tfrac12-\mathcal K'\Bigr)\psi=g_N-\sigma\,\partial_n\mathcal Gf,$$2

for interface problems (Yin et al., 22 Aug 2025).

3. Equivalent interface problems as the kernel-free mechanism

The defining numerical idea of KFBI is that each boundary or volume potential is realized as the solution of an auxiliary PDE on a simple computational domain. For the bidomain equations, the volume integral and boundary integral are replaced by two interface PDEs on the embedding box $$\Bigl(\tfrac12+\mathcal K\Bigr)\varphi=g_D-\mathcal Gf, \qquad \Bigl(\tfrac12-\mathcal K'\Bigr)\psi=g_N-\sigma\,\partial_n\mathcal Gf,$$3: a volume-integral problem with continuous solution and flux across $$\Bigl(\tfrac12+\mathcal K\Bigr)\varphi=g_D-\mathcal Gf, \qquad \Bigl(\tfrac12-\mathcal K'\Bigr)\psi=g_N-\sigma\,\partial_n\mathcal Gf,$$4, and a boundary-integral problem with a prescribed flux jump $$\Bigl(\tfrac12+\mathcal K\Bigr)\varphi=g_D-\mathcal Gf, \qquad \Bigl(\tfrac12-\mathcal K'\Bigr)\psi=g_N-\sigma\,\partial_n\mathcal Gf,$$5 on $$\Bigl(\tfrac12+\mathcal K\Bigr)\varphi=g_D-\mathcal Gf, \qquad \Bigl(\tfrac12-\mathcal K'\Bigr)\psi=g_N-\sigma\,\partial_n\mathcal Gf,$$6. The required normal derivatives at quasi-uniform interface nodes are then obtained by interpolation of the discrete Cartesian-grid solution (Gao et al., 2021).

For 3D Stokes and Navier equations, every volume potential and every double-layer potential is replaced by an interface problem on a larger cube $$\Bigl(\tfrac12+\mathcal K\Bigr)\varphi=g_D-\mathcal Gf, \qquad \Bigl(\tfrac12-\mathcal K'\Bigr)\psi=g_N-\sigma\,\partial_n\mathcal Gf,$$7. In the Stokes case, the unified auxiliary problem has the form

$$\Bigl(\tfrac12+\mathcal K\Bigr)\varphi=g_D-\mathcal Gf, \qquad \Bigl(\tfrac12-\mathcal K'\Bigr)\psi=g_N-\sigma\,\partial_n\mathcal Gf,$$8

with jump data

$$\Bigl(\tfrac12+\mathcal K\Bigr)\varphi=g_D-\mathcal Gf, \qquad \Bigl(\tfrac12-\mathcal K'\Bigr)\psi=g_N-\sigma\,\partial_n\mathcal Gf,$$9

and homogeneous Dirichlet data on $\varphi=u|_\Gamma$0. Boundary integral data are then obtained by solving on a uniform staggered Cartesian grid and interpolating the discrete solution to control points on $\varphi=u|_\Gamma$1 (Zhao et al., 2023).

The two-phase Stokes formulation makes the same correspondence explicit through three special cases of a constant-viscosity interface problem:

• volume integrals by $\varphi=u|_\Gamma$2,
• single-layer potentials by $\varphi=u|_\Gamma$3,
• double-layer potentials by $\varphi=u|_\Gamma$4 (Dong et al., 2023).

On surfaces, the same idea survives after pullback to a planar parameter domain. A surface PDE

$\varphi=u|_\Gamma$5

with jumps on a curve $\varphi=u|_\Gamma$6 is rewritten on $\varphi=u|_\Gamma$7 through a parametrization $\varphi=u|_\Gamma$8, producing a self-adjoint interface operator in local coordinates. The trace $\varphi=u|_\Gamma$9 and the conormal derivative $\psi=\sigma\partial_nu|_\Gamma$0 are then recovered by local interpolation of the grid solution (Yin et al., 22 Aug 2025).

This kernel-free mechanism removes the need for singular or hypersingular quadrature at the implementation level. In the 3D Stokes/Navier formulation, no extra work or special quadratures are required to deal with singular or hyper-singular boundary integrals, and the dependence on the analytical expressions of Green’s functions for the integral kernels is completely eliminated (Zhao et al., 2023).

4. Discretization, correction functions, and fast solvers

Most KFBI implementations use structured Cartesian discretizations with local interface corrections. For incompressible flow problems, the standard choice is a modified marker-and-cell scheme on staggered grids. In two-phase Stokes with discontinuous viscosity, the constant-viscosity auxiliary interface problems are discretized by a modified MAC finite-difference scheme, with sparse jump corrections supported only on irregular nodes adjacent to the interface. Because the coefficient matrix is identical to that of a uniform Stokes solve, FFT-based preconditioned CG solvers remain applicable, with cost $\psi=\sigma\partial_nu|_\Gamma$1 per MAC solve (Dong et al., 2023). In 3D Stokes/Navier, an augmented saddle-point system is reduced to an SPD Schur-complement problem for pressure, solved by preconditioned CG, with each iteration requiring FFT Poisson solves and total cost $\psi=\sigma\partial_nu|_\Gamma$2 per interface-problem solve (Zhao et al., 2023).

For scalar elliptic problems, correction-function variants generalize the local-interface treatment. One formulation defines a correction function

$\psi=\sigma\partial_nu|_\Gamma$3

inside a narrow band around the interface and derives a local Cauchy problem

$\psi=\sigma\partial_nu|_\Gamma$4

A mesh-free collocation method is then used in local patches to compute $\psi=\sigma\partial_nu|_\Gamma$5, and the method explicitly avoids complicated derivation for derivative jumps of the solution and is easy to implement, especially for the fourth-order method in three space dimensions (Zhou et al., 2023).

The Brinkman CF-KFBI formulation uses the same philosophy in a vector setting. Within a narrow band around the interface, a local correction function represents the solution jump and leads to a local Cauchy problem solved by collocation. In the 2D polynomial model described there, $\psi=\sigma\partial_nu|_\Gamma$6 is approximated in $\psi=\sigma\partial_nu|_\Gamma$7, $\psi=\sigma\partial_nu|_\Gamma$8 in $\psi=\sigma\partial_nu|_\Gamma$9, and there are $$u(p)=\int_\Omega G(q;p)\,f(q)\,dq-\int_{\partial\Omega}G(q;p)\,\psi(q)\,ds_q.$$0 unknowns. A minimal set of $$u(p)=\int_\Omega G(q;p)\,f(q)\,dq-\int_{\partial\Omega}G(q;p)\,\psi(q)\,ds_q.$$1 collocation points is specified, and solvability is proved by a perturbation-of-Vandermonde argument (Zhou et al., 16 Apr 2026).

On surfaces, second-order corrected finite differences are built for the pulled-back self-adjoint operator, and the correction in each small ball is represented in the basis $$u(p)=\int_\Omega G(q;p)\,f(q)\,dq-\int_{\partial\Omega}G(q;p)\,\psi(q)\,ds_q.$$2. A $$u(p)=\int_\Omega G(q;p)\,f(q)\,dq-\int_{\partial\Omega}G(q;p)\,\psi(q)\,ds_q.$$3 local system is solved and then patched by a partition of unity. The resulting global linear system is sparse, symmetric positive-definite, and solved by geometric multigrid with full-multigrid initialization (Yin et al., 22 Aug 2025).

At the boundary-integral level, KFBI methods typically use matrix-free GMRES; the bidomain work also describes simple Richardson iteration, while the correction-function variants continue to use GMRES with fast PDE-based operator application (Gao et al., 2021, Zhou et al., 2023).

5. Reported accuracy, robustness, and application domains

Problem setting	Reported numerical behavior	Source
Bidomain equations in 2D and 3D	Second-order spatial accuracy; GMRES $$u(p)=\int_\Omega G(q;p)\,f(q)\,dq-\int_{\partial\Omega}G(q;p)\,\psi(q)\,ds_q.$$4–$$u(p)=\int_\Omega G(q;p)\,f(q)\,dq-\int_{\partial\Omega}G(q;p)\,\psi(q)\,ds_q.$$5 in 2D; Richardson $$u(p)=\int_\Omega G(q;p)\,f(q)\,dq-\int_{\partial\Omega}G(q;p)\,\psi(q)\,ds_q.$$6–$$u(p)=\int_\Omega G(q;p)\,f(q)\,dq-\int_{\partial\Omega}G(q;p)\,\psi(q)\,ds_q.$$7 in 3D; real left-ventricle simulation in a $$u(p)=\int_\Omega G(q;p)\,f(q)\,dq-\int_{\partial\Omega}G(q;p)\,\psi(q)\,ds_q.$$8 box	(Gao et al., 2021)
3D Stokes and Navier on irregular domains	Second-order convergence in $$u(p)=\int_\Omega G(q;p)\,f(q)\,dq-\int_{\partial\Omega}G(q;p)\,\psi(q)\,ds_q.$$9 and, for Stokes, in $\psi$0; GMRES and CG counts essentially constant; total CPU grows like $\psi$1	(Zhao et al., 2023)
Two-phase Stokes with discontinuous viscosity	Second order in discrete $\psi$2-norm for velocity, pressure, and velocity gradient, and second order in maximum norm for velocity and its gradient, even for high contrast viscosity	(Dong et al., 2023)
Elliptic PDEs with implicitly defined interfaces	2D Dirichlet example with measured $\psi$3 and $\psi$4 error near $\psi$5; 3D examples with fourth-order accuracy; high-contrast, close-interface, and heterogeneous cases included	(Zhou et al., 2023)
Brinkman and Stokes interface problems	Second-order accuracy in velocity and pressure and nearly second-order in velocity-gradient; fixed- and moving-interface tests; mesh-independent or mildly coefficient-dependent GMRES counts	(Zhou et al., 16 Apr 2026)
Elliptic interface problems on surfaces	Second-order accuracy; 8–20 GMRES iterations in representative tests; CPU time $\psi$6	(Yin et al., 22 Aug 2025)

These results show that KFBI is not restricted to a single operator class. The method has been validated on bidomain electrophysiology, incompressible flow, linear elasticity, coefficient-jump interface problems, and surface PDEs. The reported examples include smooth model geometries such as circles, spheres, tori, helicoids, saddles, elliptic paraboloids, and spheroids, but also more application-driven geometries such as a real human left ventricle and biomolecular-style high-contrast interface configurations (Gao et al., 2021, Zhou et al., 2023).

The numerical evidence also emphasizes a recurring trade-off. KFBI retains the favorable conditioning of second-kind boundary integral equations while shifting the dominant computational work to structured-grid solves. This is why the reported iteration counts are typically weakly dependent on mesh refinement, while the principal asymptotic cost comes from FFT or multigrid applications inside each operator evaluation.

6. Variants, distinctions, and current directions

A recurrent misconception is that “kernel-free” means the kernel disappears from both analysis and implementation. The KFBI literature does not support that reading. Green’s functions and layer potentials remain the analytical basis of the boundary integral reformulation; what is removed is the need to know or numerically integrate the explicit analytical kernel when solving the boundary integral equations (Yin et al., 22 Aug 2025). The method is therefore best understood as a hybridization of boundary integral formulations with regular-grid interface solvers.

A second distinction concerns neighboring nonsingular surface-integral methods. In field-only PEC scattering, three scalar Helmholtz surface integral equations are solved directly for the components of the electric field, the divergence-free condition is imposed through a Robin-type boundary condition involving the mean curvature, and a subtraction identity with auxiliary Helmholtz solutions $\psi$7 and $\psi$8 removes the $\psi$9 singularity. The resulting integrands remain finite, so standard Gaussian quadrature on quadratic six-node triangular elements is sufficient; the method is free of singular kernels and does not exhibit the low-frequency breakdown seen in current-based formulations (Sun et al., 2019). This is closely related in spirit but technically different from the PDE-solve-based KFBI paradigm.

Recent work also explores data-driven accelerations. A hybrid KFBI method with operator learning trains a network to approximate the solution operator that maps parameters, inhomogeneous terms, and boundary information to boundary density functions. The trained model can directly infer the boundary density function with satisfactory precision, eliminating iterative solution of the boundary integral equation for coarse-accuracy use; alternatively, its inference can be used as an initial value, retaining the inherent second-order accuracy of the KFBI method while reducing about $\phi_s$0 of the iterations in the traditional KFBI approach (Ling et al., 2024).

Current extension paths in the literature are explicit. The correction-function Brinkman framework is stated to generalize to three dimensions, Darcy–Brinkman–Stokes coupled systems via domain decomposition, time-dependent Navier–Stokes interface problems, and electrophysiology or biophysical multi-physics (Zhou et al., 16 Apr 2026). Surface KFBI extends the same operator-evaluation philosophy to PDEs on embedded manifolds (Yin et al., 22 Aug 2025). Taken together, these developments indicate that KFBI has evolved into a broad computational pattern: second-kind boundary integral formulations provide the operator structure, while Cartesian-grid or parameter-domain solvers provide the numerical realization without explicit singular-kernel quadrature.