利用者:Aozoramat/sandbox

微分方程式モデルの基本知識

最も簡単な微分方程式

簡単な１階常微分方程式は、

$\operatorname {d} \!y/\operatorname {d} \!t=ay$

は基本中の基本と考えられる微分方程式である。変数記号 $t$ は時間を意味し、 環境が良い場合の生物の増殖 (人口増加に関するマルサスの法則)、 放射性元素の崩壊、 単純化した状況下での物体の冷却など、様々な現象を記述する。

この解は次のようになる。

$y=Ce^{ax}$

この方程式の解き方

$y=0$ のときは、この方程式は満足されるのは明らかである。

$y\neq 0$ のときは、まず両辺を $y$ で除算する。

${\frac {1}{y}}$${\frac {dy}{dt}}$$=$$a$

${\frac {1}{y}}$$dy$$=$$a$$dt$

$\int {\frac {1}{y}}dy=\int adt$

$\log |y|+C_{1}=at+C_{2}$

$\log |y|=\pm at+C_{3}$

$y=e^{\pm at+C_{3}}=e^{\pm C_{3}}\times e^{at}$

$e^{\pm C_{3}}=C$ と置き換えると

$y=Ce^{at}\qquad (C\neq 0)$

$C$：任意定数

練習問題１

方程式　${dx \over dt}=rx(1-{x \over K})$ の解析解を求める

${dx \over dt}=-{r \over K}x(x-K)$

$\int {dx \over x(x-K)}={r \over K}\int xdt$

$-{1 \over K}\int {\Bigl (}{1 \over x}-{1 \over x-K}{\Bigl )}dx={r \over K}\int xdt$

この場合、$x\geqq 0\cap K>0$なので

$x\geqq K$ のときは、

$\log {\Bigl (}{x \over x-K}{\Bigl )}=-rx+C_{1}$

${x \over x-K}=e^{-rx+C_{1}}$

$x={e^{rt+C_{1}}K \over e^{rt+C_{1}}-1}={Ce^{rt}K \over Ce^{rt}-1}$　

$x(0)={CK \over C-1}$

$C={x(0) \over x(0)-K}$

これを使うと、

$x(t)={Ce^{rt}K \over Ce^{rt}-1}$

${\begin{aligned}x&={x(0)e^{rt}K \over x(0)e^{rt}-(x(0)-K)}\\&={Kx(0)e^{rt} \over K+x(0)(e^{rt}-1)}\end{aligned}}$

$0\leqq x<K$ のときは、

$\log {\Bigl (}{x \over -x+K}{\Bigl )}=-rx+C_{1}$

${x \over x-K}=-e^{-rx+C_{1}}$

$x={e^{rt+C_{1}}K \over e^{rt+C_{1}}+1}={Ce^{rt}K \over Ce^{rt}+1}$

$x(0)={CK \over C+1}$

$C={-x(0) \over x(0)-K}$

これを使うと、

$x(t)={Ce^{rt}K \over Ce^{rt}-1}$

${\begin{aligned}x&={-x(0)e^{rt}K \over x(0)e^{rt}-(x(0)-K)}\\&={-Kx(0)e^{rt} \over K+x(0)(e^{rt}-1)}\end{aligned}}$

練習問題２

方程式　${dx \over dt}=Rx(1-x)$ の解析解を求める

$\int {dx \over x(1-x)}=R\int dt$

$\int {\Bigl (}{1 \over x}+{1 \over 1-x}{\Bigl )}dx=R\int dt$

$\int {\Bigl (}{1 \over x}-{1 \over x-1}{\Bigl )}dx=R\int dt$

$1\leqq x$ のときは、

$\log(x)-\log(x-1)=\log {x \over x-1}=Rt+C_{1}$

${x \over x-1}=e^{Rt+C_{1}}$

$x(R,t)={e^{Rt+C_{1}} \over e^{Rt+C_{1}}-1}={Ce^{Rt} \over Ce^{Rt}-1}$

$x<1$ のときは、

