Integracja Verleta

[vɛʁˈlɛ] Całkowanie Verleta ( wymowa francuska: <a i=3>[ ) to metoda numeryczna używana do całkowania równań ruchu Newtona . Jest często używany do obliczania trajektorii cząstek w symulacjach dynamiki molekularnej i grafice komputerowej . Algorytm został po raz pierwszy użyty w 1791 roku przez Jeana Baptiste Delambre i od tego czasu był wielokrotnie odkrywany na nowo, ostatnio przez Loupa Verleta w latach 60. XX wieku do wykorzystania w dynamice molekularnej. Używany był również przez PH Cowell i ACC Crommelin w 1909 r. do obliczenia orbity komety Halleya , a przez Carla Størmera w 1907 r. do zbadania trajektorii cząstek elektrycznych w polu magnetycznym (stąd nazywana jest również metodą Störmera ). Integrator Verleta zapewnia dobrą stabilność numeryczną , a także inne właściwości, które są ważne w układach fizycznych , takie jak odwracalność w czasie i zachowanie postaci symplektycznej w przestrzeni fazowej , bez znaczących dodatkowych kosztów obliczeniowych w stosunku do prostej metody Eulera .

Podstawowy Störmer-Verlet

Dla równania różniczkowego drugiego rzędu typu ${\ Displaystyle {\ ddot {\ mathbf {x}}} (t) = \ mathbf {A} {\ bigl ( } \ mathbf {x} (t) {\ bigr )}}$ z warunkami początkowymi ${\ Displaystyle \ mathbf {x} (t_ {0}) = \ mathbf {x} _ {0}}$ i ${\ Displaystyle {\ kropka {\ mathbf {x}}} (t_ {0}) = \ mathbf {v} _ {0}}$ przybliżone rozwiązanie numeryczne ${\ Displaystyle \ mathbf { x} _ {n} \ około \ mathbf {x} (t_ {n})}$ w czasach ${\ Displaystyle t_ {n} = t_ {0} + n \, \ Delta t }$ z wielkością kroku można uzyskać następującą metodą: ${\ displaystyle \ Delta t> 0}$

1. zestaw ${\ textstyle \ mathbf {x} _ {1} = \ mathbf {x} _ {0} + \ mathbf {v} _ {0 }\,\Delta t+{\tfrac {1}{2}}\mathbf {A} (\mathbf {x} _{0})\,\Delta t^{2}}$ ,
2. dla n = 1, 2, ... iteracja
   $${\ Displaystyle \ mathbf {x} _ {n + 1} = 2 \ mathbf {x} _ {n} - \ mathbf {x} _ {n-1} + \ mathbf {A} (\ mathbf {x} _ {n})\,\Delta t^{2}.}$$

Równania ruchu

Równanie ruchu Newtona dla konserwatywnych układów fizycznych to

${\ Displaystyle {\ boldsymbol {M}} {\ ddot {\ mathbf {x}}} (t) = F {\bigl (}\mathbf {x} (t){\bigr )}=-\nabla V{\bigl (}\mathbf {x} (t){\bigr )},}$

lub indywidualnie

${\ Displaystyle m_ {k} {\ ddot {\ mathbf {x}}} _ { k}(t)=F_{k}{\bigl (}\mathbf {x} (t){\bigr )}=-\nabla _{\mathbf {x} _{k}}V\left(\mathbf {x} (t)\prawo),}$

Gdzie

• ${\ displaystyle t}$ to czas,
• ${\ Displaystyle \ mathbf {x} (t) = {\ bigl (} \ mathbf {x} _ {1} (t), \ldots,\mathbf {x} _ {N} (t) {\ bigr )}}$ to zespół wektora pozycji obiektów , ${\ displaystyle N}$
• jest funkcją potencjału skalarnego, ${\ displaystyle V}$
• ${\ displaystyle F}$ to ujemny gradient potencjału , dający zespół sił działających na cząstki, fa {\ displaystyle F}
• $jest$ macierzą mas , zwykle ukośną z blokami o $każdej$ cząstki

To równanie, dla różnych wyborów potencjalnej funkcji , może być użyte do opisania ewolucji różnych układów fizycznych, od ruchu $cząsteczek$ po orbity planet .

Po przekształceniu mającym na celu przeniesienie masy na prawą stronę i zapomnieniu o strukturze wielu cząstek, równanie można uprościć do

${\ Displaystyle {\ ddot {\ mathbf {x}}} (t) = \ mathbf {A} {\ bigl (} \ mathbf {x} (t) { \bigr )}}$

z odpowiednią funkcją o wartościach wektorowych reprezentującą przyspieszenie zależne od położenia. ${\ Displaystyle \ mathbf {A} (\ mathbf {x})} .$$początkowa$ pozycja początkowa $x$ prędkość również podane.

Całkowanie Verleta (bez prędkości)

Aby zdyskretyzować i numerycznie rozwiązać ten problem z wartością początkową , wybiera się krok czasowy i sekwencję punktów próbkowania $t$$= n \ ,\Delta t}$ brane pod uwagę. Zadanie polega na skonstruowaniu sekwencji punktów $(t_ {n})}$$x$ \ na trajektorii rozwiązania dokładnego.

