数理モデル

同義語:mathematical model数理モデリングmathematical modeling

概要

数理モデル(mathematical model)とは、現象を表すために選んだ状態・入力・出力の変数、パラメータ、それらを結ぶ方程式、仮定と適用範囲の組である。量の出入りの収支式に構成則を加えると微分方程式が得られ、例えば 1 コンパートメントモデル $x'=-kx+r$ の解は定常値 $r/k$ に近づく。$x'=f(x)$ の平衡点 $x^*$ は $f'(x^*)<0$ なら漸近安定、$f'(x^*)>0$ なら不安定であり、一定量の収穫があるロジスティックモデルでは収穫が $rK/4$ を越えると平衡点が消える。薬の投与、個体群、振動、流出などの現象の予測と、仮定の検証に使う。

$$\newcommand{C}[0]{\mathbb{C}} \newcommand{div}[0]{\mathbin{÷}} \newcommand{N}[0]{\mathbb{N}} \newcommand{Q}[0]{\mathbb{Q}} \newcommand{R}[0]{\mathbb{R}} \newcommand{Z}[0]{\mathbb{Z}} $$

前提知識: 常微分方程式, 変数分離形の微分方程式, 平均値の定理

薬を点滴で入れ続けると、血液中の量はどこまで増えるのか。漁獲を毎年一定量ずつ続けると、魚の数は保たれるのか。どちらの問いも、現象をそのまま観察するだけでは答えにくい。そこで、量の変わり方についての仮定を式に書き、その式から結論を導く。このように現象を数学の対象に置き換えたものを数理モデルという。
数理モデルの価値は、仮定から結論までの道筋がすべて見えることにある。点滴の例では「薬が体から消える速さは、体内の量に比例する」と仮定すれば、量は一定の値に近づき、その値は点滴の速さに比例し、消える速さの係数に反比例することが証明できる。結論が観測と合わなければ、どの仮定が間違っていたかを調べられる。この記事では、数理モデルを定義し、収支の考え方から方程式を立てる方法、平衡点の安定性という最も基本的な解析の道具、パラメータの決め方を、証明と数値の例とともに述べる。

定義と構成要素

数理モデル

数理モデルとは、現象を表すために選んだ次のものの組である。
(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 コンパートメントモデルと呼ぶ。コンパートメントとは、中が十分よく混ざっていて一つの量で代表できる区画のことである。

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_-$ を下回れば絶滅へ向かう。パラメータを連続的に変えたときに平衡点の数が変わる現象を分岐といい、この例は二つの平衡点が衝突して消える型である。

計算の手順

数理モデルを作って使う手順を、番号付きで述べる。

  1. 問いを決める。 何を予測・説明したいかを一文で書く。問いによって状態変数と適用範囲が決まる。
  2. 変数とパラメータを列挙する。 状態・入力・出力・パラメータを分け、それぞれの単位を書く。次元解析 で、関係しうる量の組み合わせを先に絞れることがある。
  3. 収支式を立てる。 量の出入りを数え上げ、prop-mm-balance の形で書く。
  4. 構成則を選ぶ。 流出・反応・力を状態で表す仮定を選ぶ。最初は比例のような最も簡単な仮定から始める。
  5. 解析する。 解の存在と一意性、平衡点と安定性(thm-mm-stability)、パラメータを変えたときの分岐を調べる。閉じた形で解けるなら解く。
  6. パラメータを決める。 観測からパラメータを求める。観測の数がパラメータの数より少ないと、パラメータは一つに決まらない。
  7. 検証し、直す。 予測を、パラメータを決めるのに使っていないデータと比べる。合わなければ手順 4 に戻り、構成則を取り替える。
    手順 5 の分岐は、表で場合分けすると見通しがよい。収穫のモデルでは次のとおりである。
    場合平衡点長い時間での振る舞い
    $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$ に達する
    表の 1 行目と 3 行目の「近づく」は、次の理由で成り立つ。1 次元の自励系では、解は平衡点を通り越せない(平衡点も解なので、一意性から二つの解のグラフは交わらない)。そのため解は単調で、有界なら極限をもち、極限は平衡点でなければならない。極限 $\ell$ が平衡点でなければ $x'\to f(\ell)\ne0$ となり、解が有界であることに反するからである。2 行目と 4 行目の「有限の時間で $0$ に達する」は、prop-mm-harvest の (iii) と同じく、$x'$ が負の定数以下に抑えられることから従う。

例題

細菌の増加

指数モデル $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アソシエイト)の紹介料で運営されています。 支援について / 寄付する