使用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)

因此,我可以u(t)通过快速傅里叶变换(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(?).积分几乎是逆操作,但如果有DC组件,它可能会增加一个冲动:?udt <=> U(?) / j? + ?*U(0)*?(?)

要使用采样序列逼近导数,必须按照采样率缩放离散时间角频率ω(范围从0到2π,或-π到π)fs.例如,假设离散时间频率是π/ 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/2DFT分量是实值的,具有0或π弧度的相位,即a*cos(?*fs*t + {0, ?}).给定后者,fs/2组件的导数-a*?*fs*cos(?*fs*t + {?/2, 3*?/2})为,样本时间为0 t == n/fs.

(关于后者,DFT的标准三角插值使用余弦,在这种情况下,导数在采样点处将为零.如果您同时对信号及其导数进行采样,则不一定正确.采样失去相位fs/2相对于信号的分量,但不是相对于它的导数.根据你开始采样的时间,fs/2组分及其导数在样本点可能都不为零.如果幸运的话,其中一个在样本中为0时代,另一个将处于巅峰状态,因为它们是?/2弧形的.)

鉴于fs/2DFT分量始终是实值信号的实数值,当您j在计算导数或积分的过程中将其乘以时,会在结果中引入虚值分量.有一个简单的解决方法.如果N是偶数,则只需将fs/2索引处的组件清零N/2.另一个问题是除以j?积分时除以零.这可以通过向?向量的索引0添加一个小值来解决(例如,finfo(float64).tiny是最小的双精度浮点数).

鉴于? = fs * ?,问题中显示的系统在频域中具有以下形式:

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

这是一个单极点低通滤波器.您派生的解决方案有2个问题.

1. 你是不是缩放频率可变w的fs.
2. 您使用fftfreq,使用范围-0.5到0.5.你需要-π到π.实际上你只需要0到π,因为它i(t)是实值的.在这种情况下,您可以使用rfft和irfft实际值信号,跳过计算负频率.

所有这些,你可能仍然对结果感到失望,因为DFT使用信号的周期性扩展.

例子

首先,这是一个简单的例子,1个正弦曲线(以红色绘制),以1024个样本/秒采样2秒,其衍生物通过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,以及中心和边缘的导数的不连续性.理想情况下,导数应该是方波,但计算完全需要无限带宽.相反,吉布斯在不连续性周围响起:

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.5秒; 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(红色)对从差分方程(蓝色)的计算值.在两种情况下,它们大致相等.我将odeintto 的初始值设置u_dft[0]为显示它返回相同初始值的DFT解决方案.不同之处在于odeint解决方案将继续衰减为零,但DFT解决方案是周期性的,具有2秒的周期.如果更多的i(t)被采样,则DFT结果在这种情况下看起来会更好,因为i(t)从零开始并衰减到零.

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

John P. Boyd讨论了Chebyshev 和Fourier光谱方法中非周期函数的Chebyshev近似方法(在他的密歇根大学页面免费在线阅读).

如果你在信号处理或数学堆栈交换中询问,你可能会得到更多关于这个问题的帮助.