このサイトは学部では早稲田で物理を, 修士では東大で数学を専攻し, 今も非アカデミックの立場で数学や物理と向き合っている一市民の奮闘の記録です. 運営者情報および運営理念についてはこちらをご覧ください.

理系のための総合語学・リベラルアーツの視点から数学・物理・プログラミング・語学 (特に英語) の情報を発信しています. コンテンツアーカイブに見やすくまとめているのでぜひご覧ください.

こちらに PDF が置いてあります.
サイトでは見づらい方は PDF をご覧ください.

マルサスの『人口論』と微分方程式

あなたが高校生以上なら社会の授業で次のような言葉を聞いたことがあるかもしれません.

人口は制限されなければ幾何級数的に増加するが生活資源は算術級数的にしか増加しないので生活資源は必ず不足する.

これは有名なマルサスの『人口論』で議論された話です.
幾何級数はいわゆる等比数列の和で,
算術級数は等差数列の和です.
マルサスは経済学者なので経済の話です.
経済の話でこういう数学ネタが出てくるわけですね.

ちなみに大学の経済学だと経済の話で微分積分が本当に出てきます.
ときどき経済学部の入試で数学が受験科目に入っていることがあります.
実際に大学に入ってから使うからです.
それも高校の理系のレベルを越えた数学です.
統計データをいろいろいじらないといけないので統計学が必要で,
そっちでも割といろいろ数学が必要です.

話を元に戻しましょう.
マルサスの人口論に合わせて人口の増え方を考えます.
ここではもっと一般的に生物の集団だと思いましょう.
まずは生物の種類が 1 種類だとして考えていきます.
生物っぽく人の数というよりも一般に個体数と呼ぶことにして,
時刻 $t$ での個体数を $x(t)$ とします¹.
個体数の増加率を出生率の死亡率の差として定義します.
微分方程式的にいうとこれは $x'(t) / x(t)$ です.
増加率が一定値 $\alpha$ に等しいとすると $x'(t) / x(t) = \alpha$ だから次の微分方程式が出てきます.

\begin{align}
\frac{d x(t)}{dt}
=
\alpha x(t).
\end{align}

この微分方程式をマルサスモデルと言います.
これを解くと $x(t) = x(0) e^{\alpha t}$ です.
これで指数関数が出てきました.
等比数列は一般項が $a_n = 2^{n}$ のように指数関数で書けている数列だから,
これを足し上げていけばまさに幾何級数になりますね.
放射性物質の崩壊とは指数の肩の符号が変わっているだけで,
それ以外は同じ式です.
放射性物質の崩壊と生物の個体数の変化を追う方程式が同じ形をしているわけです.

シミュレーション

シミュレーションの結果は放射性物質の崩壊のときと基本的に同じです.
いちおうきちんとやっておきましょうか.
まず微分方程式を次のように近似します.

\begin{align}\label{math_refuge00005}
\frac{x_{n+1} – x_{n}}{h}
=
\alpha x_n, \quad
x_{n+1}
=
\alpha x_n + h x_n.
\end{align}

これにしたがって解いてみましょう.
プログラムと結果の図はこのページの後半を確認してください.

少し現実的な設定にしてみる

もう少し現実的な設定にしてみます.
個体数 $x(t)$ が大きくなると食物が足りなかったり,
衛生状態の悪化で病気にかかりやすくなったりして増加率が減るはずです.
そこで増加率 $\alpha$ を $\alpha – \beta x$ ($\alpha$, $\beta > 0$) に変えてみます.
\begin{align}
\frac{dx}{dt}
=
(\alpha – \beta x) x.
\end{align}
これにはロジスティック方程式という名前がついています.
悪化を $\beta x$ と書くことに必然性はありません.
とりあえずやってみただけです.

ロジスティック方程式の厳密解

