Numeryczne rozwiązywanie równań różniczkowych – metody Eulera i Rungego – Kutty

Mgr inż.

Sebastian Dziarmaga-Działyński

Wirtualna kancelaria korepetytorska i konsultacyjna – wszechwiedza.pl

Równania różniczkowe opisują zmiany zachodzące w czasie lub przestrzeni. Występują między innymi w mechanice, elektrotechnice, termodynamice, biologii, ekonomii, automatyce i teorii sterowania. Pozwalają modelować ruch ciał, przepływ ciepła, drgania, rozwój populacji, przebieg reakcji chemicznych oraz działanie obwodów elektrycznych.

Nie każde równanie różniczkowe można jednak rozwiązać analitycznie. Nawet jeżeli rozwiązanie dokładne istnieje, jego wyznaczenie może być bardzo trudne albo otrzymana postać może być mało użyteczna w praktycznych obliczeniach.

W takich sytuacjach stosuje się numeryczne metody rozwiązywania równań różniczkowych. Zamiast szukać jawnego wzoru funkcji, wyznaczamy jej przybliżone wartości w kolejnych punktach siatki.

W artykule omówimy między innymi:

• jawną i niejawną metodę Eulera,
• ulepszoną metodę Eulera, nazywaną metodą Heuna,
• metodę punktu środkowego,
• klasyczną metodę Rungego – Kutty czwartego rzędu,
• metody adaptacyjne,
• metody wielokrokowe,
• stabilność i równania sztywne.

Podstawowe znaczenie algorytmu numerycznego, zbieżności, stabilności i błędu zostało przedstawione w artykułach Metody numeryczne – czym są i do czego służą? oraz Błędy obliczeń numerycznych – rodzaje, źródła i ocena dokładności.

Zagadnienie początkowe

Rozważmy równanie różniczkowe pierwszego rzędu:

\[ y'(t)=f(t,y). \]

Samo równanie może mieć całą rodzinę rozwiązań. Aby wybrać jedną konkretną funkcję, podajemy warunek początkowy:

\[ y(t_0)=y_0. \]

Otrzymujemy w ten sposób zagadnienie początkowe, nazywane także zagadnieniem Cauchy’ego:

\[ \begin{cases} y'(t)=f(t,y),\\[2mm] y(t_0)=y_0. \end{cases} \]

Funkcja \(f(t,y)\) określa nachylenie rozwiązania w punkcie \((t,y)\). Jeżeli znamy aktualny stan układu, równanie różniczkowe mówi nam, w jakim kierunku rozwiązanie powinno się lokalnie zmieniać.

Rozwiązanie analityczne a rozwiązanie numeryczne

W rozwiązaniu analitycznym próbujemy wyznaczyć funkcję:

\[ y=y(t), \]

która spełnia równanie różniczkowe i warunek początkowy na całym rozpatrywanym przedziale.

Przykładowo dla zagadnienia:

\[ y'=-2y, \qquad y(0)=1 \]

rozwiązaniem dokładnym jest funkcja:

\[ y(t)=e^{-2t}. \]

Rozwiązanie numeryczne ma zwykle inną postać. Wybieramy punkty:

\[ t_0,\ t_1,\ t_2,\ldots,t_N \]

i obliczamy wartości:

\[ y_0,\ y_1,\ y_2,\ldots,y_N, \]

które powinny przybliżać dokładne rozwiązanie:

\[ y_n\approx y(t_n). \]

Podstawowym wynikiem obliczeń jest więc zbiór punktów:

\[ (t_0,y_0),\ (t_1,y_1),\ldots,(t_N,y_N). \]

Punkty można następnie połączyć odcinkami albo zastosować interpolację, aby otrzymać przybliżony wykres ciągły. Nadal jednak nie jest to na ogół symboliczny wzór funkcji, lecz jej numeryczna reprezentacja.

Siatka obliczeniowa i krok metody

W najprostszym przypadku stosujemy stały krok:

\[ h>0. \]

Kolejne punkty siatki wyznaczamy ze wzoru:

\[ t_{n+1}=t_n+h. \]

Stąd:

\[ t_n=t_0+nh. \]

Jeżeli chcemy obliczać rozwiązanie do chwili \(T\), liczba kroków jest w przybliżeniu równa:

\[ N=\frac{T-t_0}{h}. \]

Im mniejszy jest krok \(h\), tym gęściej rozmieszczone są punkty rozwiązania i zazwyczaj tym większa jest dokładność. Jednocześnie rośnie jednak liczba działań, koszt obliczeń oraz wpływ błędów zaokrągleń.

W metodach adaptacyjnych krok nie jest stały. Zamiast \(h\) stosuje się kolejne wartości:

\[ h_0,\ h_1,\ h_2,\ldots, \]

dostosowywane do lokalnego zachowania rozwiązania.

Całkowa interpretacja jednego kroku

Z równania:

\[ y'(t)=f(t,y(t)) \]

wynika:

\[ y(t_{n+1})-y(t_n) = \int_{t_n}^{t_{n+1}} f(t,y(t))\,dt. \]

A zatem:

\[ y(t_{n+1}) = y(t_n) + \int_{t_n}^{t_{n+1}} f(t,y(t))\,dt. \]

Wiele metod rozwiązywania równań różniczkowych można rozumieć jako różne sposoby przybliżania tej całki.

• Euler jawny wykorzystuje wartość funkcji na początku kroku.
• Euler niejawny wykorzystuje wartość na końcu kroku.
• Metoda trapezowa uśrednia wartości na początku i końcu.
• Metody Rungego – Kutty obliczają kilka pośrednich ocen nachylenia.

Związek ten łączy numeryczne rozwiązywanie równań różniczkowych z całkowaniem numerycznym.

Pole kierunków

W każdym punkcie płaszczyzny \((t,y)\) wartość:

\[ f(t,y) \]

określa nachylenie krzywej rozwiązania. Jeżeli w wielu punktach narysujemy krótkie odcinki o takim nachyleniu, otrzymamy pole kierunków.

Krzywe rozwiązań powinny w każdym swoim punkcie przebiegać zgodnie z lokalnym kierunkiem pola. Metody numeryczne można interpretować jako wędrówkę po tym polu: od aktualnego punktu przechodzimy do następnego, korzystając z jednej lub kilku ocen kierunku.

Jawna metoda Eulera

Jawna metoda Eulera jest najprostszą metodą jednokrokową. W aktualnym punkcie:

\[ (t_n,y_n) \]

obliczamy nachylenie:

\[ f(t_n,y_n). \]

Następnie przesuwamy się o krok \(h\) wzdłuż prostej o tym nachyleniu:

\[ y_{n+1} = y_n+h f(t_n,y_n). \]

Metodę tę można otrzymać z rozwinięcia Taylora:

\[ y(t_n+h) = y(t_n) + hy'(t_n) + \frac{h^2}{2}y''(\xi_n), \]

gdzie \(\xi_n\in(t_n,t_{n+1})\). Po pominięciu składnika zawierającego \(h^2\) oraz wykorzystaniu równania \(y'=f(t,y)\) otrzymujemy wzór Eulera.

Algorytm jawnej metody Eulera

1. Ustal punkt początkowy \((t_0,y_0)\).
2. Wybierz długość kroku \(h\).
3. Oblicz nachylenie \(f(t_n,y_n)\).
4. Wyznacz:
   \[
   y_{n+1}=y_n+h f(t_n,y_n).
   \]
5. Przyjmij:
   \[
   t_{n+1}=t_n+h.
   \]
6. Powtarzaj obliczenia do osiągnięcia końca przedziału.

Przykład jawnej metody Eulera

Rozważmy zagadnienie:

\[ y'=-2y, \qquad y(0)=1. \]

Przyjmujemy krok:

\[ h=0{,}25. \]

Wzór Eulera przyjmuje postać:

\[ y_{n+1} = y_n+0{,}25(-2y_n) = 0{,}5y_n. \]

Kolejne wartości wynoszą więc:

\[ y_0=1, \qquad y_1=0{,}5, \qquad y_2=0{,}25, \qquad y_3=0{,}125, \qquad y_4=0{,}0625. \]

Dokładna wartość w chwili \(t=1\) wynosi:

\[ y(1)=e^{-2}\approx0{,}1353352832. \]

Euler jawny daje:

\[ y_4=0{,}0625, \]

więc błąd jest wyraźnie widoczny. Zmniejszenie kroku poprawiłoby dokładność.

Dokładność jawnej metody Eulera

Lokalny błąd obcięcia, czyli błąd powstający podczas jednego kroku rozpoczętego od dokładnej wartości, jest rzędu:

\[ O(h^2). \]

Błędy poszczególnych kroków kumulują się. Globalny błąd rozwiązania na skończonym przedziale jest rzędu:

\[ O(h). \]

Dwukrotne zmniejszenie kroku powinno więc w typowej sytuacji zmniejszyć główny składnik błędu globalnego około dwa razy.

Niejawna metoda Eulera

W niejawnej metodzie Eulera wykorzystujemy nachylenie w punkcie końcowym kroku:

\[ y_{n+1} = y_n+ h f(t_{n+1},y_{n+1}). \]

Nowa wartość \(y_{n+1}\) występuje po obu stronach równania. Metoda nie podaje więc wyniku bezpośrednio i w każdym kroku trzeba rozwiązać równanie algebraiczne.

Dla funkcji liniowej może być to proste. Dla równania:

\[ y'=-2y \]

otrzymujemy:

\[ y_{n+1} = y_n-2h y_{n+1}. \]

Stąd:

\[ (1+2h)y_{n+1}=y_n, \]

czyli:

\[ y_{n+1} = \frac{y_n}{1+2h}. \]

Dla \(h=0{,}25\):

\[ y_{n+1} = \frac23y_n. \]

Kolejne przybliżenia wynoszą:

\[ 1,\quad 0{,}6666666667,\quad 0{,}4444444444,\quad 0{,}2962962963,\quad 0{,}1975308642. \]

Niejawna metoda Eulera ma, podobnie jak wariant jawny, globalny błąd rzędu \(O(h)\). Jej najważniejszą zaletą nie jest wyższy rząd dokładności, lecz znacznie lepsza stabilność.

Rozwiązywanie równania wewnątrz kroku

Dla nieliniowej funkcji \(f\) zależność:

\[ y_{n+1} - y_n - h f(t_{n+1},y_{n+1}) = 0 \]

jest równaniem nieliniowym względem \(y_{n+1}\). Można je rozwiązywać metodą Newtona lub innym algorytmem opisanym w artykule Numeryczne rozwiązywanie równań nieliniowych.

W takim przypadku pojedynczy krok metody niejawnej zawiera wewnętrzną iterację z własnym kryterium stopu. Jest to istotna różnica w porównaniu z metodą jawną.

Ulepszona metoda Eulera – metoda Heuna

Jawna metoda Eulera korzysta tylko z nachylenia na początku kroku. Metoda Heuna oblicza dodatkowo przewidywane nachylenie na jego końcu, a następnie uśrednia oba kierunki.

Jest to jawny schemat typu predyktor – korektor.

Predykcja

Najpierw wykonujemy zwykły krok Eulera:

\[ \widetilde y_{n+1} = y_n+ h f(t_n,y_n). \]

Wartość \(\widetilde y_{n+1}\) jest przewidywanym położeniem na końcu kroku.

Korekcja

Obliczamy nachylenie w przewidzianym punkcie:

\[ f(t_{n+1},\widetilde y_{n+1}), \]

a następnie stosujemy średnią nachyleń:

\[ y_{n+1} = y_n+ \frac h2 \left[ f(t_n,y_n) + f(t_{n+1},\widetilde y_{n+1}) \right]. \]

Przykład metody Heuna

Dla równania:

\[ y'=-2y, \qquad y_0=1, \qquad h=0{,}25 \]

predykcja daje:

\[ \widetilde y_1 = 1+0{,}25(-2) = 0{,}5. \]

Nachylenia na początku i w punkcie przewidywanym wynoszą:

\[ k_1=-2, \qquad k_2=-2\cdot0{,}5=-1. \]

Po korekcji:

\[ y_1 = 1+ \frac{0{,}25}{2} (-2-1) = 0{,}625. \]

Kolejne przybliżenia są równe:

\[ 1,\quad 0{,}625,\quad 0{,}390625,\quad 0{,}244140625,\quad 0{,}1525878906. \]

Metoda Heuna ma lokalny błąd rzędu:

\[ O(h^3), \]

oraz globalny błąd rzędu:

\[ O(h^2). \]

Jawna metoda punktu środkowego

Inną metodą drugiego rzędu jest jawna metoda punktu środkowego. Najpierw obliczamy przybliżenie rozwiązania w połowie kroku:

\[ y_{n+\frac12} = y_n+ \frac h2f(t_n,y_n). \]

Następnie wykorzystujemy nachylenie w tym punkcie:

\[ y_{n+1} = y_n+ h f\left( t_n+\frac h2, y_{n+\frac12} \right). \]

Metoda Heuna i metoda punktu środkowego są różnymi metodami Rungego – Kutty drugiego rzędu. W obu wykonuje się dwie oceny funkcji w jednym kroku, lecz w innych punktach.

Metody Rungego – Kutty

Metody Rungego – Kutty zwiększają dokładność bez jawnego obliczania wyższych pochodnych rozwiązania. W jednym kroku wykonują kilka ocen funkcji \(f(t,y)\), a następnie łączą je w odpowiednio dobraną średnią ważoną.

Metoda Eulera może być traktowana jako metoda Rungego – Kutty pierwszego rzędu. Heun i metoda punktu środkowego należą do metod drugiego rzędu.

Klasyczna metoda Rungego – Kutty czwartego rzędu – RK4

Najbardziej znanym schematem jest klasyczna metoda Rungego – Kutty czwartego rzędu, oznaczana jako RK4.

W każdym kroku obliczamy cztery nachylenia.

Nachylenie na początku kroku:

\[ k_1=f(t_n,y_n). \]

Pierwsza ocena w połowie kroku:

\[ k_2= f\left( t_n+\frac h2, y_n+\frac h2k_1 \right). \]

Druga ocena w połowie kroku:

\[ k_3= f\left( t_n+\frac h2, y_n+\frac h2k_2 \right). \]

Nachylenie na końcu kroku:

\[ k_4= f(t_n+h,y_n+hk_3). \]

Nową wartość obliczamy ze wzoru:

\[ y_{n+1} = y_n+ \frac h6 \left( k_1+2k_2+2k_3+k_4 \right). \]

Nachylenia \(k_2\) i \(k_3\), obliczane w środku kroku, otrzymują dwukrotnie większe wagi niż wartości na początku i końcu.

Przykład jednego kroku RK4

Dla równania:

\[ y'=-2y, \qquad y_0=1, \qquad h=0{,}25 \]

otrzymujemy:

\[ k_1=-2, \]

\[ k_2= -2\left(1+\frac{0{,}25}{2}(-2)\right) =-1{,}5, \]

\[ k_3= -2\left(1+\frac{0{,}25}{2}(-1{,}5)\right) =-1{,}625, \]

\[ k_4= -2\left(1+0{,}25(-1{,}625)\right) =-1{,}1875. \]

Stąd:

\[ y_1 = 1+ \frac{0{,}25}{6} \left[ -2+2(-1{,}5)+2(-1{,}625)-1{,}1875 \right] \]

\[ y_1 \approx 0{,}6067708333. \]

Dokładna wartość wynosi:

\[ y(0{,}25)=e^{-0{,}5} \approx 0{,}6065306597. \]

Już po jednym kroku wynik jest bardzo dokładny.

Dokładność RK4

Lokalny błąd obcięcia klasycznej metody RK4 jest rzędu:

\[ O(h^5), \]

a błąd globalny:

\[ O(h^4). \]

Dwukrotne zmniejszenie kroku powinno więc w typowej sytuacji zmniejszyć główny składnik błędu globalnego około:

\[ 2^4=16 \]

razy.

Porównanie metod na wspólnym przykładzie

Porównajmy wyniki dla zagadnienia:

\[ y'=-2y, \qquad y(0)=1, \]

przy kroku:

\[ h=0{,}25. \]

Błędy bezwzględne w chwili \(t=1\) wynoszą w przybliżeniu:

RK4 wykonuje więcej obliczeń w jednym kroku, ale uzyskuje znacznie większą dokładność.

Błąd lokalny i globalny

Przy analizie metod rozwiązywania równań różniczkowych należy rozróżnić dwa podstawowe rodzaje błędów.

Lokalny błąd obcięcia

Jest to błąd wykonany w pojedynczym kroku przy założeniu, że rozpoczynamy go od dokładnej wartości:

\[ y_n=y(t_n). \]

Błąd globalny

Jest to różnica pomiędzy dokładnym rozwiązaniem i wartością otrzymaną po wielu krokach:

\[ e_n=y(t_n)-y_n. \]

Błąd globalny zawiera wpływ błędów powstałych we wszystkich wcześniejszych krokach.

Czy istnieje kryterium stopu?

W metodach rozwiązywania równania:

\[ f(x)=0 \]

kolejne iteracje prowadzi się zazwyczaj do chwili, gdy przybliżenie pierwiastka spełni zadane kryterium dokładności.

Przy numerycznym rozwiązywaniu równania różniczkowego sytuacja jest inna. Najczęściej z góry ustalamy przedział:

\[ [t_0,T] \]

i obliczamy rozwiązanie aż do osiągnięcia chwili końcowej:

\[ t_N=T. \]

Nie próbujemy uzyskać jednego punktu granicznego przez nieskończone poprawianie przybliżenia. Każdy krok wyznacza kolejny fragment trajektorii.

Obliczenia mogą zakończyć się również wtedy, gdy:

• rozwiązanie osiągnie określoną wartość,
• zajdzie zdarzenie, na przykład \(y(t)=0\),
• trajektoria opuści dopuszczalny obszar,
• zostanie osiągnięta maksymalna liczba kroków,
• rozwiązanie zbliży się do stanu ustalonego,
• długość kroku spadnie poniżej dopuszczalnej wartości.

W metodach adaptacyjnych tolerancja błędu służy przede wszystkim do przyjmowania lub odrzucania pojedynczego kroku oraz doboru jego długości, a nie do określania końca całej trajektorii.

Metody adaptacyjne Rungego – Kutty

Stały krok nie zawsze jest efektywny. Na fragmentach, gdzie rozwiązanie zmienia się powoli, można stosować duży krok. W miejscach gwałtownej zmiany potrzebny jest krok mniejszy.

Popularne metody adaptacyjne wykorzystują dwie formuły Rungego – Kutty różnych rzędów, obliczane przy użyciu wspólnych wartości pośrednich. Przykładami są:

• metoda Rungego – Kutty – Fehlberga,
• metoda Dormanda – Prince’a,
• schematy RK45.

Jeżeli w jednym kroku otrzymano dwa przybliżenia:

\[ y_{n+1}^{(p)} \qquad\text{oraz}\qquad y_{n+1}^{(p+1)}, \]

ich różnica służy jako oszacowanie błędu:

\[ e_n \approx \left| y_{n+1}^{(p+1)} - y_{n+1}^{(p)} \right|. \]

Jeżeli błąd przekracza tolerancję, krok zostaje odrzucony i wykonany ponownie z mniejszym \(h\). Jeżeli jest znacznie mniejszy od tolerancji, kolejny krok można zwiększyć.

W efekcie punkty rozwiązania mogą być rozmieszczone nierównomiernie:

\[ t_{n+1}=t_n+h_n. \]

Stabilność metody numerycznej

Wysoki rząd dokładności nie gwarantuje poprawnego zachowania metody dla dowolnego kroku. Algorytm może być dokładny dla małego \(h\), ale niestabilny po jego zwiększeniu.

Do analizy stabilności stosuje się równanie testowe:

\[ y'=\lambda y, \qquad \operatorname{Re}(\lambda)<0. \]

Rozwiązanie dokładne:

\[ y(t)=y_0e^{\lambda t} \]

zanika do zera.

Stabilność jawnej metody Eulera

Euler jawny daje:

\[ y_{n+1} = (1+h\lambda)y_n. \]

Aby kolejne wartości zanikały, potrzebny jest warunek:

\[ |1+h\lambda|<1. \]

Dla równania:

\[ y'=-15y \]

otrzymujemy:

\[ |1-15h|<1, \]

czyli:

\[ 0<h<\frac{2}{15} \approx0{,}1333. \]

Dla \(h=0{,}2\):

\[ 1-15h=-2. \]

Kolejne wartości wynoszą:

\[ 1,\ -2,\ 4,\ -8,\ 16,\ldots \]

Metoda oscyluje i rozbiega się, mimo że rozwiązanie dokładne bardzo szybko zanika.

Stabilność niejawnej metody Eulera

Euler niejawny daje:

\[ y_{n+1} = \frac{y_n}{1-h\lambda}. \]

Dla \(\lambda=-15\) oraz \(h=0{,}2\):

\[ y_{n+1} = \frac{y_n}{4}. \]

Kolejne wartości wynoszą:

\[ 1,\quad 0{,}25,\quad 0{,}0625,\quad 0{,}015625,\ldots \]

Rozwiązanie numeryczne zachowuje stabilny zanik.

Stabilność nie jest tym samym co dokładność. Metoda niejawna może zachowywać właściwy charakter rozwiązania przy dużym kroku, ale jej wartości nadal mogą wyraźnie różnić się od rozwiązania dokładnego.

Równania sztywne

Równanie lub układ równań nazywa się sztywnym, gdy w rozwiązaniu występują procesy zachodzące w bardzo różnych skalach czasowych. Niektóre składniki zanikają bardzo szybko, podczas gdy inne zmieniają się powoli.

Jawna metoda może wówczas wymagać niezwykle małego kroku wyłącznie ze względu na stabilność, mimo że dokładny przebieg rozwiązania jest już prawie gładki.

W takich problemach często stosuje się:

• niejawną metodę Eulera,
• niejawną metodę trapezową,
• metody Adamsa – Moultona,
• wsteczne formuły różnicowe BDF,
• niejawne metody Rungego – Kutty.

Niejawna metoda trapezowa

Metoda trapezowa wykorzystuje średnią nachyleń na początku i końcu kroku:

\[ y_{n+1} = y_n+ \frac h2 \left[ f(t_n,y_n) + f(t_{n+1},y_{n+1}) \right]. \]

W przeciwieństwie do metody Heuna wartość \(y_{n+1}\) wewnątrz funkcji nie jest tylko predykcją, lecz poszukiwanym wynikiem. Jest to więc metoda niejawna, wymagająca rozwiązania równania w każdym kroku.

Metoda trapezowa ma globalny błąd rzędu \(O(h^2)\) i bardzo dobre właściwości stabilnościowe.

Metody wielokrokowe

Metody Eulera i Rungego – Kutty są metodami jednokrokowymi: do wyznaczenia \(y_{n+1}\) wykorzystują przede wszystkim informacje z bieżącego kroku.

Metody wielokrokowe korzystają również z wcześniej obliczonych wartości. Pozwala to ograniczyć liczbę nowych obliczeń funkcji \(f\).

Jawna metoda Adamsa – Bashfortha

Dwukrokowa metoda Adamsa – Bashfortha ma postać:

\[ y_{n+1} = y_n+ \frac h2 \left[ 3f(t_n,y_n) - f(t_{n-1},y_{n-1}) \right]. \]

Do rozpoczęcia obliczeń potrzebna jest dodatkowa wartość \(y_1\), wyznaczona na przykład metodą Rungego – Kutty.

Metody Adamsa – Moultona

Metody Adamsa – Moultona są schematami niejawnymi. Wartość funkcji w punkcie końcowym zależy w nich od poszukiwanego \(y_{n+1}\).

Często łączy się metodę Adamsa – Bashfortha jako predyktor z metodą Adamsa – Moultona jako korektor.

Wsteczne formuły różnicowe BDF

Metody BDF przybliżają pochodną za pomocą wartości rozwiązania z kilku kolejnych chwil. Dwukrokowy wzór BDF ma postać:

\[ \frac{ 3y_{n+1}-4y_n+y_{n-1} }{ 2h } = f(t_{n+1},y_{n+1}). \]

Metody te są niejawne i szczególnie użyteczne w rozwiązywaniu układów sztywnych.

Układy równań różniczkowych

Większość metod można bezpośrednio stosować do układów równań. Zamiast jednej liczby \(y\) mamy wektor:

\[ \mathbf y(t) = \begin{bmatrix} y_1(t)\\ y_2(t)\\ \vdots\\ y_m(t) \end{bmatrix}. \]

Układ zapisujemy jako:

\[ \mathbf y'(t) = \mathbf f(t,\mathbf y). \]

Jawna metoda Eulera przyjmuje wtedy postać:

\[ \mathbf y_{n+1} = \mathbf y_n+ h\mathbf f(t_n,\mathbf y_n). \]

Wzory Rungego – Kutty pozostają analogiczne, lecz wartości \(k_i\) są wektorami.

Równania wyższego rzędu

Równanie drugiego rzędu:

\[ y''=g(t,y,y') \]

można sprowadzić do układu dwóch równań pierwszego rzędu. Wprowadzamy:

\[ y_1=y, \qquad y_2=y'. \]

Wtedy:

\[ y_1'=y_2, \]

\[ y_2'=g(t,y_1,y_2). \]

Dzięki temu tymi samymi metodami można rozwiązywać między innymi równania ruchu mechanicznego:

\[ m y''+c y'+ky=F(t). \]

Zagadnienia brzegowe

Artykuł koncentruje się na zagadnieniach początkowych, w których wszystkie warunki są podane w jednym punkcie \(t_0\).

W zagadnieniach brzegowych warunki są zadane w różnych punktach, na przykład:

\[ y(a)=\alpha, \qquad y(b)=\beta. \]

Do ich rozwiązywania stosuje się między innymi:

• metodę strzałów,
• metody różnic skończonych,
• metody elementów skończonych,
• metody kolokacji.

W metodzie strzałów zagadnienie brzegowe sprowadza się do wielokrotnego rozwiązywania zagadnienia początkowego i dobierania brakującego warunku tak, aby spełnić warunek na drugim końcu przedziału.

Wpływ błędów zaokrągleń

Zmniejszanie kroku redukuje błąd obcięcia, ale zwiększa liczbę działań. Każde działanie jest wykonywane w skończonej precyzji, dlatego bardzo duża liczba kroków może powodować kumulowanie błędów zaokrągleń.

Znaczenie ograniczonej reprezentacji liczb zostało omówione w artykule Arytmetyka zmiennoprzecinkowa – cecha, mantysa i zapis binarny.

W praktyce krok powinien być wystarczająco mały, aby ograniczyć błąd metody, ale nie mniejszy, niż wymaga tego oczekiwana dokładność.

Porównanie najważniejszych metod

Którą metodę wybrać?

• Euler jawny jest odpowiedni do nauki, prostych symulacji i wstępnego sprawdzania modeli.
• Euler niejawny jest użyteczny, gdy najważniejsza jest stabilność, szczególnie dla równań sztywnych.
• Heun i metoda punktu środkowego zapewniają dobrą równowagę między prostotą i dokładnością.
• RK4 jest dobrym wyborem dla wielu niesztwnych zagadnień z ustalonym krokiem.
• RK45 i inne metody adaptacyjne sprawdzają się, gdy rozwiązanie zmienia się z różną szybkością.
• Metody BDF są często stosowane w dużych, sztywnych układach.

Praktyczne zasady obliczeń

• Sprawdź, czy rozwiązujesz zagadnienie początkowe, czy brzegowe.
• Ustal przedział obliczeń i oczekiwaną dokładność.
• Dobierz krok odpowiedni do szybkości zmian rozwiązania.
• Porównuj wyniki dla kroków \(h\) i \(h/2\).
• Nie utożsamiaj stabilności z dokładnością.
• Dla metod niejawnych kontroluj dokładność wewnętrznego rozwiązania równania.
• W metodach adaptacyjnych ustal tolerancję bezwzględną i względną.
• Kontroluj, czy rozwiązanie nie opuszcza dziedziny modelu.
• Uwzględniaj możliwość wystąpienia zdarzeń pomiędzy punktami siatki.
• Przy układach sztywnych nie zwiększaj bezkrytycznie kroku metody jawnej.
• Sprawdzaj zachowanie wielkości, które powinny być zachowane, na przykład energii lub masy.
• Porównuj wynik z rozwiązaniem dokładnym lub inną metodą, jeżeli jest to możliwe.

Numeryczne rozwiązanie równania różniczkowego nie jest zazwyczaj wzorem funkcji, lecz obliczoną krok po kroku trajektorią jej przybliżonych wartości.

Podsumowanie

Numeryczne rozwiązywanie równań różniczkowych polega najczęściej na wyznaczaniu kolejnych wartości rozwiązania w punktach siatki czasowej. Zamiast jawnej postaci funkcji otrzymujemy dyskretny zbiór punktów opisujących jej przybliżoną trajektorię.

Jawna metoda Eulera wykorzystuje nachylenie na początku kroku. Jest bardzo prosta, lecz ma dokładność pierwszego rzędu i może być niestabilna.

Niejawna metoda Eulera korzysta z nachylenia w punkcie końcowym. Wymaga rozwiązania równania w każdym kroku, ale ma znacznie lepsze właściwości stabilnościowe.

Metoda Heuna i metoda punktu środkowego wykonują po dwie oceny nachylenia i osiągają dokładność drugiego rzędu.

Klasyczna metoda RK4 wykorzystuje cztery nachylenia w jednym kroku. Jej globalny błąd jest rzędu \(h^4\), dzięki czemu zapewnia dużą dokładność bez konieczności obliczania wyższych pochodnych.

Metody adaptacyjne samodzielnie zmieniają długość kroku na podstawie oszacowania lokalnego błędu. Metody niejawne i wielokrokowe są szczególnie istotne w rozwiązywaniu równań sztywnych.

Obliczenia prowadzi się zwykle do końca ustalonego przedziału albo do wystąpienia określonego zdarzenia. Klasyczne kryterium stopu pojawia się głównie wewnątrz metod niejawnych, gdy w pojedynczym kroku trzeba iteracyjnie rozwiązać równanie algebraiczne.

Utworzono: 17.06.2026 | Zmodyfikowano: 26.07.2026

Powiązane artykuły

• Błędy obliczeń numerycznych – rodzaje, źródła i ocena dokładności 
• Metody numeryczne – czym są i do czego służą? 

Masz problem z tym tematem?

Wszechwiedza.pl pomaga zrozumieć matematykę, statystykę, ekonometrię, badania operacyjne, analizę danych, mechanikę, rachunkowość i wiele innych przedmiotów — spokojnie, konkretnie i krok po kroku. 

Zapytaj o pomoc