如何在 python 中实现随机 ODE (SDE) 系统？答案

【问题标题】：How to implement a system of stochastic ODEs (SDEs) in python?如何在 python 中实现随机 ODE (SDE) 系统？
【发布时间】：2019-06-29 03:27:42
【问题描述】：

我有一个 ODE 系统，我试图在其中包含一个“错误”术语，以便它成为一个随机 ODE 系统。

为了解决 Python 中的 ODE 系统，我通常使用 scipy 的 odeint。

一个源自Scipy Cookbook的例子，涉及著名的僵尸启示录：

# zombie apocalypse modeling
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
plt.rcParams['figure.figsize'] = 10, 8

P = 0      # birth rate
d = 0.0001  # natural death percent (per day)
B = 0.0095  # transmission percent  (per day)
G = 0.0001  # resurect percent (per day)
A = 0.0001  # destroy percent  (per day)

# solve the system dy/dt = f(y, t)
def f(y, t):
     Si = y[0]
     Zi = y[1]
     Ri = y[2]
     # the model equations (see Munz et al. 2009)
     f0 = P - B*Si*Zi - d*Si
     f1 = B*Si*Zi + G*Ri - A*Si*Zi
     f2 = d*Si + A*Si*Zi - G*Ri
     return [f0, f1, f2]

# initial conditions
S0 = 500.              # initial population
Z0 = 0                 # initial zombie population
R0 = 0                 # initial death population
y0 = [S0, Z0, R0]     # initial condition vector
t  = np.linspace(0, 5., 1000)         # time grid

# solve the DEs
soln = odeint(f, y0, t)
S = soln[:, 0]
Z = soln[:, 1]
R = soln[:, 2]

# plot results
plt.figure()
plt.plot(t, S, label='Living')
plt.plot(t, Z, label='Zombies')
plt.xlabel('Days from outbreak')
plt.ylabel('Population')
plt.title('Zombie Apocalypse - No Init. Dead Pop.; No New Births.')
plt.legend(loc=0)
plt.show()

是否可以使用odeint 来求解随机 ODE 系统？例如，如果我想在方程式的出生率 (P) 中包含一个错误项/随机游走？

我的想法是在系统中使用一个额外的方程来定义随机游走（随机抽样死亡率（使用 random.normalvariate()）并像这样求解系统：

 f0 = P - B*Si*Zi - f3*Si
 f1 = B*Si*Zi + G*Ri - A*Si*Zi
 f2 = f3*Si + A*Si*Zi - G*Ri
 f3 = random.normalvariate(mu, sigma)
 return [f0, f1, f2]

这是解决 SDE 系统的正确方法吗？还是我必须为随机 ODE 使用不同的求解器？

【问题讨论】：

是的，您必须采用不同的求解器，因为对于时间步 dt，随机游走步的大小应为 sqrt(dt)，ODE 求解器没有针对不同时间尺度的设施。查找 Euler-Maruyama 和 Milstein 方法，SDE 的高阶方法成倍增加。
我和github.com/mattja/sdeint玩过一次，结果很好

标签： python ode stochastic

【解决方案1】：

在帮助下，ODE 系统被改写为 SDE 系统，其中出生率是一个随机过程。

使用SDEint 包是一个很好的建议。

# Zombie apocalypse SDE model

import matplotlib.pyplot as plt
import numpy as np
import sdeint

P, d, B, G, A = 0.0001, 0.0001, 0.0095, 0.0001, 0.0001

tspan = np.linspace(0, 5., 1000)
y0 = np.array([500., 0., 0., P])


def f(y, t):
    Si = y[0]
    Zi = y[1]
    Ri = y[2]

    f0 = y[3] - B * Si * Zi - d * Si
    f1 = B * Si * Zi + G * Ri - A * Si * Zi
    f2 = d * Si + A * Si * Zi - G * Ri
    f3 = 0
    return np.array([f0, f1, f2, f3])


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

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

plt.plot(result)
plt.show()

【讨论】：