解き方はさておき解は次のようになります.
あなたがもしこの関数を微分できるなら,
厳密解になることを計算して確認してみてください.
\begin{align}
x(t)
=
\frac{\alpha}{\beta} \frac{1}{1 + e^{- \alpha t}}.
\end{align}

これを近似すると次のようになります.
\begin{align}
\frac{x_{n+1} – x_{n}}{h}
=
(\alpha – \beta x_{n}) x_{n}, \quad
x_{n+1}
=
x_n + h (\alpha – \beta x_{n}) x_{n}.
\end{align}

これも厳密解と比較しながら計算してみます.
プログラムや図はこの記事の後半を参照してください.

近似の精度が気になります

これもよく近似できていることは認めてもらえるんじゃないでしょうか?
でもうるさい人は「近似の精度はどう決めるの?」みたいなことを考えているはずです.

もちろんいたって真っ当な指摘です.
ただけっこう難しい問題がたくさんあって,
なかなかスカっと綺麗に回答できません.
あとで少し考えることにしてここではこのまま進みます.

生物学として正しいの?

あなたは数値計算の精度とかそれ以前の問題として,
次のように思っているかもしれません.

• この微分方程式で実際の生物の増減にあてはめられるの?
• もっと現実的なことを考える必要があるんじゃないか?

もちろんこれも正しい指摘です.

例えばマルサスモデルでもロジスティック方程式でも右辺が気になります.
実際には成長しないと子作りできないですからね.
その成長が必要だという情報が方程式に盛り込まれていません.

生物で考えるなら食物連鎖があるわけで食べられて個体数が減ることだってあります.
天敵が増えたらその個体の数はあおりを受けて激減するはずです.

もちろんこんなことは折り込み済みで,
例えば遅延型微分方程式,
ロトカ-ボルテラ方程式と呼ばれる方程式が対応しています.
これを調べるのも面白いんですが,
今回はこのくらいにしておきましょう.

今回のポイント

簡単にまとめると,
今回は近似の精度が気になるあなたのために,
実際どのくらいよく本当の解を近似できているのかを調べました.
経済学や生物のようにあんまり数学との関係がなさそうな分野と数学の関係も紹介しています.

数学的ポイント: 数列と漸化式

長くなってきましたが最後にちょっと大事な話.
式 \eqref{math_refuge00005} で $\alpha = 1$, $h = 1$ とすると $x_{n+1} = 2 x_n$ になります.
この漸化式は高校でも出てくるやつで, もちろん公比 $2$ の等比数列の漸化式です.
底が $2$ か $e$ かの違いはあっても指数関数で書けることは同じです.
これはたまたまではなくて数学的に意味があることです.

$x_{n+1} – x_{n}$ のように数列の隣の項の差を差分と言うことがあります.
もちろん微分との関係を強く意識した言葉です.
数列の問題, もっと言えば漸化式の問題は微分方程式とも深い関係があります².
あまり深入りはしませんが数学や物理に興味があるなら覚えておくといろいろ楽しいことがあります.

アンケート

今回もアンケートがあります.
改善に繋げたいのでぜひ積極的に回答をお願いします.

• https://goo.gl/forms/SVa4aNyaqOetTzYx2

ではまた次回をお楽しみに!

プログラミングパート

マルサスモデル

基本は放射性物質の崩壊と同じだからあっさり行こう.

\begin{align}
\frac{dx}{dt}
=
\alpha u.
\end{align}
厳密解は $x(t) = C_0 e^{\alpha t}$ だ.
初期値を設定すれば $C_0$ はそこから決まる.

微分を単純に離散化すると次のようになる.

\begin{align}
\frac{x_{n+1} – x_{n}}{h}
=
\alpha x_{n}.
\end{align}

整理すると次の通り.

\begin{align}
x_{n+1}
=
x_{n} + \alpha h x_{n}.
\end{align}

オイラー法でコードに落とそう.

マルサスモデルの近似解プログラム

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

