学生向け
微分方程式を数値的に簡単に解く方法を教えてください。
と訊かれた。Runge-Kutta法を実装すれば勉強になると思いつつも、(1)モデル化して、(2)数値的に解いて、(3)検討する、というプロセスの(2)のウェイトが大きすぎるのはかわいそうだし。(計算制度が重要な時は、Runge-Kutta法ではダメなことに注意。参考:ピタゴラス3体問題)
もちろん、Mathematicaでやるのが簡単。
a = 60;
b = 4;
c = 3;
d = 150;
endTime = 0.2;
solution = NDSolve[{
r'[t] == (a - b f[t]) r[t], r[0] == 20,
f'[t] == (c r[t] - d) f[t], f[0] == 10},
{r, f},
{t, 0, endTime}
];
結果の図示で苦労している場合ではなくて、
Plot[Evaluate[{r[t], f[t]} /. solution], {t, 0, endTime}];
ParametricPlot[{r[t], f[t]} /. solution, {t, 0, endTime}];
いつもMathematicaとMaximaじゃおもしろくないから、Octaveを使ってみよう(Maximaで微分方程式を数値的に解く方法知らないし)。
お手軽に使いたいならKnoppix/Mathが一番かな(ISOイメージのダウンロードして、VMware Server上で起動すれば使える)。手元にUbuntuがあるなら「sudo apt-get install octave」のほうが簡単。
WindowsでCygwinを使っているなら、Mathの中にあるgnuplotとOctaveをインストールすればいい(startx &でXを起動して、octave。Cygwin Shellから使いたいなら、xterm上でxhost +、Cygwin Shell上でexport DISPLAY=localhost:0.0; octaveかな)。
#関数定義
function xdot=f(x,t)
a=60;
b=4;
c=3;
d=150;
xdot(1)=(a-b*x(2))*x(1);
xdot(2)=(c*x(1)-d)*x(2);
endfunction
#初期条件
x0=[20;10];
#積分範囲
t=linspace(0,0.2,100);
#積分
x=lsode("f",x0,t);
plot(x)
#xは行列になっていて、x(:,1)で第1列を取り出せるから
plot(x(:,1),x(:,2))