"Dynamische" N-dimensionale endliche Differenz in Python entlang einer Achse

Question

Ich habe eine Funktion, um die endliche Differenz eines 1d np.array zu berechnen, und ich möchte auf ein n-d-Array extrapolieren.

Die Funktion ist wie folgt:

def fpp_fourth_order_term(U): 
    """Returns the second derivative of fourth order term without the interval multiplier.""" 
    # U-slices 
    fm2 = values[ :-4] 
    fm1 = values[1:-3] 
    fc0 = values[2:-2] 
    fp1 = values[3:-1] 
    fp2 = values[4: ] 

    return -fm2 + 16*(fm1+fp1) - 30*fc0 - fp2 

Es ist die 4. Ordnung Multiplikator (1/(12*h**2)) fehlt, aber das ist ok, weil ich multiplizieren, wenn die Bedingungen zu gruppieren.

Ich würde es gerne als N-dimensional erweitern. Denn das würde ich folgende Änderungen tun:

def fpp_fourth_order_term(U, axis=0): 
    """Returns the second derivative of fourth order term along an axis without the interval multiplier.""" 
    # U-slices 

Aber hier ist das Problem

fm2 = values[ :-4] 
    fm1 = values[1:-3] 
    fc0 = values[2:-2] 
    fp1 = values[3:-1] 
    fp2 = values[4: ] 

Dieses in 1D funktioniert gut, wenn ist 2D entlang der ersten Achse zum Beispiel würde ich ändern müssen so etwas wie:

fm2 = values[:-4,:] 
    fm1 = values[1:-3,:] 
    fc0 = values[2:-2,:] 
    fp1 = values[3:-1,:] 
    fp2 = values[4:,:] 

Aber entlang der zweite Achse wäre:

fm2 = values[:,:-4] 
    fm1 = values[:,1:-3] 
    fc0 = values[:,2:-2] 
    fp1 = values[:,3:-1] 
    fp2 = values[:,4:] 

Das Gleiche gilt für 3D, hat aber 3 Möglichkeiten und geht weiter und weiter. Die Rückgabe funktioniert immer, wenn die Nachbarn richtig eingestellt sind.

return -fm2 + 16*(fm1+fp1) - 30*fc0 - fp2 

Natürlich axis kann nicht größer sein als len(U.shape)-1 (I nenne dies die Dimension, ist es eine Möglichkeit, statt diese Schnipsel zu extrahieren?

Wie kann ich einen eleganten und pythonic Ansatz für diese Codierung Problem zu tun?

gibt es einen besseren Weg, es zu tun

PS: In Bezug auf np.diff und np.gradient, diejenigen nicht funktionieren, da die erste erster Ordnung und die zweite i s zweite Ordnung, ich mache eine Näherung vierter Ordnung. In der Tat beende ich bald dieses Problem, ich werde auch den Auftrag verallgemeinern. Aber ja, ich möchte in der Lage sein, in jeder Achse wie np.gradient zu tun.

Answer 1

Eine einfache und effektive Lösung ist swapaxes ganz am Anfang und am Ende der Prozedur zu verwenden:

import numpy as np 

def f(values, axis=-1): 
    values = values.swapaxes(0, axis) 

    fm2 = values[ :-4] 
    fm1 = values[1:-3] 
    fc0 = values[2:-2] 
    fp1 = values[3:-1] 
    fp2 = values[4: ] 

    return (-fm2 + 16*(fm1+fp1) - 30*fc0 - fp2).swapaxes(0, axis) 

a = (np.arange(4*7*8)**3).reshape(4,7,8) 
res = f(a, axis=1) 
print(res) 
print(res.flags) 

Ausgang:

# [[[ 73728 78336 82944 87552 92160 96768 101376 105984] 
# [110592 115200 119808 124416 129024 133632 138240 142848] 
# [147456 152064 156672 161280 165888 170496 175104 179712]] 

# [[331776 336384 340992 345600 350208 354816 359424 364032] 
# [368640 373248 377856 382464 387072 391680 396288 400896] 
# [405504 410112 414720 419328 423936 428544 433152 437760]] 

# [[589824 594432 599040 603648 608256 612864 617472 622080] 
# [626688 631296 635904 640512 645120 649728 654336 658944] 
# [663552 668160 672768 677376 681984 686592 691200 695808]] 

# [[847872 852480 857088 861696 866304 870912 875520 880128] 
# [884736 889344 893952 898560 903168 907776 912384 916992] 
# [921600 926208 930816 935424 940032 944640 949248 953856]]] 

Das Ergebnis ist auch angrenzend.

# C_CONTIGUOUS : True 
# F_CONTIGUOUS : False 
# OWNDATA : False 
# WRITEABLE : True 
# ALIGNED : True 
# UPDATEIFCOPY : False