マルサスの人口論

アップデート 2010.4.13

マルサスの人口論

　イギリスの経済学者であり数学者でもあったマルサス(Malthus)は、一定の地域に住む人間の、時刻tにおける人口x(t)について考察した。出生数bが人口に比例すると仮定すると、出生率はb/xになり、死亡数dが人口に比例すると仮定すると、死亡率はd/xとなる。従って人口の成長率aは
a＝b/x－d/x＝(b－d)/x
となるというモデルを考えた。このとき人口xに対する微分方程式は、導関数dx/dtをx'(t)と書くことにすると、
x'＝ax
で与えられることとなる。この増殖率aをマルサス係数(Malthusian parameter)と言う。この微分方程式をマルサスの法則という。各点で微分係数x'が計算できる。微分係数x'は解曲線の接線の傾きである。各点での傾きx'を図示すると、次のベクトル場が得られる。傾きを表す直線に沿って解曲線が描かれる。微分方程式の解曲線の概形はこのベクトル場のグラフから推測できる。

　一般に
x'＝f(t)g(x)
の様にtだけの関数f(t)と、xだけの関数g(x)との積の形で、導関数x'(t)が与えられた微分方程式を、変数分離形という。g(x)が0でないとき、両辺をg(x)で割ると
dx/g(x)＝f(t)dt
となる。この式を積分すると
∫(1/g(x))dx＝∫f(t)dt
となり、変数分離形の微分方程式を解くことは、上式を積分をする問題に帰着される。
　マルサスの成長の方程式
x'＝ax
の場合
∫(1/x)dx＝∫adt
であるから、積分すると
log|x|＝at＋C₁
となる。両辺の指数をとると
|x|＝exp(at＋C₁)＝C₂e^at　ただし　C₂＝exp(C₁)
となり、絶対値をはずすと
x＝Ce^at　ただし　C＝±C₂
と表現できる。
これがマルサスの微分方程式の一般解である。定数Cを微分方程式の任意定数という。
x₀＝x(t₀)
を初期値とする解は
x(t)＝x₀e^atとなる。x₀は初期人口であるからx₀＞0である。また人口は増加すると考えているからa＞0である。したがってx＞0である。ｅ＝2.718…をネピアの定数または自然対数の底という。解曲線は指数関数である。時間の経過と共に急激に増大する。
　このマルサスの方程式をルンゲ_クッタ法で解いてみる。
　Excelのシートの内容は次の様になっている。(クリックすると別のウインドウに説明が表示される。)

A B C D E F

1
成長曲線：x'＝ax

2
変数t₀ 0

f(t,x) ＝B5D2

3
関数 x₀ 3 変数 t 関数 x 計算式メモ　'＝B5D2

4
刻み幅h 0.1

5
係数a 1

ここで数値を代入する箇所は

B2 変数tの初期値0を代入している。

B3 関数xの初期値3を代入している。

B4 刻み幅hに0.1代入している。

B5 微分方程式の係数aに１を代入している。

である。この数値aを変更することにより成長率aの異なる方程式を解いたり、刻み幅hを変えることにより解の精度を変更したり、初期条件x₀を変更することができる。
　また数式を代入する箇所はF2である。

F2 ＝B5D2
＝は、この欄は数式を扱うという記号である。
B5は、マルサスの係数aに与えられた数値1を参照している。
D2は、関数xの一時保管場所にある未知の値を参照している。
＝B5D2は、結局axを計算していることになる。
半角英数文字を使用しなければならない。

　このF2の計算式を変更することにより、x'(t)＝f(t,x)という形の1階常微分方程式を解くことができる。
　その他の欄は、計算には関係のない欄である。ただし、どんな計算をしているかが後で理解できるように、変数名とか、計算式とかを記入しておく。
　C4とD39とで囲まれる長方形の中に、計算結果が貯えられる。
　微分方程式の数値解は、表計算でもできるが、ここではルンゲ_クッタ法をマクロで記述する。そうすることにより、EXCELのシートには途中結果が表示されず、管理や保守が簡単になる。

マクロの作成方法の説明

ルンゲ_クッタ法のマクロについて、その機能の説明を以下に示す。

Sub RK_method() マクロの名前

Range("c4:d39").Clear C4からD39までの領域を空白にする

Cells(4, 3) ＝ Cells(2, 2) C4にB2(すなわちtの初期値)を代入する

Cells(4, 4) ＝ Cells(3, 2) D4にB3(すなわちxの初期値)を代入する

h ＝ Cells(4, 2) 刻み幅hをB4から読み取り、記憶する

h2 ＝ h / 2 刻み幅の半分を計算し、記憶する

For n ＝ 4 To 38 4から38までくり返す

