Numerické metody I Jan Schee FPF SÚ Opava 2. října 2018 1 Strojová reprezentace čísel Ve smyslu numerických výpočtů rozlišujeme dva číselné typy: celočíselný a s plovoucí desetinnou čárkou. V jazyce C jsou celočíselné typy implementovány v datových typech unsigned int, int, unsigned char, char. Co se týče čísel s plovoucí desetinnou čárkou tak jsou to datové typy float a double. Aritmetika číslicové techniky je postavena dvojkové soustavě. Čísla jsou ukládaná do bitových polí, která jsou organizována po osmi do jednotky byte (B). 1.1 Formát celého čísla Celé číslo zapsané do dvojkové (binární) baze bude mít tvar a = n−1 k=0 ak2k (1) 1 kde koeficienty rozvoje a k jsou jednotlivé hodnoty bitového pole, viz. tabulka ukazující zápis čísla 136 do 8-bitového pole. a7 a6 a5 a4 a3 a2 a1 a0 1 0 0 0 1 0 0 0 Uvažujme například datový typ unsigned char, ketrý je tvořen polem 8 bitů, tj. n = 8. Mějme binární číslo (10001000)2 které převedeme na číslo v desítkové soustavě a obdržíme (10001000)2 = (1 × 27 + 0 × 26 + 0 × 25 + 0 × 24 + 1 × 23 +0 × 22 + 0 × 21 + 0 × 20 )10 = (128 + 8)10 = (136)10 (2) V případě znaménkových typů je nejvyšší bit vyhrazen hodnotě znaménka a to tak, že sign(a) = (−1)a7 (3) (v případě datového typu char). 1.2 Formát čísla s plovoucí desetinou čárkou Norma IEEE 754 definuje single-precision binary floating-point format: binary32. Jeho 32 bitů je rozděleno do tří polí: • Znaménkový typ s: 1 bit • Exponent e: 8 bitů 2 • Mantisa m: 23 bitů Následující tabulka ilustruje upspořádání jednotlivých polí v datovém typu float: s e7 e6 e5 e4 e3 e2 e1 e0 m22 m21 · · · m0 Hodnotu takto zakódovaného desetinného čísla získáme ze vztahu r = (−1s ) × 1 + 23 k=1 m23−k2−k × 2e−127 . (4) Jako příklad uveďme číslo 5.125 zakódované podle normy IEEE 754. Bitové pole je zaplněno tak jak ukazuje tabulka 0 10000001 01001000000000000000000 Odkud je zřejmě exponent e = (10000001)2 = (129)10 a mantisa je m = (010 0100 0000 0000 0000 0000)2 = (2359296)10. Po dosazení do předchozího vztahu obdržíme hodnotu r = 1 × (1 + 1 × 2−2 + 1 × 2−5 ) × 2129−127 = (1 + 0.25 + 0.03125) × 22 = 1.28125 × 4 = 5.125. (5) 3 Cvičení 1. Vytvořte program, který má na vstupu celé číslo v desístkové soustavě a převeďte jej do dvojkové soustavy a zobrazte na obrazovce (hint: použijte datový typ unsigned int, výsledek uložte do pole char a[32]). 2. Napište program, který zobrazí bitové pole čísla uloženého v proměnnné o datovém typu float a určete exponent a mantisu. Výsledky zobrazte na obrazovce. 2 Typy chyb a jejich šíření Díky konečné velikosti datových typů kódujících čísla s plovoucí detetinou čárkou není většina čísel strojově reprezento- vatelná. Například číslo (0.1)10 není strojově reprezovatelné, protože má nekonečný rozklad do dvojkové baze, jak ukazuje následující roz- klad 0.110 = 0.00011001100110011001100110011001100110011001 10011001100110011001100110011001100110011001100 110011001 · · ·2 . (6) Nechť je množina strojově reprezovatelných čísel A. Pokud pro číslo x platí, že x /∈ A pak hledáme takové číslo g ∈ A pro které platí |x − g| ≤ |x − b| , ∀b ∈ A, (7) 4 tj., hledáme nejbližší strojově reprezentovatelné číslo g k číslu x. Základním mechanizmem, kterým lze toho dosáhnou je zao- krouhlování. Uvažujme číslo a < 1 a definujme následující reprezentaci a = 0.a1a2a3 · · · aiai+1 · · · , 0 ≤ ai ≤ 9 ∧ a1 = 0. (8) Z čísla a vytoříme zaokrouhlené číslo na t -platných číslic takto a = 0.a1a2a3 · · · at pro at+1 ≤ 4 0.a1a2a3 · · · (at + 1) pro at+1 ≥ 5 (9) Zaokrouhlené číslo potom píšeme ve tvaru rd(x) = sign(x) × a × 10b . (10) Pomocí této definice snadno určíme relativní chybu zaokrouhlení, tj. r = rd(x) − x x . (11) Nejprve určeme, čemu se rovná ∆ ≡ rd(x) − x: • at+1 ≤ 4 ∆ = 0.a1a2 · · · at × 10b − 0.a1a2 · · · atat+1 · · · × 10b = at+1 · at+2 · · · × 10b−(t+1) ≤ 5 × 10b−(t+1) , (12) protože je at+1 ≤ 4 a at+2 ∈ [0, 1, 2, . . . , 9]. 5 • at+1 ≥ 5 ∆ = 0.a1a2 · · · (at + 1) × 10b − 0.a1a2 · · · atat+1 · · · × 10b = bt+1 ˙bt+2 · · · × 10b−(t+1) ≤ 5 × 10b−(t+1) (13) protože zřejme platí bt+1 = 10 − at+1 ≤ 5 a současně at+1 ≥ 5. Pro relativní chybu zaokrouhlování tak dostaneme odhad r ≤ 5 × 10b−(t+1) a × 10b = 5 a × 10−t+1 ≤ 5 × 10−t (14) což je důsledkem faktu, že a ≥ 10−1 . Zavádíme označení EPS = 5×10−t a nazýváme tuto veličinu strojová přesnost. Chyba, které se při zaokrouhlování dopouštíme je tzv. zaokrouhlovací chyba (trucation error). Tato chyba se netýká pouze vstupních dat ale mohou a jsou touto chybou zatíženy i aritmetické operace , které nad čísly s plovoucí desetinou čárkou provádíme. Může se totiž stát, že pro x, y ∈ A bude x + y /∈ A. Je tedy důležité vědět, jak se při danám algoritmu (který je tvořen aritmetickými operacemi) zaokrouhlovací chyba šíří. Nechť φ je funkce reprezentující algoritmus, tj. φ : D → Rn , φ(x) =    φ1(x1, . . . , xn) ... φm(x1, . . . , xn)    (15) dále nechť je ˜x aproximace čísla x a označme ∆xi ≡ ˜xi − xi a ∆x = ˜x − x absolutní chybu a ˜xi ≡ ∆xi/xi relativní chybu 6 aproximace. Dosadíme-li na vstupu ˜x za x tak dostaneme výsledek ˜y = φ(˜x), který se liší od y = φ(x). Absolutní chyba výsledku bude (Taylorův rozvoj do prvího řádu v ∆xi) ∆yi ≡ ˜yi − yi = φi(˜x) − φi(x) = n j=1 ∆xj ∂φi ∂xj . (16) Veličina ∂φi/∂xj reprezentuje "citlivost algoritmu"na změnu ∆x. Analogicky, pro relativní chybu výsledku dostaneme vztah yi = n j=1 xj φi(x) ∂φi(x) ∂xj xj . (17) Dokažme poslední vztah. Relativní chybu definujeme vztahem yi ≡ ∆yi yi . (18) Dosazením za ∆yi obdržíme vztah yi = φi(˜x) − φi(x) φi(x) n j=1 ∆xj φi(x) ∂φi(x) ∂xj = n j=1 ∆xj φi(x) ∂φi(x) ∂xj xj xj = n j=1 xj φi(x) ∂φi(x) ∂xj xj (19) C.B.D. 7 Ilustrujme šíření zaokrouhlovací chyby na příkladu y = φ(a, b, c) = a + b + c. Pro relativní chybu výsledku obdržíme výsledek y = 3 j=1 xj φ(x) ∂φ(x) ∂xj xj = a a + b + c a+ b a + b + c b+ c a + b + c c. (20) Z výsledku vyčteme, tento problém je dobře definován pokud je součet a + b + c mnohem větší než jednotlivé hodnoty a, b, c. Podívejme se na další příklad. Nechť je tentokráte y = −p + p2 + q. Pro relativní chybu výpočtu tohoto vzorce dostaneme výraz y = − p p2 + q p + p + p2 + q 2 p2 + q q. (21) Pro q > 0 jsou oba faktory < 1. Tento algoritmus je špatně definován pro q = −p2 . Nakonec určeme relativní chyby aritmetických operací ±, ∗ a /. Obdržíme následující výsledky φ(x, y) ≡ x ∗ y → xy = x + y, (22) φ(x, y) ≡ x/y → xy = x − y, (23) φ(x, y) ≡ x ± y → xy = x x ± y x ± y x ± y y. (24) Vydíme, že násobení a dělení nejsou nebezpečné operace. Pozor si musíme dávat při odečítání, kdy dochází k amplifikaci chyby pro x y. 8 3 Hledání kořene nelineární 1D-rovnice f(x) = 0 Pokud je funkce f(x) spojitá a hladká, tak podmínkou pro existenci kořene r na intervalu [x1, x2] je • f(x1)f(x2) < 0, • ∀x ∈ [x1, x2] bude buď f (x) < 0 nebo f (x) > 0. Pokud hledáme všechny kořeny rovnice f(x) = 0 tak musíme provést na intervalu [a, b] tzv. “bracketing”, tj. musíme najít subintervaly na kterých funkce f(x) splnuje výše uvedené pod- mínky. Bracketing Nejednoduší způsob jak zajistit aby na pod-intervalech byla funkce f(x) monotónní, je rozdělit definiční obor na co nejmenší části. Následuje pseudokód ilustrující myšlenku bracketingu. 1 (∗ bracketing ∗) 2 h=(b−a )/N 3 indx=0 4 for i =0 to N step 1 5 x1=a+i ∗h ; 6 x2=x1+h ; 7 i f f ( x1 )∗ f ( x2)<0 then 8 xl [ indx ]=x1 9 xu [ indx ]=x2 10 indx=indx+1 11 end 12 end 9 Po ukončení tohoto algoritmu budou sub-intervaly na kterých existuje kořen rovnice f(x) = 0, Ji = [xl[i], xu[i]]. Jejich počet bude uložen v proměnné indx. Metoda bisekce Půlením intervalu [x1, x2] postupně udoláme kořen. Jestliže po n - tém kroku je šířka subintervalu n pak v následujícím bude zřejmě poloviční, tj. n+1 = n 2 . (25) Počáteční interval má šířku 0. Počet kroků n potřebných k dosažení šířky subintervalu n potom je n = log2 0 n . (26) Dokažme toto tvrzení. Pro tři po sobě následující kroky bisekce postupně dostáváme • 1 = 0/2 • 2 = 1/2 = 0/4 • 3 = 2/2 = 0/8 • . . . Z této posloupnosti je vidět, že pro obecné n platí 0 n = 2n ⇒ n = log2 0 n . (27) 10 Když budeme chtít najít kořen s přesností odpovidájící n = 10−12 v okolí x = 1 bude potřeba n 40 kroků. Následující pseudokód ilustruje implementaci metody bisekce při hledání kořene 1D rovnice f(x) = 0. 1 (∗ Bisekce ∗) 2 x1=a 3 x2=b 4 m=(x1+x2 )/2 5 while abs ( x1−x2)>eps and abs ( f (m))> eps do 6 i f f (m)∗ f ( x2)<0 then 7 x1=m 8 else 9 x2=m 10 end 11 m=(x1+x2 )/2 12 end Na konci algoritmu bisekce bude přibližný kořen rovnice (určený s přesností eps) uložen v proměnné m. Metoda sečen Tato metoda je založena na lineární aproximací fuknce f(x) na daném intervalu. Označme interval monotónosti I ≡ [a, b] . Přímka aproximující f(x) na daném intervalu bude předepsána formulí y(x) = fb − fa b − a (x − a) + fa (28) kde je fa ≡ f(a) a fb ≡ f(b). Odhad kořene rovnice f(x) = 0 potom plyne jako řešení rce y(x) = 0, což nám dá výsledek xc = b − fa b − a fb − fa . (29) 11 Pokud je f(xc) < tak jsme našli kořen rovnice f(x) = 0. Pokud zmíněná podmínka neplatí, pak je interval [a, b] rozdělen na dva podintervaly [a, xc] a [xc, b] a skutečný kořen leží v jednom z nich. Pokud leží v [a, xc] tak položíme b = xc v opačném případě položíme a = xc . Celý postup opakujeme dokud není dosaženo požadované přesnosti . Metodu sečen shrnuje následující pseu- dokód: 1 a i=a 2 bi=b 3 f c=f ( ( a+b )/2) 4 while abs ( ai−bi)>eps do 5 fa=f ( a i ) 6 fb=f ( bi ) 7 xc=b−fa ∗(b−a )/( fb−fa ) 8 f c=f ( xc ) 9 i f abs ( f c )eps and ( f ( c )!=0 and f ( x2 )!=0) do 6 i f f0 != f1 and f2 != f1 then 7 c=f1 ∗ f2 ∗x0 /(( f0−f1 )∗( f0−f2 ) ) 13 8 c=c+f0 ∗ f2 ∗x1 /(( f1−f0 )∗( f1−f2 ) ) 9 c=c+f0 ∗ f1 ∗x2 /(( f2−f0 )∗( f2−f1 ) ) 10 else 11 c=x2−f1 ∗( x2−x0 )/( f2−f0 ) 12 end 13 f c=f ( c ) 14 i f x1 4 tak neznáme explicitní vyjádření pro kořeny polynomiálních rovnic Pn(x) = 0. Použijeme účinnou numerickou metodu založenou na Laugerre algoritmu. Nechť jsou x1, x2, . . . , xn kořeny rovnice Pn(x) = 0. Pak lze Pn(x) pomocí těchto kořenů zapsat vztahem Pn(x) = (x − x1)(x − x2) · · · (x − xn). (33) Nyní na poslední rovnici aplikujeme přirozený logaritmus a ob- držíme ln(Pn(x)) = ln |x − x1| + ln|x − x2| + · · · + ln |x − xn|. (34) Nakonec, pro poslední rovici určíme první a druhou derivaci d ln Pn(x) dx = 1 x − x1 + 1 x − x2 +· · ·+ 1 x − xn = Pn(x) Pn(x) ≡ G (35) a − d2 ln Pn(x) dx2 = 1 (x − x1)2 + 1 (x − x2 2) + · · · + 1 (x − xn)2 = Pn(x) Pn(x) 2 − Pn (x) Pn(x) . (36) 15 Hledaný kořen x1 nechť je vzdálen od aktuálního odhadu o délku a a ostatní kořeny o délku b, tj. x1 − x = a a xi − x = b, ∀ i = 2, 3, . . . , n. (37) Po dosazení do (35) a (36) obdržíme rovnice 1 a + n − 1 b = G, (38) 1 a2 + n − 1 b2 = H. (39) Po několika algebraických úpravách posledních dvou rovnic najdeme pro délku a výraz a = n G ± (n − 1)(nH − G2) . (40) Kořen x1 teď snadno najdeme iterativně. Pro odhad kořene x určíme, ze vztahu (40), parametr a a získáme nový odhad x := x + a (operátor := zde reprezetuje operaci přiřazení!). Takto pokračujeme dokud není |a| < nebo |f(x)| < . Pro určení zbývajících kořenů nejprve zredukujeme původní polynom Pn(x) = (x − x1)Qn−1(x) ⇒ Qn−1(x) = Pn(x)/(x − x1) (41) a aplikujeme celou předchozí iteraci na polynom stupně n − 1, Qn−1(x), atd. dokud neurčíme všechny kořeny. Redukce polynomu n -tého stupně na polynom stupně n − 1 si zaslouží detailní diskuzi. Nechť je z komplexní, pak jsou-li koeficienty ai polynomu Pn(z) reálná čísla pak komplexním sdružením tohoto polynomu dostaneme [Pn(z)]∗ = Pn(z∗). Dokažme toto tvrzení. 16 Počítejme [Pn(z)]∗ = (anzn )∗ + (an−1zn−1 )∗ + · · · + (a1z)∗ + a0 = an(zn )∗ + an−1(zn−1 )∗ + · · · + a1z∗ + a0 = an(z∗ )n + an−1(z∗ )n−1 + · · · + a1z∗ + a0 = Pn(z∗ ).(42) Zde jsme využili následující identity. Nechť v = a + ib a w = c + id jsou komplexní, pak plati: • (v + w)∗ = (a + ib + c + id)∗ = (a + c + i(b + d))∗ = a + c − i(b + d) = v∗ + w∗ (43) • (wn )∗ = (|w|einφ )∗ = |w|n (einφ )∗ = |w|n e−inφ = (w∗ )n (44) C.B.D. Nyní, pokud je z kořenem rovnice Pn(z) = 0 pak je kořenem i příslušné komplexně sdružené číslo z∗ ,tj. platí Pn(z∗) = 0. To znamená, že v algoritmu pro redukci polynomu vystačíme s reálnými poli a[] a v[] neboť zřejmě platí Pn(z) = (z −z1)(z −z∗ 1 )Qn−2(z) = (z2 +2az +a2 +b2 )Qn−2(z), (45) kde a a b určují kořen z1 = a+ib a zjevně jsou a a b reálná čísla. Polynom n − 2 stupně, Qn−2(z), má reálné koeficienty. Dále si, na základě této diskuze, musíme uvědomit, že kořeny rovnice 17 Pn(z) = 0 jsou bud reálná čísla, nebo komplexní. A pokud je nalezený kořen komplexní tak řešením této rovnice je určitě i kořen komplexně sdružený. Následující algoritmy ilustrují realizace Laguerrovy iterace, dělení polynomu a vyčíslení polynomu a jeho prvních dvou de- rivací. 1 (∗ Laguerrova i t e r a c e ∗) 2 l a g u e r r e ( c , n , x , eps ) 3 imax =100 4 z=0+I ∗0 5 while ( i<=imax ) 6 poly ( z , c , n , p) 7 i f abs (p [0]) < eps break 8 F=(n−1) 9 F=F∗(( n−1)∗p [ 1 ] ∗ p[1] −n∗p [ 0 ] ∗ p [ 2 ] ) 10 G1=p [1]+ s q r t (F) 11 G2=p[1] − s q r t (F) 12 i f ( abs (G1)>abs (G2) ) then G=G1 13 else G=G2 14 z=z−n∗p [ 0 ] /G 15 end while 16 end Redukci polynomu n- tého stupně na n − 1 stupeň provede algoritmus: 1 (∗ Deleni polynomu ∗) 2 poldiv ( a , n , z1 , v ) 3 p=a [ n ] 4 for i=n−1 to 0 step −1 5 v [ i ]=p 6 p=p∗z1+a [ i ] 7 end for 8 end Dále je potřeba vyčíslit polynom Pn v z a jeho první a druhou derivaci. To zajistí algoritmus: 18 1 (∗ Vycisleni P, P’ a P’ ’ ∗) 2 poly ( z , a , n , p) 3 //polynom 4 p [0]= a [ n ] 5 for i=n−1 to 0 step −1 6 p [0]= p [ 0 ] ∗ z 7 p [0]= p [0]+ a [ i ] 8 end for 9 // prvni d e r i v a c e polynomu 10 p [1]= n∗a [ n ] 11 for i=n−1 to 1 step −1 12 p [1]= p [ 1 ] ∗ z 13 p [1]= p [1]+ i ∗a [ i ] 14 end for 15 // druha d e r i v a c e polynomu 16 p [2]= n∗(n−1)∗a [ n ] 17 for i=n−1 to 2 step −1 18 p [2]= p [ 2 ] ∗ z 19 p [2]= p [2]+ i ∗( i −1)∗a [ i ] 20 end for 21 end 5 Lagrangova interpolace Mějme N bodů (x1, y1), (x2, y2), ..., (xN , yN ), kde je yi = f(xi) pro i = 1, 2, ..., N. Těmito body vede jedinečný polynom N − 1 stupně určený Lagrangeovou formulí PN−1(x) N i=1 yi N−1 j=1,j=i x − xj xi − xj . (46) Například pro polynom 2. stupně dostaneme výraz P2(x) = (x − x2)(x − x3) (x1 − x2)(x1 − x3) y1 + (x − x1)(x − x3) (x2 − x1)(x2 − x3) y2 19 + (x − x1)(x − x2) (x3 − x1)(x3 − x2) y3. (47) Pro určení chyby, které se dopouštíme při Lagrangově interpolaci definujeme funkci F(z) F(z) ≡ f(z) − y(z) − [f(x) − y(x)] pn(z) pn(x) (48) kde je pn(x) = n i=1 (x − xi) (49) Funkce F(z) má n + 1 uzlových bodů. Když n-krát zderivujeme funkci F(z) podle z dostaneme výraz F(n) (z) = f(n) (z) − y(n) (z) − f(x) − y(x) pn(x) n!. (50) Tato fukce má na intervalu, které určují nulové body x1, x2,. . . a bod x aspoň jeden nulový bod (n-krát jsme derivovali funkci, která jich má n+1). Označme tento nulový bod z = ξ a dosadíme jej do vztahu pro F(n) (z) a obdržíme rovnici 0 = F(n) (ξ) = f(n) (ξ) − y(n) (ξ) − f(x) − y(x) pn(x) n! (51) Všimněme si, že f(n) (z) = 0, protože f je polynom n−1 stupně. Označíme-li chybu interpolace E(x) = f(x) − y(x) tak z předchozí rovnice dostaneme výsledek E(x) = pn(x) n! f(n) (ξ). (52) 20 K vyčíslení Lagrangova polynomu se využívá Nevillův algoritmus. Pošme jej pomocí příkladu pro N = 4. Strom Nevillova pak vypadá následovně: x1: y1 = P1 P12 x2: y2 = P2 P123 P23 P1234 x3: y3 = P3 P234 P34 x4: y4 = P4 kde jsou P1, . . . , P4 polynomy 0-tého stupně odpovídající funkčním hodnotám f(xi). P12 je polynom 1. řádu procházející body (x1, y1), (x2, y2) a analogicky ostatní Pij. P123 je polynom druhého stupně procházející body (x1, y1), (x2, y2) a (x3, y3). Nakonec P1234 je požadovaný výstup, polynom třetího stupně procházející všemi čtyřmi body. Vztah mezi rodičovskými a dceřinými polynomy vyjadřuje výraz Pi(i+1)···(i+m) = 1 xi − xi+m (x − xi+m)Pi(i+1)···(i+m−1)+ (xi − x)P(i+1)(i+2)···(i+m) . (53) Následující pseudokód implementuje Nevillův algoritmus: 1 lagrange ( xx , x [ ] , y [ ] , n) 2 for m =1 to n 3 A[m ]=y [m ] 4 end for 5 for m =1 to n−1 21 6 for j=1 to n−m 7 B[ j ]=(( xx−x [ j+m ] ) ∗A[ j ] 8 −(xx−x [ j ] ) ∗A[ j +1])/( x [ j ]−x [ j+m ] ) 9 end for 10 for j=1 to n−m 11 A[ j ]=B[ j ] 12 end for 13 end for 14 return B [ 1 ] 15 end 6 Spline interpolace Máme N + 1 bodů (x0, y0), . . . , (xN , yN ), kde je yi = f(xi). Na intervalech [xi, xi+1] budeme aproximovat funkci f(x) lineární funkcí gi(x) = aix + bi. (54) Mluvíme pak o lineárním spline. Podmínkou pro konstrukci lineárního spline je aby byl spojitý, tj. aby platilo gi(xi+1) = gi+1(xi+1) (55) kde gi(x) je lineární funkce na intervalu [xi, xi+1]. Lineární spline má potom tvar gi(x) = (yi+1 − yi) (xi+1 − xi) (x − xi) + yi, ∀ i = 0, . . . , N − 1. (56) Snadno ověříme, že platí výše zavedená podmínka spojitosti. Platí totiž gi(xi+1) = yi+1 − yi + yi = yi+1 = gi+1(xi+1). (57) 22 Lineární spline g(x) zapisujeme ve tvaru g(x) =    g0(x) pro x ∈ [x0, x1] g1(x) pro x ∈ [x1, x2] ... gN−1(x) pro x ∈ [xN−1, xN ] (58) Chceme-li aby se lineární spline g(x) co nejvíce přiblížil funkci f(x), pak musíme mít k dispozici co nejvíce bodů (xi, yi). Lepší aproximace dosáhneme volbou kubické funkce na intervalech [xi, xi+1]. V případě kubického spline nahradíme lineární funkce (54) kubickou funkcí Si(x) = aix3 + bix2 + cix + di. (59) Každá funkce Sije určena čtyřmi parametry (ai, bi, ci, di). To je celkem 4N neznámých parametrů. Následující podmínky, jejichž splnění u kubického spline uplatňujeme nám dá potřebných 4N rovnic. 1. Požadujeme, aby náš spline interpoloval všech N +1 bodů, dostáváme tak podmínky Si(xi) = yi pro 0 ≤ i ≤ N − 1 a SN−1(xN ) = yN . (60) Tím jsme získali n + 1 podmínek. 2. Požadujeme spojitost slpline, jeho prví a druhé derivace. To nám dává následující podmínky Si−1(xi) = Si(xi) pro 1 ≤ i ≤ N − 1 dává N − 1 podmínek 23 Si−1(xi) = Si(xi) pro 1 ≤ i ≤ N − 1 dává N − 1 podmínek Si−1(xi) = Si (xi) pro 1 ≤ i ≤ N − 1 dává N − 1 podmínek Máme celkem (N + 1) + 3(N − 1) = 4N − 2 podmínek. Chybí nám ještě dvě. V případě tzv. přirozeného spline se požaduje aby byly navíc splněny podmínky S (x0) = S (xN ) = 0 (61) tj. v krajních bodech interpolované funkce požadujeme nulovost druhých derivací kubického spline. Uplatněním výše uvedených podmínek určíme kubický spline následovně. Nechť pro všechna 0 ≤ i ≤ N platí zi ≡ Si (xi) (62) a dále nechť je hi ≡ xi+1 − xi. (63) Protože je Si(x) kubickým polynomem an intervalu [xi, xi+1], tak Si (x) je na tomto intervalu lineární, tj. Si (x) = zi + mi (x − xi) (64) kde je mi = zi+1 − zi hi (65) Rovnici (64) snadno převedeme na tvar Si (x) = zi+1 hi (x − xi) + zi hi (xi+1 − x) . (66) 24 Ze znalosti zijsme schopni určit druhou derivaci kubického spline. Nás ale zajímá funkce Si(x) takže budeme funkci (66) dvakrát integrovat a dostaneme Si(x) = zi+1 6hi (x − xi) 3 + zi 6hi (xi+1 − x) 3 +Ci (x − xi)+Di (xi+1 − x) (67) kde jsou Cia Di integrační konstanty, které určíme z podmínek Si(xi) = yi a Si(xi+1) = yi+1. (68) Takto obdržíme pro Cia Di vztahy Ci = yi+1 hi − hi 6 zi+1, Di = yi hi − hi 6 zi. Ke konečnému výpočtu funkce (67) musíme určit samotné zi. K jejich určení použijeme poslední podmínku spojisti prvních derivací kubického spline, tj. Si−1(xi) = Si(xi). (69) Derivací rovnice (67) dostaneme následující dvě rovnice Si−1(xi) = hi−1 6 zi−1 + hi−1 3 zi + bi−1, Si(xi) = − hi 6 zi+1 − hi 3 zi + bi, kde je bi ≡ (yi+1 −yi)/hi. Dosazením posledních dvou rovnic do (69) dostáváme rovnici pro 1 ≤ i ≤ N − 1 hi−1zi−1 + 2 (hi−1 + hi) zi + hizi+1 = 6 (bi − bi−1) . (70) 25 Nezapomeňme, že pro přirozený spline platí z0 = zn = 0. Dostali jsme tak tridiagonální lineární systém          u1 h1 h1 u2 h2 h2 u3 h3 ... ... ... hN−3 uN−2 hN−2 hN−2 uN−1          ·          z1 z2 z3 ... zN−2 zN−1          =          p1 p2 p3 ... pN−2 pN−1          . (71) ui = 2(hi + hi−1) (72) a pi = 6(bi − bi−1). (73) Tridiagonální systém K vyřešení tridiagonálního lineárního systému použijeme následující algoritmus. Hledáme řešení systému, který má tvar       b0 c0 0 . . . a1 b1 c1 . . . . . . . . . aN−2 bN−2 cN−2 . . . 0 aN−1 bN−1       ·       u0 u1 · · · uN−2 uN−1       =       r0 r1 · · · rN−2 rN−1       (74) Řešení spočíva v postupné eliminaci prvního nenulového prvku v řádcích 1 až N −1. Tímto určíme hodnotu N-té složky vektoru 26 u, tj. uN−1. Postupnou zpětnou substitucí od řádku č.N − 2 po řádek 0 získáme zbylé složky vektoru u. V maticové notaci bude soustava vypadat následovně M · u = r. (75) Budeme ekvavalentními operaceni upravovat matici (M|r) (76) nebo ve složkovém zápise      M00 M01 M02 . . . M0N−1 M10 M11 M12 . . . M1N−1 ... MN−10 MN−11 MN−12 . . . MN−1N−1 M0N = r0 M1N = r1 ... MN−1N = rN−1      . (77) Protože máme tridiagonální systém, tak stačí eliminovat na itém řádku prvek Mii−1. Provedeme to tak, že i-tý řádek vynásobíme zlomkem −Mi−1i−1/Mii−1 a přičteme řádek i − 1, tzn. pro j-tý sloupec i-tého řádku provedeme následující výpočet ˜Mij = − Mi−1i−1 Mii−1 Mij + Mi−1j. (78) Tento vztah je aplikován na řádky matice Mij pro i = 1, . . . , N− 1. Novou matici soustavy ˜M budeme opět značit M. Řešení soustavy (77) najdeme zpětnou substitucí, tj. pro i = N − 1 bude uN−1 = MN−1N /MN−1N−1. (79) 27 Zbývající složky vektoru u vycházejí ze vztahu ui = 1 Mii  MiN − N−1 j=i−1 ujMij   (80) Tento postup je detailně ukázách na následujícím pseudokódu: 1 (∗ transformace t r i d i a g systemu ∗) 2 for i =1 to N−1 step 1 3 MM=M[ i −1][ i −1]/M[ i ] [ i −1]; 4 for j=0 to N step 1 5 M[ i ] [ j ]=−MM∗M[ i ] [ j ]+M[ i −1][ j ] 6 end 7 end 8 (∗ vypocet korenu u_i ∗) 9 u [N−1]=M[N−1][N] /M[N−1][N−1] 10 for i=N−2 to 0 step −1 11 u [ i ]=M[ i ] [ N] 12 for j=i −1 to N−1 step 1 13 u [ i ]=u [ i ]−u [ j ] ∗M[ i ] [ j ] 14 end 15 end Spline je definován na intervalech [xi, xi+1] a je tedy potřeba, k určení jeho hodnoty v bodě x, znát do kterého intervalu x náleží. Pokud jsou jsou body xi ekvidistantní se šířkou intervalu h, pak index i, který identifikuje levý, krajní bod intervalu do kterého x náleží, určíme takto x = x0 + i ∗ h ⇒ i = Chop x − x0 h . (81) Funkce Chop odtrhává od reálného čísla jeho desetinnou část, takže nám zůstane celé číslo. 28 Shrnutí Interval [x0, xN ] rozdělíme na N+1 částí. Na každé z nich určíme funkci Si(x) = aix3 + bix2 + cix + di (82) a požadujeme aby tyto funkce na jednotlivých intervalech splňovali podmínky Si(xi+1) = Si+1(xi+1) (83) a Si (xi+1) = Si+1(xi+1). (84) Postupnými algebraickými úpravami dostaneme Si(x) = zi+1 6hi (x−xi)3 + zi 6hi (xi+1−x)3 +Ci(x−xi)+Di(xi+1−x) (85) kde je Ci = yi+1 hi − hi 6 zi+1, (86) Di = yi hi − hi 6 zi, (87) bi = yi+1 − yi hi , (88) yi = Si(xi), (89) zi = Si (xi), (90) hi = xi+1 − xi. (91) K určení splinu je tedy nutné určit pouze zi. Hodnoty zi jsou řešením tridiagonálního systému (71). 29 7 Numerická integrace Neprve se seznámíme s integrací s fixním krokem kdy budeme aproximovat určitý integrál I = b a f(x)dx (92) pomocí kvadratury, tj. sumou I = h N i−0 αifi (93) kde je N dělení intervalu [a, b], h ≡ (b − a)/N, fi ≡ f(xi) a xi = a + i · h pro i = 0, · · · , N. Koeficienty αi jsou váhy součtu. Aproximace pomocí kvadratur je založena na myšlence, že na funkci f(x) na intervalu [a, b] nahradíme Legengrovým interpolačním polynomem, který jsme studovali výše. Napišme jej ve tvaru PN (x) ≡ N i=0 fiLi(x) (94) kde zřejmě je Li(x) = N k=0,k=i x − xk xi − xk . (95) Integrál I pak aproximujeme takto b a f(x)dx b a PN (x)dx = N i=0 fi b a Li(x)dx. (96) 30 Když dále zavedeme substituci x = a + t · h (97) bude mít funkce Li(x) tvar Li(x) = φi(t) = N k=0,k=i a + th − a − kh a + ih − a − kh = N k=0,k=i t − k i − k (98) a diferenciál dx = hdt. Výsledný vzorec pro aproximaci integrálu I kvadraturou bude mít tvar I h N i=0 fi N 0 φi(t)dt = h N i=0 fiαi. (99) Koeficienty αi jsou potom určeny vztahem αi ≡ N 0 φi(t)dt. (100) Pro názornost odvodíme lichoběžníkové a Simpsnovo pravidlo. Lichoběžníkové pravidlo K aproximaci použijeme pouze dva krajní body intervalu [a, b], tzn. N = 1. Určíme koeficienty alpha0 a α1 : α0 = 1 0 t − 1 −1 dt = − t2 2 − t 1 0 = −(1/2 − 1) = 1/2, (101) 31 a α1 = 1 0 t − 0 1 − 0 dt t2 2 1 0 = 1/2. (102) Výsledkem bude integrační pravidlo I h 2 (fa + fb). (103) Simpsonovo pravidlo V tomto případě aproximace se použijí tři body intervalu [a, b], oba krajní a jeden prostřední, tzn. N = 2. Opět určíme koeficienty αi , tentokráte pro i = 0, 1 a 2. Dostaneme α0 = 2 0 (t − 1)(t − 2) (0 − 1)(0 − 2) dt = 1 2 2 0 (t2 − 3t + 2)dt = 1 2 t3 3 − 3 t2 3 + 2t 2 0 = 1 3 (104) α1 = 2 0 (t − 0)(t − 2) (1 − 0)(1 − 2) dt = − 2 0 (t2 − 2t)dt = − t3 3 − t2 2 0 = 4 3 (105) α2 = 2 0 (t − 0)(t − 1) (2 − 0)(2 − 1) dt = 1 2 2 0 (t2 − t)dt = 1 2 t3 3 − t2 2 2 0 = 1 3 . (106) Výsledné integrační pravidlo pak píšeme ve tvaru I h 3 (f0 + 4f1 + f2). (107) Vytvořme pseudokód pro implementaci Simpsonova pravidla. Mějme reálnou funkci f(x), na intervalu [a, b], který rozdělíme na M podintervalů, přičemž M je sudé číslo. Pro ilustraci se podívejte na dělení intervalu [a, b] pro M = 8 na obrázku. 32 Hodnoty Si představují dílčí aproximace integrálu nad příslušným podintervalem. V následujícím pseudókdu je ukázána impementace Simpsonova pravidla. 1 simpson ( f , a , b ,M) 2 h=(b−a )/M 3 S=0 4 for i =0 to M/2 5 x0=a+2∗ i ∗h 6 x1=a+(2∗ i +1)∗h 7 x2=a+(2∗ i +2)∗h 8 f0=f ( x0 ) 9 f1=f ( x1 ) 10 f2=f ( x2 ) 11 S = S + h∗( f0 +4∗ f1+f2 )/3 12 end for 13 return S 14 end Gauss-Legendrova kvadratura V předchozích podkapitolách byl určitý integrál funkce aproximován součtem funkčních hodnot na množině ekvidistantních bodů xi, vynásobené vhodně zvoleným váhovým parametrem. Gaussova kvadratura, nám umožˇuje vybrat jak váhové parametry tak souřadnice (abscissa). Vhodnou volbou vah a souřadnic lze zajistit, to aby pro integrandy typu polynom × W(x) 33 (W(x) je známá fukce) byl určitý integrál exaktní. Funkci W(x) můžeme zvolit tak, abychom odstraili integrovatelné singularity z požadovaného ntegrálu. Pro danou funkci W(x) a dané N najdeme množinu vah wj a souřadnic xj tak, že aproximace b a W(x)f(x)dx N j=1 wjf(xj) (108) je exaktní pokud je f(x) polynom. Určení vah a souřadnic je postaveno na existenci ortogonálních polynomů. Polynomy f(x) a g(x) jsou ortogonální, pokud je jejich skalární součin < f|g >≡ b a W(x)f(x)g(x)dx (109) roven nule. Funkce f(x) je normalizovaná pokud platí < f|f >= 1. (110) Množina polynomů splňujících obě podmínky je ortonormální možina polynomů. Lze najít množiny polynomů které splňují: • zahrnují přesně jeden polynom j-tého řádu Pj(x) ∀ j = 0, 1, 2, . . .. • všechny jsou vzájemně ortogonální přes danou váhovou funkci W(x). 34 Procedura pro nalezení takovéto množiny polynomů je určena rekurentním vztahem P−1(x) ≡ 0, (111) P0(x) ≡ 1, (112) Pj+1(x) = (x − aj)xPj(x) − bjPj−1(x) (113) kde je aj ≡ < xPj|Pj > < Pj|Pj > pro j = 0, 1, . . . (114) bj ≡ < Pj|Pj > < Pj−1|Pj−1 > pro j = 1, 2, . . . (115) Polynomy Pj(x) mají přesně j kořenů na intervalu (a, b). Souřadnice(abscisy) N-bodové Gaussovy kvadratury s váhovou funkcí W(x) na intervalu (a, b) jsou kořeny ortogoválního polynomu PN (x) pro stejný interval a stejnou váhovou funkci. Soustřeďme pozrnost na Gauss-Legendrovu kvadraturu, tj. na případ kdy je váhová funkce W(x) = 1. Ortonormální polynomy na intervalu −1 ≤ x ≤ 1 jsou dány rekurencí (j + 1)Pj+1 = (2j + 1)xPj − jPj−1. (116) Po nalezení souřadnic xi zbývá určit váhové koeficienty wj, které vypočítáme podle vztahu wj = < PN−1|PN−1 > PN−1(xj)PN (xj) . (117) 35 8 Numerické řešení obyčejných diferenciálních rovnic 8.1 Eulerova metoda Přímý útok na řešení problému y (t) = f(t, y) (118) je Eulerova metoda, která spočívá na hledání řešení ve tvaru Taylorova rozvoje do prvního řádu, tj. y(t + h) = y(t) + y (t)h + O(h2 ) (119) nebo ve tvaru y(t + h) = y(t) + f(t, y)h. (120) Pro fixní krok h převedeme vztah pro numerické výpočty na tvar yn+1 = yn + hf(tn, yn). (121) 8.2 Runge-Kutta (klasická RK4) Řešíme problém s počátečními podmínkami ve tvaru ˙y = f(t, y) a y(t0) = y0. (122) Pro fixní krokh > 0 definujeme yn+1 = yn + 1 6 h(k1 + 2k2 + 2k3 + k4) (123) 36 kde je k1 = f(tn, yn), (124) k2 = f(tn + h/2, yn + hk1/2), (125) k3 = f(tn + h/2, yn + hk2/2), (126) k4 = f(tn + h, yn + hk3). (127) 8.3 Runge-Kutta druhého řádu Pro jednoduchost bude odvozeny vzorce Rungeho-Kuttovy metody druhého řádu. Případě RK metod vyšších řádů je postup analogický ale mnohem pracnější. Uvažujme 1-D problém y (t) = f(t, y). (128) Řešení y(t) hledáme ve tvaru Taylorova rozvoje , tj. y(t + h) = y(t) + hy (t) + h2 2 y (t) + O(h3 ) (129) kde je y (t) = f(t, y) a y obdržíme diferenciací y , tj. y (t) = ∂f ∂t + ∂f ∂y y = ft(t, y) + fy(t, y)f(t, y). (130) Dostáváme tak řešení ve tvaru y(t + h) = y(t) + hf(t, y) + h2 2 [ft(t, y) + fy(t, y)f(t, y)] , (131) 37 kterou po úpravě převedeme na výraz y(t+h) = y(t)+ h 2 f(t, y)+ h 2 [f(t, y) + hft(t, y) + hfy(t, y)f(t, y)] . (132) Taylorův rozvoj funkce dvou proměnných do prvního řáde je f(t + h, y + k) = f(t, y) + hft(t, y) + kfy(t, y) + · · · . (133) Srovnáním s výrazem v závorce ve vztahu (132) dostáváme f(t+h, y+hf(t, y)) = f(t, y)+ h 2 ft(t, y)+ h 2 f(t+h, y+hf(t, y)) (134) a tak bude mít vztah pro RK2 krok tvar y(t + h) = y(t) + h 2 f(t, y) + h 2 f(t + h, y + hf(t, y)). (135) Zapsáno pro účely numerických výpočtů máme výsledek yn+1 = yn + h 1 2 k1 + 1 2 k2 , k1 = f(tn, yn), k2 = f(tn + h, yn + hk1). 8.4 Chyby (Demonstrováno na Eulerově metotě integrace) • Local Truncation Error (LTE) - rozdíl mezi numerickým řešením po ukončení jednoho kroku y1a exaktním řešením 38 v čase t1 = t0 + h, tj. y1 = y0 + hf(y0, t0) (numerické řešení), y(t0 + h) = y(t0) + hy (t0) + 1 2 h2 y (t0) + O(h3 ) (exaktní řešení rozvedené do Taylorova rozvoje), LTE = y(t0 + h) − y1 = 1 2 h2 y (t0) + O(h3 ). • Global Truncation Error (GTE) - je chyba v určitém časovém okomžiku t dosženého po provedené n kroků potřebných k tomu aby metoda dospěla od počátečního okamžiku t0k okamžiku t. Pro fixní krok h počet kroků bude n = t − t0 h . (136) Tento druh chyby je kumulativní efekt LTE který je v případě Eulerovy integrace ∼ h2 . To znamená, že po n krocích bude GTE ∼ nh2 ∼ h . 8.5 Adaptivní krok a RK4 Dosud byl integrační krop pevně daný. Pro získání výsledku s předepsanou přesností eps bylo nutné vhodně nastavit integrační krok h tak aby maximální chyba výsledku byla právě eps. Nevýhoda spočívá v tom, že při fixním kroku jsou určité části řešení určovány se zbytečně malým krokem (zbytečně velkou přesností) za cenu velkého počtu integračních kroků. 39 Lepší je v každém integračním kroku určit chybu výsledku a pokud bude větší než požadovaná přesnost tak zmenšit integrační krok a v opačném případě integrační krok zmenšit. Přímočarý algoritmus spočívá v "doubling"integračního kroku. Nechť je aproximace řešení ODE v bodě t + 2h reprezentováno veličinou y1. Určeme aproximaci řešení ODE v tomtéž bodě ale použitím dvou po sobě následujících RK4 kroků s krokem h. Tuto aproximaci označme y2. Chybu dvojnásobného kroku určuje rozdíl ∆ = |y2 − y1|. (137) Bude-li ∆ > EPS pak integrační krok zmenšíme na polovinu, tj. h → h/2, a zopakujeme výpočet y1 a y2. Pokud je ∆ < EPS tak další krok zdojnásobíme, tj. h → 2h, a určíme y1 a y2 v následujícím bodě. Následující pseudokód inlustruje myšlenku adaptivního kroku. 1 h=1e−6 2 t=0 3 y=y0 4 while teps then 19 h=h/2 40 20 y1=y 21 y2=y 22 else 23 h=2∗h 24 y=y1 25 break 26 end 27 end while 28 29 p r i n t t , y 30 31 end while 8.6 Adaptivní krok a embedované RK metody Cílem embedovaných RK metod je co nejefektivnějí vypočítat aproximaci řešení v p-tém a p + 1 - řádu. Nechť je ˜xn+1 přesné řešení diferenciální rovnice dx dt = f(t, x) (138) v n + 1 integračním kroku. Vyjádřeme jej pomocí aproximativních řešení x1 a x2 (první aproximace je určena krokem h a druhá krokem h/2) ˜xn+1 = x1 + Chp+1 + o(hp+2 ), (139) ˜xn+1 = x2 + 2C h 2 h p+1 + o(hp+2 ). (140) Odtud vidíme, že pro rozdíl x1 − x2 platí |x1 − x2| = Chp+1 (1 − 1 2p ) ⇒ C = |x1 − x2| hp+1(1 − 1/2p) . (141) 41 Označme = |x1 − x2| 1 − 2−p (142) a dostáváme výraz pro C ve tvaru C = hp+1 . (143) Odhadněme nyní největší optimální krok pro RK metodu při požadované toleranci 0. Jestliže bude krok hnew pak vztah mezi aproximací a skutečným řešením je ˜xn+1 = xn+1 + Chp+1 new. (144) Hledáme aproximaci s tolerancí 0, tj. ˜xn+1 = ˜xn+1 + 0. (145) To ale znamená, že musí platit Chnew ≤ 0 ⇒ hnew = 0 1 1+p h. (146) V praktických aplikacích se ukazuje, že je užitečné použít tzn. safety faktor β 1 který modifikuje odhad optimálního kroku takto hnew = β 0 1 1+p h. (147) Volba 0 závisí na daném pro problému, který řešíme. Často se volí 0 úměrné h. Tady se musíme mít na pozoru, protože ve fázi kdy snižujeme krok h z velkých hodnot pak nový krok hnew 42 určený exponentem 1/(1+p) nevede k požadované přesnosti. Pro > 0 je lepší použít exponent 1/p. Určení nového optimálního kroku hnew tak shrnuje vztah hnew = βh 0 1/(1+p) pro < 0 βh 0 1/(p) pro ≥ 0 (148) Pro embedovanou RK metodu čtvrtého řádu tak dostaneme ( p = 4 ) výraz pro nový krok ve tvaru hnew = βh 0 0.2 pro < 0 βh 0 0.25 pro ≥ 0 (149) Algoritmus pro RK s adaptivním krokem pak vypadá takto 1 VSTUP: t0 , x0 , eps0 , h , j =0, beta =0.8 2 3 KROK1: URCI x ( t+h , h ) , x (h+h , h/2) a eps 4 5 KROK2: POKUD JE eps