$\log(x)-\log(1-x)=\log {x \over 1-x}=Rt+C_{2}$

${x \over 1-x}=e^{Rt+C_{2}}$

$x(R,t)={e^{-Rt+C_{2}} \over e^{-Rt+C_{2}}+1}={Ce^{-Rt} \over Ce^{-Rt}+1}$

　

ppppppppppppppppppppppppppppppppppppppppppppp

$x(0)={CK \over C-1}$

$C={x(0) \over x(0)-K}$

これを使うと、

$x(t)={Ce^{rt}K \over Ce^{rt}-1}$

${\begin{aligned}x&={x(0)e^{rt}K \over x(0)e^{rt}-(x(0)-K)}\\&={Kx(0)e^{rt} \over K+x(0)(e^{rt}-1)}\end{aligned}}$

$0\leqq x<K$ のときは、

$\log {\Bigl (}{x \over -x+K}{\Bigl )}=-rx+C_{1}$

${x \over x-K}=-e^{-rx+C_{1}}$

$x={e^{rt+C_{1}}K \over e^{rt+C_{1}}+1}={Ce^{rt}K \over Ce^{rt}+1}$

$x(0)={CK \over C+1}$

$C={-x(0) \over x(0)-K}$

これを使うと、

$x(t)={Ce^{rt}K \over Ce^{rt}-1}$

${\begin{aligned}x&={-x(0)e^{rt}K \over x(0)e^{rt}-(x(0)-K)}\\&={-Kx(0)e^{rt} \over K+x(0)(e^{rt}-1)}\end{aligned}}$

$\log {\Bigl (}{x \over x-K}{\Bigl )}=-rx+C_{1}$

$dx=-{r \over K}x(x-a)(x-K)dt$

$dx=-{r \over K}x(x-a)(x-K)dt$

アリー効果の個体群動態方程式

方程式　${dx \over dt}=rx(x-a)(1-{x \over K})$ の解析解を求める。

これは直接積分法で解くことができる。

$dx=-{r \over K}x(x-a)(x-K)dt$

$x=0,x=a,x=K$ のいずれかのときは、与式は明らかに満足される。 したがって、${\textstyle x\nsubseteqq \{0,a,K\}}$ のときについて解く。

${dx \over x(x-a)(x-K)}=-{r \over K}dt$ $\int {1 \over x(x-a)(x-K)}dx=-{r \over K}\int dt$

左辺の分数を部分分数に分解する

${u \over x}+{v \over x-a}+{w \over x-K}$$={u(x-a)(x-K)+vx(x-K)+wx(x-a) \over x(x-a)(x-K)}$

上式の分子を$P$とおき、これを整理する。

${\begin{aligned}P&=u(x-a)(x-K)+vx(x-K)+wx(x-a)\\&=ux^{2}-uax-uKx+uaK+vx^{2}-vKx+wx^{2}-wax\\&=(u+v+w)x^{2}-(ua+uK+vk+wa)x+uaK\end{aligned}}$

$P=1$　であるので、次の連立方程式を得る。

連立方程式

$u+v+w=0\ldots \ldots \ldots (1)$

$ua+uK+vK+wa=0\qquad (2)$

$uaK=1\qquad \qquad \qquad (3)$

$(3)$より、$u={1 \over aK}$

$(2)$に代入して、

${1 \over K}+{1 \over a}+vK+wa=0$

これにより、$w=-{1 \over aK}-{1 \over a^{2}}-{vK \over a}$ を得る。これを$(1)$に代入

${1 \over aK}+v-{1 \over aK}-{1 \over a^{2}}-{vK \over a}=0$

$v(1-{K \over a})=-{1 \over aK}+{1 \over aK}+{1 \over a^{2}}={1 \over a^{2}}$

${(a-K) \over a}v={1 \over a^{2}}$

$v={1 \over a(a-K)}$

ここから

${\begin{aligned}w=-u-v&=-{1 \over aK}-{1 \over a(a-K)}\\&=-{a-K \over aK(a-K)}-{K \over aK(a-K)}\\&={-a+K-K \over a(a-K)K}\\&={1 \over K(K-a)}\\\end{aligned}}$

