Getting Started

Data must be an AbstractVector of point vectors — each point needs an inferrable dimension (e.g., SVector{2,Float64} from StaticArrays.jl).

using RadialBasisFunctions
using StaticArrays

Interpolation

Suppose we have a set of data where  , and we want to interpolate a function   

f(x) = 2*x[1]^2 + 3*x[2]
x = rand(SVector{2,Float64}, 1000)
y = f.(x)

and now we can build the interpolator

interp = Interpolator(x, y)

Interpolator
├─Input type: StaticArraysCore.SVector{2, Float64}
├─Output type: Float64
├─Number of points: 1000
└─Basis: Polyharmonic spline (r³) with degree 2 Monomial

and evaluate it at a new point

x_new = rand(SVector{2,Float64}, 5)
y_new = interp(x_new)
y_true = f.(x_new)

and compare the error

abs.(y_true .- y_new)

5-element Vector{Float64}:
 4.440892098500626e-16
 2.220446049250313e-16
 4.440892098500626e-16
 2.220446049250313e-16
 4.440892098500626e-16

The error is numerically zero because the default basis — PHS(3; poly_deg=2) — includes quadratic polynomial augmentation, which can represent our 2nd-order polynomial f exactly. Reducing the polynomial degree shows the effect:

interp = Interpolator(x, y, PHS(3; poly_deg=1))
y_new = interp(x_new)
abs.(y_true .- y_new)

5-element Vector{Float64}:
 3.0824116148409075e-8
 5.375260103690493e-7
 8.036068588346268e-9
 2.4177835289407312e-8
 1.5450165591346376e-9

Operators

Operators compute RBF-FD weights for differentiation on scattered data. Weights are built at construction and cached; invalidate the cache to trigger recomputation.

Partial Derivative

df_x_rbf = partial(x, 1, 1)

# define exact
df_x(x) = 4*x[1]

# error
all(abs.(df_x.(x) .- df_x_rbf(y)) .< 1e-10)

true

Laplacian

lap_rbf = laplacian(x)

# define exact
lap(x) = 4

# error
all(abs.(lap.(x) .- lap_rbf(y)) .< 1e-8)

true

Gradient / Jacobian

The jacobian function computes all partial derivatives. For scalar fields, this is the gradient. The gradient function is a convenience alias for jacobian.

op = jacobian(x)  # or equivalently: gradient(x)
result = op(y)    # Matrix of size (N, dim)

# define exacts
df_x(x) = 4*x[1]
df_y(x) = 3

# error - access columns for each partial derivative
all(df_x.(x) .≈ result[:, 1])

true

all(df_y.(x) .≈ result[:, 2])

true

Directional Derivative

Compute derivatives in any direction using directional. The direction can be constant or vary spatially:

using LinearAlgebra: normalize

# Constant direction (same for all points)
v = normalize([1.0, 1.0])
dir_op = directional(x, v)
result = dir_op(y)
typeof(result)

Vector{Float64} (alias for Array{Float64, 1})

The direction can also vary per-point, useful for computing normal derivatives:

# Spatially-varying direction (e.g., radial directions)
normals = map(normalize, x)
normal_deriv = directional(x, normals)
typeof(normal_deriv(y))

Vector{Float64} (alias for Array{Float64, 1})

For the common case of differentiating along outward normals, normal_derivative wraps directional and normalizes the vectors you pass it.

Custom & PDE Operators

Beyond the built-ins, the @operator macro lets you write PDE operators in mathematical notation — Helmholtz, diffusion, advection-diffusion — and call them directly with data points. See Building PDE Operators for the recipes, and Operators & Type Hierarchy for an in-depth guide to the operator system.

Regridding

Interpolate field values from one set of points to another using regrid:

# Target points (fine grid, different from original x)
x_fine = rand(SVector{2,Float64}, 500)

# Build regridding operator from x to x_fine
rg = regrid(x, x_fine)
y_fine = rg(y)
length(y_fine)

500

Operator Algebra

Operators can be combined using + and -:

# Create individual operators
∂x = partial(x, 1, 1)
∂y = partial(x, 1, 2)

# Combine them: ∂f/∂x + ∂f/∂y
combined = ∂x + ∂y
result = combined(y)
typeof(result)

Vector{Float64} (alias for Array{Float64, 1})

Enforcing Boundary Conditions

For PDE applications, operators support Hermite interpolation to enforce Dirichlet, Neumann, or Robin conditions at boundary nodes. See Boundary Conditions for the condition types and the hermite keyword.

Where to Next

• Operators & Type Hierarchy — the operator system, rank semantics, and virtual operators

• Building PDE Operators — assemble Helmholtz, diffusion, and advection-diffusion operators

• Quick Reference — data formats, basis options, and operator constructors at a glance

• Convergence & Parameter Selection — how to pick a basis, poly_deg, and stencil size

Current Limitations

2. Global interpolation: Interpolator currently uses all points globally. Local collocation support (like the operators use) is planned for future releases.

3. GPU weight computation: weight computation (stencil assembly and solve) currently runs on CPU only; a GPU-compatible dense solver is needed for full GPU support (#88). Built operators can be moved to GPU for evaluation — see GPU Evaluation in the Quick Reference.