KernelInterpolation.jl API

KernelInterpolation

KernelInterpolation.jl is a Julia package that implements methods for multivariate interpolation in arbitrary dimension based on symmetric (conditionally) positive-definite kernels with a focus on radial basis functions. It can be used for classical interpolation of scattered data, as well as for generalized (Hermite-Birkhoff) interpolation by using a meshfree collocation approach. This can be used to solve partial differential equations both stationary ones and time-dependent ones by using some time integration method from OrdinaryDiffEq.jl.

See also: KernelInterpolation.jl

AbstractKernel

An abstract supertype of kernels.

dim(kernel)

Return the dimension of a kernel, i.e. the size of the input vector.

get_name(kernel::AbstractKernel)

Returns the canonical, human-readable name for the given system of equations.

GaussKernel{Dim}(; shape_parameter = 1.0)

Gaussian kernel function with

\[ \phi(r) = \exp(-(\varepsilon r)^2),\]

where $\varepsilon$ is the shape parameter. The Gaussian kernel is always positive definite. See Wendland (2004), p. 74.

See also RadialSymmetricKernel.

• Holger Wendland (2004) Scattered Data Approximation Cambridge University Press DOI: 10.1017/CBO9780511617539

InverseMultiquadricKernel{Dim}(beta = 0.5; shape_parameter = 1.0)

Inverse multiquadric kernel function with

\[ \phi(r) = (1 + (\varepsilon r)^2)^{-\beta},\]

where $\varepsilon$ is the shape parameter. The inverse multiquadric kernel is always positive definite. See Wendland (2004), p. 76 and p. 95.

See also RadialSymmetricKernel.

• Holger Wendland (2004) Scattered Data Approximation Cambridge University Press DOI: 10.1017/CBO9780511617539

Matern12Kernel{Dim}(; shape_parameter = 1.0)

Matern kernel with $\nu = 1/2$, i.e.,

\[ \phi(r) = \exp(-\varepsilon r),\]

where $\varepsilon$ is the shape parameter. The Matern kernel is positive definite.

See Wikipedia and Fasshauer (2007), p. 41.

See also MaternKernel, RadialSymmetricKernel.

• Gregory Fasshauer (2007) Meshfree Approximation Methods with MATLAB World Scientific DOI: 10.1142/6437

Matern32Kernel{Dim}(; shape_parameter = 1.0)

Matern kernel with $\nu = 3/2$, i.e.,

\[ \phi(r) = (1 + \sqrt{3}\varepsilon r)\exp(-\sqrt{3}\varepsilon r),\]

where $\varepsilon$ is the shape parameter. The Matern kernel is positive definite.

See Wikipedia and Fasshauer (2007), p. 41.

See also MaternKernel, RadialSymmetricKernel.

• Gregory Fasshauer (2007) Meshfree Approximation Methods with MATLAB World Scientific DOI: 10.1142/6437

Matern52Kernel{Dim}(; shape_parameter = 1.0)

Matern kernel with $\nu = 5/2$, i.e.,

\[ \phi(r) = (1 + \sqrt{5}\varepsilon r + 5\cdot(\varepsilon r)^2/3)\exp(-\sqrt{5}\varepsilon r),\]

where $\varepsilon$ is the shape parameter. The Matern kernel is positive definite.

See Wikipedia and Fasshauer (2007), p. 41.

See also MaternKernel, RadialSymmetricKernel.

• Gregory Fasshauer (2007) Meshfree Approximation Methods with MATLAB World Scientific DOI: 10.1142/6437

Matern72Kernel{Dim}(; shape_parameter = 1.0)

Matern kernel with $\nu = 7/2$, i.e.,

\[ \phi(r) = (1 + \sqrt{7}\varepsilon r + 12\cdot(\varepsilon r)^2/5 + 7\cdot(\varepsilon r)^3/15)\exp(-\sqrt{7}\varepsilon r),\]

where $\varepsilon$ is the shape parameter. The Matern kernel is positive definite.

See Wikipedia and Fasshauer (2007), p. 41.

See also MaternKernel, RadialSymmetricKernel.

• Gregory Fasshauer (2007) Meshfree Approximation Methods with MATLAB World Scientific DOI: 10.1142/6437

MaternKernel{Dim}(nu = 1.5; shape_parameter = 1.0)

Matern kernel with

\[ \phi_\nu(r) = \frac{2^{1-\nu}}{\Gamma(\nu)}\big(\sqrt{2\nu}\varepsilon r\big)^\nu K_\nu\big(\sqrt{2\nu}\varepsilon r\big),\]

where $\varepsilon$ is the shape parameter. The Matern kernel is positive definite.

See Wikipedia and Fasshauer (2007), p. 41.

See also RadialSymmetricKernel.

• Gregory Fasshauer (2007) Meshfree Approximation Methods with MATLAB World Scientific DOI: 10.1142/6437

MultiquadricKernel{Dim}(beta = 0.5; shape_parameter = 1.0)

Multiquadric kernel function with

\[ \phi(r) = (1 + (\varepsilon r)^2)^\beta,\]

where $\varepsilon$ is the shape parameter. The multiquadric kernel is conditionally positive definite of order $m = \lceil\beta \rceil$. See Wendland (2004), p. 109.

See also RadialSymmetricKernel.

• Holger Wendland (2004) Scattered Data Approximation Cambridge University Press DOI: 10.1017/CBO9780511617539

PolyharmonicSplineKernel{Dim}(k)

Polyharmonic spline kernel function with

\[ \phi_k(r) = \begin{cases} r^k, &\text{ if } k \text{ odd}\\ r^k\log(r), &\text{ if } k \text{ even} \end{cases}.\]

The polyharmonic spline is conditionally positive definite of order $m = \lceil k/2\rceil$ for odd k and order $m = k/2 + 1$ for even k. See Wendland (2004), pp. 111–112.

See also RadialSymmetricKernel.

• Holger Wendland (2004) Scattered Data Approximation Cambridge University Press DOI: 10.1017/CBO9780511617539

RadialCharacteristicKernel{Dim}(beta = 2.0; shape_parameter = 1.0)

Radial characteristic function (or also called truncated power or Askey) kernel function with

\[ \phi(r) = (1 - \varepsilon r)^\beta_+,\]

