Jawny schemat Eulera Czy bd cakowity maleje gdy

Jawny schemat Eulera Czy błąd całkowity maleje gdy Dt maleje ? Czy maleje do zera? eksperyment numeryczny problem początkowy: u’=lu, u(0)=1 l=-100 Dt=0. 001 dokładny jawny Euler e (błąd globalny) = dokładny - numeryczny z rozwiązaniem dokładnym u(t)=exp(lt)

rząd zbieżności (dokładności) określa jakość metody: jak szybko błąd globalny dla chwili T zmierza do zera w funkcji Dt jawna metoda Eulera = pierwszy rząd dokładności O(Dt) jest to minimalny rząd dokładności dla użytecznej metody zmniejszajmy krok czasowy, jaki wynik w chwili t=0. 01 ? [1/e=0. 3678794] n 10 102 103 104 Dt 10 -3 10 -4 10 -5 10 -6 un 0. 34867 0. 36603 0. 36769 0. 36784 exp(-1)-un 1. 920 10 -2 1. 847 10 -3 1. 840 10 -4 1. 839 10 -5

Definicja: Metody różnicowa jest zbieżna jeśli błąd globalny znika do zera w chwili T gdy z Dt do 0 zbieżność a błędy zaokrągleń (skończona dokładność arytmetyki zmiennoprzecinkowej)

błędy zaokrągleń a zbieżność do tej pory zakładaliśmy, że błędy zaokrągleń nie ma (że arytmetyka dokładna) arytmetyka zmiennoprzecinkowa nie jest dokładna. pojedyncza precyzja: 32 bity podwójna : 64 bity arytmetyka 21 – bitowa błąd minimalny zmniejszanie kroku czasowego nie poprawi już wyniku

błędy zaokrągleń a metody różnicowe rozwiązanie równania różniczkowego w chwili tn rozwiązanie równania różnicowego z dokładną arytmetyką rozwiązanie uzyskane z arytmetyką skończonej dokładności błąd całkowity błąd globalny (jak wcześniej zdefiniowano) błąd zaokrąglenia oszacowanie od góry błędu całkowitego

błędy zaokrągleń a metody różnicowe e ani ó dg ry uc d ę bł eg t i ow ałk o o błąd cow a z os błąd globalny dla schematu Eulera błąd zaokrągleń rzędu liczby wykonanych kroków, czyli 1/dt dt optymalny krok czasowy błędy zaokrągleń dają o sobie znać gdy wykonamy zbyt wiele kroków remedium: używać się schematów o wyższym rzędzie zbieżności niż pierwszy

błędy zaokrągleń a metody różnicowe oszacowanie od góry błędu całkowitego błąd globalny dla schematu Eulera błąd całkowity błąd zaokrągleń błąd globalny błąd zaokrąglenia optymalny krok czasowy Definicja: Metody różnicowa jest zbieżna jeśli błąd globalny znika do zera w chwili T gdy z Dt do 0 uwaga: definicja zbieżności dotyczy błędu globalnego a nie całkowitego

Wróćmy do eksperymentu numerycznego i zwiększmy krok czasowy do Dt=0. 05 problem początkowy: u’=-100 u, u(0)=1 z rozwiązaniem dokładnym u(t)=exp(-100 t) iteracja się rozbiega tn 0 0. 05 0. 15 0. 25 0. 3 un 1 -4 8 -16 256 -1024 4096 wniosek: wyniki metody zbieżnej mogą eksplodować dla zbyt dużego kroku czasowego

bezwzględna stabilność schematu różnicowego* schemat różnicowy dla du/dt = f (dla danego f) i dla danego kroku czasowego jest bezwzględnie stabilny jeśli kolejne generowane przez niego wartości pozostają skończone. stabilność bezwzględna: określana zazwyczaj dla autonomicznego problemu liniowego bezwzględna stabilność jawnej metody Eulera wsp. wzmocnienia wyniki będą skończone, jeśli

z= l. Dt region bezwzględnej stabilności metody: zbiór par : krok czasowy Dt oraz zespolone wartości l dla których metoda jest bezwzględnie stabilna region bezwzględnej stabilności jawnej metody Eulera: zbiór punktów odległych od (-1, 0) o nie więcej niż 1 Dt Im (l) koło o środku w (-1, 0) i promieniu 1 -2 -1 1 -1 Dt Re(l) niestabilność bezwzględna metody dla prawej połowy p. Gaussa = nic dziwnego rozwiązania dokładne eksploduje do nieskończoności gdy t .

