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 )

This post is similar - but I don't want to do any numerical evaluation: stackoverflow.com/questions/24718313/… — D Adams
scipy.interpolate.UnivariateSpline has methods for derivatives and integrals — Stelios
@Stelios -> To provide motivation for a sympy solution: I agree that scipy can do it numerically -> but I want to construct expressions for the derivatives, integrals, and other manipulation of the spline. Each of those operations you do numerically, stack up to slow down execution when you plug numbers in for a final result. If I create 100 splines, take their integrals, add them all up. Now plug in a number. Scipy will be orders of magnitude slower than plugging that number into a an analytical expression generated by sympy. — D Adams
I will hack something together using docs.sympy.org/dev/modules/functions/elementary.html#piecewise unless someone has a better idea — D Adams
Is this what you're looking for? docs.sympy.org/latest/modules/functions/special.html#b-splines — asmeurer

Unknown Unknown · Accepted Answer · 2017-12-19T01:24:26

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.

Construct a symbolic interpolating spline through given points using SymPy

1 Answers

Old answer