where $\varepsilon$ is the shape parameter. The radial characteristic function is positive definite if $\beta\ge (d + 1)/2$. It is compactly supported. See Wendland (2004), p. 80, Iske (2018), p. 281.

See also RadialSymmetricKernel.

• Holger Wendland (2004) Scattered Data Approximation Cambridge University Press DOI: 10.1017/CBO9780511617539
• Armin Iske (2018) Approximation Theory and Algorithms for Data Analysis Texts in Applied Mathematics (Springer) DOI: 10.1007/978-3-030-05228-7

RadialSymmetricKernel

An abstract supertype of radial symmetric kernels. Radial symmetric kernels are generated by an even and continuous function $\Phi: \mathbb{R}^d\to\mathbb{R}$, which is radial-symmetric meaning that there exists a $\phi:[0,\infty)\to\mathbb{R}$ such that

\[ \Phi(x) = \phi(\Vert x\Vert).\]

The kernel is then defined by

\[ K(x, y) = \Phi(x - y).\]

A RadialSymmetricKernel can be evaluated at two points x and y by calling kernel(x, y) or at a single point x by calling kernel(x), which implicitly sets y to zero.

RieszKernel{Dim}(beta; shape_parameter = 1.0)

Riesz kernel with

\[ \phi(r) = -(\varepsilon r)^\beta,\]

where $\varepsilon$ is the shape parameter and $\beta\in (0,2)$. The Riesz kernel is conditionally positive definite of order 1. See Hertrich et al. (2023).

See also RadialSymmetricKernel.

• Johannes Hertrich, Christian Wald, Fabian Altekrüger, Paul Hagemann (2023) Generative Sliced MMD Flows with Riesz Kernels arXiv: 2305.11463

ThinPlateSplineKernel{Dim}()

Thin plate spline kernel function with

\[ \phi(r) = r^2\log(r),\]

i.e., PolyharmonicSplineKernel with $k = 2$. The thin plate spline is conditionally positive definite of order $m = 2$. See Wendland (2004), p. 112.

See also RadialSymmetricKernel.

• Holger Wendland (2004) Scattered Data Approximation Cambridge University Press DOI: 10.1017/CBO9780511617539

WendlandKernel{Dim}(k; shape_parameter = 1.0, d = Dim)

Wendland kernel with

\[ \phi_{d,k}(r) = \begin{cases} p_{d,k}(\varepsilon r), \text{ if } 0\le \varepsilon r\le 1\\ 0, \text{ if } \varepsilon r > 1 \end{cases},\]

where $\varepsilon$ is the shape parameter and $p$ is a polynomial with minimal degree. The Wendland kernel is positive definite for $d\le$ Dim and compactly supported. See Wendland (2004), p. 129 or Fasshauer (2007), pp. 87–88.

See also RadialSymmetricKernel.

• Holger Wendland (2004) Scattered Data Approximation Cambridge University Press DOI: 10.1017/CBO9780511617539
• Gregory Fasshauer (2007) Meshfree Approximation Methods with MATLAB World Scientific DOI: 10.1142/6437

WuKernel{Dim}(l, k; shape_parameter = 1.0)

Wu kernel with

\[ \phi_{l,k}(r) = \begin{cases} p_{l,k}(\varepsilon r), \text{ if } 0\le \varepsilon r\le 1\\ 0, \text{ if } \varepsilon r > 1 \end{cases},\]

where $\varepsilon$ is the shape parameter, $k\le l$, and $p$ is a polynomial of degree $4l - 2k + 1$. The Wu kernel is positive definite for Dim $\le 2k + 1$ and compactly supported. See Fasshauer (2007), pp. 88–90 and Wu (1995).

See also RadialSymmetricKernel.

• Gregory Fasshauer (2007) Meshfree Approximation Methods with MATLAB World Scientific DOI: 10.1142/6437
• Zongmin Wu (1995) Compactly supported positive definite radial functions Advances in Computational Mathematics DOI: 10.1007/BF03177517

Phi(kernel, x)

For a RadialSymmetricKernel kernel return value of the multivariate function $\Phi$ defined by $\Phi(x) = \phi(\Vert x\Vert)$.

See also RadialSymmetricKernel, phi.

order(kernel)

Return order of kernel.

phi(kernel, r)

For a RadialSymmetricKernel kernel return value of the univariate function $\phi$ defining the kernel.

See also RadialSymmetricKernel, Phi.

ProductKernel{Dim}(kernels)

Given a vector of kernels, construct a new kernel that multiplies the results of the component kernels, i.e., the new kernel $K$ is given by

\[ K(x, y) = \prod_{i = 1}^n K_i(x, y),\]

where $K_i$ are the component kernels and $n$ the number of kernels. Note that all component kernels need to have the same dim. A ProductKernel can also be constructed using the * operator.

SumKernel{Dim}(kernels)

Given a vector of kernels, construct a new kernel that sums the results of the component kernels, i.e., the new kernel $K$ is given by

\[ K(x, y) = \sum_{i = 1}^n K_i(x, y),\]

where $K_i$ are the component kernels and $n$ the number of kernels. Note that all component kernels need to have the same dim. A SumKernel can also be constructed using the + operator.

TransformationKernel{Dim}(kernel, transformation)

Given a base kernel and a bijective transformation function, construct a new kernel that applies the transformation to both arguments $x$ and $y$, i.e., the new kernel $K_T$ is given by

\[ K_T(x, y) = K(Tx, Ty),\]

where $K$ is the base kernel and $T$ the transformation, i.e. if $K$ is a kernel of dimension $d$, $T$ is a function from dimension Dim to $d$, where Dim is the dimension of the new kernel.

NodeSet(nodes)

Set of interpolation nodes.

distance_matrix(nodeset1::NodeSet, nodeset2::NodeSet)

Compute the distance matrix between two NodeSets, which is a matrix $D$ with $D_{ij} = \|x_i - \xi_j\|$ for all $i$ and $j$, where $x_i$ are the nodes in nodeset1 and $\xi_j$ are the nodes on nodeset2.

empty_nodeset(Dim, RealT = Float64)

Create an empty NodeSet.

fill_distance(nodeset, reference)

