python - 使用scipy fft和ifft数值求解常微分方程

Question

我在时域中有一个ordinary differential equation，如下所示:

C*du/dt = -g*u + I

其中 I = A*t/tau*exp^(1-t/tau)
在freq域中:

u(w) = I(w)/(g*(1+C/g*j*w))

j是复数 sqrt(-1)
因此，我可以通过使用 fast Fourier transform( u(t))进入freq域，然后再使用 fft来获取 ifft。

代码:

t = np.linspace(0.,499.9,5000)
I = q*t*np.exp(1-t/Tau_ca)/Tau_ca
u1 = np.fft.ifft(np.fft.fft(I)/(G_l*(1.+1.j*(C_m/G_l)*np.fft.fftfreq(t.shape[-1]))))

但是，当我将此获得的 u(t)与其他方法(例如微分方程的数值积分或其解析形式)进行比较时，这是不正确的。我已经尝试过并且未能弄清楚我的错误在哪里。

请指教。

Answer 1

正弦波或复指数的导数与频率成正比，并通过π/2相移。对于复杂的指数，相移等效于j的乘积。例如，d/dt exp(j*Ω*t) == j*Ω * exp(j*Ω*t) == Ω * exp(j*π/2) * exp(j*Ω*t) == Ω * exp(j*(Ω*t + π/2))。因此，如果您具有傅立叶变换对u(t) <=> U(Ω)，则为du/dt <=> jΩ * U(Ω)。积分几乎是逆运算，但如果有直流分量，则可能会增加一个脉冲:∫udt <=> U(Ω) / jΩ + π*U(0)*δ(Ω)
为了使用采样序列近似导数，您必须按采样率fs缩放离散时间角频率ω(范围为0到2π，或-π到π)。例如，假设离散时间频率为π/2弧度/样本，例如序列[1, 0, -1, 0, ...]。在原始信号中，这对应于fs/4。派生词是d/dt cos(2*π*fs/4*t) == d/dt cos(π/2*fs*t) == -π/2*fs*sin(π/2*fs*t) == π/2*fs*cos(π/2*fs*t + π/2)。

您必须以大于信号带宽两倍的fs进行采样。确切的fs/2处的采样组件不可靠。具体而言，每个周期只有2个样本，fs/2分量的振幅会交替显示序列中第一个样本的符号。因此，对于实值信号，fs/2 DFT分量是实数值，其相位为0或π弧度，即a*cos(π*fs*t + {0, π})。给定后者，fs/2分量的导数为-a*π*fs*cos(π*fs*t + {π/2, 3*π/2})，对于采样时间t == n/fs，该值为0。

(关于后者，DFT的标准三角插值使用余弦，在这种情况下，导数在采样点处为零。如果同时对信号及其导数进行采样，则不一定是正确的。采样会丢失相位相对于信号，而不是相对于其导数，fs/2分量的比率。取决于您开始采样的时间，fs/2分量及其导数在采样点上可能都非零。采样时间，另一个将达到峰值，因为它们的π/2弧度分开。)

鉴于fs/2 DFT分量始终是实值信号的实值，因此在计算导数或积分的过程中将其乘以j时，这会在结果中引入虚值分量。有一个简单的解决方法。如果N是偶数，则将索引fs/2的N/2组件置零。另一个问题是用jΩ除以进行积分时除以零。这可以通过向Ω向量的索引0添加一个较小的值来解决(例如finfo(float64).tiny是最小的 double 浮点数)。

给定Ω = fs * ω，问题中显示的系统在频域中具有以下形式:

H(Ω) = 1 / (g + j*Ω*C)
U(Ω) = I(Ω) * H(Ω)

这是一个单极点低通滤波器。您得出的解决方案有2个问题。

您不会通过w缩放频率变量fs。

您使用fftfreq，它的范围是-0.5到0.5。您需要-π到π。
实际上，您只需要0到π，因为i(t)是实值。在这个
在这种情况下，您可以将rfft和irfft用于实值信号，
跳过计算负频率。

