Równania różniczkowe

Po metodzie Eulera (która jest metodą pierwszego rzędu), w analizie numerycznej równań różniczkowych zwyczajnych (ang. ODE Ordinary Differential Equations) stosuje się bardziej zaawansowane algorytmy. Zapewniają one znacznie wyższą dokładność i stabilność bez konieczności drastycznego zmniejszania odstępu \(h\).

Oto krótkie zestawienie metod rozwiązywania równań różniczkowych:

Metody jawne wyższych rzędów (Jednokrokowe)

To bezpośredni, naturalni następcy metody Eulera. Też liczą wartość w punkcie \(x + h\) na podstawie punktu poprzedniego, ale robią to z większą precyzją.

• Zmodyfikowana metoda Eulera (metoda Heuna): Metoda II rzędu dokładności. Najpierw wylicza się punkt próbny metodą Eulera, a następnie koryguje wynik za pomocą średniej arytmetycznej z pochodnych na początku i na końcu kroku.
• Ulepszona metoda Eulera (metoda punktu środkowego): Również II rząd. Wyznacza wartość pochodnej dokładnie w połowie kroku \({(x + \frac{h}{2})}\) i na tej podstawie wykonuje ostateczny skok.
• Metoda Rungego-Kutty IV rzędu (RK4): Absolutny klasyk i złoty standard w fizyce komputerowej oraz inżynierii. Liczy cztery współczynniki próbne \({k_1, k_2, k_3, k_4}\) wewnątrz jednego kroku i wyciąga z nich średnią ważoną. Lokalny błąd obcięcia wynosi aż \(h^{5}\), a błąd globalny \(h^{4}\).

Metody adaptacyjne (Z automatycznym doborem kroku)

• W tych algorytmach komputer sam decyduje, jak duży powinien być krok \(h\). Jeśli funkcja mocno skręca, krok maleje; jeśli jest liniowa, krok rośnie.
• Metoda Rungego-Kutty-Fehlberga (RKF45): Łączy metodę IV i V rzędu, porównując ich wyniki, aby oszacować bieżący błąd i na tej podstawie dopasować wielkość kroku \(h\).
• Metoda Casha-Karpa / Dormanda-Prince'a (RK54): Nowocześniejsze i bardziej wydajne odmiany metody adaptacyjnej, powszechnie stosowane w profesjonalnych pakietach obliczeniowych.

Metody niejawne (Dla równań sztywnych)

Niezbędne do rozwiązywania układów, w których procesy zachodzą w skrajnie różnych skalach czasowych (tzw. stiff equations), gdzie klasyczny Euler czy RK4 błyskawicznie "wylatują w nieskończoność".

• Niejawna metoda Eulera (wsteczna): Wzór wymaga wyznaczenia wartości \({f(x+h, y_{nowe})}\), co prowadzi do konieczności rozwiązywania równania algebraicznego w każdym kroku. Oferuje jednak niemal absolutną stabilność numeryczną.
• Metoda trapezów: Niejawna metoda II rzędu, bardzo stabilna, często stosowana w symulacjach obwodów elektrycznych.

Metody wielokrokowe (Pamiętające historię)

Zamiast liczyć wszystko od nowa w każdym kroku przy użyciu punktów próbnych (jak przy metodzie Rungego-Kutty), te metody wykorzystują wartości wyliczone w kilku poprzednich krokach.

• Metody Adamsa-Bashfortha: Jawne metody wielokrokowe.
• Metody Adamsa-Moultona: Niejawne metody wielokrokowe o bardzo wysokiej stabilności.
• Metody wstecznego różniczkowania (BDF): Najpotężniejsze narzędzia do skrajnie sztywnych układów równań.

Jak widzisz, temat jest bardzo rozległy i wymagałby osobnego artykułu. Tutaj opiszemy klasyczną metodę Rungego-Kutty czwartego rzędu (ang. Classical fourth-order Runge-Kutta method).

Krótka historia metody