Approximate the fill distance of a NodeSet $X = \{x_1,\ldots, x_n\}$ with respect to a domain $\Omega$ represented by a reference NodeSet $\Xi$:

\[ h_{X,\Omega} \approx \max_{\xi_i \in \Xi}\,\min_{x_j \in X} \|\xi_i - x_j\|.\]

The result is an approximation on the true fill distance

\[ h_{X,\Omega} = \sup_{x \in \Omega}\,\min_{x_j \in X} \|x - x_j\|;\]

accuracy improves with the density of reference.

See also separation_distance, distance_matrix.

homogeneous_hypercube(n, x_min = ntuple(_ -> 0.0, dim), x_max = ntuple(_ -> 1.0, dim); [dim])

If n is integer, create a NodeSet with n homogeneously distributed nodes in every dimension each of dimension dim inside a hypercube defined by the bounds x_min and x_max. If n is a Tuple of length dim, then use as many nodes in each dimension as described by n. The resulting NodeSet will have $n^{\textrm{dim}}$ respectively $\prod_{j = 1}^{\textrm{dim}}n_j$ points. If the bounds are given as single values, they are applied for each dimension. If they are Tuples of size dim, the hypercube has the according bounds. If dim is not given explicitly, it is inferred by the lengths of n, x_min and x_max if possible.

homogeneous_hypercube_boundary(n, x_min = ntuple(_ -> 0.0, dim), x_max = ntuple(_ -> 1.0, dim); [dim])

If n is integer, create a NodeSet with n homogeneously distributed nodes in every dimension each of dimension dim on the boundary of a hypercube defined by the bounds x_min and x_max. If n is a Tuple of length dim, then use as many nodes in each dimension as described by n. If the bounds are given as single values, they are applied for each dimension. If they are Tuples of size dim, the hypercube has the according bounds. If dim is not given explicitly, it is inferred by the lengths of n, x_min and x_max if possible.

random_hypercube([rng], n, x_min = ntuple(_ -> 0.0, dim), x_max = ntuple(_ -> 1.0, dim); [dim])

Create a NodeSet with n random nodes each of dimension dim inside a hypercube defined by the bounds x_min and x_max. If the bounds are given as single values, they are applied for each dimension. If they are Tuples of size dim the hypercube has the according bounds. If dim is not given explicitly, it is inferred by the lengths of x_min and x_max if possible. Optionally, pass a random number generator rng.

random_hypercube_boundary([rng], n, x_min = ntuple(_ -> 0.0, dim), x_max = ntuple(_ -> 1.0, dim); [dim])

Create a NodeSet with n random nodes each of dimension dim on the boundary of a hypercube defined by the bounds x_min and x_max. If the bounds are given as single values, they are applied for each dimension. If they are Tuples of size dim the hypercube has the according bounds. If dim is not given explicitly, it is inferred by the lengths of x_min and x_max if possible. Optionally, pass a random number generator rng.

random_hypersphere([rng], n, r = 1.0, center = Tuple(zeros(dim)); [dim])

Create a NodeSet with n random nodes each of dimension dim inside a hypersphere with radius r around the center center, which is given as a tuple. If dim is not given explicitly, it is inferred by the length of center if possible. Optionally, pass a random number generator rng.

random_hypersphere_boundary([rng], n, r = 1.0, center = Tuple(zeros(dim)); [dim])

Create a NodeSet with n random nodes each of dimension dim at the boundary of a hypersphere with radius r around the center center, which is given as a tuple. If dim is not given explicitly, it is inferred by the length of center if possible. Optionally, pass a random number generator rng.

separation_distance(nodeset::NodeSet)

Return the separation distance of a NodeSet $X = \{x_1,\ldots, x_n\}$ defined by

\[ q_X = \frac{1}{2}\min_{x_i\neq x_j}\|x_i - x_j\|.\]

AbstractBasis

Abstract type for a basis of a kernel function space. Every basis represents a set of functions, which can be obtained by indexing the basis object. Every basis object holds a kernel function and a NodeSet of centers and potentially more fields depending on the concrete basis type.

AbstractRBFFDLocalBasis

Abstract type selecting the local basis representation used to compute RBF-FD weights.

See also RBFFDStandardBasis and RBFFDLagrangeBasis.

LagrangeBasis(centers, kernel, m = order(kernel))

The Lagrange (or cardinal) basis with respect to a kernel and a NodeSet of centers. This basis already includes polynomial augmentation of order m defaulting to order(kernel). The basis functions are given such that

\[ b_j(x_i) = \delta_{ij},\]

which means that the kernel_matrix of this basis is the identity matrix making it suitable if multiple interpolations with the same centers of the basis and the same kernel, but with different right-hand sides or nodesets are performed. Since the basis already includes polynomials no additional polynomial augmentation is needed for interpolation with this basis.

RBFFDBasis(nodeset, kernel, stencil_selection;
           m = order(kernel), local_basis = RBFFDLagrangeBasis())

Configuration object for RBF-FD discretizations.

RBFFDBasis stores the global node set, kernel, and stencil strategy used to build local RBF-FD weights. All stencils are precomputed at construction time. What else is precomputed is selected by local_basis:

• RBFFDLagrangeBasis (default): the Lagrange cardinal functions on each stencil are precomputed and stored. Weight computation, matrix assembly, and interpolant evaluation evaluate (operators applied to) these cardinal functions.
• RBFFDStandardBasis: the factorization of each stencil's augmented kernel/ polynomial system is precomputed and stored. The same operations are obtained by solving the local system with the appropriate right-hand side.

Both choices produce the same weights mathematically; they differ only in the numerical route and in what basis[i, k] returns (see getindex).

Use basis[i, k] to access the k-th local basis function on the stencil around center index i.

The polynomial augmentation is local to each stencil; its order is given by local_order (and the monomials by polynomial_basis). No global polynomials are added, so order of an RBFFDBasis is always 0, consistent with LagrangeBasis.

RBFFDLagrangeBasis()

Compute local RBF-FD weights from local cardinal (Lagrange) basis functions on each stencil, i.e., w_j = 𝓛 ℓ_j(x_i). The cardinal functions are built once at construction and reused for weight computation, matrix assembly, and interpolant evaluation.

See also AbstractRBFFDLocalBasis and RBFFDBasis.

