Derivadas de alta ordem de splines usando SciPy

Eu criei um spline para ajustar meus dados em python usando:

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

A equação que quero usar envolve um somatório entre n = 2 en = infinito, onde n é a ordem do diferencial em um ponto Eo. No entanto, usando;

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

para chamar o valor, não consigo obter um valor para nada além do diferencial de quarta ordem. Existe alguma outra função que as pessoas conheçam para avaliar essa função? Acima da 8ª ordem, existe um pré-multiplicador que deve definir o valor como zero, mas ainda preciso ir além da 4ª ordem.

Se você precisar de derivativos de ordem superior, não obterá bons resultados usando pontos de dados equidistantes. Se você pode testar sua função em nós arbitrários, eu recomendaria o uso de pontos Chebyshev, ou seja, para um polinômio de grau .

Você pode avaliar o polinômio de forma estável usando a interpolação baricêntrica . Observe que, como você está usando um polinômio de alto grau durante todo o intervalo, provavelmente precisará de menos pontos de dados. Observe também que isso pressupõe que seus dados podem ser representados por um polinômio, ou seja, são contínuos e suaves.

Obter derivadas de ordem superior do interpolante é um pouco complicado e está mal condicionado para derivadas altas. No entanto, isso pode ser feito usando os polinômios de Chebyshev . No Matlab / Octave (desculpe, meu Python não é bom), você pode fazer o seguinte:

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

O código para calcular a derivada pode ser reaplicado várias vezes para calcular derivadas mais altas.

Se você usa o Matlab, pode estar interessado no projeto Chebfun , que faz a maior parte disso automaticamente e de quais partes do exemplo de código acima foram tiradas. O Chebfun pode criar um interpolante a partir de literalmente qualquer função, por exemplo, contínua, descontínua, com singularidades, etc.

Examinando a documentação do UnivariateSplineobjeto, parece que você está construindo um spline de 4ª ordem, e é por isso que não é possível obter valores para derivadas maiores que a 4ª ordem. (Todos esses derivativos, como você sabe, seriam zero.)

SciPy limita o grau polinomial de splines para ser de 5ª ordem ou menos (ou seja, k <= 5). Se você precisar de um interpolador de spline polinomial de 8ª ordem, terá que encontrar uma biblioteca alternativa ou possivelmente codificá-la por conta própria.

Por padrão, você usa splines cúbicos, que são polinômios por partes de 3ª ordem. Se você tomar a quarta derivada de um polinômio de 3ª ordem, você terminará com 0. Se realmente precisar desses diferenciais de alta ordem, o uso de splines cúbicos não será bom.

No Scipy, se você tentar calcular a derivada de n-ésima ordem de um spline de k-ésima ordem, em que n> k, obterá ValueError:

ValueError: 0 <= der = 5 <= k = 4 deve manter

Você poderia escrever algo como isto:

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

PS: Como você pode ver, em vez de usar __call__, você pode simplesmente escrever

spline(e0, nu=n)