これで、連立方程式の解が求まった。

$u={1 \over aK}$

$v={1 \over a(a-K)}$

$w={1 \over K(K-a)}$

検算

$W={u \over x}+{v \over x-a}+{w \over x-K}={1 \over x(x-a)(x-K)}$ であることを確かめる

${\begin{aligned}W&={1 \over aKx}+{1 \over a(a-K)(x-a)}+{1 \over K(K-a)(x-K)}\\&={(K-a)(x-a)(x-K) \over aK(K-a)x(x-a)(x-K)}+{-Kx(x-K) \over aK(K-a)x(x-a)(x-K)}+{ax(x-a) \over aK(K-a)x(x-a)(x-K)}\end{aligned}}$

${\begin{aligned}n1=(Kx-ax-aK+a^{2})(x-K)&=Kx^{2}-ax^{2}-aKx+a^{2}x-K^{2}x+aKx+aK^{2}-a^{2}K\\&=Kx^{2}-ax^{2}+a^{2}x+-K^{2}x+aK^{2}-a^{2}K\end{aligned}}$

${\begin{aligned}n2=-Kx(x-K)=-Kx^{2}+K^{2}x\end{aligned}}$$n3=ax(x-a)=ax^{2}-a^{2}x$

${\begin{aligned}n1+n2+n3&=aK^{2}-a^{2}K\\&=aK(K-a)\end{aligned}}$

${\begin{aligned}{n1+n2+n3 \over aK(K-a)x(x-a)(x-k)}&={aK(K-a) \over aK(K-a)x(x-a)(x-k)}\\&={1 \over x(x-a)(x-K)}\end{aligned}}$

以上で検算も成功した。

微分方程式の解

$\int {dx \over aKx}+\int {dx \over a(a-K)(x-a)}+\int {dx \over K(K-a)(x-k)}=-{r \over K}\int dt$

${1 \over Ak}\int {dx \over x}+{1 \over a(a-K)}\int {dx \over x-a}+{1 \over K(K-a)}\int {dx \over (x-k)}=-{r \over K}\int dt$

$\left\{{1 \over Ak}\log(|x|)+C_{1}\right\}+\left\{{1 \over a(a-K)}\log(|x-a|)+C_{2}\right\}+\left\{{1 \over K(K-a)}\log(|x-K|)+C_{3}\right\}=-{r \over K}\int dt$

${1 \over Ak}\log(|x|)+{1 \over a(a-K)}\log(|x-a|)+{1 \over K(K-a)}\log(|x-K|)+C_{4}=-{r \over K}\int dt$

積分公式　$r\log M=logM^{r}$ $r$:実数、 $M>0$

　　　　　 $\log M+\log N=log(M\cdot N)$

より

$\log \left\{x^{1 \over ak}(x-a)^{1 \over a(a-K)}(x-K)^{1 \over K(K-a)}\right\}+C4=-{r \over K}t+C_{5}$

$\log \left\{x^{1 \over ak}(x-a)^{1 \over a(a-K)}(x-K)^{1 \over K(K-a)}\right\}=-{r \over K}t+C_{6}$

$x^{1 \over ak}(x-a)^{1 \over a(a-K)}(x-K)^{1 \over K(K-a)}=e^{{-r \over K}t+C_{6}}=Ce^{{-r \over K}t}$

これは解けない!!!!!!!!!!!!!!!!!!!

$A$

$A$

$A$

$A$

${\begin{aligned}&{\frac {dS}{dt}}=-{\frac {\beta IS}{N}},\\[6pt]&{\frac {dI}{dt}}={\frac {\beta IS}{N}}-\gamma I,\\[6pt]&{\frac {dR}{dt}}=\gamma I,\end{aligned}}$

${\begin{alignedat}{2}f(x)&=(a+b)^{2}\\&=a^{2}+2ab+b^{2}\\\end{alignedat}}$