jawna metoda Eulera jest bezwzględnie stabilna jeśli |1+l. Dt| 1 czyli dla rzeczywistego l <0 : Dt -2/ l Weźmy l = -1, u(0)=1, rozwiązanie u(t)=exp(-t), warunek stabilności bzwz: Dt 2 Przepis Eulera: un+1=un(1 -Dt) 6 u(t) 0. 4 4 dokładny Dt=0. 5 Dt=0. 9 u(t) 0. 6 Dt=1 : wszędzie 0 0. 2 Dt=1. 2 Dt=2. 5 2 0 -2 0. 0 0 2 4 t 6 8 10 -4 0 2 4 t 6 8 10 uwaga: rozwiązanie bezwzględnie stabilne (np. Dt=1 lub Dt=2) może być bardzo niedokładne lub wręcz - jakościowo złe = tutaj stałe i niemonotoniczne odpowiednio

pozbyć się ograniczenia na krok czasowy: niejawna metoda Eulera metoda jawna metoda Eulera (funkcjonuje jak podstawienie) „metoda odważna” . . niejawna funkcjonuje jak równanie nieliniowe „metoda ostrożna” zmiana u zgodna z prawą stroną w punkcie docelowym

niejawna metoda Eulera: niejawny Euler problem początkowy: u’=-100 u, u(0)=1 z rozwiązaniem dokładnym u(t)=exp(-100 t) jawny Euler tn 0 0. 05 0. 15 0. 25 0. 3 un 1 -4 8 -16 256 -1024 4096 tn 0 0. 05 0. 15 0. 25 0. 3 un 1. 166(6). 027(7). 004(629). 0007716. 0001286. 00002143 e(tn) 0 -. 15992 -. 02773 -. 00462 itd. . exp(-100 tn) gaśnie znacznie szybciej niż 1/6 n mało dokładne, ale zawsze to lepiej niż eksplodująca oscylacja jawnego Eulera

niejawna metoda Eulera: region bezwzględnej stabilności Dt Im (l) 1 Dt Re(l) -1 rejon bezwzględnej stabilności: dopełnienie pustego koła o środku w (1, 0) i promieniu 1

regiony stabilności metod Eulera Dt Im (l) -2 -1 1 1 -1 metoda Eulera jawna Dt Re(l) -1 niejawna metoda Eulera

niejawna metoda Eulera: region bezwzględnej stabilności Dt Im (l) 1 Dt Re(l) -1 Re(l)<0 : niejawny Euler bezwzględnie stabilny dla dowolnego kroku czasowego takie metody: tzw. A-stabilne dla Re(l)>0, poza kołem metoda Eulera jest bezwzględnie stabilna mimo, że rozwiązania równania różniczkowego są niestabilne (patrz wyżej) w tym obszarze metoda jest nadstabilna daje skończone wartości, mimo że rozwiązania dokładne dąży do nieskończoności bezwzględna stabilność nie oznacza dobrej dokładności. W regionie stabilności dla Re(l)>0 błędy będą rosły w nieskończoność.

jak rozwiązać, gdy nie można rozwikłać równania (f nieliniowe względem u) lub gdy f nieznane w formie wzoru 1) iteracja funkcjonalna iterować do zbieżności jeśli się zbiegnie um=um-1 i mamy rozwiązanie równania nieliniowego

iteracja funkcjonalna przykład problem początkowy: u’=-100 u, u(0)=1, dt=0. 05 z rozwiązaniem dokładnym u(t)=exp(-100 t) tn 0 0. 05 0. 15 0. 25 0. 3 un 1. 166(6). 027(7). 004(629). 0007716. 0001286. 00002143 kolejne oszacowania: 1, -4, 21, -104, 521, -2604, . . . iteracja się nie zbiega cały zysk z niejawności stracony bo nie potrafimy wykonać kroku e(tn) 0 -. 15992 -. 02773 -. 00462

iteracja funkcjonalna przykład iteracja się nie zbiega . zmniejszymy krok dt, zaczynając iterację od un-1 będziemy bliżej rozwiązania. Może się zbiegnie. dt=0. 01 dt=0. 001 (1, 0, 1, 0) 1, 0. 91, 0. 9091, 0. 909091, . . . 0. 9090909 iteracja funkcjonalna się zbiega gdy Dt u’=-100 u, max|fu(t, u)| 1 (w otoczeniu u) Dt 100 < 1 uwaga: w tej sytuacji jawny Euler jest bezwzględnie stabilny dla 2 -krotnie większego kroku! [dla jawnego Eulera Dt 100 < 2] Z iteracją funkcjonalną stosować wstecznego Eulera nie ma sensu.

