Python中沿轴的“动态”N维有限差分答案

【问题标题】："Dynamic" N-dimensional finite difference in Python along an axisPython中沿轴的“动态”N维有限差分
【发布时间】：2017-12-18 02:20:36
【问题描述】：

我有一个计算一维 np.array 的有限差分的函数，我想推断一个 n-d 数组。

函数是这样的：

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

它缺少 4 阶乘数 (1/(12*h**2))，但这没关系，因为我会在对术语进行分组时进行乘数。

我很想将其扩展为 N 维。为此，我将进行以下更改：

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

问题来了

    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（我把这个叫做维度，有没有办法提取这个sn-p？

我怎样才能为这个编码问题采取一种优雅和 Pythonic 的方法？

有没有更好的方法？

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

【问题讨论】：

优雅可能要求太多，但these links 可能会帮助您入门。
slice 模式不支持 fancy 索引作为U[:,1:-3] 的另一种方式（np.take）模糊支持，我需要建立范围（因为range(1,-1) 没有意义），np.take 也复制数据（我想避免。但谢谢。我正在进步。
fc0 = values.take(indices=range(2, values.shape[axis]-2), axis=axis) 之类的东西可以工作，但是复制数组 5 次会让我的记忆大打折扣。这些数组会变大。 :-/
实际上，简单有效的方法是在输入上使用swapaxes 将轴移动到最左边进行操作，然后在最终结果上再次使用swapaxes。
我已阅读 swapaxes，但不知道如何处理。

标签： python arrays numpy differential-equations

【解决方案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)

输出：

# [[[ 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]]]

结果甚至是连续的。

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

【讨论】：