Algorytm Rungego-Kutty
Rodzaj | Metoda numeryczna do iteracyjnego rozwiązywania równań różniczkowych zwyczajnych |
---|
![](http://upload.wikimedia.org/wikipedia/commons/thumb/0/00/Runge-kutta.svg/400px-Runge-kutta.svg.png)
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
Niech dane będzie równanie różniczkowe zwyczajne 1-go rzędu postaci
oraz warunek początkowy
gdzie:
- – zmienna niezależna,
- – szukana funkcja zmiennej
- – pochodna zmiennej zależnej względem zmiennej niezależnej
- – znane wyrażenie algebraiczne, zależne od zmiennych oraz
Metoda numeryczna znalezienia rozwiązania polega na obliczeniu dyskretnych wartości szukanej funkcji w kolejnych krokach obliczeń. Niech oznacza przyjęty do obliczeń, niewielki skok zmiennej niezależnej tj. zmienna przyjmuje w kolejnych krokach dyskretne wartości zwiększające się o wartość
Metoda RK 4. rzędu
Iterację zaczynamy zadając oraz funkcję W -tym kroku iteracji, dla obliczamy następujące wielkości:
- (1)
- (2)
- (3)
- (4)
- (5)
następnie liczymy
- (6)
i ostatecznie otrzymamy:
- (7)
Widać, że zależy od wartości wcześniej obliczonej oraz wielkości kroku
W ten sposób otrzymuje się (podobnie jak w innych iteracyjnych metodach rozwiązywania równań różniczkowych) kolejne wartości które przybliżają szukaną funkcję
Metoda RK 2. rzędu
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)
Oznaczenia takie same, jak wyżej.
Metoda RK 1. rzędu, czyli metoda Eulera
Metoda Eulera jest szczególnym przypadkiem metod Rungego-Kutty (metoda Eulera pojawiła się historycznie najpierw, w związku z tym zachowano tradycyjną nazwę):
Oznaczenia takie same jak poprzednio.
Równanie różniczkowe zwyczajne rzędu 2
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
1. Zamiana równania 2-go rzędu na układ 2 równań 1-go rzędu
Równanie drgań wahadła matematycznego ma postać[1]:
gdzie:
- – kąt odchylenia wahadła od pionu w chwili
- – stała liczba ( – przyspieszenie ziemskie, – długość nici wahadła).
Podstawiając równanie to sprowadza się do układu dwóch równań 1-go rzędu
który można rozwiązać jedną z metod Rungego-Kutty.
2. Obliczenia numeryczne – obliczenie okresu drgań wahadła metodą Eulera
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 W programie podanym niżej przejęto (Dokładność metody można sprawdzić, porównując wyniki z obliczeniami okresu drgań wahadła, opartymi na całkach eliptycznych pierwszego 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.
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ć
gdzie:
- – współczynnik tłumienia,
- – siła wymuszająca, zależna dowolnie od czasu.
Wprowadzając nową zmienną powyższe równanie sprowadza się do układu dwóch równań różniczkowych pierwszego rzędu
Następnie układ tych równań rozwiązuje się iteracyjnie, np. metodą Rungego-Kutty, co prowadzi do znalezienia dyskretnych wartości w zadanym przedziale całkowania równań.
Pokazany tu przykładowy program zakłada zerową siłę wymuszającą, tj. 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 gdzie – 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
Jeżeli mamy równanie różniczkowe zwyczajne rzędu 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
to zamienia się je na układ równań rzędu 1 poprzez wprowadzenie dodatkowych zmiennych
tj. pochodne wyższych rzędów niż 1 traktuje się jako nowe zmienne; przyjmując dodatkowo oznaczenie dla jednolitości zapisu otrzymamy układ n równań różniczkowych pierwszego rzędu w postaci:
lub bardziej zwięźle w notacji wektorowej:
gdzie:
Otrzymany układ równań rozwiązuje się stosując jedną z wyżej wymienionych metod Rungego-Kutty.
Zobacz też
Całkowanie równań różniczkowych:
Inne:
Przypisy
- ↑ Resnick i Halliday 1999 ↓, s. 361–364.
- ↑ Guter i Janpolski 1980 ↓, s. 134–139.
Bibliografia
- 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
- MyPhysicsLab – Runge-Kutta Algorithm – omówienie (w języku angielskim)
- Numerical Recipes in C – Runge-Kutta Method – algorytm w języku C
- p
- d
- e
- SNL: Runge-Kutta-metoder