2. Całkowanie Verleta

Sercem symulacji jest system cząsteczek. Zazwyczaj w implementacji systemu cząstek każda cząstka ma dwie główne zmienne: x – pozycja i v – jej prędkość. Następnie w 'czasokrokowej’ pętli, nowa pozycja x’ i prędkość v’ są obliczane według:

\textbf{x}' = \textbf{x} + \textbf{v} \cdot \Delta t
\textbf{v}' = \textbf{v} + \textbf{a} \cdot \Delta t

gdzie dt jest krokiem czasu, a – przyspieszeniem obliczanym z wykorzystaniem prawa Newtona \textbf{F} = m \cdot \textbf{a} (gdzie F to wypadkowa sił działających na cząstki). Jest to prosta integracja Eulera.

W naszym przypadku jednak wybieramy schemat w mniejszym stopniu oparty na prędkości: zamiast przechowywać pozycję i prędkość każdej cząstki, będziemy przechowywać jej aktualną pozycję x oraz wcześniejszą \textbf{x}^*. Utrzymując stały czasowy krok, pętla aktualizacyjna wygląda wtedy następująco:

\textbf{x}' = 2\textbf{x} - \textbf{x}^* + \textbf{a} \cdot \Delta t^2
\textbf{x}^* = \textbf{x}

Nazywamy to całkowaniem Verleta (patrz [((bibcite Verlet))]); jest ono intensywnie wykorzystywane podczas symulowania dynamiki molekularnej. Metoda jest dość stabilna, gdyż prędkość jest domyślnie podana, a co za tym idzie trudniej jest prędkości i pozycji być niezsynchronizowanymi. (Na marginesie warto zauważyć, że podobne podejście wykorzystuje znany efekt tworzenia fal w wodzie). Metoda działa z uwagi na fakt, że 2x - x^* = x + (x - x^*) oraz x - x^* są przybliżeniem aktualnej prędkości (tak naprawdę jest to droga przebyta od ostatniego kroku). Metoda ta nie zawsze jest wyjątkowo dokładna (energia może np. opuścić system), ale za to szybka i stabilna. Zmieniając aktualizację reguły np. na taką:

\textbf{x}' = 1,99\textbf{x} - 0,99\textbf{x}^* + \textbf{a} \cdot \Delta t^2

możemy wprowadzić niewielką ilość oporu do systemu.

Pod koniec każdego etapu aktualna pozycja x każdej cząstki zostaje przechowana w odpowiedniej zmiennej \textbf{x}^*. Zauważ, że gdy manipulujemy wieloma cząstkami, przydatną optymalizacją może być po prostu zamiana tablicy wskaźników.

Wynikowy kod będzie wyglądać mniej więcej tak (klasa Vector3 powinna zawierać odpowiednie funkcje i przeciążenia operatorów manipulacji wektorów):

// przykładowy kod symulacji fizyki
class ParticleSystem
{
Vector3   m_x[NUM_PARTICLES];            // aktualne pozycje
Vector3   m_oldx[NUM_PARTICLES];        // poprzednie pozycje
Vector3   m_a[NUM_PARTICLES];            // obliczone siły
Vector3   m_vGravity;                           // grawitacja
float   m_fTimeStep;

public:
void   TimeStep();

private:
void   Verlet();
void   SatisfyConstraints();
void   AccumulateForces();

// konstruktory, inicjalizacje itp. pominięte
};

// krok całkowania Verleta
void ParticleSystem::Verlet()
{
for (int i = 0; i < NUM_PARTICLES; i++)
{
Vector3& x = m_x[i];
Vector3 temp = x;
Vector3& oldx = m_oldx[i];
Vector3& a = m_a[i];
x += x - oldx + a*fTimeStep*fTimeStep;
oldx = temp;
}
}

// funkcja ta powinna obliczać siły dla każdej cząstki
void ParticleSystem::AccumulateForces()
{
// wszystkie cząstki podlegają grawitacji
for (int i = 0; i < NUM_PARTICLES; i++)
m_a[i] = m_vGravity;
}

// tu uwzględniamy ograniczenia
void ParticleSystem::SatisfyConstraints()
{
// ta funkcja jest na razie nieistotna
}

void ParticleSystem::TimeStep()
{
AccumulateForces();
Verlet();
SatisfyConstraints();
}

Powyższy kod został napisany dla przejrzystości, a nie szybkości działania. Jedną z optymalizacji byłoby użycie tablic typu float zamiast Vector3 dla reprezentacji stanów. Mogłoby to również ułatwić implementację na procesorach wektorowych.

Zaprezentowana metoda prawdopodobnie nie wygląda jeszcze bardzo przełomowo. Jednakże korzyści powinny być wyraźne, gdy zaczniemy używać ograniczeń i przełączymy się na ciała sztywne. Zauważymy wtedy, jak powyższa metoda prowadzi do zwiększenia stabilności i zmniejszenia ilości obliczeń w porównaniu do innych metod.

Spróbuj np. ustawić \textbf{a} = (0,0,1) i użyć jako warunku początkowego \textbf{x} = (0,0,0) i \textbf{x}^* = (0,0,0), a następnie wykonać ręcznie kilka powtórzeń, aby zobaczyć, co się dzieje.

Dodaj komentarz