Klasyczna metoda Rungego-Kutty nie została stworzona przez jednego uczonego, lecz jest efektem pracy dwóch wybitnych niemieckich matematyków na przełomie XIX i XX wieku. Podwaliny pod algorytm położył w 1895 roku Carl Runge. Zaproponował on rozbudowanie prostych metod numerycznych (takich jak metoda Eulera) o dodatkowe punkty próbne pomiędzy krokami obliczeniowymi, co miało drastycznie zwiększyć dokładność otrzymywanych wyników. Pomysł ten w 1901 roku udoskonalił i uogólnił Martin Kutta, formułując ostateczny, genialny wariant czwartego rzędu (znany dziś na całym świecie jako RK4).Opracowany przez nich algorytm okazał się tak potężny i stabilny, że ponad sto lat później — w erze superkomputerów — wciąż stanowi złoty standard i absolutną podstawę w fizyce komputerowej, inżynierii oraz symulacjach mechaniki nieba.

Rdzeń matematyczny RK4

W metodzie Rungego-Kutty zamiast jednej stycznej (jak w metodzie Eulera) liczymy cztery próbne współczynniki nachylenia \({(k_1, k_2, k_3, k_4)}\), a nowe \(y\) wyznaczamy jako średnią ważoną z tych czterech wartości.

Dla równania \({\frac{dy}{dx} = f(x,y)}\) i odstępu \(h\), algorytm przejścia do kolejnego punktu wygląda następująco:

1. Współczynnik \(k_1\) – nachylenie na początku kroku (to jest dokładnie to, co robi metoda Eulera):

2. Współczynnik \(k_2\) – nachylenie w połowie kroku, przy użyciu wyliczonego \(k_1\):

3. Współczynnik \(k_3\) – ponowne nachylenie w połowie kroku, ale przy użyciu dokładniejszego współczynnika \(k_2\):

4. Współczynnik \(k_4\) – nachylenie na końcu kroku, przy użyciu \(k_3\):

5. Ostateczne wyznaczenie nowego \(y\):

6. Przesunięcie \(x\) na kolejną pozycję:

Algorytm metody RK4

Specyfikacja problemu

Dane wejściowe

Dane wyjściowe

Zmienne pomocnicze

Lista kroków

Poniższe programy rozwiązują równanie różniczkowe \({\frac{dy}{dx} = x + y}\) z warunkiem  początkowym \({y(0) = 1}\) na przedziale od \({x = 0}\) do \({x = 2{,}5}\). Jako odstęp \(h\) wybieramy \(0{,}125\) (taka liczba jest przez komputer pamiętana dokładnie bez błędu zaokrągleń zmiennoprzecinkowych). Wyjaśnienia są w komentarzach.

Dla porównania programy wyświetlają również wartość dokładną funkcji \(y(x)\). Rozwiązaniem analitycznym (dokładnym) równania różniczkowego \({\frac{dy}{dx} = x + y}\) z warunkiem początkowym \({y(0) = 1}\) jest funkcja:

Porównaj wyniki z poprzednią metodą Eulera.

// Równania różniczkowe
// metoda Rungego-Kutty
// Metody numeryczne 0086
// --------------------------
// (C)2026 mgr Jerzy Wałaszek

#include <windows.h>
#include <iostream>
#include <iomanip>
#include <cmath>

using namespace std;

// Stałe
//------
// Przedział (start,koniec)
const double x_s = 0.000;
const double x_k = 2.500;
// Warunek początkowy
const double y_s = 1.000;
// Odstęp na osi OX
const double h   = 0.125;

// Funkcja prawej strony
// równania różniczkowego
//-----------------------
double f(double x,double y)
{
  return x + y;
}

// Funkcja wynikowa
// obliczona analitycznie
//-----------------------
double rf(double x)
{
  return 2 * exp(x) - x - 1;
}

// Wypisuje wyniki
//----------------
void pisz(int c, double x, double y)
{
  cout << setw(3) << c
  << ": x = " << setw(7) << x
  << "  y = " << setw(8) << y
  << "| dokładnie y = "
  << setw(8) << rf(x)
  << endl;
}

// Program główny
//---------------