Tam, gdzie metoda Eulera wykorzystuje przybliżenie różnicy w przód do pierwszej pochodnej w równaniach różniczkowych rzędu pierwszego, całkowanie Verleta można postrzegać jako wykorzystujące przybliżenie różnicy centralnej do drugiej pochodnej:

${\ Displaystyle {\ rozpocząć {wyrównane}} {\ Frac {\ Delta ^ {2} \ mathbf {x} _ {n}} {\ Delta t ^ {2}}} & = {\ Frac {{\ Frac {\ mathbf {x} _{n+1}-\mathbf {x} _{n}}{\Delta t}}-{\frac {\mathbf {x} _{n}-\mathbf {x} _{n -1}}{\Delta t}}}{\Delta t}}\\[6pt]&={\frac {\mathbf {x} _{n+1}-2\mathbf {x} _{n} +\mathbf {x} _{n-1}}{\Delta t^{2}}}=\mathbf {a} _{n}=\mathbf {A} (\mathbf {x} _{n}) .\end{wyrównane}}}$

Całkowanie Verleta w postaci używanej jako metoda Störmera wykorzystuje to równanie do uzyskania następnego wektora położenia z dwóch poprzednich bez użycia prędkości jako

${\ Displaystyle {\ rozpocząć {wyrównane} \ mathbf {x} _ {n + 1} & = 2 \ mathbf {x} _ {n} - \ mathbf {x} _ {n-1} + \ mathbf {a} _{n}\,\Delta t^{2},\\[6pt]\mathbf {a} _{n}&=\mathbf {A} (\mathbf {x} _{n}).\end{ wyrównany}}}$

Błąd dyskretyzacji

$stopnia$ poziom lokalnych błędów wprowadzanych do całkowania przez dyskretyzację poprzez usunięcie wszystkich wyrazów nieparzystych stopni, tutaj wyrazów trzeciego . Błąd lokalny jest kwantyfikowany przez wstawienie dokładnych wartości ${\ Displaystyle \ mathbf {x} (t_ {n-1}), \ mathbf {x} (t_ {n}), \ mathbf {x} (t_ {n + 1})} do iteracji i$ obliczania Taylora $t = t_ {n}}$ Displaystyle wektora pozycji ${\ Displaystyle \ mathbf {x} (t ​​\ pm \ Delta t)}$ w różnych kierunkach czasowych: t = t n {\ displaystyle t = t_ {n } }

${\ Displaystyle {\ rozpocząć {wyrównane} \ mathbf {x} (t ​​+ \ Delta t) & = \ mathbf {x} (t) + \ mathbf {v} (t) \ Delta t + {\ Frac {\ mathbf {a } (t)\Delta t^{2}}{2}}+{\frac {\mathbf {b} (t)\Delta t^{3}}{6}}+{\mathcal {O}}\ left(\Delta t^{4}\right)\\\mathbf {x} (t-\Delta t)&=\mathbf {x} (t)-\mathbf {v} (t)\Delta t+{\ frac {\mathbf {a} (t)\Delta t^{2}}{2}}-{\frac {\mathbf {b} (t)\Delta t^{3}}{6}}+{\ mathcal {O}}\left(\Delta t^{4}\right),\end{wyrównane}}}$

