用python 写出ｄＳ／ｄｔ＝－β１ＳＩ－β２ＳＥ－φＳ， ｄＶ／ｄｔ＝φＳ－ １－（ ）σ （ ） β１Ｉ＋β２Ｅ Ｖ， ｄＥ／ｄｔ＝β１ＳＩ＋β２ＳＥ＋ １－（ ）σ （ ） β１Ｉ＋β２Ｅ Ｖ－αＥ， ｄＩ／ｄｔ＝αＥ－γＩ， ｄＲ／ｄｔ＝γＩ ■ ■

import numpy as np
from scipy.integrate import odeint

# 定义ODE方程
def model(z, t, beta1, beta2, phi, sigma, alpha):
S, V, E, I, R = z
dSdt = -beta1*S*I - beta2*S*E - phi*S
dVdt = phi*S - (1-sigma*V)*(beta1*I+beta2*E)*V
dEdt = beta1*S*I + beta2*S*E + (1-sigma*V)*(beta1*I+beta2*E)*V - alpha*E
dIdt = alpha*E - gamma*I
dRdt = gamma*I
return [dSdt, dVdt, dEdt, dIdt, dRdt]

# 初始值
S0 = 0.99
V0 = 0.01
E0 = 0
I0 = 0
R0 = 0

# 参数设定
beta1 = 0.4
beta2 = 0.3
phi = 0.1
sigma = 0.1
alpha = 0.2
gamma = 0.1

# 时间点
t = np.linspace(0, 100, 1000)

# 解ODE方程
z = odeint(model, [S0, V0, E0, I0, R0], t, args=(beta1, beta2, phi, sigma, alpha))

# 画图
import matplotlib.pyplot as plt
plt.plot(t, z[:,0], 'r', label='S')
plt.plot(t, z[:,1], 'g', label='V')
plt.plot(t, z[:,2], 'b', label='E')
plt.plot(t, z[:,3], 'y', label='I')
plt.plot(t, z[:,4], 'm', label='R')
plt.legend(loc='best')
plt.xlabel('t')
plt.ylabel('Population')
plt.show()