W tym rozdziale przedstawimy teorię zbieżności schematów jednokrokowych i wielokrokowych.
Rozpatrzymy tylko przypadek skalarny, ale teoria dla układów równań jest analogiczna. Wystarczy tylko moduł zastąpić przez jakąś normę w np. normę euklidesową.
Teoria zbieżności schematów jednokrokowych jest teorią odrębną od teorii dla schematów wielokrokowych liniowych.
Okaże się, że kluczowym pojęciem jest tu zgodność schematu jednokrokowego - inaczej konsystentność, którą definiujemy następująco:
Schemat jednokrokowy (4.10) jest zgodny (konsystentny), (ang. consistent) jeśli:
jest ciągłą ze względu na wszystkie zmienne
dla wszystkich
.
jest lipschitzowska ze względu na zmienne
i
tzn. istnieje
takie, że dla wszystkich
![]() |
Jeśli rozwiązanie zagadnienia początkowego (3.1) , schemat jednokrokowy jest zgodny i jest rzędu
, to
ten schemat jest zbieżny z rzędem
.
Dowód zostanie tutaj przedstawiony
dla prostoty tylko dla schematów otwartych z rzędem tj.
.
Dowód w całej ogólności można znaleźć np. w [23].
Oznaczmy przez , czyli błąd pomiędzy obliczonym schematem przybliżeniem
rozwiązania dla czasu
, a dokładną wartością rozwiązania
. Niech
,
czyli
to lokalny błąd schematu dla czasu
.
Wtedy otrzymujemy, że
![]() |
a stąd, korzystając ze zgodności schematu, a dokładniej lipschitzowskości funkcji ,
por. Definicja 5.1, otrzymujemy
![]() |
Dalej, poprzez indukcję matematyczną otrzymujemy:
![]() |
Korzystając z tego, że ()
![]() |
dla takich, że
widzimy, że
![]() |
Zauważmy, że . Widzimy też, że
![]() |
Ponieważ schemat ma rząd i
, to
zatem
![]() |
W szczególności zbieżność z rzędem oznacza dla ustalonego
i
, że
![]() |
Zgodność otwartego schematu Eulera.
W zasadzie sprawdzenie konsystentności (zgodności) schematu jest w tym przypadku oczywiste,
ponieważ funkcja , czyli spełnia założenia zgodności.
Dodatkowo wiemy, że schemat ma rząd jeden, więc jeśli
rozwiązanie zadania początkowego (3.1) należy do
, to
dla
, co potwierdzają wyniki eksperymentów z tabeli 4.1 w przypadku skalarnym dla równania
z
.
Teoria zbieżności dla schematów wielokrokowych liniowych różni się od teorii zbieżności dla schematów jednokrokowych. Precyzyjniej pisząc - używamy tak samo jak w przypadku schematów jednokrokowych pojęcie rzędu schematu, ale równie ważne jest pojęcie stabilności schematu. Używając nieformalnego języka - stabilność schematu oznacza to, że błędy z poprzednich kroków się nie kumulują, a wręcz zanikają.
Na początku rozważmy dowolny schemat liniowy wielokrokowy (4.8) zastosowany do równania z polem wektorowym równym zero tzn. do zagadnienia początkowego:
![]() |
którego jedynym rozwiązaniem jest rozwiązania stałe .
Nasz schemat staje się wtedy
równaniem różnicowym liniowym jednorodnym o stałych współczynnikach:
![]() |
(5.1) |
Wprowadzając wielomian
![]() |
(5.2) |
otrzymujemy, że jeśli jest zerem (pierwiastkiem) wielomianu
, to
ciąg
![]() |
jest rozwiązaniem równania różnicowego. Jeśli zero jest dwukrotne , to otrzymujemy nowe rozwiązanie związane z tym pierwiastkiem tzn.:
![]() |
i tak dalej - jeśli zero jest trzykrotne, to kolejne rozwiązanie:
![]() |
W każdym razie - co ważne - jeśli jakiś pierwiastek ma moduł
lub
, ale krotność pierwiastka jest większa od jeden, to istnieją rozwiązania nieograniczone, a dokładnie, zachodzi dla pewnych rozwiązań:
![]() |
Jeśli pierwiastek to zawsze wszystkie rozwiązania z nim związane spełniają
![]() |
Podsumowując: jeśli jakiś pierwiastek ma moduł mniejszy od jeden, albo ma moduł jeden i krotność jeden, to rozwiązania z nim związane są ograniczone. Wydaje się, że wymaganie żeby wszystkie rozwiązania schematu zastosowanego do równania z prawą stroną równą zero były ograniczone jest jak najbardziej uzasadnione.
Dlatego w ten sposób definiujemy stabilność (ang. stability) schematu liniowego wielokrokowego (4.8):
Schemat (4.8) jest stabilny, jeśli każdy pierwiastek wielomianu
spełnia
![]() |
a w przypadku jeśli , to krotność pierwiastka
wynosi jeden.
Dodatkowo wprowadzamy pojęcie silnej stabilności schematu (ang. strong stability):
Schemat (4.8) jest silnie stabilny, jeśli spełnia Definicję 5.2 stabilności, i
jeśli jest zerem wielomianu
takim, że
to
.
Pojęcie silnej stabilności nie wpływa na samą teorię zbieżności schematów, ale można się spodziewać, że schematy silnie stabilne zachowują się lepiej, por. rozdział 5.2.2.
W praktycznych obliczeniach wszystkie używane schematy liniowe wielokrokowe są silnie stabilne.
Możemy oczekiwać, że jednym z rozwiązań naszego równania różnicowego (5.1) będzie rozwiązanie stałe:
, czyli oczekujemy żeby
było pierwiastkiem wielomianu
, tzn.
. Ten warunek nazwiemy prezgodnością (ang. preconsistency) schematu.
Dodatkowo rozpatrzmy drugie proste równanie z warunkiem początkowym:
![]() |
którego rozwiązaniem jest .
Jeśli zastosujemy schemat do tego równania (zawsze
), to otrzymujemy następujące liniowe równanie różnicowe o stałych współczynnikach:
![]() |
Dodatkowo wprowadzimy wielomian:
![]() |
Zakładając, że schemat jest dokładny w dla tego zagadnienia początkowego tzn.
wstawiając
otrzymujemy:
![]() |
Jeśli założymy, że schemat jest prezgodny, to powyższe równanie daje nam:
![]() |
(5.3) |
Wprowadzamy pojęcie zgodności (konsystentności) schematu liniowego wielokrokowego:
Schemat liniowy wielokrokowy (4.8) jest zgodny (konsystentny), (ang.consistent) jeśli zachodzi (5.3).
Nie jest trudno pokazać, że ten warunek jest równoważny temu, że rząd schematu jest co najmniej jeden, tzn. prawdziwe jest następujące stwierdzenie:
Schemat liniowy wielokrokowy (4.8) jest zgodny wtedy i tylko wtedy, gdy rząd schematu wynosi co najmniej jeden.
Dowód pozostawiamy jako zadanie, por. ćwiczenie 5.14.
Jeśli rozwiązanie zagadnienia początkowego , schemat liniowy
-krokowy jest rzędu
dla
i jest stabilny, oraz wartości startowe
dla
spełniają nierówność:
![]() |
dla pewnej stałej nieujemnej niezależnej od
,
to ten
schemat jest zbieżny z rzędem
.
Zbieżność otwartego schematu Eulera jako schematu liniowego wielokrokowego.
Wiemy już, że rząd otwartego schematu Eulera wynosi jeden. Wystarczy więc zbadać stabilność schematu.
Wielomian zatem ma jeden pierwiastek
więc schemat jest stabilny, a nawet silnie stabilny.
Jeśli rozwiązanie jest klasy
, to zbieżność zachodzi z rzędem jeden, co potwierdzają wyniki eksperymentów z tabeli 4.1 w przypadku skalarnym dla równania
z
.
Rozpatrzmy schemat kroku środkowego (4.3), który jest dwukrokowym schematem liniowym rzędu dwa (zadanie). Nietrudno sprawdzić, że jest to schemat stabilny, ale nie silnie stabilny.
Rozpatrzmy nasze modelowe równanie liniowe:
![]() |
którego rozwiązanie jest równe
. Zastosujmy schemat midpoint biorąc za
dokładną wartość rozwiązania
. Narysujmy wykres rozwiązania przybliżonego otrzymanego tym schematem i schematem Eulera otwartym dla
na odcinku
, por. rysunek 5.1, i potem na
, por. rysunek 5.2.
W tym drugim przypadku rozwiązanie przybliżone otrzymane schematem midpoint od pewnego momentu przestaje zachowywać się jak rozwiązanie dokładne
, tzn. widzimy coraz większe oscylacje.
Weźmy mniejsze . Wtedy na
wszystko jest w porządku, ale na odcinku
znów rozwiązanie otrzymane schematem midpoint zaczyna od pewnego momentu oscylować z coraz większą amplitudą, zamiast maleć do zera. Im
jest mniejsze, tym odcinek na którym rozwiązanie otrzymane schematem midpoint dobrze aproksymuje rozwiązanie równania wyjściowego jest większy, ale zawsze w pewnym momencie pojawi się zaburzenie numeryczne, por. rysunki 5.3 i 5.4.
Numeryczne zaburzenie wynika właśnie z tego, że schemat nie jest silnie stabilny. Można wypisać wzór na rozwiązania równania różnicowego zadanego przez schemat midpoint dla tego równania i wykazać, że po jakimś czasie zawsze pojawią się oscylacje.
Wyniki tu otrzymane nie stoją w jakiekolwiek sprzeczności z teorią zbieżności schematów liniowych wielokrokowych.
Można sprawdzić, że schemat midpoint spełnia założenia tej teorii dla tego równania, i że
na ustalonym odcinku zachodzi zbieżność z rzędem dwa.
W praktyce nie należy stosować schematów, które nie są silnie stabilne, ponieważ zawsze mogą pojawić się oscylacje po a priori nie znanym czasie, szczególnie kiedy chcemy rozwiązać równanie na długim odcinku czasu.
W tym miejscu warto zwrócić uwagę, jak w praktycznych obliczeniach implementować schematy liniowe wielokrokowe.
Aby przy pomocy danego schematu -krokowego rzędu
obliczyć
, należy znać poprzednie
przybliżenia
. Czyli aby wystartować schemat musimy znaleźć odpowiednio dobre (tzn. takie, że
, por. Twierdzenie 5.2)
startowe przybliżenia
. Wartość
jest zadana jako warunek początkowy, tzn.
możemy przyjąć dany warunek początkowy
. Natomiast musimy wyliczyć
.
W praktyce do wyliczenia tych wartości możemy np. zastosować
razy schemat jednokrokowy tego rzędu co najmniej
.
Dla schematu Adamsa-Bashfortha dwukrokowego rzędu dwa za możemy przyjąć
obliczone jednym krokiem schematu Heuna.
Tu warto zauważyć że musimy zastosować tylko jeden krok schematu Heuna. Uwzględniając to widzimy, że
można by też obliczyć stosując dowolny schemat jednokrokowy rzędu jeden np. otwarty schemat Eulera.
Ogólniej do obliczenia startowych wartości
dla schematu
-krokowego rzędu
można zastosować schemat jednokrokowy rzędu
lub większego.
![]() |
![]() |
Rząd zbieżności schematu dla czasu można badać eksperymentalnie dla ustalonego zagadnienia początkowego
z warunkiem początkowym
i ze znanym rozwiązaniem
określonym na odcinku
.
Możemy wtedy przy pomocy naszego schematu, np. schematu Heuna, zmodyfikowanego schematu Eulera i otwartego schematu Eulera
obliczać przybliżone wartości
dla ustalonego
, z połowionym krokiem
:
dla ustalonego
. Będziemy badać jak zmienia się błąd w kolejnych krokach, a dokładniej, jak zmienia się stosunek błędów dla kolejnych połowionych
.
Jeśli błąd dla ustalonego zachowywałby się jak
, a dokładniej, gdyby
, to stosunek błędu powinien się zachowywać jak:
![]() |
dla dostatecznie małych ,
czyli dla schematów Heuna i zmodyfikowanego schematu Eulera jak cztery, a dla otwartego schematu Eulera jak dwa, a dla schematu rzędu cztery jak
- czyli szesnaście.
W tabeli 5.1 widzimy wyniki eksperymentu dla dla schematów: otwartego schematu
Eulera, Heuna i zmodyfikowanego schematu Eulera zastosowanych do zagadnienia początkowego
z
, którego rozwiązaniem jest
, czyli
.
W tabeli 5.2 przedstawiliśmy wyniki tego samego eksperymentu, ale dla schematu Rungego-Kutty rzędu cztery dla zadania początkowego z
i dla
, czyli z rozwiązaniem
. Skoro schemat jest rzędu cztery,
możemy oczekiwać, że błąd będzie jak
dla ustalonego
.
Wyniki uzyskane w tabeli 5.2 potwierdzają nasze przypuszczenie.
Dla kroku o połowę mniejszego błąd maleje około razy dla dostatecznie małych
(im
mniejsze, tym ten stosunek bliższy jest szesnastu). Przy okazji zauważmy, jak ogromna jest różnica w błędzie dla
schematu rzędu cztery (Rungego-Kutty), a dla schematu Heuna rzędu dwa.
W przypadku tego ostatniego -
błąd bezwzględny dla
jest rzędu
, a dla tego pierwszego błąd jest rzędu
.
Z kolei patrząc na błędy względne, w przypadku schematu Heuna błąd jest na poziomie
,
a dla schematu Rungego
, czyli poziom błędu w arytmetyce podwójnej precyzji praktycznie jest
wystarczająco dokładny.
Błędem bezwzględnym nazywamy , a względnym
.
W przypadku gdy
jest bardzo duże lub bardzo małe należy rozpatrywać błąd względny choćby z powodu
własności arytmetyki zmiennopozycyjnej.
Zbadaj rząd zbieżności zamkniętego schematu Eulera korzystając z teorii zbieżności dla schematów jednokrokowych.
Zbadaj rząd zbieżności zamkniętego schematu Eulera korzystając z teorii zbieżności dla schematów wielokrokowych liniowych. Czy schemat jest silnie stabilny?
Zbadaj rząd zbieżności schematu trapezów dany wzorem (4.7).
korzystając z teorii zbieżności dla schematów jednokrokowych,
korzystając z teorii zbieżności dla schematów wielokrokowych liniowych.
Czy schemat jest silnie stabilny?
Zbadaj rząd zbieżności zmodyfikowanego Eulera i schematu Heuna.
Zbadaj rząd zbieżności schematu punktu środkowego (ang. midpoint). Czy schemat jest silnie stabilny?
Zaimplementuj schematy: otwarty i zamknięty Eulera w octave i przetestuj rzędy zbieżności eksperymentalnie metodą połowionego kroku, jak opisano w rozdziale 5.4 dla równania
skalarnego oraz dla równania drugiego rzędu
z
.
Do rozwiązywania nieliniowego równania przy implementacji schematu zamkniętego w każdym kroku możesz użyć funkcji fsolve()
.
Zaimplementuj schemat trapezów, por. (4.7), w octave i przetestuj rząd zbieżności eksperymentalnie metodą połowionego kroku, jak opisano w rozdziale 5.4 dla równania
skalarnego oraz dla równania drugiego rzędu
z
.
Do rozwiązywania nieliniowego równania w każdym kroku możesz użyć funkcji fsolve()
.
Zaimplementuj zmodyfikowany schemat Eulera oraz schemat Heuna w octave i przetestuj rzędy zbieżności eksperymentalnie metodą połowionego kroku, jak opisano w rozdziale 5.4.
Zaimplementuj schemat Rungego rzędu cztery, dany wzorem (4.15) w octave i przetestuj rząd zbieżności schematu eksperymentalnie metodą połowionego kroku, jak opisano w rozdziale 5.4.
Zaimplementuj schemat (ang. midpoint) w octave i przetestuj rząd zbieżności schematu eksperymentalnie metodą połowionego kroku,
jak opisano w rozdziale 5.4, przyjmując, że
, czyli jest równe dokładnemu rozwiązaniu. Dla równania
z
zastosuj schemat na długim odcinku czasu dla ustalonego
. Czy rozwiązanie zachowuje się tak jak tego oczekujemy? Zmniejsz
i zobacz czy sytuacja się poprawia?
Zbadaj rząd, stabilność i rząd zbieżności otwartego dwukrokowego schematu Adamsa, por. Ćwiczenie 4.7.
Zaimplementuj w octave otwarty dwukrokowy schemat Adamsa, por. Ćwiczenie 4.7.
Następnie przetestuj eksperymentalnie rząd zbieżności tego schematu metodą połowionego kroku,
jak opisano w rozdziale 5.4 dla z
. Biorąc za
:
czyli rozwiązanie dokładne dla
.
obliczone schematem Heuna czyli schematem Rungego-Kutty rzędu dwa.
obliczone schematem Eulera otwartym czyli schematem rzędu jeden.
Znajdź dokładne rozwiązanie schematu różnicowego, który jest schematem punktu środkowego zastosowanym dla równania
. Pokaż, że jeśli
, a
różne od jednej konkretnej wartości (co np. odpowiada warunkom startowym
), to
, czyli pojawią się niepotrzebne kumulujące się zaburzenia numeryczne.
Udowodnij stwierdzenie 5.1.
Z tego, że rząd schematu jest co najmniej jeden, wynika od razu warunek zgodności schematu (wystarczy rozważyć oba równania z
i
z
). Natomiast aby pokazać, że ze zgodności schematu wynika to, że rząd schematu jest większy od jeden,
należy zastosować wzór Taylora.
Czy schemat Adamsa (por. rozdział 4.1.2) może nie być stabilny, ewentualnie nie być silnie stabilny? Uzasadnij podając ewentualnie kontrprzykład niestabilnego (czy nie będącego silnie stabilnym) schematu Adamsa.
Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.
strona główna | webmaster | o portalu | pomoc
© Wydział Matematyki, Informatyki i
Mechaniki UW, 2009-2010. Niniejsze materiały są udostępnione bezpłatnie na licencji Creative Commons Uznanie autorstwa-Użycie niekomercyjne-Bez utworów zależnych 3.0 Polska.
Projekt współfinansowany przez Unię Europejską w ramach
Europejskiego Funduszu Społecznego.
Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.