Zajmiemy się równaniem
y(x
0)=y
0 (1)
Podobnie jak przy konstrukcji metody Eulera postawmy sobie następujące zadanie: znaleźć wartość funkcji y=y(x), w x=x
1=x
0+h
Załóżmy, że rozwiązanie zadania (1) jest parabolą i skorzystajmy z następującej własności paraboli:
- współczynnik kierunkowy siecznej (A0A1) jest średnią
arytmetyczną współczynnika kierunkowego stycznej (A0S0)
i stycznej (A1S1).
Współczynnik kierunkowy stycznej (A
0S
0)
umiemy wyznaczyć: jest to f(x
0,y
0). Wyznaczmy przybliżoną
wartość funkcji y, w x=x
1 tak, jak to robiliśmy w metodzie Eulera,
czyli
Teraz możemy wyznaczyć przybliżoną wartość współczynnika kierunkowego stycznej (A
1S
1), wynosi on
Następnie wyznaczmy współczynnik kierunkowej siecznej (A
0A
1),
jest on równy
Jest to przybliżona wartość tego współczynnika ze względu na sposób liczenia y
1E. Zamiast siecznej (A
0A
1)
mamy sieczną (A
0A
1'). Z trójkąta A
0PA
1'
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
2).
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
1,w
2,w
3,... takich, by wartość funkcji y(x)
dla x=x
1, czyli y
1 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
0+ah,y
0+
ak
1)
jej rozwinięciem w szereg Taylora, czyli
(5)
gdzie
Z (5) po uporządkowaniu i wstawieniu k
1=hf mamy
(6)
Załóżmy, że y
0*=y(x
0)
i y
1*=y(x
1) 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
2 z obu rozwinięć. Otrzymujemy
w
1+w
2=1 ,
w
2a=1/2 ,
w
2a=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=1/2, w
2=1/2, są to wzory (4).
Wartość funkcji y w x=x
2=x
1+h=x
0+2h możemy
policzyć stosując (4), o ile (x
1, y(x
1)=y
1)
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
|
k 1=hf(x i,y i)
k2=hf(xi+h,yi+k1)
|
Heuna
|
O(h3)
|
O(h2)
|
3
|
|
k 1=hf(x i,y i)
k 2=hf(x i+h/2,y i+k 1/2)
k3=hf(xi+h,yi-k1+2k2)
|
pokrewna metodzie
Simsona
|
O(h4)
|
O(h3)
|
4
|
|
k 1=hf(x i,y i)
k 2=hf(x i+h/2,y i+k 1/2)
k 3=hf(x i+h/2,y i+k 2/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
0+(i+1)h).
Dla i=0 liczymy wartości y
1=y(x
1)=y(x
0+h)=y(0.1).
Najpierw policzymy
dalej
Analogicznie obliczenia przeprowadzamy dla i=1 (czyli dla x
2=x
0+2h),
policzymy y
2=y(x
2). Policzymy wartości rozwiązania ścisłego
dla x=x
1 i x=x
2 i wyniki tych obliczeń zestawimy w tabeli
i+1
|
x
|
y
|
k=0.1(x+y)
|
Y(x)=2ex-x-1
|
0
|
0
|
1
|
|
|
|
|
|
|
|
1
|
0.1
|
1.11034
|
|
1.1103420
|
|
|
|
k 1=0.121034
k 2=0.13220857
k 3=0.132638285
k4=0.1442978285
|
|
2
|
0.2
|
1.242803
|
|
1.2428055
|