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

Question

我目前有一个具有随时间变化的常量的颂歌系统。例如

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"s 的列表，用于每个时间步，我想在每个时间步插入，当使用 scipy ode 求解器时......这可能吗？

谢谢！

Answer 1

不，字面意义上的那是不可能的

"I currently have a list of "a"s for every time-step which I would like to insert at every time-step"

因为求解器具有自适应步长控制，也就是说，它将使用您无法控制的内部时间步长，并且每个时间步长使用函数的多个评估。因此求解器时间步长和数据时间步长之间没有联系。

然而，在给定数据定义分段常数阶跃函数的扩展意义上，有几种方法可以获得解决方案。

您可以从一个跳跃点积分到另一个跳跃点，使用 ODE 函数和这个时间段的常量参数。之后使用像 concatenate 这样的 numpy 数组操作到 assemble 完整的解决方案。
您可以使用 numpy.interp 或 scipy.interpolate.interp1d 等插值函数。第一个给出分段线性插值，这里可能不需要。第二个returns一个函数对象，可以配置为“零阶保持”，这是一个分段常数阶跃函数。
您可以实现自己的逻辑，从时间 t 到这些参数的正确值。这主要适用于数据具有某种结构的情况，例如，如果它们具有 f(int(t/h)).
的形式

请注意，数值积分的近似阶数不仅受RK（solve_ivp）或多步（odeint）方法阶数的限制，而且受（部分）的可微阶数限制的）微分方程。如果 ODE 远不如方法的阶数平滑，则违反了步长控制机制的隐式假设，这可能导致非常小的步长需要大量积分步骤。

Answer 2

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

希望对您有所帮助。

Answer 3

我也遇到了类似的问题。在我的例子中，参数 a, b, 和 c 不直接随时间变化，而是由当时的 x, y, 和 z 决定的。所以我必须在时间 t 得到 x, y, z，并在 t+dt 计算 x, y, z 的积分计算 a, b, c。事实证明，如果我改变dt值，整个积分结果会发生巨大变化，甚至变得不合理。

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

Solving a system of odes (with changing constant!) using scipy.integrate.odeint?

python

scipy

ode