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

谢谢!

Answer 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()

我希望这能帮到您.