Fix pchip endpoints

[egavazzi]

Fix #415.

According to Cleve Moler, Numerical Computing with MATLAB, Chap 3.6, the endpoints need special treatment even when the first/last two slopes don't have the same sign.

Case 2 from the issue (different slopes at the endpoints + extrapolation) ⬇️

Code for the figure

using DataInterpolations
using CondaPkg
using PythonCall
using MATLAB

## Interpolate
u = [0.9, 0.8, 0.86, 0.65, 0.44, 0.76, 0.73, 0.8]
t = [-2, -1, 0.0022, 0.68, 1.41, 2.22, 2.46, 2.76]
t_eval = -4:0.01:3.5

# DataInterpolations.jl
A_datainterpolations = DataInterpolations.PCHIPInterpolation(u, t; extrapolation = ExtrapolationType.Extension)(t_eval)

# SciPy
# CondaPkg.add("scipy")  # to install scipy in your local environment the first time you run the code
pyinterpolate = pyimport("scipy.interpolate")
A_python = pyconvert.(Array,pyinterpolate.PchipInterpolator(t, u)(t_eval))

# Matlab
using MATLAB
@mput t
@mput u
@mput t_eval
mat"A_matlab = interp1(t, u, t_eval,'pchip')"
@mget A_matlab

# Plot
using GLMakie
fig = Figure(size = (1600, 1000))
ax = Axis(fig[1:2, 1]; title = "PCHIP")
lines!(ax, t_eval, A_datainterpolations; linestyle = :dash, linewidth = 3, label = "DataInterpolations.jl")
lines!(ax, t_eval, A_python; linestyle = :solid, linewidth = 3, label = "SciPy")
lines!(ax, t_eval, A_matlab; linestyle = :dashdot, linewidth = 3, label = "MATLAB")
scatter!(ax, t, u, marker = :utriangle, markersize = 15, color = :black, label = "Data")
vlines!(ax, [t[1], t[end]], color = :black, linewidth = 2, linestyle = :dash)
axislegend()

ax = Axis(fig[1, 2]; title = "Difference")
scatterlines!(ax, t_eval, A_python .- A_datainterpolations; linewidth = 3, label = "Scipy - DataInterpolations.jl")
vlines!(ax, [t[1], t[end]], color = :black, linewidth = 2, linestyle = :dash)
axislegend()

ax = Axis(fig[2, 2]; title = "Difference", yticks = [-2e-16, 0, 2e-16])
scatterlines!(ax, t_eval, A_python .- A_matlab; linewidth = 3, label = "Scipy - MATLAB")
vlines!(ax, [t[1], t[end]], color = :black, linewidth = 2, linestyle = :dash)
axislegend()

display(fig)

[SouthEndMusic]

This is probably OK if it matches the reference implementations. Is there some condition that the boundary derivatives satisfy?

[egavazzi]

My understanding is that if the signs of the two nearest slopes to the boundary are the same, then the derivative satisfies d = ((2 * h₁ + h₂) * δ₁ - h₁ * δ₂) / (h₁ + h₂), which is what was already implemented (note: subscript ₁ indicates the nearest slope). In the case where the signs are different, then the derivative should be zero like the interior points (also already implemented), except in the special case abs(d) > 3 * abs(δ₁) for which the derivative should be 3 * δ₁. It is that last special case that was missing and that is being added here.

[SouthEndMusic]

It does look like extrapolation is better behaved 👍