problem początkowy: u’=-100 u, u(0)=1, dt=0. 05 z rozwiązaniem dokładnym u(t)=exp(-100 t) 1, -4, 21, -104, 521, -2604, . . . oscylująca rozbieżność - stłumimy ją: iteracja funkcjonalna – zapewniamy zbieżność modyfikując przepis iteracyjny zamiast: „mieszając” nowe i stare rozwiązania z wagą w, 0 w 1 Zabieg podobny do “podrelaksacji” jeśli się zbiegnie – to do rozwiązania schematu niejawnego

problem początkowy: u’=-100 u, u(0)=1, dt=0. 05 z rozwiązaniem dokładnym u(t)=exp(-100 t) iterujemy u(dt) w=0. 1 w=0. 3 w=0. 2 wybierając w odpowiedni sposób wagę w: potrafimy ustabilizować iterację i doprowadzić ją do zbieżności.

dt=0. 01 w=0 (1, 0, 1, 0) w=. 2 (optymalne dla dt=0. 05) 0. 8, 0. 608, 0. 5648, 0. 53888, 0. 5233, 0. 51399, 0. 50839, 0. 50503, 0. 5018, 0. 5010, 0. 50065, 0. 50039, . . . , 1/2 tutaj optymalne byłoby w=1/2 (zbieżność w jednej iteracji) dla w=. 7 0. 3, 0. 58, 0. 468, 0. 512, 0. 4948, 0. 5003, 0. 4998 dt=0. 001 w=1 1, 0. 91, 0. 9091, 0. 909091, . . . 0. 9090909 w=0. 2

dt=0. 01 w=0 (1, 0, 1, 0) w=. 2 (optymalne dla dt=0. 05) 0. 8, 0. 608, 0. 5648, 0. 53888, 0. 5233, 0. 51399, 0. 50839, 0. 50503, 0. 5018, 0. 5010, 0. 50065, 0. 50039, . . . , 1/2 tutaj optymalne byłoby w=. 5 (zbieżność w jednej iteracji natychmiastowa) dt=0. 001 w=1 1, 0. 91, 0. 9091, 0. 909091, . . . 0. 9090909 w=0. 2 dla w=. 7 0. 3, 0. 58, 0. 468, 0. 512, 0. 4948, 0. 5003, 0. 4998 Problem: 1) w trzeba odpowiednio dobrać (mniejszy krok czasowy, w bliższe 1) 2) dla źle dobranego w iteracja może być wolnozbieżna Proces optymalizacji (np. dynamicznej) w może być kłopotliwy.

Na szczęście nie jesteśmy skazani na iterację funkcjonalną 2) metoda Newtona-Raphsona szukamy zera równania nieliniowego F(x) F(xn+Dx)=F(xn)+Dx F’(xn) F(xn+Dx)=0 xn+1=xn-F(xn)/F’(xn)

2) metoda Newtona-Raphsona szukamy zera równania nieliniowego F(xn+Dx)=F(xn)+Dx F’(xn) F(xn+Dx)=0 xn+1=xn-F(xn)/F’(xn)

niejawny schemat Eulera z metodą Newtona-Raphsona, zastosowanie problem początkowy: u’=-100 u, u(0)=1 z rozwiązaniem dokładnym u(t)=exp(-100 t) kolejne przybliżenia: Dt=0. 05 (jawny Euler stabilny bezwzględnie dla Dt <0. 02) 1, 0. 1666677 zbieżność w jednej iteracji - f jest liniowa w u zakres zbieżności: w praktyce Dt znacznie większy niż w iteracji funkcjonalnej ale: niedostępne proste oszacowane przedziału zbieżności w praktyce iteracja Newtona – szybsza i szerzej zbieżna niż iteracja funkcjonalna

niejawny schemat Eulera z metodą Newtona-Raphsona, zastosowanie dla problemu nieliniowego problem początkowy: dla równania: u’=u(u-1)

czerwone niejawny Euler z krokiem Dt=0. 1 u(0)=0. 8 iteracja dla u(Dt) ze startem w u(0): 0. 80000 0. 78297 0. 78300