RBFFDStandardBasis()

Compute local RBF-FD weights by solving the local kernel (or kernel+polynomial) system. The per-stencil augmented system matrix is factorized once at construction and reused for weight computation, matrix assembly, and interpolant evaluation.

See also AbstractRBFFDLocalBasis and RBFFDBasis.

StandardBasis(centers, kernel)

The standard basis for a function space defined by a kernel and a NodeSet of centers. The basis functions are given by

\[ b_j(x) = K(x, x_j)\]

where K is the kernel and x_j are the nodes in centers.

getindex(basis::RBFFDBasis, i::Integer, k::Integer)

Return the k-th local stencil basis function associated with center index i. The returned object is callable.

In both cases k ranges over the stencil size. For RBFFDLagrangeBasis it is the k-th Lagrange cardinal function (with polynomial augmentation baked in, like the global LagrangeBasis); for RBFFDStandardBasis it is the k-th kernel translate x -> kernel(x, x_{i_k}), matching the global StandardBasis. The polynomial part is not exposed through indexing (see polynomial_basis); only for the Lagrange basis do the returned functions carry the cardinal (nodal-value) semantics.

centers(basis)

Return the centers from a basis object.

interpolation_kernel(basis)

Return the kernel from a basis.

local_order(basis::RBFFDBasis)

Return the order m of the polynomial basis applied separately on each RBF-FD stencil, of an RBFFDBasis, i.e., the polynomial degree plus 1 (0 if no separate polynomials are used). This is distinct from order, which is the (always zero) global augmentation.

Like the global StandardBasis/LagrangeBasis distinction, this depends on the local-basis policy: for RBFFDStandardBasis the polynomials are an explicit part of the local system, so local_order is their order m; for RBFFDLagrangeBasis they are baked into the cardinal functions, so there is no separate polynomial basis and local_order is 0.

order(basis)

Return the order $m$ of the polynomial, which is needed by this basis for the interpolation, i.e., the polynomial degree plus 1. If $m = 0$, no polynomial is added.

polynomial_basis(basis::RBFFDBasis)

Return the (stencil-independent) polynomial basis functions used for the local polynomial augmentation of the RBF-FD stencils. Empty if local_order(basis) == 0.

Interpolation

Interpolation object that can be evaluated at a node and represents a kernel interpolation of the form

\[ s(x) = \sum_{j = 1}^N c_jb_j(x) + \sum_{k = 1}^Q d_kp_k(x),\]

where $b_j$ are the basis functions and $p_k$ is a basis of the Q-dimensional space of multivariate polynomials of order order. The additional conditions

\[ \sum_{j = 1}^N c_jp_k(x_j) = 0, \quad k = 1,\ldots, Q\]

are enforced.

See also interpolate.

MultiscaleInterpolation

Container for a multiscale interpolation composed of several single-scale Interpolation objects. Evaluation is the sum of the evaluations of the individual scales.

See also multiscale_interpolate.

RBFFDInterpolation

Type alias for Interpolation{<:RBFFDBasis}.

The coefficient vector c holds nodal values $u(x_i)$ for all nodes $x_1,\ldots,x_N$ (inner nodes first, then boundary nodes), ordered consistently with merge(nodeset_inner, nodeset_boundary).

Evaluation at a point x uses the local stencil $S(j)$ of the nearest center $x_j$:

\[u_h(x) = \sum_{k \in S(j)} c_k \, w_k(x;\, S(j)),\]

where $w_k(x)$ are the local cardinal weights mapping the nodal values to the value of the local interpolant at x (i.e. local_weights(basis, j, x, Identity())). This formula holds regardless of the weight computation algorithm: for RBFFDLagrangeBasis the $w_k(x) = \ell_k(x)$ are the precomputed cardinal functions, while for RBFFDStandardBasis they are obtained by solving the cached local system with the evaluation right-hand side. Both give the same $w_k$ mathematically.

TemporalInterpolation(ode_sol::ODESolution)

Temporal interpolation of an ODE solution. The result can be evaluated at a time t and a spatial point x. Evaluating the interpolation at a time t returns an Interpolation object that can be evaluated at a spatial point x.

basis(itp)

Return the basis from an Interpolation object.

centers(itp)

Return the centers from the basis of an Interpolation object.

coefficients(itp::Interpolation)

Obtain all the coefficients of the linear combination for the Interpolation, i.e., both the coefficients for the kernel part and for the polynomial part.

See also kernel_coefficients and polynomial_coefficients.

dim(itp::Interpolation)

Return the dimension of the input variables of the Interpolation.

interpolate(basis, values, nodeset = centers(basis); m = order(basis),
            regularization = NoRegularization(), factorization_method = nothing,
            linsolve = nothing)
interpolate(centers, [nodeset,] values, kernel = GaussKernel{dim(nodeset)}();
            m = order(kernel), regularization = NoRegularization(),
            factorization_method = nothing, linsolve = nothing)

Interpolate the values evaluated at the nodes in the nodeset to a function using the kernel kernel and polynomials up to a order m (i.e. degree - 1), i.e., determine the coefficients $c_j$ and $d_k$ in the expansion

\[ s(x) = \sum_{j = 1}^N c_jb_j(x) + \sum_{k = 1}^Q d_kp_k(x),\]

where $b_j$ are the basis functions in the basis and $s(x)$ the interpolant $s(x_j) = f(x_j)$, where $f(x_j)$ are given by values, $x_j$ are the nodes in the nodeset, and $p_k$ is a basis of the $Q$-dimensional space of multivariate polynomials with maximum degree of m - 1. If m = 0, no polynomial is added. The additional conditions

\[ \sum_{j = 1}^N c_jp_k(x_j) = 0, \quad k = 1,\ldots, Q = \begin{pmatrix}m - 1 + d\\d\end{pmatrix}\]

are enforced. Returns an Interpolation object.

If nodeset is provided, the interpolant is a least squares approximation with a different set of nodes as the centers used for the basis. Otherwise, nodeset is set to centers(basis) or centers.

