Webページでは数式を表すためにLaTeX表記が使えるMathJaxを利用しています。 WebブラウザにはSafari/Chrome/Firefoxを使って下さい（IEでは表示できないようです。）

Euler法とRunge-Kutta法の直接比較

Lorenz方程式では、原点付近から出発した解軌道は無限遠には離れないことがわかっているために、刻み時間幅 $\Delta t$ をある程度小さくして原点付近から出発させた数値解がいきなり無限編に発散してしまうようなことはない。

では、Runge-Kutta法はEuler法に比べてどれほど改良されているのかは直ちに明らかではない。 比較検討のためには、厳密解が分かっている場合について　調べてみるとよい。

演習（調和振子）：

質量 $m$ の質点がバネ定数 $K=k^2 (k>0)$ のバネでつながれて直線上を運動する調和振子の微分方程式は、位置を $x$ および運動量を $p$ として次のようになる。

\begin{align*} \frac{dx}{dt} &= \frac{p}{m}\\ \frac{dp}{dt} &= -k^2 x \end{align*}

次の関数 $E(t)$ \begin{equation} E(t)=\frac{1}{2m}p(t)^2+\frac{k^2}{2}x(t)^2 \end{equation} を微分すると直ぐに確かめられるように$dE/dt=0$ つまり $E$ は一定で（エネルギー保存則）で、次式。 \begin{align*} x(t) &= \frac{\sqrt{2E}}{k} \sin\left( \frac{k}{\sqrt{m}}t + \theta_0 \right)\\ p(t) &=\sqrt{2mE}\cos\left( \frac{k}{\sqrt{m}}t + \theta_0 \right) \end{align*} は解になっていて、$(x,p)$ は楕円上を周期 $\displaystyle\frac{2\sqrt{m}}{k}\pi$ で周回する。

簡単のために、$m=1$, $k=\sqrt{2}$, $E=1$, $\theta_0=0$、つまり $x(0)=0$, $p(0)=\sqrt{2}$ を初期条件とする運動を考える（周期 $\sqrt{2}\pi$ ）。

  1. Euler法とRunge-Kutta法とで求めた数値解 $(x_n,p_n)$ を（csvファイルに書き出して）それぞれプロットして、比較してみなさい。
  2. Euler法とRunge-Kutta法とで求めた数値解 $(x_n,p_n)$ の式(1)のエネルギー値 $E_n = \frac{1}{2m}p_n^2 + \frac{k^2}{2}x_n^2$ を計算し、その値を比較しなさい。