niejawny schemat Eulera z metodą Newtona-Raphsona gdy przepis funkcyjny nieznany (np. programujemy metodę dla dowolnego f ) można szacować z ilorazu różnicowego (poniżej centralny = dokładnie różniczkuje parabole) cena zastąpienia dokładnej pochodnej ilorazem różnicowym? przy osiągniętej zbieżności - nie zmieni rozwiązania! może tylko spowolnić iterację! dla naszego przykładu u’=u(u-1) centralny iloraz różnicowy zadziała dokładnie dla dowolnego du żeby przykład był ciekawszy: policzmy pochodną z wstecznego ilorazu różnicowego

u(0)=0. 8, pierwszy krok t=Dt: u’=f(u)=u(u-1) metoda Newtona dla pochodnej f liczonej numerycznie w każdej iteracji: dokładna pochodna iloraz wsteczny du=u/10 iloraz wsteczny du=u/2 0. 80000 0. 78297 0. 78300 0. 80000 0. 78312 0. 78301 0. 78300 0. 80000 0. 78367 0. 78303 0. 78301 0. 78300 przybliżenie w liczeniu pochodnej nie zmienia wyniku do którego iteracja zbiega bo: xn+1=xn-F(xn)/F’(xn) nieco spowalnia iterację numeryczne liczenie pochodnych w każdej iteracji może być kosztowne w praktyce można np. wstawić tutaj un-1 można również używać oszacowania pochodnej w wielu kolejnych iteracjach odnawiać pochodną gdy iteracja zwalnia

zamiast iteracja Newtona z pochodną liczoną w poprzednim kroku (nieiterowaną) dla naszego przykładu: f(u)=u(u-1) z dt=0. 1 dt=0. 5 pochodna brana z z punktu tn-1, un-1 iterowan poprzednieg a o kroku 0. 8 0. 6857 0. 80000 0. 6937 0. 6857 0. 78297 0. 6937 0. 6950 0. 78300 0. 6935 0. 78300 0. 6938 bez różnicy! wolni 0. 6937 ej w mianowniku: 1 -dt(2 u-1) stara: 0. 7 stara: 0. 94 , doiterowana 0. 9434 doiterowana: 0. 806 iterowana pochodna dt=1 iterowan stara a 0. 8 0. 4 0. 8 0. 5333 0. 4 0. 5523 0. 8 0. 5527 brak zbieżności stara 0. 4 doiterowana 0. 89

jawny schemat Eulera [globalny błąd O(Dt)] u(t) niejawny schemat Eulera [globalny błąd O(Dt)] f(t, u) u(t) [t, u(t)] dokładne f(t, u) [t+ Dt, u(t+Dt)] Dt t Dt u(t) t przesunięcie wyliczane na podstawie średniej arytmetycznej z chwil t i t+Dt [wzór trapezów] f(t, u) [t, u(t)] dokładne [t+ Dt, u(t+Dt)] Dt t

dokładność wzóru trapezów a jawnego schematu Eulera: eksperyment Równanie: Warunek początkowy: u 1=u(t 1=0)=1 Rozwiązanie: Punkt t 2=0. 5 u 2 = ? [dokładnie: 1. 1331] u 2 Euler jawny jeden krok: u 2=u 1+Dt t 1 u 1= 1 wzór trapezów u 2=u 1+(t 1 u 1+t 2 u 2) Dt /2 = u 1+t 2 u 2 Dt /2 jawny Euler u 2: = u 1+u 2 /8 iteracja funkcjonalna wynik 8/7 niejawny Euler: 4/3=1. 333 (pokazać) nieco gorzej niż jawny TR wygląda na bardziej dokładny od E:

Oszacować błąd lokalny wzoru trapezów 1. rozw. Taylora wstecz 2. dla dowolnej funkcji ciągłej f(t)=f(t+Dt)+O(Dt) (wstawimy, rząd błędu pozostanie trzeci) 3. Rozwiązać na u(t+Dt) Euler miał rząd obcięcia Dt 2 pozbyć się go.

3. Rozwiązać na u(t+Dt) [przepisane] 4. Uśrednić z rozwinięciem Taylora do przodu 5. Wynik 6. Korzystamy z równania jawny i niejawny Euler wzór trapezów – lokalny błąd rzędu drugiego (rząd dokładności 1) – lokalny błąd rzędu trzeciego (rząd dokładności 2)