A regularization can be applied to the kernel matrix using the regularization argument, cf. regularize!. In addition, the factorization_method can be specified to determine how the system matrix is factorized. By default, the system matrix is just wrapped as a Symmetric matrix for interpolation and no factorization is applied for a least squares solution, but you can, e.g., also explicitly use cholesky, lu, or qr factorization. If linsolve is provided, the linear system is solved with LinearSolve.jl and any LinearSolve.jl algorithm can be passed there. If linsolve is not provided, the linear system is solved with the backslash operator, which will automatically use the factorization if factorization_method is provided.

interpolation_kernel(itp)

Return the kernel from an Interpolation object.

kernel_coefficients(itp::Interpolation)

Obtain the coefficients of the kernel part of the linear combination for the Interpolation.

See also coefficients and polynomial_coefficients.

kernel_inner_product(itp1, itp2)

Inner product of the native space for two Interpolation objects itp1 and itp2 with the same kernel. The inner product is defined as

\[ \langle f, g\rangle_K = \sum_{i = 1}^N\sum_{j = 1}^Mc_i^fc_j^gK(x_i, \xi_j)\]

for the interpolants $f(x) = \sum_{i = 1}^Nc_i^fK(x, x_i)$ and $g(x) = \sum_{j = 1}^Mc_j^gK(x, \xi_j)$.

See also kernel_norm.

kernel_norm(itp)

Norm of the native space defined by the kernel of the Interpolation itp. The norm is defined as

\[ \|f\|_K^2 = \sum_{i,j=1}^Nc_ic_jK(x_i, x_j)\]

for the interpolant $f(x) = \sum_{j = 1}^nc_jK(x, x_j)$.

See also kernel_inner_product.

multiscale_interpolate(nodesets, valuesets, kernels; kwargs...)

Construct a multiscale interpolation by fitting successive interpolants with the given kernels to the residual. The nodesets and valuesets must be provided as vectors of the same length as kernels, allowing grids and data to grow between scales. Each scale is constructed by calling interpolate on the corresponding nodeset and residual values.

• Armin Iske (2018) Multiresolution Methods in Scattered Data Modelling Lecture Notes in Computational Science and Engineering (Springer) DOI: 10.1007/978-3-642-18754-4

nodeset(itp)

Return the node set from an Interpolation object.

order(itp)

Return the order $m$ of the polynomial used for the Interpolation, i.e., the polynomial degree plus 1. If $m = 0$, no polynomial is added.

polynomial_basis(itp::Interpolation)

Return a vector of the polynomial basis functions used for the Interpolation.

See also polyvars.

polynomial_coefficients(itp::Interpolation)

Obtain the coefficients of the polynomial part of the linear combination for the Interpolation.

See also coefficients and kernel_coefficients.

polyvars(itp::Interpolation)

Return a vector of the polynomial variables for the Interpolation.

See also polynomial_basis.

system_matrix(itp::Interpolation)

Return the system matrix, i.e., the matrix $A$ in the linear system

\[ Ac = f,\]

where $c$ are the coefficients of the Interpolation and $f$ the vector of known values. The exact form of $A$ differs depending on which method is used.

AbstractRegularization

An abstract supertype of regularizations. A regularization implements a function regularize! that takes a matrix and returns a regularized version of it.

L2Regularization(regularization_parameter::Real)

A regularization that adds a multiple of the identity matrix to the input matrix.

NoRegularization()

A regularization that does nothing.

regularize!(A, reg::AbstractRegularization)

Apply the regularization reg to the matrix A in place.

EllipticOperator(A, b, c)

Linear second-order elliptic operator with matrix $A(x)\in\mathbb{R}^{d\times d}$, vector $b(x)\in\mathbb{R}^d$, and scalar $c(x)$. The operator is defined as

\[ \mathcal{L}u = -\sum_{i,j = 1}^d a_{ij}(x)\partial_{x_i,x_j}^2u + \sum_{i = 1}^db_i(x)\partial_{x_i}u + c(x)u.\]

A, b and c are space-dependent functions returning a matrix, a vector, and a scalar, respectively. The matrix A should be symmetric and positive definite for any input x. The operator can be called with a RadialSymmetricKernel and points x and y to evaluate the operator of the kernel at x - y. It can also be called with an Interpolation object and a point x to evaluate the elliptic operator of the interpolation at x.

Gradient()

The gradient operator. It can be called with a RadialSymmetricKernel and points x and y to evaluate the gradient of the kernel at x - y. It can also be called with an Interpolation object and a point x to evaluate the gradient of the interpolation at x.

Identity()

The identity operator, i.e. the differential operator of order zero. Applied to a function it returns its value, and like the other operators it can be called with a RadialSymmetricKernel and points x and y to evaluate the kernel at x - y, or with an Interpolation object and a point x to evaluate the interpolation at x. This lets plain evaluation reuse the same interface as the differential operators (for example in the RBF-FD weight computation, see local_weights).

Laplacian()

The Laplacian operator. It can be called with a RadialSymmetricKernel and points x and y to evaluate the Laplacian of the kernel at x - y. It can also be called with an Interpolation object and a point x to evaluate the Laplacian of the interpolation at x.

PartialDerivative(i)

Partial derivative operator with respect to the i-th component. The operator can be called with a RadialSymmetricKernel and points x and y to evaluate the derivative of the kernel at x - y. It can also be called with an Interpolation object and a point x to evaluate the first partial derivative of the interpolation at x in the i-th direction.

AdvectionDiffusionEquation(diffusivity, advection_velocity, f)

Advection-diffusion equation with diffusivity diffusivity and advection velocity advection_velocity. The advection-diffusion equation is defined as

\[ \partial_t u + \mathbf{a}\cdot\nabla u = \kappa\Delta u + f,\]

where $\mathbf{a}$ is the advection velocity, $\kappa$ is the diffusivity, and $f$ is the right-hand side, which can be a time- and space-dependent function or a vector.

AdvectionEquation(advection_velocity)

Advection equation with advection velocity advection_velocity. The advection equation is defined as

\[ \partial_t u + \mathbf{a}\cdot\nabla u = f,\]

where $\mathbf{a}$ is the advection velocity and $f$ a source term.

EllipticEquation(A, b, c, f)

Libear second-order elliptic equation with matrix A, vector b, and scalar c and right-hand side f. The elliptic equation is defined as

\[ \mathcal{L}u = -\sum_{i,j = 1}^d a_{ij}(x)\partial_{x_i,x_j}^2u + \sum_{i = 1}^db_i(x)\partial_{x_i}u + c(x)u = f,\]