def malthus_euler(nt, init = 10):
    dt = 2 / (nt - 1)
    # 初期条件設定
    x = np.zeros(nt)
    x[0] = init

    for i in range(1, nt):
        x[i] = x[i-1] + c * dt * x[i-1]

    return x


c = 1
init = 1
nt = 101
x_approx = malthus_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)


t = np.linspace(0, 2, nt)
x_exact = init * np.exp(c * t)
plt.plot(t, x_exact)


plt.legend(['approximation', 'exact'])

plt.show()

RESULT

放射性物質の崩壊の時との違い

t が大きくなると厳密解の方が少し大きくなる.
指数の肩の符号が違うだけなので,
放射性物質の崩壊で形式的に時間を負の方に伸ばせば同じようにズレが出てくるはずだ.
私の今の感覚だとこのズレが今回のように大きくなるように出るのかどうかまではわからない.

何はともあれ区間の分割数 nt を大きくしてみたのが次の結果:
具体的には 101 から 1001 にした.

パラメータを修正して再計算

ライブラリが読み込まれていたり関数 malthus_euler が定義されている前提のコード片なので,
これだけで実行できるわけではない.
IPython か Jupyter なら順にセルを読み込んでいけばそのまま実行できる.

c = 1
init = 1
nt = 1001
x_approx = malthus_euler(nt, init)
plt.plot(np.linspace(0, 2, nt), x_approx)


t = np.linspace(0, 2, nt)
x_exact = init * np.exp(c * t)
plt.plot(t, x_exact)


plt.legend(['approximation', 'exact'])

plt.show()

RESULT

注意

$nt = 101$ よりも近似の精度が良くなった.
この辺は単純に振る舞ってくれるようだ.

調和振動子だとオイラー法自体がうまく動かないから注意したい.

ロジスティック方程式

方程式は次の通り.

\begin{align}
\frac{dx}{dt}
=
(\alpha – \beta x) x.
\end{align}

厳密解は次の通り.

\begin{align}
x(t)
=
\frac{\alpha}{\beta} \frac{1}{1 + e^{- \alpha t}}.
\end{align}

近似すると次の通り.

\begin{align}
\frac{x_{n+1} – x_{n}}{h}
=
(\alpha – \beta x_{n}) x_{n}, \quad
x_{n+1}
=
x_n + h (\alpha – \beta x_{n}) x_{n}.
\end{align}

では数値的に解いてみよう.
とりあえず $\alpha = 2$, $\beta = 1$ で計算している.

ロジスティック方程式の近似解プログラム

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

def logistics_euler(nt, init = 10):
    dt = T / (nt - 1)
    # 初期条件設定
    x = np.zeros(nt)
    x[0] = init

    # ベクトル計算で一気に計算したい
    for i in range(1, nt):
        x[i] = x[i-1] + dt * (alpha - beta * x[i-1]) * x[i-1]

    return x


alpha = 2
beta = 1
init = 1
nt = 101
T = 5 # 時間変化を見る最大値
x_approx = logistics_euler(nt, init)
plt.plot(np.linspace(0, T, nt), x_approx)


t = np.linspace(0, T, nt)
x_exact = (alpha / beta) / (1 + np.exp(- alpha * t))
plt.plot(t, x_exact)


plt.legend(['approximation', 'exact'])

plt.show()

RESULT

時間の刻みを変えてみる

近似解と厳密解にちょっとズレが見える.
そこで nt = 101 から nt = 1001 にしてみた.

パラメータを修正して再計算

alpha = 2
beta = 1
init = 1
nt = 1001
T = 5 # 時間変化を見る最大値
x_approx = logistics_euler(nt, init)
plt.plot(np.linspace(0, T, nt), x_approx)


t = np.linspace(0, T, nt)
x_exact = (alpha / beta) / (1 + np.exp(- alpha * t))
plt.plot(t, x_exact)


plt.legend(['approximation', 'exact'])

plt.show()

RESULT