在 python 中使用正向欧拉法计算微分方程时出错

Question

我正在做一个项目，我使用 SEIRDV 模型预测 COVID 的行为。我已经创建了微分方程并设置了初始参数。我正在使用正向欧拉方法求解这些 ODE。当我创建一个程序以在 180 天的每一天查找 S、E、I、R、D、V 并使用 matplotlib 显示它时。我收到一个错误，如下所示。我相信这是由于以下操作。

I[i]/N[i]

我对 NumPy 有点陌生，因此不熟悉出现的错误。谁能告诉我如何解决这个问题？此外，我曾尝试使用 matplotlib 绘制数据，但数据点无处不在，而不是 SEIRDV 模型表示应有的样子。任何有关如何解决此问题的见解都将不胜感激。

import matplotlib.pyplot as plt
import numpy as np

ndays = 180 # time length for prediction
dt = 1 # time step in days

S = np.zeros(ndays) # susceptibles
E = np.zeros(ndays) # exposed
I = np.zeros(ndays) # infective
R = np.zeros(ndays) # recovered / removed
D = np.zeros(ndays) # fatalities
V = np.zeros(ndays) # vaccinated but have had the vaccine for less than 21 days
N = np.zeros(ndays) # total population
t = np.arange(ndays)*dt

S[0] = 45167146 # initial S
E[0] = 20000 # initial E
I[0] = 3085 # initial I
R[0] = 19744004 # initial R
V[0] = 65000 # initial V
D[0] = 0 # initial D
N[0] = 60500000 # initial N

beta = (0.75) # infection rate
gam = (0.1) # recovery rate
Clamb = (0.00005) # birth rate
mu = (0.00003) # natural death rate
epsi = (0.07) # rate at which someone goes from E to I
alpha = (0.1) # rate at which Infectives die from the virus
Ldelta = (0.0006) # rate at which people get vaccinated
zeta = (0.07) # rate at which someone goes from V to R

for i in range(ndays-1):
    S[i+1] = S[i] + (Clamb - mu*S[i] - beta*S[i]*(I[i]/N[i]) - Ldelta*S[i])*dt
    E[i+1] = E[i] + (beta*S[i]*(I[i]/N[i]) + beta*V[i]*(I[i]/N[i]) - mu*E[i] - epsi*E[i])*dt
    I[i+1] = I[i] + (epsi*E[i] - gam*I[i] - mu*I[i] - alpha*I[i])*dt
    R[i+1] = R[i] + (gam*I[i] + zeta*V[i] - mu*R[i])*dt
    D[i+1] = D[i] + (alpha*I[i])*dt
    V[i+1] = V[i] + (Ldelta*S[i] - beta*V[i]*(I[i]/N[i]) - zeta*V[i])*dt

display (S)
fig = plt.figure(1) 
plt.plot(t, S, "r", lw=2, label="Susceptible")
plt.plot(t, E, "y", lw=2, label="Exposed")
plt.plot(t, I, "g", lw=2, label="Infective")
plt.plot(t, R, "p", lw=2, label="Removed")
plt.plot(t, D, "b", lw=2, label="Fatalities")
plt.plot(t, V, "o", lw=2, label="Vaccinated before immunity")
fig.legend(); plt.xlabel("Days"); plt.ylabel("Amount of people")

这是我得到的错误

C:\Users\jared\AppData\Local\Temp/ipykernel_18008/1585353571.py:34: RuntimeWarning: divide by zero encountered in double_scalars
  S[i+1] = S[i] + (Clamb - mu*S[i] - beta*S[i]*(I[i]/N[i]) - Ldelta*S[i])*dt
C:\Users\jared\AppData\Local\Temp/ipykernel_18008/1585353571.py:35: RuntimeWarning: divide by zero encountered in double_scalars
  E[i+1] = E[i] + (beta*S[i]*(I[i]/N[i]) + beta*V[i]*(I[i]/N[i]) - mu*E[i] - epsi*E[i])*dt
C:\Users\jared\AppData\Local\Temp/ipykernel_18008/1585353571.py:39: RuntimeWarning: divide by zero encountered in double_scalars
  V[i+1] = V[i] + (Ldelta*S[i] - beta*V[i]*(I[i]/N[i]) - zeta*V[i])*dt
C:\Users\jared\AppData\Local\Temp/ipykernel_18008/1585353571.py:34: RuntimeWarning: invalid value encountered in double_scalars
  S[i+1] = S[i] + (Clamb - mu*S[i] - beta*S[i]*(I[i]/N[i]) - Ldelta*S[i])*dt
C:\Users\jared\AppData\Local\Temp/ipykernel_18008/1585353571.py:35: RuntimeWarning: invalid value encountered in double_scalars
  E[i+1] = E[i] + (beta*S[i]*(I[i]/N[i]) + beta*V[i]*(I[i]/N[i]) - mu*E[i] - epsi*E[i])*dt
C:\Users\jared\AppData\Local\Temp/ipykernel_18008/1585353571.py:39: RuntimeWarning: invalid value encountered in double_scalars
  V[i+1] = V[i] + (Ldelta*S[i] - beta*V[i]*(I[i]/N[i]) - zeta*V[i])*dt

Answer 1

您永远不会更新 N[i]，所以它是在第一步之后 0 .. 它可能应该在每次迭代时设置，也许 filled 从

开始

编辑：作为，直接设置 N[i] 并非易事 “添加此行” 并且需要做更多工作

获取前一个D[i-1]（可能从索引 1 开始）
切换N为累积生活
切换 D 立即死亡

i=0 S[i]=45167146.0 I[i]=3085.0 N[i]=60500000.0 dt=1
i=1 S[i]=45136963.33469715 I[i]=3867.90745 N[i]=0.0 dt=1

您可以嵌入 print() 或使用 PDB 来非常轻松地调试小程序

% python3 -m pdb test_70976944.py
> test_70976944.py(1)<module>()
-> import matplotlib.pyplot as plt
(Pdb) b 34
Breakpoint 1 at test_70976944.py:34
(Pdb) c
> test_70976944.py(34)<module>()
-> S[i+1] = S[i] + (Clamb - mu*S[i] - beta*S[i]*(I[i]/N[i]) - Ldelta*S[i])*dt
(Pdb) print(f"i={i} S[i]={S[i]} I[i]={I[i]} N[i]={N[i]} dt={dt}")
i=0 S[i]=45167146.0 I[i]=3085.0 N[i]=60500000.0 dt=1
(Pdb) c
> test_70976944.py(34)<module>()
-> S[i+1] = S[i] + (Clamb - mu*S[i] - beta*S[i]*(I[i]/N[i]) - Ldelta*S[i])*dt
(Pdb) print(f"i={i} S[i]={S[i]} I[i]={I[i]} N[i]={N[i]} dt={dt}")
i=1 S[i]=45136963.33469715 I[i]=3867.90745 N[i]=0.0 dt=1

在 python 中使用正向欧拉法计算微分方程时出错

Error in computing differential equations using the forward euler method in python

python

matplotlib

numerical-methods