“动态" Python中沿轴的N维有限差分 [英] "Dynamic" N-dimensional finite difference in Python along an axis
问题描述
我有一个函数来计算一维np.array的有限差分,我想外推到n维数组.
I have a function to compute the finite difference of a 1d np.array and I want to extrapolate to a n-d array.
函数是这样的:
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)
),但这没关系,因为在对术语进行分组时我会乘以
It is missing the 4th order multiplier (1/(12*h**2)
), but that is ok, because I will multiply when grouping the terms.
我想将其扩展为N维.为此,我将进行以下更改:
I would love to extend it as a N-dimensional. For that I would do the following changes:
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: ]
例如在2D沿第一轴的情况下,这在1D中效果很好,例如,我将不得不更改为以下内容:
This works fine in 1D, if is 2D along first axis for example I would have to change for something like:
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种可能性,并且还在不断发展.如果邻居设置正确,则返回始终有效.
The same applies to 3d, but has 3 possibilities and goes on and on. The return always works if the neighbors are correctly set.
return -fm2 + 16*(fm1+fp1) - 30*fc0 - fp2
当然axis
不能大于len(U.shape)-1
(我称其为尺寸,有什么方法可以代替此代码段吗?
Of course axis
cannot be larger than len(U.shape)-1
(I call this the dimension, is there any way to extract instead this snippet?
如何为这种编码问题做一个优雅而pythonic的方法?
How can I do a elegant and pythonic approach for this coding problem?
有更好的方法吗?
PS:关于np.diff
和np.gradient
,它们不起作用,因为第一个是一阶,第二个是二阶,我正在做四阶近似.实际上,我很快会解决此问题,因此我还将概括该顺序.但是,是的,我希望能够像np.gradient
一样在任何轴上进行操作.
PS: Regarding np.diff
and np.gradient
, those do not work since the first one is first order and the second one is second order, I'm doing a fourth order approximation. In fact soon I finish this problem I will also generalize the order. But yes, I want to be able to do in any axis as np.gradient
do.
推荐答案
一种简单有效的解决方案是在过程的开始和结束时使用swapaxes
:
A simple and effective solution is to use swapaxes
at the very beginning and end of your procedure:
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
这篇关于“动态" Python中沿轴的N维有限差分的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!