Metoda Dormanda–Prince’a

Istota. Metoda Dormanda–Prince’a 5(4) (DP5) jest adaptacyjną metodą Rungego–Kutty o rzędach 5 i 4, podobnie jak RKF45, lecz z innym doborem współczynników. Wariant Dormanda–Prince’a został zaprojektowany tak, aby rozwiązanie rzędu 5 (główne) miało szczególnie mały błąd globalny, a rozwiązanie rzędu 4 służyło jedynie do kontroli błędu lokalnego. Dla układu \( \dot{\mathbf{x}} = f(\mathbf{x}, t, \theta) \) metoda w jednym kroku oblicza wektory nachyleń \( \mathbf{k}_1, \dots, \mathbf{k}_7 \) według tablic Butchera, a następnie na ich podstawie wyznacza przybliżenie rzędu 5, oznaczane zwykle jako \( \mathbf{y}_5 \), oraz konstruuje „towarzyszące” przybliżenie rzędu 4 przy użyciu innych wag. Różnica między tymi dwoma przybliżeniami jest estymatą lokalnego błędu w danym kroku. Na podstawie tej estymaty dobierana jest nowa długość kroku czasowego \( \Delta t \), podobnie jak w metodzie RKF45: jeśli błąd jest mniejszy niż zadana tolerancja, krok jest akceptowany i można krok wydłużyć; jeśli błąd przekracza tolerancję, krok jest odrzucany, a metoda automatycznie zmniejsza krok czasowy.

Zastosowanie. Metoda Dormanda–Prince’a jest jednym z najczęściej stosowanych schematów adaptacyjnych w praktycznych solverach ODE (rodzina metod spokrewnionych z popularnym „ode45” w MATLAB-ie). W porównaniu z RKF45, DP5 kładzie większy nacisk na to, aby główne rozwiązanie rzędu 5 miało jak najmniejszy błąd globalny, co często przekłada się na lepszą jakość trajektorii przy podobnej liczbie kroków. Dzięki adaptacyjnemu doborowi kroku czasowego metoda dobrze sprawdza się przy wizualizacji złożonych trajektorii atraktorów, takich jak Lorenz czy Rössler, gdzie gęstość próbkowania musi automatycznie dostosowywać się do lokalnej dynamiki.

Implementacja. W dodatku Attractor Builder każdy krok metody Dormanda–Prince’a 5(4) realizuje funkcja step_dp5(rhs_func, p: Vector, params: dict, current_dt: float, tolerance: float). Metoda DP5, podobnie jak RKF45, wykorzystuje tzw. tablice Butchera — układ współczynników określających przesunięcia czasowe (macierz \( A \)) oraz wagi kombinacji liniowych (\( B \), \( B^* \)), z których otrzymywane są dwa przybliżenia rozwiązania: rząd 5 (główny) oraz wbudowany rząd 4. Pełną tablicę współczynników metody Dormanda–Prince’a 5(4) można znaleźć w (Hairer, Nørsett, Wanner 1993, Table 5.2, s. 178), natomiast ogólną definicję tablic Butchera opisano w Butcher (2008, rozdz. 5, s. 96). Pierwszym etapem jest obliczenie wektorów nachyleń \( \mathbf{k}_1, \dots, \mathbf{k}_6 \) zgodnie z macierzą \( A \): k1 = rhs_func(p, params) k2 = rhs_func(p + current_dt * A[1][0] * k1, params) k3 = rhs_func(p + current_dt * (A[2][0]*k1 + A[2][1]*k2), params) k4 = rhs_func(p + current_dt * (A[3][0]*k1 + A[3][1]*k2 + A[3][2]*k3), params) k5 = rhs_func(p + current_dt * (A[4][0]*k1 + A[4][1]*k2 + A[4][2]*k3 + A[4][3]*k4), params) k6 = rhs_func(p + current_dt * (A[5][0]*k1 + A[5][1]*k2 + A[5][2]*k3 + A[5][3]*k4 + A[5][4]*k5), params) Następnie z użyciem wag \( B \) wyznaczane jest przybliżenie rzędu 5: y5 = p + current_dt * (B[0]*k1 + B[1]*k2 + B[2]*k3 + B[3]*k4 + B[4]*k5 + B[5]*k6) Dodatkowo obliczane jest nachylenie k7 = rhs_func(y5, params), wykorzystywane do estymaty błędu przy użyciu wag \( B^* \). Wektor błędu lokalnego ma postać: all_k = [k1, k2, k3, k4, k5, k6, k7] error_vec = current_dt * sum([(b - bs) * k for b, bs, k in zip(B, B_star, all_k)], Vector()) error = error_vec.length co odpowiada różnicy pomiędzy głównym rozwiązaniem rzędu 5 a wbudowanym rozwiązaniem rzędu 4. Na tej podstawie obliczany jest optymalny krok czasowy: optimal_dt = 0.9 * current_dt * (tolerance / error)**0.2 Jeżeli error <= tolerance, krok jest akceptowany (nowy punkt trajektorii to y5). Jeśli błąd jest zbyt duży, krok zostaje odrzucony i wykonywana jest ponowna próba z krokiem optimal_dt. Podobnie jak w przypadku RKF45, funkcja step_dp5 zwraca trójkę: (nowy_punkt, proponowany_krok_dt, flaga_akceptacji).

Przykładem zastosowania może być układ Lorenza:

Równania:
dx/dt = a * (y - x)
dy/dt = x * (b - z) - y
dz/dt = x * y - c * z
Parametry:
| a = 10 | b = 28 | c = 8/3 |
Ustawienia symulacji:
Stan początkowy: x₀ = 0.01, y₀ = 0.01, z₀ = 0.01
Metoda: DP5
Tolerancja: 0.000001
min_dt: 0.000001
max_dt: 0.05
Zapisane kroki: 3
Burn-in: 1000
Skala: 1

Korzystając z funkcji Raw Data → Export w dodatku Attractor Builder otrzymujemy surowe dane w formacie CSV z kolumnami steps, dt, x, y, z. Dla powyższych ustawień pierwsze cztery wiersze pliku CSV wyglądają następująco:

Pierwszy wiersz odpowiada warunkowi początkowemu \( (x_0, y_0, z_0) = (0.01, 0.01, 0.01) \) i ma przypisany czasowy krok \( dt = 0 \). W kolejnych wierszach widać, jak metoda Dormanda–Prince’a adaptuje krok czasowy: po pierwszym kroku dt rośnie z 0.01 do około 0.029 (3 krok), co oznacza, że w tym fragmencie trajektorii ruch jest na tyle regularny, że można używać większego kroku bez utraty zadanej dokładności.

Źródła:

Butcher, J. C. (2008). Numerical Methods for Ordinary Differential Equations. Wiley.

Hairer, E., Nørsett, S., Wanner, G. (1993). Solving Ordinary Differential Equations I. Springer.