where A, b and c are space-dependent functions returning a matrix, a vector, and a scalar, respectively.

See also EllipticOperator.

HeatEquation(diffusivity, f)

Heat equation with thermal diffusivity diffusivity. The heat equation is defined as

\[ \partial_t u = \kappa\Delta u + f,\]

where $\kappa$ is the thermal diffusivity and $f$ is the right-hand side, which can be a time- and space-dependent function or a vector.

PoissonEquation(f)

Poisson equation with right-hand side f, which can be a space-dependent function or a vector. The Poisson equation is defined as

\[ -\Delta u = f\]

See also Laplacian.

AbstractSpatialMethod

Abstract type tagging the spatial discretization strategy.

Collocation()

Global collocation strategy (Kansa method).

See also AbstractSpatialMethod and RBFFD.

RBFFD()

Local radial basis function finite difference strategy.

The local_basis keyword of SpatialDiscretization (or the local_basis field of RBFFDBasis) selects the algorithm used consistently for weight computation, matrix assembly, and interpolant evaluation:

• RBFFDLagrangeBasis() (default): precompute the local cardinal functions and evaluate w_k = 𝓛ℓ_k(x_i).
• RBFFDStandardBasis(): cache the factorization of each local kernel/polynomial system and solve M w = rhs, rhs_k = 𝓛K(x_i, x_k) (plus polynomial rows).

Both give the same weights up to numerical precision.

See also AbstractSpatialMethod and Collocation.

Semidiscretization(spatial_discretization, initial_condition)
Semidiscretization(equations, nodeset_inner, boundary_condition, nodeset_boundary, [centers,]
                   initial_condition, kernel = GaussKernel{dim(nodeset_inner)}())

Semidiscretization of a partial differential equation with Dirichlet boundary conditions and initial condition initial_condition. The boundary_condition function can be time- and space-dependent. The initial_condition function is time- and space-dependent to be able to reuse it as analytical solution if available. If no analytical solution is available, the time variable can be ignored in the initial_condition function. If a basis is passed via SpatialDiscretization, this basis is used consistently in the semidiscretization. For convenience constructors with centers and kernel, a StandardBasis is used and by default centers is set to merge(nodeset_inner, nodeset_boundary). The number of basis functions must be equal to the number of inner and boundary nodes because OrdinaryDiffEq.jl does not support DAEs with rectangular mass matrices.

If method in spatial_discretization is Collocation, the semidiscretization is built by collocating the PDE at the inner nodes and the boundary conditions at the boundary nodes, leading to a system of ODEs for the coefficients of the kernel interpolant. If method is RBFFD, the semidiscretization is built by applying the RBF-FD operator matrix to the vector of nodal values, leading to a system of ODEs for the nodal values themselves.

See also SpatialDiscretization, semidiscretize.

SpatialDiscretization(equations, nodeset_inner, boundary_condition, nodeset_boundary, basis)
SpatialDiscretization(equations, nodeset_inner, boundary_condition, nodeset_boundary, method, basis)
SpatialDiscretization(equations, nodeset_inner, boundary_condition, nodeset_boundary,
                      [centers,] kernel = GaussKernel{dim(nodeset_inner)}())
SpatialDiscretization(equations, nodeset_inner, boundary_condition, nodeset_boundary,
                      RBFFD(), kernel = GaussKernel{dim(nodeset_inner)}();
                      stencil_selection, m = order(kernel),
                      local_basis = RBFFDLagrangeBasis())

Spatial discretization of a partial differential equation with Dirichlet boundary conditions. The nodeset_inner are the nodes in the domain and nodeset_boundary are the nodes on the boundary. The boundary_condition is a function describing the Dirichlet boundary conditions. The centers are the centers of the kernel functions. By default, centers is set to merge(nodeset_inner, nodeset_boundary). Otherwise, a least squares problem is solved.

Uses method to select the spatial discretization strategy. If method is not given, it is selected based on the type of basis. If basis is an RBFFDBasis, method is set to RBFFD. Otherwise, method is set to Collocation.

See also Semidiscretization, solve_stationary.

semidiscetize(semi::Semidiscretization, tspan)

Wrap a Semidiscretization object into an ODEProblem object from SciMLBase.jl with time span tspan.

solve_stationary(spatial_discretization; linsolve = nothing)

Solve a stationary partial differential equation discretized as spatial_discretization with Dirichlet boundary conditions by non-symmetric collocation (Kansa method). Returns an Interpolation object. The linsolve keyword argument can be used to specify a linear solver from LinearSolve.jl for the linear system. If linsolve = nothing, the default backslash operator is used.

For kernel collocation, the kernel space is augmented with polynomials up to order(basis), so that conditionally positive definite kernels are augmented automatically, while strictly positive definite kernels (order == 0) use no polynomials. RBF-FD discretizations handle polynomial augmentation locally via the stencils.

See also SpatialDiscretization.

AbstractStencilSelection

Abstract type for stencil selection strategies in RBF-FD. Defines how neighboring points are selected for each local finite difference stencil.

KNearestNeighbors(k::Int)

Stencil selection strategy using k-nearest neighbors. Each interior point uses its k nearest neighbors (by Euclidean distance) to form the local FD stencil. This ensures uniform sparsity.

See also AbstractStencilSelection and RadiusSearch.

RadiusSearch(radius)

Stencil selection strategy using fixed radius search with search radius. Each interior point uses all neighbors within a given Euclidean distance radius to form the local FD stencil. This allows variable stencil sizes but adapts to local point density.

See also AbstractStencilSelection and KNearestNeighbors.

select_neighbors(i, nodeset, stencil_selection)

Selects neighbors of nodeset[i] from nodeset using stencil_selection. Returns a tuple of (neighbor_indices, neighbor_nodes).

local_weights(basis::RBFFDBasis, i, point, op)

Compute the local RBF-FD weights for stencil i, i.e. the coefficients that map the nodal values on that stencil to op applied to the local interpolant, evaluated at point. Use Identity for plain evaluation, or any differential operator/equation otherwise.