gdzie ${\ Displaystyle \ mathbf {x}}$ to pozycja, ${\ Displaystyle \ mathbf {v} = {\ kropka {\ mathbf {x}}}}$ prędkość, ${\ Displaystyle \mathbf {a} = {\ ddot {\ mathbf {x}}}}$ przyspieszenie i $}}}$${\ Displaystyle \ mathbf {b} = {\ kropka {\ mathbf {a}}} = {\ overset {\ kropki }{\ mathbf {x$ } szarpnięcie (trzecia pochodna pozycji względem czasu).

Dodanie tych dwóch rozszerzeń daje

${\ Displaystyle \ mathbf {x} (t ​​+ \ Delta t) = 2 \ mathbf {x} (t) - \ mathbf {x} (t- \ Delta t) + \ mathbf {a} (t) \ Delta t ^ {2}+{\mathcal {O}}\left(\Delta t^{4}\right).}$

Widzimy, że wyrazy pierwszego i trzeciego rzędu z rozwinięcia Taylora znoszą się, czyniąc w ten sposób integrator Verleta rzędem dokładniejszym niż całkowanie przez samo proste rozwinięcie Taylora.

Należy zwrócić uwagę na fakt, że przyspieszenie tutaj jest obliczane z dokładnego rozwiązania, ${\ Displaystyle \ mathbf {a} (t) = \ mathbf {A} {\ bigl (} \ mathbf {x} (t) {\ bigr )}}$ , podczas gdy w iteracji jest obliczany w centralnym punkcie iteracji, ${\ Displaystyle \ mathbf {a} _ {n} =\mathbf {A} (\mathbf {x} _ {n})}$ . Przy obliczaniu błędu globalnego, czyli odległości między rozwiązaniem dokładnym a ciągiem aproksymacyjnym, te dwa wyrazy nie znoszą się dokładnie, wpływając na kolejność błędu globalnego.

Prosty przykład

Aby uzyskać wgląd w relacje między błędami lokalnymi i globalnymi, pomocne jest zbadanie prostych przykładów, w których dokładne rozwiązanie, jak również rozwiązanie przybliżone, można wyrazić za pomocą jawnych wzorów. Standardowym przykładem tego zadania jest funkcja wykładnicza .

Rozważ liniowe równanie różniczkowe ${\ Displaystyle {\ ddot {x}} (t) = w ^ {2} x (t)}$ ze stałą ${\ Displaystyle w}$ . Jego dokładne rozwiązania bazowe to ${\ displaystyle e ^ {wt}}$ i ${\ displaystyle e ^ {-wt}}$ .

Metoda Störmera zastosowana do tego równania różniczkowego prowadzi do liniowej relacji powtarzalności

${\ Displaystyle x_ {n + 1} -2x_ {n} + x_ {n-1} = h ^ {2} w ^ {2}x_{n},}$

Lub

${\ Displaystyle x_ {n + 1} -2 \ lewo (1 + {\ tfrac {1} 2}}(biała)^{2}\prawa)x_{n}+x_{n-1}=0.}$

Można to rozwiązać, znajdując pierwiastki jego charakterystycznego wielomianu ${\ Displaystyle q ^ {2} -2 \ lewo (1 + {\ tfrac { 1}{2}}(wh)^{2}\right)q+1=0}$ . To są

${\ Displaystyle q _ {\ pm} = 1 + {\ tfrac {1} {2}} (wh) ^ {2} \ pm wh {\ sqrt {1 + {\ tfrac {1} {4}} (wh) ^{2}}}.}$

Podstawowymi rozwiązaniami powtarzalności liniowej są ${\ Displaystyle x_ {n} = q_ {+} ^ {n}}$ i ${\ Displaystyle x_ {n} = q_ {-} ^{n}}$ . Aby porównać je z dokładnymi rozwiązaniami, obliczane są rozwinięcia Taylora:

${\ Displaystyle {\ rozpocząć {wyrównane} q_ {+} & = 1 + {\ tfrac {1} {2}} (wh) ^ {2} + wh \ lewo (1 + {\ tfrac {1} {8} }(wh)^{2}-{\tfrac {3}{128}}(wh)^{4}+{\mathcal {O}}\left(h^{6}\right)\right)\\ &=1+(biały)+{\tfrac {1}{2}}(biały)^{2}+{\tfrac {1}{8}}(biały)^{3}-{\tfrac {3} {128}}(wh)^{5}+{\mathcal {O}}\left(h^{7}\right).\end{wyrównane}}}$

Iloraz tego szeregu z ilorazem wykładniczym $Displaystyle 1- {\ tfrac {1 }{24}}(wh)^{3}+{\mathcal {O}}\left(h^{5}\right)}$$się$ 1 , więc

${\ Displaystyle {\ rozpocząć {wyrównane} q_ {+} & = \ lewo (1- {\ tfrac {1} {24}} (wh) ^ {3} + {\ mathcal {O}} \ lewo (h ^ {5}\right)\right)e^{wh}\\&=e^{-{\frac {1}{24}}(wh)^{3}+{\mathcal {O}}\left( h^{5}\right)}\,e^{wh}.\end{wyrównane}}}$

Stąd wynika, że ​​dla pierwszego rozwiązania bazowego błąd można obliczyć jako

${\ Displaystyle {\ rozpocząć {wyrównane} x_ {n} = q_ {+} ^ {n} & = e ^ {- {\ Frac {1} {24}} (wh) ^ {2} \, wt_ {n }+{\mathcal {O}}\left(h^{4}\right)}\,e^{wt_{n}}\\&=e^{wt_{n}}\left(1-{\ tfrac {1}{24}}(wh)^{2}\,wt_{n}+{\mathcal {O}}(h^{4})\right)\\&=e^{wt_{n} }+{\mathcal {O}}\left(h^{2}t_{n}e^{wt_{n}}\right).\end{wyrównane}}}$

Oznacza to, że chociaż lokalny błąd dyskretyzacji jest rzędu 4, ze względu na drugi rząd równania różniczkowego błąd globalny jest rzędu 2, ze stałą rosnącą wykładniczo w czasie.

Rozpoczęcie iteracji

$Displaystyle t = t_ {1} = \ Delta t}$ na początku iteracji Verleta w kroku $t$ czas , obliczanie ${\ displaystyle \ mathbf {x} _ {2}}$ , potrzebny jest już wektor pozycji ${\ displaystyle \ mathbf {x} _ {1}}$ w czasie ${\ displaystyle t = t_ {1}}$ . Na pierwszy rzut oka może to sprawiać problemy, ponieważ warunki początkowe są znane tylko w czasie początkowym ${\ displaystyle t_ {0}=0}$ . Jednak $z$ przyspieszenie odpowiednie położenie w pierwszym kroku czasowym można uzyskać za pomocą wielomianu Taylora stopnia drugiego:

${\ Displaystyle \ mathbf {x} _ {1} = \ mathbf {x} _ {0} + \ mathbf {v} _ {0} \ Delta t + {\ tfrac {1} {2}} \ mathbf {a} _{0}\Delta t^{2}\około \mathbf {x} (\Delta t)+{\mathcal {O}}\left(\Delta t^{3}\right).}$

Błąd w pierwszym kroku czasowym jest więc rzędu ${\ Displaystyle {\ mathcal {O}} \ lewo (\ Delta t ^ {3} \ prawej)}$ . Nie jest to uważane za problem, ponieważ w symulacji obejmującej dużą liczbę kroków czasowych błąd w pierwszym kroku czasowym stanowi tylko pomijalnie małą część całkowitego błędu, który w czasie wynosi t n {\ displaystyle t_ { $}$ rzędu ${\ Displaystyle {\ mathcal {O}} \ lewo (e ^ {Lt_ {n}} \ Delta t ^ {2} \ prawej)} , zarówno dla odległości wektorów$ pozycji ${\ Displaystyle \ mathbf {x } _ {n}}$ do ${\ Displaystyle \ mathbf {x} (t_ {n})}$ co do odległości podzielonych różnic ${\ Displaystyle {\ tfrac {\mathbf {x} _{n+1}-\mathbf {x} _{n}}{\Delta t}}}$ do ${\ Displaystyle {\ tfrac {\ mathbf {x} (t_ {n + 1}) - \ mathbf {x} (t_ {n})} {\ Delta t}}}$ . Co więcej, aby uzyskać błąd globalny drugiego rzędu, błąd początkowy musi być co najmniej trzeciego rzędu.

Niestałe różnice czasu

Wadą metody Störmera-Verleta jest to, że jeśli zmienia się krok czasowy ( $metoda$ nie przybliża rozwiązania równania różniczkowego. Można to poprawić za pomocą wzoru

${\ Displaystyle \ mathbf {x} _ {i + 1} = \ mathbf {x} _ {i} + \ lewo (\ mathbf {x} _ {i} - \ mathbf {x} _ {i-1} \ dobrze){\frac {\Delta t_{i}}{\Delta t_{i-1}}}+\mathbf {a} _{i}\Delta t_{i}^{2}.}$

Dokładniejsze wyprowadzenie wykorzystuje szereg Taylora (do drugiego rzędu) $= t_ {i } + \ Delta t_ {i}}$$+$ Δ i ${\ Displaystyle t_ {i-1} = t_ {i} - \ Delta t_ {i-1}}$ uzyskać po eliminacja ${\ Displaystyle \ mathbf {v} _ {i}}$

${\ Displaystyle {\ Frac {\ mathbf {x} _ {i+1}-\mathbf {x} _{i}}{\Delta t_{i}}}+{\frac {\mathbf {x} _{i-1}-\mathbf {x} _{i }}{\Delta t_{i-1}}}=\mathbf {a} _{i}\,{\frac {\Delta t_{i}+\Delta t_{i-1}}{2}}, }$

tak, że formuła iteracji staje się

${\ Displaystyle \ mathbf {x} _ {i + 1} = \ mathbf {x} _ {i} + (\ mathbf {x} _ {i} - \ mathbf {x} _ {i-1}) {\ frac {\Delta t_{i}}{\Delta t_{i-1}}}+\mathbf {a} _{i}\,{\frac {\Delta t_{i}+\Delta t_{i-1 }}{2}}\,\Delta t_{i}.}$

Obliczanie prędkości – metoda Störmera-Verleta

Prędkości nie są podane wprost w podstawowym równaniu Störmera, ale często są niezbędne do obliczenia pewnych wielkości fizycznych, takich jak energia kinetyczna. Może to stwarzać techniczne wyzwania w dynamiki molekularnej , ponieważ energia kinetyczna i $chwilowe$ temperatury w czasie nie mogą ${\ displaystyle t + \ Delta t}$ obliczone dla układu, dopóki pozycje nie będą znane w czasie . Z tym niedoborem można sobie poradzić za pomocą prędkości Verlet algorytm lub przez oszacowanie prędkości przy użyciu warunków pozycji i twierdzenia o wartości średniej :

${\ Displaystyle \ mathbf {v} (t) = {\ Frac {\ mathbf {x} (t ​​+ \ Delta t) - \ mathbf {x} (t- \ Delta t)} {2 \ Delta t}} + { \mathcal {O}}\left(\Delta t^{2}\right).}$

$,$${\ Displaystyle \ mathbf {v} _ {n} = {\ tfrac {\ mathbf {x} _ {n + 1} - \ mathbf {x} _ {n-1}} {2\Delta t}}}$$n$ położenia, ponieważ dotyczy prędkości w czasie , nie , co oznacza że jest przybliżeniem drugiego rzędu do ${\ Displaystyle \ mathbf {v} (t_ {n})}$ . Z tym samym argumentem, ale o połowę w czasie, ${\ Displaystyle \ mathbf {v} _ {n + {\ Frac {1} {2}}} = { \ tfrac {\ mathbf {x} _ {n + 1} - \ mathbf {x} _ {n}} {\ Delta t}}} jest przybliżeniem drugiego rzędu do v$ ( $Displaystyle \mathbf {v} \left(t_{n+{\frac {1}{2}}}\right)}$ , gdzie ${\ Displaystyle t_ {n + {\ Frac {1} {2}}} = t_ {n} + {\ tfrac {1} {2}} \ Delta t}$ .

Można skrócić przedział, aby przybliżyć prędkość w czasie kosztem dokładności: ${\ displaystyle t + \ Delta t}$

${\ Displaystyle \ mathbf {v} (t + \ Delta t) = {\ Frac {\ mathbf {x} (t ​​+ \ Delta t) - \ mathbf {x} (t)} {\ Delta t}} + {\ mathcal {O}}(\Delta t).}$

Prędkość Verleta

Powiązanym i częściej używanym algorytmem jest algorytm Verleta dotyczący prędkości , podobny do metody przeskakującej , z tym wyjątkiem, że prędkość i położenie są obliczane przy tej samej wartości zmiennej czasowej (żaba przeskakująca nie, jak sugeruje nazwa). Wykorzystuje to podobne podejście, ale wyraźnie uwzględnia prędkość, rozwiązując problem pierwszego kroku czasowego w podstawowym algorytmie Verleta:

${\ Displaystyle {\ rozpocząć {wyrównane} \ mathbf {x} (t ​​+ \ Delta t) & = \ mathbf {x} (t) + \ mathbf {v} (t) \ \ Delta t + {\ tfrac {1} {2}}\,\mathbf {a} (t)\Delta t^{2},\\[6pt]\mathbf {v} (t+\Delta t)&=\mathbf {v} (t)+{ \frac {\mathbf {a} (t)+\mathbf {a} (t+\Delta t)}{2}}\Delta t.\end{wyrównane}}}$

Można wykazać, że błąd w prędkości Verleta jest tego samego rzędu, co w podstawowym Verlecie. Zauważ, że algorytm prędkości niekoniecznie wymaga więcej pamięci, ponieważ w podstawowym Verlecie śledzimy dwa wektory położenia, podczas gdy w Verlecie prędkości śledzimy jeden wektor położenia i jeden wektor prędkości. Standardowy schemat implementacji tego algorytmu to:

1. Oblicz ${\ Displaystyle \ mathbf {v} \ lewo (t + {\ tfrac {1} {2}} \, \ Delta t\right)=\mathbf {v} (t)+{\tfrac {1}{2}}\,\mathbf {a} (t)\,\Delta t}$ .
2. Oblicz ${\ Displaystyle \ mathbf {x} (t ​​+ \ Delta t) = \ mathbf {x} (t) + \ mathbf {v} \left(t+{\tfrac {1}{2}}\,\Delta t\right)\,\Delta t}$ .
3. za $\ Delta t)}$ z potencjału interakcji za pomocą ${\ Displaystyle \ mathbf {x} (t ​​+ \ Delta t)}$ .
4. Oblicz ${\ Displaystyle \ mathbf {v} (t + \ Delta t) = \ mathbf {v} \left(t+{\tfrac {1}{2}}\,\Delta t\right)+{\tfrac {1}{2}}\,\mathbf {a} (t+\Delta t)\Delta t}$ .

