Metody RungegoKutty

Metody RungegoKutty (RK) ciesz¡ si¦ du»¡ popularno±ci¡. W odró»nieniu od metod wielokrokowych s¡ samostartuj¡ce, za± w porównaniu z integrato-rami Taylora cechuj¡ si¦ prostot¡ i uniwersalno±ci¡. Metody RK rz¦du wy»-szego ni» pierwszy otrzymujemy dzi¦ki u»yciu pomocniczych w¦zªów przy-bli»onych w przedziale od x0 do x0+ h. W pewnym sensie mo»na uzna¢, »e metody RK polegaj¡ na tworzeniu ±redniej wa»onej z wyników otrzymywa-nych metod¡ Eulera. Je±li wspomniaotrzymywa-nych pomocniczych w¦zªów (wliczaj¡c równie» dokªadny x0) jest s, to mówimy o s-etapowej metodzie RK.

Dowoln¡ s-etapow¡ metod¦ RK deniujemy podaj¡c jej tablic¦ Butchera ci ai,j

bj

, i, j = 1, 2, . . . , s. (4.23) Tablica Butchera zawiera: s² wspóªczynników ai,j, s w¦zªów ci, oraz s wag b_i. Elementy tablicy nie s¡ niezale»ne, gdy» obowi¡zuj¡ dwa ograniczenia

∑s j=1

a_i,j = c_i,

∑s j=1

b_j = 1. (4.24)

Je±li znamy tablic¦ Butchera, to równania deniuj¡ce s-etapow¡ metod¦ RK maj¡ posta¢