Dispatches on basis.local_basis: for RBFFDLagrangeBasis the (operator applied to the) precomputed cardinal functions are evaluated at point; for RBFFDStandardBasis the cached local system is solved with the corresponding right-hand side.

Returns a weight vector (scalar operator) or matrix with one column per component (vector operator).

rbf_fd_weights(diff_op_or_pde, i, basis::RBFFDBasis)

Compute RBF-FD weights for node index i using the precomputed stencil data from basis. The algorithm is selected by basis.local_basis (see RBFFDLagrangeBasis and RBFFDStandardBasis); both yield the same weights mathematically.

Returns the weight vector (scalar operator) or matrix (vector operator).

See also local_weights.

pde_matrix(diff_op_or_pde, nodeset, basis::RBFFDBasis)

Assemble the sparse RBF-FD operator matrix. Each row j corresponds to one node in nodeset and contains the local stencil weights of the differential operator (or PDE) diff_op_or_pde, evaluated at nodeset[j] using the nearest center's stencil in basis. The result is an |nodeset| × |basis| sparse matrix, so nodeset may be any set of evaluation points and need not equal the centers of the RBFFDBasis basis.

See also differentiation_matrix, pde_boundary_matrix.

differentiation_matrix(diff_op_or_pde, basis, nodeset = centers(basis))
differentiation_matrix(diff_op_or_pde, centers, kernel, nodeset = centers)

Compute the differentiation matrix $D$ discretizing the differential operator (or PDE) $\mathcal{L}$. It maps the vector of nodal values to the values of $\mathcal{L}u$ evaluated at the nodes in nodeset,

\[ D = A_\mathcal{L} A^{-1},\]

where $A_\mathcal{L}$ is the pde_matrix of $\mathcal{L}$ evaluated at nodeset and $A$ the kernel_matrix of the basis. In contrast to operator_matrix, no boundary rows are added, so $D$ is the plain matrix of the operator and nodeset may be any set of evaluation points (it defaults to the centers of the basis). If a node set of centers and a kernel are given, the matrix is computed for the StandardBasis.

Polynomials up to order m (i.e. degree m - 1) are added to the kernel space; m defaults to order(basis). With polynomial augmentation the matrix is exact for the full kernel + polynomial space; e.g. with m >= 1 it reproduces constants, so D * ones(N) vanishes for operators that annihilate constants (such as the Laplacian).

For an RBFFDBasis the matrix is sparse and assembled directly from the local stencil weights instead (the polynomial order is taken from the basis).

See also operator_matrix, pde_matrix, and kernel_matrix.

interpolation_matrix(centers, kernel, ps, regularization = NoRegularization(); factorization_method = Symmetric)
interpolation_matrix(basis, ps, regularization; factorization_method = Symmetric)

Return the interpolation matrix for the basis, polynomials ps, and regularization. For the StandardBasis, the interpolation matrix is defined as

\[ A = \begin{pmatrix}K & P\\P^T & 0\end{pmatrix},\]

where $K$ is the regularize!d kernel_matrix and $P$ the polynomial_matrix. If a node set of centers and a kernel are given, the interpolation matrix is computed for the StandardBasis. Additionally, you can specify a factorization_method to use for the system matrix. By default, the system matrix is just wrapped as a Symmetric matrix, but you can, e.g., also explicitly use cholesky, lu, or qr factorization.

kernel_matrix(basis, nodeset = centers(basis))
kernel_matrix(nodeset1[, nodeset2], kernel)

Return the kernel matrix for the nodes and kernel. The kernel matrix is defined as

\[ A_{ij} = b_j(x_i),\]

where $b_i$ are the basis function in the basis and x_i are the nodes in the nodeset. If two nodesets and a kernel are given, the kernel matrix is computed for the StandardBasis meaning

\[ A_{ij} = K(\xi_j, x_i),\]

where $\xi_j$ are the nodes/centers in nodeset1, $x_i$ are the nodes in nodeset2, and K is the kernel. If nodeset2 is not given, it defaults to nodeset1.

least_squares_matrix(basis, nodeset, ps, regularization = NoRegularization();
                     factorization_method = Matrix)
least_squares_matrix(centers, nodeset, kernel, ps, regularization = NoRegularization();
                     factorization_method = Matrix)

Return the least squares matrix for the basis, nodeset, polynomials ps, and regularization. For the StandardBasis, the least squares matrix is defined as

\[ A = \begin{pmatrix}K & P_1\\P_2^T & 0\end{pmatrix},\]

where $K$ is the regularize!d kernel_matrix, $P_1$ the polynomial_matrix for the nodeset and $P_2$ the polynomial_matrix for the centers. If a nodeset and kernel are given, the least squares matrix is computed for the StandardBasis. Additionally, you can specify a factorization_method to use for the system matrix. By default, the system matrix is not factorized, but you can, e.g., also explicitly use the qr factorization.

nearest_node_index(y, nodeset)

Return the index of the node in nodeset that is closest to y in Euclidean norm.

operator_matrix(diff_op_or_pde, nodeset_inner, nodeset_boundary, basis)
operator_matrix(diff_op_or_pde, nodeset_inner, nodeset_boundary, kernel)

Compute the operator matrix $L$ discretizing $\mathcal{L}$ for a given kernel. The operator matrix is defined as

\[ L = A_\mathcal{L} A^{-1},\]

where $A_\mathcal{L}$ is the matrix of the differential operator (defined by the equations), and $A$ the kernel matrix. If a basis is passed, A is built via kernel_matrix for that basis on merge(nodeset_inner, nodeset_boundary) and $A_\mathcal{L}$ via pde_boundary_matrix for the same basis.

For an RBFFDBasis, the operator matrix is assembled directly from the local stencil weights instead.

In contrast to differentiation_matrix, the boundary rows act as the identity (the nodal values are imposed there), so this is the system operator of the boundary value problem rather than the plain matrix of $\mathcal{L}$.

Polynomials up to order m (i.e. degree m - 1) are added to the kernel space; m defaults to order(basis), so conditionally positive definite kernels are augmented automatically, while strictly positive definite kernels (order == 0) use no polynomials.

See also differentiation_matrix, pde_boundary_matrix, and kernel_matrix.

pde_boundary_matrix(diff_op_or_pde, nodeset_inner, nodeset_boundary, [centers,] kernel)
pde_boundary_matrix(diff_op_or_pde, nodeset_inner, nodeset_boundary, basis)