Algorytm ten działa również ze zmiennymi krokami czasowymi i jest identyczny z formą integracji metody „skok-drift-kick” .

Eliminując prędkość półkrokową, algorytm ten można skrócić do

1. Oblicz ${\ Displaystyle \ mathbf {x} (t ​​+ \ Delta t) = \ mathbf {x} (t)+\mathbf {v} (t)\,\Delta t+{\tfrac {1}{2}}\,\mathbf {a} (t)\,\Delta t^{2}}$ .
2. za $\ Delta t)}$ z potencjału interakcji za pomocą ${\ Displaystyle \ mathbf {x} (t ​​+ \ Delta t)}$ .
3. Oblicz ${\ Displaystyle \ mathbf {v} (t + \ Delta t) = \ mathbf { v} (t)+{\tfrac {1}{2}}\,{\bigl (}\mathbf {a} (t)+\mathbf {a} (t+\Delta t){\bigr )}\Delta t}$ .

Należy jednak zauważyć, że algorytm ten zakłada, że ​​przyspieszenie $Displaystyle \ mathbf {a}$$\ Delta t)$ od położenia. za i nie zależy od prędkości ${\ Displaystyle \ mathbf {v} (t + \ Delta t)}$ .

Można zauważyć, że długoterminowe wyniki prędkości Verleta i podobnie przeskakującej żaby są o jeden rząd lepsze niż półukryta metoda Eulera . Algorytmy są prawie identyczne aż do przesunięcia o pół kroku czasu w prędkości. Można to udowodnić, obracając powyższą pętlę, aby rozpocząć od kroku 3, a następnie zauważając, że składnik przyspieszenia w kroku 1 można wyeliminować, łącząc kroki 2 i 4. Jedyna różnica polega na tym, że prędkość punktu środkowego w prędkości Verlet jest uważana za prędkość końcową w półukrytej metodzie Eulera.