The dynamics of an epidemic, for example, the flu, are often much faster than the dynamics of birth and death, therefore, birth and death are often omitted in simple compartmental models. The SIR system without so-called vital dynamics (birth and death, sometimes called demography) described above can be expressed by the following set of ordinary differential equations:

${\begin{aligned}&{\frac {dS}{dt}}=-{\frac {\beta IS}{N}},\\[6pt]&{\frac {dI}{dt}}={\frac {\beta IS}{N}}-\gamma I,\\[6pt]&{\frac {dR}{dt}}=\gamma I,\end{aligned}}$

where $S$ is the stock of susceptible population, $I$ is the stock of infected, $R$ is the stock of removed population (either by death or recovery), and $N$ is the sum of these three.

This model was for the first time proposed by William Ogilvy Kermack and Anderson Gray McKendrick as a special case of what we now call Kermack–McKendrick theory, and followed work McKendrick had done with Ronald Ross.

This system is non-linear, however it is possible to derive its analytic solution in implicit form. Other numerical tools include Monte Carlo methods, such as the Gillespie algorithm.

Firstly note that from:

${\frac {dS}{dt}}+{\frac {dI}{dt}}+{\frac {dR}{dt}}=0,$

it follows that:

$S(t)+I(t)+R(t)={\text{constant}}=N,$

expressing in mathematical terms the constancy of population $N$. Note that the above relationship implies that one need only study the equation for two of the three variables.

Secondly, we note that the dynamics of the infectious class depends on the following ratio:

$R_{0}={\frac {\beta }{\gamma }},$

the so-called basic reproduction number (also called basic reproduction ratio). This ratio is derived as the expected number of new infections (these new infections are sometimes called secondary infections) from a single infection in a population where all subjects are susceptible. This idea can probably be more readily seen if we say that the typical time between contacts is $T_{c}=\beta ^{-1}$, and the typical time until removal is $T_{r}=\gamma ^{-1}$. From here it follows that, on average, the number of contacts by an infectious individual with others before the infectious has been removed is: $T_{r}/T_{c}.$

By dividing the first differential equation by the third, separating the variables and integrating we get

$S(t)=S(0)e^{-R_{0}(R(t)-R(0))/N},$

where $S(0)$ and $R(0)$ are the initial numbers of, respectively, susceptible and removed subjects. Writing $s_{0}=S(0)/N$ for the initial proportion of susceptible individuals, and $s_{\infty }=S(\infty )/N$ and $r_{\infty }=R(\infty )/N$ for the proportion of susceptible and removed individuals respectively in the limit $t\to \infty ,$ one has

$s_{\infty }=1-r_{\infty }=s_{0}e^{-R_{0}(r_{\infty }-r_{0})}$

(note that the infectious compartment empties in this limit). This transcendental equation has a solution in terms of the Lambert W function, namely

$s_{\infty }=1-r_{\infty }=-R_{0}^{-1}\,W(-s_{0}R_{0}e^{-R_{0}(1-r_{0})}).$

This shows that at the end of an epidemic, unless $s_{0}=0$, not all individuals of the population have been removed, so some must remain susceptible. This means that the end of an epidemic is caused by the decline in the number of infectious individuals rather than an absolute lack of susceptible subjects.

The role of the basic reproduction number is extremely important. In fact, upon rewriting the equation for infectious individuals as follows:

${\frac {dI}{dt}}=\left(R_{0}{\frac {S}{N}}-1\right)\gamma I,$

it yields that if:

$R_{0}>{\frac {N}{S(0)}},$

then:

${\frac {dI}{dt}}(0)>0,$

i.e., there will be a proper epidemic outbreak with an increase of the number of the infectious (which can reach a considerable fraction of the population). On the contrary, if

$R_{0}<{\frac {N}{S(0)}},$

then

${\frac {dI}{dt}}(0)<0,$

i.e., independently from the initial size of the susceptible population the disease can never cause a proper epidemic outbreak. As a consequence, it is clear that the basic reproduction number is extremely important.

The force of infection

Note that in the above model the function:

$F=\beta I,$