Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Construct a symbolic interpolating spline through given points using SymPy

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 )

like image 286
D Adams Avatar asked Sep 06 '25 09:09

D Adams


1 Answers

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:

splines

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.

Donate For Us

If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!