$ \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} $

What the package can do?

Tutorial on CubicHermiteSpline.jl

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 at each point, $k_i$, are known in advance.

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.4) 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

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 cubic polynomial is plotted.

1
2
3
4
5
6
a = 0.0
b = 2.5
x_exact = range(a, b, step=0.01);
y_exact = f.(x_exact);

plot(x_exact, y_exact)
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.

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.

Example 2: interpolating smooth functions

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

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
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);

plot(x_exact, y_exact)
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, label="input data", legend=:topright)
plot!(x_exact, y_exact, label="exact")
plot!(x_interp, y_interp, label="interpolated")
Interpolation on even grid.

It can be seen from above plot that the interpolation 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, label="input data", legend=:topright)
plot!(x_exact, y_exact, label="exact")
plot!(x_interp, y_interp, label="interpolated")
Interpolation 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.

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).
  • The captured output text blocks in this post is automatically generated by using Literate.jl.