python中的SDE求解器:手动确定积分器步长(dt)

计算科学 Python 随机颂歌
2021-12-01 18:50:44

目标:我正在尝试解决 SDE 系统,同时使用SDEintpython 3.x 中的包。它是由Munz 于 2009 年出版的 Zombie Apocalypse 模型改编并受到启发的 SDE 系统。我试图在代码中添加对术语的解释。

问题:我想解决的问题出现在下面的示例中(在 double_scalar 中出现溢出)。请注意,噪声项当前为 0,因此它基本上是 ODE 系统,而不是 SDE 系统。

import matplotlib.pyplot as plt
import numpy as np
import sdeint
from scipy.integrate import odeint

p, f, e, k = 0.7, 0.0002, 0.05, 0.02    # human birth rate, confilt occurence rate, conflict end rate, zombie killing rate

tspan = np.linspace(0, 20., 100)
y0 = np.array([100000., 0., 0., p])

def ff(y, t):
    Hi = y[0]
    Ci = y[1]
    Zi = y[2]

    # Human are created (exponential)                       --> y[3] * Hi
    # Humans and zombies engage in conflict at rate (f)     --> - f * Hi * Zi
    # the conficts end at rate (e) without a victor         --> + e * Ci
    f0 = y[3] * Hi - f * Hi * Zi + e * Ci

    # Humans and zombies engage in conflict at rate (f)     --> f * Hi * Zi
    # the conficts end at rate (e) without a victor         --> - e * Ci
    # and zombies kill humans in conflict at rate (k)       --> - k * Ci
    f1 = f * Hi * Zi - e * Ci - k * Ci

    # Zombies arise from the earth at a fixed rate          --> 10**7
    # they arise as victor from a conflict                  --> + k * Ci
    # enter conflicts                                       --> - f * Hi * Zi
    # and exit conflicts without killing                    --> + e * Ci
    f2 = 10**7 + k * Ci - f * Hi * Zi + e * Ci

    f3 = 0
    print(f0, f1 ,f2, f3)
    return np.array([f0, f1, f2, f3])

def GG(y, t):
    return np.diag([0, 0, 0, 0])

result = sdeint.itoint(ff, GG, y0, tspan)

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.set_yscale('log')
ax.set_ylim(bottom=1, top=10**13)
plt.plot(result)
plt.show()

我曾尝试使用 Scipy 解决相同的 ODE 系统(无噪声)odeint,由于可以向odeint.

result = odeint(ff, y0, tspan, mxstep=10 ** 9, rtol=10 ** (-3),
                atol=10 ** (-3), hmax=1)

所以我认为可能SDEint使用了太大的集成步骤,但是不可能手动添加额外的参数到sdeint.

潜在解决方案(?):所以我想知道,如果一个解决方案可能是手动设置集成的集合大小。np.linspace我试图通过增加例如步骤的数量来做到这一点

tspan = np.linspace(0, 20., 100000) # 100 --> 100000

但是,这似乎无法正常工作。

问题:由于我对 ODE/SDE 不是很熟悉,我不确定我是否正确地解决了这个问题。这是解决集成步骤过大问题的常用方法吗?解决方案可能是增加步长还是我完全弄错了?

任何帮助将不胜感激!

0个回答
没有发现任何回复~