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).

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.}$

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.}$

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 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.

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ń.

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.

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 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;
}

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()

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.

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

• metoda Eulera
• metoda Eulera-Maruyamy
• równanie sztywne - problem stabilności rozwiązań numerycznych

Inne:

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

• 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.

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