tn ＝ Cells(n, 3) t_nにCn(すなわちtの値)を代入する

xn ＝ Cells(n, 4) x_nにDn(すなわちxの値)を代入する

Cells(2, 3) ＝ tn 作業領域C2にt_nを代入する

Cells(2, 4) ＝ xn 作業領域D2にx_nを代入する

k1 ＝ h * Cells(2, 6) 計算結果がF2にあるので、k₁を計算し、記憶する

Cells(2, 3) ＝ tn + h2 作業領域C2に、新しいtを代入する

Cells(2, 4) ＝ xn + k1 / 2 作業領域D2に、新しいxを代入する

k2 ＝ h * Cells(2, 6) 計算結果がF2にあるので、k₂を計算し、記憶する

Cells(2, 4) ＝ xn + k2 / 2 作業領域D2に、新しいxを代入する

k3 ＝ h * Cells(2, 6) 計算結果がF2にあるので、k₃を計算し、記憶する

tnp ＝ tn + h t_n+1を計算する

Cells(2, 3) ＝ tnp 作業領域C2に、新しいtを代入する

Cells(2, 4) ＝ xn + k3 作業領域D2に、新しいxを代入する

k4 ＝ h * Cells(2, 6) 計算結果がF2にあるので、k₄を計算し、記憶する

Cells(n + 1, 3) ＝ tnp 新しいtをCnpに代入する

Cells(n+1,4)＝xn+(k1+2k2+2k3+k4)/6 新しいDをCnpに代入する

Next ループの終わり

End Sub マクロの終わり

　EXCELの基本操作3で説明した、メニューバー：「表示」→「ツールバー」→「フォーム」→「コマンドボタン」の操作では、Sub RK_method()としたマクロ名が「ボタン1.Click()」となる。

　完成したグラフを次に表示する。グラフは散布図を使って描く。

　個体数が2倍になるまでの時間Tを求めてみよう。x(T)＝2x₀であるから
　2x₀＝x₀e^aT
を解けばよい。すると
　log2＝aTであるから
　T＝(log2)/a
となる。成長率が年率2％であるとき、EXCELで　＝LN(2)/0.02　と打ち込んで計算すると、34.6573590279973年になることが分かる。
　EXCELでは、自然対数はLN、常用対数はLOG10、任意の底を用いる対数はLOGである。
　例えば＝LN(EXP(2))は、指数関数と自然対数関数が逆関数であることから２となる。
　また＝LOG10(10000)は10の4乗の常用対数であるから４となる。
　さらに＝LOG(8,2)は８が２の3乗であることから、log₂8＝log₂2³＝3となる。

　元金がx₀円である預金に対して、年利Rの利息が付くとする。このとき1年後の元利合計はx(1)＝x₀＋x₀･Ｒ＝x₀(1+R)円となる。
　2年後の元利合計はx(2)＝x₀(1+R)＋x₀(1+R)･R＝x₀(1+R)²となる。
　n年後の元利合計はx(n)＝x₀(1+R)ⁿとなる。
　この計算は単利の場合であるが、年2回の複利（1回の利率は半分になるが、年2回利息が付く）の場合、n年後の元利合計はx₂(n)＝x₀(1+R/2)²ⁿとなる。
　年k回の複利の場合は、n年後の元利合計はx_k(n)＝x₀(1+R/k)^knとなる。
　複利の回数を無限回に増やした場合n年後の元利合計は
limx_k(n)＝x₀lim(1+R/k)^kn＝x₀e^Rn＝x₀e^Reⁿ
k→∞　　　k→∞
となる。この値は微分方程式
x'(t)＝a･x(t)　ただし　a＝R
x(0)＝x₀
の解である。このときの1年後の元利合計はx₀e^R≒x₀(1+R+R²/2)＞x₀(1＋R)であり、離散型の真の解より少し大きくなる。

預金利息の問題は、微分方程式で近似することができる。

問題1：lim(1+1/x)^x＝eであることを用いて、R>0のとき
　　　　x→∞
　　　　lim(1+R/k)^kn＝e^Rnであることを証明せよ。
　　　　k→∞
問題2：年利2％の1年複利と年利2％の半年複利(半年の利率は1％)の預金では、現金100万円に対して、3年後にはどの位の差ができるか。
問題3：年率6.33％の利息がつくとき、複利で運用すれば、2倍になるのには何年かかるか。(平成2年の定額貯金利息)
問題4：年率0.5％の利息がつくとき、複利で運用すれば、2倍になるのには何年かかるか。(2000年3月の普通預金利息)
問題5：2000年3月発行の国債の利率は年当たり1.9％の利息がつく。これを複利で運用したとすると、2倍になるのにはほぼ何年かかるか。
問題6：初期値x₀を変えると、解曲線のグラフはどのように変わるか。
問題7：マルサスの係数aを変えると、解曲線のグラフはどのように変わるか。図と文章で説明せよ。
問題8：さびによる腐食の方程式はx'(t)=ae^-btx(t),a>0,b>0 である。解曲線を図示せよ。また、厳密解を求めよ。
問題9：自分で問題を考えて、自分で解け。

	A	B	C	D	E	F
