从参数值列表求解的批量绘制 ODE 曲线

Question

我的代码是一个 ODE 系统，我可以为其绘制曲线。但是，我希望生成 100 个不同的解决方案，然后我可以将所有 100 个作为曲线绘制到一个图形上。我的代码是：

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

N = 1000

I0, R0 = 1, 0

S0 = N - I0 - R0
J0 = I0

beta, gamma = 2/7, 1/7

t = np.linspace(0, 160, 160+1)


def deriv(y, t, N, beta, gamma):
    S, I, R, J = y
    dS = ((-beta * S * I) / N)
    dI = ((beta * S * I) / N) - (gamma * I)
    dR = (gamma * I)
    dJ = ((beta * S * I) / N)
    return dS, dI, dR, dJ


solve = odeint(deriv, (S0, I0, R0, J0), t, args=(N, beta, gamma))
S, I, R, J = solve.T

J_diff = np.diff(J)

fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)
ax.plot(t[1:], J_diff, 'blue', alpha=1, lw=2, label='Daily incidence')
ax.set_xlabel('Time in days')
ax.set_ylabel('Number (1000s)')
ax.grid(b=True, which='major', c='w', lw=2, ls='-')
legend = ax.legend()
legend.get_frame().set_alpha(0.5)
plt.show()

这会产生这样的曲线

但是，我希望为随机 100 个 beta 值的列表计算 ODE 方程 J，所有其他参数保持固定，我可以绘制这些参数以获得 100 条曲线。我尝试过的是使用类似的东西：

empty = []
for i in range(100):
    empty.append(random.uniform(1.5, 2.5)*gamma)

对于新的 beta 值列表，然后做

solns =  []
for empt in empty:
    ces = odeint(deriv, (S0, I0, R0, J0), t, args=(N, empty, gamma))
    solns.append(ces)

但这是不正确的。我不知道我是否在考虑错误的方式或语法不正确，而不是逻辑。是否可以获取列表empty 中的每个 beta 值，为其计算 ODE 系统，并将每个解绘制到一个数字上？

Answer 1

你的方法是正确的。问题是由于args 参数中的一个小错字，您传递的是一个beta 值列表而不是单个标量。应该是

solns =  []
for empt in empty:
    ces = odeint(deriv, (S0, I0, R0, J0), t, args=(N, empt, gamma))
    solns.append(ces)

然后，你可以像你的代码 sn-p 一样继续，即

J_diffs = []
for sol in solns:
    S, I, R, J = sol.T
    J_diffs.append(np.diff(J))

# plot all the solutions
fig = plt.figure(facecolor='w')
ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)
ax.set_xlabel('Time in days')
ax.set_ylabel('Number (1000s)')
ax.grid(b=True, which='major', c='w', lw=2, ls='-')

for J_diff in J_diffs:
    ax.plot(t[1:], J_diff, 'blue', alpha=1, lw=2, label='Daily incidence')

# plot without legend
plt.show()

Answer 2

您可以通过将循环制作成列表生成器来缩短代码，

empty = [ random.uniform(1.5, 2.5)*gamma  for _ in range(100) ]

或者使用numpy方法

empty = np.random.uniform(1.5, 2.5, size=100)*gamma

和

solns =  [ odeint(deriv, (S0, I0, R0, J0), t, args=(N, bet, gamma))  for bet in empty ]

如果您设置自定义容差或其他参数，将这些常量参数封装到一个 mysolver 函数中可能是有意义的，该函数只有可能更改的参数作为参数，beta，可能还有 gamma 和 t .