6 連立微分方程式は怖くない

前節まで、一つの未知関数に関する微分方程式だけを扱ってきたが、 目標とする Kepler 問題では、 未知関数が複数ある連立微分方程式を解く必要がある。

解の存在を問題にしたり、数値解法を用いたりする場合、 未知関数が複数個あることで本質的な難しさは生じない。 Mathematica でも扱いは難しくない。

逆に、微分の階数が高い場合、未知関数を増やすことで、 1階の微分方程式に帰着して考えるというテクニックがしばしば使われるくらいである。 これを説明した上で、Mathematica で連立微分方程式を解いてみよう。

単振り子の方程式

$$\theta''(t)=-\sin\theta(t),\quad \theta(0)=\theta_0,\quad \theta'(0)=b_0$$

を例に取る。

$$x(t):=\theta(t),\quad y(t):=\theta'(t)$$

とおくと、

$$x'(t)=y(t),\quad y'(t)=-\sin x(t),\quad x(0)=\theta_0,\quad y(0)=b_0$$

となる (実際 $x'(t)=\theta'(t)=y(t)$ , $y'(t)=\theta''(t)=-\sin\theta(t)=-\sin x(t)$ , $x(0)=\theta(0)=\theta_0$ , $y(0)=\theta'(0)=b_0$ )。

$\theta_0=0$ , $b_0=1.6$ の場合、この問題を解く Mathematica のコマンドは、 次のようになる。

  sol=NDSolve[{x'[t]==y[t], y'[t]==-Sin[x[t]],
             x[0]==0, y[0]==1.6},
             {x,y}, {t,0,20}]
  g1=Plot[x[t] /. sol, {t,0,10}]
  g2=ParametricPlot[{x[t],y[t]} /. sol, {t,0,10}]

結果はもちろん前節のものと同じである。

複数の初期値に対する、曲線 $(x(t),y(t))$ を描いてみよう (微分方程式の本に良く載っている)。

sol = Table[{x[t], y[t]} /. 
       NDSolve[{x'[t] == y[t], y'[t] == -Sin[x[t]], x[0]==a*Pi,	y[0]==b},
        {x, y}, {t, -20, 20}], {b, -2.4, 2.4, 0.2}, {a, -2, 2, 2}];
g = ParametricPlot[sol, {t, -20, 20}, PlotRange -> {{-3 Pi, 3 Pi}, {-Pi,Pi}}]

図 6: 振り子の運動の相図 (横軸は角度、縦軸は角度の時間変化率)