4. Metody dla równań różniczkowych zwyczajnych - rząd schematów
4.3. Schematy jednokrokowe
W tym podrozdziale wprowadzimy pojęcie schematu jednokrokowego:
Definicja 4.6. Dla zadania początkowego (3.1) schematem jednokrokowym dla stałego kroku h = T −t0N nazywamy równanie różnicowe:
xn+1= xn+ h φ(h, tn, xn, xn+1) n = 0, . . . , N (4.10) gdzie tj = t0+j h a φ jest funkcją ciągłą określoną na [0, H)×[t0, T )×Ux0×Ux0 dla Ux0 otoczenia x0. Dodatkowo, jeśli φ nie zależy od xn+1, to schemat jednokrokowy nazywamy otwartym, a w przeciwnym wypadku mówimy o schemacie zamkniętym.
W przypadku schematów otwartych możemy wyliczyć xn+1 znając xn, natomiast w przy-padku schematów zamkniętych musimy rozwiązać liniowy, bądź nieliniowy układ równań. Do tej pory poznaliśmy dwa schematy jednokrokowe (które zarazem są schematami liniowymi wie-lokrokowymi) - czyli oba schematy Eulera i schemat trapezów.
Analogicznie do przypadku schematów liniowych wielokrokowych, zgodnie z Definicją 4.4, schemat jednokrokowy ma rząd p 1, jeśli dla x ∈ Cp+1([t0, T )) rozwiązania zagadnienia (3.1) dla 1 ¬ p, h = T −t0N i t ∈ [t0, T ]) lokalny błąd schematu spełnia:
eh(t) := |x(t + h) − x(t) − h φ(h, t,dx dt(t),
dx
dt(t + h))| ¬ C h p+1,
ze stałą C niezależną od t, czyli eh = O(hp+1) dla dostatecznie gładkiego rozwiązania. 4.3.1. Schematy Rungego-Kutty
Podstawową klasą schematów jednokrokowych są tzw. schematy Rungego-Kutty lub - mó-wiąc krótko - schematy Rungego. Idea ich jest prosta.
Załóżmy, że znamy xn, i chcemy wyliczyć wartość xn+1 ze wzoru uwzględniającego wartość pola wektorowego nie tylko w xn, ale również w dodatkowym punkcie ˜x. Wtedy
xn+1 = F (h, tn, xn, ˜x).
Biorąc schemat otwarty Eulera z krokiem ˜h otrzymujemy punkt ˜
4.3. Schematy jednokrokowe 41 który, jak wiemy, przybliża x(t + ˜h), ale niedokładnie. Możemy policzyć wartość pola wekto-rowego f w tym punkcie i następnie, wykorzystując wartość yn = f (tn, xn) i ˜y = f (tn+ ˜h, ˜x), znaleźć lepsze przybliżenie x(tn+1) - czyli np. za przybliżenie pola wektorowego wziąć ważoną średnią obu wartości b yn+ c˜y dla pewnych ustalonych wag b, c. Możliwości jest wiele. Pojawia się pytanie: jak oceniać różne konstrukcje F ? Można tak dobierać F , aby rząd schematu był możliwie duży.
Załóżmy, że ˜h = a h. Wtedy szukamy schematu postaci:
xn+1= xn+ b h fn+ c h f (tn+ a h, xn+ a h f (tn, xn)) (4.11) tak, aby schemat miał maksymalny rząd.
Rozwijamy rozwiązanie x w szereg Taylora: x(t + h) − x(t) = hdx
dt(t) + 0.5 h
2 d2x
dt2(t) + O(h3) i rozwijając ostatni z członów (4.11) w punkcie (t, x) otrzymujemy:
f (t + a h, x + a h f (t, x)) = f + a h fxf + a h ft = f + a h(fxf + ft) = dx dt + a h d2x dt2.
Skorzystaliśmy z tego, że ddt2x2 = dtdf (t, x(t)) = fxf + ft. Zatem, wstawiając dwa ostatnie rów-nania do (4.11) otrzymujemy warunki na to, aby schemat był rzędu dwa:
b + c = 1 c a = 0.5
Rysunek 4.5. Schemat Heuna w porównaniu ze schematem otwartym Eulera na [0,1] z h = 0.1.
Rysunek 4.6. Zmodyfikowany schemat Euler w porównaniu ze schematem otwartym Eulera na [0,1] z h = 0.1.
1. Zmodyfikowany schemat Eulera
xn+1= xn+ h f tn+h 2, xn+h 2 fn (4.12) dla c = 1, b = 0, a = 0.5, 2. Schemat Heuna xn+1 = xn+h 2 (fn+ f (tn+1, xn+ h fn)) , (4.13) dla b = c = 0.5; a = 1.
Warto zauważyć, że w niektórych publikacjach wszystkie schematy otwarte Rungego-Kutty rzędu dwa nazywane są zmodyfikowanym schematem Eulera.
Na rysunkach4.6i4.5pokazano rozwiązania uzyskane tymi dwoma schematami dla zadania dx/dt = x z x(0) = 1 na [0, 1]. Widać, że wykresy się pokrywają z wykresem rozwiązania. W porównaniu do otwartego schematu Eulera widzimy znaczącą poprawę. Zobaczmy, co się dzieje na dłuższym odcinku czasu w przypadku schematu Heuna, por. rysunek4.7.
Graficzne wytłumaczenie zmodyfikowanego schematu Eulera
Na rysunku4.8 zawarliśmy graficzne wytłumaczenie jednego kroku zmodyfikowanego sche-matu Eulera.
Punkt A - oznacza (t, x), czyli wcześniej obliczone przybliżenie dla czasu t. Chcemy wyzna-czyć przybliżenie dla czasu t + h. Otwarty schemat Eulera w jednym kroku przyjmuje za różnicę między kolejnymi punktami h f (t, x), czyli idzie w kierunku pola wektorowego w punkcie t, czyli na naszym rysunku daje to punkt C. Z kolei w zmodyfikowanym schemacie Eulera przyjmujemy za różnicę h pomnożone przez kierunek pola wyznaczonego w dodatkowym pomocniczym punk-cie ˆx = x + 0.5 fn oznaczonym jako B, tzn. idziemy w kierunku f (t + 0.5 h, ˆx) i otrzymujemy w efekcie punkt D.
Na rysunku widać, że jeśli nachylenie pola wektorowego mocno się zmienia, to pole w punkcie B powinno mieć lepszy kierunek niż w A, czy w E. Jest to oczywiście argument heurystyczny.
4.3. Schematy jednokrokowe 43
Rysunek 4.7. Schemat Heuna na [0,100] i exp(t). Wykres dla t bliskich 100.
Rysunek 4.8. Graficzne wytłumaczenie zmodyfikowanego schematu Eulera.
Analogicznie konstruuje się schematy Rungego wyższych rzędów poprzez wprowadzenie więk-szej ilości kroków pośrednich, jak również schematy zamknięte Rungego - dopuszczając wartość fn+1 w schemacie.
Rysunek 4.9. Schematy Rungego rzędu cztery i schemat Heuna rzędu dwa na [0, 1].
Rysunek 4.10. Schematy Rungego rzędu cztery i schemat Heuna rzędu dwa na [0, 50] w skali pół-logarytmicznej.
Podamy kilka wzorów na powszechnie używany otwarty schemat Rungego-Kutty czwartego rzędu.
4.4. Zadania 45 Najpierw definiujemy cztery wartości:
K1 = f (tn, xn) K2 = f (tn+h2, xn+h2K1) K3 = f (tn+h2, xn+h2K2) K4 = f (tn+ h, xn+ h K3) (4.14) i otrzymujemy ostateczny wzór: xn+1= xn+h 6(K1+ 2 K2+ 2 K3+ K4) (4.15) schematu rzędu cztery (co oczywiście wymaga dowodu).
Istnieje oczywiście cała rodzina otwartych schematów Rungego-Kutty rzędu cztery. Tu po-daliśmy tylko przykładowy schemat z tej rodziny.
Popatrzmy jak działa ten schemat w porównaniu ze schematem Heuna dla naszego modelo-wego zagadnienia dx/dt = x z x(0) = 1.
Na odcinku [0, 1] dla h = 0.1 oba schematy praktycznie pokrywają się z rozwiązaniem, por. rysunek4.9. Na rysunku4.10widać wykresy rozwiązań na [0, 50] (w skali pół-logarytmicznej),
Rysunek 4.11. Schematy Rungego rzędu cztery i schemat Heuna rzędu dwa na [44, 50].
a na rysunku4.11wykres na odcinku [44, 50] - tu już widać, że schemat rzędu cztery jest jednak znacznie lepszy.
W rozdziale5.4przedstawimy metodę eksperymentalną badania rzędu schematu. Przy okazji zobaczymy ogromną różnicę w dokładności obu schematów, por. Tabela 5.2.
4.4. Zadania
Ćwiczenie 4.1. Udowodnij, że wzory (4.1) i (4.2) są prawdziwe dla funkcji dostatecznie gład-kich.
Ćwiczenie 4.2. Pokaż, że rząd schematów Eulera wynosi jeden, o ile rozwiązanie zadania Cauchy’ego jest klasy C2.
Ćwiczenie 4.3. Pokaż, że rząd schematu kroku środkowego wynosi k, o ile rozwiązanie zadania Cauchy’ego jest klasy Ck+1 dla k = 1, 2.
Ćwiczenie 4.4. Znajdź rząd schematu Taylora (4.4) dla rozwiązania dostatecznie gładkiego. Ćwiczenie 4.5. Znajdź wzór na schemat Taylora rzędu trzy.
Ćwiczenie 4.6. Rozpatrzmy rodzinę schematów:
xn+1 = xn+ h (a fn+ b fn+1)
Określ rząd schematu w zależności od wartości parametrów a, b. Dla jakich a, b rząd jest naj-większy? Dla jakich wartości parametrów schemat będzie zamknięty, a dla jakich otwarty? Ćwiczenie 4.7. Wyprowadź otwarty dwukrokowy schemat Adamsa bazujący na wielomianie interpolacyjnym stopnia jeden (żeby policzyć xn+2 potrzebujemy h, tn, xn+1, fn i fn+1, tak jak opisano w rozdziale 4.1.2). Zbadaj rząd schematu.
Rozwiązanie. Zgodnie z zasadą konstrukcji schematów Adamsa musimy scałkować na odcinku [tn, tn+ h] wielomian liniowy interpolacyjny P1(s) = fn+−h1 (fn− fn−1)(s − tn). Otrzymujemy schemat xn+1= xn+ 0.5 h (3 fn− fn−1).
Ćwiczenie 4.8. Wyprowadź otwarty trzykrokowy schemat Adamsa bazujący na wielomianie interpolacyjnym stopnia dwa (żeby policzyć xn+3 potrzebujemy h, tn, xn+2 i fn, fn+1, fn+2, tak jak opisano w rozdziale4.1.2). Zbadaj rząd schematu.
Ćwiczenie 4.9. Wyprowadź zamknięty dwukrokowy schemat Adamsa bazujący na wielomianie interpolacyjnym stopnia dwa (żeby policzyć xn+2 potrzebujemy h, tn, xn+1 i fn, fn+1, fn+2, tak jak opisano w rozdziale4.1.2). Zbadaj rząd schematu.
Ćwiczenie 4.10 (średnio trudne). Rozpatrzmy rodzinę schematów: xn+2= xn+1+ h (a fn+ b fn+1+ c fn+2).
Określ rząd schematu w zależności od wartości parametrów a, b, c. Dla jakich wartości rząd jest największy? Dla jakich wartości parametrów schemat będzie zamknięty, a dla jakich otwarty? Ćwiczenie 4.11 (trudne). Udowodnij, że schemat (4.15) ma rząd cztery.
Ćwiczenie 4.12 (laboratoryjne). Zbadaj eksperymentalnie metodą połowienia kroków rząd lokalnego błędu schematu dla schematów:
— otwartego schematu Eulera (3.4), — zamkniętego schematu Eulera (3.5), — schematu midpoint (4.3),
— schematu Taylora (4.4), — schematu Heuna (4.13), — schematu trapezów (4.7),
— zmodyfikowanego schematu Eulera (4.12), — schematu Rungego rzędu cztery (4.15).
zastosowanych do modelowego zadania dxdt = 1 + x2 z x(0) = 0 z rozwiązaniem x(t) = tan(t). Tzn. dla ustalonego t, np. t = 1 i kolejnych połowionych kroków hk = hk−1
2 = 2−kh0 z h0 = 1e − 1 liczymy lokalny błąd schematu ehk(t), por. (4.9), i następnie stosunek eehk(t)
hk+1(t). Jeśliby ten stosunek wynosił w przybliżeniu 2p+1, to lokalny błąd schematu zachowuje się jak O(hp+1), oznacza to, że schemat posiada rząd p przynajmniej dla tego zadania początkowego.