python - scipy.integrate.odeint 返回错误结果

Question

我尝试使用 python 3.5 和 scipy.integrate.odeint 函数积分方波，但结果没有任何意义，并且随着所选时间点数组的不同而变化很大。

方波的周期为 10 秒，模拟运行 100 秒。由于时间点数组的大小为 500，因此方波的每个周期上将有 50 个时间点，但这似乎没有发生。使用可选参数 hmax=0.02 可以修复该问题，但不应该自动推断吗？

代码如下:

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate

# dx/dt = f(t),  where f(t) is a square wave
def f(x, t):
    return float(t % 10.0 < 5.0) * 0.3

T = 100
tt = np.linspace(0, T, 500)
xx = integrate.odeint(f, 0, tt, hmax=0.2)

plt.figure()
plt.subplot(2,1,1)
plt.plot(tt, xx)
plt.axis([0,T,0,16])

plt.subplot(2,1,2)
plt.plot(tt, [f(None,t) for t in tt])
plt.axis([0, T, 0, 1])
plt.show()

我希望有人能够解释一下这里发生的事情。尝试将 T 更改为 80 到 100(模拟时间)。

Answer 1

我认为你的问题是 odeint 函数采用连续常微分方程，而方波则不是。

我首先将方波函数重新定义为:

def g(t):
    return float(t % 10.0 < 5.0) * 0.3

然后定义一个函数来逐步计算积分:

def get_integral(tt):
    intarray = np.zeros_like(tt)
    step_size = tt[1] -tt[0]
    for i,t in enumerate(tt):
        intarray[i] = intarray[i-1] + g(t)*step_size
    return intarray

然后:

xx = get_integral(tt)

应该给你你正在寻找的结果。