Numerické výpočty – I Jiří Zelinka Tento učební text vznikl za přispění Evropského sociálního fondu a státního rozpočtu ČR prostřednictvím Operačního programu Vzdělávání pro konkurenceschopnost v rámci projektu Univerzitní výuka matematiky v měnícím se světě (CZ.1.07/2.2.00/15.0203). Obsah Úvod 5 Kapitola 1. První pokusy 7 1. Praktické výpočty 7 2. Kombinatorické výpočty 14 3. Diferenční rovnice 18 4. Pravděpodobnost 21 5. Geometrie v rovině 24 Kapitola 2. Základy lineární algebry 31 1. Matice a vektory v Matlabu 31 2. Matice a vektory v Sage 39 3. Příklady 43 Kapitola 3. Lineární modely a maticový počet 53 1. Populační modely 53 2. Markovovy řetězce 57 3. Maticové rozklady a spol 58 Kapitola 4. Geometrie 67 3 Úvod Pane profesore, co jsou to vlastně čísla? Abych pravdu řekl, Kropáčku, tak to úplně přesně nevím. Opravdu nevíte? Ne, vlastně si ani nejsem jistý, jestli vůbec existují. Ale jinak toho o nich vím poměrně hodně. Pokud se zabýváme čísly, zpravidla se nevyhneme tomu, že někdy potřebujeme i něco spočítat. Vlastně jsme se učili nejdřív něco spočítat a až následně se občas najde nějaký šťoura, který se ptá, co to vlastně čísla jsou. V tomto textu se budeme věnovat spíš praktickému počítání, přičemž teoretické předpoklady budou někde v pozadí, ze kterého je v případě potřeby vytáhneme. Počítání nám velmi usnadňuje technika. Zpravidla si bereme na pomoc kalkulačku i při sčítání poměrně malých čísel, i když s trochou snahy bychom snad ještě uměli ručně i dělit. Nejefektivnějším pomocníkem může být samozřejmě počítač, se kterým kromě rippování CD, kódování filmů a úpravy digitálních fotek je možné pomocí vhodných prostředků i něco spočítat. Pro potřeby tohoto kursu budeme využívat hlavně dva softwarové prostředky - jsou to matematické programy Matlab a Sage. První z nich komerční a je orientován hlavně numericky, druhý je volně dostupný a zvládá i symbolické výpočty. Takže jejich vhodnou kombinací dokážeme pokrýt prakticky všechny výpočty, které budeme potřebovat. Je ovšem třeba brát v úvahu i vývoj tohoto software, proto některé některé zde uváděné informace mohou být časem zastaralé nebo neplatné. Nechť mi to laskavý čtenář promine. O dalších výpočetních prostředcích se zmíníme jen okrajově. 5 KAPITOLA 1 První pokusy 1. Praktické výpočty Kropáčku, co se s tím tak páráte, vy to ještě nemáte? Mám, pane profesore, už jsem to vypočítal. A kolik vám to vyšlo? Čtyřicet dva1. Čtyřicet dva? To je divné. Co jste to vlastně počítal? To se právě snažím zjistit. Učitelé matematiky na všech úrovních se s tím setkávají poměrně často. Zadají příklad a student nebo žák okamžitě začne počítat za použití nějakých vzorečků. Netvrdím, že výpočet musí být vždy špatně, ale zamyslet se nad problémem ještě předtím, než se na něj vrhnu se vzorečky a s kalkulačkou, je zpravidla prospěšné. Často nám to může pomoci odhalit chyby, kterých se při výpočtu dopustíme. Například když nám vyjde objem nějakého tělesa jak komplexní číslo nebo teplota materiálu −358.16◦ C. Při používání techniky je ovšem třeba mít na paměti ještě další věc, a tou jsou její omezené možnosti. Při převážné většině výpočtů se k hranicím, které výpočetní technika má, nepřiblížíme, ale je dobré o nich vědět. Tyto hranice jsou dány zejména omezeným prostorem, který je pro reprezentaci čísel vymezen. Nejlépe si to ukážeme na následujícím příkladu. Příklad 1. Čísla jsou v počítači uložena ve dvojkové soustavě, takže i jednoduchá desetinná čísla, např. 0,1, jsou při standardním způsobu ukládání uložena nepřesně, jelikož ve dvojkové soustavě je číslo 0,1 zapsáno pomocí nekonečné řady: 0, 1 = 2−4 + 2−5 + 2−8 + 2−9 + 2−12 + 2−13 + 2−16 + 2−17 + · · · a tuto řadu musíme někde useknout. V Matlabu se dá snadno zjistit, velikost chyby, která se v tomto případě do výpočtu může vloudit: 1 Viz též Douglas Adams: Stopařův průvodce po Galaxii 7 8 1. PRVNÍ POKUSY >> c=0; >> for k=1:100, c=c+0.1; end >> 10-c ans = 1.9540e-14 >> Nebo elegantněji >> c=0.1*ones(1,100);10-sum(c) ans = 1.9540e-14 >> Chyba je sice nepatrná, ale je tam. Často je možné tuto chybu zmenšit správným postupem při výpočtu. V uvedeném příkladu stačí celkový součet rozdělit na několik částečných součtů. To můžeme udělat tak, že stoprvkový vektor c přeskládáme do čtvercové matice řádu 10 (pokud ještě nevíte, co je to matice, tak si ji můžete představit jako tabulku, v našem případě o deseti řádcích a deseti sloupečcích), pak sečteme nejprve jednotlivé sloupce a následně získané mezivýsledky: >> c=0.1*ones(1,100);A=reshape(c,10,10); >> 10-sum(sum(A)) ans = 1.7764e-15 >> Vidíme, že chyba je zhruba desetkrát menší než při přímém součtu. Tuto skutečnost můžeme zdůvodnit následující úvahou. Předpokládejme, že číslo 0,1 je v počítači uloženo s chybou ε. Pokud sčítáme sto stejných hodnot, na konci výpočtu připočítáváme číslo 0,1 k číslu přibližně stokrát většímu, Protože zaokrouhlení výsledku se bere podle většího čísla, hodnota 0,1 se přičte s chybou zhruba 100ε. Je to podobné, jako bychom počítali s přesností na tři platné číslice a K číslu 111 přičítali 1,11. Celková chyba bude řádově v tisících ε, záleží na konkrétním čísle a „chytrosti software. Pokud ovšem sčítáme jen deset hodnot, je chyba řádově stokrát menší, při dalším sčítání mezivýsledků se sice zvětší, ale nedosáhne takové hodnoty jako při přímém sčítání. Při použití programu Sage si často můžeme snadno zvolit, jestli výpočet proběhne symbolicky nebo numericky. V případě symbolického výpočtu zadáváme konstanty ve formě zlomků, pro numerický výpočet použijeme desetinný zápis: 1. PRAKTICKÉ VÝPOČTY 9 sage: M=matrix(1,[1/10 for i in range(100)]) sage: 10-sum(sum(M)) 0 sage: M=matrix(1,[0.1 for i in range(100)]) sage: 10-sum(sum(M)) 1.95399252334028e-14 sage: M=matrix(10,[0.1 for i in range(100)]) sage: 10-sum(sum(M)) 1.77635683940025e-15 Také v Matlabu je možné zobrazovat hodnoty ve formě racionálních čísel, umožní to příkaz format rat. V tomto případě se ale jedná pouze o aproximace čísel uložených v počítači standardním způsobem. Takže pokud třeba potřebujeme aproximovat číslo π zlomkem, můžeme zadat >> format rat >> pi ans = 355/113 >> Matlab i Sage umí samozřejmě pracovat s komplexními čísly. V Matlabu slouží jako imaginární jednotka symbol i, na zjištění goniometrického tvaru komplexního čísla můžeme použít funkce abs a angle: >> z=1+i >> abs(z) ans = 1.4142 >> angle(z) ans = 0.7854 >> Poslední uvedená hodnota je numericky π/4. Pro zadání komplexních čísel v Sage použijeme imaginární jednotku označenou I, na zjištění goniometrického tvaru funkce argument a abs. Přitom funkci argument je potřeba použít na numerickou aproximaci komplexního čísla, a použití je následovné: sage: abs(1+I) sqrt(2) sage: numerical_approx(1+I).argument() 10 1. PRVNÍ POKUSY 0.785398163397448 sage: (1+I).numerical_approx().argument() 0.785398163397448 Patrně si pozorný čtenář už povšiml, že funkce numerical approx je možné použít dvěma způsoby: buď obvykle, kdy jak argument použijeme výraz, jehož aproximaci chceme vypočítat, nebo až za výrazem s tečkou jako oddělovačem. Tento dvojí způsob aplikace funkce je v Sage celkem obvyklý. Funkce numerical approx má navíc ještě další nepovinné parametry. Pro nás jsou zajímavé hlavně ty, které mohou specifikovat, s jak velkou přesností se má výraz zobrazit. Přirozené číslo jako další parametr určuje, kolik bitů se má pro aproximaci použít, nebo je možné specifikovat pomocí slova digits, kolik číslic se má vypsat. A pokud se někomu nechce vypisovat celý název funkce, může použit zkrácené názvy, které jsou n nebo N: sage: numerical_approx(pi) 3.14159265358979 sage: n(pi) 3.14159265358979 sage: pi.N() 3.14159265358979 sage: pi.N(100) 3.1415926535897932384626433833 sage: N(pi,100) 3.1415926535897932384626433833 sage: n(pi,digits=20) 3.1415926535897932385 sage: pi.n(digits=20) 3.1415926535897932385 Příklad 2. Podíváme se, jak si softwarové prostředky poradí s řešením rovnice x3 + x2 − 2x − 1 = 0 z příkladu 1.2. Jedná se o hledání kořenů polynomu třetího stupně, neboli kubického polynomu. V Matlabu se polynomy zadávají pomocí koeficientů od nejvyšší mocniny, všechny koeficienty tak tvoří vektor, který je v Matlabu definován pomocí hranatých závorek. Kořeny polynomu pak zjistíme pomocí funkce roots: >> p=[1 1 -2 -1]; >> roots(p) ans = -1.8019 1.2470 1. PRAKTICKÉ VÝPOČTY 11 -0.4450 >> V Sage slouží pro řešení rovnic funkce solve, ale příkaz sage: x=var(’x’) sage: solve(x^3+x^2-2*x-1,x) dá výsledek x = − 1 2 7 18 i √ 3 + 7 54 (1 3 ) i √ 3 + 1 − −7i √ 3 + 7 18 7 18 i √ 3 + 7 54 (1 3 ) − 1 3 , x = − 1 2 −i √ 3 + 1 7 18 i √ 3 + 7 54 (1 3 ) − 7i √ 3 + 7 18 7 18 i √ 3 + 7 54 (1 3 ) − 1 3 , x = 7 18 i √ 3 + 7 54 (1 3 ) + 7 9 7 18 i √ 3 + 7 54 (1 3 ) − 1 3 Je zřejmé, že z tohoto výsledku se nedá při nejlepší vůli usoudit, že imaginární část je ve skutečnosti nulová. K získání numerické hodnoty řešení můžeme použít funkci numerical approx(): sage: x=var(’x’) sage: S=solve(x^3+x^2-2*x-1.0,x) sage: x0=S[0].right() sage: x0.numerical_approx() -0.445041867912629 takže všechna řešení dostaneme např. takto: sage: x=var(’x’) sage: S=solve(x^3+x^2-2*x-1.0,x) sage: x0=S[0].right().numerical_approx() sage: x1=S[1].right().numerical_approx() sage: x2=S[2].right().numerical_approx() sage: x0;x1;x2 -0.445041867912629 -1.80193773580484 1.24697960371747 + 5.55111512312578e-17*i 12 1. PRVNÍ POKUSY Zdálo by se, že na hledání kořenů polynomů bude Matlab lepším prostředkem než Sage, ale stačí zkusit najít kořeny polynomu p(x) = (x − 1)4 = x4 − 4x3 + 6x2 − 4x + 1, nastavit věší přesnost zobrazování výsledků a budeme rychle vyvedeni z omylu: >> format long >> p=[1 -4 6 -4 1]; >> roots(p) ans = 1.000217162531685 0.999999969335974 + 0.000217131858872i 0.999999969335974 - 0.000217131858872i 0.999782898796371 >> Jak se dá čekat, Sage si s tímto problémem poradí podstatně lépe. sage: solve(x^4-4*x^3+6*x^2-4*x+1==0,x) [x == 1] > V posledním příkladu v této kapitole si ukážeme typický problém při numerických výpočtech, kterým je ztráta přesnosti při počítání s různě velkými čísly. Příklad 3. Mějme polynom P(x) = (x − 3)3 = x3 − 9x2 + 27x − 27. Předpokládejme, že potřebujeme hodnoty polynomu v malém okolí bodu 3, délka okolí bude 0,00001. Výsledek si zobrazíme graficky >> P=[1 -9 27 -27]; >> x=linspace(2.99995,3.00005,201); >> y=polyval(P,x); >> plot(x,y) 1. PRAKTICKÉ VÝPOČTY 13 2.9999 3 3 3 3 3 3.0001 −2 −1.5 −1 −0.5 0 0.5 1 1.5 x 10 −13 Výsledkem na daném intervalu by měla být čísla řádově 10−15 , protože při výpočtu používáme čísla řádově v desítkách, je chyba zhruba desetkrát větší. V tomto případě i Sage počítá numericky, protože příkaz sage: g=plot(x^3-9*x^2+27*x-27,2.99995,3.00005);g vyprodukuje následující obrázek: 3e 3e 3e 3e 3e -1e-13 -5e-14 0 5e-14 1e-13 Chybu ovšem můžeme výrazně snížit, pokud použijeme první vztah pro výpočet polynomu: >> x=linspace(2.99995,3.00005,201); >> y=(x-3).^3; >> plot(x,y) 14 1. PRVNÍ POKUSY 2.9999 3 3 3 3 3 3.0001 −2 −1.5 −1 −0.5 0 0.5 1 1.5 x 10 −13 Závěrem z této kapitoly nechť je vědomí, že ani velmi kvalitní výpočetní prostředek (zatím?) nenahradí člověka a vždycky je dobré se alespoň trochu zamyslet nad tím, jaké nám to vlastně počítač dává výsledky, co s nimi chceme dále provádět a jestli by to nešlo spočítat trochu přesněji. 2. Kombinatorické výpočty Základním výpočtem v kombinatorice je výpočet faktoriálu, který udává počet permutací a vyskytuje se v dalších kombinatorických vztazích. V Matlabu se pro tento účel hodí funkce prod, která vrací jako výsledek součin prvků vektoru, který je jejím argumentem. Takže 10! zjistíme příkazem >> prod(1:10) ans = 3628800 >> Tento postup funguje dobře i pro 0!: >> prod(1:0) ans = 1 >> Tímto standardním způsobem lze v Matlabu na 32 bitovém procesoru získat maximálně faktoriál 170: >> prod(1:170) ans = 2. KOMBINATORICKÉ VÝPOČTY 15 7.2574e+306 >> prod(1:171) ans = Inf >> Pokud bychom potřebovali pracovat s faktoriály větších čísel, museli bychom použít speciální knihovny. V tomto ohledu je Sage daleko před Matlabem. Můžeme zde použít pro výpočet funkci factorial: sage: factorial(10) 3628800 sage: factorial(100) 933262154439441526816992388562667004 907159682643816214685929638952175999 932299156089414639761565182862536979 208272237582511852109168640000000000 00000000000000 Sage si bez nesnází poradí s faktoriálem 10 000, ale pokud byste chtěli vědět, jak velký je výsledek, použijte místo počítání číslic raději numerickou aproximaci: sage: factorial(10000).numerical_approx() 2.84625968091705e35659 Počet číslic výsledku by nám prozradil i dekadický logaritmus: sage: log(factorial(10000),10).numerical_approx() 35659.4542745208 Když hovoříme o faktoriálech, nesmíme opominout gamma funkci Γ. Ta je definována pomocí integrálu a pro naše účely stačí vědět, že pro přirozená čísla platí Γ(n + 1) = n!. V Matlabu je označena jako gamma, v Sage taktéž gamma. Zajímavé je, že v Sage dává hodnoty této funkce i funkce factorial, nesmíme ovšem zapomenout na posun argumentu o 1. Můžete si vyzkoušet příkazy sage: gamma(1.5) sage: factorial(1.5) sage: factorial(0.5) Pokud se týká výpočtu kombinačních čísel, podívejme se nejdříve, jak bychom postupovali při ručním počítání, třeba v Matlabu. Výpočet 16 1. PRVNÍ POKUSY podle vztahu n k = n! k!(n − k)! je ovšem neefektivní. Lepší postup je určit m jako menší z čísel k a n − k a pak n k = n(n − 1) · · · (n − m + 1) m! . Matlabovská funkce pro výpočet kombinačního čísla pak může vypadat třeba takto: function kc=komb_c(n,k) m=min([k,n-k]); kc=prod(n-m+1:n)/prod(1:m); Pomocí ní můžeme spočítat hodnotu např. 300 100 , ale pro 300 150 už dostaneme výsledek mimo číselný rozsah Matlabu: >> komb_c(300,150) ans = Inf >> V tomto případě nám může pomoci beta funkce označovaná B, která podobně jako gamma funkce je definována pomocí integrálu, a platí pro ni vztah B(x, y) = Γ(x)Γ(y) Γ(x + y) . Snadno nahlédneme. že n k = n! k!(n − k)! = Γ(n + 1) Γ(k + 1)Γ(n − k + 1) = = Γ(n + 1)Γ(n + 2) Γ(k + 1)Γ(n − k + 1)Γ(n + 2) = Γ(n + 1) Γ(n + 2) 1 Γ(k+1)Γ(n−k+1) Γ(n+2) = = 1 (n + 1)B(k + 1, n − k + 1) S použitím tohoto vztahu dokáže Matlab spočítat daleko větší kombinační čísla: >> n=1000;k=500; >> 1/((n+1)*beta(k+1,n-k+1)) ans = 2.7029e+299 >> 2. KOMBINATORICKÉ VÝPOČTY 17 Matlabovská standardní funkce pro výpočet kombinačních čísel je nchoosek. Kdy se podíváme na její zdrojak, zjistíme, že je výpočet je prováděn trochu chytřeji, než ve výše uvedené námi vytvořené funkci: if k > n/2, k = n-k; end nums = (n-k+1):n; dens = 1:k; nums = nums./dens; c = round(prod(nums)); Při výpočtu samozřejmě vznikají numerické chyby při jednotlivých děleních, jsou ale dostatečně malé, aby se zaokrouhlením odstranily. Pro výpočet kombinačních čísel se v Sage používá funkce binomial. Zvládá spočítat hodně velká kombinační čísla, např. sage: binomial(10000,100); 65208469245472575695415972927215718683781335425416 74337221024717286920652077017898892751029134055299 08478530306159470981182823719823927054792711952961 27415562705948429404753632271959046657595132854990 606768967505457396473467998111950929802400 a když si asi minutu počkáte, můžete zjistit i kolik je třeba 1 000 000 500 000 . Kromě toho, abychom věděli, kolik kombinací kombinací k-tého stupně z n prvků existuje, se nám taky může hodit všechny tyto kombinace určit. V Matlabu to umí uvedená funkce nchoosek, když jako první parametr zadáme vektor čísel, z nichž se mají kombinace vybírat. Podobně funguje i funkce combnk, z výstupu je ovšem vidět, že každá pracuje jinak: >> nchoosek(1:5,3) ans = 1 2 3 1 2 4 1 2 5 1 3 4 1 3 5 1 4 5 2 3 4 2 3 5 2 4 5 3 4 5 >> combnk(1:5,3) ans = 18 1. PRVNÍ POKUSY 3 4 5 2 4 5 2 3 5 2 3 4 1 4 5 1 3 5 1 3 4 1 2 5 1 2 4 1 2 3 >> Sage pro tento účel používá funkci combinations: sage: mset = [1,2,3,4,5,6] sage: combinations(mset,2) [[1, 2], [1, 3], [1, 4], [1, 5], [1, 6], [2, 3], [2, 4], [2, 5], [2, 6], [3, 4], [3, 5], [3, 6], [4, 5], [4, 6], [5, 6]] Je také možné si nechat vypsat všechny permutace, v Matlabu použijeme funkci perms, v Sage funkci permutations. Co se týče dalších kombinatorických pojmů jako např. variace či kombinace s opakováním, pevně věřím, že určení jejich počtu pomocí některého v tomto textu používaného software by již čtenáři nečinilo žádný problém. 3. Diferenční rovnice Ukažme si na několika příkladech, jak je možné pracovat s diferenčními rovnicemi, aniž bychom se snažili určit explicitní řešení. Příklad 4. Vezměme si splácení auta z příkladu 1.20. Výchozí částka je 300 000 korun, úroková míra m = 6% = 0, 06, úroky se počítají měsíčně. Je-li měsíční splátka S a zůstatek po k měsících označíme dk, platí d0 = 300000, dk = dk−1 + m 12 dk−1 − S. Pokud chceme v Matlabu spočítat a zobrazit, jak rychle bude dluh klesat během tří let splácení při platbě 5 000 měsíčně, můžeme postupovat třeba takto: >> mes=36; >> d_k=300000;m=0.06;S=5000;D=d_k; >> for k=1:mes, d_k=d_k*(1+m/12)-S;... D=[D,d_k];end >> plot(0:mes,D,’*’); 3. DIFERENČNÍ ROVNICE 19 0 5 10 15 20 25 30 35 40 1.6 1.8 2 2.2 2.4 2.6 2.8 3 x 10 5 Může nás ale taky zajímat, jak poroste rozdíl mezi tím, kolik jsme zaplatili a tím, kolik nám ubylo z dluhu: >> mes=36; >> d_k=300000;m=0.06;S=5000;R=0; >> for k=1:mes, d_k=d_k*(1+m/12)-S;... R=[R,k*S-(300000-d_k)];end >> plot(0:mes,R,’*’); 0 5 10 15 20 25 30 35 40 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x 10 4 Případně si můžeme zobrazit vývoj dluhu pro různě velké splátky: >> mes=36; >> d_k=300000;m=0.06;S=5000;D1=d_k; >> for k=1:mes, d_k=d_k*(1+m/12)-S; D1=[D1,d_k];end >> d_k=300000;m=0.06;S=7000;D2=d_k; 20 1. PRVNÍ POKUSY >> for k=1:mes, d_k=d_k*(1+m/12)-S; D2=[D2,d_k];end >> d_k=300000;m=0.06;S=9000;D3=d_k; >> for k=1:mes, d_k=d_k*(1+m/12)-S; D3=[D3,d_k];end >> plot(0:mes,D1,’*’,0:mes,D2,’o’,0:mes,D3,’+’); 0 5 10 15 20 25 30 35 40 0 0.5 1 1.5 2 2.5 3 x 10 5 Další možností by mohlo být zobrazit vývoj dluhu při různých úrokových mírách, ale to už by jistě čtenář zvládl sám. Příklad 5. Uvažujme super banku, která nabízí při roční výpovědní lhůtě úrok 100 %. Takže při ročním úročení dostaneme za každou vloženou korunu koruny dvě. Co kdyby se ale úročilo pololetně: >> n=2;d_k=1;m=1; >> for k=1:n,d_k=(1+m/n)*d_k;end >> d_k d_k = 2.2500 A co třeba čtvrtletně: >> n=4;d_k=1;m=1; >> for k=1:n,d_k=(1+m/n)*d_k;end >> d_k d_k = 2.4414 Při měsíčním úročení dostaneme: >> n=12;d_k=1;m=1; >> for k=1:n,d_k=(1+m/n)*d_k;end >> d_k 4. PRAVDĚPODOBNOST 21 d_k = 2.6130 V případě, že by banka úročila týdně, či dokonce denně, vyšlo by d k = 2.6926, resp. d k = 2.7146. K čemu ty hodnoty asi konvergují? 4. Pravděpodobnost Tak mi řekněte, jaká je pravděpodobnost jevu A. Je to jedna polovina, pane profesore. Jak jste na to, proboha, přišel, Kropáčku? No, buď to vyjde, nebo ne. Pro nejrůznější pravděpodobnostní příklady se nám může hodit generování náhodných výsledků, které můžeme použít při simulacích náhodných jevů. Kromě funkcí zmíněných v kombinatorické části budeme používat i generátor náhodných čísel, i když ve skutečnosti se zpravidla jedná o čísla pseudonáhodná. V Matlabu je základním generátorem funkce rand, která dává čísla mezi 0 a 1, každé z těchto čísel by mělo mít stejnou pravděpodobnost, že bude zvoleno, v Sage je to random. Existuji ovšem funkce pro generování náhodných čísel podle zvolených kritérií, nejčastěji to bývají celá čísla v určeném rozsahu. V obou balících je pro tento účel funkce randint, takže stejného výsledku se můžeme dobrat více způsoby. V Matlabu se ovšem v nápovědě dočteme, že tato funkce již je zastaralá a novějších verzích Matlabu bude patrně scházet. Proto budeme používat funkci randi, která má jen trochu odlišnou syntaxi. V Sage navíc existuje funkce random matrix, která je použitelná pro generování náhodných čísel z daného číselného okruhu, takže s jeho pomocí můžeme generovat matice reálné, racionální, celočíselné i komplexní. Příklad 6. Vytvořme matici o dvou řádcích a čtyřech sloupcích s náhodnými celočíselnými prvky mezi 1 a 8. Nejprve možná řešení v Matlabu: >> % pomocí zaokrouhlení >> A=round(7*rand(2,4))+1 A = 1 3 5 4 6 4 8 8 >> % pomocí zaokrouhlení dolů >> A=floor(8*rand(2,4))+1 A = 2 3 4 8 22 1. PRVNÍ POKUSY 1 2 3 8 >> % pomocí zaokrouhlení nahoru >> A=ceil(8*rand(2,4)) A = 5 6 3 1 4 6 4 8 >> % přímo >> A=randi([1,8],2,4) A = 8 1 8 2 2 5 6 3 >> % o něco jednodušeji >> A=randi(8,2,4) A = 3 2 3 3 7 8 2 5 A řešení v Sage: sage: A=matrix(2,[randint(1,8) for i in range(8)]);A [6 5 7 6] [5 8 1 6] sage: A=random_matrix(ZZ,2,4,x=1,y=9);A [2 7 1 8] [4 5 7 8] V posledně uvedeném příkazu se horní mez nedosahuje, proto je potřeba jako její hodnotu vzít 9. Generování náhodných posloupností nebo matic je velmi užitečné kromě simulování některých náhodných jevů také při metodách souhrnně označovaných jako Monte Carlo. Princip spočívá v tom, že si vygenerujete množství simulovaných dat a z nich vybereme ty, které odpovídají nějakému jevu. Tímto způsobem se dají přibližně počítat nejen pravděpodobnosti jevů, ale také třeba plošné obsahy apod. Demonstrujeme si to na dvou příkladech. Příklad 7. Určíme přibližnou hodnotu čísla π (viz též učebnice). K tomu využijeme následující skutečnost: pravděpodobnost toho, že dvě náhodně zvolená nezávislá reálná čísla z intervalu [0,1] mají součet druhých mocnin menší nebo roven jedné, je rovna ploše čtvrtiny jednotkového kruhu, což je π/4. V Matlabu tedy prvky náhodné matice o dvou řádcích umocníme po jednotlivých prvcích na druhou, zjistíme, kolik z 4. PRAVDĚPODOBNOST 23 nich je menších nebo rovno jedné, vydělíme celkovým množstvím dat, vynásobíme 4 a odečteme od π: >> % 100 hodnot >> n=100;A=rand(2,n);A2=A.^2; >> pi-4*sum(sum(A2)<=1)/n ans = -0.3784 >> % 1000 hodnot >> n=1000;A=rand(2,n);A2=A.^2; >> pi-4*sum(sum(A2)<=1)/n ans = 0.0256 >> % 1 000 000 hodnot >> n=1000000;A=rand(2,n);A2=A.^2; >> pi-4*sum(sum(A2)<=1)/n ans = 0.0016 Vidíme, že s přesností to není žádná sláva. Pro miliardu hodnot už Matlab odmítá matici vytvořit. Bylo by samozřejmě možné celý výpočet provádět v cyklu, čímž bychom mohli aproximace čísla π takto spočítat ještě přesněji. Případný zájemce by jistě zvládl si příslušný prográmek vyrobit sám. Příklad 8. V dalším příkladu trochu předběhneme integrační počet. Zkusíme vypočítat plochu pod funkcí sinus na intervalu [0, π]. Budeme postupovat podobně jako v předchozím příkladě, vygenerujeme náhodné dvojice čísel na obdélníku [0, π] × [0, 1] a určíme poměrnou část těch, které leží pod grafem. Ta by měla být přibližně rovna poměru plochy pod grafem funkce k ploše celého obdélníka. Výsledek můžeme porovnat s přesnou hodnotou, která, jak se dozvíme v blízké budoucnosti, je rovna 2. >> % 100 hodnot >> n=100;x=pi*rand(1,n);y=rand(1,n); >> 2-pi*sum(y<=sin(x))/n ans = -0.1049 >> % 1000 hodnot >> n=1000;x=pi*rand(1,n);y=rand(1,n); >> 2-pi*sum(y<=sin(x))/n ans = -0.0389 24 1. PRVNÍ POKUSY >> % 1 000 000 hodnot >> n=1000000;x=pi*rand(1,n);y=rand(1,n); >> 2-pi*sum(y<=sin(x))/n ans = -9.1492e-04 Na závěr této části bych jen rád podotkl, že pro Matlab existuje speciální knihovna (toolbox) pro statistiku, který obsahuje velké kvantum nejrůznějších pravděpodobnostních a statistických funkcí. Pro ty, kteří se budou potřebovat těmito záležitostmi zabývat do větší hloubky, uvádím odkaz na dokumentaci k tomuto toolboxu, která je skutečně velmi obsáhlá: http://www.mathworks.com/help/pdf doc/stats/stats.pdf 5. Geometrie v rovině V této kapitole si ukážeme základní způsoby, jak je možné graficky znázorňovat některé rovinné geometrické objekty. Základní objektem je úsečka. Pro její zobrazení se v Matlabu i v Sage používá funkce line, syntaxe v je jednotlivých softwarech je, jak se dá očekávat, odlišná. V Matlabu při nejjednodušším použití funkce line zadáváme dva parametry, první z nich je vektor x-ových souřadnic bodů, které se mají spojit úsečkami, druhým parametrem je vektor y-ových souřadnic těchto bodů. Jednoduchou úsečku tedy získáme příkazem >> line([1,2],[2,4]) Pokud použijeme větší než dvouprvkové vektory, získáme lomenou čáru. Pro lepší přehled si také můžeme zobrazit mřížku. >> line([1,2,3,5,4],[2,4,0,1,5]) >> grid on 5. GEOMETRIE V ROVINĚ 25 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 V Sage je jako parametr funkce line seznam bodů, které se mají spojit úsečkami. Takže podobný obrázek, jako je předchozí, získáme v Sage příkazem sage: line([(1,2), (2,4), (3,0), (5,1), (4,5)]) 1 2 3 4 5 0 1 2 3 4 5 Pokud se nám zdá, že druhý obrázek je oproti prvním poněkud zdeformován, je to způsobeno tím, že měřítko není v obou osách stejné. Obdélníčky mřížky u prvního obrázku by ve skutečnosti měly být čtverečky. Obrázky Matlabu se přizpůsobují oknu ve kterém se zobrazují, v Sage je možné poměr os nastavit pomocí aspect ratio, jak si ukážeme za chvíli. V Matlabu pro zobrazování v rovině existuje univerzální funkce plot, která v základním tvaru má dva parametry podobně jako funkce line. Kromě nich lze použít další nepovinné parametry, které ovlivňují barvu obrázku, styl čáry a podobně. Například čárkovanou sinusovku získáme použitím příkazů 26 1. PRVNÍ POKUSY >> x=linspace(0,2*pi);y=sin(x); >> plot(x,y,’--’) 0 1 2 3 4 5 6 7 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Nevýhodou tohoto postupu je, že i jednoduché geometrické objekty si musíme nejdříve spočítat, přesněji řečeno, musíme spočítat body, které je tvoří a to na dostatečně husté síti. Třeba elipsu můžeme nakreslit takto: >> t=linspace(0,2*pi);x=3*cos(t);y=2*sin(t); >> plot(x,y) −3 −2 −1 0 1 2 3 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 V Sage existuje pro zobrazování objektů v rovně funkcí několik. Např. pro kružnice je tu zvláštní funkce circle, která má dva povinné parametry – střed a poloměr, nicméně obrázek získaný příkazem 5. GEOMETRIE V ROVINĚ 27 sage: circle((1,2),2) by na první pohled připomínal elipsu. Proto uděláme malinko složitější konstrukci, současně nastavíme barvu kružnice na zelenou. sage: c=circle((1,2),2, rgbcolor=(0,1,0)) sage: c.show(aspect_ratio=1) -1 1 2 3 1 2 3 4 Také v Sage existuje funkce plot. Její nejjednodušší použití vypadá takto: sage: plot(cos, (-5,5)) Předpokládám, že si každý dokáže představit obrázek, který vznikne, takže ho nemusím uvádět. Pro složitější funkce můžeme použít třeba něco jako sage: x=var(’x’) sage: plot(sqrt(x)*sin(x), (x,0,2*pi)) 1 2 3 4 5 6 -2 -1.5 -1 -0.5 0.5 1 28 1. PRVNÍ POKUSY Kromě plot je v Sage také funkce parametric plot, která pracuje podobně. Pomocí následujících příkazů získáme křivku zvanou epicykloida: sage: t=var(’t’) sage: x=2*cos(t)-cos(2*t) sage: y=2*sin(t)-sin(2*t) sage: parametric_plot((x,y),(t,0,2*pi)) -3 -2 -1 1 -2 -1 1 2 Čtenář si nyní může sám rozmyslet, jak by bylo možné nakreslit například čtverec, pro různě formulovaná zadání. Na závěr této části si předvedeme, jak lze pracovat s afinními transformacemi v rovině. Zkusme otočit a posunout výše nakreslenou elipsu. Nejprve určíme body na elipse podobně jako předtím. >> t=linspace(0,2*pi);x=3*cos(t);y=2*sin(t); Pak spočítáme matici otočení o daný úhel, v našem případě to bude 30◦ , tedy π/6. Maticí vynásobíme souřadnice bodů na elipse. Pokud souřadnice umístíme do matice tak, že každý sloupec bude tvořit souřadnice jednoho bodu, pro transformaci rotace celé elipsy nám stačí jedno maticové násobení. Výsledek pak můžeme zpátky rozdělit na xové a y-ové souřadnice. 5. GEOMETRIE V ROVINĚ 29 >> phi=pi/6; >> R=[cos(phi), -sin(phi);sin(phi), cos(phi)]; >> xyR=R*[x;y]; >> xR=xyR(1,:); >> yR=xyR(2,:); Nakonec všechny body posuneme o potřebnou délku podle souřadnic středu elipsy, v uvedeném příkladě to je [2,3]. Při tom můžeme využít konvenci Matlabu, že při přičítání skaláru k vektoru nebo k matici se přičte příslušná hodnota ke každé složce. >> xP=xR+2; >> yP=yR+3; >> plot(xP,yP) >> grid on −1 0 1 2 3 4 5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 A stejná konstrukce v Sage s tím, že maticové násobení vyjádříme ex- plicitně. sage: t=var(’t’) sage: phi=pi/6 sage: x=3*cos(phi)*cos(t)-2*sin(phi)*sin(t) sage: y=3*sin(phi)*cos(t)+2*cos(phi)*sin(t) sage: parametric_plot((x+2,y+3),(t,0,2*pi)) 30 1. PRVNÍ POKUSY 0 1 2 3 4 1 2 3 4 5 Soudím, že další afinní transformace by již čtenář zvládl bez větších problémů. KAPITOLA 2 Základy lineární algebry Kropáčku, vypadáte dnes nějak unaveně. Čím to je? Zkoušel jsem zjistit něco o maticích, jejich použití a historii, jak jste včera chtěl, pane profesore. To mě těší. A jak jste postupoval? Zadal jsem do googlu heslo Matrix a pak se probíral výsledky. A co jste zjistil? Je to fakt hrubý. Zvládl jsem za noc všechny tři díly. Z hlediska výpočetní techniky je lineární algebra jedním z nejlépe softwarově obhospodařovaným matematickým odvětvím. Existuje obrovská škála volně šiřitelných i komerčních programů a knihoven, které je možné použít pro řešení nejrůznějších typů obecných i specializovaných problémů, na něž můžeme v lineární algebře narazit. Je to dáno i tím, že nelezení řešení soustavy lineárních rovnic, což je typická úloha lineární algebry, patří mezi nejčastější inženýrské i ekonomické úlohy. Matlab patří v tomto směru ke špičce mezi numericky orientovaným software. Je to dáno jeho základní orientací od začátku jeho vývoje. Název Matlab je původně zkratka od Matrix Laboratory a v jeho dřevních dobách byly matice jediným datovým typem, přičemž číslo se považovalo za matici o jednom řádku a jednom sloupci. I třeba problém hledání kořenů polynomu se převádí na problém nalezení vlastních čísel matice, což, jak jsme viděli, ne vždy dává plně uspokojivé výsledky. K práci s maticemi a vektory jsme si již trochu přičichli v předchozí kapitole, nyní se pokusíme dané téma obsáhnout podrobněji. Trochu změníme způsob práce, nejprve se budeme věnovat Matlabu a až následně Sage. 1. Matice a vektory v Matlabu U následujících příkladů nebudu většinou uvádět jejich výsledek, jelikož předpokládám, že si je čtenář bude zkoušet sám a i v případě, že tak neučiní, dokáže si výsledek snadno představit. 31 32 2. ZÁKLADY LINEÁRNÍ ALGEBRY Matice v Matlabu můžeme zadat výčtem jejich prvků v hranatých závorkách, přičemž jednotlivé prvky na řádku oddělujeme mezerou nebo čárkou, jednotlivé řádky matice oddělujeme středníkem. Další funkcí středníku je, že pokud jej vložíme za příkaz, výsledek příkazu, např. přiřazení, se nezobrazí. >> A=[1 -1 2 -3; 3 0 4 5; 3.2, 5, -6 12] >> u=[1 2 3 4] >> v=[-1;-2;-3] Prázdnou matici vytvoříme příkazem >> Emp=[] Matice lze vytvářet pomoci už definovaných proměnných, je potřeba kontrolovat, aby souhlasily typy jednotlivých proměnných. >> x=pi/4 >> y=[x, 2*x, 3*x] >> B=[A; u] >> C=[A, v] Pro vytvoření matic větších rozměrů je možné použít některých funkci MATLABu. Příkaz >> Z=zeros(2,5) vytvoří nulovou matici příslušného typu, příkazem >> O=ones(3,4) vyrobíme matici že samých jedniček. Příkaz >> I=eye(5,8) nám dá jednotkovou matici, přičemž jedničky probíhají hlavní diagonálu. Všechny výše uvedené příkazy je možné žádat i s jedním parametrem, v tom případě je výsledkem čtvercová matice příslušného řadu. Funguje taky příkaz >> O=ones(0,5) Vektor, jehož prvky tvoří aritmetickou posloupnost, je možné vytvořit následujícím způsobem: >> x=1:10 vytvoří vektor postupně obsahující čísla 1 až 10. Pro jiný rozdíl mezi jednotlivými prvky než 1 lze použít např. příkaz 1. MATICE A VEKTORY V MATLABU 33 >> y=0:2:12 pro y se sudými čísly od 0 do 12 nebo >> z=0:0.01:2 pro z s čísly od 0 do 2 s krokem 0.01. Je možno použít i záporný krok: >> x1=10:-1:1 Kromě toho taky můžeme použít příkaz linspace >> x=linspace(a,b,n) kde a, b jsou první a poslední prvek vektoru x a n je počet prvků. Třetí parametr je nepovinný a pokud není uveden, bere se roven 100. Podobně je taky možné definovat vektor s desítkovou logaritmickou škálou: >> x=logspace(a,b,n) který je ekvivalentní příkazu >> x=10.^linspace(a,b,n) Zde je ovšem implicitně n=50. Rozměr matice zjistíme příkazem size. Na jednotlivé prvky matice je možné se odkázat pomoci kulatých závorek - tj. A(3,2) je prvek matice A ve 3. řádku a 2. sloupci. Toto vyjádření ovšem lze použít obecněji, kdy první parametr je vektor obsahující indexy vybraných řádku a druhý parametr je vektor sloupcových indexu. Pak >> A([1 3],[5 2]) je vybrána submatice z A tvořena 1. a 3. řádkem a 5. a 2. sloupcem. Je možno taky provádět přiřazení do takto vybraných prvků při zachování správných rozměrů: >> A([1 3],[5 2])=eye(2) Pokud chceme vybrat celý řádek nebo sloupec použijeme symbol ’:’, např. B(3,:) je 3. řádek matice B. Příkazem B(3,:)=[] vypustíme z matice B tento řádek. Zajímavé je, že nefunguje přiřazení >> o=[]; >> B(3,:)=o Matic stejných rozměrů lze sčítat a odčítat (+, −), matice vhodných rozměrů se dají násobit (*). Je také možné násobit konstantou nebo 34 2. ZÁKLADY LINEÁRNÍ ALGEBRY přičíst konstantu k matici, pak se přičte ke všem prvkům. Dělení jsou v MATLABu dvě – pravé a levé: ’\’, ’/’, přičemž výraz >> X=B\C dá řešení maticové rovnice B*X=C a výraz >> X=B/C je řešením maticové rovnice X*C=B. Symbolem ” ’ ” dostaneme hermitovskou transpozici matice, tj. matici transponovanou a komplexně sdruženou. Operátor ˆ slouží k umocňování čtvercových matic. Uvedené operátory mají také tzv. tečkové varianty, kdy se operace provádějí na jednotlivé prvky matice, např. >> A=B.*C vynásobí odpovídající prvky matic B a C, příkazem >> A.^2 se umocní všechny prvky matice A. Operátor ” .’ ” se používá pro transpozici matic. Na matice lze aplikovat taky všechny běžné funkce (sin, cos,. . . ), které se aplikují člen po členu. Pokud máme matici A typu m × n, pak pokus o určení hodnoty A(k,l), kde k> m nebo l> n, skončí chybou. Ovšem příkaz A(k,l)=1 přiřazení provede, přičemž chybějící členy jsou doplněny nulami. >> A=eye(2,3); >> A(4,5)=1 A = 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 Příkazem reshape je možné přeskládat matici do jiného typu. Za- dáním >> A=randi(9,3,4) A = 6 5 2 8 8 1 8 2 9 1 2 9 >> B=reshape(A,6,2) B = 6 2 1. MATICE A VEKTORY V MATLABU 35 8 8 9 2 5 8 1 2 1 9 vytvoříme z prvků matice A matici B o 6 řádcích a 2 sloupcích, přičemž prvky se při přeskládání berou po sloupcích. Pokud bychom chtěli matici A přeskládat po řádcích, museli bychom nejdřív matici A transponovat a pak opět transponovat výsledek, dostaneme tak ovšem matici o dvou řádcích a šesti sloupcích: >> B=reshape(A’,6,2)’ B = 6 5 2 8 8 1 8 2 9 1 2 9 Překlopit matici, tj. obrátit pořadí řádku nebo sloupců je možné provést příkazem flipud resp. fliplr. Také je možné matici „otočit o 90 stupňů příkazem rot90, přičemž další nepovinný parametr udává, kolikrát se má rotace provést. Všechny operátory pro práci s maticemi uvedené dříve mají své funkci ekvivalenty: funkce význam operátor plus plus + uplus unární plus + minus minus − uminus unární minus − mtimes maticové násobení * times násobení po prvcích .* mpower maticová mocnina ˆ power mocnina po prvcích .ˆ mldivide levé maticové dělení \ mrdivide právě maticové dělení / ldivide levé dělení po prvcích .\ rdivide právě dělení po prvcích ./ Příkaz a:b případně a:d:b lze nahradit příkazem colon(a,b) nebo colon(a,d,b). Dvojí úlohu má funkce diag. Jeho aplikaci na vektor získáme diagonální matici s argumentem na hlavní diagonále. Pokud ji použijeme na matici, funkce diag vybere z matice hlavní diagonálu a umístí ji do vektoru. Pokud chceme pracovat s jinou diagonálou než s hlavní, 36 2. ZÁKLADY LINEÁRNÍ ALGEBRY můžeme použít jako druhý parametr funkce diag číslo diagonály, přičemž kladná čísla se použijí nad hlavní diagonálou a záporná pod ní. Příklad: \begin{Matlab} >> A=randi(10,5,5) A = 3 3 2 10 1 1 8 7 10 5 2 5 7 5 4 6 5 3 7 10 4 9 5 9 3 >> x=diag(A) x = 3 8 7 7 3 >> B=diag(x,2) B = 0 0 3 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0 7 0 0 0 0 0 0 0 7 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 >> C=diag(pi,-3) C = 0 0 0 0 0 0 0 0 0 0 0 0 3.1416 0 0 0 >> d=diag(diag(A)) d = 3 0 0 0 0 0 8 0 0 0 0 0 7 0 0 0 0 0 7 0 0 0 0 0 3 1. MATICE A VEKTORY V MATLABU 37 Pomoci funkci tril a triu vyrobíme z daně matice dolní nebo horní trojúhelníkovou matici, přičemž je možné podobně jako u diag použít další nepovinný parametr. Příkaz max najde maximální prvek ve vektoru, přičemž pokud na výstupu uvedeme i druhý výstupní parametr, uloží se do něj index tohoto prvku. Při použití příkazu max na matici se hledají maximální prvky v jednotlivých sloupcích. Podobně funguje příkaz min: >> x=rand(1,5) x = 0.2785 0.5469 0.9575 0.9649 0.1576 >> M=max(x) M = 0.9649 >> [m,k]=min(x) m = 0.1576 k = 5 >> A=rand(4) A = 0.9706 0.1419 0.9595 0.9340 0.9572 0.4218 0.6557 0.6787 0.4854 0.9157 0.0357 0.7577 0.8003 0.7922 0.8491 0.7431 >> [m,k]=min(A) m = 0.4854 0.1419 0.0357 0.6787 k = 3 1 3 2 U komplexních matic se maximum či minimum hledá mezi absolutními hodnotami prvků. K setřídění vektoru (vzestupně) se používá příkaz sort, u komplexních hodnot se opět provádí třídění podle absolutních hodnot. Do druhého nepovinného výstupního parametrů je umístěna třídicí permutace indexu. Tj. pokud zadáme [y,ind]=sort(x) pak platí y=x(ind). Při použití příkazu sort na matice se třídí jednotlivé sloupce. Příkazy sum a prod sečtou nebo vynásobí všechny prvky vektoru. Aplikovaný na matici provedou příslušnou operaci po jednotlivých sloupcích. Např. faktoriál čísla n zjistíme snadno pomoci příkazu f=prod(1:n), přičemž příkaz funguje správné i pro n rovno 0, protože součin přes prázdnou matici je roven 1. 38 2. ZÁKLADY LINEÁRNÍ ALGEBRY Pokud bychom chtěli zjistit součet všech prvků matice, musíme příkaz sum použít dvakrát: >> sum(sum(A)) Příkazem size zjistíme rozměry matice, přičemž je možné jej použít ve formě s=size(a) nebo [m,n]=size(a). Příkaz length dává délku vektoru, přičemž pokud jej aplikujeme na matici, dostaneme maximální z rozměrů, tj. length(A) dává stejný výsledek jako max(size(A)). V MATLABu je implementováno velké množství funkci a algoritmu lineární algebry. Podrobný výpis lze získat pomoci nápovědy help matfun. Zde vyjmenujeme jen některé základní: det determinant matice rank hodnost matice null nulový prostor (jádro) matice norm maticová nebo vektorová norma trace stopa matice (součet diagonálních prvků) inv inverzní matice pinv pseudoinverzní matice lu LU rozklad qr QR rozklad svd singulární rozklad eig vlastní hodnoty a vektory matice poly charakteristický polynom matice, resp. vytvoření polynomu s danými kořeny Syntaxi a popis jednotlivých funkci zjistíme nejlépe pomoci nápovědy. Blíže se budeme věnovat jen funkci eig pro učení vlastních čísel a vlastních vektorů matice. Tato funkce má dva základní způsoby použití. Pokud je výstupní proměnná pouze jedna, vrátí funkce vektor vlastních hodnot dané ma- tice >> A=randi(10,3,3) A = 10 2 6 2 2 6 3 9 2 >> D=eig(A) D = 14.2792 5.1549 -5.4342 2. MATICE A VEKTORY V SAGE 39 V případě, že použijeme dvě výstupní proměnné, získáme i matici, jejíž sloupce jsou tvořeny vlastními vektory původní matice. Vlastní čísla už netvoří vektor ale hlavní diagonálu v diagonální matici: >> [V,D]=eig(A) V = 0.8126 0.7685 -0.2316 0.3573 -0.4241 -0.5725 0.4604 -0.4792 0.7865 D = 14.2792 0 0 0 5.1549 0 0 0 -5.4342 Pro takto uspořádané vlastní vektory a vlastní čísla by mělo platit A*V=V*D: >> A*V ans = 11.6032 3.9613 1.2584 5.1025 -2.1863 3.1110 6.5747 -2.4701 -4.2741 >> V*D ans = 11.6032 3.9613 1.2584 5.1025 -2.1863 3.1110 6.5747 -2.4701 -4.2741 Ve skutečnosti je mezi výsledky rozdíl řádově 10−14 . 2. Matice a vektory v Sage V Sage je pro vytváření vektorů určna funkce vector. Její parametrem je seznam prvků. sage: v=vector([1,2,3]);v (1, 2, 3) Takto se vytvoří vektor řádkový, v případě, že potřebujeme vektor sloupcový, musíme jej transponovat příkazem transpose, nebo ho vytvořit coby sloupcovou matici, což je popsáno o něco níže. Standardně se jedná o vektor celočíselný nebo racionální podle typu výrazů, kterými je vytvořen. Je možné ale při definici vnutit typ jiný pomocí prvního argumentu před samotným výčtem prvků. 40 2. ZÁKLADY LINEÁRNÍ ALGEBRY sage: v=vector([1,2,3]);v (1, 2, 3) sage: w=vector(RR,[1,2,3]);w (1.00000000000000, 2.00000000000000, 3.00000000000000) K jednotlivým prvkům vektoru se přistupuje pomocí hranatých závorek, prvky jsou indexovány od nuly, nikoliv od jedničky. sage: w[1] 2.00000000000000 Při vytváření vektorů je potřeba dát pozor na to, abychom nepoužili matlabovskou konstrukci sage: u=[2,3,-1];u [2, 3, -1] Výsledek vypadá stejně jako u vektoru, ale nejedná se o vektor, nýbrž seznam prvků, jak prozradí příkaz type: sage: type(v) sage: type(u) K vytváření matic používáme funkci Matrix, funguje i s počátečním písmenem malým. Nejjednodušší je vytvoření čtvercové matice, k čemuž stačí jeden parametr: sage: A=Matrix(3);A [0 0 0] [0 0 0] [0 0 0] Opět je možné vnutit vytvářené matici jiný typ: sage: B=Matrix(RR,2);B [0.000000000000000 0.000000000000000] [0.000000000000000 0.000000000000000] Pomocí dalších parametrů můžeme definovat obdélníkové matice a její prvky, nebo lze vytvořit matici pomocí seznamu jejich prvků: sage: B=Matrix(2,3);B [0 0 0] [0 0 0] 2. MATICE A VEKTORY V SAGE 41 sage: B=Matrix(2,3,[k/6 for k in range(6)]);B [ 0 1/6 1/3] [1/2 2/3 5/6] sage: B=Matrix(2,[k/6 for k in range(6)]);B [ 0 1/6 1/3] [1/2 2/3 5/6] sage: C=Matrix([[1,2,3],[4,5,6],[7,8,9]]);C [1 2 3] [4 5 6] [7 8 9] Povšimněme si, že druhý a třetí příkaz dávají stejný výsledek, pokud bychom ve druhém příkazu použili jiný počet sloupců než 3, došlo by k chybě. K jednotlivým prvkům přistupujeme opět pomocí hranatých závorek, přitom je možný dvojí způsob syntaxe. Nesmíme zapomenout na indexování od nuly: sage: B[1,1] 2/3 sage: B[1][2] 5/6 Rozměry matice A zjistíme pomocí A.nrows(), resp. A.ncols(). Jednotkovou matici získamé pomocí příkazu identity matrix, jejíž parametr udává rozměr matice (pouze čtvercové). Zdálo by se, že vektory můžeme ekvivalentně vytvořit jako matice o jednom řádku. Což je sice možné, ale operace pak fungují trochu jinak. V Sage je totiž možné násobit matici zprava řádkovým vektorem. Výsledkem je opět řádkový vektor, který má stejné prvky, které bychom dostali při obvyklém násobení sloupcovým vektorem: sage: A=Matrix([[1,2,3],[3,2,1],[1,1,1]]) sage: v=vector([1,2,3]) sage: A*transpose(v) [14] [10] [ 6] sage: A*v (14, 10, 6) Na první pohled to je sice trochu nezvyklé, ale po nějaké době používání nám to už zvláštní nepřijde. S řádkovou maticí ovšem tato operace nefunguje: 42 2. ZÁKLADY LINEÁRNÍ ALGEBRY sage: A=Matrix([[1,2,3],[3,2,1],[1,1,1]]) sage: B=Matrix([1,2,3]) sage: A*B ------------------------------------------------ TypeError atd. Pro řešení systému lineárních rovnic je možné použít levé dělení podobně jako v Matlabu, pravé dělení nefunguje. Alternativní metodou je použití funkce solve right nebo solve left: sage: A=Matrix([[1,2,3],[3,2,1],[1,-1,1]]) sage: b=vector([1,2,3]) sage: x=A\b;x (19/16, -9/8, 11/16) sage: x=A.solve_right(b);x (19/16, -9/8, 11/16) Funkce det a rank fungují tak, jak bychom čekali. Charakteristický polynom matice lze určit takto: sage: A=Matrix([[1,2,3],[3,2,1],[1,-1,1]]) sage: characteristic_polynomial(A) x^3 - 4*x^2 - 3*x + 16 Inverzní matici a stopu matice můžeme vyjádřit pouze následujícím způsobem: sage: A=Matrix([[1,2,3],[3,2,1],[1,-1,1]]) sage: A.inverse() [-3/16 5/16 1/4] [ 1/8 1/8 -1/2] [ 5/16 -3/16 1/4] sage: A.trace() 4 Podobně můžeme určit vlastní čísla a vlastní vektory matice: sage: A=Matrix([[1,2],[3,2]]) sage: A.eigenvalues() [4, -1] sage: A.eigenvectors_right() [(4, [ (1, 3/2) ], 1), (-1, [ 3. PŘÍKLADY 43 (1, -1) ], 1)] Možná, že poslední výsledek vypadá na první pohled poněkud nesrozumitelně, ale není to tak složité. Jedná se o seznam trojic, kde na první místě je vlastní číslo, pak následuje vlastní vektor a třetí v pořadí je násobnost. (Myslím, že přehlednější zformátování výsledku by neškodilo.) V Sage taky existuje funkce kernel, u níž bychom podle jména očekávali, že s její pomocí lze zjistit jádro matice jakožto lineárního zobrazení, přesněji řečeno jeho bázi. To sice možné je, ale tato funkce pracuje po řádcích: sage: A=Matrix([[1,2,3],[3,2,1],[1,1,1]]) sage: A.kernel() Free module of degree 3 and rank 1 over Integer Ring Echelon basis matrix: [ 1 1 -4] Vidíme, že první a druhý řádek spolu skutečně dávají čtyřnásobek třetího. Takže pro nalezení jádra lineárního zobrazení daného maticí bychom ji museli nejdříve transponovat. Závěrem této části můžete hádat či zkusit, co dostaneme jako výsledek příkazů A.rows() a A.columns(). 3. Příklady Příklad 9. Zkusme si nejprve demonstrovat v Matlabu úpravu matice na schodovitý tvar. Vezměme rozšířenou matici soustavy z příkladu 2.1: >> A=[0.5,0.125,0.2 270;0.25,0.75,0.2,270;... 0.25,0.125,0.6,270] A = 0.5000 0.1250 0.2000 270.0000 0.2500 0.7500 0.2000 270.0000 0.2500 0.1250 0.6000 270.0000 Vynásobení řádku konstantou odpovídá násobení zleva diagonální maticí, kde na příslušném řádku je potřebná konstanta. >> D=diag([2,4,4]) D = 2 0 0 0 4 0 0 0 4 44 2. ZÁKLADY LINEÁRNÍ ALGEBRY >> A1=D*A A1 = 1.0e+03 * 0.0010 0.0003 0.0004 0.5400 0.0010 0.0030 0.0008 1.0800 0.0010 0.0005 0.0024 1.0800 Všimněte si konstanty 1.0e+03, která je vytknuta před matici. Pokud bychom radši viděli zobrazení bez této konstanty, použijeme >> format short g >> A1 A1 = 1 0.25 0.4 540 1 3 0.8 1080 1 0.5 2.4 1080 První řádek vynásobíme -1 a přičteme ho ke druhému a třetímu. Tato úprava odpovídá násobení zleva maticí, která vychází z jednotkové a má navíc v prvním sloupci -1 ve druhém a třetím řádku >> T1=eye(3);T1(2,1)=-1;T1(3,1)=-1 T1 = 1 0 0 -1 1 0 -1 0 1 >> A2=T1*A1 A2 = 1 0.25 0.4 540 0 2.75 0.4 540 0 0.25 2 540 Druhý a třetí řádek opět vynásobíme 4 >> A3=diag([1,4,4])*A2 A3 = 1 0.25 0.4 540 0 11 1.6 2160 0 1 8 2160 přehodíme 2. a 3. řádek >> T3=eye(3);T3=T3([1,3,2],:) T3 = 1 0 0 0 0 1 3. PŘÍKLADY 45 0 1 0 >> A4=T3*A3 A4 = 1 0.25 0.4 540 0 1 8 2160 0 11 1.6 2160 a nakonec druhý řádek vynásobíme -11 a přičteme ho ke třetímu >> T4=eye(3);T4(3,2)=-11;A5=T4*A4 A5 = 1 0.25 0.4 540 0 1 8 2160 0 0 -86.4 -21600 Nakonec bychom ještě mohli vydělit pokrátit řádek, kvůli jednomu řádku je asi zbytečné vyrábět diagonální matici: >> A5(3,:)=A5(3,:)/A5(3,3) A5 = 1 0.25 0.4 540 0 1 8 2160 0 0 1 250 Pokud by se nyní někomu nechtělo ručně počítat jednotlivé neznáme, může se matici dál upravovat až na jednotkovou (přesněji její první tři sloupce), v tom případě mu v posledním sloupci vyjde řešení. Pokud bychom tuto soustavu řešili v Sage a zadávali matici podobně jako v Matlabu, dostali bychom desetinná čísla s patnácti číslicemi za desetinou tečkou, takže by obsahovala spoustu zbytečných nul. Proto lepe bude zadat matici jako racionální: sage: A=matrix(QQ,[[0.5,0.125,0.2,270], [0.25,0.75,0.2,270],[0.25,0.125,0.6,270]]);A [1/2 1/8 1/5 270] [1/4 3/4 1/5 270] [1/4 1/8 3/5 270] Sage ovšem umí provést uvedené úpravy sám, dokonce až do úplného konce sage: A.echelon_form() [ 1 0 0 400] [ 0 1 0 160] [ 0 0 1 250] 46 2. ZÁKLADY LINEÁRNÍ ALGEBRY takže hned vidíme řešení soustavy. Matlab tohle umí taky, příslušná funkce se jmenuje rref: >> rref(A) ans = 1 0 0 400 0 1 0 160 0 0 1 250 Ukazovat příklady na výpočet determinantu, hodnosti nebo inverzní matice pomocí Matlabu a Sage je patrně zbytečné. Názvy příslušných funkcí najdeme výše a jejich použití je snadné, bez nějakých záludností. Pevně věřím, že čtenář by dokonce našel několik způsobu (pšt nebo šest) jak vypočítat inverzní matici. Dobrat se řešení soustavy lineárních rovnic je také možné několika způsoby: sage: A=matrix(QQ,[[0.5,0.125,0.2],[0.25,0.75,0.2], [0.25,0.125,0.6]]);A [1/2 1/8 1/5] [1/4 3/4 1/5] [1/4 1/8 3/5] sage: b=vector([270,270,270]);b (270, 270, 270) sage: A.inverse()*b (400, 160, 250) sage: A\b (400, 160, 250) sage: A^(-1)*b (400, 160, 250) sage: A.solve_right(b) (400, 160, 250) V Matlabu lze použít podobně všechny způsoby kromě posledního, jediným drobným rozdílem je inv(A) místo A.inverse(). Oproti Sage ale v Matlabu existuje funkce linsolve, o níž je možné se více dozvědět v nápovědě. Příklad 10. Určit lineární závislost či nezávislost je poměrně jednoduché. Stačí příslušné vektory poskládat do matice a pak určit její hodnost. V Sage můžeme postupovat například takto: 3. PŘÍKLADY 47 sage: u1=vector([1,2,3]) sage: u2=vector([4,5,6]) sage: u3=vector([7,8,9]) sage: A=matrix([u1,u2,u3]);A [1 2 3] [4 5 6] [7 8 9] sage: rank(A) 2 Řešení v Matlabu je opět velmi podobné, takže ho nebudu uvádět. Problém nalezení jednoho z lineárně závislých vektorů jakožto lineární kombinace ostatních je problém hledání řešení soustavy lineárních rovnic, přičemž toto řešení nemusí být určeno jednoznačně. Zmíníme se o tom o něco níže. Nejdříve se ještě podívejme, jak nalézt bázi nějakého podprostoru, což je problém, který s lineární závislostí či nezávislostí úzce souvisí. Při hledání báze, můžeme vektory poskládat do matice a tuto pak upravovat na schodovitý tvar, čímž vynulujeme lineárně závislé vektory, nenulové řádky pak tvoří bázi. S úspěchem lze tedy použít funkce rref nebo echelon form. V Sage navíc existuje funkce basis, která je pro tento účel určena: sage: V = VectorSpace(QQ,3) sage: S = V.subspace([[1,2,0],[2,2,-1],[-1,0,1]]) sage: basis(S) [ (1, 0, -1), (0, 1, 1/2) ] Porovnáme-li výsledek s postupem uvedeným předtím, sage: u1=vector([1,2,0]) sage: u2=vector([2,2,-1]) sage: u3=vector([-1,0,1]) sage: A=matrix([u1,u2,u3]) sage: A.echelon_form() [ 1 0 -1] [ 0 2 1] [ 0 0 0] 48 2. ZÁKLADY LINEÁRNÍ ALGEBRY vidíme, že jediný rozdíl je v násobení druhého bázového vektoru tak, aby měl první koeficient roven jedné. V Matlabu lze použít funkci orth, která hledá ortonormální bázi prostoru generovaného sloupci matice. >> A=reshape(1:9,3,3) A = 1 4 7 2 5 8 3 6 9 >> Q=orth(A) Q = -0.4797 0.7767 -0.5724 0.0757 -0.6651 -0.6253 >> Q’*Q ans = 1.0000 -0.0000 -0.0000 1.0000 Tato funkce však evidentně nepracuje Grammovou–Schmidtovou ortogonalizací, to by musel být první sloupce matice Q násobek prvního sloupce matice A. Když jsme se dostali k pojmu ortonormální, měli bychom se vrátit ke skalárnímu součinu. V Sage je pro skalární součin určena operace ’*’: sage: u=vector([1,2,3]) sage: v=vector([-1,2,1]) sage: s1=u*v;s1 6 Pokud bychom chtěli použít maticové násobení, je to možné, ale výsledek má jiný typ: sage: s2=u*transpose(v);s2 (6) sage: type(s1) sage: type(s2) 3. PŘÍKLADY 49 V Matlabu se pro skalární součin používá funkce dot, ale maticové násobení dává stejný výsledek: >> u=[1,2,3];v=[-1,2,1]; >> dot(u,v) ans = 6 >> u*v’ ans = 6 Příklad 11. Zaměřme se na postup při hledání ortogonálního doplňku nějakého podprostoru, který je určen svou bází, Také pro ortogonální doplněk stačí nalézt jeho bázi. Hledáme tedy nezávislá řešení homogenní soustavy rovnic A*x=0, kde matice A je tvořena řádky báze našeho podprostoru. Takže hledaný ortogonální doplněk je vlastně jádrem lineárního zobrazení daného maticí A. Teď už tedy problém vyřešíme snadno. Příklad řešení v Sage: sage: u=vector([1,2,-1,4]) sage: v=vector([-1,2,0,3]) sage: A=matrix([u,v]) sage: B=A.transpose().kernel().basis() sage: B [ (1, 2, 1, -1), (0, 3, -2, -2) ] A pak ještě můžeme ověřit, že odpovídající vektory jsou opravdu kolmé: sage: C=matrix([B[0],B[1]]);C [ 1 2 1 -1] [ 0 3 -2 -2] sage: A*transpose(C) [0 0] [0 0] A řešení stejné úlohy v Matlabu: >> u=[1,2,-1,4]; >> v=[-1,2,0,3]; >> A=[u;v]; 50 2. ZÁKLADY LINEÁRNÍ ALGEBRY >> B=null(A) B = 0.3707 -0.2597 -0.0686 -0.8403 0.9106 -0.0456 0.1693 0.4737 >> B’*B ans = 1.0000 0.0000 0.0000 1.0000 >> A*B ans = 1.0e-15 * 0 0.2220 0 0 Vidíme, že nalezená báze je ortonormální a ve výsledku je malá numerická chyba. V podobných případech je zpravidla na uživateli, aby posoudil, jestli výsledek má být ve skutečnosti nulový a případně provedl nějaké dodatečné úpravy. Příklad 12. Na závěr této části si ukážeme, jak si poradit se systémem, který má více řešení. (Případ, kdy systém nemá řešení poznáme lehce, když porovnáme hodnost matice soustavy a hodnost matice rozšířené soustavy.) Můžeme samozřejmě matici upravit na schodovitý tvar a pak hledat řešení ručně: sage: A=matrix([[1,2,3,1],[4,5,6,1],[7,8,9,1]]);A [1 2 3 1] [4 5 6 1] [7 8 9 1] sage: A.echelon_form() [1 2 3 1] [0 3 6 3] [0 0 0 0] Z druhého řádku dostávám y + 2z = 1, můžeme najít jedno řešení dosazením hodnoty za y nebo z, případně obecné řešení v parametrickém tvaru, pokud dosadíme třeba y = t. Toto řešení umí Sage najít i sám, bohužel se mu musí zadat všechny tři rovnice explicitně: sage: x,y,z=var(’x,y,z’) sage: solve([x+2*y+3*z==1,4*x+5*y+6*z==1,7*x+8*y+9*z==1], 3. PŘÍKLADY 51 x,y,z) [[x == r1 - 1, y == -2*r1 + 1, z == r1]] Najít jedno (tzv. partikulární) řešení můžeme v Sage nejméně dvěma způsoby: sage: A=matrix([[1,2,3],[4,5,6],[7,8,9]]);A [1 2 3] [4 5 6] [7 8 9] sage: b=vector([1,1,1]);b (1, 1, 1) sage: A.solve_right(b) (-1, 1, 0) sage: A\b (-1, 1, 0) Pokud ještě určíme bázi jádra matice A, sage: kernel(A.transpose()).basis() [ (1, -2, 1) ] tak snadno určíme libovolné řešení v parametrickém tvaru, pakliže víme, že každé takové řešení dostaneme z partikulárního přičtením lineární kombinace báze jádra matice A. V Matlabu se dá určit partikulární řešení také nejméně dvěma způ- soby: >> A=[1,2,3;4,5,6;7,8,9],b=[1;1;1], A = 1 2 3 4 5 6 7 8 9 b = 1 1 1 >> A\b Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.541976e-18. ans = -2.5 4 52 2. ZÁKLADY LINEÁRNÍ ALGEBRY -1.5 >> pinv(A)*b ans = -0.5 6.9389e-18 0.5 Všimněme si jednak varování u prvního výpočtu, jednak malé druhé složky řešení u druhého výsledku, které by měla být ve skutečnosti nulová a do třetice toho, že výsledky nejsou stejné. První výsledek se počítá pomocí Gaussovy eliminace (viz help mldivide), druhý pomocí tzv. pseudoinverzní matice, čímž dostaneme řešení s nejmenší normou (velikostí). K určení báze jádra matice slouží funkce null, o které už jsme se zmiňovali. Standardně bázi normuje na velikost jednotlivých vektorů rovnu jedné, pokud chceme racionální složky bázových vektorů, můžeme použít další parametr ’r’: >> A=[1,2,3;4,5,6;7,8,9]; >> null(A) ans = -0.40825 0.8165 -0.40825 >> null(A,’r’) ans = 1 -2 1 Obecné řešení tedy opět získáme jako ve tvaru u0 +t[1, −2, 1], kde u0 je partikulární řešení. Současně jsme ověřili známou věc, že totiž obecné řešení lze zapsat vícerým způsobem. Také si čtenář může rozmyslet, pro jaké hodnoty t dostaneme z jednoho z uvedených partikulárních řešení to druhé. KAPITOLA 3 Lineární modely a maticový počet Předpokládám, že co se týče grafické prezentace diferenčních rovnic, plně stačí to,co bylo uvedeno v úvodu. Také zkoušet si počítat obecné či partikulární řešení podle vztahů odvozených v učebnici není z hlediska tohoto textu příliš zajímavé. Podívejme se tedy na to jakým způsobem můžeme pracovat s lineárními modely za pomocí vhodného softwaru. Budeme používat Matlab, kde je práce s maticemi jednodušší, ale všechny uvedené výpočty by samozřejmě šly provádět také v Sage či jiném software, který zvládá násobení matic. 1. Populační modely Podívejme se na příklad 3.15 s rybami v rybníce. Rybáře samozřejmě může zajímat, za jakou dobu po nasazení ryb do rybníka se jejich populace ustálí. Můžeme předpokládat, určité počáteční množství plůdků, na kterém příliš nezáleží, můžeme třeba začít s číslem 1 a chápat to jako jednu populaci. Nebo začít hodnotou 100 a chápat výsledky v procentech oproti výchozímu množství. Přitom použijeme hodnotu τ = 0.1, která zajišťuje stagnaci populace. >> tau=0.1; >> M=[0 3 3; 0.2 0 0;0 0.6 tau]; >> p=[100;0;0]; >> p=M*p p = 0 20 0 >> p=M*p p = 60 0 12 Vidíme, že po dvou obdobích máme už jen 60% plůdků oproti původnímu množství, přičemž do dospělosti přežilo jen 12 % ryb. Na první 53 54 3. LINEÁRNÍ MODELY A MATICOVÝ POČET pohled by se mohlo zdát, že populace vymře. Podívejme se však na další vývoj. >> p=M*p p = 36.0000 12.0000 1.2000 >> p=M*p p = 39.6000 7.2000 7.3200 Počet plůdků ještě jednou klesl, ale pak se zvýšil, počet dospělých ryb také kolísá. Pojďme dále: p = 43.5600 7.9200 5.0520 >> p=M*p p = 38.9160 8.7120 5.2572 Vypadá to, že kolísání se zmenšuje. Po dalších čtyřech iteracích dostaneme hodnoty p = 40.8802 8.1217 5.5534 a ty se pak už prakticky nemění. Velmi dobrou aproximaci limitního výsledku dostaneme např. jako >> p=M^100*p p = 40.9091 8.1818 5.4545 1. POPULAČNÍ MODELY 55 Mohli bychom také zkoumat, jak rychle se ryby přemnoží, pokud bychom nenasadili štiky redukující jejich počet (τ = 1) případně nasadili menší množství štik (0.1 < τ < 1). Věřím, že příslušné výpočty by čtenář jistě zvládl sám. Podívejme se radši na trochu modifikovaný model. Předpoklad, že ryby by bez štik přežívaly pořád, je jistě nereálný. Přidejme tedy ještě další dvě populační stadia: všechny dospělé ryby se dostanou do stadia starých ryb, které mají omezenou plodnost, z nich se polovina dostane do stadia velmi starých ryb, které už jsou neplodné a všechny vymírají. Model takového procesu je popsán maticí >> M=[0 3 3 1 0; 0.2 0 0 0 0;0 0.6 0 0 0;... 0 0 0.5 0 0;0 0 0 0.5 0] M = 0 3.0000 3.0000 1.0000 0 0.2000 0 0 0 0 0 0.6000 0 0 0 0 0 1.0000 0 0 0 0 0 0.5000 0 a opět se můžeme podívat na chování tohoto modelu: >> p=[100;0;0;0;0]; >> p=M*p p = 0 20 0 0 0 >> p=M*p p = 60 0 12 0 0 >> p=M*p p = 36 12 0 12 56 3. LINEÁRNÍ MODELY A MATICOVÝ POČET 0 >> p=M*p p = 48.0000 7.2000 7.2000 0 6.0000 >> p=M*p p = 43.2000 9.6000 4.3200 7.2000 0 Zatím je těžko usuzovat na limitní chování. Po dalších deseti etapách máme >> p=M^10*p p = 61.8795 12.0144 6.9872 6.7915 3.2846 a po dalších 10 etapách >> p=M^10*p p = 83.7110 16.2442 9.4566 9.1753 4.4511 Vidíme, že počet ryb se zvětšuje, ale zdaleka ne tak dramaticky, jak u původního modelu. Nyní je možné do modelu zakomponovat nasazení dravců (štik), kteří by redukovali stejnou (nebo i různou) mírou jednotlivé populace. Teoretické i praktické zkoumání takového modelu přenecháme čtenáři k samostatné práci. 2. MARKOVOVY ŘETĚZCE 57 2. Markovovy řetězce U Markovových řetězců je vždy potřeba si rozmyslet, jaké mohou být všechny stavy, do nichž se daný proces může dostat. Tyto stavy musí být alespoň 2 (např. u součástky funguje/nefunguje), přičemž výpočet pravděpodobnosti přechodu mezi jednotlivými stavy může být různě složitý, ale tahle otázka patří do jiné kapitoly. Zkusme se podívat na příklad 3.24 se studenty na přednášce. Zajímá nás, jak se počty studentů v jednotlivých skupinách mění s časem. Na začátku uvažujme stočlenou skupinu vnímajících studentů. >> M=[0.5 0.1 0;0.4 0.5 0;0.1 0.4 1] M = 0.5000 0.1000 0 0.4000 0.5000 0 0.1000 0.4000 1.0000 >> p=[100;0;0]; >> p=M^2*p p = 29.0000 40.0000 31.0000 Po dvou hodinách je na přednášce skoro 70 studentů, to není marné. >> p=M^2*p p = 12.4100 23.2000 64.3900 Po dalších dvou hodinách už je přítomno jen o něco více než třetina studentů, nejvyšší čas přednášku ukončit. Někdy se také můžeme setkat s procesy, které ve skutečnosti nejsou Markovovy, ačkoliv tak na první pohled vypadají. Základním požadavkem na Markovův proces je, že pravděpodobnost přechodu z jednoho stavu do jiného nezávisí na tom, jak se proces do daného stavu dostal. Kdyby toto neplatilo, není možné konstruovat přechodové matice. Ukažme si to na příkladu s ruletou (3.25). Častou strategií hazardních hráčů je sázka na některou barvu (např. černou) a v případě výhry nechá hráč sázku včetně výhry na stejné barvě, přičemž tento postup může zopakovat vícekrát. Předpokládejme tedy sázku ve výši jednoho žetonu, hráč končí sérii, pokud získá 8 žetonů, poté začíná znovu s jedním žetonem. V případě prohry začíná 58 3. LINEÁRNÍ MODELY A MATICOVÝ POČET taktéž znovu s jedním žetonem, pokud ještě nějaký má. Nechť dále má hráč na začátku n žetonů, hru končí v okamžiku, kdy množství alespoň zdvojnásobí nebo všechno prohraje. Pro n = 1 dostáváme zpočátku obdobu příkladu 3.25. Pro větší n je ovšem množina možných stavů větší něž v uvedeném příkladu, protože je možné dosáhnout libovolné částky od 0 do 2n. Označíme-li Sk stav, že hráč má k žetonů (k = 0 . . . 2n, přičemž S2n je stav kdy hráč má alespoň 2n žetonů a končí hru), můžeme se ptát např. jaká je pravděpodobnost přechodu od stavu Sk do stavu Sk+1. Tato pravděpodobnost je 18/37, pokud hráč předtím prohrál a sází jeden žeton, ale je nulová, pokud hráč předtím vyhrál dva žetony a znovu je sází, protože v tom případě se může dostat pouze do stavů Sk+2 nebo Sk−2. Naskýtá se otázka, jestli a jak se dá počítat úspěšnost takové strategie. Jenou z možností je celý proces simulovat, tedy generovat náhodně výsledky rulety a počítat počet relativní četnost úspěšnosti mezi všemi simulacemi. Pro velký počet simulací se tato relativní četnost blíží pravděpodobnosti úspěchu. Při tom je možné taky sledovat další veličiny, např. průměrnou délku celého procesu do výhry či do prohry, její závislost na n atd. 3. Maticové rozklady a spol Co se týče maticových rozkladů, je na tom Matlab daleko lépe než Sage, kde je nabídka prakticky nulová. Vzhledem k tomu, že Sage se neustále vyvíjí, dá se očekávat, že tato mezera bude v budoucnu dopl- něna. Nejprve si ukažme, jak je možné spočítat LU rozklad. Ten se v softwarových balících skutečně používá, a to např. při výpočtu inverzní matice nebo determinantu. Ukažme se nejprve pro jednoduchost tento rozklad bez výměny řádků. Hledáme tedy dolní trojúhelníkovou matici L a horní trojúhelníkovou matici U takové, že platí A = L · U. Při Gaussově eliminaci nulujeme prvky matice pod hlavní diagonálou, jejím výsledkem je hledaná matice U. Bez výměny řádků provádíme prakticky jedinou elementární operaci, tou je přičtení násobku jednoho řádku k jinému. Nechť tedy mám zpracováno prvních k − 1 sloupců, tím jsme z původní matice A získali upravenou matici A(k) , která již má v prvních k−1 sloupcích nuly pod hlavní diagonálou. Další postup je tedy takový, že k-tý řádek násobíme příslušnými vhodnými koeficienty a přičítáme k následujícím. Dá se snadno zjistit že ony koeficienty mají hodnoty −a (k) i,k /a (k) k,k, pro i > k. To ovšem odpovídá násobení 3. MATICOVÉ ROZKLADY A SPOL 59 matice A(k) zleva maticí                  1 · · · 0 ... 1 0 1 ... ... ... 0 1 −li,k 1 0 1 ... ... 0 0 1                  pro li,k = a (k) i,k /a (k) k,k a vynulování celé příslušné části k-tého sloupce maticí Lk =           1 · · · 0 ... ... 1 −lk+1,k ... ... ... 0 −ln,k 1           . Celou úpravu matice tak můžeme vyjádřit jako U = Ln−1 · · · · · L2 · L1 · A odkud dostáváme A = L−1 1 · · · · · L−1 n−1 · U = L · U. Není ovšem problém ověřit, že L−1 k =           1 · · · 0 ... ... 1 lk+1,k ... ... ... 0 ln,k 1           60 3. LINEÁRNÍ MODELY A MATICOVÝ POČET a L =           1 · · · 0 l2,1 ... ... 1 ... ... ... ln,1 . . . ln,n−1 1 ,           kde tedy li,k = a (k) i,k /a (k) k,k, přičemž a (k) i,k je prvek na dané pozici v již částečně upravené matici A. Praktický výpočet provádíme většinou tak, že si matici L předvyplníme jedničkami na hlavní diagonále a nulami nad ní a pod hlavní diagonálu postupně doplňujeme hodnoty koeficientů, které používáme při násobení daného řádku, při jeho odečítání od řádku jiného. Pokud potřebujeme při Gaussově eliminaci vyměnit řádky, děláme ve skutečnosti LU rozklad matice P · A, kde P je takzvaná permutační matice která výměnu řádků zprostředkuje. Je to vlastně jednotková matice, v níž přeházíme řádky podle potřeby. Problémem je, že dopředu nevíme, jakým způsobem máme řádky vyměnit. Naštěstí je možné informaci o pořadí řádků vytvářet a měnit i během vlastního výpočtu. Přitom nepotřebujeme znát celou matici P, stačí permutace, která příslušné pořadí řádků určuje. Na začátku Gaussovy eliminace si vytvoříme pomocný vektor p = (1, . . . , n) a při každé výměně řádků ve zpracovávané matici vyměníme i prvky v p na stejných pozicích, jako jsou vyměňované řádky. Po ukončení eliminace nám prvek pi dává číslo řádku původní matice, který je po výměně řádků na pozici i-té. Tedy např. když vyjde p = (2, 4, 1, 3) máme LU rozklad matice, která je postupně tvořena druhým, čtvrtým, prvním a třetím řádkem původní matice A. Výměna řádků při Gaussově eliminace se ve skutečnosti dělá téměř vždy. Je to nejen kvůli tomu, abychom se ve zlomcích a (k) i,k /a (k) k,k vyhnuli dělení nulou, ale hlavně kvůli tomu, že dělení malým číslem je z numerického hlediska méně přesné než dělení číslem velkým. Pokud si tedy můžeme vybrat, snažíme se ve jmenovateli použít číslo co možná největší. A při Gaussově eliminaci si vybrat můžeme. V upravované matici můžeme najít mezi k-tým až n-tým řádkem takový, pro nějž je a (k) i,k v absolutní hodnotě největší a ten pak vyměníme s řádkem k-tým. Tomu největšímu prvku se říká vedoucí nebo hlavní prvek (anglicky pivot) a tento postup se nazývá částečný výběr hlavního prvku nebo 3. MATICOVÉ ROZKLADY A SPOL 61 částečná pivotáž. Existuje i úplná pivotáž, kdy vedoucí prvek nehledáme v jednom sloupci ale v celém zbytku upravované matice, v tom případě ovšem se nám např. při řešení systému lineárních rovnic změní i pořadí hledaných hodnot řešení. V Matlabu se pro LU rozklad používá funkce lu: >> A=randi(10,4,4) A = 5 6 8 1 4 10 4 6 9 3 6 8 6 8 1 10 >> [L,U]=lu(A) L = 0.5556 0.5000 1.0000 0 0.4444 1.0000 0 0 1.0000 0 0 0 0.6667 0.6923 -0.9808 1.0000 U = 9.0000 3.0000 6.0000 8.0000 0 8.6667 1.3333 2.4444 0 0 4.0000 -4.6667 0 0 0 -1.6026 Vidíme, že L není dolní trojúhelníková, potřebovali bychom ještě vyměnit první řádek se třetím. Naštěstí je možné požadovat ještě třetí výstupní parametr obsahující informaci o výměně řádků a je možné ho získat jako permutační matici nebo jen ve formě vektoru: >> [L,U,P]=lu(A) L = 1.0000 0 0 0 0.4444 1.0000 0 0 0.5556 0.5000 1.0000 0 0.6667 0.6923 -0.9808 1.0000 U = 9.0000 3.0000 6.0000 8.0000 0 8.6667 1.3333 2.4444 0 0 4.0000 -4.6667 0 0 0 -1.6026 P = 0 0 1 0 0 1 0 0 62 3. LINEÁRNÍ MODELY A MATICOVÝ POČET 1 0 0 0 0 0 0 1 >> [L,U,p]=lu(A,’vector’) L = 1.0000 0 0 0 0.4444 1.0000 0 0 0.5556 0.5000 1.0000 0 0.6667 0.6923 -0.9808 1.0000 U = 9.0000 3.0000 6.0000 8.0000 0 8.6667 1.3333 2.4444 0 0 4.0000 -4.6667 0 0 0 -1.6026 p = 3 2 1 4 Alespoň stručně se zmíníme o Choleského rozkladu. Ten se používá pro symetrické matice a v tomto případně se snažíme reálnou matici A vyjádřit ve tvaru součinu A = TT · T, kde T je horní trojúhelníková matice. Jestliže si označíme tij prvky matice T, formálně provedeme součin TT · T a výsledek porovnáme s maticí A, dostaneme následující vztahy: t11 = √ a11 t1j = a1j t11 , j = 2, . . . , n tii = aii − i−1 l=1 t2 li, i = 2, . . . , n tij = 1 tii aij − i−1 l=1 tlitlj pro j > i tij = 0 pro i > j. Při výpočtu je ovšem možné se setkat s odmocninami ze záporných čísel, takže matice T může být komplexní. Choleského rozklad existuje pro matice, jejichž všechny hlavní minory (determinanty tvořené prvními k sloupci a řádky matice pro k = 1, . . . , n) jsou nenulové. Pokud jsou dokonce kladné, je matice A pozitivně definitní a T je pak reálná. V Matlabu se pro Choleského rozklad používá funkce chol, která ovšem funguje jen pro pozitivně definitní matice. 3. MATICOVÉ ROZKLADY A SPOL 63 Další používané rozklady jsou QR rozklad a singulární rozklad, pro které jsou v Matlabu funkce qr a svd. Jejich použití je celkem jednoduché a čtenář je jistě snadno pochopí z nápovědy. Singulární rozklad se používá pro výpočet pseudoinverzní matice, u které se trochu zdržíme. Tato matice, často také nazývaná Moorova– Penrosova pseudoinverze, se zpravidla definuje pomocí Penrosových axiomů: Matice X se nazývá pseudoinverzní maticí matice A, jestliže platí A · X · A = A X · A · X = X (A · X)∗ = A · X (X · A)∗ = X · A Z těchto axiomů se dá dovodit, že taková matice může být nejvýše jedna a ze singulárního rozkladu zase plyne její existence. Označení pseudoinverzní matice není ustálené, kromě A(−1) se často používá A+ . Pseudoinverzní matice se zpravidla počítá s použitím singulárního rozkladu, pro matice, jejichž hodnost je rovna počtu sloupců ale platí A(−1) = (A∗ · A)−1 · A∗ a podobně pro matice s hodnosti rovnou počtu řádků A(−1) = A∗ ·(A·A∗ )−1 , což se dá snadno ověřit přímým výpočtem. Pseudoinverzní matice hraje důležitou roli při hledání přibližných řešení systémů lineárních rovnic, u nichž přesné řešení neexistuje. Nejlepší přibližné řešení (někdy též řešení ve smyslu nejmenších čtverců) soustavy A·x = b je takový vektor x který minimalizuje normu rozdílu obou stran rovnice, tedy A·x−b . Řešení ve smyslu nejmenších čtverců není určeno jednoznačně, ale dá se ukázat, že je řešením tzv. systému normálních rovnic pro danou soustavu, který má tvar A∗ · A · x = A∗ · b Takový systém rovnic je řešitelný vždy a jediné řešení má právě tehdy, když jsou sloupce matice A lineárně nezávislé, tedy její hodnost je rovna počtu sloupců. Pak ale je matice A∗ · A regulární a řešení ve smyslu nejmenších čtverců můžeme zapsat jako ˆx = (A∗ · A)−1 · A∗ · b = A(−1) · b. Platí ale ještě silnější tvrzení: Vektor ˆx = A(−1) · b je vždy řešením ve smyslu nejmenších čtverců soustavy rovnic A · x = b a pokud existuje jiné řešení ˜x ve smyslu nejmenších čtverců této soustavy, pak ˆx < ˜x . Podívejme se, jak bychom v Matlabu mohli řešit příklad na hledání přibližného řešení z učebnice. Jedná se o parabolu s vrcholem 64 3. LINEÁRNÍ MODELY A MATICOVÝ POČET v počátku, která nejlépe prokládá naměřená data. Nejprve si zadáme matici soustavy a pravou stranu: >> xx=1:10; >> A=[xx;xx.^2]’ A = 1 1 2 4 3 9 4 16 5 25 6 36 7 49 8 64 9 81 10 100 >> b=[1.44 10.64 4.48 14.56 31.12 39.20 54.88... 71.28 85.92 104.16]’; Protože x používáme pro nezávislou proměnou, označme hledané koeficienty y. Spočítat je můžeme pomocí pseudoinverzní matice, pro niž se v Matlabu používá funkce pinv: >> y=pinv(A)*b y = 0.6112 0.9981 Protože se vlastně jedná o koeficienty polynomu, můžeme si jej vytvořit příkazem >> P=[y(2) y(1) 0] P = 0.9981 0.6112 0 a nakonec zobrazit původní data i proloženou parabolu: >> x=0:0.01:10; >> yP=polyval(P,x); >> plot(xx,b’,’*r’,x,yP,’b’) 3. MATICOVÉ ROZKLADY A SPOL 65 0 1 2 3 4 5 6 7 8 9 10 0 20 40 60 80 100 120 Pro prokládání dat polynomem se dá v Matlabu použít také funkce polyfit. Tato funkce hledá koeficienty polynomu, který nejlépe aproximuje zadaná data. Použití je jednoduché, na vstupu jsou zadaná data a stupeň polynomu: >> P1=polyfit(xx,b’,2) P1 = 0.9900 0.7129 -0.2680 >> yP1=polyval(P1,x); >> plot(xx,b’,’*r’,x,yP1,’b’) 0 1 2 3 4 5 6 7 8 9 10 −20 0 20 40 60 80 100 120 Všimněte si, že vrchol paraboly v tomto případě neleží v počátku souřadného systému, nicméně aproximace dat by měla být lepší než v prvním případě. KAPITOLA 4 Geometrie Soudím, že zkoušet si geometrické výpočty v Matlabu nebo v Sage celkem zbytečné, protože se jedná o aplikování již probraných postupů. A jelikož o znázorňování objektů v rovině už také byla řeč, řekneme si něco o kreslení prostorových objektů. Základem je zobrazení grafu funkce dvou proměnných. V Matlabu je k tomu možné použít funkci mesh nebo surf. jejich použití je podobné, důležité je ovšem umět správně spočítat hodnoty grafu. Zkusme například zobrazit graf funkce cos(x2 +y2 ) pro x i y z intervalu [−1, 1]. Hodnoty funkce počítáme na husté čtvercové síti, přičemž k vytvoření této sítě můžeme použít funkce meshgrid: >> x=-1:0.05:1; >> y=-1:0.05:1; >> [X,Y]=meshgrid(x,y); >> X(1:3,1:5) ans = -1.0000 -0.9500 -0.9000 -0.8500 -0.8000 -1.0000 -0.9500 -0.9000 -0.8500 -0.8000 -1.0000 -0.9500 -0.9000 -0.8500 -0.8000 >> Y(1:3,1:5) ans = -1.0000 -1.0000 -1.0000 -1.0000 -1.0000 -0.9500 -0.9500 -0.9500 -0.9500 -0.9500 -0.9000 -0.9000 -0.9000 -0.9000 -0.9000 V proměnné X je každý řádek x zopakovaný tolikrát, jaký je rozměr y, pro Y to platí obdobně po sloupcích. Nyní můžeme snadno vypočítat potřebné funkční hodnoty a graf nakreslit pomocí funkce mesh >> z=cos(X.^2+Y.^2); >> mesh(x,y,z) 67 68 4. GEOMETRIE −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 −0.5 0 0.5 1 nebo pomocí funkce surf >> surf(x,y,z) −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 −0.5 0 0.5 1 Taky je možné použít funkci waterfall, která se hodí v případě, že chceme vidět, jak se mění graf funkce jedné proměnné v závislosti na nějakém parametru. Mějme třeba funkci sin x+p·cos 2x, kde p se mění od 0 do 1. Dobrý náhled na chování funkce v závislosti na parametru p získáme v Matlabu takto: >> x=0:0.05:2*pi; >> yy=[]; >> for p=pp, y=sin(x)+p*cos(2*x); yy=[yy;y]; end >> waterfall(x,pp,yy) 4. GEOMETRIE 69 0 1 2 3 4 5 6 7 0 0.5 1 −2 −1.5 −1 −0.5 0 0.5 1 1.5 Šlo by samozřejmě výpočet provést bez cyklu s použitím funkce meshgrid podobně jako v předchozím příkladu. V Sage je pro tento účel možné použít funkci plot3d sage: x, y = var(’x, y’) sage: plot3d(cos(x^2 + y^2), (-2,2), (-2,2)) nebo parametric plot3d pro spírálu sage: parametric_plot3d((sin(u),cos(u),u/10),(u,0,20)) 70 4. GEOMETRIE nebo zakroucenou pneumatiku: sage: u, v = var(’u,v’) sage: fx = (3+sin(v)+cos(u))*cos(2*v) sage: fy = (3+sin(v)+cos(u))*sin(2*v) sage: fz = sin(u)+2*cos(v) sage: parametric_plot3d([fx,fy,fz],(u,0,2*pi),(v,0,2*pi), color="red")