stabilność bezwzględną wzoru trapezów problem modelowy: WP: u(t=0)=1. rozwiązanie u=exp(lt) zbiór punktów na p. Gaussa, które są nie dalej od (-2, 0) niż od (0, 2)

region bzwz. stabilności wzoru trapezów Dt Im (l) Dt Re(l) Wniosek: dla l<0 wzór trapezów bezwzględnie stabilny dla dowolnego kroku czasowego ! A-stabilny druga bariera Dahlquista: maksymalny rząd dokładności metody A-stabilnej =2 schemat trapezów jest najdokładniejszą metodą A-stabilną spośród liniowych metod wielokrokowych Implementowana np. w SPICE.

region bzwz. stabilności Eulera: koło o promieniu 1 Dt Im (l) i środku (-1, 0) Dt Re(l) region bzwz. stabilności wzoru trapezów Dt Im (l) niejawna metoda Eulera: region bezwzględnej stabilności Dt Re(l) Dt Im (l) między metodami można przechodzić w sposób ciągły 1 Dt Re(l) -1 q=0, 1, 1/2 – Euler jawny, niejawny i wzór trapezów odpowiednio w wykładzie na temat niejawnych formuł RK zobaczymy, że dokładność rzędu 2 uzyskana tylko dla q=1/2 region stabilności ?

iteracja funkcjonalna a wzór trapezów problem początkowy: u’=-100 u, u(0)=1 z rozwiązaniem dokładnym u(t)=exp(-100 t) Dt=0. 01 = graniczny dla zbieżności IF dla niejawnego Eulera 1, 0, 0. 5, 0. 25, 0. 375, 0. 3125, 0. 34375, 0. 328125, 0. 33593, 0. 33203, 0. 333984, . . . , 0. 333333 wzór trapezów = używa prawej strony z poprzedniego kroku czasowego z wagą 0. 5 – co nieco stabilizuje iterację. niestety iteracja funkcjonalna dla Dt=0. 02 już przestaje być zbieżna (+1, -1, itd. . ) wzór trapezów zwiększa zakres zbieżności iteracji dwukrotnie (wyraz podkreślony stabilizuje iteracje) ale to wciąż mało metoda Newtona-Raphsona pozostaje

poznane metody: 1) 2) 3) ) Poznane metody: jednokrokowe (1 -3), jawna (1) i niejawne (2 -3), pierwszego (1 -2) i drugiego (3) rzędu dokładności Metody (2 -3) A stabilne, metoda (2) nadstabilna jawne metody różnicowe wysokiej dokładności ? ?

jawne metody jednokrokowe wyższego rzędu dokładności niż jawny Euler u’=f(t, u), u(0)=u 0 rozwinięcie Taylora ponownie: liczymy pochodne: z RR. RR różniczkujemy po czasie czyli podobnie Zależnie od tego gdzie się zatrzymamy uzyskamy błąd lokalny zadanego rzędu

Zależnie od tego gdzie się zatrzymamy uzyskamy błąd lokalny zadanego rzędu np. pomysł: mało przydatny w praktyce ze względu na konieczność analitycznego wyliczenia pochodnych cząstkowych f. Dla metod ogólnych: nie powinniśmy liczyć, że f jest dane wzorem

podejście alternatywne: inspirowane całkowaniem prawa strona = funkcja tylko t z rozwiązaniem: jeśli zastąpimy całkę kwadraturą prostokątów z wywołaniem funkcji w lewym końcu przedziału tn-1 tn u(tn)=u(tn-1)+Dt f(tn-1) + O(Dt 2) - rozpoznajemy jawny schemat Eulera kwadratura prostokątów z wywołaniem funkcji w prawym końcu przedziału u(tn)=u(tn-1)+Dt f(tn) + O(Dt 2) - rozpoznajemy niejawny schemat Eulera tn-1 tn kwadratura trapezów u(tn)=u(tn-1)+Dt f(tn)/2+ Dt f(tn-1)/2 + O(Dt 3) - rozpoznajemy niejawny schemat trapezów tn-1 tn

reguła punktu środkowego wzór prostokątów z wywołaniem funkcji w środku przedziału (dokładny dla funkcji liniowej, znoszenie błędów) tn-1/2 tn uogólniony wzór na równanie równania u’=f(t, u) np. ze schematu Eulera: ale - skąd rozwiązanie w środku przedziału? błąd lokalny Eulera O(Dt 2), czy reguła punktu środkowego zachowa trzeci rząd błędu lokalnego?

