ラムゼーモデルの位相図が出たよ！（出たけどさあ……）

あらすじ

ラムゼーモデルを理解したニートは、ラムゼーモデルの収束の様子を図にすべく
講義ページのプログラムを動かそうとした。
だがそこに待ち受けていたのは数々の罠だった。

図

結論から言うと図は出せた。

↑これね。教科書と同じやつ。
$K_1 = 0.5 K_*$ から出発して定常状態に落ち着く様子と、
$K_1 = 2 K_*$ から出発して定常状態に落ち着く様子。

どうやったか

結局nleqslvが動かない問題は解決せず。
講義のページで紹介されていたOctaveというソフトをインストールして、
そっちのソフトで非線形連立方程式を解いた。
R言語とmatlab言語（？）で記法が違うので大変だった。

次回以降もR言語のプログラムをmatlab言語に書き直すのは手間がかかる（げんなり

参考

教科書と同じ図が出ているので、これでたぶんmatlabで非線形連立方程式が解けてる、はず。
ググりながら見様見真似でやったので、動いてはいるけどたぶん書き方とかメチャクチャなはず。
プログラムとかに詳しい人いたら、正しいやり方教えて下さい。

global a
global b
global d
global At
global Kt
global Kss
global Css
global T
a = 0.3
b = 0.99
d = 0.25
At = 1.0

Kt = ((1/b+d-1)/a/At)^(1/(a-1))
function res = FCt(K)
	res = At*K^a - d*K
end

Kss = ((1/b+d-1)/a/At)^(1/(a-1))
Css = At*Kss^a - d*Kss

function ret = objfun0(X, n, X0, Xss)
global a
global b
global d
global At
global Kt
global Kss
global Css
global T
	C = ones(n+1, 1);
	K = ones(n+1, 1);
	C(1:n) = X(1:n);
	C(n+1) = Xss;
	K(1)   = X0;
	K(2:n+1) = X(n+1 : 2*n);
	
	ret = ones(n, 2);
	for t = 1:n
		ret(t, 1) = C(t+1)/C(t)-b*(a*At*K(t+1)^(a-1)-d+1);
		ret(t, 2) = K(t+1)-At*K(t)^a-(1-d)*K(t)+C(t);
	end
end
function ret = obj1(X)
	global Kss
	global Css
	global T
	ret = objfun0(X, T, Kss*2, Css);
end

T = 30
Xinit = ones(2*T, 1);
X = fsolve(@obj1, Xinit)

これの出力をテキストファイルに保存して、以下のやり方でR言語に読み込ませた。
講義ページに書いてあったプログラムは、

rslt1 <- nleqslv(Xinit, objfun1)
rslt1x <- matrix(rslt1$x, maxT, 2)

の部分を消して、

rslt1 <- scan("sim1.txt")
rslt1x <- matrix(rslt1, maxT, 2)

としたら動いた。

ぐぐったもの

このあたり。