4

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 )

D A
  • 3,130
  • 4
  • 25
  • 41
  • This post is similar - but I don't want to do any numerical evaluation: http://stackoverflow.com/questions/24718313/define-numerical-evaluation-of-a-derivative-of-a-sympy-function – D A Mar 04 '17 at 23:26
  • [`scipy.interpolate.UnivariateSpline`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.UnivariateSpline.html#scipy-interpolate-univariatespline) has methods for derivatives and integrals – Stelios Mar 04 '17 at 23:33
  • @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 A Mar 06 '17 at 18:57
  • I will hack something together using http://docs.sympy.org/dev/modules/functions/elementary.html#piecewise unless someone has a better idea – D A Mar 13 '17 at 18:49
  • Is this what you're looking for? http://docs.sympy.org/latest/modules/functions/special.html#b-splines – asmeurer Apr 12 '17 at 18:27
  • If I am completely honest - when I started this rabbit hole I did not understand what bsplines actually were - and now I have a fuzzy start on them thanks to @FTP below. I do conjecture however, that the "interpolating_spline" solution provides something much more user-friendly without having to dig into the bsplines. – D A Feb 28 '18 at 03:19

1 Answers1

7

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.