3

我正在尝试计算一些数据的导数，并且我正在尝试比较有限差分的输出和光谱方法的输出。但是结果非常不同，我无法弄清楚为什么。

考虑下面的示例代码

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 关键字来校正不同的幅度，但我无法确定哪个是正确的周期（我认为它应该是period=1，但那不起作用）。

此外，如何使用 numpy 进行自己的光谱区分？

4

2 回答 2

5

该函数scipy.fftpack.diff计算导数，但它假设输入是周期性的。参数给出输入序列的period周期（即x间隔的总长度）。

就您而言，这len(x)*dx就是dx = x[1] - x[0].

这是一些绘制简单（居中）有限差分（蓝色）和diff使用period参数的结果（红色）的代码。变量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>

情节2

于 2017-02-13T07:23:13.500 回答
0

1 rad 对于差异近似来说是一个相当粗略的步幅，您应该希望数据集中有整数个周期

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 例程中周期/标准化的确切原因

于 2017-02-13T07:16:24.880 回答