15.9 ロジスティックモデル(コロナウィルスの拡散)

ここでは，感染症(コロナウイルスを想定している)の拡散状況を考えます．
都市の人口を$a$とし，感染症の罹患数を$x$とする．$x$が時間とともにどのように 増加するかについて調べる．
$x$の増加数に注目する．$x$が2倍，3倍にならば，$x$の増加量も2倍，3倍になるから， $x$の増加量は$x$に比例する．
したがって
$\dfrac{dx}{dt}=kx \cdots (\sharp)$
と表せる．ここで，$k$を増殖率と呼ぶ．
増殖率$k$が一定なら簡単に解決しますが，そうではありません．実際には，
最初気が付かないうちは一定値で推移しますが，この感染症に気づき右肩下がりになり， その後学校の閉鎖，外出禁止，都市封鎖によりさらに下がり続け，都市住民のほとんどが罹患して $k=0$になる．
となるのでしょう．専門家はデータを取り寄せ，増殖率を調べるのでしょうが， ここではもっと単純にします．
都市の人口を$n$とし，増殖率$k$は
$x=0$のとき$k=a$
$x=n$のとき$k=0$
かつ１次関数的変化をすると仮定する．すなわち
$k=a-\dfrac{a}{n}x=a(1-\dfrac{x}{n})$
とする．
$(\sharp)$より$x(t)$は微分方程式
$\dfrac{dx}{dt}=a(1-\dfrac{x}{n})x\cdots (\flat)$
を満たすことが分かる．
微分方程式$(\flat)$は後で解くとして，まず漸化式の立場から$(\flat)$を眺めよう．
$b=\dfrac{1}{n}$
とおき，$(\flat)$を
$\dfrac{dx}{dt}=a(1-bx)x$
とする．
$\dfrac{dx}{dt}=\dfrac{x(t+h)-x(t)}{h}$
と近似し
$\dfrac{x(t+h)-x(t)}{h}=a(1-bx(t))x(t)$
より
$x(t+h)=(1+ah)x(t)(1-\dfrac{abh}{1+ah}x(t))$
両辺に$\frac{abh}{1+ah}$をかけて
$\dfrac{abh}{1+ah}x(t+h)=(1+ah)\dfrac{abh}{1+ah}x(t)(1-\dfrac{abh}{1+ah}x(t))$
ここで
$y(t)=\dfrac{abh}{1+ah}x(t),A=1+ah$
とおくと
$y(t+h)=Ay(t)(1-y(t))$
を得る．
$h$として単位時間をとり
$a_n=y(n)$
により数列$\{a_n\}$を定めれば，
$a_{n+1}=Aa_n(1-a_n)\cdots (\natural)$
となる．これがロジスティックモデルの漸化式です．
この式の変形を見てわかるように$x(t)$が
$\dfrac{dx}{dt}=f(x)$
のとき(右辺が$t$を含まないとき)，そのモデルの漸化式が作れることに注意しましょう．
漸化式$(\natural)$は多くの研究がなされており，すでに紹介したネット上に優れた解説があります．

