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 \rightarrow \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}:
 1.9788455979766866
 2.9056068768448933
 2.969871299314206
 2.780167775143936
 3.3365164304658177
 2.4841717588644947
 3.8556006321068583
 3.2944912584949115
 0.8615173945037004
 4.688624353530846
 ⋮
 2.9169949794563705
 0.2812060093062602
 4.127689010892639
 1.19486141281526
 3.059807515220657
 4.1286097202765735
 4.422494726299439
 3.2787612203593257
 1.009917877623197

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}:
 2.866215709104487
 1.5164948276797192
 1.5801699460466934
 1.03183235153639
 2.567400442027525

and compare the error

abs.(y_true .- y_new)

5-element Vector{Float64}:
 3.862532960141607e-9
 6.336995752365127e-6
 0.00018309988925668819
 2.4883883833481946e-6
 1.1906733377031742e-9

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}:
 3.385523135168711e-6
 6.3370505549720235e-6
 0.0001830961695457045
 2.2612928878285032e-6
 1.0959286811385027e-9

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

We can also retrieve the gradient. This is really just a convenience wrapper around Partial.

grad = gradient(x)

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

# error
all(df_x.(x) .≈ grad(y)[1])

true

all(df_y.(x) .≈ grad(y)[2])

true