Python / SciPy: how to get cubic spline equations from CubicSpline
I am generating a cubic spline plot through a given set of data points:
import matplotlib.pyplot as plt
import numpy as np
import scipy
x = np.array([1, 2, 4, 5]) # sort data points by increasing x value
y = np.array([2, 1, 4, 3])
arr = np.arange(np.amin(x), np.amax(x), 0.01)
s = interpolate.CubicSpline(x, y)
plt.plot(x, y, 'bo', label='Data Point')
plt.plot(arr, s(arr), 'r-', label='Cubic Spline')
plt.legend()
plt.show()
How do I get the spline equations from CubicSpline
? I need equations in the form:
I've tried different methods to get the coefficients, but they all use data that was obtained using different data other than data points.
source to share
From the documentation :
c (ndarray, shape (4, n-1, ...)) Coefficients of polynomials on each segment. The rear dimensions are the same dimensions
y
, excluding the axle. For example, ify
equals 1-d, thenc[k, i]
is the coefficient for(x-x[i])**(3-k)
on the segment betweenx[i]
andx[i+1]
.
So, in your example, the coefficients for the first segment [x 1, x 2] will be in column 0:
- y 1 will be
s.c[3, 0]
- b 1 will be
s.c[2, 0]
- c 1 will be
s.c[1, 0]
- d 1will be
s.c[0, 0]
.
Then for the second segment [x 2, x 3] You are s.c[3, 1]
, s.c[2, 1]
, s.c[1, 1]
and s.c[0, 1]
for y 2, b 2, c 2, d 2etc. etc.
For example:
x = np.array([1, 2, 4, 5]) # sort data points by increasing x value
y = np.array([2, 1, 4, 3])
arr = np.arange(np.amin(x), np.amax(x), 0.01)
s = interpolate.CubicSpline(x, y)
fig, ax = plt.subplots(1, 1)
ax.hold(True)
ax.plot(x, y, 'bo', label='Data Point')
ax.plot(arr, s(arr), 'k-', label='Cubic Spline', lw=1)
for i in range(x.shape[0] - 1):
segment_x = np.linspace(x[i], x[i + 1], 100)
# A (4, 100) array, where the rows contain (x-x[i])**3, (x-x[i])**2 etc.
exp_x = (segment_x - x[i])[None, :] ** np.arange(4)[::-1, None]
# Sum over the rows of exp_x weighted by coefficients in the ith column of s.c
segment_y = s.c[:, i].dot(exp_x)
ax.plot(segment_x, segment_y, label='Segment {}'.format(i), ls='--', lw=3)
ax.legend()
plt.show()
source to share
Edit: Since the coefficients are available as an attribute (see @ali_m's answer), the approach described here is unnecessarily indirect. I leave it online in case anyone is using a more opaque library ever stumbling over this question.
One way is to evaluate in the node we are interested in:
coeffs = [s(2)] + [s.derivative(i)(2) / misc.factorial(i) for i in range(1,4)]
s(2.5)
# -> array(1.59375)
sum(coeffs[i]*(2.5-2)**i for i in range(4))
# -> 1.59375
Strictly speaking, higher derivatives don't exist in nodes, but scipy seems to return the correct one-way derivative, so it works, although it shouldn't.
source to share