sprawdźmy to rozważając bardziej ogólny schemat: obliczone na początku kroku obliczone gdzieś w środku przedziału (tn-1, tn) z odpowiednio oszacowanym rozwiązaniem u dla tego t (wzór typu Eulera) jest to jawny dwustopniowy schemat Rungego-Kutty. potencjalna wyższa dokładność od jawnego Eulera kosztem dwóch wywołań f (podobnie jak we wzorze trapezów, ale RK: jawny) b 1, b 2, a, c – parametry metody –jakie muszą być aby RK 2 (2 = rząd dokładności) reguła punktu środkowego: należy do tej klasy z b 1=0, b 2=1, c=1/2, a =1/2

Jawne metody Rungego-Kutty dwustopniowe: wybór parametrów (*) jak dobrać b 1, b 2, c, a ? – metodą brutalnej siły tak aby rozwinięcie Taylora metody zgadzało się z rozwinięciem Taylora dokładnego równania różniczkowego do wyrazów tak wysokiego rzędu jak to tylko możliwe u’=f(t, u) przypominamy: rozwinięcie Taylora dla funkcji dwóch zmiennych wstawiamy rozwiązanie dokładne u(tn), u(tn-1) do (*) i rozwijamy względem tn-1, un-1

to trzeba rozwinąć wstawmy k 2 do rozwinięcia. Zachowajmy człony do Dt 2: (wszystko liczone w tn-1, un-1) rozwinięcie Taylora rozwiązania dokładnego uzyskaliśmy kilka slajdów wcześniej czyli: rząd Dt: b 1+b 2=1, rząd Dt 2: b 2 c=b 2 a=1/2 czyli reguła punktu środkowego: b 1=0, b 2=1, c=1/2, a =1/2 ma błąd lokalny rzędu O(Dt 3) mamy metodę równie dokładną co wzór trapezów – ale jawną (co ma swoje zalety i wady) Wyższy rząd błędu do uzyskania tylko w metodach o większej niż 2 liczbie stopni

cztery parametry i trzy równania b 1+b 2=1 b 2 c=b 2 a=1/2 - pozostaje swoboda w wyborze parametrów reguła punktu środkowego RK 2 b 1=0, b 2=1, c=1/2, a =1/2 dwa zastosowania jawnego schematu Eulera oszacowanie wstępne w punkcie pośrednim (błąd lokalny rzędu drugiego) oszacowanie docelowe (błąd lokalny oszacowania: rzędu trzeciego) albo (przesunięty indeks) u(t) [t, u(t)] dokładne 1) Szacujemy metodą Eulera punkt środkowy [t+Dt/2, u(t+Dt/2)] korzystając z f(t, u) w lewym końcu przedziału 2) Wykorzystujemy wartość f w tym punkcie do wyliczenia zmiany y na całym przedziale Dt [t+Dt/2, y(t+Dt/2)] Dt t

RK punktu środkowego: b 1+b 2=1, b 2 c=b 2 a=1/2 inny wybór: b 1=b 2=1/2, wtedy musi a=c=1 metoda podobna do wzoru trapezów (ale jawna) u(t) 1) Szacujemy metodą Eulera punkt końcowy [t+Dt, u(t+Dt)] korzystając z f(t, u) w lewym końcu przedziału 2) krok z t do t+Dt wykonujemy biorąc średnią arytmetyczną z f na początku i końcu [t, u(t)] dokładne Dt metoda RK 2 trapezów t

dla błędu lokalnego O(Dt 3) potrzeba aby, rząd Dt: b 1+b 2=1, rząd Dt 2: b 2 c=b 2 a=1/2 punkt środkowy b 2=1, b 1=0 [b 1+b 2]=1 czy ma sens b 1=1, b 2=0 ?

Metody Rungego-Kutty, forma ogólna są to metody jednokrokowe, czyli można zapisać: metoda RK w s-odsłonach (stage) (unikamy słowa „krok”) z wzory przedstawiane w formie tabel Butchera c. A b

Metody Rungego-Kutty, forma ogólna czasem zapisywane w postaci: tutaj Ui – przybliżone rozwiązanie w chwili tn-1+ci. Dt zazwyczaj niższej dokładności niż rozwiązanie końcowe

