High Order derivatives of splines using SciPy

Question

I have created a spline to fit my data in python using:

spline=scipy.interpolate.UnivariateSpline(energy, fpp, k=4)

The equation I want to use involves a summation between n=2 and n=infinity, where n is the order of the differential at a point Eo. However, using;

UnivariateSpline.__call__(spline, e0, nu=n)

to call in the value, I am unable to get a value for anything past the 4th order differential. Is there any other function that people know of for evaluating this function? Above about the 8th order there is a pre-multiplier which should set the value to zero but I still need to go higher than a 4th order.

Answer 1

If you need higher-order derivatives, you won't get good results using equidistant data points. If you can sample your function at arbitrary nodes, I would recommend using Chebyshev points, i.e. $$x_k = \cos\left( \pi \frac{k}{n} \right), \quad k=0\dots n$$ for a polynomial of degree $n$.

You can evaluate the polynomial stably using Barycentric Interpolation. Note that since you're using a polynomial of high degree over the entire interval, you will probably need less data points. Note also that this assumes that your data can be represented by a polynomial, i.e. it is continuous and smooth.

Getting higher-order derivatives from the interpolant is a bit tricky and is ill-conditioned for high derivatives. It can be done, however, using Chebyshev polynomials. In Matlab/Octave (sorry, my Python is not good at all), though, you could do the following:

% We will use the sine function as a test case
f = @(x) sin( 4*pi*x );

% Set the number of points and the interval
N = 40;
a = 0; b = 1;

% Create a Vandermonde-like matrix for the interpolation using the
% three-term recurrence relation for the Chebyshev polynomials.
x = cos( pi*[0:N-1]/(N-1) )';
V = ones( N ); V(:,2) = x;
for k=3:N, V(:,k) = 2*x.*V(:,k-1) - V(:,k-2); end;

% Compute the Chebyshev coefficients of the interpolation. Note that we
% map the points x to the interval [a,b]. Note also that the matrix inverse
% can be either computed explicitly or evaluated using a discrete cosine transform.
c = V \ f( (a+b)/2 + (b-a)/2*x );

% Compute the derivative: this is a bit trickier and relies on the relationship
% between Chebyshev polynomials of the first and second kind.
temp = [ 0 ; 0 ; 2*(N-1:-1:1)'.*c(end:-1:2) ];
cdiff = zeros( N+1 , 1 );
cdiff(1:2:end) = cumsum( temp(1:2:end) );
cdiff(2:2:end) = cumsum( temp(2:2:end) );
cdiff(end) = 0.5*cdiff(end);
cdiff = cdiff(end:-1:3);

% Evaluate the derivative fp at the nodes x. This is useful if you want
% to use Barycentric Interpolation to evaluate it anywhere in the interval.
fp = V(:,1:n-1) * cdiff;

% Evaluate the polynomial and its derivative at a set of points and plot them.
xx = linspace(-1,1,200)';
Vxx = ones( length(xx) , N ); Vxx(:,2) = xx;
for k=3:N, Vxx(:,k) = 2*xx.*Vxx(:,k-1) - Vxx(:,k-2); end;
plot( (a+b)/2 + (b-a)/2*xx , [ Vxx*c , Vxx(:,1:N-1)*cdiff ] );

The code to compute the derivative can be re-applied several times to compute higher derivatives.

If you use Matlab, you may be interested in the Chebfun project, which does most of this automatically and from which parts of the code example above were taken. Chebfun can create an interpolant from literally any function, e.g. continuous, discontinuous, with singularities, etc... and then compute its integral, derivatives, use it to solve ODEs, etc...

Answer 2

Examining the documentation for the UnivariateSpline object, it seems as though you're constructing a 4th order spline, which is why you can't get values for derivatives of greater than 4th order. (All of those derivatives, as you're aware, would be zero.)

SciPy limits the polynomial degree of splines to be 5th order or less (i.e., k <= 5). If you need an 8th order polynomial spline interpolant, you'd have to find an alternate library, or possibly code it up yourself.

Answer 3

By default you use cubic splines, which are 3rd order piecewise polynomials. If you take the fourth derivative of a 3rd order polynomial, you will end up with 0. If you really need these high order differentials, using cubic splines will do you no good.

Answer 4

In Scipy, if you try to compute the n-th order derivative of a k-th order spline, where n > k, then you get a ValueError:

ValueError: 0<=der=5<=k=4 must hold

You could write something like this:

def spline_der(spline, x, nu=0):
    sd = 0
    try:
        sd = spline(x, nu=nu)
    except ValueError:
        pass
    return sd

P.S.: As you can see, instead of using __call__, you can simply write

spline(e0, nu=n)