Błąd globalny wszystkich metod Eulera jest rzędu pierwszego, podczas gdy błąd globalny tej metody jest, podobnie jak metody punktu środkowego , rzędu drugiego. Dodatkowo, jeśli przyspieszenie rzeczywiście wynika z sił w konserwatywnym układzie mechanicznym lub hamiltonowskim , energia przybliżenia zasadniczo oscyluje wokół stałej energii dokładnie rozwiązanego układu, z błędem globalnym ponownie związanym rzędu pierwszego dla półjawnego Eulera i zamów dwa dla Verlet-leapfrog. To samo dotyczy wszystkich innych zachowanych wielkości układu, takich jak pęd liniowy lub kątowy, które są zawsze zachowane lub prawie zachowane w integrator symplektyczny .

${\ tfrac {1} {2}}}$ szczególnym przypadkiem metody Newmark-beta $z$ i γ .

Reprezentacja algorytmiczna

Ponieważ prędkość Verlet jest ogólnie użytecznym algorytmem w aplikacjach 3D, ogólne rozwiązanie napisane w C++ mogłoby wyglądać jak poniżej. Uproszczona siła oporu służy do zademonstrowania zmiany przyspieszenia, jednak jest potrzebna tylko wtedy, gdy przyspieszenie nie jest stałe.

 

          
           
           
        
        struct  Body  {  Vec3d  poz  {  0.0  ,  0.0  ,  0.0  };  Wersja  Vec3d  {  2.0  ,  0.0  ,  0.0  };  // 2 m/s wzdłuż osi x  Vec3d  acc  {  0.0  ,  0.0  ,  0.0  };  // brak przyspieszenia przy pierwszej  podwójnej  masie  =  1,0  ;  //  podwójne  przeciągnięcie 1 kg   =  0,1  ;  

    



      
    
               
            // rho*C*Area – uproszczone przeciąganie dla tego przykładu  /**  * Zaktualizuj pos i vel za pomocą integracji „Velocity Verlet”  * @param dt DeltaTime / time step [np.: 0.01]  */  void  update  (  double  dt  )  {  Vec3d  new_pos  =  poz  +  vel  *  dt  +  acc  *  (  dt  *  dt  *  0.5  );  Vec3d  new_acc  =  zastosuj_siły  ();  
             
          
          
          
    

      
    
             // potrzebne tylko wtedy, gdy przyspieszenie nie jest stałe  Vec3d  new_vel  =  vel  +  (  ac  +  new_acc  )  *  (  dt  *  0.5  );  pozycja  =  nowa_pozycja  ;  poziom  =  nowy_poziom  ;  acc  =  nowy_akc  ;  }  Vec3d  apply_forces  ()  const  {  Vec3d  grav_acc  =  Vec3d  {  0.0  ,  0.0    
                  
              
           
    
 ,  -9,81  };  // 9,81 m/s² w dół w osi Z  Vec3d  drag_force  =  0,5  *  drag  *  (  vel  *  vel  );  // D = 0,5 * (rho * C * Powierzchnia * vel^2)  Vec3d  drag_acc  =  siła_przeciągania  /  masa  ;  // a = F/m  return  grav_acc  -  drag_acc  ;  }  }; 