jawne metody Rungego-Kutty jawne: aij=0 dla j i obcięte sumowanie: odsłona i-ta wyliczana na podstawie tylko wcześniejszych odsłon historycznie wszystkie RK były jawne, uogólnienie okazało się przydatne dla problemów sztywnych

Wyprowadzanie formuł RK (a, b, c) 1) 2) 3) Rozwijamy rozwiązanie dokładne w szereg Taylora względem tn-1 Podstawiamy rozwiązanie dokładne do ogólnej formy RK i rozwijamy względem tn-1 Wartości parametrów a, b, c uzyskujemy z porównania. zazwyczaj w sposób niejednoznaczny najbardziej popularne: jawne formuły 4 -etapowe RK 4: o 4 -tym stopniu zbieżności (4 -tym rzędzie dokładności) i 5 -tym rzędzie błędu lokalnego ogólna tabela Butchera: dla jawnych RK 4 c 1=0 (dla każdej jawnej RK, zaczynamy – k 1 od wyliczenia prawej strony w kroku początkowym)

klasyczna formuła RK 4: k 2 u’ u(t) k 1 k 3 4 wywołania f na krok, błąd lokalny O(Dt 5) gdy f tylko funkcja czasu RK 4 redukuje się do formuły Simpsona (dokładnie całkuje wielomiany trzeciego stopnia): k 4

Jawne schematy RK dla układu równań różniczkowych 2 zmienne zależne u 1, u 2, 2 prawe strony f 1, f 2 2 równania, s-odsłon (i=1, 2, . . . , s) zapis wektorowy un-1, un , f, U 1, U 2, . . . UN są wektorami o 2 składowych

Tabela Butchera dla klasycznej jawnej RK 4 0 0 0 1/2 0 0 1 0 1/6 1/3 1/6

Dlaczego RK 4 najbardziej popularna: Liczba kroków a rząd zbieżności jawnych metod RK: rząd 1 2 3 4 5 6 7 8 minimalna liczba odsłon 1 2 3 4 6 7 9 11 RK 4 – wyjątkowo opłacalna RK 1 – metoda RK w jednej odsłonie b 1+b 2=1, przy b 1=1, b 2=0 dostaniemy jawnego Eulera warunek a*b 2=c*b 2 =1/2 nie będzie spełniony jawny schemat Eulera to jawna metoda RK 1

jawny Euler tabela Butchera RK 2 trapezów RK 2 punktu środkowego 0 0 0 1/2 0 0 1 b 1=b 2=1/2, a=c=1

Szacowanie błędu lokalnego w metodach jednokrokowych Po co? 1) W rachunkach numerycznych musimy znać oszacowanie błędu 2) Aby ustawić krok czasowy tak, aby błąd był akceptowalny 3) Gdy oszacowanie jest w miarę dokładne: można poprawić wynik

Oszacowanie błędu lokalnego w metodach jednokrokowych W każdym kroku generujemy nowy błąd w rachunkach. Znamy jego rząd. dla RK: wstawialiśmy rozwiązanie dokładne do schematu i je rozwijaliśmy w szereg T. wybór b 1=0, b 2=1, c=1/2, a =1/2 dawał RK 2 punktu środkowego Rozwijając do jednego rzędu wyżej z Dt uzyskamy oszacowanie błędu lokalnego dn= u(tn) - un [przy założeniu, że u(tn-1) = un-1] świetny wzór choć mało praktyczny

Oszacowanie błędu (lokalnego) w metodach jednokrokowych metodą rzędu p z chwili tn-1 wykonujemy krok do tn może zależeć od tn-1 oraz un-1, ale nie zależy od Dt folia wcześniej szacowanie błędu: 1) ekstrapolacja Richardsona (step doubling) 2) osadzanie (embedding)

ekstrapolacja Richardsona dwa kroki Dt: dostaniemy lepsze oszacowanie u(tn+1) tn-1 Dt tn tn+1 Dt jeden krok 2 Dt: dostaniemy gorsze oszacowanie u(tn+1) tn-1 2 Dt tn+1 szacujemy Cn z porównania obydwu rozwiązań

ekstrapolacja Richardsona błąd lokalny u(tn)-un=dn jest: wykonujemy krok następny od tn do tn+1 odchylenie wyniku numerycznego od dokładnego u(tn+1)-un+1= g dn +dn+1 1) zakładamy, że krok jest na tyle mały, że stała błędu się nie zmienia Cn Cn+1 (lub, ze w jednym kroku zmienia się o O(Dt)] wtedy błąd lokalny popełniony w chwili tn+1 jest dn+1 dn. 2) gdy krok mały: współczynnik wzmocnienia błędu g 1 [np. du/dt=-lu jawny E. (1 -l. Dt)] (błąd popełniony w kroku pierwszym nie jest istotnie wzmacniany) Przy tym założeniu: błąd po drugim kroku – suma błędów g dn+dn+1 2 dn,

