Kansa Method for PDE Collocation

• Kansa Method is a strong-form, kernel-based collocation approach that approximates solutions of differential equations using translated kernels and direct operator evaluations at collocation points.
• It employs asymmetric collocation, least-squares, and oversampling strategies to achieve optimal H2 convergence and enhance numerical stability for various engineering applications.
• The method adapts to diverse problem settings—including diffusion, Hamiltonian waves, and fractional transport—by integrating advanced numerical solvers and structure-preserving constraints.

The Kansa method is the standard strong-form asymmetric kernel-based collocation method for differential equations. In its canonical form, one approximates an unknown field in a trial space generated by translated kernels and enforces the governing PDE and boundary conditions directly at collocation points, which yields an unsymmetric linear system because interior differential-operator rows and boundary-operator rows are different. The method is easy to implement and has therefore been widely used for engineering problems and partial differential equations; a substantial strand of later work reframes it through least-squares, oversampled asymmetric collocation, manifold-based formulations, and structure-preserving constrained solvers, thereby supplying convergence, stability, and application-specific analyses that the classical formulation long lacked (Cheung et al., 2018).

1. Canonical formulation and asymmetric collocation structure

For a second-order linear elliptic problem with Dirichlet data,

$$\begin{cases} L u = f & \text{in }\Omega,\ u = g & \text{on }\Gamma:=\partial\Omega, \end{cases}$$

the Kansa method seeks an approximation in the kernel trial space

$$U_Z=\operatorname{Span}\{\Phi(\cdot-z_j): z_j\in Z\},$$

and enforces

$$L u(z_i)=f(z_i),\quad z_i\in X,\qquad u(y_i)=g(y_i),\quad y_i\in Y.$$

Here $Z$ denotes trial centers, $X\subset\Omega$ interior collocation points, and $Y\subset\Gamma$ boundary collocation points. The defining feature is strong-form enforcement: the PDE operator and boundary conditions are collocated directly on the kernel expansion rather than being imposed in weak form (Cheung et al., 2018).

This asymmetric structure persists across later variants. On spheres and manifolds, one likewise chooses kernel centers $X$ or $Z$, forms a finite-dimensional kernel space such as $S_X(\Phi)$ or $\mathcal U_Z$, and evaluates the differential operator pointwise on that space. In time-dependent settings, the same idea is embedded inside operator splitting or time-stepping. For example, in the surface advection-diffusion setting, a Kansa-type method is applied to the diffusion subproblem along characteristic trajectories, while in Hamiltonian wave equations a least-squares Kansa discretization is coupled to Crank–Nicolson or Crank–Nicolson/Adams–Bashforth time stepping (Li et al., 26 Jan 2026).

The trial functions are radial or kernel-based and are usually global. Representative kernels explicitly mentioned in the literature summarized here include Whittle–Matérn–Sobolev kernels, Wendland kernels, the Matérn Sobolev kernel on manifolds, Hardy’s multiquadric basis, and Gaussian RBFs. The method is therefore both meshfree and kernel-dependent: geometry enters through scattered points and kernel evaluations rather than through a background mesh (Cheung et al., 2018).

2. Functional-analytic setting and $$U_Z=\operatorname{Span}\{\Phi(\cdot-z_j): z_j\in Z\},$$0-convergent least-squares variants

A major theoretical development places the Kansa method in a Sobolev approximation framework for general second-order linear elliptic operators

$$U_Z=\operatorname{Span}\{\Phi(\cdot-z_j): z_j\in Z\},$$1

with $$U_Z=\operatorname{Span}\{\Phi(\cdot-z_j): z_j\in Z\},$$2 strongly elliptic, coefficients in $$U_Z=\operatorname{Span}\{\Phi(\cdot-z_j): z_j\in Z\},$$3, solution $$U_Z=\operatorname{Span}\{\Phi(\cdot-z_j): z_j\in Z\},$$4, and a domain that is Lipschitz, satisfies an interior cone condition, and has piecewise $$U_Z=\operatorname{Span}\{\Phi(\cdot-z_j): z_j\in Z\},$$5 boundary. The kernel $$U_Z=\operatorname{Span}\{\Phi(\cdot-z_j): z_j\in Z\},$$6 is assumed to reproduce $$U_Z=\operatorname{Span}\{\Phi(\cdot-z_j): z_j\in Z\},$$7 through the Fourier decay condition