Terminy błędów

Główny artykuł: błąd obcięcia (całka numeryczna)

Globalny błąd obcięcia metody Verleta wynosi ${\ Displaystyle {\ mathcal {O}} \ lewo (\ Delta t ^ {2} \ prawej)}$ , zarówno dla położenia, jak i prędkości. Kontrastuje to z faktem, że lokalny błąd $tylko$ powyżej Różnica wynika z akumulacji lokalnego błędu obcięcia we wszystkich iteracjach.

Globalny błąd można wyprowadzić, zauważając, co następuje:

${\ Displaystyle \ operatorname {błąd} {\ bigl (} x (t_ {0} + \ Delta t) {\ bigr)} = {\ mathcal {O}}\left(\Delta t^{4}\right)}$

I

${\ Displaystyle x (t_ {0} + 2 \ delta t) = 2x (t_ {0} + \ delta t) -x (t_ {0}) + \ delta t ^ {2} {\ ddot {x}} (t_{0}+\Delta t)+{\mathcal {O}}\left(\Delta t^{4}\right).}$

Dlatego

${\ Displaystyle \ nazwa operatora {błąd} {\ bigl (} x (t_ {0} + 2 \ Delta t) {\ bigr)} = 2 \ cdot \ nazwa operatora {błąd} {\ bigl (} x (t_ {0} +\Delta t){\bigr )}+{\mathcal {O}}\left(\Delta t^{4}\right)=3\,{\mathcal {O}}\left(\Delta t^{ 4}\prawo).}$

Podobnie:

${\ Displaystyle {\ rozpocząć {wyrównane} \ operatorname {błąd} {\ bigl (} x (t_ {0} + 3 \ Delta t) {\ bigl)} i = 6 \,{\mathcal {O}}\left(\Delta t^{4}\right),\\[6px]\operatorname {błąd} {\bigl (}x(t_{0}+4\Delta t) {\bigl )}&=10\,{\mathcal {O}}\left(\Delta t^{4}\right),\\[6px]\operatorname {błąd} {\bigl (}x(t_{ 0}+5\Delta t){\bigl )}&=15\,{\mathcal {O}}\left(\Delta t^{4}\right),\end{wyrównane}}}$

co można uogólnić na (można to wykazać przez indukcję, ale jest tu podane bez dowodu):

${\ Displaystyle \ nazwa operatora {błąd} {\ bigl (} x (t_ {0} + n \ Delta t) {\ bigr)} = {\ Frac {n (n + 1)}} {2}} \, {\ mathcal {O}}\left(\Delta t^{4}\right).}$

$n\Delta t}$$Δ$$\$ błąd $między$ i gdzie , jasne jest, że

${\ Displaystyle \ nazwa operatora {błąd} {\ bigl (} x (t_ {0} +T){\bigr )}=\left({\frac {T^{2}}{2\Delta t^{2}}}+{\frac {T}{2\Delta t}}\right) {\mathcal {O}}\left(\Delta t^{4}\right),}$

a zatem globalny (skumulowany) błąd w stałym przedziale czasu jest określony przez

${\ Displaystyle \ nazwa operatora {błąd} {\ bigr (} x (t_ {0} + T) {\ bigl)} = {\ mathcal {O}} \ lewo (\ Delta t ^ {2} \ prawo).}$

Ponieważ prędkość jest określana w sposób niekumulacyjny na podstawie pozycji w integratorze Verleta, globalny błąd prędkości również wynosi ${\ Displaystyle {\ mathcal {O}} \ lewo (\ Delta t ^ { 2}\prawo)}$ .

W symulacjach dynamiki molekularnej błąd globalny jest zwykle znacznie ważniejszy niż błąd lokalny, dlatego integrator Verleta jest nazywany integratorem drugiego rzędu.

Ograniczenia

Główny artykuł: Algorytm ograniczeń

Układy wielu cząstek z ograniczeniami są łatwiejsze do rozwiązania za pomocą integracji Verleta niż metodami Eulera. Ograniczeniami między punktami mogą być na przykład potencjały ograniczające je do określonej odległości lub siły przyciągania. Można je modelować jako sprężyny łączące cząstki. Wykorzystując sprężyny o nieskończonej sztywności, model można następnie rozwiązać za pomocą algorytmu Verleta.

