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 = [SVector{2}(rand(2)) for _ in 1:1000]
y = f.(x)

1000-element Vector{Float64}:
 3.9717591125729315
 1.610067780790475
 1.8519526444552061
 0.39339965716380576
 3.4536975283290383
 4.148781785884552
 4.3741486475240965
 2.7993612345776
 2.4649904115225034
 1.828733760706511
 ⋮
 2.962191828801693
 1.342206190360624
 0.6483395520247152
 3.913924982933284
 4.1201228963632675
 1.9619814862506186
 1.582994258828927
 2.671072494148113
 1.6131820710914793

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 = [SVector{2}(rand(2)) for _ in 1:5]
y_new = interp(x_new)
y_true = f.(x_new)

5-element Vector{Float64}:
 2.75033945084214
 2.7180217589112434
 2.0698756902753277
 1.5821555513268122
 1.7830170345726586

and compare the error

abs.(y_true .- y_new)

5-element Vector{Float64}:
 4.440892098500626e-16
 4.440892098500626e-16
 8.881784197001252e-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}:
 8.268190887505966e-8
 2.2570407764987976e-9
 9.435517256406456e-9
 1.1529829958334403e-5
 7.116978117949202e-10

Operators

Operators compute RBF-FD weights for differentiation on scattered data. They're lazy — weights are built on first evaluation and cached.

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 = [normalize(collect(p)) for p in x]
normal_deriv = directional(x, normals)
typeof(normal_deriv(y))

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

Custom Operators

Define your own differential operators using custom. The function should accept a basis and return a callable (x, xc) -> value. Here's an example that creates an interpolation-like operator:

# Custom operator that evaluates the basis function
op = custom(x, basis -> (x, xc) -> basis(x, xc))
typeof(op)

RadialBasisOperator{Custom{Main.var"#12#13"}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{StaticArraysCore.SVector{2, Float64}}, Vector{StaticArraysCore.SVector{2, Float64}}, Vector{Vector{Int64}}, PHS3{Int64}, KernelAbstractions.CPU}

For more complex differential operators, use operator algebra (see below) to combine built-in operators.

Regridding

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

# Target points (fine grid, different from original x)
x_fine = [SVector{2}(rand(2)) for _ in 1: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})

Boundary Conditions (Hermite Interpolation)

For PDE applications, operators support Hermite interpolation with boundary conditions. This is useful when you need to enforce Dirichlet, Neumann, or Robin conditions at boundary nodes.

Boundary Condition Types

• Dirichlet() - Value specified:  

• Neumann() - Normal derivative specified:  

• Robin(α, β) - Mixed condition:   

• Internal() - Interior point (no boundary condition)

Example with Hermite Interpolation

using LinearAlgebra: norm

# Define boundary information
is_boundary = [norm(p) > 0.9 for p in x]  # Points near unit circle boundary
boundary_indices = findall(is_boundary)

# Create boundary conditions (Dirichlet on boundary)
bcs = [Dirichlet() for _ in boundary_indices]

# Normal vectors at boundary points
normals = [normalize(collect(x[i])) for i in boundary_indices]

# Build operator with Hermite interpolation
lap_hermite = laplacian(x; hermite=(
    is_boundary=is_boundary,
    bc=bcs,
    normals=normals
))
typeof(lap_hermite)

RadialBasisOperator{Laplacian, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{StaticArraysCore.SVector{2, Float64}}, Vector{StaticArraysCore.SVector{2, Float64}}, Vector{Vector{Int64}}, PHS3{Int64}, KernelAbstractions.CPU}

Advanced: Virtual Operators

Virtual operators (∂virtual) use finite difference formulas on interpolated values at offset points. This can be useful for certain numerical schemes:

# Virtual partial derivative in x-direction with spacing Δ=0.01
virtual_dx = ∂virtual(x, 1, 0.01)
result = virtual_dx(y)
typeof(result)

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

Current Limitations

1. Data format: The package requires AbstractVector{<:AbstractVector} input (not matrices). Each point must have inferrable dimension, e.g., SVector{2,Float64} from StaticArrays.jl. Matrix input support is planned.

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: Operators and interpolators can be moved to GPU via Adapt.jl (e.g., cu(operator)) for evaluation, but weight computation (stencil assembly and solve) currently runs on CPU only. A GPU-compatible dense solver is needed for full GPU support (#88).