Leapfrog-Verfahren

Das Leapfrog-Verfahren ist eine einfache Methode zur numerischen Integration einer gewöhnlichen Differentialgleichung vom Typ

${\displaystyle {\ddot {x}}={\dot {v}}=a(x)}$.

beziehungsweise allgemeiner von konservativen Systemen ${\displaystyle {\ddot {\mathbf {x} }}=-{\frac {\mathrm {grad} \,V(\mathbf {x} )}{M}}}$ die dem 2. Newtonschen Axiom der klassischen Dynamik folgen und beispielsweise die Bewegung eines oder mehrerer Objekte in einem Potentialfeld ${\displaystyle V}$ beschreiben:

${\displaystyle {\ddot {\vec {x}}}_{k}=-{\frac {1}{m_{k}}}\nabla _{{\vec {x}}_{k}}V\left({\vec {x}}_{1},\dots ,{\vec {x}}_{N}\right),\quad k=1,\dots ,N.}$

Die Leapfrog-Integration ist eine Methode zweiter Ordnung und liefert deshalb im Allgemeinen genauere Ergebnisse als zum Beispiel das eulersche Polygonzugverfahren, das nur von erster Ordnung ist. Außerdem ist es invariant unter Zeitumkehr und erhält in physikalischen Problemstellungen Größen wie den Impuls und Drehimpuls, die auch Erhaltungsgrößen des Originalsystems sind, exakt. Des Weiteren wird eine gestörte Energiefunktion in Ordnung 3 erhalten, während das Verfahren die globale Konvergenzordnung 2 hat.

Die Leapfrog-Integration berechnet abwechselnd die Positionen ${\displaystyle x}$ und die Geschwindigkeiten ${\displaystyle v}$ zu unterschiedlichen Zeitpunkten, ähnlich wie beim Bockspringen (engl. leapfrog). Die Schritt-Gleichungen für das Verfahren lauten:

${\displaystyle {\begin{aligned}x_{i+1}&=x_{i}+v_{i+1/2}\,\Delta t\\[.3em]v_{i+3/2}&=v_{i+1/2}+a(x_{i+1})\,\Delta t\end{aligned}}}$

mit den Startwerten ${\displaystyle x_{0}}$ und ${\displaystyle v_{1/2}}$.

Durch lineare Interpolation von Zwischenwerten kann das Leapfrog-Verfahren als Kombination der zwei Varianten des symplektischen Eulerverfahrens betrachtet werden:

${\displaystyle {\begin{aligned}{\text{(SE1)}}\quad &\left\{{\begin{array}{rcl}v_{i+1/2}&=&v_{i}+a(x_{i}){\frac {\Delta t}{2}}\\[.5em]x_{i+1/2}&=&x_{i}+v_{i+1/2}{\frac {\Delta t}{2}}\end{array}}\right.\\[.7em]{\text{(SE2)}}\quad &\left\{{\begin{array}{rcl}x_{i+1}&=&x_{i+1/2}+v_{i+1/2}{\frac {\Delta t}{2}}\\[.5em]v_{i+1}&=&v_{i+1/2}+a(x_{i+1}){\frac {\Delta t}{2}}\end{array}}\right.\\\end{aligned}}}$

Jeder einzelne Schritt und damit auch die Zusammensetzung ist eine symplektische Transformation und erhält daher Volumina im Phasenraum. Daraus folgt auch die exakte Erhaltung von Impuls und Winkelgeschwindigkeit, soweit das exakte System diese erhält. Einsetzen von (SE1) und (SE2) ineinander führt zu:

${\displaystyle {\begin{array}{rcl}x_{i+1}&=&x_{i}+v_{i}\Delta t+a(x_{i}){\frac {\Delta t^{2}}{2}}\\v_{i+1}&=&v_{i}+\left(a(x_{i+1})+a(x_{i})\right){\frac {\Delta t}{2}}\end{array}}}$

D. h., dass der neue Ort durch Taylorentwicklung bis zur zweiten Ordnung ausgehend vom alten Ort im Phasenraum projiziert wird und die neue Geschwindigkeit bis zur ersten Ordnung; allerdings mit einer modifizierten Beschleunigung, die der Mittelwert der beiden Beschleunigungen zu den Zeitpunkten ${\displaystyle t_{i}}$ und ${\displaystyle t_{i+1}}$ ist.

Eliminiert man aus der Leapfrog-Version die Geschwindigkeitsberechnungen, so ergibt sich

${\displaystyle {\begin{aligned}(x_{i+1}-x_{i})-(x_{i}-x_{i-1})&=(v_{i+1/2}-v_{i-1/2})\,\Delta t\\[.3em]x_{i+1}-2x_{i}+x_{i-1}&=a(x_{i})\,\Delta t^{2}\\[.7em]\implies \quad x_{i+1}=2x_{i}-x_{i-1}&+a(x_{i})\,\Delta t^{2},{\text{ mit Startwerten }}x_{0},x_{1}=x_{0}+v_{0}\,\Delta t+{\tfrac {1}{2}}a(x_{0})\,\Delta t^{2}\end{aligned}}}$

das Verlet-Verfahren das als symmetrische Diskretisierung von ${\displaystyle {\ddot {x}}=a(x)}$ auch direkt hergeleitet werden kann. Diese Diskretisierung hat einen lokalen Fehler von ${\displaystyle O(\Delta t^{4})}$ und damit (wegen der zweifachen Integration) einen globalen Fehler der Größe ${\displaystyle O(t_{N}^{2}\Delta t^{2})}$ für die Differenz ${\displaystyle x(t_{N})-x_{N}}$ zwischen exakter und Näherungslösung zur Endzeit ${\displaystyle t_{N}=t_{0}+N\Delta t}$. Eine weitere Variante ist unter dem Namen Velocity-Verlet bekannt^[1].

Eine erste Beschreibung dieses Verfahrens wurde von Richard Feynman in Isaac Newtons Principia von 1687 in einem Argument zur Herleitung der Keplerschen Gesetze aus den Bewegungsgleichungen gefunden. Neben anderen wurden Varianten dieses Verfahrens 1860 von J. F. Encke und 1907 von C. Störmer verwendet.^[2]

Betrachtet man die Schwingungsgleichung ${\displaystyle {\ddot {x}}+x=0}$ mit der exakten Lösung ${\displaystyle x(t)=C\,\cos(t+\phi )}$, so erhält man für ${\displaystyle {\dot {x}}=v,\;{\dot {v}}=-x}$ in der Einschrittformulierung den Übergang

${\displaystyle {\binom {x_{k+1}}{v_{k+1}}}={\begin{pmatrix}1-{\tfrac {1}{2}}\Delta t^{2}&\Delta t\\-\Delta t(1-{\tfrac {1}{4}}\Delta t^{2})&1-{\tfrac {1}{2}}\Delta t^{2}\end{pmatrix}}{\binom {x_{k}}{v_{k}}}.}$

Es stellt sich heraus, dass das modifizierte Energiefunktional ${\displaystyle E_{\Delta t}(v,x)={\tfrac {1}{2}}\left(v^{2}+(1-{\tfrac {1}{4}}\Delta t^{2})x^{2}\right)}$ exakt erhalten wird. Genauer gilt ${\displaystyle E_{\Delta t}(v_{n},x_{n})=E(v_{0},x_{0}).}$ Damit hat der Fehler in der Energie eine globale Schranke der Ordnung ${\displaystyle O(\Delta t^{2})}$. Die Näherungslösung verläuft für alle Zeiten auf den durch das konstante modifizierte Energieniveau definierten Ellipsen.

1. ↑ Michael Griebel, Stefan Knapek, Gerhard Zumbusch, Attila Caglar: Numerische Simulation in der Moleküldynamik. Springer, 2004. 
2. ↑ Ernst Hairer, Christian Lubich, Gerhard Wanner: Geometric numerical integration illustrated by the Störmer/Verlet method. In: Acta Numerica. 12. Jahrgang, 2003, S. 399–450, doi:10.1017/S0962492902000144 (psu.edu).