Algorytm skokowy

Algorytm skokowy (algorytm żabiego skoku, ang. leapfrog integration) – metoda numeryczna służąca do całkowania równań w postaci:

${\ddot {x}}=F(x),$

lub w równoważnej formie:

${\dot {v}}=F(x),\;{\dot {x}}\equiv v,$

występującej często w mechanice klasycznej przy opisie układów dynamicznych.

Opis algorytmu

Problemy typu ${\ddot {x}}=F(x)$ często przybierają postać:

${\ddot {x}}=-\nabla V(x),$

z funkcją energii zadaną wzorem:

$E(x,v)={\tfrac {1}{2}}|v|^{2}+V(x),$

gdzie $V$ jest energią potencjalną układu. Całkowanie metodą skokową (leapfrog) jest równoważne aktualizacji pozycji $x(t)$ i prędkości $v(t)={\dot {x}}(t)$ w naprzemiennych chwilach, rozłożonych w ten sposób, że „przeskakują one nad sobą” – stąd nazwa algorytmu (ang. leapfrog – skok przez plecy). Na przykład położenie jest aktualizowane dla całkowitych wielokrotności kroku czasowego, zaś wartość prędkości jest uaktualniana w chwilach czasu przesuniętych o $\Delta t/2.$

Algorytm skokowy jest metodą drugiego rzędu, w przeciwieństwie do metody Eulera, która jest metodą pierwszego rzędu, jednak liczba wywołań funkcji w każdym kroku jest taka sama. W przeciwieństwie do całkowania metodą Eulera algorytm jest stabilny w przypadku ruchu oscylacyjnego z częstością $\omega ,$ dopóki krok czasowy $\Delta t$ jest stały oraz $\Delta t\leqslant 2/\omega$.

W algorytmie skokowym równania opisujące położenie i prędkość w kolejnych chwilach czasu są następujące:

${\begin{aligned}x_{i}&=x_{i-1}+v_{i-1/2}\,\Delta t,\\[0.4em]a_{i}&=F(x_{i})\\[0.4em]v_{i+1/2}&=v_{i-1/2}+a_{i}\,\Delta t,\end{aligned}}$

gdzie $x_{i}$ jest położeniem w kroku $i,$ $v_{i+1/2\,}$ jest prędkością, lub pierwszą pochodną $x,$ w kroku $i+1/2,$ $a_{i}=F(x_{i})$ jest przyspieszeniem, lub druga pochodną $x,$ w kroku $i$ oraz $\Delta t$ rozmiarem każdego kroku. Równania te mogą być zapisane również w postaci, w której prędkość jest również wyznaczana dla całkowitych wielokrotności kroku czasowego. Jednak nawet w tej zsynchronizowanej postaci, krok czasowy $\Delta t$ musi być stały, aby zachować stabilność.

${\begin{aligned}x_{i+1}&=x_{i}+v_{i}\,\Delta t+{\tfrac {1}{2}}\,a_{i}\,\Delta t^{\,2},\\[0.4em]v_{i+1}&=v_{i}+{\tfrac {1}{2}}\,(a_{i}+a_{i+1})\,\Delta t.\end{aligned}}$

Jednym z częstszych zastosowań algorytmu skokowego są symulacje oddziaływań grawitacyjnych, gdzie przyspieszenie zależy tylko od pozycji. Jednakowoż nawet w tym przypadku metody wyższego rzędu (np. metody Rungeggo-Kutty) są używane znacznie częściej.

Algorytm skokowy ma dwie znaczące zalety:

• odwracalność w czasie – oznacza, że po obliczeniu n kroków możliwe jest odwrócenie kierunku całkowania i dojście do tej samej pozycji wyjściowej.
• symplektyczność – oznacza, że algorytm zachowuje (nieznacznie zmodyfikowaną) całkę energii systemów dynamicznych (tzn. całkowita energia układu oscyluje wokół pewnej stałej wartości bliskiej początkowej energii całkowitej). Jest to szczególnie przydatne podczas obliczania orbit w układach dynamicznych, gdzie pozostałe schematy całkowania, takie jak metody Rungego-Kutty, nie zachowują energii układu i pozwalają układowi na znaczne płynięcie w czasie.

Ze względu na powyższe dwie cechy, algorytm skokowy jest również stosowany w przypadku metod Monte Carlo.

Przykładowa implementacja

Poniżej znajduje się kod napisany w środowisku MATLAB realizujący algorytm żabiego skoku w zsynchronizowanej postaci. Program rozwiązuje równanie: ${\frac {d^{t}{\bar {r}}}{dt^{2}}}={\frac {\bar {r}}{r^{3}}}$ opisujące oddziaływanie grawitacyjne dwóch ciał. Przy czym dobrano taki układ odniesienia, w którym suma dwóch oddziałujących ze sobą mas wynosi 1 oraz stała grawitacji jest również równa 1.

close all
clear all
clc
% rozwiązanie równania d^r/dt^2 = - r/r^3 metodą leapfrog
dt = 0.001; % krok czasowy
% warunki początkowe
r = [0; 0.5; 0.75];
v = [0.25; 0; 1];
r2 = r'*r;
a = -r/(r2*sqrt(r2));
% czas symulacji
T = 10;
res = NaN(round(T/dt)+1,6);
idx = 0;
for t=0:dt:T
    idx = idx + 1;
    v = v + 0.5*a*dt;
    r = r + v*dt;
    r2 = r'*r;
    a = -r/(r2*sqrt(r2));
    v = v + 0.5*a*dt;
    res(idx,:) = [r' v'];
end
comet3(gca,res(:,1),res(:,2),res(:,3))

Zobacz też

• metoda Eulera
• algorytm Rungego-Kutty
• algorytm Verleta

Linki zewnętrzne

• Opis algorytmu. einstein.drexel.edu. [zarchiwizowane z tego adresu (2016-03-06)]., Drexel University