佐藤總夫「自然の数理と社会の数理 II 微分方程式で解析する」（日本評論社）という名著があった。ホントはI, II, IIIと３分冊で出るはずが，I, IIが出た後は音沙汰なしという，当時学部生だったワシが残念に思っていたシリーズである。今回のコロナウィルス騒ぎで散々出てくるあの「流行曲線」という奴，本書に詳しく出てたので，折角だから知識と本の虫干しがてら，Python+SciPyで描いてみることにしよう。つーことで，以下の数学的な記述は全部佐藤著の「第７話 伝染病の伝播」に基づくものなので，詳しく知りたい人は図書館にGoである。理工系大学ならあるんじゃないのかな。
感染症モデルの中でも一番簡単なものを取り上げる。前提として
- 感染者を[latex]y = y(u)[/latex], 感受性者（これから感染する可能性のある人）を\(x(u) = x\)とする。
- 全人口を\(n+1\)とし，増減はないものとする（「死者は出ない」の意）。つまり常に\(x + y = n + 1\)が保持される。
- 感染者の発生率は\(xy\)に比例する（人々の均等な接触を表現）。
- 最初の感染者は\(y(0) = 1\)とする。
とする。この時，常微分方程式の初期値問題としては
\[ \frac{dx}{du} = -xy = -x(n+1 – x), x(0) = n \]
となる。新たにどの程度の感染者が現れるかを表現したもの，つまり感染スピードを表わす関数が\(w(u)\)である。上記を解いて\(x(u)\), \(y(u)\)を求め，そこから\(w(u)\)を計算してグラフ化したものが下記の図である。小うるさい人にすれば，人数が単位の\(x(u)\), \(y(u)\)と人/時間が単位の\(w(u)\)を同じ縦軸に描くとは何事だと言いそうであるが，これらの関数の関係を表したものとご理解頂きたい・・・今気がついたが，右軸も使えば良かったんだな，面倒くさいから，やりたい人向けの課題としておく。
このグラフから，感覚的には次のことはすぐに分かる。
- 感染を妨げるものはないので，最終的には全員に感染する。
- 感染が止まるのは全員に感染するからであって，感染するスピードは山なりに変化する。
初期感染者は今更ゼロにはできないし，できることといえば，感染者を見つけたら素早く感受性者から引き離して\(xy\) を減らすしかない。つまりは初期感染スピード\(w(0)\)を極力減らすしかない訳で，実際，\(w(0) = 10, 20, 30, 40\)まで\(10\)刻みで変えてみると下記のようなグラフを得る。
このグラフが感染スピードの変化を示しており，当然ながら，初期感染スピードが遅ければ遅いほど感染率極力小さく抑えられることが分かる。現在の厚労省が取っている対策もこれに基づくもので，このグラフはまさしくこのグラフが描いている曲線ということになる。
つーことで，上のグラフを描くためのPython + SciPyコード例を下記に挙げておく。ちなみに，この常微分方程式は解析解が手計算で求められるので，ピーク時はどこになるとかといった議論も可能である。詳細は佐藤著をご覧頂きたい。
# ode_epidemic.py: 常微分方程式の初期値問題
# numpy
import numpy as np
# integrateパッケージ
import scipy.integrate as spint
# matplotlib
import matplotlib.pyplot as plt
num_div_t = 100
windex = 0
w = [[0.0] * num_div_t for i in range(4)]
#print('w[0] = ', w[0])
# 人口
all_num_people = [11, 21, 31, 41]
for num_people in all_num_people:
# 初期感染者
num_infected_people = 1
# 初期感受性者
num_noninfected_people = num_people - num_infected_people
# 陽的形式の右辺
# y' = func(t, x) = -x(n+1-x)
def func(t, y):
# return [-y[0] * (num_people - y[0])]
return np.array([-y[0] * (num_people - y[0]), y[1] * (num_people - y[1])])
# 初期値
# y(0) = [n, 1]
#y0 = [num_noninfected_people]
y0 = [num_noninfected_people, num_infected_people]
# t = [0, 1]
t_interval = [0.0, 1.0]
#print(t_interval)
# 常微分方程式を解く
#ret = spint.solve_ivp(func, t_interval, y0)
ret = spint.solve_ivp(func, t_interval, y0, t_eval = np.linspace(t_interval[0], t_interval[1], num_div_t))
# 感染者の増加率を計算
w[windex][0] = num_noninfected_people # 初期感染者発生率
for index in range(1, ret.t.size):
w[windex][index] = (ret.y[1, index] - ret.y[1, index - 1]) / (ret.t[index] - ret.t[index - 1])
windex += 1
# 結果を表示
#print(ret)
# yを表示
#print(ret.y)
# wを表示
#print('w[0] = ', w[0])
# t-yグラフを描画
# グラフ初期化
figure, axis = plt.subplots()
# 値をセット
#axis.plot(ret.t, ret.y[0, :], label='x(u)')
#axis.plot(ret.t, ret.y[1, :], label='y(u)')
axis.plot(ret.t, w[0], label='w(0) = '+ str(all_num_people[0] - 1))
axis.plot(ret.t, w[1], label='w(0) = '+ str(all_num_people[1] - 1))
axis.plot(ret.t, w[2], label='w(0) = '+ str(all_num_people[2] - 1))
axis.plot(ret.t, w[3], label='w(0) = '+ str(all_num_people[3] - 1))
# x軸，y軸，グラフタイトルをセット
axis.set(xlabel = 't', ylabel = 'w', title = 'Epidemic model')
# グリッドを描画
axis.grid()
# 凡例
axis.legend()
# グラフ保存ファイル名
figure.savefig('ode_epidemic_w.png')
# グラフを画面に描画
plt.show()