int main()
{
  SetConsoleOutputCP(CP_UTF8);
  SetConsoleCP(CP_UTF8);

  cout << setprecision(5)
       << fixed;

  // inicjujemy zmienne
  double x = x_s;
  double y = y_s;
  double k1,k2,k3,k4;
  int    c = 0;

  cout << "Rozwiązywanie równania "
          "różniczkowego metodą "
          "Rungego-Kutty\n"
          "-----------------------"
          "---------------------"
          "-------------\n\n";
  pisz(c,x,y);
  // Metoda Rungego-Kutta
  while(x < x_k)
  {
    c++;
    k1 = f(x,y);
    k2 = f(x+h/2,y+h/2*k1);
    k3 = f(x+h/2,y+h/2*k2);
    k4 = f(x+h,y+h*k3);
    y  = y+h/6*(k1+2*k2+2*k3+k4);
    x = x_s + c * h;
    pisz(c,x,y);
  }
  cout << endl;
  system("pause");
  return 0;
}

Rozwiązywanie równania różniczkowego metodą Rungego-Kutty
---------------------------------------------------------

  0: x = 0.00000  y =  1.00000| dokładnie y =  1.00000
  1: x = 0.12500  y =  1.14130| dokładnie y =  1.14130
  2: x = 0.25000  y =  1.31805| dokładnie y =  1.31805
  3: x = 0.37500  y =  1.53498| dokładnie y =  1.53498
  4: x = 0.50000  y =  1.79744| dokładnie y =  1.79744
  5: x = 0.62500  y =  2.11149| dokładnie y =  2.11149
  6: x = 0.75000  y =  2.48399| dokładnie y =  2.48400
  7: x = 0.87500  y =  2.92274| dokładnie y =  2.92275
  8: x = 1.00000  y =  3.43655| dokładnie y =  3.43656
  9: x = 1.12500  y =  4.03542| dokładnie y =  4.03543
 10: x = 1.25000  y =  4.73067| dokładnie y =  4.73069
 11: x = 1.37500  y =  5.53513| dokładnie y =  5.53515
 12: x = 1.50000  y =  6.46335| dokładnie y =  6.46338
 13: x = 1.62500  y =  7.53181| dokładnie y =  7.53184
 14: x = 1.75000  y =  8.75917| dokładnie y =  8.75921
 15: x = 1.87500  y = 10.16659| dokładnie y = 10.16664
 16: x = 2.00000  y = 11.77806| dokładnie y = 11.77811
 17: x = 2.12500  y = 13.62073| dokładnie y = 13.62079
 18: x = 2.25000  y = 15.72539| dokładnie y = 15.72547
 19: x = 2.37500  y = 18.12693| dokładnie y = 18.12703
 20: x = 2.50000  y = 20.86488| dokładnie y = 20.86499

# Równania różniczkowe
# metoda Rungego-Kutty
# Metody numeryczne 0086
# --------------------------
# (C)2026 mgr Jerzy Wałaszek

from math import exp

# Stałe
#------
# Przedział (start,koniec)
x_s = 0.0
x_k = 2.5
# Warunek początkowy
y_s = 1.0
# Odstęp na osi OX
h   = 0.125

# Funkcja prawej strony
# równania różniczkowego
#-----------------------
def f(x, y):
    return x + y

# Funkcja wynikowa
# obliczona analitycznie
#-----------------------
def rf(x):
    return 2 * exp(x) - x - 1

# Wypisuje wyniki
#----------------
def pisz(c, x, y):
    print(f"{c:3d}: x = {x:7.5f}    "
          f"y = {y:8.5f}| "
          f"dokładnie y = {rf(x):8.5f}")

# Program główny
#---------------

# inicjujemy zmienne
x = x_s
y = y_s
c = 0
print("Rozwiązywanie równania "
      "różniczkowego metodą "
      "Rungego-Kutty\n"
      "-----------------------"
      "---------------------"
      "-------------\n")
pisz(c, x, y)

# Metoda Rungego-Kutty
while x < x_k:
    c += 1
    k1 = f(x, y)
    k2 = f(x + h/2, y + h/2 * k1)
    k3 = f(x + h/2, y + h/2 * k2)
    k4 = f(x + h, y + h * k3)
    y  = y + h/6 * (k1 + 2*k2 + 2*k3 + k4)
    x = x_s + c * h
    pisz(c, x, y)

print()
input("Naciśnij Enter...")