Opiszemy tu jedynie niektóre metody różnicowe, to jest takie, które rów-nanie różniczkwe zastępują pewnym rówrów-naniem różnicowym. Weźmy pod uwagę zagadnienie Cauchy’ego
(6.22) d
dtu(t) = f (t, u(t)),
(6.23) u(t0) = u0.
Załóżmy, że rozwiązanie u istnieje i jest jednoznaczne w przedziale [t0, t0+ a] i że ma w tym przedziale tyle pochodnych ile będzie nam potrzeba. Na tym przedziale zbudujemy siatkę punktów, dla uproszczenia, równoodległych
t0 < t1 < t2· · · < tN, tj = t0+ jh, h = a
N, h > 0.
Liczbę h będziemy nazywali krokiem siatki lub krokiem całkowania. W prak-tyce często potrzebne są siatki ze zmiennym krokiem, my jednak ograni-czymy się tu do siatek o stałym kroku. Rozwiązanie u zagadnienia Cau-chy’ego (6.22),(6.23) rozwiniemy przy pomocy wzoru Taylora dla t = tk+1, w punkcie tk u(tk+1) = u(tk) + hu0(tk) + h 2 2 u 00(tk) + h 3 6u 000(tk) + · · · ,
lub
u(tk+1) = u(tk) + hf (tku(tk)) + O(h2).
Odrzucając wyrazy zawierające h w potędze 2 i wyższej, otrzymamy
równa-nie różnicowe
(6.24) uk+1 = uk+ hf (tk, uk), u0 = u0.
Rozwiązaniem tego równania jest ciąg {uj}, j = 0, 1, · · ·. Element ciągu uk
odpowiada wartości rozwiązania u(tk), zagadnienia Cauchy,ego (6.22),(6.23). Można postąpić inaczej; rozwinąć u(tk) w punkcie tk+1
u(tk) = u(tk+1) + hu0(tk+1) + h 2 2 u
00(tk+1) + O(h3).
Podobnie jak poprzednio, odrzucając wyrazy zawierające h w potędze 2 i wyższej otrzymamy inne równanie różnicowe
(6.25) uk+1 = uk + hf (tk+1, uk+1), u0 = u0.
Wygodnie będzie dalej oznaczać
fk = f (tk, uk).
Równania (6.24) i (6.25) noszą nazwę Schematów Eulera (schemat = me-toda).
(6.24) uk+1 = uk+ hfk Schemat otwarty Eulera,
(6.25) uk+1 = uk + hfk+1 Schemat zamknięty Eulera.
Schematem otwartym Eulera łatwo się posługiwać. Znając warunek począt-kowy u0, drogą podstawiania do wzoru (6.24) obliczymy uk dla każdego in-teresującego nas k.
Zupełnie inaczej jest ze schematem zamkniętym (6.25). Aby wyliczyć uk+1 znając uk, trzeba rozwiązać równanie nieliniowe (układ m-równań!)
uk+1 = uk+ hf (tk+1, uk+1).
Chwilowo nie potrafimy nic powiedzieć na temat zależności ciągów o elementach uk, oraz u(tk), k = 0, 1, 2, · · ·. Chcielibyśmy, aby spełniony był
Warunek zbieżności schematu. Przypuśćmy, że dla dowolnego ustalonego t ∈ [t0, t0+ a], siatka punktów {tk} została tak dobrana, że t = tk = t0+ kh.
Będziemy uważać rozważany schemat za zbieżny, jeśli warunek
(6.26) uk → u(t), gdy h → 0 (stąd k = t − t0
h , k → ∞).
zachodzi
• dla dowolnego rozwiązania u dowolnego równania (6.1) z warunkiem początkowym u(t0) = u0, należącego do klasy równań spełniających za-łożenia Twierdzenia Picard’a - Lindel¨of ’a,
• dla dowolnego rozwiązania {uk}, k = 0, 1, · · · rozważanego schematu, dla którego wartość u0 = u0(h) spełnia warunek
u0 → u0, gdy h → 0.
Dalej będą nas interesowały jedynie te schematy, które są zbieżne w powyż-szym sensie. Zobaczymy też jak odróżniać schematy zbieżne od niezbieżnych.
Chwilowo powróćmy do schematu zamkniętego Eulera
uk+1 = uk + fk+1.
Aby obliczyć uk+1 trzeba rozwiązać zadanie na punkt stały
x = Φ(x),
gdzie x = uk+1, Φ(x) = uk+ hf (tk+1, x). Sprobujmy zastosować Twierdzenie
Banacha o punkcie stałym. Zbudujmy ciąg wektorów {xl} l = 0, 1, 2, · · ·, x0 -dowolny element,
xl+1 = Φ(xl).
Ciąg będzie zbiegał do punktu stałego x = uk+1, jeśli Φ spełmia warunek
Lipschitza ze stałą L1, 0 ≤ L1 < 1. Przypuśćmy, że funkcja f (prawa
strona równania (6.1)) spełnia założenia Twierdzenie Picard’a - Lindel¨of’a ze stałą Lipschitza L. Wtedy dla dowolnych x i y takich, że (t, x) i (t, y) należą do dziedziny funkcji f
Widzimy, że L1 < 1, gdy
(6.27) 0 < h < 1
L.
Zatem iteracja
(6.28) ulk+1+1 = uk + hf (tk+1, ulk+1) zbiega do uk+1 dla dowolnego punktu startowego u0
k+1, jeśli
h < 1 L.
Warunek (6.27) nie jest bardzo ograniczający, jeśli stała Lipschitza L funk-cji f nie jest zbyt duża. W przypadku wielkich wartości L lepiej stosować
iterację Newtona. Zauważmy, że koszt algorytmu wykorzystującego schemat
zamknięty skupia się głównie w wyliczaniu wartości funkcji f . Zatem należy wyliczać wartości f jak najmniej razy. Liczba iteracji zależy od tego jak
dobrze dobrany został punkt startowy u0
k+1. Dobry start iteracji za-pewnia przyjęcie jako u0
k+1 wartości uk+1 uzyskanej z zastosowania schematu otwartego Eulera.
W ten sposób doszliśmy do tak zwanej METODY PREDICTOR -
COR-RECTOR opartej na schematach Eulera.
• PREDICTOR, to schemat otwarty podający punkt startowy dla
ite-racji - stosowany 1 raz na krok.
• CORRECTOR, to schemat zamknięty służący do iterowania.
Iteru-jemy małą liczbę razy, gdyż punkt startowy jest blizko rozwiązania. Metodę PREDICTOR - CORRECTOR w taki sam sposób można bu-dować w oparciu o inne pary schematów 16, złożone ze schematu otwartego (PREDICTOR) i zamknkniętego (CORRECTOR).
Narzuca sie pytanie: po co stosować skomplikowane w użyciu
sche-maty zamknięte, skoro dysponujemy bardzo wygodnymi schema-tami otwartymi? Okazuje się, że pewne cechy stawiają metodę zamkniętą
zdecydowanie wyżej od metody otwartej. Są zadania, których nie daje się 16Schematy takie poznamy w dalszej części tego wykładu.
wogóle policzyć metodą otwartą, a którym daje radę metoda zamknięta. To co odróżnia schemat Eulera otwarty od zamkniętego, to na pewno nie jest
rząd.
Co to jest rząd schematu?
Niech u(t) będzie rozwiązaniem zagadnienia Cauchy’ego (6.1), (6.2), o którym zakładamy, że ma p + 1 pochodnych ciągłych. Oznaczmy przez
(6.29) S({ul}, l = 0, 1, 2, · · ·) = 0 nasz schemat różnicowy.
Mówimy, że schemat (6.29) jest rzędu p, jeśli podstawiając do (6.29) ciąg {u(tj)}, j = 0, 1, 2, · · ·
zamiast ciągu {uj} j = 0, 1, 2, · · · , otrzymamy S({u(tj)}, j = 0, 1, 2, · · ·) = R,
gdzie reszta R spełnia warunek
R = O(hp+1),
zaś istnieje takie zadanie Cauchy’ego spełniające powyższe warunki, dla któ-rego R 6= O(hp+2).
Biorąc pod uwagę sposób w jaki otrzymaliśmy oba schematy Eulera wi-dzimy, że oba są rzędu 1.
Zanim przejdziemy, do wyjaśnienia na czym polega wyższość schematu zamkniętego nad otwartym, przyjrzyjmy się jeszcze innym schematom różni-cowym. Niech u ∈ C3. Mamy
(6.29) u(tk+1) = u(tk) + hu0(tk) + h 2 2!u 00(tk) + h 3 3!u 000(tk) + · · · , (6.30) u(tk+1) = u(tk) + hu0(tk+1) −h 2 2!u 00(tk+1) + h 3 3!u 000(tk+1) + · · · . Zauważmy jeszcze, że u00(tk+1) = u00(tk)+hu000(tk)+O(h2). Dodajmy stronami wzory (6.29) i (6.30) uwzględniając powyższą uwagę. Otrzymamy tak zwany
schemat trapezów
(6.31) uk+1 = uk +h
Jest to schemat zamknięty, rzędu 2.
Zadanie. Zapisz iteracje Banacha i Newtona dla schematu trapezów.
Wszystkie trzy schematy, które dotychczas poznaliśmy są jednokrokowe, to znaczy, że mając do dyspozycji jedynie uk, możemy wyliczyć uk+1.
Zadanie. Schematy Taylora. Używając rozwinięcia Taylora dla
rozwiąza-nia u(t) zagadnierozwiąza-nia początkowego (6.1), (6.2), uwzględrozwiąza-niając drugie i ewen-tualnie wyższe pochodne u zbuduj schematy jednokrokowe rzędu
wyż-szaego niż 1.
Wskazówka. Zauważ, że u00(t) = ∂
∂tf (t, u(t)) + ∂
∂uf (t, u(t))f (t, u(t)).
Podobnie dla wyższych pochodnych.
Odwołajmy się jeszcze raz do wzoru Taylora. Podobnie jak poprzednio
u(tk+1) = u(tk) + hu0(tk) + h 2 2 u(tk) + h3 6 u 000(tk) + dots, u(tk−1) = u(tk) − hu0(tk) + h 2 2 u(tk) − h 3 6 u 000(tk) + · · · . Odejmijmy stronami te równości. Otrzymamy schemat ”Midpoint”
uk+2 = uk+ 2hfk+1.
Schemat Midpoint, to schemat otwarty. Nie jest on schematem jednokroko-wym, gdyż uk+2 możemy wyliczyć tylko jeśli dysponujemy dwoma warto-ściami uk i uk+1. Aby schemat mógł wystartować potrzebne są dwa warunki
początkowe uo i u1. Mówimy, że taki schemat nie jest samostartujący. Jeśli dysponujemy warunkiem początkowym u0, to aby uruchomić schemat Midpo-int musimy dodatkowo doliczyć wartość u1. Można to zrobić używając jakiejś metody jednokrokowej. Nie jest jednak obojętne jakiej metody użyjemy. Ze sposobu konstrukcji schematu Midpoint wynika, że jest on rzędu 2 (reszta od-rzucona jest rzędu O(h3)). Zatem dla zachowania rzędu powinniśmy zadbać o to, aby u1 wyliczyć również schematem rzędu 2.
Okazuje się, że schemat Midpoint, mimo że ma rząd 2, zawodzi w pewnych przypadkach z którymi schemat otwarty Eulera (który jest rzędu 1) radzi sobie całkiem dobrze.
Zadanie. Napisz program rozwiązujący zagadnienie Cauchy’ego d
dtu(t) = −λu(t), λ > 0. u(0) = 1.
Użyj schematu otwartego Eulera i schematu Midpoint dla tego samego zada-nia. Porównaj zachowanie się schematów gdy wykonujesz dużą liczbę kroków przy jednakowej wartości kroku h i stałej λ > 0. Porównaj co się dzieje dla różnych wartości h i λ.
Schematy liniowe wielokrokowe. Przykładem takiego schematu jest
sche-mat Midpoint. Schesche-mat liniowy q - krokowy jest równaniem różnicowym, na
ogól nieliniowym, postaci
(6.30) q X j=0 αjuk+j = h q X j=0 βjfk+j,
gdzie jak poprzednio fl = f (tl, ul).
Aby wystartować, taki schemat potrzebuje q warunków początkowych
u0, u1, · · · , uq−1, które trzeba doliczyć schematem jednokrokowym
odpowied-nio wysokiego rzędu. Współczynniki αj, βj, j = 0, 1, · · · , q można
wyzna-czyć tak, aby rząd schematu był odpowiednio wysoki, oraz żeby posiadał on jeszcze inne cechy, o których powiemy póżniej. Zauważmy teraz, że schemat (6.30) jest
• otwarty, gdy βq= 0,
• zamknięty, gdy βq 6= 0.
Zadanie. Zbuduj schemat postaci (6.30) dlla q=1, który ma najwyższy
Powróćmy jeszcze do schematów jednokrokowych. Specjalną klasę takich schematów stanowią schematy Runge - Kutty. Schemat Runge - Kutty q - poziomowy jest postaci
(6.31) uk+1 = uk+ h[c1K1 + c2K2+ · · · + cqKq], gdzie (6.32) Kj = f (tk + haj, uk + h q X l=1 bj,lKl) j = 1, 2, · · · , q. Współczynniki (6.32) c1 c2 c3 · · · cq a1 a2 a3 · · · aq b1,1 b1,2 b1,3 · · · b1,q · · · · · · · · · · bq,1 bq,2 bq,3 · · · bq,q
wyznacza się tak, aby uzyskać możliwie wysoki rząd, oraz jeszcze inne cechy schematu. Taką cechą może być na przykład jego otwartość. Schemat będzie otwarty, jeśli zażądamy, aby
bj,l= 0 dla l ≥ j.
Schemat zamknięty wymaga rozwiązania na każdym kroku układu qm rów-nań dla wyznaczenia K1, K2, · · · , Kq. Współczynniki (6.32) dla różnych sche-matów są znane od wielu dziesiątek lat.
Przytoczymy tu dwa przykłady schematów Runge - Kutty.
Schemat 4- poziomowy otwarty, rzędu 4. uk+1 = uk +h 6[K1+ 2K2+ 2K3 + K4], K1 = f (tk, uk), (6.33) K2 = f (tk +h 2, uk + h 2K1), K3 = f (tk +h 2, uk + h 2K2),
K4 = f (tk + h, uk + hK3), Jest to bardzo często używany schemat.
Schemat 2- poziomowy zamknięty, rzędu 4. uk+1 = uk+ h 2[K1+ K2], (6.34) K1 = f (tk+ (1 2+ √ 3 6 )h, uk + h 4K1+ ( 1 2 + √ 3 6 )hK2), K2 = f (tk + (1 2 − √ 3 6 )h, uk+ ( 1 2− √ 3 6 )K1+ h 4K2).
Udowodniono, że można zbudować schematy otwarte Runge - Kutty, dla których liczba poziomów oraz rząd spełniają następujące zależności
liczba poziomów rząd 1 1 2 2 3 3 4 4 5 4 6 5 7 6 8 6 9 7 q ≤ 10 q − 2
Z tej tabelki widać, że schematy otwarte 4 poziomowe są optymalne w tym sensie, że osiągają maksymalny rząd przy minimalnej liczbie poziomów. Dla schematów zamkniętych q - poziomowych, można zawsze osiągnąć rząd 2q.
Koszt schematu determinowany jest liczbą obliczeń wartości prawej strony równania f na każdym kroku całkowania. Zatem widzimy, że schematy Runge
Zadanie. Wyprowadź wzory dla dwupoziomowego schematu Runge - Kutty,
otwartego.
uk+1 = uk+ h[c1K1 + c2K2]. Ile takich schematów rzędu 2 można zbudować?
Dotychczas, mówiąc o schematach różnicowych, podawaliśmy jako istotną ich cechę rząd. Pamiętamy jednak, że najważniejszą cechą schematu jest jego
zbieżność. Jakie znaczenie dla funkcjonowania schematu ma jego rząd
wyja-śnia teoria zbieżności schematów różnicowych. Podstawowe fakty z tej teorii, dla przypadku schematów jednokrokowych przytoczymy poniżej.
Teoria zbieżności schematów jednokrokowych.
Nierówność Gronwall’a. Niech ciąg liczb nieujemnych {vk} k = 0, 1, · · · , spełnia nierówność0 ≤ vk+1 ≤ Avk + B, k = 0, 1, · · ·
gdzie A, B > 0, to wtedy dla każdego k = 0, 1, · · ·
(6.35) 0 ≤ vk ≤ Akv0+
(Ak−1
A−1 B gdy A 6= 1 kB gdy A = 1.
Zadanie. Udowodnij nierówność Gronwall’a. Wskazówka: zastosuj indukcję
względem k.
Teraz będziemy rozważać schematy jednokrokowe otwarte 17 postaci (6.35) uk+1 = uk+ hΦ(h, tk, uk), h > 0
Zapis ten obejmuje wszystkie rozważane przez nas schematy jednokrokowe otwarte.
17Schemat zamknięty, jesli jest stosowalny, musi dać się rozwikłać przynajmiej lokalnie. Otrzymamy wtedy jego lokalny odpowiednik otwarty.
Konsystentność. Mówimy, że schemat (6.35) jest konsystentny, jeśli • funkcja Φ jest ciągła (względem wszystkich swoich argumentów) w całej
swojej dziedzinie,
• Φ spełnia warunek Lipschitza względem zmiennej u:
istnieje stała L, taka że dla wszystkich (h, t, u1), (h, t, u2) z dziedziny Φ
|Φ(h, t, u1) − Φ(h, t, u2)| ≤ L|u1− u2|, gdzie |.| oznacza normę w Rm,
• φ(0, t, u) = f(t, u), gdzie rozpatrywane przez nas równanie ma postać d
dtu(t) = f (t, u(t)).
Rozpatrujemy zagadnienie Cauchy’ego
d
dtu(t) = f (t, u(t)), u(t0) = u0,
oraz schemat jednokrokowy dla tego zagadnienia:
uk+1 = uk+ hΦ(h, tk, uk), u0= u0.
Twierdzenie o zbieżności z rzędem schematu jednokrokowego. Jeśli rozwiązanie u zagadnienia Cauchy’ego jest klasy Cp+1, p > 0 w przedziale
[t0, t0 + α] α > 0, w którym jest określone, i schemat jest konsystentny
oraz rzędu p, to schemat jest zbieżny i ponadto dla każdego ustalonego t = tk ∈ [t0, t0+ α]
|u(tk) − uk| ≤ Khp
, gdy h → 0, (h = tk − tk 0, k → ∞) gdzie K jest stałą niezależną od h.
Dowód. Podstawiając rozwiązanie zagadnienia Cauchy’ego u do schematu
różnicowego otrzymamy
u(tk+1) = u(tk) + hΦ(h, tk, u(tk)) + rk, uk+1 = uk + hΦ(h, tk, uk). Odejmując, otrzymamy
ek+1 = u(tk+1) − uk+1 = ek+ h[Φ(h, tk, u(tk)) − Φ(h, tk, uk)] + rk.
Ze względu na rząd schematu
|rk| ≤ Khp+1.
Ze względu na warunek Lipschitza otrzymamy:
|ek+1| ≤ (1 + hL)|ek| + Khp+1.
Zastosujmy teraz Nierówność Gronwalla dla A = 1 + hL i B = Khp+1. Otrzymamy |ek| ≤ (1 + hL)k|e0| + ((1+hL)k −1 hL Khp+1 dla L 6= 0, kKhp+1 dla L = 0. Ale 1 + hL ≤ ehL, i stąd (1 + hL)k ≤ ekhL ≤ eαL oraz kKhp+1 = khKhp ≤ αKhp. Ponadto przyjęliśmy, że u0 = u0, więc e0 = 0. Ostatecznie
|ek| ≤ ( eαLL Khp, gdy L > 0, αKhp gdy L = 0 = O(h p). 2
Zadanie. Udowodnij, że z samego założenia konsystentności wynika już
zbieżność schematu. Jednak nie otrzymujemy oszacowania błędu ek.
Z udowodnionego twierdzenia widać, jaką rolę odgrywa rząd schematu: jeśli
rozwiązanie u, które aproksymujemy jest dostatecznie gładkie (u ∈ Cp+1), oraz jeśli rząd schematu jest równy p, to |ek| ≤ Khp, gdy h → 0.