{ y₁ = y₀+ h^∑s_n=1 b_nk_n, kn = f

(

y₀+ h^∑s_j=1an,jkj, x0+ cnh )

. (4.25)

Nale»y wyra¹nie rozdzieli¢ dwie sytuacje powstaj¡ce przy wyliczaniu ko-lejnych wektorów pomocniczych k1, k2 itd. Albo wyliczenie kn wymaga je-dynie znajomo±ci wcze±niej wyliczonych kn−1, kn−2itd. albo równania (4.25) maj¡ charakter przest¦pny. W pierwszej sytuacji mamy do czynienia z otwar-tymi (jawnymi) metodami RK a w drugiej  z metodami zamkni¦otwar-tymi (uwi-kªanymi).

Metody otwarte RK s¡ najprostsze w u»yciu. Rozpoznajemy je po wy-gl¡dzie tablicy Butchera; poniewa» warunkiem aby metoda byªa otwarta jest

a_i,j = 0, dlaj­ i, (4.26)

to gªówna cz¦±¢ tablicy ma posta¢ trójk¡tn¡ doln¡

0 0 0 0 . . . 0 c2 a2,1 0 0 . . . 0 c₃ a_3,1 a_3,2 0 . . . 0 . . . . . . . . . . . . . . . . . .

cs as,1 as,2 as,3 . . . 0 b1 b2 b3 . . . bs

Zauwa»my, »e dla ka»dej metody otwartej RK wektor pomocniczy k1 = f (y₀, x0).

Je±li spróbujemy skonstruowa¢ jednoetapow¡ (s = 1) metod¦ otwart¡

RungegoKutty, to okazuje si¦, »e natychmiast otrzymujemy metod¦ Eu-lera (4.18) jako jedyne mo»liwe rozwi¡zanie. Istotnie: a1,1 = 0 z wªa±ciwo±ci (4.26), za± c1= 0 i b1 = 1 wynikaj¡ z (4.24). Mamy wtedy

y₁ = y₀+ h k₁,

k₁ = f (y₀, x₀), (4.27)

zgodnie z ogólnym wzorem (4.25).

Je±li chcemy skonstruowa¢ otwart¡ metod¦ dwuetapow¡, to w zagadnie-niu pojawi¡ si¦ dwa parametry a i b. W ±wietle poznanych wªasno±ci tablic Butchera, mo»emy napisa¢ dla s = 2

0 0 0

a a 0

b 1− b a wi¦c

y₁ = y₀+ h b k₁+ h (1− b) k2, k₁ = f (y₀, x₀),

k2 = f (y₀+ h a k1, x0+ a h).

Spróbujmy teraz uzgodni¢ powy»sze wzory z wielomianem Taylora (4.22).

W tym celu rozwiniemy k2 w szereg pot¦g kroku h k₂ = f (y₀, x₀) + a h

[∂f

∂x ]

x0

+ a h [Df ]_x

0 k₁+ O(h²). (4.28) Po wstawieniu tak rozwini¦tego k2 oraz k1= f₀do wyra»enia dla y₁, otrzy-mujemy

y₁ = y₀+ f₀h + (1− b) a [∂f

∂x + (Df )f ]

x0

h²+ O(h³).

Powy»szy wzór stanie si¦ równowa»ny wielomianowi ekstrapolacyjnemu Tay-lora drugiego stopnia (4.22), je»eli tylko przyjmiemy

(1− b) a = 1

2. (4.29)

Jak wida¢, istnieje niesko«czenie wiele mo»liwo±ci wyboru wspóªczynników a i b, które zdeniuj¡ nam dwuetapowo¡ metod¦ otwart¡ RK tak, aby byªa ona metod¡ drugiego rz¦du. W praktyce spotykamy jednak tylko dwa warianty zwane metod¡ trapezów (a = 1, b = 1/2) i metod¡ punktu ±rodkowego (a = 1/2, b = 0).

Metoda trapezów ma tablic¦ Butchera

0 0 0

a wi¦c przyjmuje posta¢ algorytmu k₁ = f (y₀, x₀),

k2 = f (y₀+ h k1, x0+ h).

y₁ = y₀+ (h/2) (k1+ k2) . (4.30) Natomiast metoda punktu ±rodkowego z tablic¡ Butchera

0 0 0

W podobny do przedstawionego powy»ej sposób otrzymuje si¦ metody RK o dowolnej liczbie etapów. Najbardziej znana jest czteroetapowa metoda RK4

Mo»na j¡ zapisa¢ w postaci k1 = f (y₀, x0),

k2 = f (y₀+ (h/2) k1, x0+ h/2).

k₃ = f (y₀+ (h/2) k₂, x₀+ h/2).

k₄ = f (y₀+ h k₃, x₀+ h).

y₁ = y₀+ (h/6) (k1+ 2 k2+ 2 k3+ k4) . (4.32) Jest to metoda czwartego rz¦du, z bª¦dem lokalnym proporcjonalnym do h⁵. Stosuje si¦ j¡ do±¢ powszechnie do generowania startowych w¦zªów metod wielokrokowych.

Zajmijmy si¦ teraz uwikªanymi (zamkni¦tymi) metodami RungegoKutty (IRK  skrót od implicit RungeKutta). Rozpatrzmy pozornie trywialny przykªad metody jednoetapowej. Dla metody otwartej z s = 1 nie mieli±my

»adnych parametrów swobodnych (por. wyprowadzenie (4.27)). Dla metody zamkni¦tej tablica Butchera dopuszcza jeden nieokre±lony parametr a

a a 1 i ogólny wzór (4.25) przyjmuje posta¢

k1 = f (y₀+ h a k1, x0+ h a),

y₁ = y₀+ hk1 = y₀+ h f (y₀+ h a k1, x0+ h a). (4.33) Rozwi«my teraz wzór (4.33) w szereg pot¦g h, aby sprawdzi¢, czy mo»na go uzgodni¢ z wielomianem Taylora (4.22).

y₁ = y₀+ h f (y₀, x0) + h²a [∂f

∂x ]

x0

+ h²a [Df ]_x0 k1+ O(h³).

Poniewa» k1− f0= O(h), mo»emy przepisa¢ powy»szy wzór jako

y₁ = y₀+ h f₀+ h²a ([∂f

∂x ]

x0

+ [Df ]_x

0 f₀ )

+ O(h³).

Wystarczy teraz przyj¡¢ a = ¹₂, aby nasze wyra»enie dla y₁ ró»niªo si¦ od wielomianu Taylora (4.22) tylko wyrazami z trzeci¡ pot¦g¡ h.

W ten sposób zdeniowali±my jednoetapow¡ metod¦ IRK drugiego rz¦du.

Ma ona tablic¦ Butchera

1/2 1/2 1

i opisana jest wzorami

k1 = f (y₀+ (h/2) k1, x0+ h/2),

y₁ = y₀+ h k1. (4.34)

Metoda ta nazywana jest uwikªan¡ metod¡ punktu ±rodkowego (ang. implicit midpoint) i mo»e te» zosta¢ przeksztaªcona do postaci

y₁= y₀+ h f ((y₀+ y₁)/2, x₀+ h/2) . (4.35) Przedstawiona tu metoda zachowuje si¦ lepiej nie tylko od metody Eulera lecz równie» od obydwu dwuetapowych metod otwartych RK, mimo i» for-malnie ma ten sam rz¡d co te ostatnie. Za dobre wªasno±ci tej metody pªa-cimy konieczno±ci¡ rozwi¡zywania ukªadu równa« przest¦pnych (4.34) lub (4.35) na ka»dym kroku caªkowania. W praktyce mo»na stosowa¢ metod¦

iteracji prostych: po przyj¦ciu wst¦pnego przybli»enia k^[0]₁ = f (y₀, x0),

iterujemy równanie

k^[n+1]₁ = f (y₀+ (h/2) k^[n]₁ , x0+ h/2), n = 0, 1, . . .

a» do momentu, gdy ró»nica |k^[n+1]₁ −k^[n]₁ | jest mniejsza od przyj¦tej dokªad-no±ci oblicze«.

Jaka jest górna granica wydajno±ci uwikªanych metod s-etapowych ? Okazuje si¦, »e metoda IRK mo»e mo»e osi¡gn¡¢ rz¡d 2s. Metody IRK rz¦du 2s nazywamy metodami optymalnymi RungegoKutty lub metodami GaussaLegendre'a. Ta druga nazwa wynika st¡d, »e w¦zªy ci metody opty-malnej musz¡ pokrywa¢ si¦ z pierwiastkami wielomianu Legendre'a stopnia s

P_s(2 c_i− 1) = 0.

Istotnie, dla s = 1 mamy P1(x) = x i w¦zeª c1 uwikªanej metody punktu

±rodkowego jest pierwiastkiem równania 2 c1 − 1 = 0. Metoda (4.34) jest najprostsz¡ z metod GaussaLegendre'a. Dwuetapowa (s = 2) metoda z tej rodziny ma rz¡d 4 i tablic¦ Butchera

1

2 −^√₆^{3 1}₄ ¹₄− ^√₆³

1

2 +^√₆^{3 1}₄ +^√₆^{3 1}₄

1 2

1 2

Metody optymalne pozwalaj¡ na znaczn¡ popraw¦ dokªadno±ci caªko-wania w porównaniu z klasycznymi metodami otwartymi. Na przykªad dla s = 10metoda otwarta jest rz¦du 8, natomiast metoda optymalna jest rz¦du 20. Nie mo»emy jednak ignorowa¢ faktu, »e tam, gdzie w metodzie otwartej wektor pomocniczy kn obliczany jest tylko raz na jeden krok integratora, w metodzie optymalnej musimy conajmniej kilka razy iterowa¢ równania przest¦pne. Ust¦pstwem na rzecz zmniejszenia kosztu obliczeniowego s¡ tak zwane metody GaussaRadau. S¡ to równie» metody uwikªane Rungego

Kutty, ale wyprowadzane przy zaªo»eniu warunku c₁= a_1,1 = . . . = a_1,s= 0,

tak, aby pierwszy wektor pomocniczy k1 = f (y₀, x₀) nie wymagaª iteracji.

Takie wymaganie oznacza, »e nie mo»na osi¡gn¡¢ rz¦du równie wysokiego jak w metodach optymalnych, lecz strata jest niewielka, gdy» s-etapowe metody GaussaRadau zapewniaj¡ rz¡d m = 2s − 1. Na przykªad o±mioetapowy integrator RA15 opracowany przez E. Everharta jest metod¡ pi¦tnastego rz¦du.