常微分方程式の数値解法 2

今回の講義は常微分方程式の数値解法のベクトル編です（従属変数がスカラーではない場合です；すなわち，前回の記号でいえば，\(d\) が \(2\) 以上の場合です）．

例題として，\(d=2\) の場合の次の微分方程式を考えましょう： \begin{equation} \frac{\mathrm d}{\mathrm dt} \begin{bmatrix} q \cr p \end{bmatrix}= \begin{bmatrix} p \cr -q\end{bmatrix}, \quad \begin{bmatrix} q(0) \cr p(0) \end{bmatrix} = \begin{bmatrix} 1 \cr 0 \end{bmatrix} . \end{equation} 前回の表記との対応は \(y = [q,p]^\top\) です． この方程式は，調和振動子という微分方程式の特殊ケースで，厳密解は \begin{equation} \begin{bmatrix} q(t) \cr p(t) \end{bmatrix} = \begin{bmatrix} \cos (t) \cr -\sin(t) \end{bmatrix} \end{equation} です．すなわち，\((q,p)\)-平面において単位円を時計回りに運動します．

今回は，前回の陽的 Euler 法と Runge-Kutta 法に加えて，symplectic Euler 法とよばれる解法を実装してみましょう． Symplectic Euler 法とは次のような解法です： \begin{align} q_1 &= q_0 + h p_0, \cr p_1 &= p_0 - h q_1. \end{align} \( p\) についての更新式の右辺の \( q_1\) がもし \(q_0 \) であれば陽的 Euler 法と一致します（ここだけが陽的 Euler 法と異なります）． なお，一般には \begin{align} q_{n+1} &= q_{n} + h p_{n}, \cr p_{n+1} &= p_{n} - h q_{n+1} \end{align} です．

では，プログラムを見て課題を行っていきましょう．

サンプルファイルはこちら