落体の運動とEuler法

\documentclass{jsarticle}
\begin{document}

\subsection*{落体の運動方程式をEuler法で解く}

鉛直下方向を鉛直方向の位置の変数$y$の正の方向とする。
時刻$t_0$で位置$y_0$、速度$v_0$であるような自由落下する質点の運動方程式は、
\[
m\frac{\mathrm{d}^2y}{\mathrm{d}t^2} = mg
\]
である。ここで、$m$は質点の質量、$g$は重力加速度である。

これは位置$y$と速度$v$の連立一階常微分方程式
\begin{eqnarray*}
\frac{\mathrm{d}y}{\mathrm{d}t} &=& v \\
\frac{\mathrm{d}v}{\mathrm{d}t} &=& g
\end{eqnarray*}
と書き直せるので、Euler法により、
\begin{eqnarray*}
y_{n+1} &=& y_n + hv_n \\
v_{n+1} &=& v_n + hg
\end{eqnarray*}
のように差分方程式で近似できる。

これを解くと、$n \ge 1$で、
\begin{eqnarray*}
y_n &=& y_0 + h\left(nv_0 + \frac{n(n-1)}{2}hg\right) = y_0 + nhv_0 + \frac{1}{2}gnh(n-1)h \\
v_n &=& v_0 + nhg
\end{eqnarray*}
となる。定義から$nh = t_n - t_0$であるので、
\begin{eqnarray*}
y_n &=& y_0 + v_0(t_n - t_0) + \frac{1}{2}g(t_n - t_0)(t_{n-1} - t_0) \\
v_n &=& v_0 + g(t_n - t_0)
\end{eqnarray*}

さて、この落体の運動方程式はもちろん簡単に解けて、
\begin{eqnarray*}
y &=& y_0 + v_0(t - t_0) + \frac{1}{2}g(t - t_0)^2 \\
v &=& v_0 + g(t - t_0)
\end{eqnarray*}
したがって、、Euler法による$y_n$の近似解の誤差は、
\begin{eqnarray*}
& & \frac{1}{2}g(t_n - t_0)(t_{n-1} - t_0) - \frac{1}{2}g(t_n - t_0)^2 \\
&=& \frac{1}{2}g(t_n - t_0)(t_{n-1} - t_n) \\
&=& - \frac{1}{2}gh(t_n - t_0)
\end{eqnarray*}
であり、$h(t_n - t_0)$に比例する。
すなわち、初期時刻$t_0$からの時間経過とともに誤差は線形に増えていき、
また、刻み幅$h$を一桁小さくすれば同時刻における誤差も一桁小さくなる。

$h(t_n - t_0) = nh^2$であるが、誤差が$h^2$に比例するわけではないことに注意する。
$h$が半分になれば同じ時刻であるためには$n$が倍でなければならない。
したがって$h(t_n - t_0)$で表現されるように同時刻の誤差は$h$に比例する。

\end{document}