Search code examples
pythonscipyspline

Curvature when extrapolated using Natural Cubic Spline

One of the assumptions behind the Natural Cubic Spline is that at the endpoints of the interval of interpolation, the second derivative of the spline polynomials is set to be equal to 0. I tried to show that using the Natural Cubic Spline via from scipy.interpolate import CubicSplines in the example (code below).

from scipy.interpolate import CubicSpline
from numpy import linspace

import matplotlib.pyplot as plt

runge_f = lambda x: 1 / (1 + 25*x**2)

x = linspace(-2, 2, 11)
y = runge_f(x)
cs = CubicSpline(x, y, bc_type = "natural")
t = linspace(-5, 5, 1000)

plt.plot(x, y, "p", color="red")
plt.plot(t, runge_f(t), color="black")
plt.plot(t, cs(t), color="lightblue")

plt.show()

Natural Cubic Spline extrapolation

In the presented example, the extrapolated points' curvature is not equal to zero - shouldn't the extrapolation outside the interval be linear in the Natural Cubic Spline?

Solution

  • The curvature (second derivative) of the spline at the end points is indeed 0, as you can check by running this

    print(cs(x[0],2), cs(x[-1], 2))

    which calculates second derivatives at both ends of your x interpolation interval. However that does not mean the spline is flat beyond the limits -- it continues on as a cubic polynomial. If you want to extrapolate linearly outside the range, you have to do it yourself. It is easy enough to extrapolate flat: replace your cs=.. line with

    from scipy.interpolate import interp1d
cs = interp1d(x,y,fill_value = (y[0],y[-1]), bounds_error = False)

    to get something like this: extrapolate flat

    but a bit more work to extrapolate linearly. I am not sure there is a Python function for that (or rather I am sure there is one I just don't know what it is)