如此说来，您可能仍然会对结果感到失望，因为DFT使用信号的周期性扩展。

例子

首先，这是一个简单的示例，它以1024个样本/秒的速度采样2秒钟的1 Hz正弦波(以红色绘制)，并通过DFT计算出其导数(以蓝色绘制):

from pylab import *

fs = 1024
t = arange(2*fs, dtype=float64) / fs
N = len(t) // 2 + 1    # non-negative frequencies
w = 2 * pi / len(t) * arange(N)
Omega = w * fs

x0 = cos(2*pi*t)    # w0 == 2*pi/fs
X0 = rfft(x0);
# Zero the fs/2 component. It's zero here anyway.
X0[-1] = 0  
dx0 = irfft(X0 * 1j*Omega)
plot(t, x0, 'r', t, dx0, 'b')
show()

这是一个简单的情况-具有有限带宽的周期性信号。只需确保以足够高的速率对整数周期进行采样即可避免混叠。

下一个示例是一个三角波，其斜率为1和-1，并且导数在中心和边缘处不连续。理想情况下，导数应为方波，但完美地进行计算将需要无限带宽。而是在不连续性周围有 Gibbs ringing:

t2 = t[:fs]
m = len(t) // (2*len(t2))
x1 = hstack([t2, 1.0 - t2] * m)
X1 = rfft(x1)
X1[-1] = 0
dx1 = irfft(X1 * 1j*Omega)
plot(t, x1, 'r', t, dx1, 'b')
show()

如果您要求解非周期性系统，则DFT的隐式周期性扩展会出现问题。这是同时使用 odeint和DFT的相关系统的解决方案( tau设置为0.5s； g和 C设置为1 Hz转折频率):

from scipy.integrate import odeint

a = 1.0; tau = 0.5
g = 1.0; C = 1.0 / (2 * pi)
i = lambda t: a * t / tau * exp(1 - t / tau)
f = lambda u, t: (-g*u + i(t)) / C

H = 1 / (g + 1j*Omega*C)  # system function
I = rfft(i(t))
I[-1] = 0
U_DFT = I * H
u_dft = irfft(U_DFT)
u_ode = odeint(f, u_dft[0], t)[:,0]

td = t[:-1] + 0.5/fs
subplot('221'); title('odeint u(t)');
plot(t, u_ode)
subplot('222'); title('DFT u(t)');
plot(t, u_dft)
subplot('223'); title('odeint du/dt')
plot(td, diff(u_ode)*fs, 'r',
     t, (-g*u_ode + i(t)) / C, 'b')           
subplot('224'); title('DFT du/dt')
plot(td, diff(u_dft)*fs, 'r',
     t, (-g*u_dft + i(t)) / C, 'b')
show()

du/dt图覆盖了由 diff估计的导数(红色)与根据微分方程计算得到的值(蓝色)的重叠。它们在两种情况下大致相等。我将 odeint的初始值设置为 u_dft[0]，以显示它为相同的初始值返回DFT解决方案。不同之处在于 odeint解决方案将继续衰减为零，但DFT解决方案是周期性的，周期为2s。在这种情况下，如果对i(t)进行更多采样，DFT结果将看起来更好，因为i(t)从零开始衰减到零。

DFT使用零填充来执行线性卷积。特别是在这种情况下，输入的零填充将有助于将零状态响应的 transient 与其稳态隔离开。但是，更常见的是使用Laplace或z转换系统函数来分析ZSR/ZIR。 scipy.signal 具有分析LTI系统的工具，包括以多项式，零极点增益和状态空间形式映射连续时间到离散时间。

约翰·博伊德(John P. Boyd)讨论了 Chebyshev and Fourier Spectral Methods中非周期性函数的Chebyshev逼近方法(可在他的密歇根大学页面免费在线获得)。

如果您在 Signal Processing或 Mathematics堆栈交换上询问，您可能会在诸如此类的问题上获得更多帮助。