$$U_Z=\operatorname{Span}\{\Phi(\cdot-z_j): z_j\in Z\},$$8

so that its native space is norm-equivalent to $$U_Z=\operatorname{Span}\{\Phi(\cdot-z_j): z_j\in Z\},$$9 (Cheung et al., 2018).

Point-set geometry is encoded through fill distance, separation distance, and mesh ratio: $$L u(z_i)=f(z_i),\quad z_i\in X,\qquad u(y_i)=g(y_i),\quad y_i\in Y.$$0 Analogous quantities are used for $$L u(z_i)=f(z_i),\quad z_i\in X,\qquad u(y_i)=g(y_i),\quad y_i\in Y.$$1 and $$L u(z_i)=f(z_i),\quad z_i\in X,\qquad u(y_i)=g(y_i),\quad y_i\in Y.$$2. The analysis assumes asymptotic quasi-uniformity,

$$L u(z_i)=f(z_i),\quad z_i\in X,\qquad u(y_i)=g(y_i),\quad y_i\in Y.$$3

together with explicit denseness conditions relating the trial centers to interior and boundary sampling. These conditions are the basis for discrete stability in $$L u(z_i)=f(z_i),\quad z_i\in X,\qquad u(y_i)=g(y_i),\quad y_i\in Y.$$4 (Cheung et al., 2018).

Within this framework, the constrained least-squares formulation

$$L u(z_i)=f(z_i),\quad z_i\in X,\qquad u(y_i)=g(y_i),\quad y_i\in Y.$$5

enforces boundary data exactly and minimizes only the interior PDE residual. Under the stated smoothness and density assumptions, it satisfies

$$L u(z_i)=f(z_i),\quad z_i\in X,\qquad u(y_i)=g(y_i),\quad y_i\in Y.$$6

which is the optimal $$L u(z_i)=f(z_i),\quad z_i\in X,\qquad u(y_i)=g(y_i),\quad y_i\in Y.$$7 rate in the sense that the decay matches the underlying interpolation benchmark (Cheung et al., 2018).

Weighted least squares introduces boundary residuals as penalized terms: $$L u(z_i)=f(z_i),\quad z_i\in X,\qquad u(y_i)=g(y_i),\quad y_i\in Y.$$8 In the reported theory, $$L u(z_i)=f(z_i),\quad z_i\in X,\qquad u(y_i)=g(y_i),\quad y_i\in Y.$$9 maximizes the theoretical convergence rate, and for this choice the weighted formulation matches the constrained least-squares rate. The smaller trial space $Z$0, which does not include boundary points, remains convergent but carries extra boundary terms in the estimates; numerically it is described as computationally cheaper, with moderate weights such as $Z$1 or $Z$2 effective and stable (Cheung et al., 2018).

Three ingredients organize the convergence proof: sampling inequalities on $Z$3 and $Z$4, an inverse/Bernstein-type inequality on the trial space,

$Z$5

and a stability-plus-consistency argument comparing the least-squares minimizer to a kernel interpolant. This places the Kansa method, when reformulated as least-squares kernel collocation, on a rigorous approximation-theoretic footing (Cheung et al., 2018).

3. Oversampling, norming sets, and rectangular Kansa-like matrices

A distinct line of analysis extends the classical square Kansa system into an oversampled asymmetric collocation method on spheres and compact manifolds. In the classical form for an elliptic PDE $Z$6, one seeks

$Z$7

such that

$Z$8

The associated matrix,

$Z$9

is generally not symmetric positive definite, and in the square case $X\subset\Omega$0 it can be badly conditioned or singular for certain point sets (Hangelbroek et al., 20 Jul 2025).

The oversampled formulation replaces the usual test set by a larger set $X\subset\Omega$1, often chosen to be a norming set for a differentiated trial space, and it permits alternative bases $X\subset\Omega$2, especially Lagrange or local Lagrange bases. The resulting “Kansa-like” matrix is

$X\subset\Omega$3

The method thereby becomes overdetermined rather than square (Hangelbroek et al., 20 Jul 2025).

The key analytical notion is the norming inequality. For a finite-dimensional space $X\subset\Omega$4,

$X\subset\Omega$5

and in the $X\subset\Omega$6 case,

$X\subset\Omega$7

When $X\subset\Omega$8 is norming for the image of the trial space under $X\subset\Omega$9, the matrix acquires a lower bound of the form

