我目前有一个带有时间相关常数的颂歌系统。例如
def fun(u, t, a, b, c):
x = u[0]
y = u[1]
z = u[2]
dx_dt = a * x + y * z
dy_dt = b * (y-z)
dz_dt = -x*y+c*y-z
return [dx_dt, dy_dt, dz_dt]
常量是“a”、“b”和“c”。我目前有每个时间步长的“a”列表，我想在每个时间步长插入，当使用 scipy ode 求解器时……这可能吗？
谢谢！
【问题讨论】：
是的，这是可能的。在a 是常量的情况下，我猜你叫scipy.integrate.odeint(fun, u0, t, args) 其中fun 在你的问题中被定义，u0 = [x0, y0, z0] 是初始条件，t 是要解决的时间点序列ODE 和args = (a, b, c) 是要传递给fun 的额外参数。
在a 取决于时间的情况下，您只需将a 重新考虑为一个函数，例如（给定一个常量a0）：
def a(t):
return a0 * t
然后您将不得不修改 fun，它计算每个时间步的导数，以将先前的更改考虑在内：
def fun(u, t, a, b, c):
x = u[0]
y = u[1]
z = u[2]
dx_dt = a(t) * x + y * z # A change on this line: a -> a(t)
dy_dt = b * (y - z)
dz_dt = - x * y + c * y - z
return [dx_dt, dy_dt, dz_dt]
最后，请注意u0、t 和args 保持不变，您可以再次调用scipy.integrate.odeint(fun, u0, t, args)。
关于这种方法的正确性的一句话。数值积分近似的性能受到影响，我不知道具体如何（没有理论上的保证），但这是一个有效的简单示例：
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
import scipy.integrate
tmax = 10.0
def a(t):
if t < tmax / 2.0:
return ((tmax / 2.0) - t) / (tmax / 2.0)
else:
return 1.0
def func(x, t, a):
return - (x - a(t))
x0 = 0.8
t = np.linspace(0.0, tmax, 1000)
args = (a,)
y = sp.integrate.odeint(func, x0, t, args)
fig = plt.figure()
ax = fig.add_subplot(111)
h1, = ax.plot(t, y)
h2, = ax.plot(t, [a(s) for s in t])
ax.legend([h1, h2], ["y", "a"])
ax.set_xlabel("t")
ax.grid()
plt.show()
我希望这会对你有所帮助。
【讨论】：
不，这在字面意义上是不可能的
“我目前有一个每个时间步的“a”列表，我想在每个时间步插入”
由于求解器具有自适应步长控制，也就是说，它将使用您无法控制的内部时间步长，并且每个时间步长使用函数的多个评估。因此求解器时间步长和数据时间步长之间没有联系。
然而，在给定数据定义分段常数阶跃函数的扩展意义上，有几种方法可以得到解决方案。
您可以使用 ODE 函数和该时间段的常量参数从一个跳转点到另一个跳转点进行积分。之后使用像 concatenate 这样的 numpy 数组操作来组装完整的解决方案。
您可以使用numpy.interp 或scipy.interpolate.interp1d 等插值函数。第一个给出分段线性插值，这里可能不需要。第二个返回一个函数对象，可以配置为“零阶保持”，即分段常数阶跃函数。
您可以实现自己的逻辑，从时间t 到这些参数的正确值。这主要适用于数据具有某种结构的情况，例如，如果它们具有f(int(t/h)) 的形式。
请注意，数值积分的近似阶不仅受限于 RK (solve_ivp) 或多步 (odeint) 方法的阶，还受限于微分方程（部分）的可微阶。如果 ODE 的平滑度远低于方法的阶数，则违反了步长控制机制的隐含假设，这可能导致步长非常小，需要大量的积分步骤。
【讨论】：
scipy.integrate.odeint 在内部使用 Runge-Kutta。我认为它更像是线性多步方法（Adams 方法和后向微分公式方法）等方法的混合，它们比 Runge-Kutta 方法更强大。
numpy.interp 进行立即分段插值，scipy.interpolation.interp1d 返回可以在各种模式下生成的插值函数，包括零阶保持，即，基于提供的时间值数据的分段常数函数。
我也遇到过类似的问题。就我而言，参数a, b, 和c 与时间没有直接关系，而是由当时的x, y, 和z 确定。所以我必须在时间t 得到x, y, z，并计算a, b, c 用于x, y, z 在t+dt 的积分计算。事实证明，如果我更改dt 的值，整个集成结果会发生巨大变化，甚至会变得不合理。
【讨论】：