Algorytm skokowy

Z Wikipedii, wolnej encyklopedii
Skocz do: nawigacja, szukaj

Algorytm skokowy (algorytm żabiego skoku, ang. leapfrog integration) to 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 szczególnie często w mechanice klasycznej przy opisie układów dynamicznych.

Opis algorytmu[edytuj | edytuj kod]

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

\ddot x=-\nabla V(x),

z funkcją energii zadaną wzorem:

E(x,v)=\tfrac12|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 czasu, rozłożonych w ten sposób, że "przeskakują one nad sobą" – stąd nazwa algorytmu (ang. leapfrogskok 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 \leq 2/\omega[1].

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

\begin{align}
  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{align}

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[2]. Jednak nawet w tej zsynchronizowanej postaci, krok czasowy  \Delta t musi być stały, aby zachować stabilność[3].

\begin{align}
  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{align}

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[4].

Przykładowa implementacja[edytuj | edytuj kod]

Poniżej znajduję 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 soba 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 poczatkowe
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ż[edytuj | edytuj kod]

Przypisy

  1. [C. K. Birdsall and A. B. Langdon, Plasma Physics via Computer Simulations, McGraw-Hill Book Company, 1985, p. 56]
  2. 4.1 Two Ways to Write the Leapfrog
  3. [Skeel, R. D., "Variable Step Size Detabilizes the Stömer/Leapfrog/Verlet Method," BIT Numerical Mathematics, Vol. 33, 1993, pp. 172-175.]
  4. Bishop Christopher: Pattern Recognition and Machine Learning. Nowy Jork: Springer-Verlag, 2006, s. 548–554. ISBN 978-0-387-31073-2.

Linki zewnętrzne[edytuj | edytuj kod]