Recently I had a chance to test two interpolation packages, Interpolations.jl and Dierckx.jl. I found a difference between the two packages in the way they treat the boundary points when evaluating derivatives. I think this is an important difference that the users need to be aware of. And yet, they are not clearly discussed in their documentations, so here I report it.
First off, here in the version information of the two packages:
(@v1.6) pkg> st Interpolations Status `~/.julia/environments/v1.6/Project.toml` [a98d9a8b] Interpolations v0.13.2 (@v1.6) pkg> st Dierckx Status `~/.julia/environments/v1.6/Project.toml` [39dd38d3] Dierckx v0.5.1
And here is the test code to interpolate the data points sampled from
f(x) = x^3 in
-1 ≤ x ≤ 1:
using Interpolations using Dierckx xmin, xmax = -1, 1 f(x) = x^3 xs = range(xmin, xmax, length=21) ys = f.(xs) # Interpolators fi = CubicSplineInterpolation(xs, ys, extrapolation_bc=Flat()) # use Interpolations fd = Spline1D(xs, ys, bc="nearest") # use Dierckx # First derivatives of interpolators fi⁽¹⁾(x) = only(Interpolations.gradient(fi, x)) fd⁽¹⁾(x) = derivative(fd, x) # Second derivatives of interpolators fi⁽²⁾(x) = only(Interpolations.hessian(fi, x)) fd⁽²⁾(x) = derivative(fd, x, nu=2)
Note that the boundary conditions
Dierckx.jl are equivalent boundary conditions that extrapolate the estimate function with the values of the interpolator evaluated at the boundary points.
Then, I plot the interpolators, their first derivatives, and their second derivatives at the sampling points. I also plot the exact curve y = x³ from which the sampling data are sampled, and its first and second derivatives y = 3x² and y = 6x. Here is the result (x-coordinates of the dots in the plots correspond to the sampling points):
In the first panel of the figure, we observe that the interpolators produced by the two packages are visually indistinguishable and perfectly aligned with the exact curve y = x³. (This is expected, because cubic spline passes though knots.)
However, the two packages start deviating when evaluating derivatives. In the second panel of the figure, we observe that the two packages produce different first derivatives around the boundary points. Specifically,
Dierckx.jl produces first derivatives perfectly aligned with the exact curve y = 3x², but
Interpolations.jl doesn’t near the boundary points
The deviation between the two packages becomes very noticeable when evaluating second derivatives. In the third panel of the figure,
Dierckx.jl still produces second derivatives perfectly aligned with the the exact curve y = 6x, but
Interpolations.jl produces very different values around the boundary points of the sampling range.
Though I didn’t have time to take a close look at the code, it seems that
Interpolations.jl takes the extrapolated values outside the sampling range into account when evaluating the derivatives. For example, in the third panel of the figure, the second derivatives produced by
Interpolations.jl at the boundary points
x = ±1 are exactly zero. This is probably because the second derivatives at the boundary points were evaluated with the values outside the sampling range, where the interpolator produces a constant value due to the boundary condition
On the other hand,
Dierckx.jl seems to use the interpolator without extrapolated points when evaluating derivatives, even at the boundary points. This means the derivatives at the boundary points are calculated with one-sided limit within the sampling range. I confirm that regardless of the values of
Dierckx.jl maintains this behavior, i.e., the values of the derivatives within the sampling range, including the boundary points, does not change with the boundary condition.
Depending on the situation, it could be desirable or not to incorporate the extrapolated points outside the sampling range when evaluating derivatives at boundary points, so I don’t think one package is superior to the other. In any case, I think the users of these interpolation packages should know the difference before they evaluate derivatives.