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 `extrapolation_bc=Flat()`

for `Interpolations.jl`

and `bc="nearest"`

for `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 `extrapolate_bc=Flat()`

.

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 `bc`

used, `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.