Numerické výpočty Jiří Zelinka Obsah 1 Polynomy 4 1.1 Homérovo schéma a základní úkony s polynomy........ 4 1.2 Interpolace............................. 9 1.3 Kořeny polynomů......................... 16 1.4 Splajny .............................. 24 1.5 Bersteinovy polynomy a Bézierovy křivky............ 30 2 Derivace 37 2.1 Limity............................... 37 2.2 Symbolické derivování...................... 40 2.3 Taylorův rozvoj.......... ................ 41 2.4 Numerické derivování....................... 43 3 Integrály 45 3.1 Symbolické integrování...................... 45 3.2 Numerické integrování...................... 47 4 Řady 52 4.1 Symbolické součty řad...................... 52 4.2 Fourierovy řady.......................... 54 5 Řešení nelineárních rovnic 57 5.1 Motivační úloha.......................... 57 5.2 Základní metody řešení nelineární rovnice ........... 61 5.3 Rychlost konvergence....................... 73 5.4 Řešení systému nelineárních rovnic............... 77 1 6 Funkce více proměnných 79 6.1 Zobrazovnání grafů........................ 79 6.2 Výpočet parciáních derivací................... 83 6.3 Integrování ve více proměnných................. 86 7 Řešení diferenciálních rovnic 89 7.1 Analytické řešení......................... 89 7.2 Numerické řešení......................... 92 8 Statistické metody 95 8.1 Intevalové odhady......................... 97 8.2 Testování hypotéz......................... 100 8.3 Lineární regrese.......................... 102 2 Úvod Tento text je zaměřen na základy numerické aproximace a řešení nelineárních rovnic, Bližší informace o teoretickém základu, na němž tento text staví může čtenář najít zejména ve skriptech [1] a [2], dále také v knihách [3] nebo [4]. 3 Kapitola 1 Polynomy Nebudeme se zda zabývat příliš teorií a základními pojmy jako např. stupeň polynomu apod. Literatura v tomto směru je velmi obsáhlá a čtenáři jistě nebude činit pražádné potíže si nějakou vyhledat. Nás budou zajímat především výpočty. První věc, na kterou bych rád upozornil, je nejednotnost zápisu polynomů v různé literatuře. Zpravidla se setkáme se zápisem P{x) = anxn + an_\xn~x + • • • + a\x + ao ale občas narazíme na opačné pořadí indexů, tedy na polynom ve tvaru P{x) = aoxn + aix™-1 + • • • + an_\x + an. Pro konkrétní polynom, kdy koeficienty nejsou vyjádřeny obecně ale čísly, je to samozřejmě jedno, ale pokud chceme pro polynom použít nějaké tvrzení, je potřeba si dávat pozor, v jakém pořadí ta či ona věta uvádí pořadí koeficientů. V této příručce se budu držet prvního způsobu, tedy indexy u koeficientů budou stejné jako příslušná mocnina x. Pokud bude potřeba zvýraznit stupeň polynomu, budeme používat značení Pn. 1.1 Hornerovo schéma a základní úkony s polynomy Na polynom budeme nahlížet jako na funkci - vrazíte do ní x a vypadne funkční hodnota podle daného vztahu. Základní věc, kterou tedy s polynomem potřebujeme dělat, je výpočet funkční hodnoty. Nejrychlejší metodou 4 k tomu je tzv. Homérovo schéma. Nebudu uvádět jeho podrobné odvození, které by pro čtenáře bylo snadným cvičením. Homérovo schéma je založeno na dělení polynomu P polynomem prvního stupně P±(x) = x — c, tedy P{x) = {x - c) ■ Q{x) + A, (1.1) kde Q je podíl (polynom stupně o jedna menší než P) a A je zbytek po dělení, který musí mít stupeň menší než dělitel, tedy je to konstanta. Z uvedeného zápisu je jasné, že P(c) = A, tedy tento zbytek po dělení je současně hodnota polynomu P v bodě c. Z (1.1) se porovnáním koeficientů u stejných mocnin dají odvodit vztahy pro koeficienty polynomu Q, Předpokládáme-li polynom Q ve tvaru Q(x) = fe^ii""1 H-----h bľx + b0 dostáváme postupně bn-i = 0>n b„-2 = an-i + c ■ bn_x fofc_i = ak + c-bk b0 = ai + c-bi A = a0 + c-b0. Hornerovo schéma zpravidla v literatuře nalezname coby následuj íccí tabulku: Q-n-2 ■ a2 ai a0 c bn-1 bn_2 b„-3 ■ ■ h bo A Pravidlo pro zapamatování výpočetního algoritmu je: „První číslo opíšeme (bn-i = an), každé další dostaneme tak, že k číslu nad čarou připočteme c násobek předchozího čísla (bk_i = ak + c • bk)LL. Hornerovo schéma se při ručních výpočtech nejčastěji používá k testování, zda dané číslo c je kořenem polynomu, v tom případě vyjde A = 0. Ale vidíme, že kromě hodnoty polynomu v bodě c dostáváme i podíl - polynom Q. Někoho možná napadne, co by se stalo, kdybychom Hornerovo schéma použili se stejnou hodnotou c znovu, tentokrát na polynom Q, pak znovu na výsledný podíl a tak dál. Pro tento účel si označíme polynom Q jako Qi a 5 hodnotu A jakožto Ao, v dalším kroku dostaneme podíl Q2 a hodnotu A±, a tak dál, takže obecně máme Qk(x) = {x - c) ■ Qk+i(x) + Ak. Homérovo schéma pak (symbolicky zkráceno) vypadá takto: p c Qi A0 c Q2 A1 c Q3 A2 c An Pro polynom P pak dostáváme P {x) = (x - c)Q!(x) + Ao = = (x - c)((x - c)Q2(x) + AJ + Ao = = (x - c)2Q2{x) + Ax(x - c) + A0 = = (x- c)2((x - c)Q3(x) + A2) + Ax(x -c)+A0 = = (x- c)3Q3(x) + A2(x - c)2 + Ax(x -c) + Ao = --- = = An(x - c)n + - c)™"1 + • • • + A^x -c) + A0 Hodnoty An,..., Ao jsou tedy koeficienty polynomu P posunutého do bodu c. Toto vyjádření se dá také velmi dobře použít pro výpočet derivací polynomu P v bodě c, ale to bychom předbíhali. V Matlabu definujeme polynom pomocí vektoru jeho koeficientů od nejvyšší mocniny, takže příkazem » P=[l 3 -2 0 5] definujeme polynom Pix) = x4 + 3x3 — 2x2 + 5, jak se dá snadno ověřit převodem polynomu na symbolický objekt: » poly2sym(P) ans = x"4 + 3*x"3 - 2*x"2 + 5 Pro výpočet hodnoty polynomu v bodě použijeme funkci polyval a můžeme ji aplikovat nejen na jednu hodnotu ale na vektor či matici hodnot, přičemž hodnoty polynomu se počítají pro každou služku zvlášť. Takže graf polynomu na intervalu můžeme získat třeba pomocí příkazů 6 » P=[l 3 -2 0 5]; » x=-3:0.01:3; >> y=polyval(P,x); >> plot(x,y) Operace s polynomy Mezi základní úkony, které lze s polynomy provádět, patří sčítání, násobení skalárem, vzájemné násobení a dělení se zbytek. První dvě operace patrně nepotřebují další komentáře, jen je potřeba vědět, že pro sčítání polynomů v Matlabu musí mít příslušné vektory koeficientů stejnou délku, takže je nutné je doplnit v případě potřeby nulami. Pro násobení polynomů se používá funkce conv, kde jako vstupní parametry zadáme polynomy, které chceme násobit. Dělení polynomu se zbytkem se provádí podobně, jako je tomu u celých čísel. Nejznámějším použitím dělení polynomu se zbytkem je Eukleidův algoritmus pro hledání největšího společného dělitele dvou polynomů. Algerotmus je všeobecně známý, proto jej zde nebudeme uvádět. Pro dělení se zbytkem v Matlabu je možné použít funkci deconv, která má dva vstupní a dva výstupní parametry. Na výstupu dostáváme podíl a zbytek po dělení. Například při dělení polynomu x4 + 3x3 — 4x + 1 polynomem x2 + 1 dostáváme » p=[l 3 0-41]; » Q=[l 0 1] ; » [D,R] =deconv(P,Q) D = 1 3 -1 R = 0 0 0-7 2 Zbytek po dělení R má formálně stejný stupeň jako dělenec P, aby platila rovnost P=D*Q+R. Tuhle skutečnost je potřeba brát v potaz při některých výpočtech v Matlabu, například právě u zmíněného Eukleidova algoritmu. 7 1.1.1 Derivace polynomu V této části trochu předběhneme učivo matematického kursu, protože derivace polynomu se nám bude hodit v dalším výkladu a navíc pro polynomy se počítá poměrně jednoduše. Obecně lze říci, že derivace funkce udává, jak rychle se funkce mění. Tedy pokud funkce hodně rychle roste, má derivace velké hodnoty, pokud se funkce na nějakém intervalu nemění (je konstantní), je tam její derivace nulová a tak podobně. Derivace je kladná tam, kde funkce roste, záporná, když funkce klesá, v bodech, kde má funkce lokální minimum nebo maximum, je derivace nulová. Pro polynom P{x) = anxn + an_\xn~x + • • • + a\x + ao se jeho derivace vypočítá podle vztahu P\x) = n ■ anxn~x + (n — 1) • an_\xn~2 + • • • + 2 • a> y=polyval(P,x); » dy=polyval(DP,x); >> z=zeros(size(x)); >> plot(x,y,x,dy,x,z,'k—') 8 Modře je zobrazen polynom, zeleně jeho derivace. 1.2 Interpolace Problém interpolace patří do teorie aproximace. Zpravidla se snažíme aproximovat nějakou funkci, z níž známe jenom hodnoty v diskrétních bodech, někdy je dokonce můžeme mít nepřesné. Lagrangeova interpolace Předpokládejme, že jsou dány body xt, i = 0,1,... ,n, xt ^ Xk pro í ^ k a hodnoty funkce / v těchto bodech: /(xj) = f i, i = 0,1,... ,n. Hledáme polynom Pn stupně nejvýše n takový, že Pn{Xi) = f i, í = 0,l,...,n. Body Xi, i = 0,1, ... ,n, Xi ^ Xk pro í ^ k, budeme nazývat uzly, polynom Pn interpolační polynom. Platí následující tvrzení, které se dá poměrně snadno dokázat: Pro in + 1) daných dvojic čísel (xi,fi), í = 0,l,...,n, Xi^Xkproi ^ k, existuje nejvýše jeden polynom Pn stupně menšího nebo rovného n takový, že Pn{Xi) = f i, í = 0,l,...,n. 9 Pro důkaz existence interpolačního polynomu použijeme Lagrangeovu konstrukci, podle níž se interpolační polynom nazývá Lagrangeův: Položme ^ ^ ^ (x ^o) • • • (x 1)(^ ^ž+l) • • • (x •En) "j^j" (x Xj) i^Xi Xq^ • • • (^ž — l)(^ž ^ž+l) • • • (^ž ^n) j=o (Xí Xj) Je zřejmé, že Zj je polynom stupně n a J II pro i = j. Definujme nyní Lagrangeův interpolační polynom Pn vztahem: n Pn(x) = lo(x)fo + h(x)fl + • • • + l„(x)f„ = J2lÁX)fc- i=0 Snadno se ověří, že tento polynom splňuje dané interpolační podmínky a je polynomem stupně nejvýše n. Polynomy Zj se nazývají Lagrangeovy fundamentální polynomy a tvoří bázi prostoru polynomů stupně nejvýše n. Podívejme se, jaké je výpočetní složitost Konstrukce Lagrangeova interpolačního polynomu. Při konstrukci jednoho fundamentálního polynomu začínáme polynomem prvního stupně x — xo, který postupně násobíme členy x—Xj, stupeň polynomu se postupně zvětšuje a pro jeho násobení potřebujeme řádově tolik operaci, jaký je jeho stupeň. Celkový počet operací pro konstrukci Zj tedy dostaneme jako součet konečné aritmetické řady, což je řádově n2 operací. Abychom sestrojili všechny fundamentální polynomy potřebujeme tedy řádově n3 operací. Při ruční konstrukci interpolačního polynomu ovšem zjistíme, že spoustu operací provádím opakovaně, takže při vhodném postupu by bylo možné konstrukci urychlit. A skutečně, optimální konstrukce Lagrangeova interpolačního polynomu je založena na funkci un+i(x) = ix — xq) . .. (x — xn) = n™=o(x ~ xi). Pro určení funkce uJn+i potřebujeme řádově také n2 operací, ale tuto funkci konstruujeme jen jedenkrát. Čitatel fundamentálního polynomu Zj získáme dělením ojn+i(x) : (x — xt), k čemuž se dá využít Hornerovo schéma, které je co se týče počtu operací lineární. 10 Dále se dá odvodit, že jmenovatel fundamentálního polynomu Zj je roven derivaci funkce con+i v bodě xt. Pro derivování polynomu je potřeba také řádově n operací a pro výpočet hodnoty derivace v bodě jakbysmet. Jednoduššeji se dá ovšem jmenovatel spočítat na základě faktu, že je roven hodnotě čitatele v uzlu xt. Určení hodnoty polynomu v bodě ale také dokážeme určit v linerání závislsti na stupni polynomu. Celkově tedy pro jeden fundamentální polynom potřebujeme jen lineární množství operací, pro všechny fundamentální polynomy je to pak řádově n2 operací. Matlabovská funkce pro konstrukci Lagrangeova interpolačního polynomu pak může vypadat následovně: function [lip,lfp] = laintpol2(uzly,ff) % function [lip,lfp] = laintpol(uzly,ff) 7o Lagrangeuv interpolační polynom % uzly,ff - radkove vektory 7o lip _ La.gr. interpol, polynom 1 lfp - matice fundamentalnich polynomu om=poly(uzly); % funkce omega omd=polyder(om); % omega' n=length(uzly)-1; lip=zeros(l,n+l); ifp=[]; for ii=0:n il=ii+l; xi=uzly(il); citatel=deconv(om,[1,-xi]); jmenovatel=polyval(omd,xi); li=l/jmenovatel*citatel; % fundamentálni polynom lfp=[lfp;li]; end % konec cyklu lip=ff*lfp; end % konec funkce Jako nepovinný druhý výstupní parametr funkce vrací fundamentální polynomy řádcích matice. Pro vyzkoušení programu zkusíme interpolovat hodnoty funkce sin. Nejprve definujeme uzly a určíme funkční hodnoty: 11 » uzly=[l 2 3 4 6]; >> ff=sin(uzly); » [lip,lfp] = laintpol2(uzly,ff); Protože máme spočtené i Lagrangeovy fundamentální polynomy, můžeme si je prohlédnout: » ll=lfp(l,:); » 12=lfp(l,:); » 12=lfp(2,:); » 13=lfp(3,:); » 14=lfp(4,:); » 15=lfp(5,:); » xx=0:0.01:7; » Ll=polyval(ll,xx); » L2=polyval(12,xx); » L3=polyval(13,xx); » L4=polyval(14,xx); >> L5=polyval(15,xx); » plot(xx,[L1;L2;L3;L4;L5]) >> grid on » axis([0,7,-5,5]) 0 1 2 3 4 5 6 7 12 Nakonec si necháme vykreslit i interpolační polynom spolu s funkcí sin pro porovnání: » Px=polyval(lip,xx); » fx=sin(xx); » plot(xx,Px,xx,fx,'— -',uzly,ff,'*') Vidíme, že polynom funkci aproximuje poměrně dobře, chyba je větší v oblasti, kde jsou uzly dále od sebe. Chyba se samozřejmě zvětšuje vně intervalu obsahující uzly. Hlavní nevýhodou Lagrangeovy konstrukce je ten fakt, že v případě, když potřebujeme přidat uzel, musíme přepočítat všechny fundamentální polynomy. Proto je v některých případech výhodnější Newtonova konstrukce, kterou si teď odvodíme. Je několik způsobů, jak to udělat, my použijeme tzv. iterovanou interpolaci. Iterovaná a Newtonova interpolace Označme Pi:i+i:...:i+k interpolační polynom na uzlech xl,... ,xl+k, tedy polynom stupně k. Pro k = 0 dostáváme polynom nultého stupně Pt, tedy se jedná o konstantní polynom, který je všude roven hodnotě fl. Polynom Pi,i+i,...,i+k sestrojíme z polynomů nižších stupňů následujícím způsobem: Nechť Pi:...,i+k-i a Pi+i,...,i+k Jsou interpolační polynomy na příslušných uz- 13 lech, oba stupně k — 1. Pak polynom Pi:i+i:...:i+k získáme ze vztahu Pi+l,...,i+k{x) Pi,...,i+k-l{x) _ \' _ v * v * _ v * = -í3—[Pí+iv..)í+fc(x) • (x - Xi) - Pív..)í+fc_i(x) • (x - xl+fc)]. xi-\-k xi Dosazením jednotlivých uzlů do této formule snadno zjistíme, že výsledný polynom skutečně splňuje interpolační podmínky Uvedený vztah trochu upravíme: Pi,...,i+k{x) = -^-[Pív..,j+fc_l(x) • {{Xl+k - Xi) - (X - Xi)) + + Pí+iv..)í+fc(x) • (X - Xi)] = P,..,+k-,{x) + ^--^Oe) -WiW (x _ gi) xi-\-k xi Rozdíl polynomů v čitateli zlomku je samozřejmě polynom stupně k — 1, označme jej například fž. Protože jak polynom Pjv..)í+fc_i tak i Pi+i:...:i+k mají stejné funkční hodnoty v uzlech xl+lr . . ,xl+k_i, jsou tyto uzly kořeny polynomu R. Protože těchto uzlů je A; — 1 lze polynom R zapsat ve tvaru R(x) = a ■ (x - xl+1) ■■■{x- xl+fc_i), kde a je koeficient u nejvyšší mocniny x, tedy u xfe_1. Tento koeficient je ale roven rozdílu koeficientů u nejvyšší mocniny polynomů Pi+i:...:i+k a Pi,...,i+k-i-Celkem dostáváme Pi,...,i+k{x) = Pi:...:i+k-i(x) H--^——(x -Xi)---(x- xl+fc_i). (1.2) xi-\-k xi Interpolační polynom na uzlech xl1... ,xl+k tedy dostaneme z interpolačního polynomu na uzlech xlv .. ,xl+k-i přidáním dalšího členu, který je roven fc-i součinu kořenových činitelů n (x—xi+j) vynásobenému koeficientem u nejvyšší 3=0 mocniny. Označme tento koeficient jako f[xl1... ,xl+k\. Je jasné, že f[xj\ = fl a dále ze vztahu (1.2) dostáváme rr n _ f[xi+li ■ ■ ■ i xi+k] — f[xi, • • • i xi+k-l] J lxi i • • • i xi+k\ • xiJi-k xi 14 Pí, .,i+k X Způsob výpočtu těchto koeficientů jim dává také název - označujeme je jako poměrné diference. Jejich výpočet pro všechny uzly lze provést pomocí následujícího schématu: fo > f[xo,xi] fi > f[xo,xi,x2] > f[xi,x2] Í2 f[x0,...,Xn] ^ f [Xn—2 i Xn—i, f n Postupným použitím vztahu (1.2) pak dostáváme Po,...,nix) = flxo\ + f[xo,xi]{x ~ xo)+ + f[xo, xi, x2](x — x0)(x — xi)+ (1.3) + • • • + f[xo, ■ ■ ■ ,xn](x - x0) ■ ■ ■ (x - X„_i). Tato konstrukce se nazývá Newtonův interpolační polynom. Výsledek samozřejmě musí být totožný s Lagrangeovým interpolačním polynomem. Hlavní výhodou Newtonovy konstrukce ale je snadné přidávání dalších uzlů - stačí vypočítat další řádek v tabulce pro výpočet poměrných diferencí a k předchozímu interpolačnímu polynomu přidat součin příslušných kořenových činitelů násobený výslednou poměrnou diferencí. Pro interpolaci v Matlabu je možné použít také jeho vlastní funkci polyf it, která je ale trochu obecnější a umožňuje najít také aproximaci dat polynomem nižšího stupně pomocí metody nejmenších čtverců. Intrpolovat lze samozřejmě take v Sage. Abychom dostali jako výsledek racionální polynom, použijeme pro příklad jiná data než výše uvedená: sage: P = PolynomialRing(QQ, 'x') sage: P.lagrange_polynomial([(0,1),(2,2),(3,-2),(-4,9)]) -23/84*x"3 - ll/84*x"2 + 13/7*x + 1 Dají se zde určit i poměrné diference pro Newtonův interpolační polynom: xo X\ x2 f[xo] f[xi] f[x2] f[Xn] = 15 ságe: P.divided_difference([(O,1),(2,2),(3,-2),(-4,9)]) [1, 1/2, -3/2, -23/84] Předpokádám, že čtenáři nečiní žádné potíže si ověřit, že z vypočítaných nodnot pomocí Newtonovy konstrukce získáme stejný polynom jako v prvním případě. 1.3 Kořeny polynomů Najít kořen daného polynomu je jednou ze základních matematických úloh. K ověřování, zda dané číslo je nebo není kořenem polynomu zpravidla používáme Hornerovo schéma, existují ale další nástroje, které nám mohou hledání kořenů usnadnit. Například pokud máme polynom s celočíselnými koeficienty a zajímají nás racionální kořeny, tedy kořeny ve tvaru |, kde p je celé číslo, q je přirozené číslo a p a q jsou nesoudělná, tak se dá lehce zjistit, že koeficient u nejvyšší mocniny polynomu (tedy an) musí být dělitelný číslem q, zatímco absolutní člen (t.j. a0) musí být dělitelný p. Další pomůckou může být stanovení oblasti, kde všechny kořeny leží. K tomu nám může pomoci například tvrzení, které lze najít i s důkazem v [1]. Nechí A = max (|a0|, • • • , , B = max (|ai|,..., \an\), kde cik, k = 0,1,..., n, aoQ™ 7^ 0, jsou koeficienty polynomu P, Pak pro všechny kořeny xk, k = 0, 1, ..., n, polynomu P platí 1 , , A , , <\xk\ 1, pak tento kořen je i kořenem derivace polynomu P s násobností rij — 1. Odtud plyne, že snížit násobnost všech kořenů na 1 lze tak, že určíme největší společný dělitel D polynomu P a jeho derivace P'. Kořeny tohoto největšího společného dělitele jsou právě ty kořeny původního polynomu, které mají násobnost větší než jedna a jejich násobnost v děliteli D je o jedna menší než v polynomu P. Vydělením polynomu P dělitelem D získáme tedy polynom, který má stejné kořeny jako P, ale všechny jsou násobnosti 1. Z numerického hlediska je potřeba poznamenat, že tento postup v praxi funguje dobře pro polynomy se celočíselnými koeficienty, které se pohybují v „rozumných" mezích. Ale i v jednoduchých případech je potřeba postupovat obezřetně: >> format long » P=poly([l 111]) p — r 1-4 6 -4 1 » DP=polyder(P) DP = 4 -12 12 -4 » roots(P) ans = 1.000217162530750 0.999999969335793 + 0.000217131857745i 0.999999969335793 - 0.000217131857745i 0.999782898797672 » roots(DP) ans = 1.000004930958320 + 0.000008540746793i 1.000004930958320 - 0.000008540746793i 0.999990138083364 Vidíme, že pokud bychom posuzovali shodnost kořenů u uvedeného polynomu a jeho derivace, přičemž ani nepožadujeme příliš velkou přesnost, tak nejenže ke shodě nedochází, ale kořeny ani nejsou násobné. Nicméně pokud bychom hledali největší společný dělitel, postup povede ke správnému výsledku: 1» [Q,R]=deconv(P,DP) 18 Q = O.250000000000000 -0.250000000000000 R = 0 0 0 0 0 >> D=Q % D je nejv.spol.dělitel D = 0.250000000000000 -0.250000000000000 » roots(D) ans = 1 Vidíme, že výsledek je přesný, a to i navzdory numerickým nepřesnostem během výpočtu. Ukažme si ještě jeden příklad, kde budou dva násobné kořeny blízko sebe: » P=poly([0.9 0.9 1.1 1.1 1.1]) P = 1.0000 -5.1000 10.3800 -10.5380 5.3361 -1 0781 » DP=polyder(P) DP = 5.0000 -20.4000 31.1400 -21.0760 5.3361 » [Ql,Rl]=deconv(P,DP) Ql = 0.2000 -0.2040 Rl = 0 0 -0.0096 0.0298 -0 0306 0 0105 >> Rl(l:2)=[] % odstránime pocatecni nuly Rl = -0.0096 0.0298 -0.0306 0 0105 » [Q2,R2]=deconv(DP,Rl) Q2 = -520.8333 510.4167 R2 = 1.0e-ll * 0 0 -0.1766 0 2515 -0 0769 V tomto místě je nutné usoudit, že zbytek po dělení je ve skutečnosti nulový, tudíž nějvětší společný dělitel je předchozí zbytek, tedy polynom Rl. 19 » D=R1 D = -0.0096 0.0298 -0 0306 0.0105 » P0=deconv(P,D) PO = -104.1667 208.3333 -103 1250 » roots(PO) ans = 1.1000 0.9000 >> format long » roots(P0) ans = 1.100000000001784 0.899999999998279 Vidíme, že výsledek je opět poměrně přesný. Kořeny původního polynomu Matlab spočítá s daleko větší chybou: » roots(P) ans = 1.100021372535598 + 0 000037011183255i 1.100021372535598 - 0 000037011183255i 1.099957254925198 0.900000434848828 0.899999565154781 1.3.1 Separace reálných kořenů Ukážeme si algoritmus, s jehož pomocí je možné přesně určit počet reálných polynomů v daném intervalu, pokud předpokládáme, že všechny reálné kořeny jsou jednoduché. Tento algoritmus je založen na konstrukci tzv. Sturmovy posloupnosti, která se definuje následovně: Posloupnost reálných polynomů P = Po, Pl, ■ ■ ■ , Pfn se nazývá Sturmovou posloupností příslušnou polynomu P, jestliže 20 1. Všechny reálné kořeny polynomu Pq jsou jednoduché. 2. Je-li xq reálný kořen polynomu Pq, pak signPi{xo) = —sígnP^xo). 3. Jestliže Xq je reálný kořen polynomu Pt, platí Pl+1{x0)Pl_1(x0) < O, pro i = 1, 2,. .., m — 1. 4- Poslední polynom Pm nemá reálné kořeny. Sturmovu posloupnost lze zkonstruovat poměrně snadno, Pokud předpokládáme, že výchozí polynom P = Pq nemá reálné kořeny, stačí položit Pi = —Pq, abychom zaručili splnění druhé vlastnosti. Třetí vlastnost z definice dostaneme, pokud polynom Pl+i vezmeme jako záporně vzatý zbytek po dělení : Pl5tedy Pi-i = P% - Q — Pi+i Pro kořen xq polynomu Pj tedy máme P1_i(xq) = —P1+i(xq), přičemž daný výraz nemůže být nulový, jinak bychom indukcí dostali, že xq je kořenem Po i Pi, což je spor s tím, že kořeny Po Jsou jednoduché. Konstrukce založená na dělení se zbytkem zaručuje, že stupně polynomů v posloupnosti postupně klesají. Poslední polynom Pm je zpravidla konstantní, ale je možné konstrukci ukončit i dříve, pokud víme, že polynom nemá reálné kořeny, např. pro polynom x2 + 1. Počet kořenů v daném intervalu pak lze zjistit pomocí Sturmovy věty: Počet reálných kořenů polynomu P v intervalu [a, b] je roven Wib) — W(a), kde W(x) je počet znaménkových změn ve Sturmově posloupnosti Pq(x), ... , Pm(x) v bodě x (z níž jsou vyškrtnuty nuly). Předpokládám, že pojem počet znaménkových změn je natolik intuitivně jasný, že nepotřebuje definici. Důkaz Sturmovy věty je založen na faktu, že k navýšení či snížení počtu znaménkových změn může dojít jen při přechodu přes kořen některého z polynomů ve Sturmově posloupnosti. Ze druhé vlastnosti v definici Sturmovy posloupnosti skutečně plyne, že pokud x roste, pak při přechodu přes kořen P = P0 jsou dvě možnosti: buď přecházíme ze záporných hodnot polynomu do kladných, v tom případě polynom Po roste, jeho derivace je tedy v okolí kořene kladná, tudíž Pi je tam záporný a přibude 21 jedna znaménková změna. Nebo Po v kořenu klesá, takže Pi je kladný a také přibude jedna znaménková změna. x — h x x + h Po Pl - 0 + W(x) 0 0 1 x — h x x + h Po Pi + 0 -+ + + W(x) 0 0 1 Při přechodu přes kořen polynomu Pl pro i > 0 jsou celkem čtyři možnosti, při žádné z nich však podle třetí vlastnosti z definice nedochází ke změně počtu znaménkových změn: x — h x x + h P-1 Pl Pl+1 - 0 + + + + W(x) 1 1 1 x — h x x + h P-i P Pi+i + 0 -+ + + W(x) 1 1 1 x — h x x + h P-i Pi Pi+1 + + + - 0 + W(x) 1 1 1 x — h x x + h P-i P Pi+1 + + + + 0 - W(x) 1 1 1 K navýšení počtu znaménkových změn tak dochází jenom při přechodu přes kořen polynomu Pq, odkud bezprostředně plyne tvrzení věty. Příklad 1. Zjistíme počet kladných a záporných kořenů polynomu x4 — 4x2 + 8x — 2. Nejprve sestrojím Sturmovu posloupnost: » P0=[1 -4 0 8 -2] ; » Pl=-polyder(P0) Pl = -4 12 0 Protože nás zajímají jen znaménka a jejich změny, je možné polynom vydělit nebo vynásobit kladným číslem. To se hodí hlavně při ručním počítání, když se chceme vyhnout složitějším zlomkům. Pak dál pokračujeme v konstrukci Sturmovy posloupnosti: provedeme dělení se zbytkem, odstraníme nulové koeficienty a polynom opět můžeme vydělit kladným číslem: » Pl=Pl/4 Pl = 22 -1 3 0-2 » [Q,R] =deconv(P0,Pl) Q = -1 1 R = 0 0-360 » P2=-R/3; » P2(l: 2) = [] P2 = 1 -2 0 Celou činnost opakujeme, dokud nedostaneme konstantní polynom: » [Q,R] =deconv(Pl,P2) Q = -1 1 R = 0 0 2-2 » P3=-R/2; » P3(l: 2) = [] P3 = -1 1 » [Q,R] =deconv(P2,P3) Q = -1 1 R = 0 0 -1 » P4=l; Podle předchozího textu je horní hranice pro absolutní hodnotu všech kořenů rovna 9, takže všechny kořeny leží v intervalu [—9,9]. Při ručním počítání je jednodušší určit znaménko limity hodnoty v ±oo podle parity nejvyšší mocniny a znaménka jejího koeficientu. Počet znaménkových změn v Matlabu určíme pomocí příkazů: » x=-9; » [polyval(PO,x),polyval(Pl,x),polyval(P2,x),... polyval(P3,x),polyval(P4,x)] 23 ans = 9403 970 99 10 1 » x=0; » [polyval(PO,x),polyval(Pl,x),polyval(P2,x),... polyval(P3,x),polyval(P4,x)] ans = -2-2011 » x=9; » [polyval(PO,x),polyval(Pl,x),polyval(P2,x),... polyval(P3,x),polyval(P4,x)] ans = 3715 -488 63 -8 1 Při ručním počítání stačí určovat znaménka a výsledky můžeme přehledně umístit do tabulky: X P0(x) Pl(x) P2{x) P3(x) PA{x) W(x) —oo + + + + + 0 0 — - 0 + + 1 oo + - + . + 4 Celkem tedy máme 4 kořeny, z toho je jeden záporný a tři kladné. 1.4 Splajny Teď zase trochu předběhneme, budeme totiž opět potřebovat derivace, o kterých budeme podrobněji mluvit později, ovšem výklad o splajnech dobře navazuje na interpolaci, takže jsem se rozhodl jej zařadit už sem. Snad to nebude příliš vadit, protože tato část bude spíše jen doplňková a informativní. V češtině už se pro funkce, o nichž teď budeme mluvit, zaběhl počeštěný výraz splajn. Ve starší literatuře je možné narazit na anglický termín splíne, výslovnost je ovšem stejná jako u počeštěného výrazu. Splajny patří mezi interpolační techniky, protože většinou procházejí přesně zadanými hodnotami. Existují ale i tzv. vyhlazovací splajny, o těch tady v tuto chvíli taktně pomlčíme. Splajn se zpravidla definuje jako po částech polynomiální funkce, která je spojitá do určitého řádu. Nejjednodušším příkladem splajnu je spojitá po částech lineární funkce, která prochází zadanými body v rovině. Přesněji 24 řečeno předpokládejme opět, že jsou dány body xt, i = 0,1,.. ., n, xt ^ Xk pro i 7^ k a hodnoty funkce / v těchto bodech: fixi) = ft, i = 0,1, ... ,n. Lineární spojitý splajn je funkce s, která je spojitá, lineární na každém intervalu [xl, i = 0,..., n — l a platí s(xt) = fl, i = 0,..., n. Na následujícím obrázku můžeme vidět příklad takové funkce: 1 2 3 4 5 6 7 Nejčastěji se setkáme s kubickými splajny. Jedná se o funkce, které jsou po částech polynomy třetího stupně a jsou spojité do řádu dva, jsou tedy spojité a mají spojité i první a druhé derivace, Uvedeme radši opět přesnou definici: Nechť jsou dány body xt, i = 0, 1, ..., n, xt ^ Xk pro í ^ k (zvané uzly) a hodnoty funkce / v těchto bodech: /(xj) = /l5 i = 0,1,..., n. Kubický splajn je funkce s, která je spojitá včetně první a druhé derivace a je kubickým polynomem na každém intervalu [xl,xl+i], i = 0,... ,n — 1 a platí s(xi) = fi, i = 0, ... ,n. Podmínky spojitosti musíme uplatnit v uzlech, jelikož polynomy mají spojité derivace všech řádů, takže uvnitř intervalů není se spojitostí problém. Máme celkem n intervalů, na každém potřebujeme sestrojit polynom 3. stupně, který má čtyři neznámé koeficienty, dohromady tedy máme An neznámých koeficientů. Ovšem pro každý z n polynomů máme zadána dvě funkční hodnoty (jednu v každém krajním bodě intervalu [xl,xl+i]), tím se nám počet neznámých koeficientů redukuje na polovinu a současně zajistíme spojitost kubického splajnu. Ve vnitřních bodech (x±,. .., x„_i) dále máme podmínky na spojitost první a druhé derivace, čímž se nám počet neznámých koeficientů zredukuje ne 2. Další podmínky už ale nemáme, takže kubický splan 25 není funkčními hodnotami jednoznačně zadán a je potřeba dvě podmínky přidat,aby bylo možné jej sestrojit. Nejčastěji se zadávají hodnoty první derivace v krajních bodech xq a sn, pak hovoříme o úplném kubickém splajnu, nebo hodnoty druhé derivace v těchto bodech. Pokud mají být tyto hodnoty nulové, splajn se nazývá přirozený kubický splan. Nebudeme tady uvádět přesnou konstrukci kubických splajnů, zájemce najde podrobné informace třeba v [2]. Ukážeme si, jak zkonstruovat splajn za pomocí software. V Matlabu se pro konstrukci splajnu dá použít základní funkce spline. V nejjednodušší verzi počítá koeficienty polynomů, jako výsledek pak vrací složitější strukturu, která obsahuje víc informací. Pro data z předchozího obrázku dostaneme: » x=[l 2 3 4 6 7] ; » f=[3 1 0 2 1 1] ; » plot(x,f,'-*') >> print -depsc splajnl.eps >> s=spline(x,f) form: 'pp' breaks: [12 3 4 6 7] coefs: [5x4 double] pieces: 5 order: 4 dim: 1 Označení 'pp' ukazuje, že se jedná o po částech polynomiální funkci, Dále následují uzly, tedy hraniční body pro jednotlivé části polynomu. Pak můžeme zjisti koeficienty na jednotlivých intervalech, počet intervalů a stupeň polynomu (ve skutečnosti spíše počet koeficientů na každém intervalu). Význam poslední položky by měl jasný. Pokud bychom tedy chtěli znát koeficienty splajnu na jednotlivých intervalech stačí zadat >> s.coefs ans = s = 0.6978 0.6978 -1.4891 -1.5935 0.5000 2.5935 -1.1044 -2.1978 0.8956 3.0000 1.0000 0 26 0.4081 -1.8738 1.6153 2.0000 0.4081 0.5748 -0.9829 1.0000 Nás budou ovšem spíš zajímat hodnoty splanu, abychom si jej mohli případně nakreslit. K tomu lze použít funkci ppval: » xx=l:0.01:7; >> ss=ppval(s,xx); >> plot(xx,ss,x,f,'r*') V případě, že nás nezajímají koeficienty na jednotlivých intervalech, ale jen výsledné funkční hodnoty, stačí zadat >> ss=spline(x,f,xx); Zvídavého čtenáře teď možná napadne, že jsme při konstrukci zadávali jen funkční hodnoty a žádné okrajové podmínky. Právem se tedy může ptát, jaké okrajové podmínky byly použity. Tuto otázku jistě uspokojivě zodpoví dokumentace Mat labu. V Sage se splajn definuje pro body v rovině, jimiž má procházet, takže pro stejná data jako výše můžeme použít následující postup: 27 sage: values = [[1,3],[2,1],[3,0], [4,2], [6,1], [7,1]] sage: S = spline(values) sage: S [[1, 3] , [2, 1] , [3, 0] , [4, 2] , [6, 1] , [7, 1]] Vidíme, že samotný obsah proměnné S nám moc informací nedá. Příkazem sage: type(S) sage.gsl.interpolation.Spline aspoň zjistíme, o jaký datový typ se jedná, ale koeficienty na jednotlivých intervalech nezjistíme. A pokud vím, u této základní funkce se tyto koeficienty přímo zjistit nedají. Pokud by je člověk opravdu potřeboval, musí si stáhnout další knihovny, které obsahují poněkud více sofistikované funkce pro práci se splajny. Pro naše účely ovšem základní funkce zcela postačuje, například hodnotu polynomu v bodě zjistíme snadno: sage: S(l) 3.0 sage: S(1.5) 1.9917763157894737 sage: S(sqrt(2)) 2.1640477491461865 Nechat si splajn vykreslit (případně i s uzly a příslušnými hodnotami) taky není problém: plot(S,(1,7))+list_plot(values,color='red') 28 V Sage je velmi jednoduché změnit hodnotu v některém uzlu: sage: S[0] [1, 3] sage: S[0] = [l,2] sage: S [[1, 2], [2, 1], [3, 0], [4, 2], [6, 1], [7, 1]] sage: plot(S,(1,7))+list_plot(values,color='reď) 2.5 29 První hodnota samozřejmě neleží na splajnu, jelikož jsem ji nepředefinovali. Splajn ale prochází danými body. Jednoduše se dá taky bod ubrat nebo naopak přidat: sage del S[l] sage S [[1, 2] , [3, 0] , [4, 2] , [6, 1], [7, 1]] sage S.append([5,2]) sage S [[1, 2] , [3, 0] , [4, 2] , [6, 1], [7, 1], [5, 2]] Když si necháme splajn zobrazit, vidíme, že nazáleží na pořadí uzlů. 1.5 Bersteinovy polynomy a Bézierovy křivky Sergej Bernstein byl ruský matematik, který se během svého dlouhého života dosáhl významných výsledků v několika matematických odvětvích. Jeho jméno by se mělo vyslovovat rusky, tedy „bernštejn11 a někdy se takto i v češtině píše. Často se ovšem setkáme s výslovností německou („bernštajn11) nebo s anglickou („bernstain11). Nás bude zajímat jeho příspěvek k teorií aproximací. Bernsteinovy polynomy totiž mají tu důležitou vlastnost, že s libovolnou přesností aproximují spojitou funkci na intervalu [0,1] Jednoduchou substitucí pak jimi můžeme 30 aproximovat spojitou funkci na libovolném uzavřeném intervalu. Jejich konstrukce je navíc velmi jednoduchá a tedy i snadno naprogramovatelná. Nechť tedy n je pevně dané přirozené číslo a k je přirozená číslo nebo nula, k < n. Bersteinův bázový polynom bn^ definujeme vztahem bn,k{x)=(fyxk{\-x)n-k. Dále nechť / je funkce definovaná na intervalu [0,1], položme fk = /(^) pro k = 0,.. ., n. Bernsteinův polynom stupně n funkce / definujeme výrazem n B„j(x) = Y^fk- b„,k(x)- i1-5) k=0 Snadno se ověří, že n Yl h ■ b„,k(x) = 1, k=0 odkud pak bezprostředně plyne, že pro funkci / konstantní na [0,1] platí fix) = Bnj(x) pro každé x E [0,1]. Totéž platí i v případě, že / je polynomem prvního stupně, jak se dá též poměrně snadno spočítat. Pro polynomy vyšších stupňů to už ale neplatí, jak dokládá následující obrázek, kde čárkovaně je zobrazena funkce x2, křivky, které se k ní blíží jsou postupně Bernsteinovy polynomy druhého, třetího a čtvrtého stupně. 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 31 Tento obrázek byl získán pomocí matlabovského programu pro výpočet Bersteinova polynomu, přesněji řečeno jeho koeficientů. Jako vedlejší produkt dostáváme i Bernsteinovy bázové polynomy uspořádané v matici. Při jejich konstrukci se využívá toho, že daný bázový polynom má kořeny 0 a 1 násobností k a n — k. function [BP,B] = bern_pol(fk) % function BP = bern_pol(fk) % Bernsteinuv polynom pro hodnoty fk 7o B - matice Bern. bazových polynomu n=length(fk)-l; B= [] ; 7o matice bazových Bernsteinovych polynomu for k=0:n bk=nchoosek(n,k)*poly( [zeros(1,k),ones(1,n-k)]); bk=bk*(-l)~(n-k); B=[B;bk] ; end BP=fk*B; end Uvedený obrázek pak dostaneme například pomocí následujícího postupu: » f=inline('x."2'); » n=2; » k=0:n; » xk=k/n; » f2=f(xk); » B2=bern_pol(f2); » n=3; » k=0:n; » xk=k/n; » f3=f(xk); » B3=bern_pol(f3); » n=4; » k=0:n; » xk=k/n; » f4=f(xk); 32 » B4=bern_pol(f4); » xx=0:0.01:l; » Bx2=polyval(B2,xx); >> Bx3=polyval(B3,xx); » Bx4=polyval(B4,xx); » plot(xx,[Bx2;Bx3;Bx4],xx,xx."2,'—') Nej důležitější vlastností Bernsteinových polynomů, jak už bylo naznačeno, je, že pro funkcí / spojitou na [0.1] konverguje posloupnost Bnj stejnoměrně Kromě teoretického významu se Bernsteinovy polynomy používají i prakticky, totiž ke konstrukci Bězierových křivek. „Legenda" praví, že tyto křivky objevili nezávisle na sobě dva pracovníci konstrukčních kancelářích dvou renomovaných francouzských automobilek. První objevitel ovšem přišel o prvenství tím, že kvůli konkurenčnímu boji jeho objev podléhal utajení. V dnešní době jsou Bézierovy křivky běžnou součástí téměř každého kreslícího software, často bez uvedení, o jaký typ křivek se jedná. Nejčastěji můžeme narazit na kubické Bézierovy křivky, jejichž konstrukci si zde ukážeme Zobecnění pro vyšší stupně je zcela přímé a čtenář by je jistě bez problémů zvládl. Bézierovy křivky jsou zadány parametricky, konstrukce je navíc v podstatě totožná v rovině i v trojrozměrném prostoru. Mějme tedy 4 body v rovině, označme je Po,Pi,P2,P3, Pt = [xl,yl], tyto body se nazývají kontrolní nebo řídící. Bézierova křivka je množina všech bodů Q = [x, y) pro něž platí kde proměnná t probíhá interval [0,1] Pokud bychom uvažovali body v prostoru, t.j. Pj = [xl,yl,zl], přibyla by nám rovnice pro třetí souřadnici bodu x = x(t) xo ■ h,o{t) + xi ■ bz,i{ť) + x2 ■ k^2{t) + x3 ■ b3;i(t) x0(l - tf + 3xiŕ(l - ŕ)2 + 3x2ŕ2(l - t) + x3t3 Vo ■ h,o{t) + yi ■ h,i{t) + 2/2 • h,2{t) + ž/3 • h;í{t) y0(l - tf + 3yií(l - ŕ)2 + 3y2ŕ2(l - t) + y3ŕ, y = y{t) Q z = z(ť) zo ■ b3,o{t) + zi ■ 63)i(ŕ) + z2 ■ bzfliť) + z3 ■ b3;i(t) zq(1 - tf + 3ziŕ(l - ŕ)2 + 3z2ŕ2(l - t) + z3ŕ 33 Souřadnice bodu Q na Bézierově křivce jsou tedy hodnoty Bersteinova polynomu pro příslušné souřadnice kontrolních bodů. Stručně se body na Bézierově křivce dají vyjádřit jako Q = Q(t) = P0(l - tf + 3Pií(l - ŕ)2 + 3P2ŕ2(l - t) + P3í3 Z vyjádření navíc plyne, že pro t = 0 dostaneme výchozí kontrolní bod Pq, pro t = 1 koncový kontrolní bod P3. Sestrojit Bézierovu křivku v Matlabu je velmi jednoduché, není ani potřeba vytvářet pro tento účel zvláštní funkci: » P0=[0;1]; Pl=[2;3] ;P2=[4;2] ;P3= [4;0]; » t=0:0.01:1; » Q=P0*((l-t).~3)+3*Pl*(t.*(l-t). ~2) + .. . 3*P2*(t."2.*(l-t))+P3*(t."3); +» plot(Q(l, :),Q(2, :)) Pokud chceme zobrazit i kontrolní body, nejjednodušší je poskládat je do matice, kterou zobrazíme po řádcích: >> hold on » PP=[P0,P1,P2,P3]; » plot(PP(l,:),PP(2,:),'r*') 3i-1-1-1-»- 2.5 - )l-1-1-1-1-1-1-1-J, 0 0.5 1 1.5 2 2.5 3 3.5 4 34 Pro konstrukci Bézierových křivek se dá použít také algoritmus de Casteljau pojmenovaný po jednou z objevitelů Bézierových křivek. Spočívá v rozdělení spojnic mezi kontrolními body v daném poměru. Tím získáme tři body, jejichž spojnice opět rozdělíme ve stejném poměru, Stejně tak rozdělíme i spojnici takto získaných dvou bodů a tím získáme bod, který leží na Bézierově křivce. Následující obrázek ukazuje tento proces pro dělící poměr 1:1 a 1:2, to znamená, že všechny spojnice dělíme na poloviny, respektive na třetinu a dvě třetiny. 0.5 - 0I-1-1-1-1-1-1-1-J, 0 0.5 1 1.5 2 2.5 3 3.5 4 To, že uvedený postup skutečně dává body na Bézierově křivce, snadno dokážeme následujícím výpočtem. Dělící poměr pro rozdělení spojnice dvou bodů je určen hodnotou a 6 (0,1), např. bod C na spojnici bodů A a B se dá vyjádřit jako C = a ■ A + (1 — a) ■ B, přičemž výsledný bod je blízko bodu A pro a blízké jedné. Délky ůdeček AC a C B jsou pak v poměru (1 — t) : t. Položme tedy R0 = a-P0 + (l-a)-P1, Ri = a ■ P1 + (1 - a) ■ P2, R2 = a-P2 + (l-a)-P3 čímž získáme první trojici bodů. Pokračujme dále: 5*0 = a ■ Rq + (1 — a) ■ Ri, Si = a ■ Ri + (1 — a) ■ i?2 35 a konečně T = a ■ Sq + (1 — a) ■ S±. Zpětným dosazováním dostáváme T = a ■ Sq + (1 — a) ■ S± = a ■ (a ■ R0 + (1 - a) ■ Rx) + (1 - a) ■ (a ■ Rx + (1 - a) ■ R2) = a2 ■ R0 + 2a(l - a) ■ i?i + (1 - a)2 ■ R2 = a2 ■ (a ■ P0 + (1 - a) ■ + 2a(l - a) ■ (a ■ P1 + (1 - a) ■ P2) + + (1 - a)2 ■ (a ■ P2 + (1 - a) ■ P3) = a3 • P0 + 3a2(l - a) ■ P1 + 3a(l - a)2 ■ P2 + (1 - a)3 ■ Pz. Vidíme, že pro t = 1 — a jsem dostali přesně vyjádření bodu na Bézierově křivce. 36 Kapitola 2 Derivace 2.1 Limity Ještě před samotnými derivacemi se aspoň krátce zmíníme o limitách. V Sage to není problém, stačí zadat například sage: n=var('n' ) sage: limit((1+1/n)~n,n=oo) e sage: limit((1-1/n)~n,n=oo) e~(-l) sage: limit((1+10/n)~n,n=oo) e"10 sage: x=var('x') sage: limit((1+x/n)~n,n=oo) e"x sage: limit(sin(x)/x,x=0) 1 V Matlabu symbolické limity nenajdeme, nicméně některé limity pro n —> oo se dají spočítat nebo odhadnout tak, že daný výraz se pokoušíme spočítat pro hodně velká čísla. Má to ale svoje úskalí. Pokud bychom se pomocí výše uvedené limity pokoušeli určit Eulerovo číslo s přesností na 6 desetinných míst, dostáváme: >> formát long » n=l; 37 » eO=(l+l/n)~n eO = 2 » n=10; » el=(l+l/n)"n el = 2.593742460100002 » while abs(el-eO)>0.000001, eO=el; n=10*n;... el=(l+l/n)~n, end el = 2.704813829421528 el = 2.716923932235594 el = 2.718145926824926 el = 2.718268237192297 el = 2.718280469095753 el = 2.718281694132082 el = 2.718281798347358 Pokud bychom ale chtěli tímto způsobem určit Eulerovo číslo s přesností na 15 desetinných čísel, což je zhruba přesnost, s jakou jsou běžně v počítači uložena čísla řádově v jednotkách, vyjde » while abs(el-e0)>0.000000000000001, e0=el; n=10*n;... el=(l+l/n)"n, end el = 2.704813829421528 el = 2.716923932235594 el = 2.718145926824926 el = 2.718268237192297 38 el = 2.718280469095753 el = 2.718281694132082 el = 2.718281798347358 el = 2.718282052011560 el = 2.718282053234788 el = 2.718282053357110 el = 2.718523496037238 el = 2.716110034086901 el = 2.716110034087023 el = 3.035035206549262 el = 1 el = 1 Pro velká n je totiž hodnota 1/n menší než přesnost, s jakou je uložena jednička, takže při sčítání 1 + 1/n dostáváme výsledek 1. Pokud bychom chtěli přesto s pomocí Matlabu Eulerovo číslo určit jako limitu nějakého nekonečného procesu, můžeme využít vztah oo -i n=0 X- Součet nekonečné řady je limita částečných součtů, které můžeme snadno vyjádřit: » n=0;a=l;s=l; » while a>0.000000000000001, n=n+l;a=a/n;s=s+a; end » s 39 s = 2.718281828459046 » exp(l) ans = 2.718281828459046 » n n = 18 Ve výpisu vidíme, že dosažená hodnota je v rámci přesnosti stejná jako hodnota, se kterou Matlab pracuje. Na určení stačilo sečíst 18 členů řady, protože 1/18! je řádově 10~16, tedy další členy by se ve výsledku již neprojevily. 2.2 Symbolické derivování Symbolické derivování je poměrně snadné v Matlabu i v Ságe. V Matlabu si nejdříve definujeme symbolický objekt a pak s ním můžeme provádět symbolické operace, tedy i derivování: >> fl=sym('sin(x)*cos(x)') fl = cos(x)*sin(x) » diff(fl) ans = cos(x)"2 - sin(x)"2 » f2=sym('log(t"2)'); » diff(f2) ans = 2/t Vidíme, že pokud zadaná funkce obsahuje jen jednu proměnnou, derivuje se automaticky podle ní a výsledek je současně upraven. Pokud máme funkci více proměnných, je potřeba si zvolit, podle které chceme derivovat, pokud to není zrovna x. » f3=sym('cos(x"2)*sin(y)'); » diff(f3) ans = 40 -2*x*sin(x~2)*sin(y) » diff(f3,'y') ans = cos(x"2)*cos(y) Je taky možné počítat vyšší derivace: » diff(f2,3) ans = 4/t"3 » diff(f3,'y',2) ans = -cos(x"2)*sin(y) Derivování v Sage je podobně jednoduché, jen nejdříve musíme definovat proměnné. Máme dokonce na výběr, který příkaz pro derivování použijeme a samozřejmostí je počítání vyšších derivací. sage: x=var('x') sage: t=var('V) sage: diff(exp(x)*cos(x),x) -e"x*sin(x) + e"x*cos(x) sage: derivative(exp(x)*cos(x),x) -e"x*sin(x) + e"x*cos(x) sage: diff(log(t)*exp(x),t) e"x/t sage: diff(log(t)*exp(x),x) e"x*log(t) sage: diff(log(t)*exp(x),t,3) 2*e"x/t"3 2.3 Taylorův rozvoj Určit Taylorův rozvoj funkce, když dokážeme počítat derivace, je už snadné. Spočítali bychom derivace funkce v daném bodě, určili tak příslušné koeficienty, věřím, že by s s tímto problémem čtenář bez problémů poradil. V 41 Matlabu ale existuje hezký prográmek taylortool, ve kterém si můžete zadat, po jakou funkci chcete Taylorův rozvoj spočítat, ve kterém bodě a do jakého stupně a jako výsledek uvidíte něco jako na dalším obrázku. Taylor Series Approximation Jak by se dalo čekat, určit Taylorův rozvoj funkce v Sage je jednoduché. Stačí definovat proměnnou, funkci a pak použít funkci taylor: sage: x=var('x') sage: fl=sqrt(x+1) sage: f1.taylor(x,0,6) -21/1024*x~6 + 7/256*x"5 - 5/128*x"4 + l/16*x"3 - l/8*x"2 + l/2*x + 1 sage: f2=sqrt(x) sage: f2.taylor(x,1,6) -21/1024*(x - 1)"6 + 7/256*(x - 1)"5 - 5/128*(x - 1)"4 + l/16*(x - 1)"3 - l/8*(x - 1)"2 + l/2*x + 1/2 Jde to i bez definování funkce: sage: taylor(sin(x)*cos(x),x,pi,6) -pi - 2/15*(pi - x)"5 + 2/3*(pi - x)"3 + X Možná by stálo za to trochu prozkoumat, podle čeho se určuje pořadí členů na výstupu. 42 2.4 Numerické derivování Může se stát, že známe hodnoty funkce jen v diskrétních uzlech, nicméně potřebujeme spočítat alespoň přibližně hodnotu její derivace v nějakém bodě. Zpravidla se jedná o některý z uzlů, ale nemusí to tak být vždy. V takové situaci máme dvě možnosti: sestrojíme interpolační polynom a ten zderivujeme, nebo si vypomůžeme Taylorovým rozvojem. Oba přístupy dávají stejné výsledky, u Taylorova rozvoje dostaneme navíc odhad chyby, výsledky dosažené derivací interpolačního polynomu zase dostaneme bez velkého přemýšlení. Začneme tedy s interpolačním polynomem. Zde se hodí Lagrangeův tvar n Pn{x)=YJhk{x)1 i=0 kde li jsou Lagrangeovy fundamentální polynomy. Derivováním této formule dostaneme n i=0 Podrobnosti o přesnosti a dalších vlastnostech této formule se čtenář může dočíst v [1], my se zde blíže podíváme jen na některé speciální případy. Odvodíme formuli pro přibližný výpočet derivace v prostředním ze tří ekvi-distantních uzlů. Kvůli symetrii formulí si je označíme trochu jinak než obvykle: Xo, Xi, Xj = Xq + íh, í = ±1. Nejprve sestrojíme Lagrangeův interpolační polynom: X—i Xq X\ fi f-i fo fl _ (x-x0)(x-x1) (x - - Xi) ť2{x) - J-i--w-rJ-i + Jo 7-r;-r + {X_i — Xq){X_i — Xi) {Xq — X_i){Xq — Xi) + ^ {X ~ X_1){X - Xq) {Xi — X_i) {Xi — Xq) ' Výraz zderivujeme 43 pro Xj = xq + ih, i = ±1. Pokud dosadíme x = xq dostaneme /'(xo)«^(x0) = ^(/i (2.1) Všimněme si ještě geometrického významu této formule. Výraz (fi — f_i)/2h je směrnice sečny, která je určena body /_i) a (xi, /i). Podobně můžeme odvodit i formule pro výpočet derivace v dalších uzlových bodech x_i, xi. /'(*-0 « ^(-3/_i + 4/o-A) (2.2) í\xi) « ^(/-i-4/0 + 3/0 (2.3) Tyto formule se nazývají tříbodové. Podívejme se ještě, jak se dají stejné formule odvodit z Taylorova rozvoje. Pro uzly x±i = xq ± h dostáváme 2 o 24 Odečtením obou rovností a vydělením 2h dostáváme /,K) = ^(/i-/-0-^£l^----, (2-4) takže máme stejnou formuli jako v předchozím případě, navíc vidíme, že chyba je řádově h2, takže pokud například hodnotu h zmenšíme na polovinu, chyba klesne přibližně na čtvrtinu. Pokud bychom pro přibližné vyjádření derivace v bodě xq použili jen jeden Taylorův rozvoj, dostaneme f\xo) = lifi-fo)-^h-..., tedy výraz s chybou o řád horší než vztah předchozí. Uvedené Taylorovy rozvoje ale můžeme i sečíst, v tom případě dostaneme vztah pro přibližný výpočet druhé derivace: /"(zo) = ^(/-i " 2/o + A) " ^p±h2 (2.5) Vidíme, že chyba uvedeného vztahu je řádově opět rovna h2. 44 Kapitola 3 Integrály 3.1 Symbolické integrování Co se týče počítání symbolického počítání neurčitých i určitých integrálů, zvládají to oba námi používané softwarové prostředky bez velkých problémů. Třeba v Matlabu pro výpočet neurčitého nebo určitého integrálu používáme funkci int, rozdíl je jenom v zadaných parametrech: » f=sym('sin(x)*cos(2*x)'); » int(f) ans = cos(x) - (2*cos(x)~3)/3 » int(f,0,pi) ans = -2/3 Je možné i integrovat podle jedné z více proměnných: » g=sym('exp(-t)*sin(2*x*t)'); » int(g,'t') ans = -(sin(2*t*x) + 2*x*cos(2*t*x))/(exp(t)*(4*x"2 + 1)) » int(g,'x',0,pi/2) ans = -(cos(pi*t) - l)/(2*t*exp(t)) 45 Dají se dokonce vyjádřit integrály z funkcí, jejichž primitivní funkce sice existují, ale nedají se vyjádřit v rozumném konečném tvaru: » G=sym('exp(-x~2)'); » int(G) ans = (pi~(l/2)*erf(x))/2 » I=int(G,0,l) I = (pi~(l/2)*erf(l))/2 » eval(I) ans = 0.7468 » int(G,0,inf) ans = pi~(l/2)/2 Výrazem erf v Matlabu získáme primitivní funkci ke1, jak zjistíme z nápovědy, a Matlab umí spočítat hodnoty této funkce: >> help erf erf Error function. Y = erf(X) is the error function for each element of X. X must be real. The error function is defined as: erf(x) = 2/sqrt(pi) * integral from 0 to x of exp(-t"2) dt. See also erfc, erfcx, erfinv, erfcinv. Overloaded methods: sym/erf Reference page in Help browser doc erf » erf(l) ans = 0.8427 46 V Sage je to velmi podobné: sage: x=var('x') sage: integral(sin(x)*x,x) -x*cos(x) + sin(x) sage: integral(sin(x)*x,x,0,pi) Pi sage: integrál(sin(x)*x,x,0,pi/2) 1 sage: t=var('ť) sage: integrál(sin(x),t) t*sin(x) sage: integrál(1/sqrt(x),x,0,1) 2 sage: integral(l/x~2,x,0,1) Traceback (click to tne left of this block for traceback) ValueError: Integrál is divergent. I primitivní funkce k e x se vyjadřuje podobně: sage: integrál(exp(-x"2),x) l/2*sqrt(pi)*erf(x) sage: integrál(exp(-x"2),x,0,oo) l/2*sqrt(pi) sage: IO=integral(exp(-x"2),x,0,1);10 l/2*sqrt(pi)*erf(1) sage: 10.N() 0.746824132812427 3.2 Numerické integrování Při numerickém integrování se podobně jako při derivování nesnažíme určit primitivní funkci z funkčních hodnot v zadaných diskrétních uzlech, ale snažíme se na základě těchto dat určit přibližně určitý integrál přes nějaký interval. 47 Základem získaných formulí jsou opět interpolační polynomy. Při integrování Lagrangeova tvaru dostáváme: b b b / f{x)dx ~ / Pn(x)dx = J2fi kix)dx = J^fiA 1=0 „ 1=0 b pro Al = J li{x)dx. Podobný typ výrazů dostáváme také při odvozování Rie- a mannova integrálu, kde aproximaci určitého integrálu dostáváme jako součet funkčních hodnot násobených koeficienty zahrnujícími délky subintervalů při dělení původního intervalu na menší dílky. Těmto výrazům, tedy ° i=0 říkáme kvadraturní formule. Kromě odhadu chyby nás u kvadraturních formulí zajímá také její stupeň přesnosti, který udává, pro polynomy jakých stupňů je tato formule přesná. Přesněji, pokud je stupeň přesnosti formule roven n pak je tato formule je přesná pro polynomy do stupně n včetně. Dá se snadno ukázat (viz [1]), že pro n + 1 uzlů je stupeň přesnosti roven maximálně 2n + 1. Na druhou stranu je zřejmé, že pokud kvadraturní formuli získáme postupem ukázaným výše, tedy integrací interpolačního polynomu, musí být stupeň přesnosti roven alespoň n. Konstrukce kvadraturních formulí je možné ještě poněkud zobecnit tím, že do integrálu přidáme ještě tzv. váhovou funkci. Do ní můžeme zahrnout například singularity nebo společnou část pro nějakou třídu funkcí. V tom případě kvadraturní formule nahrazují integrál J^fiA ~ / w(x(-f(x)dx. i=0 Koeficienty kvadraturní formule pak získáme ze vztahu b Al = J w(x) ■ li{x)dx. a Výhodou takového postupu je, že váhová funkce není zahrnuta ve funkčních hodnotách funkce /. V tomto textu se ale budeme zabývat pouze jednoduchou situací s váhovou funkcí rovnou jedné. 48 Odvodíme si některé nej používanější formule. Situací s jedním uzlem se zabývat nebudeme, tak je až příliš jednoduchá. Pro dva uzly většinou máme a = xq < x\ = b. V tomto případě má interpolační polynom tvar například ft(l) = /(o) + MzM(l_o). b — a Jeho integrováním přes interval [a, b] vyjde kvadratur ní formule a což je v podstatě plocha lichoběžníka (postaveného na bok),jehož strany mají délky rovny funkčním hodnotám funkce / v krajních bodech intervalu a jehož výška je rovna b — a. Proto se taky této formuli říká lichoběžníkové pravidlo. V případě, že pudeme funkci aproximovat polynomem druhého stupně, tedy parabolou, dostaneme pro uzly xq = a, x\ = x2 = b formuli / f(x)dx ~ ^ (f(a) + 4/(^) + /(&)) , která se nazývá Simpsonovo pravidlo. 49 (a+b)/2 b Z těchto dvou formulí se odvozují patrně nejpoužívanější kvadraturní formule. Jejich myšlenka je jednoduchá: pokud máme více uzlů, pak nám rozdělují interval [a, b] na více subintervalů, na každém tomto subintervalu použijeme lichoběžníkové pravidlo nebo Simpsonovo pravidlo a výsledek pak sečteme. Tím dostaneme složené lichoběžníkové nebo složené Simponovo pravidlo. U složeného lichoběžníkovécho pravidla v podstatě počítáme integrál z aproximace funkce lomenou čárou, respektive lineárním splajnem. U Složeného Simpsonova pravidla zase vyžadujeme, aby celkový počet uzlů byl lichý, na což by jistě čtenář přišel sám. pro ekvidistantní uzly kde vzdálenost sousedních dvou je rovna h dostáváme složené lichoběžníkové pravidlo ve tvaru r h j f(x)dx « - (/o + 2/i + 2/2 + • • • + 2fn_x + /„), a kde fi = /(x,). 50 51 Kapitola 4 Rady Řad už jsme částečné dotkli, když jsme se zmínili o Taylorově rozvoji v souvislosti s derivacemi. Teď si zkusíme ukázat, jak se dají počítat některé součty řad a nesmíme zapomenout taktéž na řady Fourierovy. 4.1 Symbolické součty řad Jak se dá čekat, v Matlabu je se symbolickým sčítáním řad potíž. Můžeme zde maximálně odhadnout číselný součet řady, jestliže připočítáváme k součtu členy, které jsou už z numerického hlediska zanedbatelné, jak jsme to ukázali při výpočtu Eulerova čísla. Proto se v dalším výkladu soustředíme na Sage. Nejprve si zkusíme sečíst jednoduchou geometrickou řadu: var('x','k','n') sum(x~k,k,0,oo) Traceback (click to tne left of this block for traceback) Is abs(x)-l positive, negative, or zero? Získali jsme pouze chybové hlášení. Chybí totiž bližší informace o proměnné x a víme, že uvedená řada má součet jen v případě, že |x| < 1. Proto tento požadavek zadáme do Sage: sage: assume(abs(x)<1) sage: sum(x~k,k,0,oo) -l/(x - 1) 52 Vidíme, že tohle funguje a můžeme zkoušet další součty: sage: sum(k*x~k,k,1,oo) x/(x~2 - 2*x + 1) sage: sum(l/k~4, k, 1, oo) l/90*pi~4 sage: sum(x"k/factorial(k), k, 0, oo) e"x Zatím to vypadá všechno bezproblémově. V některých situacích si ovšem ani Sage neporadí. Ukážeme si to na poměrně jednoduchém příkladu. Nejprve si spočítáme Taylorův rozvoj funkce \/l — x. sage: f=sqrt(l-x) sage: f.taylor(x,0,6) -21/1024*x"6 - 7/256*x"5 - 5/128*x"4 - l/16*x"3 - l/8*x"2 - l/2*x + 1 Čísla, která se zde objevují, jsou binomické koeficienty, které jsou pro reálný argument a definované vztahem ía\ a(a-l)---(a-k + l) \k) = k\ V našem případě jsou to koeficienty pro a = 1/2 (až na znaménko), jak snadno ověříme: [binomial(l/2,k) for k in range(7)] [1, 1/2, -1/8, 1/16, -5/128, 7/256, -21/1024] Zkusíme tedy řadu zpětně sečíst. Problémem bude trochu první člen, ale ten by se měl projevit jen posunem výsledku o konstantu: sum(x"k*(-l)~k*binomial(l/2,k),k,0,oo) Traceback (click to tne left of this block for traceback) TypeError: Either m or x-m must be an integer Vidíme, že na tomto Sage ztroskotal. 53 4.2 Fourierovy řady Spočítat koeficienty Fourierovy řady na základě toho, co už známe, jistě nebude pro čtenáře problém. Například pokud bychom chtěli v Matlabu určit prvních pár členů Fourierovy řady pro funkci x2, můžeme postupovat třeba takto: » f=sym('x~2'); » aO=int(f,-pi,pi)/(2*pi) aO = pi"2/3 >> al=int('cos(x)'*f,-pi,pi)/pi al = -4 » a2=int('cos(2*x)'*f,-pi,pi)/pi a2 = 1 >> a3=int('cos(3*x)'*f,-pi,pi)/pi a3 = -4/9 Koeficienty u sinových členů jsou nulové, neboť x2 je sudá funkce. Sečteme tedy první čtyři členy řady a výsledek porovnáme s původní funkcí: » x=-pi:0.01:pi; » S=eval(a0+al*cos(x)+a2*cos(2*x)+a3*cos(3*x)); » plot(x,S,x,x. "2,' —') Pro výpočet hodnot součtů v bodech daných vektorem x musíme použít funkci e val, protože koeficienty aO,.. . ,a3 jsou symbolické objekty. 54 Pokud chceme určit Fourierovu řadu v Sage, je potřeba si funkci definovat jen na příslušném intervalu, tedy např. na [—7r,7r]: sage: x=var('x') sage: f=Piecewise( [ [(-pi,pi),x~2]]) Stejné koeficienty jako v Matlabu určíme pomocí příkazů sage: f.fourier_series_cosine_coefficient(0,pi) 2/3*pi~2 sage: f.fourier_series_cosine_coefficient(1,pi) -4 sage: f.fourier_series_cosine_coefficient(2,pi) 1 sage: f.fourier_series_cosine_coefficient(3,pi) -4/9 První vstupní parametr říká, který koeficient určujeme, druhý je polovina intervalu, na němž řadu počítáme. Částečný součet řady ale můžeme získat i přímo a porovnat ho s původní funkcí: sage: g=f.fourier_series_partial_sum(3,pi) sage: pl=plot( [g,x"2],-pi,pi) 55 56 Kapitola 5 Řešení nelineárních rovnic 5.1 Motivační úloha Začneme geometrickou úlohou s jednoduchým zadáním. Mějme kruh v rovině o zadaném poloměru r. Máme z hranice kruhu udělat kruhový oblouk o neznámém poloměru x tak, aby část ohraničená oběma křivkami měla plochu, která je rovna polovině plochy původního kruhu. Situaci si můžeme znázornit na obrázku: 57 Tato úloha je matematickým vyjádřením úlohy ze zemědělství. Sedlák má kruhový pozemek, na kterém se pase koza. Protože sedlák chce, aby jí tráva na pozemku vystačila na dva dny, uváže ji ke kůlu na okraji pozemku tak dlouhým provazem, aby za první den spásla polovinu trávy. Druhý den ji nechá k dispozici celý pozemek, kde může spást zbylou trávu. Pro řešení úlohy si do obrázku doplníme několik údajů. 58 Spojíme střed kružnice i oblouku s průsečíky na kružnici. Průvodiče vycházející od kraje kružnice svírají úhel a, průvodiče vycházející ze středu kružnice tedy svírají úhel 2a. Oblast, která nás zajímá, se skládá ze dvou kruhových úsečí, první z nich má obsah 5*1 a je dána poloměrem x a úhlem a, druhá má obsah 5*2, poloměr r a úhel 2tt — 2a. Označíme-li obsah celého kruhu S, dostáváme Si + S2 = -S. Obsah kruhové úseče určíme tak, že od obsahu kruhové výseče odečte obsah trojúhelníka o stranách rovných poloměru výsečejež svírají úhel, jež určuje výseč. Tedy Si = ^x2(a - sin a), S2 = ^2((2tt - 2a) - sin(27r - 2a)) Ještě potřebujeme vztah mezi r a x. K tomu můžeme použít pravoúhlý trojúhelník, který dostaneme, když spojíme střed kružnice se středem tětivy délky x podle následujícího obrázku: 59 V uvedeném obrázku platí x 12 a - = cos —, r 2 tedy a x = 2r cos —. 2 Můžeme tedy namísto x hledat úhel a, z předchozího vztahu pak x snadno vypočítáme. Nyní už máme všechny potřebné vztahy pohromadě a použitím goniometrických vzorečků (toto cvičení jistě čtenář snadno zvládne) získáme rovnici a = tan a 2 cos a Přesné řešení této rovnice ovšem získat neumíme. Nezbývá, než použít nějakou numerickou metodu, kterýmžto se budou věnovat následující řádky. 60 5.2 Základní metody řešení nelineární rovnice Budeme se zabývat hledáním řešení rovnice f{x) = 0, xel=[a,b] (5.1) pro spojitou funkci /. To že řešení na uvedeném intervalu existuje nám zaručí například podmínka f (a) ■ f(b) < 0, tedy v krajních bodech intervalu I má funkce / opačná znaménka. Řešení ovšem může být více. Označme řešení rovnice jako £, nazýváme jej také kořenem funkce /. Všechny metody a tvrzení, které zde budou vedeny, může čtenář nalézt v detailním znění ve skriptech [1]. 5.2.1 Půlení intervalu Nejjednodušší metodou pro nalezení řešení je metoda půlení intervalu, zvané též metoda bisekce. Její myšlenka je velmi jednoduchá: jestliže máme splněnu podmínku f (a) ■ f(b) < 0, rozpůlíme interval [a, b] a jako nový interval vezmeme ten, pro který platí, že v jeho krajních bodech má funkce / opět opačná znaménka, tedy tento poloviční interval opět obsahuje řešení rovnice. Pokračujeme s půlením tak dlouho, dokud nedostaneme interval dostatečné malý (stále obsahující řešení), tedy máme přibližné řešení s požadovanou přesností. Jako přibližné řešení £ stanovíme střed konečného intervalu. U metody půlení intervalu je možné dokonce předem určit, kolik iterací budeme potřebovat. Jestliže například požadujeme, aby chyba přibližného řešení byla menší než ó, dostáváme pro počet kroků k nerovnost h-a<8 2k+i odkud logaritmování při základu 2 dostáváme b — a k > log2 —. Metoda půlení intervalu má jednu základní výhodu, že totiž za uvedených předpokladů vždy konverguje. Její nevýhodou je, že konverguje poměrně pomalu. A jelikož jsme v podstatě vyčerpali všechny důležité informace o ní, budeme se věnovat metodám dalším. 61 5.2.2 Prostá iterační metoda U této metody se zaměříme na rovnici v jiném tvaru: x = g(x), hledáme tedy hodnotu £, kterou funkce g zobrazí samu na sebe. Proto se bod £ nazývá pevný bod funkce g. Pro funkci g samozřejmě požadujeme spojitost na intervalu I = [a,b]. Pro existenci pevného bodu na tomto intervalu pak postačuje, aby se interval I zobrazil sám do sebe, tedy Vx £ I : g (x) £ I. Tato podmínka ovšem nezaručí jednoznačnost pevného bodu. Pro ni potřebujeme, aby funkce g byla kontrakcí na intervalu I, to jest musí existovat konstanta L, 0 < L < 1, že pro každé x,y £ I platí \g{x) -g(y)\ > while abs(x-cos(x))>presnost, x=cos(x), end x = 0.5403 x = 0.8576 x = 0.6543 x = 0.7935 x = 0.7014 x = 0.7392 x = 0.7390 x = 0.7391 » Na dosažení požadované přesnosti je potřeba něco přes dvacet iterací. Pokud bychom chtěli použít prostou iterační metodu na řešení úvodního příkladu pro rovnici 7T x — tčin x ^ 2 cos x tedy pro funkci gix) = tanx— 2cqSXi můžeme si vypomoci obrázkem, abychom mohli lépe odhadnout interval, kde se pevný bod nachází. Při tom využijeme toho, že musí ležet někde v intervalu [7r/2,7r], který je třeba o něco zkrátit, jelikož v levém krajním bodě jde funkce g do nekonečna. 64 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Zjistíme, že pevný bod leží v intervalu [1.8, 2]. Na tomto intervalu ale použít prostou iterační metodu není možné, neboť absolutní hodnota derivace funkce (/je zde větší než jedna. 5.2.3 Newtonova metoda Vrátíme se nyní zpět k rovnici f(x) = 0. Na řešení této rovnice můžeme nahlížet jakožto na průsečík grafu funkce / s osou x. Myšlenka Newtonovy metody je prostá. Zvolíme počáteční iteraci xq a další iterací je průsečík tečny ke grafu s osou x. V podstatě hledáme přibližné řešení na lineární aproximaci funkce /. Tento postup opakujeme tak dlouho, dokud nemáme přibližné řešení s dostatečnou přesností. Metodu si můžeme dobře ilustrovat graficky: 65 Je zřejmé, že v tomto případě kromě spojitosti funkce / potřebujeme také spojitost její derivace. Vyjádříme-li si z rovnice tečny k / v bodě [xk-f(xk)] její průsečík z osou x, dostáváme iterační vztah Xk+1-Xk~iwr Stejný vztah bychom obdrželi použitím Taylorova rozvoje, v němž bychom zanedbali členy od druhé derivace včetně. Máme vlastně opět prostou iterační metodu, tentokrát pro funkci g ve speciálním tvaru 9{x)=x~wv můžeme tedy opět zkoumat podmínky, za nichž je tato funkce kontrakcí. Obecné vyjádření je ovšem složité, platí ale tvrzení, že pokud je derivace funkce / nenulová na intervalu obsahujícím řešení rovnice, pak existuje okolí tohoto řešení, na němž je funkce g kontrakcí, tedy pro libovolnou počáteční aproximaci z tohoto okolí Newtonova metoda konverguje. Toto kriterium vypadá v praxi jen špatně použitelné, pokud se týče hledání počáteční aproximace. Naštěstí existují jiné podmínky, snadněji ověřitelné, umožňující její volbu, jsou to například Fourierovy podmínky: 1. Funkce / má na intervalu I spojité druhé derivace a v krajních bodech I opačná znaménka. 66 2. Derivace funkce / je nenulová na celém intervalu, tedy ani nemění znaménko. 3. Druhá derivace funkce / nemění na I znaménko. Při splnění těchto podmínek Newtonova metoda konverguje na /, a to dokonce monotónně, pakliže jako počáteční iteraci zvolíme ten krajní bod intervalu, v němž má funkční hodnota stejné znaménko jako je znaménko druhé derivace. Uvedené podmínky vlastně říkají, že funkce / je na intervalu I ryze monotónní (rostoucí nebo klesající) a to buď konvexní nebo konkávni. Celkem jsou tedy čtyři případy, následující obrázek ukazuje situaci, kdy je / rostoucí a konkávni. V tom případě volíme jako počáteční aproximaci ten krajní bod, v němž je funkční hodnota záporná. Další tři možnosti by si čtenář jistě dokázal snadno představit. Příklad 3. Zkusme požít Newtonovu metodu na úvodní příklad, tedy na rovnici 7T x — tčin x . 2 cosx Nejprve všechny členy převedeme na jednu stranu, abychom dostali funkce /: x — tan x H--= 0 2 cosx 67 Pak zapojíme Matlab, přičemž použijeme toho, že z už známe interval, kde řešení leží. » f=sym('x-tan(x)+pi/(2*cos(x))') f = x - tan(x) + pi/(2*cos(x)) » df=diff(f) df = (pi*sin(x))/(2*cos(x)"2) - tan(x)"2 » x=1.8; » x=x-eval(f)/eval(df) x = 1.8735 » x=x-eval(f)/eval(df) x = 1.9028 » x=x-eval(f)/eval(df) x = 1.9057 » x=x-eval(f)/eval(df) x = 1.9057 » Vidíme, že aproximaci řešení na čtyři desetinná místa jsem nalezli velmi rychle. Mohli bychom samozřejmě nastavit větší přesnost a počítat dál, ale i tak bychom dosáhli během několika dalších iterací nej lepší možné aproximace v rámci počítačové přesnosti. Vypadá to tedy, že Newtonova metoda je velmi rychlá, i při jejím použití ovšem nesmíme zapomínat na jistou míru bedlivost, zejména co se týče volby počáteční aproximace. To si můžeme dobře demonstrovat na jednoduchém příkladě - totiž na funkci fix) = arctanx. Zvolíme-li jako počáteční aproximaci hodnotu 1.5, iterační posloupnost nekonverguje, jak naznačuje obrázek. Volbou xq = 1 by vše bylo v pořádku. 68 1.5 r 1 " 0.5 - -0.5 - -1 - -1.5 -6 -4 -2 4 5.2.4 Metoda sečen a regula falsi Další metody, o nichž si něco řekneme, jsou podobné. Jejich hlavní myšlenka spočívá v náhradě derivace v bodě xk. To je nutné v případech, že derivaci nelze vyjádřit, například když hodnoty funkce / nejsou dány nějakým vzorcem, ale jsou to třeba výsledky fyzikálního měření. Derivaci pak nahradíme poměrnou diferencí, potřebujeme ovšem funkční hodnoty ve dvou bodech: Použijeme-li tuto náhradu v iteračním vztahu pro Newtonovu metodu, dostáváme Metoda daná tímto iteračním vztahem se nazývá metoda sečen. Tento název je odvozen z toho, že tečná z Newtonovy metody je nahrazena sečnou protínající graf funkce / ve dvou bodech. Z iteračního vztahu metody sečen také plyne, že na jejím začátku potřebujeme dvě iterace Také v každém kroku používáme dvě funkční hodnoty, stačí ale počítat jen jednu a tu druhou si pamatovat od minule. Vyzkoušejme si opět metodu sečen v Matlabu na úvodním příkladu: f\xk) f{xk) ~ f{xk-i) %k %k — l 69 » f=inline('x-tan(x)+pi/(2*cos(x))'); » x0=1.8; » xl=2; » fO=f(xO) fO = -0.8274 » fl=f(xl) fl = 0.4104 » x2=xl-fl*(xl-xO)/(f1-fO) x2 = 1.9337 » f2=f(x2) f2 = 0.1422 » x3=x2-f2*(x2-xl)/(f2-f1) x3 = 1.8985 » f3=f(x3) f3 = -0.0401 » X4=x3-f3*(x3-x2)/(f3-f2) x4 = 1.9063 » f4=f(x4) f4 = 0.0031 » X5=x4-f4*(x4-x3)/(f4-f3) x5 = 1.9057 » f5=f(x5) f5 = 6.1196e-05 » X6=x5-f5*(x5-x4)/(f5-f4) x6 = 1.9057 » f6=f(x6) 70 f6 = -9.5082e-08 » Počítat tímto způsobem je jistě poněkud nepraktické, zde to slouží jen jako ukázka. Čtenář by jistě zvládl celý problém řešit elegantněji, například pomocí cyklu. Grafické znázornění metody sečen ukazuje následující obrázek: 0.51- -1'-'-1-1-1-1-1-1-1-1-1 1.8 1.82 1.84 1.86 1.88 1.9 1.92 1.94 1.96 1.98 2 Není ovšem potřeba volit počáteční aproximace tak, aby funkční hodnoty měly opačná znaménka. Naopak, pokud budou blízko sebe, bude příslušná sečna bližší tečně z Newtonovy metody: 71 0.5 r -1 '-'-'-'-'-'-'-'-'-'-1 1.8 1.82 1.84 1.86 1.88 1.9 1.92 1.94 1.96 1.98 2 U metody regula falsi naopak počáteční aproximace z opačnými funkčními hodnotami požadujeme a tento požadavek pak platí pro všechny další iterace. Jedná se o podobnou myšlenku jako u metody půlení intervalu, iterace ale počítáme podobně jako u metody sečen: r f \ Xfc Xs Xk+1 = Xk~f[Xk)' f{Xk)-f{Xsy kde s je největší index menší než k takový, že f(xs) ■ fixk) < 0. Pro tuto metodu je typické, že pokud je funkce / na intervalu I konvexní nebo konkávni, jeden z bodů počátečních aproximací je tzv. fixní, tedy index s je pořád stejný, jak vidíme i z následujícího obrázku: 72 -1.2'-1-1-1-1-1-1-1-1-1-1 0123456789 10 5.3 Rychlost konvergence Z uvedených příkladů vidíme, že rychlost konvergence se liší u různých metod. Proto si definujme řád konvergence iterační metody, který udává rychlost její konvergence: Nechť posloupnost (x^ je iterační posloupnost získaná nějakou metodou, Xk —t* £. Označme chybu v k-tém kroku, t.j. = Xk — £. Řekneme, že iterační metoda je řádu p > 1, jestliže platí lim = C^0. Rád konvergence vlastně říká, jak musíme umocnit chybu v k-tém kroku, abychom limitně dostali chybu v kroku k + 1 (až na násobení konstatnou). Je jasné, že čím vyšší je řád metody, tím je konvergence rychlejší. Například pro p = 1 (říkáme také, že metoda je lineárního řádu) je přechod mezi jednotlivými chybami jen přibližně násobení konstantou (ta tedy musí být menší než 1), kdežto pro metodu řádu 2 (kvadratického řádu) platí, že pokud je chyba v některém kroku v desetinách, v dalším bude v setinách, v následujícím už v desetitisícinách atd. Pro metody, které jsem si ukazovali je řád znám. Bisekce je lineární, Newtonova metoda kvadratická, metoda sečen má řád roven zlatému řezu 1+2V^, metoda regula falši je opět jen lineární. Pro prostou iterační metodu s funkcí 73 g platí, že pokud g má spojité derivace do řádu p včetně, £ je pevný bod funkce a g'(£) = ... = = 0, <7^(£) 7^ 0, pak metoda je řádu p. Je známa metoda, pomocí níž se dá rychlost konvergence urychlit. Jedná se o Aitkenovu ó2 metodu, která je založena na tom, že pokud by iterační posloupnost konvergovala přesně lineárně, je možné pomocí jednoduchého výpočtu ze tří členů určit limitu. Mějme tedy posloupnost (xk). Z ní sestrojíme posloupnost (xk) s použitím vztahu ~ _ {xk+i — Xk) Xk Xk ~ ' . xk+2 — ^xk+l + Xk Platí, že pokud původní posloupnost konverguje k £, konverguje k této hodnotě i nová posloupnost. Dále pokud je řád konvergence původní posloupnosti lineární, je nová posloupnost řádu alespoň kvadratického, je li řád původní posloupnosti p > 1, je řád nové posloupnosti 2p—l. Tento postup tedy urychlí i Newtonovu metodu. 5.3.1 Steffensenova metoda Aitkenovu metodu nelze použít tak, jak je uvedena, jelikož bychom potřebovali znát celou původní posloupnost, dokázali bychom tedy odhadnout i její limitu, nebylo by tedy potřeba počítat ji znovu, byť rychleji. Nicméně kombinací Aitkenovy metody s prostou iterační metodou získáme Steffensenovu metodu, která zpravidla dává rychle konvergující posloupnost. Myšlenka Ste-ffensenovy metody je prostá: určíme tři iterace prostou iterační metodou (počáteční iteraci a dvě další) a pak použijeme Aitkenovu metodu na urychlení výpočtu. Pak určíme další tři iterace prostou iterační metodou a zase provedeme urychlení, což opakujeme, dokud nedosáhneme požadované přesnosti. Celý postup můžeme zapsat třeba takto (pro počáteční iteraci xo): /x /x (yo-xo)2 yo = 9{x0), z0 = g(y0), x1 = x0 zo ~ 2y0 + x0 a obecný krok je pak {Vk ~ Xk)2 yk = g(xk), zk = g(yk), xk+i = xk zk ~ 2yk + xk Steffensenova metoda konverguje pro počáteční iteraci z nějakého okolí pevného bodu £ funkce g, jestliže g'(£) 7^ 1. 74 Příklad 4. Vyzkoušíme Steffensenovu metodu na řešení úvodního příkladu pro funkci g{x) = tana" — 2Zsx' Ja^° počáteční iteraci zvolíme xq = 1.8 a výpočet provedeme v Matlabu: » g=inline('tan(x)-pi/(2*cos(x))'); » x=1.8; » y=g(x),z=g(y) y = 2.6274 z = 1.2392 » x=x-(x-y)"2/(z-2*y+x) x = 2.1090 » y=g(x),z=g(y) y = 1.3894 z = -3.2542 » x=x-(x-y)"2/(z-2*y+x) x = 2.2410 » y=g(x),z=g(y) y = 1.2672 z = -2.0623 » x=x-(x-y)"2/(z-2*y+x) x = 2.6435 » y=g(x),z=g(y) y = 1.2442 z = -1.9440 » x=x-(x-y)~2/(z-2*y+x) x = 3.7379 75 » Vidíme, že nevhodně zvolená počáteční iterace ke konvergenci nevede. Zkusíme počáteční iteraci poněkud zpřesnit: » x=1.85; » y=g(x),z=g(y) y = 2.2117 z = 1.2865 » x=x-(x-y)"2/(z-2*y+x) x = 1.9517 » y=g(x),z=g(y) y = 1.7283 z = 3.7177 » x=x-(x-y)"2/(z-2*y+x) x = 1.9291 » y=g(x),z=g(y) y = 1.8087 z = 2.5417 » x=x-(x-y)"2/(z-2*y+x) x = 1.9121 » y=g(x),z=g(y) y = 1.8775 z = 2.0449 » x=x-(x-y)"2/(z-2*y+x) x = 1.9062 76 » y=g(x),z=g(y) y = 1.9034 z = 1.9159 » x=x-(x-y)"2/(z-2*y+x) x = 1.9057 » y=g(x),z=g(y) y = 1.9057 z = 1.9058 » x=x-(x-y)"2/(z-2*y+x) x = 1.9057 » 5.4 Řešení systému nelineárních rovnic Mějme systém nelineárních rovnic /i(xi,... ,xn) = 0 fn{xi, ■ ■ ■ ,x„) = 0 Můžeme jej stručně zapsat ve tvaru F(x) = 0, kde F i x bereme jako vektory. Nejpožívanější metodou je Newtonova metoda, kterou dostaneme tak, že vezmeme Taylorův rozvoj funkce více proměnných do prvního řádu. Označme JF matici derivací vektoru funkcí F, tedy 77 Pak Newtonovu metodu můžeme vyjádřit podobným způsobem jako je to pro funkci jedné proměnné: xfc+i = xfc - J^x(xk) • F(xk). Protože ale práci s funkcemi více proměnných jsme zatím v Matlabu ani v Sage nezkoušeli, odložíme prozatím praktické vyzkoušení Newtonovy metody. 78 Kapitola 6 Funkce více proměnných 6.1 Zobrazovnání grafů První věcí, kterou si ukážeme, či spíše zopakujeme, pro funkce více proměnných je nakreslení jejich grafu. To samozřejmě dokážeme jen pro funkce dvou proměnných. V Sage můžeme postupovat třeba takto: sage: x,y = var('x,y') sage: plot3d(sin(x+y)*exp(-(x~2+y~2)/2) , (-pi,pi) , (-pi,p|i)) Tím dostaneme následující obrázek: 3.1 -3.1 Získat podobný obrázek v Matlabu dá ovšem o něco víc práce. Nejprve musíme definovat síť bodů, v nichž se má funkce vyčíslit. K tomu se hodí matlabovská funkce meshgrid. Ukážeme její použití i výsledek, který dává: 79 >> x=linspace(-pi, pi,31); >> y=linspace(-pi, pi,31); » [X,Y]=meshgrid( x,y); » X(l:5,l:5) ans = -3.1416 -2.9322 -2. 7227 -2 5133 -2 3038 -3.1416 -2.9322 -2. 7227 -2 5133 -2 3038 -3.1416 -2.9322 -2. 7227 -2 5133 -2 3038 -3.1416 -2.9322 -2. 7227 -2 5133 -2 3038 -3.1416 -2.9322 -2. 7227 -2 5133 -2 3038 » Y(l:5,l:5) ans = -3.1416 -3.1416 -3. 1416 -3 1416 -3 1416 -2.9322 -2.9322 -2. 9322 -2 9322 -2 9322 -2.7227 -2.7227 -2. 7227 -2 7227 -2 7227 -2.5133 -2.5133 -2. 5133 -2 5133 -2 5133 -2.3038 -2.3038 -2. 3038 -2 3038 -2 3038 » size(X) ans = 31 31 » Vidíme, že ve výsledku se opakují vstupní vektory, v jedné výstupní matici po sloupcích, ve druhé po řádcích. Tyto výstupní matice se pak dají použít pro výpočet funkčních hodnot a výsledek si pak zobrazíme: » f=sin(X+Y).*exp(-(X."2+Y."2)/2); >> mesh(x,y,f) » 80 -4 Graf funkce je zobrazen z jiné strany, pomocí ovládacích prostředků Matlabu si jej ovšem můžeme natočit do stejné polohy jako je to na předchozím obrázku. 4 -4 Kromě funkce mesh můžeme také pro zobrazení použít funkci surf, která graf funkce / zobrazuje pomocí plošek: 81 Další možností je zobrazit vrstevnice grafu, k tomu slouží funkce contour: 3^ -3-2-1 0 1 2 3 Můžeme samozřejmě graf funkce / zobrazit na jemnější síti, abychom dostali podobný výsledek jako u Sage. Je také možné použít jiné barevné zobrazení, k tomu slouží funkce colormap, a nastavit přechody mezi jednotlivými ploškami. Další podrobnosti je možné získat v dokumentaci: >> x=linspace(-pi,pi,51); >> y=linspace(-pi,pi,51); 82 » [X,Y]=meshgrid(x,y); » f=sin(X+Y).*exp(-(X.~2+Y.~2)/2); >> surf(x,y,f) >> colormap(cool) >> shading interp » 4 -4 6.2 Výpočet parciáních derivací Výpočet parciálních derivací jev obou našich používaných softwarech velmi jednoduchý. V Sage použijeme následující postup: sage: x,y=var('x,y') sage: f=sin(x+y)*exp(x-y) sage: diff(f,x) e" (x - - y)*sin(x + y) + e"(x - y)*cos(x + y) sage: diff(f,y) -e" (x - y)*sin(x + y) + e"(x - y)*cos(x + y) sage: V Matlabu je postup obdobný: 83 >> f=sym('sin(x+y) f = *exp(x-y)') J. exp(x - y)*sin(x + y) » diff(f,'x') ans = exp(x - y)*cos(x + y) + exp(x - y)*sin(x + y) » diff(f,'y') ans = exp(x - y)*cos(x + y) - exp(x - y)*sin(x + y) » Příklad 5. Nyní si můžeme v Matlabu ukázat, jak bychom řešili systém dvou nelineárních rovnic pomocí Newtonovy metody. Uvažujme rovnice fi(x,y) = x2 — 2x + y2 = 0 h{x,y) = x-y2-l = 0 Hledáme vlastně průsečík kružnice a paraboly, i 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 84 Najít v Matlabu přibližné řešení za pomocí Newtonovy metody můžeme třeba takto: » fl=sym('x~2-2*x+y~2'); » f2=sym('x-y~2-l'); » flx=diff(f1,'x') fix = 2*x - 2 » fly=diff(f1,'y') fly = 2*y » f2x=diff(f2,'x') f2x = 1 » f2y=diff(f2,'y') f2y = -2*y » J=[flx,fIy;f2x,f2y] J = [ 2*x - 2, 2*y] [ 1, -2*y] » F=[fl;f2] F = x"2 - 2*x + y"2 - y~2 + x - 1 >> x=2;y=l; %pocatecni aproximace >> xy= [x;y]-inv(eval(J))*eval(F) xy = 1.6667 0.8333 » x=xy(l);y=xy(2); >> xy=[x;y]-inv(eval(J))*eval(F) xy = 1.6190 0.7881 » x=xy(l);y=xy(2); >> xy= [x;y]-inv(eval(J))*eval(F) xy = 85 1.6180 0.7862 » x=xy(l);y=xy(2); >> xy=[x;y]-inv(eval(J))*eval(F) xy = 1.6180 0.7862 » Newtonova metoda i v tomto případě konverguje velmi rychle. Zde by samozřejmě bylo možné najít analytické vyjádření řešení, neboť dosazením druhé rovnice do první bychom dostali kvadratickou rovnici. Čtenář tak může učinit pro jeho porovnání s výše získaným přibližným řešením. 6.3 Integrování ve více proměnných Určit primitivní funkci podle jedné či druhé proměnné s použitím Sage nebo Matlabu je jednoduché: sage: x,y=var('x,y') sage: f=sin(x+y)*cos(x-y)"2 sage: fl=integral(f,x);f1 -l/12*cos(-3*x + y) + l/4*cos(-x + 3*y) sage: f2=integral(f1,y);f2 -l/12*sin(-3*x + y) + l/12*sin(-x + 3*y) sage: >> f=sym('sin(x+y)*cos(x- -y)"2'); » fl=int(f,'y') fl = cos(y - 3*x)/4 - cos(3*y - x)/12 - cos(x + y)/2 » f2=int(f1,'x') f2 = sin(3*x - y)/12 - sin(x - - 3*y)/12 - sin(x + y)/2 » l/2*cos(x + 3 - l/2*sin(x + ) y) 86 Pak také není problémem spočítat určitý integrál přes obdélníkovou oblast, například pro uvedenou funkci přes [0, tt] x [0,7r/2], přitom by nemělo záležet na pořadí integrování: ságe: Il=integral(f,x,0,pi);II -l/2*cos(3*y) + 7/6*cos(y) sage: Ix=integral(f,x,0,pi);Ix -l/2*cos(3*y) + 7/6*cos(y) sage: II=integral(Ix,y,0,pi/2);II 4/3 sage: » Iy=int(f,'y',0,pi/2) iy = (4*cos(x))/3 + sin(x)/3 - cos(x)"3 + cos(x)~2*sin(x) » II=int(Iy,'x',0,pi) II = 4/3 » Jestliže ale potřebujeme spočítat integrál přes jinou oblast, nezbývá, než si meze oblasti vyjádřit ručně, neboť zdá se, že integrovat přes obecné oblasti zatím námi používaný software neumí. Takže pokus bychom potřebovali určit integrál z funkce f(x,y) = ix + y)2 přes kruh se středem v počátku a poloměrem 2, musíme si pro dané y určit, v jakých mezích se pohybuje x. To je celkem snadné meze jsou [—\/A — y2, \/A — y2]. Protože ale v mezích je obsaženo ometení pro y, musíme ho zadat i do výpočtu: sage: f=(x+y)"2 sage: assume(abs(y)<=2) sage: Il=integral(f,x,-sqrt(4-y"2),sqrt(4-y"2));I1 4/3*sqrt(-y"2 + 4)*(y"2 + 2) sage: I2=integral(II,y,-2,2);12 8*pi sage: U jednoduchých oblastí, jako je například kruh či elipsa, si můžeme nechat hranice spočítat: 87 ságe: S=solve(x~2+y~2-4,x);S [x == -sqrt(-y~2 +4), x == sqrt(-y"2 + 4)] ságe: Il=integral (f ,x,S [0] .right() ,S[1] .rightO) ;I1 4/3*sqrt(-y"2 + 4)*(y"2 + 2) ságe: I2=integral(II,y,-2,2);12 8*pi ságe: V Matlabu se obejdeme bez předpokladů, nutno ovšem dodat, že výpočet posledního integrálů zde trval daleko delší dobu, než tomu bylo v Sage: » f=sym('(x+y)"2'); » ylims=solve('x~2+y~2-4','y') ylims = (2 - x)~(l/2)*(x + 2)~(l/2) -(2 - x)~(l/2)*(x + 2)~(l/2) » f=sym('(x+y)"2'); » ylims=solve('x~2+y~2-4','x') ylims = (2 - y)~(l/2)*(y + 2)~(l/2) -(2 - y)~(l/2)*(y + 2)~(l/2) » Il=int(f,'x',ylims(2),ylims(1)) 11 = (4*(y~2 + 2)*(2 - y)~(l/2)*(y + 2)~(l/2))/3 » I2=int(Il,'y',-2,2) 12 = 8*pi » 88 Kapitola 7 Řešení diferenciálních rovnic Budeme se zabývat zejména hledáním řešení rovnice ve tvaru V = f(x,y) pro nějakou funkci /, s případnou počáteční podmínkou y(xo) = yo pak hovoříme o počáteční úloze. Zobecněním této úlohy jsou rovnice vyšších řádů, z nichž velkou důležitost mají okrajové úlohy - tedy rovnice druhého řádu se zadanými podmínkami v krajních bodech nějakého intervalu [a, b]: y" = f(x,y,y), y(á) = y1, y(b) = y2 7.1 Analytické řešení Matlab i Sage obsahují nástroje pro nalezení analytického řešení diferenciální rovnice. V Matlabu se používá funkce dsolve, které v nejjednodušší podobě stačí jeden vstupní parametr - textový řetězec s rovnicí, kde derivace je označena jako D: Základní typy rovnic zvládá bez problémů: » dsolve('Dy = y') ans = C4*exp(t) » dsolve('Dy = -y') ans = C2/exp(t) » dsolve('Dy = t*y') 89 ans = C6*exp(ť2/2) » dsolveCDy = y/ť) ans = C8*t » Je možné použít i rovnici s parametrem a pokud Matlab rovnici nedokáže vyřešit, slušně nám to oznámí: » dsolveCDy = a*y*sin(t) ') ans = C10/exp(a*cos(t)) » dsolveCDy = cos(t*y)') Warning: Explicit solution could not be found. > In dsolve at 161 ans = [ empty sym ] » Lze také samozřejmě označit jinak výstupní nezávislou proměnnou nebo funkci, kterou hledáme, či zadat počáteční podmínku: » dsolveCDy = y"2','x') ans = 0 -1/CC38 + x) » dsolveCDx = -s*x"2','s') ans = 0 -l/(- s"2/2 + C54) » dsolveCDy = -y"2' , 'y(l)=l' , 'x' ) ans = 1/x » Problémem nejsou ani okrajové úlohy pro rovnice druhého řádu: 90 » dsolve('D2y = -y+cos(x)','y(0)=0 , 'y(pi/2)=ľ , 'x') ans = cos(3*x)/8 - cos(x)/8 + sin(x)*(x/2 + sin(2*x)/4) - sin(x)*(pi/4 - 1) » Řešení diferenciálních rovnic v Sage je také snadné. Nejprve definujeme x jako nezávislou proměnnou, y jakožto funkci proměnné x a pro nalezení řešení rovnice použijeme funkcei desolve. Zde se implicitně předpokládá, že pokud pravá strana rovnice není uvedena, je nulová. Stejnou rovnici tedy mohu zadat dvěma způsoby: sage: x=var('x') sage: y=function('y',x) sage: desolve(diff(y,x) ==x*y, y) c*e"(l/2*x"2) sage: desolve(diff(y,x) -x*y, y) c*e"(l/2*x"2) sage: V některých případech je potřeba z výsledku funkci y vyjádřit a to i při zadání počáteční podmínky: sage: desolve(diff(y,x) == y"2, y) -l/y(x) == c + x sage: desolve(diff(y,x) == y"2, y,[l,-l]) -l/y(x) == x sage: Pro rovnice vyšších řádu se uvádí řád derivace jako třetí parametr funcke diff. Zkusíme-li vyřešit v Sage stejnou okrajovou úlohu jako v Matlabu, dostaneme o něco jednodušší vyjádření výsledku: sage: desolve(diff(y,x,2) == -y+cos(x), y,[0,0,pi/2,1]] -l/4*(pi - 4)*sin(x) + l/2*x*sin(x) sage: 91 Čtenář si jako cvičení může zkusit ověřit, zda oba výsledky jsou totožné. V případě, že Sage nějakou rovnici nedokáže vyřešit (nebo uděláme chybu v zadání), objeví se chybové hlášení, ze kterého je vidět, že Sage pro řešení diferenciálních rovnic používá software zvaný Maxima. Více informací o něm mohou zájemci zjistit na internetových stránkách http: //maxima, sourcef orge .net/. 7.2 Numerické řešení Pro nalezení numerického řešení budeme používat Matlab. Numerické řešení hledáme na nějaké síti bodů (xj)™=0, takže nelze hledat obecné řešení, ale pouze přibližné partikulární řešení nějaké počáteční úlohy. V současné době patrně nejpoužívanější metody jsou metody Runge-Kutta, které patří mezi jednokrokové metody, to znamená, že výpočet přibližných hodnot hledané funkce y v bodě Xk+i probíhá na základě hodnoty v předchozím bodě Xk- Kromě toho známe také vícekrokové metody (pro výpočet v Xk+i se použije více předchozích hodnot), ale ty jsou v základním matlabovském balíku využívány jen velmi málo (funkce ode 113). Nejpoužívanějšími metodami Runge-Kutta jsou metody čtvrtého řádu, v Matlabu je jedna z nich implementována ve funkci ode45. Nápověda k této funkci je poměrně obsáhlá a taky v ní můžeme zjisti další funkce použitelné pro numerické řešení diferenciálních rovnic. Zaměříme se jen na tuto funkci a postačí nám základní použití. V nejjednodušší verzi musíme zadat funkci / z pravé strany diferenciální rovnice, interval, na němž chceme přibližné řešení spočítat a také počáteční podmínku. Výsledkem jsou dva vektory - síť bodů, a příslušné funkční hodnoty. Výsledek si zobrazíme: » f=inline('y*sin(t)'); » l=[0,pi] ; » cond= [1]; » [t,y]=ode45(f,I,cond); » plot(t,y) » 92 1 O 0.5 1 1.5 2 2.5 3 3.5 Protože nás zajímá, s jakou chybou je řešení určeno, zkusíme je porovnat s hodnotami řešení přesného: » dsolve('Dy = y*sin(x)','y(0)=l','x') ans = exp(l)/exp(cos(x)) >> ff=inline('exp(l)./exp(cos(x))') ff = Inline function: ff(x) = exp(l)./exp(cos(x)) » yy=ff(t); >> max(abs(y-yy)) ans = 2.1830e-05 » Chyba je tedy poměrně velmi malá, na grafickém znázornění by se přesné řešení od přibližného nedalo rozlišit. Co se týče počátečních podmínek, berou se automaticky v krajním bodě zadaného intervalu, ale je možné jich zadat více najednou: 93 » cond=[0.5 1 1.5 2] ; » [t,y]=ode45(f,I,cond); » plot(t,y) » 15 10 0 0.5 1 1.5 2 2.5 3 3.5 94 Kapitola 8 Statistické metody V této kapitole se nebudeme věnovat teorii pravděpodobnosti a teoretické statistice, zkusíme se zaměřit na praktické statistické výpočty. Co se týče pravděpodobnosti a statistiky, je Matlab skvěle vybaven, zejména pokud máte k dispozici statistický toolbox. Zde jsou implementovány všechny běžné metody a algoritmy, ale také spousta metod a algoritmů speciálních. Protože popis celého statistického toolboxu je nad rámec tohoto textu, omezíme se toliko na základní použití. Pokud chceme vyzkoušet některé dostupné programy a nemáme k dispozici vhodná data, není problém si je vygenerovat, nabízí se generátor náhodných čísel, který umí vytvořit data s požadovaným rozdělením. Tak například pro data s normálním rozdělením lze použít funkci normrnd, s exonenciálním exprnd, s rovnoměrným unif rnd, s beta rozdělením betarnd atd. Podobně pro výpočet hustoty náhodné veličiny jsou k použití funkce končící „pdf " (z anglického probability density function), tedy normpdf, unif pdf, exppdf apod., pro distribuční funkce můžeme požít funkce normcdf, unif cdf, expcdf atd. (z anglického cumulative distribution function). Pro kvantily se používají funkce norinv, unif inv, expinv apod. jako hodnoty inverzní distribuční funkce. Co se týče parametrů uvedených funkcí, zaleží samozřejmě na tom, jakými parametry je určeno konkrétní rozdělení, takže nejlepší bude poradit se v každém konkrétním případě s nápovědou. Pokud bychom se třeba chtěli podívat, jak se mění hustota normální rozdělení v závislosti na směrodatné odchylce, stačí zadat: >> x=linspace(-3,3); 95 » sl=0.5;s2=l;s3=2; » yl=normpdf(x,0,sl) >> y2=normpdf(x,O,s2) >> y3=normpdf(x,O,s3) » plot(x, [yl;y2;y3]) Když už jsem se zmínili o směrodatné odchylce, pro její výpočet se použává příkaz std, pro rozptyl příkaz var. A nesmíme zapomenout na střední hodnotu a medián - mean a medián. Tyto čtyři funkce nejsou součástí statistického toolboxu ale jsou již v základní matlabovské výbavě. Jedná se samozřejmě o odhady příslušných charakteristik náhodných veličin, nikoliv o přesné hodnoty. Zkusme porovnat u vygenerovaných dat teoretické hodnoty s vypočítanými. Pokud bychom v podobném duchu udělali rozsáhlé testování, mohou nám výsledky leccos napovědět o kvalitě generátoru náhodných čísel, který Matlab používá. Vygenerujeme 100 normálně rozdělených hodnot se střední hodnotou 1 a rozptylem 0.25. Nesmíme ovšem zapomenout, že jako parametr zadáváme směrodatnou odchylku: 96 » x=normrnd(l,0.5,1,100); >> mean(x) ans = 1.0615 >> median(x) ans = 1.0477 >> var(x) ans = 0.3378 » std(x) ans = 0.5812 » Vidíme, že u rozptylu je relativně největší chyba, z jednoho vzorku ovšem nemůžeme vyvozovat žádné závěry. 8.1 Intevalové odhady Jednou ze základních statistických úloh je určení intervalu, v němž s danou pravděpodobností leží nějaký parametr náhodného rozdělení. Nejčastěji předpokládáme, že máme data s normálním rozdělením nebo s rozdělením s normálního odvozeným. Podíváme se, jak se mění interval spolehlivosti pro odhad střední hodnoty normálních dat, pokud roste jejich počet. Nejprve předpokládejme, že známe hodnotu rozptylu, pak se používají kvantily standardizovaného normálního rozdělení: » n=100; » mu=l; >> sigma=0.5; >> X=normrnd(mu,sigma,1,n); » M=mean(X) M = 0.9813 » alpha=0.05; 97 » D=M-norminv(l-alpha/2,0,1)*sigma/sqrt(n); » H=M+norminv(l-alpha/2,0,1)*sigma/sqrt(n); » I=[D,H] I = O.8833 1.0793 » n=500; >> X=normrnd(mu,sigma,1,n); » M=mean(X) M = 0.9878 » D=M-norminv(l-alpha/2,0,1)*sigma/sqrt(n); » H=M+norminv(l-alpha/2,0,l)*sigma/sqrt(n); » I=[D,H] I = 0.9439 1.0316 » n=1000; >> X=normrnd(mu,sigma,1,n); » M=mean(X) M = 1.0266 » D=M-norminv(l-alpha/2,0,1)*sigma/sqrt(n); » H=M+norminv(l-alpha/2,0,l)*sigma/sqrt(n); » I=[D,H] I = 0.9956 1.0576 » V případě, že rozptyl není známý, použijeme kvantily Studentova rozdělení: » n=100; » alpha=0.05; » mu=l; >> sigma=0.5; >> X=normrnd(mu,sigma,1,n); » M=mean(X) M = 0.9590 » S=std(X) 98 s = O.5047 » D=M-tinv(l-alpha/2,n-l)*S/sqrt(n); » H=M+tinv(l-alpha/2,n-l)*S/sqrt(n); » I=[D,H] I = 0.8589 1.0592 » n=500; >> X=normrnd(mu,sigma,1,n); » M=mean(X) M = 1.0160 » S=std(X) S = 0.4909 » D=M-tinv(l-alpha/2,n-l)*S/sqrt(n) ; » H=M+tinv(l-alpha/2,n-l)*S/sqrt(n); » I=[D,H] I = 0.9729 1.0592 » n=1000; >> X=normrnd(mu,sigma,1,n); » M=mean(X) M = 1.0172 » S=std(X) S = 0.4930 » D=M-tinv(l-alpha/2,n-l)*S/sqrt(n) ; » H=M+tinv(l-alpha/2,n-l)*S/sqrt(n); » I=[D,H] I = 0.9866 1.0478 » Je jasné, že inerval se zkracuje (délka klesá s o odmocninou n), ale rozdíl mezi oběma případy není patriw. 99 8.2 Testování hypotéz Pro testování hypotéz je v Matlabu možné využít již hotových programů. Základním z nich je ttest, který v nejjednodušší verzi testuje, zda rozptyl náhodného výběru je nulový. Hladina významnosti je implicitně 5 procent: » n=100; » mu=l; >> sigma=0.5; >> X=normrnd(mu,sigma,1,n); » ttest(X) ans = 1 Výsledek 1 znamená zamítnutí hypotézy. Pokud bychom chtěli vědět také pravděpodobnost (tzv. jo-hodnota), přidáme výstupní parametr: » [H,P]=ttest(X) H = 1 P = 1.1433e-32 » Vidíme, že střední hodnota není nulová téměř jistě. Provedeme tady drobnou úpravu: » M=mean(X) M = 0.9926 » [H,P]=ttest(X,M) H = 0 P = 1 » Ještě si otestujeme, jak dopadne skutečná střední hodnota použitá při generování náhodného výběru: 100 » [H,P]=ttest(X,l) H = O P = 0.8946 » Funkce ttest se dá také použít na testování rovnosti středních hodnot dvou náhodných výběru stejného rozsahu: >> Y=normrnd(mu,sigma,1,n); » [H,P]=ttest(X,Y) H = 0 P = 0.5040 >> Y=normrnd(mu-0.1,sigma,1,n); » [H,P]=ttest(X,Y) H = 0 P = 0.2128 » K testování rovnosti středních hodnot u dvou výběrů se dá ooužít také funkce ttest2, zdá se ale, že pracuje poněkud jinak než funkce ttest: » [H,P]=ttest(X,Y) H = 0 P = 0.5040 » [H,P]=ttest2(X,Y) H = 0 P = 0.5276 » 101 8.3 Lineární regrese Studenti si lineární regresi občas pletou s konstrukcí regresní přímky, což je ale jen speciální případ lineární regrese. Na tuto problematiku jsme už narazili v prvním semestru, když jsem hovořili o pseudoinverzní matici. V lineární regresi předpokládáme, že pozorovaná data, která zpracováváme, jsou lineárními kombinacemi funkcí předem daných hodnot, přičemž v pozorování je nějaký náhodná chyba. Pozorování je tedy náhodná veličina Y(x) = + • • • + /3fc/fc(x) + e. Měříme hodnoty Y v zadaných bodech > k, tím dostáváme celkem soustavu n rovnic pro k neznámých. Matice soustavy se ve statistice nazývá matice plánu a zpravidla se označuje X. Při řešení se snažíme minimalizovat reziduum (rozdíl mezi levou a pravou stranou), k čemuž se výborně hodí pseudoinverzní matice. Jestliže navíc předpokládáme normalitu chyby veps se známým rozptylem, dají se spočítat i intervalové odhady pro parametry f3y Příklad 6. Mějme třídící algoritmus a chceme otestovat jeho kvalitu. Dobré třídící algoritmy pro n dat potřebují Oin log n) času, horší algoritmy potřebují O in2) času. Algoritmus jsme otestovali pro náhodné výběry různých rozsahů a následná tabulka udává průměrné časy pro jednotlivé rozsahy v milisekundách: n 100 200 300 400 500 t 0.42.883 127.39 264.08 435.69 665.0 n 600 700 800 900 1000 t 920.24 1220.4 1565.8 1941.4 2348.6 Vyrobíme matici plánu, přičemž předpokládáme, že výsledná funkce může být polynom druhého stupně doplněná funkcemi log nan log n. » n n = Columns 1 through 5 100 200 300 400 500 Columns 6 through 10 600 700 800 900 1000 » t t = Columns 1 through 5 102 42.883 127.39 264.08 435.69 665.08 Columns 6 through 10 920.24 1220.4 1565.8 1941.4 2348.6 » Xl=ones(10,1); » X2=n'; » X3=(n."2)'; » X4=(log(n))'; » X5=(n.*log(n))'; » X=[X1,X2,X3,X4,X5] X = 1 100 10000 4 6052 460.52 1 200 40000 5 2983 1059.7 1 300 90000 5 7038 1711.1 1 400 1 6e+05 5 9915 2396.6 1 500 2 5e+05 6 2146 3107.3 1 600 3 6e+05 6 3969 3838.2 1 700 4 9e+05 6 5511 4585.8 1 800 6 4e+05 6 6846 5347.7 1 900 8 le+05 6 8024 6122.2 1 1000 le+06 6 9078 6907.8 » Nakonec provedeme výpočet odhadu parametrů s použitím pseudoinverzní matice: » Y=t'; » beta=pinv(X)*Y beta = -278.98 -5.6652 0.0013371 99.146 0.90803 » Zdá se že koeficient u druhé mocniny je zanedbatelný oproti ostatním koeficientům. Podíváme se ještě, jak funkce, kterou jsme získali, prochází aproximuje data. Přitom ověříme rozdíl mezi tím když uvedený koeficient zanedbáme a když nikoliv: 103 » 1111=100:1000; » fl=beta(l)+beta(2)*nn+beta(4)*log(nn)+... beta(5)*nn.*log(nn); » f2=beta(l)+beta(2)*nn+beta(3)*nn. "2+beta(4)*log(nn)--. . . beta(5)*nn.*log(nn); » plot(n,t, '*r' ,nn,f 1, 'b' ,nn,f 1, 'g—') » 2500 Vidíme, že koeficient zanedbat nemůžeme, protože je sice řádově daleko menší než ostatní koeficienty, ale vzhledem k tomu, že se jím násobí velká čísla (druhé mocniny n), nelze jej vynechat. 104 Literatura [1] Horová, I.; Zelinka, J.: Numerické metody. Brno: Masarykova univerzita, druhé vydání, 2008, ISBN 978-80-210-3317-7. [2] Kobza, J.: Splajny. Univerzita Palackého, 1993, ISBN 9788070672655. URL http://books.google.cz/books?id=-0nXYgEACAAJ [3] Ralston, A.: Základy numerické matematiky. Praha: Academia, druhé vydání, 1978. [4] Schumaker, L.: Spline Functions: Basic Theory. Cambridge University Press, ISBN 9781139463430. URL http://books.google.cz/books?id=2uZLawUhXfgC 105