Python中沿轴的“动态”nDimension有限差分

Question

我有一个计算一维有限差分的函数np.数组我想外推到n维数组。在

函数如下：

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

它缺少四阶乘法器（1/(12*h**2)），但这没问题，因为我将在分组项时乘法。在

我想把它扩展成N维。为此，我将做以下更改：

^{pr2}$

但问题是

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

这在1D中工作得很好，例如，如果是沿第一个轴的2D，我将不得不改变如下内容：

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

但是沿着第二个轴是：

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

同样的道理也适用于3d，但有3种可能，而且会一直持续下去。如果邻域设置正确，则返回始终有效。在

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

当然axis不能大于len(U.shape)-1（我称之为维度，有没有任何方法来代替这个片段的提取？在

我如何为这个编码问题提供一个优雅的python方法？在

有更好的方法吗？在

注：关于np.diff和np.gradient，它们不起作用，因为第一个是一阶，第二个是二阶，我做的是四阶近似。事实上，我很快就完成了这个问题，我也将推广这个顺序。但是是的，我希望能够像np.gradient那样在任何轴上执行。在

Answer 1

一个简单有效的解决方案是在程序的开始和结束时使用swapaxes：

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)

输出：

^{pr2}$

结果甚至是连续的。在

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