1 00000 1/2 000 0 00 10010 1/61/3 1/6 Tabela Butchera dla klasycznej jawnej RK4

2 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 RK4 – wyjątkowo opłacalna Dlaczego RK4 najbardziej popularna: RK1 – metoda RK w jednej odsłonie b1+b2=1, przy b 1 =1, b2=0 dostaniemy jawnego Eulera warunek a*b2=c*b2 =1/2 nie będzie spełniony jawny schemat Eulera to jawna metoda RK1

3 jawny Euler RK2 trapezów tabela Butchera b 1 =b 2 =1/2, a=c=1 RK2 punktu środkowego 000 1/2 0 01

4 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

5 Oszacowanie błędu lokalnego w metodach jednokrokowych wybór b 1 =0, b 2 =1, c=1/2, a =1/2 dawał RK2 punktu środkowego 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. Rozwijając do jednego rzędu wyżej z  t uzyskamy oszacowanie błędu lokalnego d n = u(t n ) - u n [przy założeniu, że u(t n-1 ) = u n-1 ] świetny wzór choć mało praktyczny

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

7 t n-1 t n t n+1 t n-1 t n+1 tt tt  t ekstrapolacja Richardsona dwa kroki  t: dostaniemy lepsze oszacowanie u(t n+1 ) jeden krok 2  t: dostaniemy gorsze oszacowanie u(t n+1 ) szacujemy C n z porównania obydwu rozwiązań

