微分方程式、それが目的ではないとき


学生向け

微分方程式を数値的に簡単に解く方法を教えてください。

と訊かれた。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))