$Y\subset\Gamma$0

and, with a Riesz basis, essentially

$Y\subset\Gamma$1

This ties stability to oversampling and basis quality rather than to square collocation alone (Hangelbroek et al., 20 Jul 2025).

Two solver paradigms are then analyzed. The first is discrete least squares,

$Y\subset\Gamma$2

with a resulting approximation $Y\subset\Gamma$3. For thin-plate splines and sufficiently smooth solutions $Y\subset\Gamma$4, the reported estimate is

$Y\subset\Gamma$5

The second uses a strong RRQR factorization to thin the rectangular system to a square one; the corresponding error estimate is

$Y\subset\Gamma$6

This analysis reinterprets the Kansa method as an oversampled asymmetric kernel method with rectangular matrices, norming-set stability, and two practical solvers: least squares and thinning (Hangelbroek et al., 20 Jul 2025).

4. Surface, wave, and other time-dependent formulations

On smooth compact manifolds, the Kansa method appears as the diffusion discretization inside a trajectory-based RBF collocation framework for the surface advection-diffusion equation

$Y\subset\Gamma$7

The method first treats advection by characteristics through surface-restricted backtracking,

$Y\subset\Gamma$8

and then applies a Kansa-type global RBF expansion

$Y\subset\Gamma$9

to the diffusion step, with trial centers $X$0 and collocation points $X$1 (Li et al., 26 Jan 2026).

After Crank–Nicolson time discretization, collocation produces the unsymmetric algebraic system

$X$2

Here the Laplace–Beltrami operator replaces the Euclidean Laplacian, and the right-hand side is evaluated at backtracked points. The paper further proves a chain of equivalences: the original surface PDE is equivalent to an embedded PDE in a narrow band, that embedded PDE is equivalent to an operator-split characteristic system, and the restriction of that system to the manifold recovers the surface dynamics with no operator splitting error (Li et al., 26 Jan 2026).

For second-order time-dependent Hamiltonian wave equations, the Kansa method is used in a least-squares strong-form discretization with oversampling,

$X$3

where $X$4. At each time level, the unconstrained least-squares Kansa formulation is

$X$5

The principal innovation is to embed energy conservation as a nonlinear quadratic constraint: $X$6 where the constraint enforces equality of the discrete energy and the prescribed initial energy $X$7. The solver uses GSVD, Lagrange multipliers, and Newton iteration with successive linearization. In the reported experiments, the resulting EC-LS Kansa methods preserve energy to near machine precision, CN and CNAB are second-order in time, CNAB is generally faster, and CPU-time reductions relative to a traditional secant method are about $X$8 to more than $X$9 (Li et al., 6 Jul 2025).

These developments broaden the meaning of “Kansa method.” In the manifold setting it functions as a surface diffusion engine inside a characteristic framework; in Hamiltonian waves it becomes a structure-preserving constrained solver whose spatial core remains asymmetric kernel collocation. The common element is still direct operator evaluation on a kernel trial space, but the surrounding algorithmic architecture has become substantially richer (Li et al., 26 Jan 2026).

5. Fractional transport and biological ODE systems

For two-dimensional spatiotemporal fractional diffusion and advection-dispersion equations, the Kansa method is used to discretize the spatial fractional derivative while time is handled analytically. The spatial approximation is the global MQ-RBF expansion

$Z$0

with $Z$1 interior nodes and $Z$2 boundary nodes. PDE collocation at interior nodes and boundary enforcement at boundary nodes lead to the semi-discrete fractional ODE system

$Z$3

After shifting the source term, one obtains

$Z$4

whose solution is expressed through the Mittag-Leffler function,

$Z$5

The fractional derivatives of the basis are evaluated by Gauss-Jacobi quadrature, which is the paper’s main mechanism for handling weakly singular fractional operators (Sun et al., 2016).

The reported numerical behavior in that setting includes super-linear convergence in one dimension, low MAEs and relative errors in two-dimensional tests, and nearly constant relative error over long final times because temporal evolution is not advanced by step-by-step marching. The method is presented as especially attractive for irregular domains, circular domains, unit disks, and anisotropic transport, and as a fast deterministic alternative to previous particle Monte-Carlo approaches for solute plumes and related anomalous diffusion problems (Sun et al., 2016).

In nonlinear biological ODE systems, the Kansa method becomes a meshfree RBF collocation procedure on the time interval. Each state variable is approximated by an RBF expansion such as

