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

【问题标题】：python - scipy.integrate.odeint returning wrong resultspython - scipy.integrate.odeint 返回错误结果
【发布时间】：2016-02-12 02:45:34
【问题描述】：

我尝试使用 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()

我希望有人能对这里发生的事情有所了解。尝试在 80 到 100（模拟时间）之间更改 T。

【问题讨论】：

标签： python scipy odeint

【解决方案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)

应该会给你想要的结果。

【讨论】：

谢谢，我不知道 odeint 只处理连续函数。你知道是否有其他可用的积分器可以集成不连续的功能？您提出的 Euler 方法在这种情况下有效，但通常是一个糟糕的解决方案。
我认为没有更通用的方法@Miguel。我担心您的问题更多是数学问题而不是编程问题。上面的 get_integral 方法可以修改为梯形规则，这在许多情况下会产生更好的结果，但在方波的情况下会产生错误。我当然不知道任何数学或计算的“通用”方法。