`numpy.diff` 和 `scipy.fftpack.diff` 在微分时给出不同的结果
`numpy.diff` and `scipy.fftpack.diff` giving different results when differentiating
我正在尝试计算一些数据的导数,并且正在尝试比较有限差分的输出和谱方法的输出。但是结果却大相径庭,我不知道为什么。
考虑下面的示例代码
import numpy as np
from scipy import fftpack as sp
from matplotlib import pyplot as plt
x = np.arange(-100,100,1)
y = np.sin(x)
plt.plot(np.diff(y)/np.diff(x))
plt.plot(sp.diff(y))
plt.show()
这将输出以下结果
橙色输出是 fftpack
输出。没关系,这只是为了举例。
那么,为什么它们如此不同?它们不应该(大致)相同吗?
我很确定可以用 fftpack.diff
的句点关键字校正不同的振幅,但我不知道哪个是正确的句点(我认为应该是 period=1
但这不起作用)。
此外,如何使用 numpy 进行自己的光谱区分?
1 弧度对于差值近似来说是一个相当粗略的步长,你应该希望数据集中的周期数为整数
x = np.arange(-200,200,1)
y = np.sin(np.pi/50*x)
plt.plot(np.diff(y)/np.diff(x))
plt.plot(sp.diff(y,order=1, period=400))
匹配得很好 - 但我不知道 fft 例程
中 period/normalization 的确切合理性
函数scipy.fftpack.diff
计算导数,但它假定输入是周期性的。 period
参数给出了输入序列的周期(即x
区间的总长度)。
在您的情况下,这是 len(x)*dx
,其中 dx = x[1] - x[0]
。
这是一些绘制简单(中心)有限差分(蓝色)和使用 period
参数(红色)的 diff
的结果的代码。变量 x
和 y
与代码中使用的变量相同:
In [115]: plt.plot(0.5*(x[1:]+x[:-1]), np.diff(y)/np.diff(x), 'b')
Out[115]: [<matplotlib.lines.Line2D at 0x1188d01d0>]
In [116]: plt.plot(x, sp.diff(y, period=len(x)*(x[1]-x[0])), 'r')
Out[116]: [<matplotlib.lines.Line2D at 0x1188fc9d0>]
In [117]: plt.xlabel('x')
Out[117]: <matplotlib.text.Text at 0x1157425d0>
请注意,如果您的输入实际上不是周期性的,则 diff
计算的导数在区间末端附近将不准确。
这是另一个示例,使用一个较短的序列,该序列仅包含区间 [0, 1] 中正弦函数的一个完整周期:
In [149]: x = np.linspace(0, 1, 20, endpoint=False)
In [150]: y = np.sin(2*np.pi*x)
In [151]: plt.plot(0.5*(x[1:]+x[:-1]), np.diff(y)/np.diff(x), 'b')
Out[151]: [<matplotlib.lines.Line2D at 0x119872d90>]
In [152]: plt.plot(x, sp.diff(y, period=len(x)*(x[1]-x[0])), 'r')
Out[152]: [<matplotlib.lines.Line2D at 0x119c49090>]
In [153]: plt.xlabel('x')
Out[153]: <matplotlib.text.Text at 0x1197823d0>
我正在尝试计算一些数据的导数,并且正在尝试比较有限差分的输出和谱方法的输出。但是结果却大相径庭,我不知道为什么。
考虑下面的示例代码
import numpy as np
from scipy import fftpack as sp
from matplotlib import pyplot as plt
x = np.arange(-100,100,1)
y = np.sin(x)
plt.plot(np.diff(y)/np.diff(x))
plt.plot(sp.diff(y))
plt.show()
这将输出以下结果
橙色输出是 fftpack
输出。没关系,这只是为了举例。
那么,为什么它们如此不同?它们不应该(大致)相同吗?
我很确定可以用 fftpack.diff
的句点关键字校正不同的振幅,但我不知道哪个是正确的句点(我认为应该是 period=1
但这不起作用)。
此外,如何使用 numpy 进行自己的光谱区分?
1 弧度对于差值近似来说是一个相当粗略的步长,你应该希望数据集中的周期数为整数
x = np.arange(-200,200,1)
y = np.sin(np.pi/50*x)
plt.plot(np.diff(y)/np.diff(x))
plt.plot(sp.diff(y,order=1, period=400))
匹配得很好 - 但我不知道 fft 例程
中 period/normalization 的确切合理性函数scipy.fftpack.diff
计算导数,但它假定输入是周期性的。 period
参数给出了输入序列的周期(即x
区间的总长度)。
在您的情况下,这是 len(x)*dx
,其中 dx = x[1] - x[0]
。
这是一些绘制简单(中心)有限差分(蓝色)和使用 period
参数(红色)的 diff
的结果的代码。变量 x
和 y
与代码中使用的变量相同:
In [115]: plt.plot(0.5*(x[1:]+x[:-1]), np.diff(y)/np.diff(x), 'b')
Out[115]: [<matplotlib.lines.Line2D at 0x1188d01d0>]
In [116]: plt.plot(x, sp.diff(y, period=len(x)*(x[1]-x[0])), 'r')
Out[116]: [<matplotlib.lines.Line2D at 0x1188fc9d0>]
In [117]: plt.xlabel('x')
Out[117]: <matplotlib.text.Text at 0x1157425d0>
请注意,如果您的输入实际上不是周期性的,则 diff
计算的导数在区间末端附近将不准确。
这是另一个示例,使用一个较短的序列,该序列仅包含区间 [0, 1] 中正弦函数的一个完整周期:
In [149]: x = np.linspace(0, 1, 20, endpoint=False)
In [150]: y = np.sin(2*np.pi*x)
In [151]: plt.plot(0.5*(x[1:]+x[:-1]), np.diff(y)/np.diff(x), 'b')
Out[151]: [<matplotlib.lines.Line2D at 0x119872d90>]
In [152]: plt.plot(x, sp.diff(y, period=len(x)*(x[1]-x[0])), 'r')
Out[152]: [<matplotlib.lines.Line2D at 0x119c49090>]
In [153]: plt.xlabel('x')
Out[153]: <matplotlib.text.Text at 0x1197823d0>