By: Yi-Xin Liu
Re-posted from: http://www.yxliu.group/2020/06/cubic-hermite-spline
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.
Table of Contents
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)
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
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")
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
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")
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
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")
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
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")
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
- Star the package on github.com.
- File an issue or make a pull request on github.com.
- Contact me via email lyx@fudan.edu.cn.
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).
Tutorial on CubicHermiteSpline.jl was originally published by Yi-Xin Liu at Yi-Xin Liu on June 19, 2020.