1.5.4 Metody typu Runge-Kutty

Zajmiemy się równaniem

               y(x₀)=y₀                                                     (1)

Podobnie jak przy konstrukcji metody Eulera postawmy sobie następujące zadanie: znaleźć wartość funkcji y=y(x), w x=x₁=x₀+h

Załóżmy, że rozwiązanie zadania (1) jest parabolą i skorzystajmy z następującej własności paraboli:

- współczynnik kierunkowy siecznej (A₀A₁) jest średnią arytmetyczną współczynnika kierunkowego stycznej (A₀S₀) i stycznej (A₁S₁).

Współczynnik kierunkowy stycznej (A₀S₀) umiemy wyznaczyć: jest to f(x₀,y₀). Wyznaczmy przybliżoną wartość funkcji y, w x=x₁ tak, jak to robiliśmy w metodzie Eulera, czyli



Teraz możemy wyznaczyć przybliżoną wartość współczynnika kierunkowego stycznej (A₁S₁), wynosi on



Następnie wyznaczmy współczynnik kierunkowej siecznej (A₀A₁), jest on równy



Jest to przybliżona wartość tego współczynnika ze względu na sposób liczenia y₁^E. Zamiast siecznej (A₀A₁) mamy sieczną (A₀A₁^'). Z trójkąta A₀PA₁^' policzymy




Uporządkujemy dotychczasowe rozwiązania wprowadzając oznaczenia

                                                       (2)

Wzory (2) są wzorami Runge-Kutty II rzędu. Można wykazać, że całkowity błąd metody jest rzędu O(h²).

Spróbujmy dojść do tych samych wzorów (2) inną drogą i wyprowadzić wzory Runge-Kutty wyższych rzędów. Zapiszmy uogólnienie algebraiczne wzorów (2)

                                                                   (3)

Będziemy szukać stałych a,b,c,... , a,b,g,... i w₁,w₂,w₃,... takich, by wartość funkcji y(x) dla x=x₁, czyli y₁ była możliwie bliska wartości dokładnej. Uzyskamy to rozwijając w szereg Taylora funkcję przybliżoną i ścisłą, żądając zgodności (do wyrazu odpowiedniego rzedu) tych rozwinięć. Wykonajmy obliczenie dla metody Runge-Kutty II rzędu, czyli wartość przybliżoną obliczamy następująco

                                                                               (4)

stąd otrzymujemy



Zastąpmy f(x₀+ah,y₀+ak₁) jej rozwinięciem w szereg Taylora, czyli

                                                     (5)

gdzie



Z (5) po uporządkowaniu i wstawieniu k₁=hf mamy

                                                 (6)

Załóżmy, że y₀^*=y(x₀) i y₁^*=y(x₁) są dokładnymi wartościami funkcji y i zapiszmy korzystając z rozwinięcia w szereg Taylora

                                             (7)

Żądamy zgodności obu (6, 7) rozwinięć do drugiego wyrazu rozwinięcia w szereg Taylora. Porównajmy więc współczynniki przy h i h² z obu rozwinięć. Otrzymujemy

w₁+w₂=1 ,

w₂a=1/2 , w₂a=1/2 ,

stąd



Mamy więc całą rodzinę wzorów Runge-Kutty II rzędu zależną od parametru a. I tak dla a=1 mamy a=1, w₁=1/2, w₂=1/2, są to wzory (4).

Wartość funkcji y w x=x₂=x₁+h=x₀+2h możemy policzyć stosując (4), o ile (x₁, y(x₁)=y₁) potraktujmy tak, jak warunek początkowy. Czyli ogólnie (4) zapiszemy



Wzory wyższych rzędów można otrzymać analogicznie. W poniższej tabeli są zapisane najczęściej używane współczynniki dla metod typu Runge-Kutty rzędu od I do IV

Rząd metody P	Stałe (wzór 3) w	Wartości współczynników k	Nazwa metody	Błąd lokalny e	Błąd całkowity E
1	w1=1	k1=hf(xi,yi)	Eulera	O(h2)	O(h)
2	w1=w2=1/2	k1=hf(xi,yi) k2=hf(xi+h,yi+k1)	Heuna	O(h3)	O(h2)
3	w1=w3=1/6 w2=2/3	k1=hf(xi,yi) k2=hf(xi+h/2,yi+k1/2) k3=hf(xi+h,yi-k1+2k2)	pokrewna metodzie Simsona	O(h4)	O(h3)
4	w1=w4=1/6 w2=w3=1/3	k1=hf(xi,yi) k2=hf(xi+h/2,yi+k1/2) k3=hf(xi+h/2,yi+k2/2) k4=hf(xi+h,yi+k3)	klasyczna metoda Runge-Kutty	O(h5)	O(h4)

Metody typu Runge-Kutty są łatwe do zaprogramowani, zmiana kroku całkowania może być dokonywana w dowolnym etapie obliczeń i nie wymaga dużego nakładu pracy, są metodami samostartującymi, tzn. znajomość warunku początkowego wystarcza, by rozpocząć obliczenia. Wymagają one jednak wielokrotnego (dla metody rzędu p p-krotnego) obliczania wartości funkcji f w każdym kroku całkowania, mogą więc być metodami kosztownymi (zwykle najbardziej czasochłonnymi, a więc i kosztownym zadaniem jest obliczanie wartości funkcji f).

Przykład

Znaleźć rozwiązanie zagadnienia początkowego

y'(x)=x+y

y(0)=1

metodą Runge-Kutty IV rzędu. Obliczenia wykonać dla xÎ[0,0.2] z krokiem h=0.1 . Porównać otrzymane wyniki z rozwiązaniem ścisłym y(x)=2e^x-x-1.

Wzory Runge-Kutty IV rzędu (patrz (3) i powyższa tabela) mają postać



gdzie



Dla przypomnienia y_i+1=y(x_i+1)=y(x₀+(i+1)h). Dla i=0 liczymy wartości y₁=y(x₁)=y(x₀+h)=y(0.1). Najpierw policzymy



dalej



Analogicznie obliczenia przeprowadzamy dla i=1 (czyli dla x₂=x₀+2h), policzymy y₂=y(x₂). Policzymy wartości rozwiązania ścisłego dla x=x₁ i x=x₂ i wyniki tych obliczeń zestawimy w tabeli

i+1	x	y	k=0.1(x+y)	Y(x)=2ex-x-1
0	0	1
	0. 0.05 0.05 0.1	1. 1.05 1.055 1.1105	k1=0.1 k2=0.11 k3=0.1105 k4=0.12105
1	0.1	1.11034		1.1103420
	0.1 0.15 0.15 0.2	1.11034 1.170857 1.17638285 1.242978285	k1=0.121034 k2=0.13220857 k3=0.132638285 k4=0.1442978285
2	0.2	1.242803		1.2428055