$ \newcommand\vu{\mathbf{u}} \newcommand\vv{\mathbf{v}} \newcommand\vw{\mathbf{w}} \newcommand\vi{\mathbf{i}} \newcommand\vj{\mathbf{j}} \newcommand\vk{\mathbf{k}} \newcommand\vo{\mathbf{o}} \newcommand\vr{\mathbf{r}} \newcommand\vt{\mathbf{t}} \newcommand\vx{\mathbf{x}} \newcommand\vy{\mathbf{y}} \newcommand\vz{\mathbf{z}} \newcommand\va{\mathbf{a}} \newcommand\vb{\mathbf{b}} \newcommand\vc{\mathbf{c}} \newcommand\ve{\mathbf{e}} \newcommand{\vE}{\mathbf{E}} \newcommand{\vS}{\mathbf{S}} \newcommand{\vk}{\mathbf{k}} \newcommand{\vq}{\mathbf{q}} \newcommand{\vzero}{\mathbf{0}} \newcommand{\vomega}{\mathbf{ω}} \newcommand{\mI}{\mathbf{I}} \newcommand{\mM}{\mathbf{M}} \newcommand{\mR}{\mathbf{R}} \newcommand{\mP}{\mathbf{P}} \newcommand{\mK}{\mathbf{K}} \newcommand{\mW}{\mathbf{W}} \DeclareMathOperator{\Tr}{Tr} \newcommand{\abs}[1]{\left\lvert {#1} \right\rvert} \newcommand{\norm}[1]{\left\lVert {#1} \right\rVert} \newcommand{\ensemble}[1]{\left\langle {#1} \right\rangle} \newcommand{\der}[2]{\frac{\partial {#1}}{\partial {#2}}} \newcommand{\Der}[2]{\frac{\delta {#1}}{\delta {#2}}} \newcommand{\lap}{\nabla^2} $

Perform cubic Hermite spline interpolation in Julia

Tutorial on CubicHermiteSpline.jl

Software · Tutorials Julia CubicHermiteSpline.jl PhaseDiagram.jl Algorithms

This is a tutorial on how to use the Julia package CubicHermiteSpline.jl, which performs a cubic Hermite spline interpolation on an array of data points, $(x_i, y_i)$, given that their associated gradients, $k_i=(dy/dx)_i$, are known in advance.

New

  • v0.3.0 can now perform bivariate cubic Hermite spline interpolation for 2D data points (regular and irregular grids are both supported).
  • v0.2.2 can now compute the 1st order derivative of the interpolated function.

Getting started

First, the package can be easily installed from the Julia REPL:

1
2

julia> using Pkg
julia> Pkg.add("CubicHermiteSpline")

Or in the package mode (typing ] in REPL):

1

(v1.5) pkg> add CubicHermiteSpline

Then you can load the package:

1
2
3
4

using CubicHermiteSpline
import Plots
using Plots: plot, plot!
using LaTeXStrings

Plotting related packages are also loaded. Let’s tweak the appearance of our plots:

1
2
3
4
5
6
7
8
9
10
11
12

Plots.default(size=(600, 370))
fntf = :Helvetica
titlefont = Plots.font(fntf, pointsize=12)
guidefont = Plots.font(fntf, pointsize=12)
tickfont = Plots.font(fntf, pointsize=9)
legendfont = Plots.font(fntf, pointsize=8)
Plots.default(fontfamily=fntf)
Plots.default(titlefont=titlefont, guidefont=guidefont, tickfont=tickfont, legendfont=legendfont)
Plots.default(minorticks=true)
Plots.default(linewidth=1.2)
Plots.default(foreground_color_legend=nothing)
Plots.default(legend=false)

Now we are ready to go.

Example 1: interpolating cubic polynomial

The cubic polynomial of the form

\[f(x) = ax^3 + bx^2 + cx + d\]

should be exactly interpolated by the cubic Hermite spline interpolation.

Below we use CubicHermiteSpline.jl to demonstrate this fact. First we define a typical cubic polynomial:

1

f(x) = x^3 - 3x^2 + 2x - 5;

Its gradient are available in an analytical form as

1

df(x) = 3x^2 - 6x + 2;

The exact cubic polynomial is evaluated at evenly spaced points.

1
2
3
4
5

a = 0.0
b = 2.5
x_exact = range(a, b, step=0.01);
y_exact = f.(x_exact);
k_exact = df.(x_exact);

The exact cubic polynomial is plotted.

1
2
3

xlabel = L"x"
ylabel = L"f(x)"
plot(x_exact, y_exact, xlabel=xlabel, ylabel=ylabel)
Cubic polynomial.

Generate a set of data points to be interpolated in the cubic polynomial curve on a evenly spaced grid:

1
2
3

x_input = range(a, b, step=0.5);
y_input = f.(x_input);
k_input = df.(x_input);

Create an interpolation instance:

1

spl = CubicHermiteSplineInterpolation(x_input, y_input, k_input);

Perform the actual cubic Hermite spline interpolation:

1
2

x_interp = range(a, b, step=0.01);
y_interp = spl(x_interp);

Plotting input data points, exact solution, and the interpolated result in a single plot shows that the interpolation is exact for cubic polynomials.

1
2
3

Plots.scatter(x_input, y_input, label="input data", legend=:topleft)
plot!(x_exact, y_exact, label="exact")
plot!(x_interp, y_interp, label="interpolated")
Interpolation on even grid.

Interpolate the gradients of the cubic polynomial

1

k_interp = spl(x_interp; grad=true);

The results are plotted.

1
2
3

Plots.scatter(x_input, k_input, xlabel=xlabel, ylabel=L"f'(x)", label="input data", legend=:topleft)
plot!(x_exact, k_exact, label="exact")
plot!(x_interp, k_interp, label="interpolated")
Interpolated gradients on even grid.

Do the interpolation again for data points on an irregular grid:

1
2
3
4
5
6
7
8
9
10
11
12

x_input = [a, 0.1, 0.6, 1.75, b];
y_input = f.(x_input);
k_input = df.(x_input);

spl = CubicHermiteSplineInterpolation(x_input, y_input, k_input);

x_interp = range(a, b, step=0.01);
y_interp = spl(x_interp);

Plots.scatter(x_input, y_input, label="input data", legend=:topleft)
plot!(x_exact, y_exact, label="exact")
plot!(x_interp, y_interp, label="interpolated")
Interpolation on uneven grid.

Also compute the 1st order derivative of the interpolated function.

1
2
3
4
5

k_interp = spl(x_interp; grad=true);

Plots.scatter(x_input, k_input, xlabel=xlabel, ylabel=L"f'(x)", label="input data", legend=:topleft)
plot!(x_exact, k_exact, label="exact")
plot!(x_interp, k_interp, label="interpolated")
Interpolated gradients on uneven grid.

Example 2: interpolating smooth functions

Here we use the spherical bessel function of the first kind:

\[j_1(x) = \frac{\sin x - x\cos x}{x^2}.\]

Below we show that the cubic Hermite spline performs impressively well for interpolating such smooth functions.

1
2
3
4
5
6
7
8
9
10
11
12

g(x) = (sin(x) - x*cos(x)) / x^2;

# dg(x) = ((x^2 - 2) * sin(x) + 2x*cos(x)) / x^3;
dg(x) = (sin(x) - 2*g(x)) / x;

a = 0.01
b = 8.01
x_exact = range(a, b, step=0.01);
y_exact = g.(x_exact);

ylabel = L"g(x)"
plot(x_exact, y_exact, xlabel=xlabel, ylabel=ylabel)
spherical bessel function of the first kind.

Perform cubic Hermite interpolation on an even grid:

1
2
3
4
5
6
7
8
9
10
11
12

x_input = range(a, b, step=0.5);
y_input = g.(x_input);
k_input = dg.(x_input);

spl = CubicHermiteSplineInterpolation(x_input, y_input, k_input);

x_interp = range(a, b, step=0.01);
y_interp = spl(x_interp);

Plots.scatter(x_input, y_input, xlabel=xlabel, ylabel=ylabel, label="input data", legend=:topright)
plot!(x_exact, y_exact, label="exact")
plot!(x_interp, y_interp, label="interpolated")
Interpolation on even grid.

Compute the interpolated gradients of the target function.

1
2
3
4
5

k_interp = spl(x_interp; grad=true);

Plots.scatter(x_input, k_input, xlabel=xlabel, ylabel=L"g'(x)", label="input data", legend=:topright)
plot!(x_exact, k_exact, label="exact")
plot!(x_interp, k_interp, label="interpolated")
Interpolated gradients on even grid.

It can be seen from above plot that interpolation on data points and gradients works perfectly for known data points on a dense even grid.

Interpolate the function on an unevenly spaced grid:

1
2
3
4
5
6
7
8
9
10
11
12

x_input = [a, 0.6, 1.5, 3.4, 4.5, 5.0, 6.2, 7.5, b];
y_input = g.(x_input);
k_input = dg.(x_input);

spl = CubicHermiteSplineInterpolation(x_input, y_input, k_input);

x_interp = range(a, b, step=0.01);
y_interp = spl(x_interp);

Plots.scatter(x_input, y_input, xlabel=xlabel, ylabel=ylabel, label="input data", legend=:topright)
plot!(x_exact, y_exact, label="exact")
plot!(x_interp, y_interp, label="interpolated")
Interpolation on uneven grid.

The gradients of the target function can be also interpolated.

1
2
3
4
5

k_interp = spl(x_interp; grad=true);

Plots.scatter(x_input, k_input, xlabel=xlabel, ylabel=L"g'(x)", label="input data", legend=:topright)
plot!(x_exact, k_exact, label="exact")
plot!(x_interp, k_interp, label="interpolated")
Interpolated gradients on uneven grid.

Note that the accuracy of the interpolation degrades near the extrema of the function if there is no input data point near that region.

Available methods

  • Interpolation
1
2
3
4
5

# `spl` is an instance of `CubicHermiteSplineInterpolation`.
# `x` can be a real number of an 1D array of real numbers.
# Following two methods are equivalent.
spl(x; grad=false)
interp(spl, x)
  • Interpolated gradients
1
2

spl(x; grad=true)
grad(spl, x)

Contribute

Acknowledgements

  • This work is partially supported by the General Program of the National Natural Science Foundation of China (No. 21873021) and Shanghai Pujiang Program (No. 18PJ1401200).