Algorytm Rungego-Kutty

Algorytm Rungego-Kutty (metoda Rungego-Kutty) – metoda numeryczna do iteracyjnego rozwiązywania równań różniczkowych zwyczajnych (opisujących procesy fizyczne, chemiczne, biologiczne itp.) Opracowana około 1900 przez niemieckich matematyków: Carla Rungego oraz Martina Kuttę.

Obecnie nazwą metody Rungego-Kutty określa się rodzinę jawnych i niejawnych metod wielokrokowych, jak również pewne ich modyfikacje. Istnieje wiele metod RK, o wielu stopniach, wielu krokach, różnych rzędach, i różniących się między sobą własnościami (jak stabilność, jawność, niejawność, metody osadzone, szybkość działania itp.).

Potocznie metodą Rungego-Kutty określa się najczęściej metodę Rungego-Kutty 4. rzędu ze współczynnikami podanymi poniżej, gdyż jest powszechnie stosowana ze względu na prostotę implementacji, relatywnie proste wzory, dużą szybkość oraz wysoki rząd metody.

Na początku artykułu zostanie omówione podejście do rozwiązywania równań różniczkowych rzędu 1-go, a następnie zostanie pokazane uogólnienie na równania rzędu 2-go i wyższych. Całość dopełnia przykład programu napisany w języku C++, rozwiązujący równanie różniczkowe 2-go rzędu nieliniowe (opis drgań wahadła).

Równanie różniczkowe zwyczajne rzędu 1[edytuj | edytuj kod]

Niech dane będzie równanie różniczkowe zwyczajne 1-go rzędu postaci