8 błąd lokalny u(t n )-u n =d n jest: wykonujemy krok następny od t n do t n+1 1)zakładamy, że krok jest na tyle mały, że stała błędu się nie zmienia C n  C n+1 (lub, ze w jednym kroku zmienia się o O(  t)] wtedy błąd lokalny popełniony w chwili t n+1 jest d n+1  d n. 2) gdy krok mały: współczynnik wzmocnienia błędu   1 [np. du/dt=- u jawny E. (1-  t)] (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  d n +d n+1  2d n, ekstrapolacja Richardsona odchylenie wyniku numerycznego od dokładnego u(t n+1 )-u n+1 =  d n +d n+1

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

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

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

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

13 Euler RK2 RK2 z odciętym błędem ekstrapolacja Richardsona

14 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

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

18 Przykład: oscylator harmoniczny używane oszacowanie błędu z RK2 Uwaga: tutaj rozwiązania nie poprawiamy przez ekstrapolacje tol=1e-2 tol=1e-3 start x2x2 V2V2 tt algorytm ustawia minimalny krok czasowy gdy zmiany prędkości lub położenia są maksymalne RK2 tolerancja błędu obcięcia tol=0.1

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

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

21 u(0)=0 proste równanie traktowane jawnym schematem Eulera

22 niech  >> 0 szybkozmienna składowa wolnozmienna  prosty problem nieco skomplikujemy

23 rozwiązanie  krok dt=0.02 jest bardzo drobny w porównaniu ze skalą zmienności składowej parabolicznej krzyżyki : t 2 kółka : t 2 + exp(-100t)

24 rozwiązanie dokładne dt=0.019 dt=0.02 dt=0.021  krok dt=0.02 okazuje się zbyt długi gdy włączyć składową szybkozmienną nawet tam, gdy znika ona z rozwiązania

25 rozwiązanie dokładne dt=0.019 dt=0.02 dt=0.021  część szybkozmienna gaśnie szybko, ale w schemacie jawnym Eulera nakłada ograniczenie na krok czasowy : u’=-  u  =100  dt

26  t Re( )  t Im ( ) 1  t Re( )  t Im ( ) -2 1 metoda Eulera jawna niejawna metoda Eulera regiony stabilności metod Eulera w metodzie niejawnej problemu ze stabilnością bezwzględna nie ma...

27 niejawna metoda Eulera: zastosowanie do problemu sztywnego dokładny dt=0.02 dt=0.04 rozwiązania są stabilne i dokładne dla dużych t nawet gdy dt duże dla małych t można wstawić mniejsze dt, potem krok zwiększyć

28 Problemy sztywne (drętwe) (stiff, stiffness) Problem jest praktyczny i ścisłej definicji, która byłaby użyteczna, nie ma. Jedna z możliwych: problem jest sztywny, gdy stosując schemat jawny musimy przyjąć krok czasowy bardzo mały w porównaniu ze skalą zmienności funkcji. RRZ jest problemem sztywnym gdy: 1.Problem jest charakteryzowany bardzo różnymi skalami czasowymi 2.Stabilność bzwz nakłada silniejsze ograniczenia na krok czasowy niż dokładność. 3.Metody jawne się nie sprawdzają. niech  >> 0 szybkozmienna składowa wolnozmienna

29 Problemy sztywne (drętwe) (stiff) Ogólna postać układu równań pierwszego rzędu wektor R n fcja R  R n  R Tylko niekiedy można podać rozwiązanie w zamkniętej formie analitycznej. Można, np. dla jednorodnego problemu liniowego problem najczęściej spotykany dla układ równań różniczkowych opisujących sprzężone procesy o bardzo różnych skalach czasowych

30 gdzie np. problem rozpadu promieniotwórczego Izotop 2 o stałej rozpadu 2 rozpada się promieniotwórczo na inny izotop 1 o stałej rozpadu 1 c j liczone z warunku początkowego y 1 (0)=0 y 2 (0)=1 dla niezdegenerowanych wartości własnych problemy sztywne wartości własne – 1, – 2 rozłożyć warunek początkowy na wektory własne

32 podobny do poprzedniego problem sztywny z liniowego równania drugiego rzędu o bliskich współczynnikach następny przykład: wartości / wektory własne: -1 / [-1,1] T -1000 / [1,-1000] T bardzo różne skale czasowe u’’+1001u’+1000u=0

33 szczególnie dotkliwy przypadek: równanie niejednorodne (bez rozwiązania analitycznego) Rozwiązanie będzie miało postać: stan przejściowy (wszystkie zgasną) stan ustalony wolnozmienny Na czym polega problem? : Rozwiązując problem numerycznie metodą jawną (Euler, RK2) musimy przyjąć krok czasowy  t < 2/| _max| aby uniknąć eksplozji rozwiązań nawet gdy wszystkie wyrazy z powyższej sumy w rozwiązaniu znikają problemy sztywne załóżmy, że wartości własne A są ujemne

34 y 2 (0)=1 y 1 (0)=0 1 =1/10 2 =1/10 000 bardzo wolno się rozpada [taka i większa rozpiętość lambd typowa również dla reakcji chemicznych spotykana również dla układów elektrycznych] y 2 – izotop matka wolno rozpadająca się na y 1 y 1 – izotop szybko rozpadający się, niejednorodność: dodatkowo pewna ilość jest w stałym tempie doprowadzana z zewnątrz przy zaniedbywalnej wielkości 2 y 1 =0.5 spełnia pierwsze równanie

35 tol=0.001 t t tt automatyczna kontrola kroku czasowego dla jawnego RK2 z krokiem czasowym ustawianym przez ekstrapolację Richardsona tol=0.00001 w obydwu przypadkach  t tylko chwilowo przekracza krytyczną wartość 2/(1/10)=20 zęby: zaakceptowane błędy y y t 1 =1/10 2 =1/10 000

36 RK4 2.78 / 1

37 Wzór trapezów stały krok, bardzo dłuuugi stały  t=200 nic złego się nie dzieje ze stabilnością w stanie „ustalonym” Zastosujmy metodę A-stabilną = wzór trapezów (p=2) t y

38 Wzór trapezów i krok automatycznie dobierany przez ekstrapolację Richardsona tol=0.01 kropki -tam gdzie postawiony krok Krok czasowy – zmienność 4 rzędów wielkości. raptem 10 kroków i załatwione! zamiast 10 4 kroków RK4 t t y y

39 trapezy (najdokładniejsza metoda A-stabilna spośród wielokrokowych) tt maksymalnie parę tysięcy metoda trapezów: jako A-stabilna radzi sobie nieźle z doborem kroku czasowego w problemach sztywnych – ale jest stosunkowo mało dokładna dokładniejsza A-stabilna pozwoliłaby stawiać jeszcze dłuższe kroki niestety = dokładniejszej A-stabilnej tej w klasie metod (liniowe wielokrokowe) nie ma dlatego : niejawne metody RK (jednokrokowe, nieliniowe) poziom jawnych RK t y t z tolerancją 0.00001

40 trapezy z tolerancją 0.00001 (najdokładniejsza metoda A-stabilna spośród wielokrokowych) niejawna dwustopniowa metoda RK (rzędu 4) z tolerancją 0.00001 (A-stabilna) tt tt maksymalnie parę tysięcy maksymalnie kilkadziesiąt tysięcy t t tt y y