Getting Started

First, let's load the package along with the StaticArrays.jl package which we use for each data point

using RadialBasisFunctions
using StaticArrays

Interpolation

Suppose we have a set of data $\mathbf{x}$ where ${\mathbf{x}}_i\in {\mathbb{R}}^2$, and we want to interpolate a function $f:{\mathbb{R}}^2\to \mathbb{R}$

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}:
 2.0006792591140914
 1.5445054348176694
 3.273748991537368
 1.6658071504763172
 1.445402970709216
 4.014623483005918
 0.6658329758830654
 2.532894729350871
 1.0404216010508114
 2.1270888118075244
 ⋮
 2.052812915133434
 4.720189219358042
 1.1569973923633432
 2.069780671719789
 2.547387320356096
 0.8037488077188456
 1.6397828809254094
 1.4588531200647892
 2.0374114018103335

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

5-element Vector{Float64}:
 1.3585462831253774
 3.0251392593462434
 2.5695487349864563
 4.055040497756922
 1.8580288682902035

and compare the error

abs.(y_true .- y_new)

5-element Vector{Float64}:
 4.440892098500626e-16
 8.881784197001252e-16
 8.881784197001252e-16
 0.0
 4.440892098500626e-16

Wow! The error is numerically zero! Well... we set ourselves up for success here. Interpolator (along with RadialBasisOperator) has an optional argument to provide the type of radial basis including the degree of polynomial augmentation. The default basis is a cubic polyharmonic spline with 2nd degree polynomial augmentation (which the constructor is PHS(3, poly_deg=2)) and given the underlying function we are interpolating is a 2nd order polynomial itself, we are able to represent it exactly (up to machine precision). Let's see what happens when we only use 1st order polynomial augmentation

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

5-element Vector{Float64}:
 4.5903730416796407e-7
 7.201775620302442e-9
 2.3064976595321696e-8
 1.37589984205988e-6
 2.6868224429765775e-6

Operators

This package also provides an API for operators. There is support for several built-in operators along with support for user-defined operators. Currently, we have implementations for

• partial derivative (1st and 2nd order)

• laplacian

• gradient

but we plan to add more in the future. Please make and issue or pull request for additional operators.

Partial Derivative

We can take the same data as above and build a partial derivative operator with a similar construction as the interpolator. For the partial we need to specify the order of differentiation we want along with the dimension for which to take the partial. We can also supply some optional arguments such as the basis and number of points in the stencil. The function inputs order is partial(x, order, dim, basis; k=5)

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

Building a laplacian operator is as easy as

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

Current Limitations

1. Data format: The package requires Vector{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.