1	成長曲線：x'＝ax
2	変数t₀	0			f(t,x)	＝B5*D2
3	関数 x₀	3	変数 t	関数 x	計算式	メモ　'＝B5*D2
4	刻み幅h	0.1
5	係数a	1

B2	変数tの初期値0を代入している。
B3	関数xの初期値3を代入している。
B4	刻み幅hに0.1代入している。
B5	微分方程式の係数aに１を代入している。

Sub RK_method()	マクロの名前
Range("c4:d39").Clear	C4からD39までの領域を空白にする
Cells(4, 3) ＝ Cells(2, 2)	C4にB2(すなわちtの初期値)を代入する
Cells(4, 4) ＝ Cells(3, 2)	D4にB3(すなわちxの初期値)を代入する
h ＝ Cells(4, 2)	刻み幅hをB4から読み取り、記憶する
h2 ＝ h / 2	刻み幅の半分を計算し、記憶する
For n ＝ 4 To 38	4から38までくり返す
tn ＝ Cells(n, 3)	t_nにCn(すなわちtの値)を代入する
xn ＝ Cells(n, 4)	x_nにDn(すなわちxの値)を代入する
Cells(2, 3) ＝ tn	作業領域C2にt_nを代入する
Cells(2, 4) ＝ xn	作業領域D2にx_nを代入する
k1 ＝ h * Cells(2, 6)	計算結果がF2にあるので、k₁を計算し、記憶する
Cells(2, 3) ＝ tn + h2	作業領域C2に、新しいtを代入する
Cells(2, 4) ＝ xn + k1 / 2	作業領域D2に、新しいxを代入する
k2 ＝ h * Cells(2, 6)	計算結果がF2にあるので、k₂を計算し、記憶する
Cells(2, 4) ＝ xn + k2 / 2	作業領域D2に、新しいxを代入する
k3 ＝ h * Cells(2, 6)	計算結果がF2にあるので、k₃を計算し、記憶する
tnp ＝ tn + h	t_n+1を計算する
Cells(2, 3) ＝ tnp	作業領域C2に、新しいtを代入する
Cells(2, 4) ＝ xn + k3	作業領域D2に、新しいxを代入する
k4 ＝ h * Cells(2, 6)	計算結果がF2にあるので、k₄を計算し、記憶する
Cells(n + 1, 3) ＝ tnp	新しいtをCnpに代入する
Cells(n+1,4)＝xn+(k1+2k2+2k3+k4)/6	新しいDをCnpに代入する
Next	ループの終わり
End Sub	マクロの終わり

マルサスの人口論

A B C D E F 1 成長曲線：x'＝ax 2 変数t0 0 f(t,x) ＝B5*D2 3 関数x0 3 変数 t 関数 x 計算式 メモ '＝B5*D2 4 刻み幅h 0.1 5 係数a 1

ここで数値を代入する箇所は

B2 変数tの初期値0を代入している。 B3 関数xの初期値3を代入している。 B4 刻み幅hに0.1代入している。 B5 微分方程式の係数aに１を代入している。

である。この数値aを変更することにより成長率aの異なる方程式を解いたり、刻み幅hを変えることにより解の精度を変更したり、初期条件x0を変更することができる。 また数式を代入する箇所はF2である。

マクロの作成方法の説明

ルンゲ_クッタ法のマクロについて、その機能の説明を以下に示す。

完成したグラフを次に表示する。グラフは散布図を使って描く。

A B C D E F

1
成長曲線：x'＝ax

2
変数t₀ 0

f(t,x) ＝B5D2

3
関数 x₀ 3 変数 t 関数 x 計算式メモ　'＝B5D2

4
刻み幅h 0.1

5
係数a 1

B2 変数tの初期値0を代入している。

B3 関数xの初期値3を代入している。

B4 刻み幅hに0.1代入している。

B5 微分方程式の係数aに１を代入している。

である。この数値aを変更することにより成長率aの異なる方程式を解いたり、刻み幅hを変えることにより解の精度を変更したり、初期条件x₀を変更することができる。
　また数式を代入する箇所はF2である。

　完成したグラフを次に表示する。グラフは散布図を使って描く。