Derivate höherer Ordnung von Splines mit SciPy

Ich habe einen Spline erstellt, um meine Daten in Python anzupassen, indem ich:

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

Die Gleichung, die ich verwenden möchte, beinhaltet eine Summation zwischen n = 2 und n = unendlich, wobei n die Ordnung des Differentials an einem Punkt Eo ist. Verwenden Sie jedoch;

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

Um den Wert aufzurufen, kann ich keinen Wert für etwas erhalten, das über das Differential 4. Ordnung hinausgeht. Gibt es eine andere Funktion, die Menschen zur Bewertung dieser Funktion kennen? Oberhalb der 8. Ordnung gibt es einen Vorvervielfacher, der den Wert auf Null setzen sollte, aber ich muss immer noch höher als eine 4. Ordnung gehen.

Wenn Sie Derivate höherer Ordnung benötigen, erzielen Sie mit äquidistanten Datenpunkten keine guten Ergebnisse. Wenn Sie Ihre Funktion an beliebigen Knoten abtasten können, würde ich die Verwendung von Chebyshev-Punkten empfehlen, dh für ein Polynom Grad .

Sie können das Polynom mithilfe der baryzentrischen Interpolation stabil auswerten . Beachten Sie, dass Sie wahrscheinlich weniger Datenpunkte benötigen, da Sie über das gesamte Intervall ein Polynom von hohem Grad verwenden. Beachten Sie auch, dass dies voraussetzt, dass Ihre Daten durch ein Polynom dargestellt werden können, dh sie sind kontinuierlich und glatt.

Ableitungen höherer Ordnung aus dem Interpolanten zu erhalten, ist etwas schwierig und für Ableitungen hoher Ordnung schlecht konditioniert. Dies kann jedoch unter Verwendung von Chebyshev-Polynomen erfolgen . In Matlab / Octave (sorry, mein Python ist überhaupt nicht gut) können Sie jedoch Folgendes tun:

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

Der Code zum Berechnen der Ableitung kann mehrmals erneut angewendet werden, um höhere Ableitungen zu berechnen.

Wenn Sie Matlab verwenden, könnte Sie das Chebfun- Projekt interessieren , das das meiste automatisch erledigt und aus dem Teile des obigen Codebeispiels stammen. Chebfun kann einen Interpolanten aus buchstäblich jeder Funktion erstellen, z. B. stetig, diskontinuierlich, mit Singularitäten usw. und dann sein Integral und seine Ableitungen berechnen, um ODEs usw. zu lösen.

Wenn Sie die Dokumentation für das UnivariateSplineObjekt untersuchen, scheint es, als würden Sie einen Spline 4. Ordnung konstruieren, weshalb Sie keine Werte für Ableitungen größer als 4. Ordnung erhalten können. (Alle diese Derivate wären, wie Sie wissen, Null.)

SciPy begrenzt den Polynomgrad von Splines auf 5. Ordnung oder weniger (dh k <= 5). Wenn Sie einen Polynom-Spline-Interpolanten 8. Ordnung benötigen, müssen Sie eine alternative Bibliothek finden oder diese möglicherweise selbst codieren.

Standardmäßig verwenden Sie kubische Splines, die stückweise Polynome 3. Ordnung sind. Wenn Sie die vierte Ableitung eines Polynoms 3. Ordnung nehmen, erhalten Sie 0. Wenn Sie diese Differentiale hoher Ordnung wirklich benötigen, hilft es Ihnen nichts, kubische Splines zu verwenden.

Wenn Sie in Scipy versuchen, die Ableitung n-ter Ordnung eines Splines k-ter Ordnung zu berechnen, wobei n> k ist, erhalten Sie a ValueError:

ValueError: 0 <= der = 5 <= k = 4 muss gelten

Sie könnten so etwas schreiben:

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

PS: Wie Sie sehen, können Sie anstelle von __call__einfach schreiben

spline(e0, nu=n)