$Z$6

and the governing equations are enforced at equidistant collocation points, with initial conditions appended as algebraic constraints. This produces nonlinear algebraic systems, solved in the cited work by Newton–Raphson. The method is applied to an HIV model and an Influenza SIRC model (Salehi et al., 2017).

A central issue in that application is the RBF shape parameter $Z$7. Because no universal analytical rule is reported for selecting an optimal value, the paper formulates shape selection as a continuous genetic algorithm problem. Fitness is defined by

$Z$8

where $Z$9 is ASN2R for HIV or ARE for Influenza. Parent selection uses the Roulette Wheel Technique, and crossover uses the pseudo-combination formula

$S_X(\Phi)$0

The reported experiments show close agreement with Runge–Kutta or Runge–Kutta–Fehlberg references, successful narrowing of the shape-parameter search interval, and residuals decreasing to very small values as the number of collocation points increases (Salehi et al., 2017).

6. Numerical behavior, conditioning, and recurring analytical themes

Across these formulations, several recurrent properties define the modern understanding of the Kansa method.

Formulation	Defining algebraic form	Reported property
Strong-form Kansa	Direct pointwise PDE and boundary collocation	Unsymmetric linear system
CLS Kansa	$S_X(\Phi)$1 with exact boundary constraint	$S_X(\Phi)$2-convergent, optimal $S_X(\Phi)$3
WLS Kansa	Interior residual plus weighted boundary residual	$S_X(\Phi)$4 maximizes theoretical rate
Oversampled Kansa-like	Rectangular $S_X(\Phi)$5, $S_X(\Phi)$6	Stability lower bounds via norming sets
TBRBF diffusion step	Kansa collocation after characteristic backtracking	No operator splitting error
EC-LS Kansa	Least squares with nonlinear quadratic energy constraint	Near machine precision energy conservation

The first recurring issue is that the Kansa matrix is generally unsymmetric and need not be symmetric positive definite. In the square case $S_X(\Phi)$7, later work explicitly notes that the matrix can be badly conditioned or even singular for certain point sets. This corrects a common simplification that treats the method as merely a meshfree interpolation scheme; its stability depends on operator action, basis choice, oversampling, and point-set geometry, not only on kernel interpolation properties (Hangelbroek et al., 20 Jul 2025).

The second recurring issue is that least-squares and oversampling are not peripheral modifications but central analytical devices. The rigorous $S_X(\Phi)$8-convergence theory for elliptic Dirichlet problems is formulated for constrained or weighted least-squares variants, while the manifold theory based on norming sets uses rectangular matrices and either least squares or RRQR-based thinning. This suggests that the most analyzable forms of the Kansa method are often overdetermined rather than square (Cheung et al., 2018).

The third issue is conditioning versus accuracy. In the elliptic least-squares study, large weights or very high smoothness can increase sensitivity to ill-conditioning, and stable Gaussian evaluation via the RBF-QR technique can dramatically improve performance. The same study states that the condition number alone does not predict accuracy or convergence rate, especially for Gaussian and multiquadric kernels. A plausible implication is that implementation quality, evaluation stability, and sampling geometry must be treated as coequal design variables alongside kernel smoothness (Cheung et al., 2018).

The fourth issue is that operator splitting and structure preservation can coexist with Kansa discretization when the surrounding framework is designed carefully. In the surface advection-diffusion setting, the characteristic formulation is proved equivalent to the original PDE with no operator splitting error. In Hamiltonian wave equations, the spatial Kansa discretization is embedded in an energy-conserving constrained solver that preserves the discrete Hamiltonian to near machine precision. These results expand the method from a collocation technique into a component of broader geometric and invariant-preserving algorithms (Li et al., 26 Jan 2026).

Taken together, the cited literature portrays the Kansa method not as a single fixed algorithm but as a family of asymmetric kernel-collocation strategies. Its core remains the same—global kernel trial spaces and direct operator collocation—but its modern forms include optimal $S_X(\Phi)$9-convergent least-squares elliptic schemes, oversampled rectangular formulations on spheres and manifolds, surface methods based on characteristics and Laplace–Beltrami operators, energy-conserving solvers for Hamiltonian waves, semi-discrete fractional transport solvers with analytic time propagation, and nonlinear ODE collocation procedures whose practical success depends strongly on shape-parameter selection (Li et al., 6 Jul 2025).