${\displaystyle y'=f(x,y)}$

oraz warunek początkowy

${\displaystyle x_{0},\,y_{0}\equiv y(x_{0}),}$

gdzie:

• ${\displaystyle x}$ – zmienna niezależna,
• ${\displaystyle y=y(x)}$ – szukana funkcja zmiennej ${\displaystyle x,}$
• ${\displaystyle y'}$ – pochodna zmiennej zależnej ${\displaystyle y}$ względem zmiennej niezależnej ${\displaystyle x,}$
• ${\displaystyle f(x,y)}$ – znane wyrażenie algebraiczne, zależne od zmiennych ${\displaystyle x}$ oraz ${\displaystyle y.}$

Metoda numeryczna znalezienia rozwiązania polega na obliczeniu dyskretnych wartości ${\displaystyle y_{1},y_{2},\dots ,y_{n}}$ szukanej funkcji w kolejnych krokach obliczeń. Niech ${\displaystyle h}$ oznacza przyjęty do obliczeń, niewielki skok zmiennej niezależnej ${\displaystyle x,}$ tj. zmienna ${\displaystyle x}$ przyjmuje w kolejnych krokach dyskretne wartości ${\displaystyle x_{n},n=0,1,\dots }$ zwiększające się o wartość ${\displaystyle h.}$

Metoda RK 4. rzędu[edytuj | edytuj kod]

Iterację zaczynamy zadając ${\displaystyle x_{0},\,y_{0}}$ oraz funkcję ${\displaystyle f(x,y).}$ W ${\displaystyle n}$-tym kroku iteracji, dla ${\displaystyle n=0,1,2,\dots ,}$ obliczamy następujące wielkości:

(1) ${\displaystyle x_{n+1}=x_{n}+h,}$

(2) ${\displaystyle k_{1}=h\cdot f\left(x_{n},y_{n}\right),}$

(3) ${\displaystyle k_{2}=h\cdot f{\Big (}x_{n}+{\frac {h}{2}},y_{n}+{\frac {k_{1}}{2}}{\Big )},}$

(4) ${\displaystyle k_{3}=h\cdot f{\Big (}x_{n}+{\frac {h}{2}},y_{n}+{\frac {k_{2}}{2}}{\Big )},}$

(5) ${\displaystyle k_{4}=h\cdot f\left(x_{n}+h,y_{n}+k_{3}\right),}$

następnie liczymy

(6) ${\displaystyle dy_{n}={\frac {1}{6}}(k_{1}+2k_{2}+2k_{3}+k_{4})}$

i ostatecznie otrzymamy:

(7) ${\displaystyle y_{n+1}=y_{n}+dy_{n}.}$

Widać, że ${\displaystyle y_{n+1}}$ zależy od wartości wcześniej obliczonej ${\displaystyle y_{n}}$ oraz wielkości kroku ${\displaystyle h.}$

W ten sposób otrzymuje się (podobnie jak w innych iteracyjnych metodach rozwiązywania równań różniczkowych) kolejne wartości ${\displaystyle y_{n},n=1,2,\dots ,}$ które przybliżają szukaną funkcję ${\displaystyle y.}$

Metoda RK 2. rzędu[edytuj | edytuj kod]

Istnieje cała klasa metod RK 2. rzędu, z czego z nazwy wyróżnia się przynajmniej dwie metody.

Jawna metoda RK 2, znana jako metoda punktu pośredniego (ang. midpoint method)

${\displaystyle y_{n+1}=y_{n}+h\cdot f\left(x_{n}+{\frac {h}{2}},y_{n}+{\frac {1}{2}}k_{1}\right).}$

Oznaczenia takie same, jak wyżej.

Metoda RK 1. rzędu, czyli metoda Eulera[edytuj | edytuj kod]

Metoda Eulera jest szczególnym przypadkiem metod Rungego-Kutty (metoda Eulera pojawiła się historycznie najpierw, w związku z tym zachowano tradycyjną nazwę):

${\displaystyle y_{n+1}=y_{n}+k_{1}.}$

Oznaczenia takie same jak poprzednio.

Równanie różniczkowe zwyczajne rzędu 2[edytuj | edytuj kod]

Jeżeli mamy równanie różniczkowe zwyczajne rzędu 2-go, to zamienia się je na układ dwóch równań rzędu 1-go i stosuje wyżej wymienione metody Rungego-Kutty, iteracyjnie rozwiązując dany układ równań.

Przykład: Równanie drgań wahadła[edytuj | edytuj kod]

1. Zamiana równania 2-go rzędu na układ 2 równań 1-go rzędu[edytuj | edytuj kod]

Równanie drgań wahadła matematycznego ma postać^[1]:

${\displaystyle \theta (t)''+k\cdot \sin \theta (t)=0,}$

gdzie:

${\displaystyle \theta (t)}$ – kąt odchylenia wahadła od pionu w chwili ${\displaystyle t,}$

${\displaystyle k=g/\ell }$ – stała liczba (${\displaystyle g}$ – przyspieszenie ziemskie, ${\displaystyle \ell }$ – długość nici wahadła).

Podstawiając ${\displaystyle \omega =\theta (t)',}$ równanie to sprowadza się do układu dwóch równań 1-go rzędu

${\displaystyle {\begin{cases}\theta (t)'=\omega ,\\\omega '+k\cdot \sin \theta =0,\end{cases}}}$

który można rozwiązać jedną z metod Rungego-Kutty.

2. Obliczenia numeryczne – obliczenie okresu drgań wahadła metodą Eulera[edytuj | edytuj kod]

Poniżej podano kod programu w C++ liczącego okres drgań wahadła – zastosowano całkowanie równania ruchu metodą Eulera. Pomimo prostoty metody uzyskuje się bardzo duże dokładności obliczeń, gdy dobierze się odpowiednio mały krok iteracji ${\displaystyle dt.}$ W programie podanym niżej przejęto ${\displaystyle dt=0.000001.}$ (Dokładność metody można sprawdzić, porównując wyniki z obliczeniami okresu drgań wahadła, opartymi o całki eliptyczne 1-go rodzaju – por. Tabela okresów drgań wahadła)

Program można testować online, korzystając np. z darmowych compilatorów języka C++.

// Metoda Eulera całkowania równania różniczkowego drgań wahadła matematycznego #include <iostream> #include <cmath> using namespace std;  const double g = 9.81, l = 0.25;//9.81;  // Funkcja liczy połowę okresu drgań wahadła matematycznego double PolOkresDrgan(double theta0, double dt) {     double theta, theta_2, v, v_2, t, PolOkres;     // Inicjalizacja zmiennych     theta = theta0; //Kąt odchylenia od pionu     v = 0.0;  //prędkość kątowa     t = 0.0;     // Obliczania Theta_2, v_2 dla kolejnych kroków czasowych     while (true) {         theta_2 = theta + dt * v;         v_2 = v - dt * (g / l) * sin(theta);         if(v * v_2 < 0) break; //wahadło osiąga makymalne wychylenie, gdy v zmiania znak - upływa wtedy pół okresu drgań         // Aktualizacja wartości         theta = theta_2;         v = v_2;         t = t + dt;     }     PolOkres = t;     return PolOkres; }  int main() {     double angle, theta0, dt, Okres;     angle= 179;   //Amplituda drgań w stopniach     dt = 0.000001; // krok iteracji      theta0 = angle * M_PI / 180.0;    // Konwersja kąta na radiany     Okres = 2.0 * PolOkresDrgan(theta0, dt);     cout << " theta0 = " << angle << " [stopnie]      T = " << Okres <<" [sekundy]"<< endl;      return 0; } 

Przykład: Obliczenia numeryczne równania rzędu 2-go metodą Rungego-Kutty rzędu 4.[edytuj | edytuj kod]

Napisanie programu, rozwiązującego równanie różniczkowe rzędu 2-go metodą RK rzędu 4. wymaga liczenia podwójnej liczby współczynników. Poniżej podano kod programu w języku python; program rozwiązuje równanie różniczkowe wahadła z tłumieniem. Współczynniki metody RK są zdefiniowane w liniach 14 – 30 programu.

Znalezienie rozwiązania analitycznego równania ruchu złożonych układów fizycznych jest w ogólnym przypadku niemożliwe. Jednak metody numeryczne pozwalają efektywnie rozwiązywać te równania ruchu. Np. równanie ruchu wahadła z tłumieniem i z siłą wymuszającą ma postać

${\displaystyle {\frac {d^{2}\theta (t)}{dt^{2}}}+2\beta {\frac {d\theta (t)}{dt}}+{\frac {g}{\ell }}\sin \theta (t)=F(t),}$

gdzie:

${\displaystyle \beta }$ – współczynnik tłumienia,

${\displaystyle F(t)}$ – siła wymuszająca, zależna dowolnie od czasu.

Wprowadzając nową zmienną ${\displaystyle \omega =\theta (t)'}$ powyższe równanie sprowadza się do układu dwóch równań różniczkowych pierwszego rzędu

${\displaystyle {\begin{cases}{\frac {d\theta (t)}{dt}}=\omega (t)\\{\frac {d\omega (t)}{dt}}=-2\beta \omega (t)-{\frac {g}{\ell }}\sin \theta (t)+F(t)\end{cases}}}$

Następnie układ tych równań rozwiązuje się iteracyjnie, np. metodą Rungego-Kutty, co prowadzi do znalezienia dyskretnych wartości ${\displaystyle \theta (t_{n}),\omega (t_{n})}$ w zadanym przedziale całkowania równań.

Pokazany tu przykładowy program zakłada zerową siłę wymuszającą, tj. ${\displaystyle F(t)=0.}$ Uwzględnienie niezerowej siły wymuszającej w kodzie programu jest nieskomplikowane – wystarczy uzupełnić linię 11 kodu, uzupełniając funkcję o dodatkowy parametr – aktualną wartość siły wymuszającej.

Dodatkowo program generuje cztery wykresy zależności ${\displaystyle T(\theta _{0},c),}$ gdzie ${\displaystyle c}$ – współczynnik tłumienia, korzystając online z biblioteki mathplotlib.pyplot.

Program można testować, korzystając np. z darmowego notatnika colab google online.

import math import matplotlib.pyplot as plt  g = 9.81  # Przyspieszenie ziemskie L = 9.81   # Długość wahadła  T_0 = 2 * math.pi * math.sqrt(L / g)  def f(theta, omega, c):     # Równania ruchu dla wahadła z tłumieniem     omega_dot = (-c * omega - (g / L) * math.sin(theta))     return omega_dot  def runge_kutta(theta, omega, h, c):     k1_theta = h * omega     k1_omega = h * f(theta, omega, c)      k2_theta = h * (omega + k1_omega / 2)     k2_omega = h * f(theta + k1_theta / 2, omega + k1_omega / 2, c)      k3_theta = h * (omega + k2_omega / 2)     k3_omega = h * f(theta + k2_theta / 2, omega + k2_omega / 2, c)      k4_theta = h * (omega + k3_omega)     k4_omega = h * f(theta + k3_theta, omega + k3_omega, c)      theta += (k1_theta + 2 * k2_theta + 2 * k3_theta + k4_theta) / 6     omega += (k1_omega + 2 * k2_omega + 2 * k3_omega + k4_omega) / 6      return theta, omega  def main():     # DANE:     theta_0 = 170  # Początkowy kąt wychylenia [w stopniach]      h = 0.00001          # Krok czasu     iterations = 4000000  # Liczba iteracji      c_values = [0, 0.2, 0.5, 1.5]  # Wartości współczynnika tłumienia c     print(c_values)      for c in c_values:         omega = 0         theta = theta_0 / 180 * math.pi         time_values = []         theta_values = []          print("Tłumienie c =", c)         time_2 = 0         for i in range(iterations):             time = i * h             time_values.append(time / T_0)             theta_values.append(theta)              omega_2 = omega             theta, omega = runge_kutta(theta, omega, h, c)             if (omega * omega_2 < 0) or (omega_2 == 0):                 print("T/2/T_0=", (time - time_2) / T_0)                 time_2 = time          plt.plot(time_values, theta_values)      # WYKRESY - RYSOWANIE:     plt.xlabel(r'Czas unormowany $t/T_0$')     plt.ylabel(r'$\theta(t)$ [stopnie]')     plt.title('Wykres kąta wychylenia wahadła z tłumieniem od czasu')     plt.grid(True)     plt.legend(['c = 0', 'c = 0.2', 'c = 0.5', 'c = 1.5'])      plt.show()  if __name__ == "__main__":     main() 

Równanie różniczkowe zwyczajne rzędu n[edytuj | edytuj kod]

Jeżeli mamy równanie różniczkowe zwyczajne rzędu ${\displaystyle n,}$ to w szczególnych wypadkach można obniżyć jego rząd. Zagadnienie obniżania stopnia równań różniczkowych omawia np. ^[2] Np. gdy dane jest równanie różniczkowe w postaci jawnej

${\displaystyle F\left(x,y,y',y'',\ \dots ,\ y^{(n-1)}\right)=y^{(n)},}$

to zamienia się je na układ ${\displaystyle n}$ równań rzędu 1 poprzez wprowadzenie dodatkowych ${\displaystyle n-1}$ zmiennych ${\displaystyle u_{i}=y^{(i-1)},\quad i=2,3,\dots ,n,}$

tj. pochodne wyższych rzędów niż 1 traktuje się jako nowe zmienne; przyjmując dodatkowo oznaczenie ${\displaystyle u_{1}\equiv y}$ dla jednolitości zapisu otrzymamy układ n równań różniczkowych pierwszego rzędu w postaci:

${\displaystyle {\begin{array}{rcl}u_{1}'&=&u_{2}\\u_{2}'&=&u_{3}\\&\vdots &\\u_{n-1}'&=&u_{n}\\u_{n}'&=&F(x,u_{1},\dots ,u_{n}).\end{array}}}$

lub bardziej zwięźle w notacji wektorowej:

${\displaystyle \mathbf {u} '=\mathbf {F} (x,\mathbf {u} )}$

gdzie:

${\displaystyle \mathbf {u} =(u_{1},\dots ,u_{n}),\quad \mathbf {F} (x,u_{1},\dots ,u_{n})=[u_{2},u_{3},\dots ,u_{n},F(x,u_{1},\dots ,u_{n})]}$

Otrzymany układ równań rozwiązuje się stosując jedną z wyżej wymienionych metod Rungego-Kutty.

Zobacz też[edytuj | edytuj kod]

Całkowanie równań różniczkowych:

• metoda Eulera
• metoda Eulera-Maruyamy

Inne:

• całkowanie numeryczne
• metody numeryczne
• metoda Monte Carlo
• metoda Monte-Carlo Metropolisa-Hastingsa

Przypisy[edytuj | edytuj kod]

Bibliografia[edytuj | edytuj kod]

• Z. Fortuna, B. Macukow, J. Wąsowski, Metody numeryczne Podręczniki akademickie Elektronika, informatyka, telekomunikacja, Warszawa: Wydawnictwa Naukowo-Techniczne, 1982, s. 285–312.
• R.S. Guter, A.R. Janpolski: Równania różniczkowe. Warszawa: Wydawnictwo Naukowe PWN, 1980.
• D. Potter, Metody obliczeniowe fizyki, fizyka komputerowa, Warszawa: PWN, 1982, s. 19–43.
• J. Szmelter, Metody komputerowe w mechanice, Warszawa: PWN, 1980, s. 150–157.

Linki zewnętrzne[edytuj | edytuj kod]

• MyPhysicsLab – Runge-Kutta Algorithm – omówienie (w języku angielskim)
• Numerical Recipes in C – Runge-Kutta Method – algorytm w języku C