ekstrapolacja Richardsona tn-1 Dt tn tn+1 Dt to chwili tn+1 dojdziemy z tn-1 w pojedynczym kroku 2 Dt tn-1 tn+1 2 Dt dostaniemy gorsze oszacowanie u(tn+1) chcemy poznać Cn (to + znajomość p da nam oszacowanie błędu): odejmujemy niebieskie wzory tak aby wyeliminować rozwiązanie dokładne (nam niedostępne)

ekstrapolacja Richardsona błąd wykonany po dwóch krokach Dt wynosi więc: pierwszy wniosek: jeśli znamy rząd metody p to potrafimy go podnieść o jeden

ekstrapolacja Richardsona podnosimy rząd dokładności metody „algorytm”

Przykład: ekstrapolacja Richardsona r. dokładne: Euler (p=1) błąd lokalny O(Dt 2) Euler po poprawce: błąd lokalny O(Dt 3) kreski: RK 2 punktu środkowego (p=2), b. lok. O(Dt 3) znając rząd dokładności możemy radykalnie poprawić dokładność metody przy natdatku (50 procent) numeryki

Euler ekstrapolacja Richardsona RK 2 z odciętym błędem

Oszacowanie błędu lokalnego w metodach jednokrokowych 1) ekstrapolacja Richardsona (step doubling) 2) osadzanie (embedding) cel: szacujemy błąd lokalny metody rzędu p przy pomocy lepszej metody, np. rzędu p+1 obydwie metody szacują rozwiązanie w tych samych chwilach czasowych co daje oszacowanie błędu gorszej metody nie nadaje się do poprawiania schematu p po cóż zresztą poprawiać gdy mamy p+1

celem szacowania błędu nie jest poprawa wyniku, (dla poprawy zawsze można Dt zmienić) lecz adaptacja Dt : stały krok zawsze może okazać się zbyt wielki albo zbyt mały. JAKI KROK CZASOWY SYMULACJI USTAWIĆ gdy coś ciekawego zdarza się tylko czasem ?

Automatyczna kontrola kroku czasowego dla metod jednokrokowych Program może sam dobierać krok czasowy w zależności od tego co dzieje się w symulacji. Chcemy utrzymać błąd na poziomie zbliżonym do parametru tol. nie większy aby zachować wymaganą dokładność, nie mniejszy aby nie tracić czasu na rachunki zbyt dokładne Szacujemy błąd lokalny E (ekstrapolacja Richardsona lub metody embedding) E=C[Dt]p+1 chcemy zmienić krok odpowiednio do naszych wymagań z Dt do Dt(nowy) tol=C[Dt(nowy)]p+1 Dt(nowy)=(S tol /E)1/(p+1) Dt Dt(nowy)=(tol/E)1/(p+1) Dt dla bezpieczeństwa S<1 wzór zwiększy zbyt mały krok i vice versa uwaga: błąd jest szacowany, zawsze warto dorzucić sztywne ograniczenia na Dt

Automatyczna kontrola kroku czasowego dla metod jednokrokowych symulacja ustawiająca krok czasowy może wyglądać np. tak: u 0= warunek początkowy t 0=0 n=1 do { jeśli E<tol { tn: =tn+Dt n: =n+1 (oznacza akceptację wyniku) } Dt: =(S tol /E)1/(p+1) Dt } while ( t<T)

Przykład: oscylator harmoniczny używane oszacowanie błędu z RK 2 Uwaga: tutaj rozwiązania nie poprawiamy przez ekstrapolacje tolerancja błędu obcięcia tol=0. 1 RK 2 start tol=1 e-2 tol=1 e-3 algorytm ustawia minimalny krok czasowy gdy zmiany prędkości lub położenia są maksymalne x 2 V 2 Dt

wyniki Konrada Rekiecia RK 4 – spirala się skręca zamiast rozkręcać E(t) tol=. 1 RK 2 RK 4 przy założonej tolerancji RK 4 wcale nie jest dokładniejsze od RK 2

. . . tylko pozwala stawiać dłuższe kroki dt