数理モデル(mathematical model)とは、現象を表すために選んだ状態・入力・出力の変数、パラメータ、それらを結ぶ方程式、仮定と適用範囲の組である。量の出入りの収支式に構成則を加えると微分方程式が得られ、例えば 1 コンパートメントモデル $x'=-kx+r$ の解は定常値 $r/k$ に近づく。$x'=f(x)$ の平衡点 $x^*$ は $f'(x^*)<0$ なら漸近安定、$f'(x^*)>0$ なら不安定であり、一定量の収穫があるロジスティックモデルでは収穫が $rK/4$ を越えると平衡点が消える。薬の投与、個体群、振動、流出などの現象の予測と、仮定の検証に使う。
前提知識: 常微分方程式, 変数分離形の微分方程式, 平均値の定理
薬を点滴で入れ続けると、血液中の量はどこまで増えるのか。漁獲を毎年一定量ずつ続けると、魚の数は保たれるのか。どちらの問いも、現象をそのまま観察するだけでは答えにくい。そこで、量の変わり方についての仮定を式に書き、その式から結論を導く。このように現象を数学の対象に置き換えたものを数理モデルという。
数理モデルの価値は、仮定から結論までの道筋がすべて見えることにある。点滴の例では「薬が体から消える速さは、体内の量に比例する」と仮定すれば、量は一定の値に近づき、その値は点滴の速さに比例し、消える速さの係数に反比例することが証明できる。結論が観測と合わなければ、どの仮定が間違っていたかを調べられる。この記事では、数理モデルを定義し、収支の考え方から方程式を立てる方法、平衡点の安定性という最も基本的な解析の道具、パラメータの決め方を、証明と数値の例とともに述べる。
数理モデルとは、現象を表すために選んだ次のものの組である。
(i) 変数:時間とともに変わる量。現象の現在の様子を決める状態変数 $x$、外から与える入力 $u$、観測できる出力 $y$ に分ける。
(ii) パラメータ:考えている時間の範囲で一定とみなす量 $\theta$。
(iii) 関係式:変数とパラメータの間の方程式。動的なモデルでは $x'=f(x,u,\theta)$ のような微分方程式や、$x_{n+1}=g(x_n,u_n,\theta)$ のような漸化式、出力の式 $y=h(x,\theta)$。
(iv) 仮定と適用範囲:関係式を選んだ理由となる仮定と、それが成り立つとみなす変数・パラメータの範囲。
(iv) を定義に含めるのは、同じ方程式でも適用範囲が違えば別のモデルだからである。例えば $x'=kx$ は、資源が十分ある間の個体数のモデルとしては適切でも、何十世代も先の予測には使えない。状態変数とは、その時刻の値さえ分かれば、入力を与えたときの未来が決まる変数の組である。ばねにつながれた物体では位置だけでは足りず、位置と速度の組が状態になる。
数理モデルにはいくつかの分け方がある。代表的なものを表にまとめる。
| 分け方 | 一方 | 他方 | 例 |
|---|---|---|---|
| 時間の扱い | 連続時間(微分方程式) | 離散時間(漸化式) | 冷却則/年ごとの個体数 |
| 変化の有無 | 動的(時間変化を追う) | 静的(つり合いだけを見る) | ばねの振動/ばねの伸び |
| 偶然の扱い | 決定論的 | 確率的 | 冷却則/待ち行列 |
| 式の形 | 線形 | 非線形 | 点滴のモデル/収穫のあるロジスティックモデル |
| 由来 | 機構的(法則から導く) | 経験的(データの形から選ぶ) | 収支式/当てはめた多項式 |
機構的なモデルのパラメータには物理的な意味があり、外挿にも一定の根拠がある。経験的なモデルはデータへの当てはまりはよくても、データの範囲の外での予測には根拠がない。実際のモデルは多くの場合その中間で、収支式の骨組みに、経験的に選んだ比例の仮定を組み合わせて作られる。
多くの機構的なモデルは、「量の変化 = 入ってくる量 - 出ていく量」という収支の考え方から立てられる。体内の薬の量、容器の中の水の量、地域の人口、物体の熱量など、増減の原因を数え上げられる量であれば同じ形になる。まず、この考え方を微分方程式に移す根拠を確かめる。
$a,b$ を区間 $I$ 上の連続関数とする。区間 $I$ 上の関数 $x$ について、次の二つは同値である。
(i) 任意の $t_1< t_2$($t_1,t_2\in I$)で $x(t_2)-x(t_1)=\displaystyle\int_{t_1}^{t_2}\bigl(a(s)-b(s)\bigr)\,ds$。
(ii) $x$ は $I$ で微分可能で、$x'(t)=a(t)-b(t)$。
(i) ならば (ii):$t_0\in I$ を一つ固定すると、(i) から $x(t)=x(t_0)+\int_{t_0}^{t}(a-b)\,ds$ が任意の $t\in I$ で成り立つ($t< t_0$ のときは (i) を $t_1=t$、$t_2=t_0$ に使う)。被積分関数 $a-b$ は連続なので、微分積分学の基本定理により右辺は微分可能で、導関数は $a(t)-b(t)$ である。
(ii) ならば (i):$x'=a-b$ は連続なので $x$ は $C^1$ 級であり、微分積分学の基本定理により $x(t_2)-x(t_1)=\int_{t_1}^{t_2}x'(s)\,ds=\int_{t_1}^{t_2}(a-b)\,ds$。$\square$
収支そのものは (i) の形で観察されることが多い。「1 時間に入った量と出た量の差が、1 時間の増加量になる」という形である。命題により、流入・流出の速さが連続なら、これは微分方程式 (ii) と同じ情報をもつ。流入が途中で急に切り替わる場合は $a$ が不連続になり、(ii) は切り替えの時刻を除いて成り立つ。
収支式だけでは方程式は閉じない。出ていく速さ $b$ を状態でどう表すかを決める仮定が要る。これを構成則と呼ぶ。最も簡単な構成則は「出ていく速さは量に比例する」で、$b=kx$($k>0$)と書く。薬の消失、放射性物質の崩壊、穴の開いた容器からの細い流れ、預金の引き出しなど、多くの現象の第一近似になる。この仮定のもとで、外から速さ $u(t)$ で入れる場合の方程式は
$$
x'=-kx+u(t)
$$
になる。これを 1 コンパートメントモデルと呼ぶ。コンパートメントとは、中が十分よく混ざっていて一つの量で代表できる区画のことである。
$k>0$ とし、$x'=-kx+u(t)$ を考える。
(i) $u$ が一定値 $r\ge0$ のとき、$x(0)=x_0$ を満たす解は
$$
x(t)=\frac rk+\Bigl(x_0-\frac rk\Bigr)e^{-kt}
$$
であり、$t\to\infty$ で $x(t)\to r/k$。特に $x_0=0$ なら、$x(t)$ が定常値 $r/k$ の $90\%$ に達する時刻は $(\log10)/k$ である。
(ii) $u=0$ とし、時刻 $0,\tau,2\tau,\dots$ の直前の量に毎回 $D>0$ を一度に加える($x(0)=D$)。$n$ 回目の投与の直後の量を $P_n$ とすると
$$
P_n=D\,\frac{1-q^n}{1-q},\qquad q=e^{-k\tau},
$$
であり、$P_n$ は単調に増えて $D/(1-q)$ に収束する。次の投与の直前の量は $qP_n$ で、これは $qD/(1-q)$ に収束する。
(i) $v(t)=e^{kt}\bigl(x(t)-r/k\bigr)$ とおくと、
$$
v'=ke^{kt}\Bigl(x-\frac rk\Bigr)+e^{kt}(-kx+r)=0
$$
なので $v$ は定数 $v(0)=x_0-r/k$ に等しい。よって解は主張の式に限られ、逆にこの式は方程式と初期条件を満たす。$k>0$ から $e^{-kt}\to0$ なので $x(t)\to r/k$。$x_0=0$ のとき $x(t)=\frac rk(1-e^{-kt})$ で、$x(t)=0.9\,r/k$ は $e^{-kt}=0.1$、つまり $t=(\log10)/k$ と同値である。
(ii) 投与と投与の間では $u=0$ なので、(i) を $r=0$ で使えば量は $e^{-k(t-t_0)}$ 倍になる。よって次の投与の直前の量は直後の量の $q$ 倍であり、$P_{n+1}=qP_n+D$、$P_1=D$ が成り立つ。$P_n=D(1+q+\dots+q^{n-1})$ を帰納法で示す。$n=1$ では正しく、$n$ で成り立てば $P_{n+1}=D(q+\dots+q^n)+D$ となる。等比数列の和の公式から主張の式を得る。$0< q<1$ なので $q^n$ は単調に減って $0$ に近づき、$P_n$ は単調に増えて $D/(1-q)$ に収束する。直前の量 $qP_n$ の極限はその $q$ 倍である。$\square$
(i) の定常値 $r/k$ は、流入 $r$ と流出 $kx$ がつり合う量である。定常値に近づく速さは流入の速さによらず $k$ だけで決まる。点滴の速さを倍にすれば定常値は倍になるが、そこに近づくまでの時間は変わらない。(ii) の $1/(1-q)$ は蓄積係数と呼ばれ、繰り返しの投与で量が 1 回分の何倍まで積み上がるかを表す。
よく使われるモデルを、仮定と使える条件とともに並べる。根拠の列は、解や性質を証明している箇所を指す。
| モデル | 方程式 | 主な仮定 | 使える条件 | 根拠 |
|---|---|---|---|---|
| 指数モデル | $x'=kx$ | 増える速さが量に比例 | 資源が十分ある範囲 | 指数モデルとロジスティックモデル |
| ロジスティックモデル | $x'=rx(1-x/K)$ | 相対的な増え方が量とともに 1 次式で減る | $r,K>0$ | 指数モデルとロジスティックモデル |
| Newton の冷却則 | $T'=-k(T-T_a)$ | 熱の出入りが温度差に比例 | 室温 $T_a$ が一定 | 応用数学 |
| 1 コンパートメント | $x'=-kx+u$ | 中がよく混ざる、消失が量に比例 | $k>0$ | thm-mm-compartment |
| 一定量の収穫 | $x'=rx(1-x/K)-H$ | ロジスティックの増え方から毎年 $H$ を取る | $x>0$ の範囲 | prop-mm-harvest |
| ばねと質量 | $mx''+cx'+kx=u$ | 復元力が変位に、抵抗が速度に比例 | 変位が小さい範囲 | ex-mm-spring |
| 容器からの流出 | $h'=-c\sqrt h$ | 流れの速さが水深の平方根に比例 | $h\ge0$ | 計算(反例の表の 2 行目) |
表の下の 3 行のように、同じ収支の骨組みに異なる構成則を入れると異なるモデルになる。どの構成則を選ぶかは数学からは決まらず、現象についての知識と、データとの比較で決める。
動的なモデルで最初に調べるのは、時間がたっても変わらない状態、つまり平衡点である。$x'=f(x)$ の平衡点とは $f(x^*)=0$ を満たす $x^*$ で、定数関数 $x(t)=x^*$ は解になる。現実には状態は少しずつ乱されるので、平衡点から少しずれたときに戻るかどうかが重要になる。
$f$ を開区間 $J$ 上の $C^1$ 級関数とし、$x^*\in J$ で $f(x^*)=0$ とする。
(i) $f'(x^*)<0$ なら、ある $\delta>0$ が存在して、$\lvert x_0-x^*\rvert\le\delta$ を満たす任意の $x_0$ について、$x'=f(x)$、$x(0)=x_0$ の解は $t\ge0$ 全体で存在し、$c=-f'(x^*)/2$ として
$$
\lvert x(t)-x^*\rvert\le\lvert x_0-x^*\rvert\,e^{-ct}\qquad(t\ge0).
$$
(ii) $f'(x^*)>0$ なら、ある $\delta>0$ が存在して、$0<\lvert x_0-x^*\rvert\le\delta$ を満たす任意の $x_0$ から出発する解は、有限の時間のうちに区間 $[x^*-\delta,x^*+\delta]$ の外に出る。
$f$ は $C^1$ 級なので局所 Lipschitz 条件を満たし、常微分方程式 の定理「Picard–Lindelöf の定理」と系「解の一意性」により解は局所的にただ一つ存在する。
(i) $f'$ は連続なので、$\delta>0$ を十分小さくとれば、$[x^*-\delta,x^*+\delta]\subset J$ かつこの区間で $f'(x)\le f'(x^*)/2=-c$ となる。区間内の $x$ について、平均値の定理により $x$ と $x^*$ の間の $\xi$ で $f(x)=f(x)-f(x^*)=f'(\xi)(x-x^*)$ となるから
$$
(x-x^*)f(x)=f'(\xi)(x-x^*)^2\le-c\,(x-x^*)^2 .
$$
解 $x(t)$ に対して $V(t)=(x(t)-x^*)^2$ とおくと、解が区間内にある間は
$$
V'(t)=2(x-x^*)f(x)\le-2cV(t).
$$
まず解が区間 $I_\delta=[x^*-\delta,x^*+\delta]$ を出ないことを示す。$x_0\ne x^*$ としてよい($x_0=x^*$ なら解は定数で、主張は明らか)。このとき、定数解 $x^*$ も解であることと一意性から、どの時刻でも $x(t)\ne x^*$、つまり $V>0$ である。解の最大存在区間の右端を $\beta$ とし、$S$ を「$[0,t]$ の全体で $x(s)\in I_\delta$」となる $t\in[0,\beta)$ の集合とする。$S$ は $0$ を含む区間で、解の連続性から $[0,\beta)$ の中で閉じている。$S$ が $[0,\beta)$ の中で開いていることを示す。$t\in S$ とすると、$[0,t]$ で $V'\le-2cV<0$なので $V$ は減少し、$t>0$ なら $V(t)< V(0)\le\delta^2$ である。よって $V(t)<\delta^2$ なら、連続性から $t$ の少し先まで $V<\delta^2$ である。$V(t)=\delta^2$ となるのは $t=0$ かつ $V(0)=\delta^2$ の場合だけで、そのときは $V'(0)<0$ から、$0$ の少し先まで $V<\delta^2$ である。どちらの場合も、$t$ の少し先までの時刻は $S$ に属する。連結性により $S=[0,\beta)$ となり、解は $I_\delta$ にとどまる。コンパクト集合 $I_\delta\subset J$ にとどまる解は、常微分方程式 の定理「爆発の二者択一」により $\beta=\infty$、つまり $t\ge0$ 全体に延長できる。
次に $W(t)=e^{2ct}V(t)$ とおくと $W'=e^{2ct}(V'+2cV)\le0$ なので $W(t)\le W(0)=V(0)$、つまり $V(t)\le V(0)e^{-2ct}$。平方根をとれば主張の評価を得る。
(ii) $c=f'(x^*)/2>0$ とし、$\delta$ を $[x^*-\delta,x^*+\delta]\subset J$ かつこの区間で $f'\ge c$ となるようにとる。同じく平均値の定理から、区間内で $(x-x^*)f(x)\ge c(x-x^*)^2$。解が区間 $I_\delta$ にとどまり続けるとすると、(i) と同じく「爆発の二者択一」により解は $t\ge0$ 全体で存在する。任意の $T>0$ について $[0,T]$ で $V'=2(x-x^*)f(x)\ge2cV$ なので $e^{-2ct}V(t)$ は増加し、$V(t)\ge V(0)e^{2ct}$。$V(0)>0$ だから、$T$ が $\frac1{2c}\log\frac{\delta^2}{V(0)}$ を越えると $V(T)>\delta^2$ となり矛盾する。よって解はそれまでに区間の外に出る。$\square$
定理の (i) の状況の $x^*$ を漸近安定な平衡点、(ii) の状況の $x^*$ を不安定な平衡点という。判定に使うのは $f'(x^*)$ の符号だけで、解の式は要らない。これが線形化の考え方である。$f'(x^*)=0$ の場合は、この定理からは何も言えない(反例の表を見よ)。
この定理を使って、漁獲のある個体群のモデルを調べる。個体数 $x$ がロジスティックモデルで増え、毎年一定量 $H\ge0$ を取るとすると、モデルは
$$
x'=f(x)=rx\Bigl(1-\frac xK\Bigr)-H\qquad(r,K>0)
$$
になる。$H=0$ のときの平衡点 $K$ が、$H$ を増やすとどう動くかが問いである。
上のモデルで $H_{\max}=rK/4$ とおく。
(i) $0\le H< H_{\max}$ なら、平衡点は $x_\pm=\dfrac K2\Bigl(1\pm\sqrt{1-H/H_{\max}}\Bigr)$ の二つで、$x_+$ は漸近安定、$x_-$ は不安定である。
(ii) $H=H_{\max}$ なら、平衡点は $x=K/2$ の一つだけである。
(iii) $H>H_{\max}$ なら平衡点はなく、$x(0)>0$ から出発する解はすべて有限の時間で $0$ に達する。
$f(x)=0$ は $\frac rKx^2-rx+H=0$、つまり $x^2-Kx+\frac{KH}r=0$ と同値で、判別式は $K^2-\frac{4KH}r=K^2(1-H/H_{\max})$ である。これが正・零・負に応じて、平衡点は (i) の二つ・(ii) の一つ・なし、となる。
(i) $f'(x)=r\bigl(1-\frac{2x}K\bigr)$ に $x_\pm$ を代入すると $f'(x_\pm)=\mp r\sqrt{1-H/H_{\max}}$。$H< H_{\max}$ では根号の中は正なので、$f'(x_+)<0$、$f'(x_-)>0$。thm-mm-stability により $x_+$ は漸近安定、$x_-$ は不安定である。
(iii) $rx(1-x/K)$ は $x=K/2$ で最大値 $rK/4=H_{\max}$ をとるので、すべての $x$ で $f(x)\le H_{\max}-H=-\varepsilon<0$。解が $0$ より大きい間は $x'(t)\le-\varepsilon$ だから、$x(t)\le x(0)-\varepsilon t$ である。解が $0$ より大きい間は $x$ は減少して $(0,x(0)]$ にとどまるので、もし解が $0$ に達しないまま最大存在区間の右端 $\beta$ まで続くなら、解はコンパクト集合 $[0,x(0)]$ にとどまり、常微分方程式 の定理「爆発の二者択一」により $\beta=\infty$ となる。これは $x(t)\le x(0)-\varepsilon t$ と矛盾する。よって解は存在したまま、遅くとも $t=x(0)/\varepsilon$ までに $0$ に達する。$\square$
$H_{\max}=rK/4$ は最大持続収穫量と呼ばれる。$H$ がこの値を少しでも越えると、平衡点が消え、個体群は必ず絶える。$H$ が $H_{\max}$ より少し小さい場合も、安定な $x_+$ と不安定な $x_-$ が近づいているので、偶然の減少で $x_-$ を下回れば絶滅へ向かう。パラメータを連続的に変えたときに平衡点の数が変わる現象を分岐といい、この例は二つの平衡点が衝突して消える型である。
数理モデルを作って使う手順を、番号付きで述べる。
| 場合 | 平衡点 | 長い時間での振る舞い |
|---|---|---|
| $0\le H< rK/4$、$x(0)>x_-$ | $x_-$(不安定)、$x_+$(安定) | $x_+$ に近づく |
| $0\le H< rK/4$、$0< x(0)< x_-$ | 同上 | 有限の時間で $0$ に達する |
| $H=rK/4$ | $K/2$ のみ | $x(0)>K/2$ なら $K/2$ に近づき、$x(0)< K/2$ なら $0$ に達する |
| $H>rK/4$ | なし | 必ず $0$ に達する |
指数モデル $P'=kP$ に従う細菌が、時刻 $0$ に $300$ 個、5 時間後に $450$ 個だった。20 時間後の数と倍加時間を求める。
使う公式:$P(t)=P(0)e^{kt}$、観測からの $k$ の決定(指数モデルとロジスティックモデル の命題「倍加時間と半減期」)。
答え:$k=\frac15\log1.5\approx0.0810930$(毎時)、$P(20)=300\cdot1.5^4=1518.75$、倍加時間は $\frac{\log2}k\approx8.54756$ 時間。
手順 1:$450=300e^{5k}$ から $e^{5k}=1.5$。手順 2:$P(20)=300(e^{5k})^4=300\cdot5.0625$。手順 3:$e^{kT}=2$ から $T=\log2/k=5\log2/\log1.5$。
消失の係数 $k=0.2$(毎時)の薬を、速さ $r=10$ mg/時 で点滴する。体内の量 $x$ の定常値と、その $90\%$ に達する時刻を求める($x(0)=0$)。
使う公式:thm-mm-compartment の (i)。
答え:定常値 $r/k=50$ mg、$90\%$ に達する時刻は $(\log10)/0.2\approx11.5129$ 時間。
$x(t)=50(1-e^{-0.2t})$。$x=45$ は $e^{-0.2t}=0.1$ と同値。点滴の速さを $20$ mg/時 にしても、定常値が $100$ mg になるだけで、$90\%$ に達する時刻は同じ $11.5129$ 時間である。
同じ薬($k=0.2$)を $D=100$ mg ずつ $\tau=6$ 時間ごとに投与する。投与直後の量の極限と、投与直前の量の極限を求める。
使う公式:thm-mm-compartment の (ii)。
答え:$q=e^{-1.2}\approx0.301194$、蓄積係数 $1/(1-q)\approx1.43101$。直後の量は $143.101$ mg に、直前の量は $43.1013$ mg に近づく。
$P_1=100$、$P_2=100+100q\approx130.119$、$P_3\approx139.191$、$P_5\approx142.747$、$P_{10}\approx143.100$。直前の量は $qP_n$ で、極限は $100q/(1-q)$。
$r=1$、$K=1$ とし、$H=0.2$ と $H=0.3$ の場合を調べる。
使う公式:prop-mm-harvest、thm-mm-stability。
答え:$H=0.2$ では平衡点 $x_\pm=\frac{1\pm\sqrt{0.2}}2$、つまり $x_+\approx0.723607$($f'\approx-0.447214$、安定)と $x_-\approx0.276393$($f'\approx0.447214$、不安定)。$x(0)=0.25$ から出発すると約 $4.30409$ で $0$ に達する。$H=0.3>0.25$ では平衡点がなく、$x(0)=1$ からでも約 $10.2883$ で $0$ に達する。
$H=0.2$:$x^2-x+0.2=0$ の解。$x(0)=0.25< x_-$ なので減少し、$0$ に達する時刻は $\int_0^{0.25}\frac{dx}{(x-x_-)(x-x_+)}=\frac1{x_+-x_-}\Bigl[\log\frac{x_+-x}{x_--x}\Bigr]_0^{0.25}$。$H=0.3$:$x'=-\bigl((x-\tfrac12)^2+0.05\bigr)$ なので、時刻は $\int_0^1\frac{dx}{(x-1/2)^2+0.05}=\frac2{\sqrt{0.05}}\arctan\frac{0.5}{\sqrt{0.05}}$。
図: 収穫量と平衡点
図は、横軸に収穫量 $H$、縦軸に平衡点をとったものである($r=K=1$)。実線が安定な $x_+$、破線が不安定な $x_-$ で、$H=0.25$ で二つが出会って消える。矢印は解の動く向きを表す。
質量 $m=2$ kg の物体が、ばね定数 $k=50$ N/m のばねと、抵抗の係数 $c=4$ N·s/m の減衰器につながれている。Newton の運動方程式から状態方程式を作り、振動の周期を求める。
使う公式:力のつり合い $mx''=-kx-cx'$、状態 $(x,v)=(x,x')$。
答え:$\begin{pmatrix}x\\v\end{pmatrix}'=\begin{pmatrix}0&1\\-25&-2\end{pmatrix}\begin{pmatrix}x\\v\end{pmatrix}$。特性根は $-1\pm2\sqrt6\,i$、振動の周期は $\frac{2\pi}{2\sqrt6}\approx1.28255$ 秒。
手順 1:$v=x'$ とおくと $v'=-(k/m)x-(c/m)v=-25x-2v$。手順 2:特性方程式 $2s^2+4s+50=0$ の根は $s=-1\pm\sqrt{-24}$。手順 3:減衰の速さは実部 $-1$、振動の角振動数は虚部 $2\sqrt6\approx4.89898$。
この例で状態が位置だけでなく速度を含むのは、位置だけでは未来が決まらないからである。同じ位置にいても、動いている向きによってその後の運動は違う。2 階の方程式を 1 階の連立方程式に直す操作は、常微分方程式 の命題「高階の方程式と 1 階の連立方程式」で一般に扱われている。
| 外す条件 | 反例 | 成り立たなくなること |
|---|---|---|
| $f'(x^*)\ne0$ | $x'=x^3$ と $x'=-x^3$($x^*=0$) | どちらも $f'(0)=0$ だが、一方は不安定、他方は安定で、線形化では判定できない |
| $f$ が Lipschitz(一意性) | 容器からの流出 $h'=-c\sqrt h$ | 空の容器を見ても、いつ空になったかが決まらない |
| 適用範囲(資源が十分) | 指数モデルを長い時間に外挿する | 予測が現実の上限を大きく越える |
| 連続時間で考えること | ロジスティックモデルを大きな刻みの漸化式にする | 平衡点に単調に近づかず、振動や不規則な動きが出る |
1 行目の反例を確かめる。$x'=x^3$、$x(0)=x_0>0$ の解は $x(t)=x_0/\sqrt{1-2x_0^2t}$ で、$t=1/(2x_0^2)$ で発散する。$x_0=0.5$ なら $t=2$ である。一方 $x'=-x^3$ の解は $x(t)=x_0/\sqrt{1+2x_0^2t}$ で $0$ に近づき、$x_0=0.5$ なら $x(100)=1/\sqrt{204}\approx0.0700140$ である。近づき方は指数的でなく $t^{-1/2}$ 程度と遅い。$f'(0)=0$ の平衡点では、高次の項が振る舞いを決める。
2 行目の反例を確かめる。断面積 $A$ の容器の底に面積 $a$ の穴があり、水深 $h$ のとき流出の速さが $\sqrt{2gh}$ に比例するとすると、収支式は $Ah'=-a\sqrt{2gh}$、つまり $h'=-c\sqrt h$($c=\frac aA\sqrt{2g}$)になる。$h(0)=h_0$ の解は $\sqrt h=\sqrt{h_0}-\frac c2t$ で、時刻 $2\sqrt{h_0}/c$ に空になり、その後は $h=0$ が続く。$c=0.2$、$h_0=1$ なら $10$ 秒で空になり、$h(5)=0.25$ である。ところが、時刻 $10$ に空だったと観測しても、時刻 $T\le10$ に空になったどの解とも矛盾しない。$\sqrt h$ は $h=0$ で Lipschitz 条件を満たさず、過去へ向かう一意性が崩れているからである。モデルから過去を推定する使い方では、一意性の条件が欠かせない。
3 行目は ex-mm-growth で見られる。5 時間で $1.5$ 倍という率をそのまま 100 時間先へ延ばすと、$300\cdot1.5^{20}\approx9.98\times10^5$ 個になる。実際の培地には栄養と空間の上限があるので、どこかで増え方は鈍る。指数モデルの結論は、仮定「資源が十分」が成り立つ時間の範囲でだけ意味をもつ。
4 行目は、指数モデルとロジスティックモデル の例「$r$ を変えて計算する」で数値が示されている。漸化式 $b_{n+1}=b_n+rb_n(1-b_n)$ は、刻み $1$ の Euler 法をロジスティック方程式に使ったものと見ることもでき、同じ記事の命題「$0< r\le1$ なら単調に $1$ に近づく」の条件を外すと、$r=2.5$ では 4 つの値をくり返し、$r=2.7$ では不規則に動く。連続時間のモデルでは起こらない振る舞いであり、離散化の刻みがモデルの結論を変えてしまう例である。
| 誤り | 正しくは | なぜ |
|---|---|---|
| 平衡点があれば、そこに落ち着く | 安定かどうかを調べる | thm-mm-stability |
| $f'(x^*)=0$ なら安定 | 線形化では判定できない | 1 行目の反例 |
| 点滴の速さを上げれば早く定常に達する | 近づく速さは $k$ だけで決まる | thm-mm-compartment |
| 収穫を少し増やしても個体数は少し減るだけ | $rK/4$ を越えると必ず絶滅する | prop-mm-harvest |
| 当てはまったモデルは外挿にも使える | 適用範囲の外では根拠がない | 3 行目の反例 |
この記事では 1 次元のモデルを中心に扱った。2 つ以上の量が互いに影響し合うモデル、例えば捕食者と被食者の個体数のモデルや、感染症の広がりのモデルでは、平衡点の安定性は線形化した行列の固有値の実部で判定する。thm-mm-stability の (i) は、その 1 次元の場合にあたる。多次元では、平衡点に落ち着かず周期的に振動し続ける解(リミットサイクル)や、不規則な動きが現れることがある。
確率的なモデルでは、変数の値そのものではなく、その分布の時間変化を追う。個体数が少ないときの絶滅の確率のように、決定論的なモデルでは捉えられない現象を扱える。
数理モデルの考え方は、制御工学の教科書でも体系的に扱われている。状態・入力・出力の分け方や、モデルが近似であることへの注意は AMF20 §3.1(pp. 3-1–3-2)に、モデルを作る方法論は同 §3.3(p. 3-22–)にある。現実の問題から方程式を立てて解き、結果を解釈する往復は LebDQ26 §0.2.3(p. 12)でも述べられている。
Mathpediaは寄付と、参考文献の書籍リンク(Amazonアソシエイト)の紹介料で運営されています。 支援について / 寄付する