つぎに，ロジスティックモデルの微分方程式
$\dfrac{dx}{dt}=a(1-bx)x$
を解きましょう．
微分方程式は，高校では教科書の最後に簡単に書いているだけです．知らない人もいるでしょうから まず基本事項から始めます．
微分方程式を解くとは積分をすることです．
例えば，
$y=x^3+2x^2$
のとき
$\dfrac{dy}{dx}=3x^2+4x$
で$\frac{dy}{dx}$が$x$だけの式なら問題はありません．$y$を求めるのは単に積分するだけです．
しかし，例えば$x^2+y^2=1$の微分は
$2x+2y\dfrac{dy}{dx}=0$
より
$\dfrac{dy}{dx}=-\dfrac{x}{y}$
となり，$\frac{dy}{dx}$が$x,y$両方用いて表される場合も少なくありません．
また，円を
$x=\cos \theta ,y=\sin \theta $
と表すなら，
$x'(\theta)=-\sin \theta ,y'(\theta)=\cos \theta$
より
$\dfrac{dy}{dx}=-\dfrac{1}{\tan \theta}$
と媒介変数表示される場合もあります．
例えば,2次導関数が
$\dfrac{d^2y}{dx^2}=x+y$
などもあります．このように$\frac{dy}{dx}$が$x,y$やさらにパラメータを 含む式で与えられている場合に$y$ を求めるというのが微分方程式を解くということです．
微分方程式の解法は沢山ありますが，その中で必要なのは変数分離法です．
変数分離法とは$\frac{dy}{dx}$が$x$だけの式と$y$だけの式の積または商で表せる場合です．
すなわち
$\dfrac{dy}{dx}=f(x)g(y),\dfrac{f(x)}{g(x)},\dfrac{1}{f(x)g(x)}$
の場合です．
解き方は
$\dfrac{dy}{dx}=f(x)g(y)$
で説明すると．
$g(y)$を左辺に移項して
$\dfrac{1}{g(y)}\dfrac{dy}{dx}=f(x)$
両辺を$x$で積分すると
$\displaystyle \int \dfrac{1}{g(y)}\dfrac{dy}{dx}dx=\int f(x)dx$
左辺は置換積分より
$\displaystyle \int \dfrac{1}{g(y)}dy=\int f(x)dx$
となり，両辺をそれぞれ積分し
$G(y)=F(x)+C$($C$は積分定数)
初期条件が与えられていれば，$C$を求め，必要であれば$y$を$x$について解けばよい．
ここで述べた一連の解法は
$\dfrac{dy}{dx}=f(x)g(y)$
を
$\dfrac{dy}{g(y)}=f(x)dx$
と変形し，両辺を積分し
$\displaystyle \int \dfrac{1}{g(y)}dy=\int f(x)dx$
としてよいことを意味しています．

図形や物理現象を微分を用いた関係式で表し，それが微分方程式で，それを解くことで問題が 解決されることがよくあります．
例を挙げましょう．
例1
各点における法線が原点を通る曲線を求めよう．
解 曲線上の点$(x,y)$における法線の方程式は
$Y-y=-\dfrac{1}{\frac{dy}{dx}}(X-x)$
原点を通るから
$-y=\dfrac{1}{\frac{dy}{dx}}x$ したがって,微分方程式
$\dfrac{dy}{dx}=-\dfrac{x}{y}$
を得る．
$\int ydy=-\int xdx$
$\dfrac{y^2}{2}=-\dfrac{x^2}{2}+C$
したがって，
$x^2+y^2=C (C > 0)$

例2
一方が固定され，もう一方に質量$m$のおもりが付いているばねが滑らかな水平面上の 置かれている．
自然の長さより一定の長さだけ伸ばし静かに離した時の質量の運動を考える．
時刻$t$におけるばねの伸びた長さを$x(t)$とすると，
ばねがおもりに加える力$F$は向きも考慮して
$F=-kx (k > 0)$,ここで$k$はばね定数．
一方運動方程式より
$F=m\cdot \dfrac{d^2x}{dt^2}$
$\therefore m\dfrac{d^2x}{dt^2}=-kx\cdots (\sharp)$
これがばねの振動の微分方程式です．
$a=\sqrt{\dfrac{k}{m}}$
とおくと$(\sharp)$は
$\dfrac{d^2x}{dt^2}=-a^2x$
この微分方程式の解は
$x_1(t)=\sin at,x_2(t)=\cos at$
がともに微分方程式を満たすから
$x=A\sin at+B\cos at =C\sin (at+\alpha)$
が解であることがわかる．
もちろんこれだけでは不十分である．解が他にないことを示す必要があるが，ばねの振動の 理解にはこれで十分である．

このように，図形や運動の性質を微分方程式で表し，それを解く．似たようなことを高校で学んでいます． そう軌跡です．軌跡では，図形の性質を図形上の点$(x,y)$の$x,y$の関係式で表しました． それを導関数$y'$や$y''$を用いただけです．