$)}}$ jednym wymiarze związek między nieograniczonymi pozycjami a $rzeczywistymi$ pozycjami ${\ Displaystyle {\ tylda {x}} _ {i} ^ {$$displaystyle$$t$ punktów w czasie , biorąc pod uwagę pożądaną odległość ograniczenia , można znaleźć za pomocą algorytmu $r}$

${\ Displaystyle {\ rozpocząć {wyrównane} d_ {1} & = x_ {2} ^ {(t)} -x_ {1} ^ {(t)}, \\ [6px] d_ {2} & = \ | d_{1}\|,\\[6px]d_{3}&={\frac {d_{2}-r}{d_{2}}},\\[6px]x_{1}^{(t+ \Delta t)}&={\tylda {x}}_{1}^{(t+\Delta t)}+{\tfrac {1}{2}}d_{1}d_{3},\\[ 6px]x_{2}^{(t+\Delta t)}&={\tylda {x}}_{2}^{(t+\Delta t)}-{\tfrac {1}{2}}d_{ 1}d_{3}.\end{wyrównane}}}$

Całkowanie Verleta jest przydatne, ponieważ bezpośrednio wiąże siłę z położeniem, zamiast rozwiązywać problem za pomocą prędkości.

Problemy pojawiają się jednak, gdy na każdą cząstkę działa wiele sił ograniczających. Jednym ze sposobów rozwiązania tego problemu jest przejście przez każdy punkt symulacji, tak aby w każdym punkcie rozluźnienie ograniczeń ostatniego punktu było już wykorzystane do przyspieszenia rozprzestrzeniania się informacji. W symulacji można to zaimplementować, stosując małe kroki czasowe do symulacji, stosując ustaloną liczbę kroków rozwiązywania ograniczeń na krok czasowy lub rozwiązując ograniczenia, dopóki nie zostaną spełnione przez określone odchylenie.

Podczas lokalnego przybliżania ograniczeń do pierwszego rzędu jest to to samo, co metoda Gaussa – Seidela . Wiadomo, że dla małych macierzy rozkład LU jest szybszy. Duże systemy można podzielić na klastry (na przykład każdy ragdoll = klaster). Wewnątrz klastrów stosowana jest metoda LU, pomiędzy klastrami metoda Gaussa–Seidela . Kod macierzowy może być ponownie użyty: zależność sił od pozycji może być lokalnie przybliżona do pierwszego rzędu, a integracja Verleta może być bardziej ukryta.

Wyrafinowane oprogramowanie, takie jak SuperLU, istnieje do rozwiązywania złożonych problemów przy użyciu rzadkich macierzy. Specyficzne techniki, takie jak stosowanie (klastrów) matryc, mogą być stosowane w celu rozwiązania konkretnego problemu, takiego jak siła rozchodząca się przez arkusz materiału bez tworzenia fali dźwiękowej .

Innym sposobem rozwiązania ograniczeń holonomicznych jest użycie algorytmów ograniczeń .

Reakcje zderzenia

Jednym ze sposobów reagowania na kolizje jest stosowanie systemu opartego na karach, który zasadniczo przykłada ustaloną siłę do punktu po zetknięciu. Problem polega na tym, że bardzo trudno jest wybrać przekazywaną siłę. Użyj zbyt dużej siły, a obiekty staną się niestabilne, zbyt słabe, a przedmioty będą się przenikać. Innym sposobem jest użycie reakcji kolizji projekcji, które biorą punkt obrażający i próbują przesunąć go na jak najkrótszą odległość, aby wysunąć go z drugiego obiektu.

W tym drugim przypadku całkowanie Verleta automatycznie obsłużyłoby prędkość nadaną przez zderzenie; należy jednak pamiętać, że nie gwarantuje się, że zrobi to w sposób zgodny z fizyką zderzeń (to znaczy, że zmiany pędu nie są gwarantowane jako realistyczne). Zamiast niejawnie zmieniać składnik prędkości, należałoby jawnie kontrolować końcowe prędkości zderzających się obiektów (poprzez zmianę zarejestrowanej pozycji z poprzedniego kroku czasowego).

Dwie najprostsze metody określania nowej prędkości to zderzenia doskonale sprężyste i niesprężyste . Nieco bardziej skomplikowana strategia, która zapewnia większą kontrolę, wymagałaby zastosowania współczynnika restytucji .

Zobacz też

• Warunek Couranta-Friedrichsa-Lewy'ego
• Dryf energii
• Integrator symplektyczny
• Integracja Leapfroga
• Algorytm Beemana

Linki zewnętrzne

• Demo i kod integracji Verlet jako aplet Java
• Zaawansowana fizyka postaci autorstwa Thomasa Jakobsena
• Teoria Symulacji Dynamiki Molekularnej – dół strony