One-dimensional interpolation and differentiation

In this tutorial, we will create a simple one-dimensional interpolation, investigate how to tune the interpolation method, and show how to apply differential operators on the resulting Interpolation object.

Define problem setup and perform interpolation

We start by defining a simple one-dimensional interpolation problem. We will interpolate the oscillatory function

\[f(x) = \exp(\sin(2x^2)) + 0.1(x - \pi/2)^2\]

between $x = -3$ and $x = 3$. For simplicity, we take 25 equidistant points in the interval $[-3, 3]$ as interpolation points.

using KernelInterpolation
f(x) = exp(sin(2*x[1]^2)) + 0.1*(x[1] - pi/2)^2
x_min = -3
x_max = 3
N = 25
nodeset = NodeSet(LinRange(x_min, x_max, N))
values = f.(nodeset)

25-element Vector{Float64}:
 2.561118346612702
 3.6010482803853834
 2.592967848165464
 1.9848324553623153
 3.9645665583549765
 1.957023058914906
 1.3192182268719703
 1.8124194938335236
 3.143477103401214
 3.003796380105595
 ⋮
 2.5151585726832555
 1.0270213304360754
 0.3767404307950323
 0.8574656301584783
 2.7079294969190593
 0.5711157612469084
 1.0221715213705675
 1.8731723209109972
 0.676162754458826

Next, we choose the kernel (radial basis function) for the interpolation. We use the Gaussian kernel with a fixed shape parameter of 0.5 and interpolate the function values.

kernel = GaussKernel{1}(shape_parameter = 0.5)
itp = interpolate(nodeset, values, kernel)

Interpolation with 25 nodes, kernel GaussKernel{1}(shape_parameter = 0.5) and polynomial of order 0.

Let's plot the interpolated function and the original function on a finer grid to see how well the interpolation works.

using Plots
many_nodes = NodeSet(LinRange(x_min, x_max, 200))
plot(many_nodes, f, label = "Original function")
plot!(many_nodes, itp)

Uhh, that doesn't look too good. What happened?

Finding a well-suited interpolation method

We used the GaussKernel with a rather small shape parameter of 0.5, which leads to an ill-conditioned linear system of equations. We can inspect the condition number of the interpolation matrix to confirm this.

using LinearAlgebra
A = system_matrix(itp)
cond(A)

6.239147369601393e17

Here, we used the system_matrix function to obtain the interpolation matrix A and calculated the condition number of the matrix. For this specific example the system matrix simply is the kernel_matrix, but for more sophisticated interpolations the system matrix contains additional parts like the polynomial augmentation. The condition number is a measure of how well-conditioned the matrix is. A large condition number indicates that the matrix is ill-conditioned, which usually leads to high numerical errors. To avoid this, we have different options. We can either increase the shape parameter of the kernel or we can use a different kernel. The GaussKernel is known to be rather ill-conditioned and other kernels like the WendlandKernel usually lead to better condition numbers. Here, we choose to increase the shape parameter of the Gaussian kernel to 1.8, which makes the interpolation more localized. Note, however, that you might need to choose another kernel if you increase the number of interpolation points.

kernel = GaussKernel{1}(shape_parameter = 1.8)
itp = interpolate(nodeset, values, kernel)
plot(many_nodes, f, label = "Original function")
plot!(many_nodes, itp)

We can see a much better agreement between the original function and the interpolated function. We still observe some undershoots, but this is expected due to the oscillatory nature of the function and the limited number of interpolation points. Let's confirm that increasing the shape parameter improved the condition number of the interpolation matrix.

A = system_matrix(itp)
cond(A)

57037.29656737549

Indeed, the condition number is much smaller than before!

Applying differential operators

Sometimes, we are not only interested in interpolating a function, but also in computing its derivatives. Remember that in the simplest case, where no polynomial augmentation is used, the interpolation itp represents a linear combination

\[s(x) = \sum_{j = 1}^N c_j\phi(\|x - x_j\|_2)\]

with $\phi$ given by the radial basis function, in this case the Gaussian. Because we know $\phi$ and its derivatives, we can compute the derivatives of $s$ by differentiating the kernel function. For a general dimension $d$, the partial derivative in the $i$-th direction, $i\in\{1,\ldots,d\}$, of the interpolation is then given by

\[\frac{\partial s}{\partial x_i}(x) = \sum_{j = 1}^N c_j\frac{\partial \phi}{\partial x_i}(\|x - x_j\|_2).\]

KernelInterpolation.jl already provides some common differential operators. For example, we can compute the first derivative of the interpolation itp at a specific point x by using the PartialDerivative operator.

d1 = PartialDerivative(1)
x = 0.0
itp_dx = d1(itp, x)

-0.3128030378004036

Let's plot the first derivative of the interpolated function and compare it to the analytical first derivative.

itp_dx_many_nodes = d1.(Ref(itp), many_nodes)
f_dx(x) = 4*exp(sin(2*x[1]^2))*x[1]*cos(2*x[1]^2) + 0.2*x[1] - pi/10
plot(many_nodes, f_dx, label = "Derivative of original function")
plot!(many_nodes, itp_dx_many_nodes, label = "Derivative of interpolated function")