代码如下：
import ode_7 as ode
import numpy as np
import matplotlib.pyplot as plt
# 采用matplotlib作图时默认设置下是无法显示中文的，凡是汉字都会显示成方块。
# 实际上，matplotlib是支持unicode编码的，不能正常显示汉字主要是没有找到合适的中文字体。
from pylab import mpl
mpl.rcParams['font.sans-serif']= ['SimHei']
# 解决负号显示问题
import matplotlib
matplotlib.rcParams['axes.unicode_minus']=False
'''
根据传染病传播的SIR模型，微分方程组可以转换为向量运算：
y = [S, I, R]
f(t, y) = [-beta*y[0]*y[1], beta*y[0]*y[1]-gamma*y[1], gamma*y[1]]
前向欧拉公式：
y[k+1] = y[k] + dt * f(t[k], y[k])
后向欧拉公式（两步迭代）：
y_p = y[k] + dt * f(t[k], y[k])
y[k+1] = y[k] + dt * f(x[k+1], y_p)
改进的欧拉法：
y_p = y[k] + dt * f(t[k], y[k])
y_c = y[k] + dt * f(t[k+1], y_p)
y[k+1] = 1/2 * (y_p + y_c)
'''
class SIR(ode.ODE):
def f(self, t, y): # Y为向量[S,I,R]，t为标量
#请在此处编程，实现f函数
########### Begin ###########
return np.array([-self.beta*y[0]*y[1], self.beta*y[0]*y[1]-self.gamma*y[1], self.gamma*y[1]])
########### End ###########
def __init__(self, beta, gamma, y0):
#请在此处编程，调用ODE类对SIR进行初始化
########### Begin ###########
self.beta = beta; self.gamma = gamma; self.y0 = y0 # 参数属性
ode.ODE.__init__(self, self.f, self.y0)
########### End ###########
N = 1e8 # 武汉总人数：1000万人
beta = 1.0/N # 1名病毒携带者平均每天感染1.0人
gamma = 1/25 # 假设肺炎患者平均25天治愈（15天潜伏 + 10天治疗）
y0 = [N-1, 1, 0] # 初始发病1人，其他人员正常 [S0, I0, R0]
t = np.arange(0, 61, 1) # 模拟60天的发展情况，单位时间为1天
simulation = SIR(beta=beta, gamma=gamma, y0=y0)
y = simulation.solve(t, 'T') # 'T' 代表改进的欧拉法
print(y)
plt.plot(t, y[:, 0], '-g', label='易感者')
plt.plot(t, y[:, 1], '-r', label='感染者')
plt.plot(t, y[:, 2], '-k', label='恢复（隔离）者')
plt.legend(loc=0)
plt.title('继承ODE类的新冠病毒传播分析（SIR模型）')
plt.xlabel('时间（天）')
plt.ylabel('人数')
Vx = [1.0e6, 5.0e6] + [i * 1.e7 for i in range(11)]
plt.yticks(Vx, ['%d'%e for e in Vx])
plt.savefig('test1/test1.png')
|