Przybliżenie linią styczną
Idea geometryczna
Równanie różniczkowe zwyczajne dy/dx = f(x,y) przypisuje nachylenie każdemu punktowi na płaszczyźnie (x,y) — pole kierunków. Prawdziwe rozwiązanie y(x) jest krzywą, która wszędzie podąża przypisanymi nachyleniami.
Metoda Eulera konwertuje ciągłe pole kierunków na dyskretny spacer:
> (xₙ, yₙ) → (xₙ + h, yₙ + h·f(xₙ, yₙ))
Z punktu (xₙ, yₙ), poruszaj się o odległość h wzdłuż linii stycznej. Osiągnij przybliżony następny punkt. Powtarzaj.
Błąd geometryczny: linia styczna w punkcie (xₙ, yₙ) ma nachylenie f(xₙ, yₙ), ale prawdziwa krzywa ma inne nachylenie w każdym punkcie wzdłuż interwału [xₙ, xₙ + h]. Krok Eulera używa nachylenia w lewym końcu punktu — 'nachylenie, które było.' Błąd na krok rośnie jako h².
Skumulowany błąd
Nad N krokami, aby osiągnąć stały punkt końcowy x = a, z h = a/N:
- Lokalny błąd obcięcia na krok: O(h²)
- Liczba kroków: N = a/h
- Błąd globalny: O(h²) × (a/h) = O(h) — dokładność pierwszego rzędu
Metoda Eulera jest pierwszego rzędu: zmniejszenie h o połowę zmniejsza błąd globalny o połowę.
Uruchamianie metody Eulera
Rozważ dy/dx = y, z warunkiem początkowym y(0) = 1. Prawdziwe rozwiązanie: y(x) = eˣ, więc y(1) = e ≈ 2,71828.
Zastosuj metodę Eulera z h = 0,5, od x = 0 do x = 1 (2 kroki):
Krok 1: y₁ = y₀ + h·f(x₀, y₀) = 1 + 0,5·(1) = 1,5. Nowy punkt: (0,5, 1,5).
Krok 2: y₂ = y₁ + h·f(x₁, y₁) = 1,5 + 0,5·(1,5) = 2,25. Nowy punkt: (1,0, 2,25).
Euler podaje 2,25 vs prawdziwą wartość 2,71828. Błąd: 0,468. Błąd względny: ~17%.
Wyprowadzanie regionu stabilności Eulera
Dla równania testowego dy/dx = λy (gdzie λ jest liczbą zespoloną), metoda Eulera daje:
> yₙ₊₁ = yₙ + h·λ·yₙ = yₙ·(1 + hλ)
Współczynnik wzmocnienia na krok: z = 1 + hλ.
Warunek stabilności: obliczone rozwiązanie pozostaje ograniczone wtedy i tylko wtedy, gdy |z| ≤ 1, tj. |1 + hλ| ≤ 1.
To jest warunek geometryczny na złożonej płaszczyźnie hλ: punkt hλ musi leżeć wewnątrz koła o promieniu 1 wyśrodkowanego w (-1, 0).
Region stabilności Eulera: { hλ ∈ ℂ : |1 + hλ| ≤ 1 }
Dla rzeczywistych, ujemnych λ (zanikającego ODE takiego jak dy/dx = -2y): hλ musi leżeć w interwale (-2, 0] na osi rzeczywistej. Z λ = -2 i h = 0,5: hλ = -1. To jest dokładnie na granicy stabilności — metoda jest marginalnie stabilna, co wyjaśnia jakościową porażkę w wcześniejszym przykładzie.
Z h = 1 i λ = -2: hλ = -2, wysyłając nas poza region stabilności. Rozwiązanie oscyluje z rosnącą amplitudą.
Znajdowanie granicy stabilności
Runge-Kutta 4 (RK4) ma większy region stabilności niż Euler, co jest jedną z przyczyn, dla którego jest preferowany dla większości problemów.
Dla rzeczywistych ujemnych λ, RK4 pozwala hλ aż do około -2,785 na osi rzeczywistej (vs limit Eulera -2).
Dla sztywnych równań z wartościami własnymi λ o bardzo różnych magnitudach — powiedzmy λ₁ = -1 i λ₂ = -1000 — stabilność wymaga, aby hλ₂ pozostało wewnątrz regionu. Dla RK4 na osi rzeczywistej: h·(-1000) ≥ -2,785, więc h ≤ 0,002785.
Ten mały rozmiar kroku, podyktowany sztywną wartością własną λ₂, czyni symulację drogą, nawet jeśli wolny komponent λ₁ mógłby używać h = 2.
Punkty stałe & baseny atrakcji
Metoda Eulera zastosowana do dy/dx = f(y) definiuje mapę dyskretną: yₙ₊₁ = g(yₙ) = yₙ + h·f(yₙ).
Punkt stały tej mapy: y taki, że g(y) = y. Dla Eulera na dy/dx = f(y), punkty stałe spełniają f(y) = 0 — równowaga ODE.
Stabilność punktu stałego: jeśli |g'(y)| < 1, pobliskie iteracje zbiegają się do y. Jeśli |g'(y*)| > 1, rozbiegają się.
g'(y) = 1 + h·f'(y). W punkcie stałym y: |1 + h·f'(y)| < 1 do stabilności.
To jest dokładnie warunek stabilności Eulera z λ = f'(y*) — linearyzacja ODE na równowadze.
Basen atrakcji: zbiór warunków początkowych, które zbiegają się do y* pod mapą Eulera. Dla systemów nieliniowych, granica basenu określa, gdzie symulacja będzie niezawodnie śledźić równowagę ODE vs rozbiec się do innego atraktora.
Pętla symulacji jest dyskretnym systemem dynamicznym. Jego jakościowe zachowanie — zbieżność, oscylacja, rozbieżność — zależy od wielkości kroku h w stosunku do geometrii pola kierunków ODE.
Łączenie geometrii z projektowaniem symulacji
Geometria symulacji numerycznej sprowadza się do trzech pytań:
1. Gdzie jest region stabilności? Dla Eulera: dysk |1 + hλ| ≤ 1. Większy dla RK4. Nieograniczony (całą lewą półpłaszczyznę) dla metod niejawnych.
2. Gdzie są wartości własne ODE? Wartości własne λ jakobianu f w każdym punkcie określają, który region stabilności musi zawierać hλ.
3. Co h utrzymuje wszystkie hλ wewnątrz regionu? Maksymalnie dozwolony h = (promień regionu stabilności) / max|λ|.
Dla sztywnych systemów: max|λ| jest ogromny, zmuszając mały h dla metod jawnych. Metody niejawne są kosztowne na krok, ale pozwalają na duży h — wymieniają koszt na krok za stabilność.
Wgląd Hamminga przetłumaczy: wybór metody numerycznej koduje zakład o geometrii spektrum wartości własnych ODE. Uczynienie tego zakładu wyraźnym jest pierwszą decyzją projektową w każdej symulacji.