质量绘制从参数值列表求解的 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)

但这不正确。我不知道我是否以错误的方式思考它或语法不正确而不是逻辑。是否有可能获取列表中的每个 beta 值 empty，为其计算 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)

然后，您可以按照与您的代码片段类似的方式继续操作，即

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.

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

Mass plotting ODE curves that are solved from a list of parameter values

python

scipy

ode