Algorytm skokowy

Algorytm skokowy (algorytm żabiego skoku, Szablon:Ang.) – 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}=-\mathrm{\nabla }V(x),$

z funkcją energii zadaną wzorem:

$E(x,v)=\frac{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 $\mathrm{\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 $\mathrm{\Delta }t$ jest stały oraz $\mathrm{\Delta }t⩽2/\omega$.

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

$\begin{array}{cc}{x}_i& =x_{i-1}+v_{i-1/2}\mathrm{\Delta }t,\\ a_i& =F(x_i)\\ v_{i+1/2}& =v_{i-1/2}+a_i\mathrm{\Delta }t,\end{array}$

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 $\mathrm{\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 $\mathrm{\Delta }t$ musi być stały, aby zachować stabilność.

$\begin{array}{cc}{x}_{i+1}& =x_i+v_i\mathrm{\Delta }t+\frac{1}{2}{a}_i\mathrm{\Delta }{t}^2,\\ v_{i+1}& =v_i+\frac{1}{2}(a_i+a_{i+1})\mathrm{\Delta }t.\end{array}$

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\overline{r}}{d{t}^2}=\frac{\overline{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

• Szablon:Cytuj stronę, Drexel University