Derivados de alto orden de splines usando SciPy

He creado una spline para ajustar mis datos en python usando:

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

La ecuación que quiero usar implica una suma entre n = 2 yn = infinito, donde n es el orden del diferencial en un punto Eo. Sin embargo, usando;

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

para invocar el valor, no puedo obtener un valor para nada más allá del diferencial de cuarto orden. ¿Hay alguna otra función que la gente conozca para evaluar esta función? Arriba, aproximadamente en el octavo orden, hay un pre-multiplicador que debería establecer el valor en cero, pero aún necesito ir más alto que un cuarto orden.

Si necesita derivados de orden superior, no obtendrá buenos resultados utilizando puntos de datos equidistantes. Si puede muestrear su función en nodos arbitrarios, recomendaría usar puntos Chebyshev, es decir,

Puede evaluar el polinomio de forma estable utilizando la interpolación baricéntrica . Tenga en cuenta que dado que está utilizando un polinomio de alto grado en todo el intervalo, probablemente necesitará menos puntos de datos. Tenga en cuenta también que esto supone que sus datos pueden ser representados por un polinomio, es decir, son continuos y suaves.

Obtener derivados de orden superior del interpolante es un poco complicado y está mal condicionado para derivados altos. Sin embargo, se puede hacer usando polinomios de Chebyshev . Sin embargo, en Matlab / Octave (lo siento, mi Python no es bueno en absoluto), puede hacer lo siguiente:

% 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 ] );

El código para calcular la derivada se puede volver a aplicar varias veces para calcular derivadas más altas.

Si usa Matlab, puede estar interesado en el proyecto Chebfun , que hace la mayor parte de esto automáticamente y de qué partes del ejemplo de código anterior se tomaron. Chebfun puede crear un interpolante a partir de literalmente cualquier función, por ejemplo, continua, discontinua, con singularidades, etc. y luego calcular su integral, derivada, usarla para resolver EDO, etc.

Al examinar la documentación del UnivariateSplineobjeto, parece que está construyendo una spline de 4º orden, por lo que no puede obtener valores para derivadas de más de 4º orden. (Todos esos derivados, como saben, serían cero).

SciPy limita el grado polinómico de splines para que sea de quinto orden o menos (es decir, k <= 5). Si necesita un interpolante de spline polinomial de octavo orden, tendría que encontrar una biblioteca alternativa, o posiblemente codificarlo usted mismo.

De manera predeterminada, utiliza splines cúbicas, que son polinomios por partes de tercer orden. Si toma la cuarta derivada de un polinomio de tercer orden, terminará en 0. Si realmente necesita estos diferenciales de alto orden, el uso de splines cúbicos no servirá de nada.

En Scipy, si intenta calcular la derivada de orden n-ésima de una spline de orden k, donde n> k, entonces obtiene un ValueError:

ValueError: 0 <= der = 5 <= k = 4 debe contener

Podrías escribir algo como esto:

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

PD: Como puede ver, en lugar de usar __call__, simplemente puede escribir

spline(e0, nu=n)