python - Construct a symbolic interpolating spline through given points using SymPy

Question

Pretend I start with some simple dataset which is defined on R2 follows:

DataPointsDomain = [0,1,2,3,4,5]
DataPointsRange =  [3,6,5,7,9,1]

With scipy I can make a lazy polynomial spline using the following:

ScipySplineObject = scipy.interpolate.InterpolatedUnivariateSpline( 
    DataPointsDomain, 
    DataPointsRange, 
    k = 1, )

What is the equivalent object in sympy??

SympySplineObject = ...???

(I want to define this object and do analytic sympy manipulation like taking integrals, derivatives, etc... on the sympy object )

Answer 1

In SymPy versions above 1.1.1, including the current development version, there is a built-in method interpolating_spline which takes four arguments: the spline degree, the variable, domain values and range values.

from sympy import *
DataPointsDomain = [0,1,2,3,4,5]
DataPointsRange =  [3,6,5,7,9,1]
x = symbols('x')
s = interpolating_spline(3, x, DataPointsDomain, DataPointsRange)

This returns

Piecewise((23*x**3/15 - 33*x**2/5 + 121*x/15 + 3, (x >= 0) & (x <= 2)),
   (-2*x**3/3 + 33*x**2/5 - 55*x/3 + 103/5, (x >= 2) & (x <= 3)), 
   (-28*x**3/15 + 87*x**2/5 - 761*x/15 + 53, (x >= 3) & (x <= 5)))

which is a "not a knot" cubic spline through the given points.

---

Old answer

An interpolating spline can be constructed with SymPy, but this takes some effort. The method bspline_basis_set returns the basis of B-splines for given x-values, but then it's up to you to find their coefficients.

First, we need the list of knots, which is not exactly the same as the list of x-values (xv below). The endpoints xv[0] and xv[-1] will appear deg+1 times where deg is the degree of the spline, because at the endpoints all the coefficients change values (from something to zero). Also, some of the x-values close to them may not appear at all, as there will be no changes of coefficients there ("not a knot" conditions). Finally, for even-degree splines (yuck) the interior knots are placed midway between data points. So we need this helper function:

from sympy import *
def knots(xv, deg):
    if deg % 2 == 1:
        j = (deg+1) // 2
        interior_knots = xv[j:-j]
    else:
        j = deg // 2
        interior_knots = [Rational(a+b, 2) for a, b in zip(xv[j:-j-1], xv[j+1:-j])]
    return [xv[0]] * (deg+1) + interior_knots + [xv[-1]] * (deg+1)

After getting b-splines from bspline_basis_set method, one has to plug in the x-values and form a linear system from which to find the coefficients coeff. At last, the spline is constructed:

xv = [0, 1, 2, 3, 4, 5]
yv =  [3, 6, 5, 7, 9, 1]
deg = 3
x = Symbol("x")
basis = bspline_basis_set(deg, knots(xv, deg), x)
A = [[b.subs(x, v) for b in basis] for v in xv]
coeff = linsolve((Matrix(A), Matrix(yv)), symbols('c0:{}'.format(len(xv))))
spline = sum([c*b for c, b in zip(list(coeff)[0], basis)])
print(spline)

This spline is a SymPy object. Here it is for degree 3:

3*Piecewise((-x**3/8 + 3*x**2/4 - 3*x/2 + 1, (x >= 0) & (x <= 2)), (0, True)) + Piecewise((x**3/8 - 9*x**2/8 + 27*x/8 - 27/8, (x >= 3) & (x <= 5)), (0, True)) + 377*Piecewise((19*x**3/72 - 5*x**2/4 + 3*x/2, (x >= 0) & (x <= 2)), (-x**3/9 + x**2 - 3*x + 3, (x >= 2) & (x <= 3)), (0, True))/45 + 547*Piecewise((x**3/9 - 2*x**2/3 + 4*x/3 - 8/9, (x >= 2) & (x <= 3)), (-19*x**3/72 + 65*x**2/24 - 211*x/24 + 665/72, (x >= 3) & (x <= 5)), (0, True))/45 + 346*Piecewise((x**3/30, (x >= 0) & (x <= 2)), (-11*x**3/45 + 5*x**2/3 - 10*x/3 + 20/9, (x >= 2) & (x <= 3)), (31*x**3/180 - 25*x**2/12 + 95*x/12 - 325/36, (x >= 3) & (x <= 5)), (0, True))/45 + 146*Piecewise((-31*x**3/180 + x**2/2, (x >= 0) & (x <= 2)), (11*x**3/45 - 2*x**2 + 5*x - 10/3, (x >= 2) & (x <= 3)), (-x**3/30 + x**2/2 - 5*x/2 + 25/6, (x >= 3) & (x <= 5)), (0, True))/45

You can differentiate it, with

spline.diff(x)

You can integrate it:

integrate(spline, (x, 0, 5))   #  197/3

You can plot it and see it indeed interpolates the given values:

plot(spline, (x, 0, 5))

I even plotted them for degrees 1,2,3 together:

Disclaimers:

• The code given above works in the development version of SymPy, and should work in 1.1.2+; there was a bug in B-spline method in previous versions.
• Some of this takes a good deal of time because Piecewise objects are slow. In my experience, the basis construction takes longest.