Verlet, Leapfrog & RK4: Wybieranie Integratora Numerycznego

Co właściwie robi integrator

Prnajmniej dla większości symulacji na tej stronie rozwiązujemy tę samą problematykę: znasz siły działające na ciało w chwili bieżącej, a chcesz znać jego położenie w przyszłej chwili. Newton daje drugorzędną równanie różniczkowe (ODE), a = F(x)/m, a integrator numeryczny konwertuje je na regułę przesuwającą (x, v) o czas h. Wybór reguły decyduje, czy planeta pozostanie w orbitie przez godzinę czasu zegara lub spadnie w słońce do trzeciej minuty.

Dwa cechy są istotne i nie są tym samym. Dokładność to jak mała jest błąd po jednym kroku (błąd lokalny) lub po ustalonej liczbie kroków (błąd globalny). Stabilność to czy błąd pozostaje ograniczony, gdy kontynuujesz przesuwanie. Metoda może być dokładna i niestabilna, lub nieprecyzyjna i pięknie stabilna — Verlet w prędkościach jest drugiego rodzaju, a w fizyce to zwykle to, czego chcesz.

Euler niejawny: ten, który kłamię

Najprostsza schemat ocenia siłę tylko raz i przesuwa wszystko o tej siłce:

// jawny (przód) Euler — nie używaj do orbity v += a(x) * h; x += v * h; // tutaj v to stara v, jeśli zapiszesz to w tym porządku Jego błąd lokalny wynosi O(h²) i globalny O(h) — pierwszego rzędu. Gorsza jest sytuacja dla oscylatora: jest niezależnie niewygodny — amplituda rośnie na każdym kroku, niezależnie od tego jak mała h jest. Symulując sprężynę lub okrągłą orbitę z użyciem jawnego Eulera, energia wzrasta monotonicznie; planeta spirali zewnątrz. Zmniejszenie kroku tylko opóźni katastrofę, ale nie zapobiegnie jej.

// explicit (forward) Euler — do not use for orbits
v += a(x) * h;
x += v * h;        // v here is the OLD v if you write it in this order

Symplektyczne, i dlaczego to lepsze od prostych metod dokładności

Systemy Hamiltona — grawitacja, sprężyny, wahadła, dynamiczna molekuła — zachowują objętość przestrzeni fazowej. Symplektyczny integrator zachowuje dyskretną odpowiednicję tej własności, a konsekwencją jest dokładne zachowanie lekko perturbowanego „cienia” energii, które pozostaje w granicach O(h²) prawdziwej energii na zawsze. Praktycznie: energia symplektycznego obliczenia oscyluje wokół poprawnej wartości zamiast odchodzić od niej. Metoda nie-symplektyczna o wyższym porządku błąd będzie mniejszy podczas jednej orbity, ale ma secularny błąd po tysiącu orbitach.

Szybkościowa metoda Verleta i skok leaping

Metoda szybkościowa Verleta jest podstawą dla dynamiki molekularnej oraz rozwiązań ciał ruchomych w każdym silniku gier. Jest symplecticzna, czasowo odwracalna, drugiego rzędu i wymaga dokładnie jednego oceny siły na krok, jeśli przechowujemy przyspieszenie:

// Metoda szybkościowa Verleta — 2. rzędu, symplecticzna, 1 ocena siły na krok x += v * h + 0.5 * aOld * h * h; const aNew = a(x); // jedyna ocena siły v += 0.5 * (aOld + aNew) * h; aOld = aNew; Metoda skok leaping to taka sama metoda, tylko w odzieży innej: pozycje żyją na całkowitych krokach, a prędkości na półkrokach, i dwie

leap

// Velocity Verlet — 2nd order, symplectic, 1 force eval/step
x += v * h + 0.5 * aOld * h * h;
const aNew = a(x);                 // the only force evaluation
v += 0.5 * (aOld + aNew) * h;
aOld = aNew;

RK4: precyzyjny, nie symplektyczny

Klasyczne czwarte rzędu Runge–Kutta próbuje pochodnej cztery razy na krok i łączy je. Błąd globalny wynosi O(h⁴), więc jest znacznie precyzyjniejsze na krok niż Verlet — ale kosztuje to cztery oceny sił, a nie jest symplektyczne. Zastosuj RK4 do zamkniętej orbity i energia maleje wolno i monotonicznie: planeta spiraliem się zbliża. Dla chaotycznego systemu integrowanego na krótki horizont (dwupięciokąt, atraktor Lorenza) to jest dobre i precyzja jest zasługiwana; dla systemu słonecznego integrowanego przez milion lat to katastrofa.

const k1 = f(t, y);
const k2 = f(t + h/2, y + h/2 * k1);
const k3 = f(t + h/2, y + h/2 * k2);
const k4 = f(t + h,   y + h   * k3);
y += h/6 * (k1 + 2*k2 + 2*k3 + k4);   // 4 evaluations, O(h⁴), not symplectic

Wybór wielkości kroku

Krok jest ograniczony szybkim elementem w systemie. Dla oscylatora z częstotliwością kątową ω, metoda jawna drugiego rzędu potrzebuje około h·ω < 2 aby pozostawać stabilną, a chcesz znacznie niżej — reguła dobrej praktyki to 20 do 50 kroków na okres najstiffniejszego trybu przed tym, jak ruch wygląda prawidłowo. Dla systemów cząsteczkowych odpowiednikiem tego ograniczenia jest warunek CFL: żadna cząsteczka nie może przekroczyć więcej niż ułamek promienia interakcji w jednym kroku.

Dwa praktyczne zwyczaje. Pierwszy, odseparuj krok fizyki od częstotliwości ramki: uruchom stały h (np. 1/240 s) w pętli akumulatora i interpoluj dla renderowania, inaczej utracony ramek cicho zmienia zachowanie integratora. Drugi, narysuj całą energię. Jest to najtańszy detektor błędu — monotoniczne wzrostu oznacza jawny Euler gdzieś tam, monotoniczne maleńkość oznacza metodę nie-symplektyczną lub za dużo tłumienia, a ograniczone drgania oznaczają, że jesteś w porządku.

Gdzie używamy tych metod

Ten artykuł omawia trzy metody symulacji: Verlet, Leapfrog i RK4.

orbits, N-body, springs, molecular dynamics  → Velocity Verlet / Leapfrog
cloth, ropes, soft bodies                    → position Verlet + constraints
SPH fluids                                   → Leapfrog with a CFL-limited step
double pendulum, Lorenz, chaotic systems     → RK4 (short horizon, high accuracy)

Często zadawane pytania

Wypróbuj na żywo

Wszystko powyżej działa bezpośrednio w Twojej przeglądarce — otwórz Spring-Mass Cloth Simulation i zmieniaj parametry podczas działania. Nic nie jest instalowane ani przesyłane na serwer, cały model działa w jednej karcie.