準備が整いました．ロジスティックモデルの微分方程式
$\dfrac{dx}{dt}=a(1-\dfrac{x}{n})x\cdots (\flat)$
を解きましょう．ちなみに，これは人口$n$（万）人，最初の増殖率が$a$で，増殖率が罹患者とともに， 1次関数的に減少するモデルです．
$(\flat)$を
$\dfrac{1}{(1-\frac{x}{n})x}dx=adt$
と変形し，積分する．左辺を部分分数に変形し
$\int (\dfrac{1}{n-x}+\dfrac{1}{x})dx=\int adt$
$-\log (n-x)+\log x=at+C$
$\log \dfrac{x}{n-x}=at+C$
$\dfrac{x}{n-x}=Ce^{at}$($e^C$を新しく$C$と置いた)
あとは$x$について解けば
$x=\dfrac{n}{1+Ce^{-at}}$
参考
微分方程式は，いろいろな解法があります．
微分方程式$(\flat)$は
$s=\dfrac{1}{x}$
と変数変換すると，少し簡単になります．
$\dfrac{ds}{dt}=\dfrac{ds}{dx}\dfrac{dx}{dt}$
$=-\dfrac{1}{x^2}a(1-\dfrac{x}{n})x$
$=-a(s-\dfrac{1}{n})$
したがって
$\displaystyle \int \dfrac{ds}{s-\dfrac{1}{n}}=-\int adt$
$\log (s-\dfrac{1}{n})=-at+C$
$s-\dfrac{1}{n}=Ce^{-at}$
$\dfrac{n-x}{nx}=Ce^{-at}$
$\dfrac{nx}{n-x}=\dfrac{e^{at}}{C}$
積分定数を変換すれば上で求めた式と一致する．

話をすすめます．初期条件を代入します．
$t=0$のとき$x=b$とし定数$C$を定めます．
$\dfrac{n}{1+C}=b$
より
$C=\dfrac{n-b}{b}$
$C$を代入し
$x=\dfrac{n}{1+\dfrac{(n-b)e^{-at}}{b}}$
$=\dfrac{bn}{b+(n-b)e^{-at}}$
この式が求めるものです．

グラフの形を調べましょう．
東京を考え
$n=1000(\text{万人}),a=1,b=0.5$
として，すなわち $t=0$のとき$5000$人，その時の増殖率を$1$としてその後の変化を図示したのが右図です．
さらに，$0\leqq t\leqq 4$では実際の値が分かりづらいので，グラフの下に$0.1$間隔で 数値を表記しました．
$y=0,y=1000$が漸近線のグラフで，$t=4$ぐらいから指数関数的に増加し，途中変曲点を迎え増加量は減少していきます．

日本では4月7日に非常事態宣言が発動されました．その時の判明している東京の感染者数が約$1000$名．実際の感染者数を $50000$人と考え，感染者数推移のグラフを作りました．
問題は増殖率ですが，資料がないので，公表されている東京の感染者数の推移より推測し， 4月7日の増殖率を$0.2$すなわち，感染者は ５日に一人うつすとしました．
増殖率が$a=0.2$(左図)とその8割減の$a=0.04$のグラフ（右図）です．












一見グラフは似ていますが縦座標のメモリを見ればわかるように 全く異なることが分かります．
実際，例えば20日後を見てみましょう（横座標が$20$のときの縦座標），増殖率が$0.2$の時は約$200$万人が感染していますが， 増殖率が$0.04$では約$10$万人です．
実際には，治っている人数を引くので，増殖率が$0.04$の時は直線的か減少に転ずるのでしょう．

専門家が助言を与えているのでしょうが，とても冷静に判断していることが分かります．

実際の増殖率は右図のように，要所要所で不連続に下がっていくと思われるので $a,b,c$で新しいグラフに置き換わるグラフになるのでしょう．

グラフや表は数式ソフトMathmatica を利用しました．入力式もシンプルにし表示しました． 興味があれば楽しんでください．