蛙飛び法

ファインマン物理学の第Ｉ巻の第９章に蛙跳び法による運動のシミュレーションが取り上げられている。バネの運動の場合，運動方程式は，$m \ddot{x} = -k x$であるが，これを$\dot{x}=v$と$\dot{v}= -k/m x = -\lambda x$としてオイラー法を適用する。$\epsilon=t_{n+1}-t_n, \ f_n=f(t_n)$などとして，
前進差分は，
\begin{equation}
\begin{aligned}
x_{n+1}=x_{n}+\epsilon v_{n}\\
v_{n+1}=v_{n}- \lambda \epsilon x_{n}
\end{aligned}
\end{equation}
後退差分は，
\begin{equation}
\begin{aligned}
x_{n+1}=x_{n}+\epsilon v_{n+1}\\
v_{n+1}=v_{n}- \lambda \epsilon x_{n+1}
\end{aligned}
\end{equation}
中心差分は，
\begin{equation}
\begin{aligned}
x_{n+1}=x_{n-1}+2\epsilon v_{n}\\
v_{n+1}=v_{n-1}- 2\lambda \epsilon x_{n}
\end{aligned}
\end{equation}
中心差分の片方をずらすと蛙跳び法の表式が得られる。
\begin{equation}
\begin{aligned}
x_{n+2}=x_{n}+2\epsilon v_{n+1}\\
v_{n+1}=v_{n-1}- 2\lambda \epsilon x_{n}
\end{aligned}
\end{equation}
一方，前進差分と後退差分は，次のように表せる。
\begin{equation}
\begin{aligned}
\begin{pmatrix}x_{n+1} \\ v_{n+1}\end{pmatrix} &=
\begin{pmatrix} 1 & \epsilon \\ - \lambda \epsilon & 1 \end{pmatrix}
\begin{pmatrix}x_{n} \\ v_{n}\end{pmatrix} \\
\begin{pmatrix} 1 & -\epsilon \\ \lambda \epsilon & 1 \end{pmatrix}
\begin{pmatrix}x_{n+1} \\ v_{n+1}\end{pmatrix}
&= \begin{pmatrix}x_{n} \\ v_{n}\end{pmatrix}
\end{aligned}
\end{equation}
陰解法となっている後退差分を解いて$ o(\epsilon^2)$まで考えると，
\begin{equation}
\begin{pmatrix}x_{n+1} \\ v_{n+1}\end{pmatrix}
=\begin{pmatrix}\dfrac{1}{1+\lambda \epsilon^2} & \epsilon \\
-\lambda \epsilon &  \dfrac{1}{1+\lambda \epsilon^2}\end{pmatrix}
\begin{pmatrix}x_{n} \\ v_{n}\end{pmatrix}
\end{equation}
前進差分と後退差分の平均をとると，中心差分に相当するものが得られる。
\begin{equation}
\begin{pmatrix}x_{n+1} \\ v_{n+1}\end{pmatrix}
=\begin{pmatrix}1 - \lambda \epsilon^2 /2 & \epsilon \\
-\lambda \epsilon &  1 -\lambda \epsilon^2 /2 \end{pmatrix}
\begin{pmatrix}x_{n} \\ v_{n}\end{pmatrix}
\end{equation}