使用 scipy.integrate.odeint 解决颂歌系统（不断变化的常数！）？

Posted 2023-03-12

【中文标题】使用 scipy.integrate.odeint 解决颂歌系统（不断变化的常数！）？

【英文标题】：Solving a system of odes (with changing constant!) using scipy.integrate.odeint?

【发布时间】：2016-01-04 09:23:18

【问题描述】：

我目前有一个带有时间相关常数的颂歌系统。例如

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 求解器时……这可能吗？

谢谢！

【参考方案1】：

是的，这是可能的。在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()

我希望这会对你有所帮助。

【讨论】：

在扭结的情况下，平均结果（对于固定步长的积分器）将具有全局误差 O(dt)，与欧拉方法相同。使用自适应步长，扭结是一个无限曲率的点，将步长减小到救助点（从而增加步数，从而增加浮点错误的累积）。一个扭结不会产生明显的影响。

【参考方案2】：

不，这在字面意义上是不可能的

“我目前有一个每个时间步的“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 方法更强大。

这似乎对应于分段常数的情况。

感谢您的回答，实际上我有 sae 问题，并且我对时间有分段线性依赖性。你提到它可以被插值。您能否提供有关该案例的更多详细信息？

@Ibrahimzawra ：您有函数 numpy.interp 进行立即分段插值，scipy.interpolation.interp1d 返回可以在各种模式下生成的插值函数，包括零阶保持，即，基于提供的时间值数据的分段常数函数。

【参考方案3】：

我也遇到过类似的问题。就我而言，参数a, b, 和c 与时间没有直接关系，而是由当时的x, y, 和z 确定。所以我必须在时间t 得到x, y, z，并计算a, b, c 用于x, y, z 在t+dt 的积分计算。事实证明，如果我更改dt 的值，整个集成结果会发生巨大变化，甚至会变得不合理。

【讨论】：

这没有提供问题的答案。一旦你有足够的reputation，你就可以comment on any post；相反，provide answers that don't require clarification from the asker。 - From Review

或者至少，就@camille 而言，如果这是一个答案，我很难理解它。听起来您更像是遇到了类似的问题，尝试了一些事情，并提供了有关您遇到问题的其他信息。我是不是误会了？如果是这样，编辑您的答案以增加清晰度会很有用。