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


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

1


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)


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


Interpolate the gradients of the cubic polynomial

1


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


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


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

1
2
3
4
5

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


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


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


Compute the interpolated gradients of the target function.

1
2
3
4
5

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


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


The gradients of the target function can be also interpolated.

1
2
3
4
5

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


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.
interp(spl, x)

1