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).
1
Obsah
1 První pokusy 5
1.1 Praktické výpočty . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Kombinatorické výpočty . . . . . . . . . . . . . . . . . . . . . 13
1.3 Diferenční rovnice . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.4 Pravděpodobnost . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.5 Geometrie v rovině . . . . . . . . . . . . . . . . . . . . . . . . 25
2 Základy lineární algebry 34
2.1 Matice a vektory v Matlabu . . . . . . . . . . . . . . . . . . . 35
2.2 Matice a vektory v Sage . . . . . . . . . . . . . . . . . . . . . 44
2.3 Příklady . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3 Lineární modely a maticový počet 60
3.1 Populační modely . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.2 Markovovy řetězce . . . . . . . . . . . . . . . . . . . . . . . . 64
3.3 Maticové rozklady a spol . . . . . . . . . . . . . . . . . . . . . 66
4 Geometrie 75
2
Ú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ě.
Pokud by čtenáře zajímaly zdroje, kde se může dozmědět více, co se týče
software je nejlepší hledat v dokumentaci (http://www.mathworks.com/help/
a http://www.sage.org/). Pokud se čtenář bude zajímat o teoretickou stránku
3
problémů, většinu z nich pokrývá libovolny kurs lineární algebry (např. [2]),
o pseudoinverzi se můžeme dozvědět více třeba v kalsické učebnici statistiky
– [1].
4
Kapitola 1
První pokusy
1.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í
1
Viz též Douglas Adams: Stopařův průvodce po Galaxii
5
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:
>> 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
6
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:
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:
7
>> 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()
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
8
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
-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
,
9
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
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:
10
>> 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)
11
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)
12
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.
1.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
>>
13
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 =
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
14
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 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)
.
15
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
>>
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
16
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 =
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],
17
[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.
1.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,’*’);
18
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,’*’);
19
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;
>> 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,’+’);
20
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:
21
>> n=12;d_k=1;m=1;
>> for k=1:n,d_k=(1+m/n)*d_k;end
>> d_k
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í?
1.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 kom-
plexní.
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:
22
>> % 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
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.
23
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 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
24
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
>> % 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
1.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
25
>> 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
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říka-
zem
sage: line([(1,2), (2,4), (3,0), (5,1), (4,5)])
26
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ů
>> x=linspace(0,2*pi);y=sin(x);
>> plot(x,y,’--’)
27
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
28
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
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))
29
1 2 3 4 5 6
-2
-1.5
-1
-0.5
0.5
1
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))
30
-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 x-ové a y-ové souřadnice.
>> phi=pi/6;
>> R=[cos(phi), -sin(phi);sin(phi), cos(phi)];
>> xyR=R*[x;y];
31
>> 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 explicitně.
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))
32
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ů.
33
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.
34
2.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.
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)
35
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
>> 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)
36
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 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ř.
37
>> 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. Zadá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
8 8
9 2
5 8
1 2
1 9
38
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í, 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:
39
\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
40
0 0 0 0 3
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.
41
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.
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é matice
42
>> 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
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
43
Ve skutečnosti je mezi výsledky rozdíl řádově 10−14
.
2.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ů.
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:
44
sage: type(v)
<type ’sage.modules.vector_integer_dense.
Vector_integer_dense’>
sage: type(u)
<type ’list’>
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]
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]
45
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:
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:
46
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ůso-
bem:
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, [
(1, -1)
], 1)]
47
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().
2.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 =
48
2 0 0
0 4 0
0 0 4
>> 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 =
49
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
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í.
50
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]
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
51
[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:
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í.
52
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]
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
53
-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)
<type ’sage.rings.integer.Integer’>
sage: type(s2)
<type ’sage.modules.vector_integer_dense.
Vector_integer_dense’>
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’
54
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];
>> B=null(A)
B =
55
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ě:
56
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],
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 =
57
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
-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
58
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é.
59
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.
3.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 popu-
lace.
>> tau=0.1;
>> M=[0 3 3; 0.2 0 0;0 0.6 tau];
>> p=[100;0;0];
>> p=M*p
60
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í 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
61
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
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:
62
>> 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
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
63
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.
3.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ů.
64
>> 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á 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é
65
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.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 doplně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.
66
To ovšem odpovídá násobení 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














67
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
částečná pivotáž. Existuje i úplná pivotáž, kdy vedoucí prvek nehledáme v
68
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
69
0 0 0 -1.6026
P =
0 0 1 0
0 1 0 0
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,
70
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.
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.
71
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 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
72
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’)
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’)
73
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ě.
74
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
75
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)
−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
76
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)
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))
77
nebo parametric plot3d pro spírálu
sage: parametric_plot3d((sin(u),cos(u),u/10),(u,0,20))
nebo zakroucenou pneumatiku:
78
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")
79
Literatura
[1] Anděl, J.: Matematická Statistika. Řada teoretické literatury, SNTL,
1978.
URL http://books.google.cz/books?id=AP mGgAACAAJ
[2] PASEKA, J.; ZLATOŠ, P.: Lineární algebra a geometrie I [online]. Elportál,
Masarykova univerzita, první vydání, 2010 [cit. 2013-10-16].
URL Dostupné z WWW <http://is.muni.cz/elportal/?id=902268>
80