Assemble the PDE + boundary matrix

\[ \begin{pmatrix}\tilde A_\mathcal{L}\\\tilde A\end{pmatrix},\]

where $\tilde A_\mathcal{L} =$pde_matrix evaluates the differential operator at the inner nodes and $\tilde A =$kernel_matrix evaluates the basis at the boundary nodes. With a kernel and optional centers (defaulting to merge(nodeset_inner, nodeset_boundary)) the entries are

\[ (\tilde A_\mathcal{L})_{ij} = \mathcal{L}K(x_i, \xi_j), \quad \tilde A_{ij} = K(x_i^b, \xi_j),\]

where $x_i$ are inner nodes, $x_i^b$ boundary nodes, and $\xi_j$ centers. With a basis the same structure holds with $b_j$ replacing $K(\cdot, \xi_j)$.

This method is used for all basis types, including RBFFDBasis, where pde_matrix and kernel_matrix use nearest-center stencil lookup. When |centers(basis)| < |nodeset_inner| + |nodeset_boundary|, the resulting system is overdetermined and solve_stationary solves it in the least-squares sense.

See also pde_matrix and kernel_matrix.

pde_boundary_matrix(diff_op_or_pde, nodeset_inner, nodeset_boundary, basis, ps)

Polynomial-augmented variant of pde_boundary_matrix. With the polynomials ps the system matrix becomes the saddle-point matrix

\[ \begin{pmatrix}\tilde A_\mathcal{L} & \mathcal{L}P_\mathrm{i}\\ \tilde A & P_\mathrm{b}\\ P^T & 0\end{pmatrix},\]

where $\mathcal{L}P_\mathrm{i}$ (pde_polynomial_matrix) and $P_\mathrm{b}$ (polynomial_matrix) are the operator/evaluation of the polynomials at the inner and boundary nodes, and $P$ the polynomial_matrix of the centers enforcing $\sum_j c_j p_k(\xi_j) = 0$. If ps is empty, this reduces to the unaugmented pde_boundary_matrix.

pde_matrix(diff_op_or_pde, nodeset1, nodeset2, kernel)
pde_matrix(diff_op_or_pde, nodeset, basis)

Compute the matrix of a partial differential equation (or differential operator) with a given kernel. The matrix is defined as

\[ (\tilde A_\mathcal{L})_{ij} = \mathcal{L}K(x_i, \xi_j),\]

where $\mathcal{L}$ is the differential operator (defined by the equations), $K$ the kernel, $x_i$ are the nodes in nodeset1 and $\xi_j$ are the nodes in nodeset2. If a basis is given, the matrix is defined as

\[ (\tilde A_\mathcal{L})_{ij} = \mathcal{L}b_j(x_i),\]

where $b_j$ are the basis functions in the basis and $x_i$ are the nodes in nodeset.

See also differentiation_matrix, pde_boundary_matrix.

pde_polynomial_matrix(diff_op_or_pde, nodeset, ps)

Return the matrix of a differential operator (or PDE) applied to the polynomials ps,

\[ (A_{\mathcal{L}}^p)_{ik} = \mathcal{L}p_k(x_i),\]

where $x_i$ are the nodes in nodeset and $p_k$ the polynomials. This is the polynomial counterpart of pde_matrix used for polynomial augmentation.

polynomial_matrix(nodeset, ps)

Return the polynomial matrix for the nodeset and polynomials. The polynomial matrix is defined as

\[ A_{ij} = p_j(x_i),\]

where $x_i$ are the nodes in the nodeset and $p_j$ the polynomials.

AliveCallback(io::IO = stdout; interval::Integer=0, dt=nothing)

Inexpensive callback showing that a simulation is still running by printing some information such as the current time to the screen every interval time steps or after a time of dt in terms of integration time by adding additional (shortened) time steps where necessary (note that this may change the solution).

SummaryCallback(io::IO = stdout)

Create and return a callback that resets the timer at the beginning of a simulation and prints the timer values at the end of the simulation.

SaveSolutionCallback(; interval::Integer=0,
                     dt=nothing,
                     save_initial_solution=true,
                     save_final_solution=true,
                     output_directory="out",
                     extra_functions=(),
                     keys=append!(["itp"], "value_" .* string.(eachindex(extra_functions))))

Save the current numerical solution in regular intervals in VTK format as a Paraview Collection (.pvd). Either pass interval to save every interval time steps or pass dt to save in intervals of dt in terms of integration time by adding additional (shortened) time steps where necessary (note that this may change the solution). The interpolation object will always be saved at the inner and boundary nodes of the corresponding Semidiscretization. You can pass extra functions (time- and space-dependent) or vectors to save at these nodes via extra_functions. The corresponding keys in the .vtu files can be specified by keys.

See also add_to_pvd, vtk_save.

add_to_pvd(filename, pvd, time, nodeset::NodeSet, functions_or_vectors...;
           keys = "value_" .* string.(eachindex(functions_or_vectors)))

Same as vtk_save, but appends the data to a Paraview collection file pvd at time time. In contrast to vtk_save, the functions are time- and space-dependent.

vtk_read(filename)

Read a set of nodes from a VTK file and return them as a NodeSet. Note that the data will always be returned as a 3D NodeSet, even if the data is 1D or 2D. The point data will be returned as a dictionary with the keys being the names of the data arrays in the VTK file.

See also vtk_save.

vtk_save(filename, nodeset::NodeSet, functions_or_vectors...;
         keys = "value_" .* string.(eachindex(functions_or_vectors)))

Save a NodeSet to a VTK file. You can optionally pass a list of space-dependent functions or vectors to save the values of the functions at the nodes. The functions can also be passed as Interpolation or directly as vectors. The optional keyword argument keys is used to specify the names of the data arrays in the VTK file.

See also add_to_pvd, vtk_read.

default_example()

Return the path to an example that can be used to quickly see KernelInterpolation.jl in action. See also examples_dir and get_examples.

examples_dir()

Return the directory where the example files provided with KernelInterpolation.jl are located.

Examples

readdir(examples_dir())

get_examples()

Return a list of all examples that are provided by KernelInterpolation.jl. See also examples_dir and default_example.