3 Soustavy lineárních rovnic, tvorba grafů

3 Soustavy lineárních rovnic, tvorba grafů Studijní cíl Třetí blok je věnován využití maticových operací při řešení soustav lineárních rovnic a dále tvorbě grafů (2D , 3D a speciálních). Též se zabývá možnostmi editace grafů a jejich exportem.

Doba nutná k nastudování

2 ‐ 3 hodiny

Průvodce studiem Předpokládá se, že čtenář je seznámen s prostředím MATLABu a jeho základ‐ ními příkazy, se základy maticové počtu a komplexními čísly. Dále se předpo‐ kládá, že čtenář zná principy a základní pojmy z objektového programování. Při studiu je vhodné mít spuštěný MATLAB a jednotlivé příkazy či funkce hned vyzkoušet a případně získat podrobnější informace využíváním nápovědy k programu. Na závěr jsou uvedeny řešené příklady na procvičení. V textu jsou použity následující typografické konvence: Calibri 11, Bold, Italic nové pojmy či informace k zapamatování Calibri 11 označení klávesy Courier New 10, Bold názvy nástrojů MATLABu Courier New 10, Bold upřesnění nápovědy (help téma) Courier New 10 názvy příkazů, funkcí a objektů Courier New 10 označení části příkazu, která se může vynechat Courier New 10, Italic názvy proměnných použitých v programu Courier New 9 příkazy příkazové řádky

3.1 Soustavy lineárních rovnic Jednou ze základních úloh lineární algebry je úloha nalézt řešení soustavy line‐ árních rovnic. Tuto úlohu můžeme zapsat buď jako samostatné rovnice

a11 x1 + a12 x2 + K + a1s xs = b1 a21 x1 + a22 x2 + K + a2 s xs = b2 M ar1 x1 + ar 2 x2 + K + ars xs = br

r

nebo

∑a x i =1

ij

j

= bj

pro j = 1K s

nebo v ekvivalentním maticovém tvaru

Ax = b

⎡ a11 ⎢a A = ⎢ 12 ⎢ M ⎢ ⎣ a r1

a12 K a1s ⎤ a22 K a2 s ⎥⎥ M O M ⎥ ⎥ ar 2 K ars ⎦

KŘP/IMSW Modelování ve výpočtových software

⎡ x1 ⎤ ⎡ b1 ⎤ ⎢x ⎥ ⎢b ⎥ 2⎥ ⎢ x= b = ⎢ 2⎥ ⎢M⎥ ⎢M⎥ ⎢ ⎥ ⎢ ⎥ ⎣ xs ⎦ ⎣bs ⎦ 3‐1 (24) 17.8.11

František Dušek KŘP FEI Univerzita Pardubice

V následujícím textu budeme používat pouze zápis v maticovém tvaru a bude‐ me se zabývat možnostmi MATLABu pro řešení soustavy lineárních rovnic.

Soustava lineárních rovnic má buď jedno, nekonečně mnoho nebo žádné řeše‐ ní. Aby řešení existovalo, musí vždy platit, že hodnost 1 matice (funkce rank) soustavy A musí být stejná jako hodnost tzv. rozšířené matice soustavy A|b tj. h(A)= h(A|b). O kterou variantu jde lze obecně určit podle hodnosti matice soustavy A. Pokud je hodnost h(A) rovna počtu neznámých, pak má soustava právě jedno řešení, pokud je hodnost h(A) menší pak má soustava nekonečně mnoho řešení a pokud je hodnost h(A) větší pak řešení neexistuje. Podle počtu rovnic r a neznámých s je možné ihned určit že: a) soustava má nekonečně mnoho řešení – nedourčená soustava r<s počet rovnic r je menší jako počet neznámých s tj. h(A) ≤ r < s b) soustava pravděpodobně nemá žádné řešení – přeurčená soustava r>s počet rovnic r je větší jako počet neznámých s tj. h(A)≤s ≠ h(A|b) ≤ s+1 V případě že počet rovnic r je stejný jako počet neznámých s (tj. r=s=n soustava n rovnic o n neznámých) pak nelze jednoduše odhadnout o kterou variantu jde h(A) = h(A|b) ≤ n.

3.1.1 Soustava n rovnic o n neznámých Základní řešení (s využitím inverzní matice) je možné zapsat jako maticovou rovnici pro nalezení všech neznámých ve formě sloupcového vektoru

x = A −1b

V případě, že soustava rovnic nemá právě jedno řešení, inverze čtvercové ma‐ tice soustavy A neexistuje (matice je singulární). >> A=[0,2,3; 1,0,3; 1,2,0]; >> b=[1;1;1]; >> rank(A) ans = 3 >> rank([A,b]) ans = 3 >> x=inv(A)*b >> x=A^(-1)*b >> x=A^-1*b x = 0.5000 0.2500 0.1667 >> A*x ans = 1 1 1

% čtvercová matice soustavy % sloupcový vektor pravé strany % hodnost matice soustavy % hodnost rozšířené matice soustavy % funkce inverze % mocnina matice % i tento zápis je správně

% kontrola řešení

1

Hodnost matice je vždy rovna nebo menší než je menší z rozměrů matice. Vyjadřuje maximální počet vzájemně nezávislých řádků či sloupců.


3‐2 (24) 17.8.11


Pro přímý výpočet jednotlivých neznámých je také možné využít i Crammerovo pravidlo 2 obecně popsané následující rovnicí (podíl determinantů 3 matic).

⎛ ⎡ a11 K a1i−1 b1 ⎜ xi = det ⎜ ⎢⎢ M O M M ⎜ ⎢a K a bn ni −1 ⎝ ⎣ n1

a1i+1 K a1n ⎤ ⎞ ⎟ M O M ⎥⎥ ⎟ det( A) ani +1 K ann ⎥⎦ ⎟⎠

V případě, že řešení neexistuje je determinat matice soustavy roven nule. >> >> >> >> x1 >> >> x2 >> >> x3

A=[0,2,3;1,0,3;1,2,0]; b=[1;1;1]; A1=[b,A(:,2:3)]; x1=det(A1)/det(A) = 0.5000 A2=[A(:,1),b,A(:,3)]; x2=det(A2)/det(A) = 0.2500 A3=[A(:,1:2),b]; x3=det(A3)/det(A) = 0.1667

% čtvercová matice soustavy % sloupcový vektor pravé strany % doplněná matice soustavy i=1 % doplněná matice soustavy i=2 % doplněná matice soustavy i=3

Obě uvedená řešení využívají standardní maticové operace a funkce. V MATLABu existuje speciální operátor \ (levé maticové dělení, též funkce mldivide), který určí řešení soustavy lineárních rovnic stejně jako postup s in‐ verzí – tj. platí x=inv(A)*b=A\b. Výhoda použití tohoto postupu je v numericky efektivnějším způsobu výpočtu. Místo numericky náročného výpočtu inverzní matice a následného maticového násobení se úloha v případě, že je matice soustavy A trojúhelníková, řeší Gaussovou eliminací 4 .

⎡ a11 0 ⎢a ⎢ 21 a22 ⎢ a31 a32 ⎢ M ⎣ M

0 0 a33 M

L⎤ ⎡ x1 ⎤ ⎡ b1 ⎤ L⎥⎥ ⎢⎢ x2 ⎥⎥ ⎢⎢b2 ⎥⎥ = L⎥ ⎢ x3 ⎥ ⎢b3 ⎥ ⎥⎢ ⎥ ⎢ ⎥ O⎦ ⎣ M ⎦ ⎣ M ⎦

x1 = b1 a11 x2 = (b2 − a21 x1 ) a22 x3 = (b3 − a32 x2 − a31 x1 ) a33 M

V případě obecné matice soustavy se tato rozloží na součin dolní a horní trojú‐ helníkové matice ([L,U]=LU(M) LU dekompozice M=L*U) a řeší se dvojná‐ sobná Gaussova eliminace.

A.x = b

{A = L.U

x = A\b

y = U \ b x = L \ y}

Následující ukázka demonstruje efektivitu výpočtu pomocí operátoru levého maticového dělení oproti použití inverzní matice na výpočtu soustavy 2000 rovnic o 2000 neznámých. Je dosaženo úspory cca 33% původního času. >> A=rand(2000); >> rank(A) ans = 2000 >> b=rand(2000,1);

% matice soustavy 2000x2000 % kontrola existence řešení % vektor pravé strany 2000x1

2

Neznámá xi je podíl determinantu doplněné matice soustavy, ve které je i‐tý sloupec nahrazen vektorem pravé strany a determinantu původní matice soustavy.

3

Funkce det MATLABu

4

Postupný výpočet neznámých s využitím již spočítaných hodnot.


3‐3 (24) 17.8.11


>> t1 t1 t1 >> t2 t2 t2

tic, = = = tic, = = =

x1=inv(A)*b; t1=toc 1.2655 1.2909 1.2203 x2=A\b; t2=toc 0.7157 0.8218 0.6723

% % % % % % % %

výpočet 1 s měřením času doba výpočtu 1.běh v sec doba výpočtu 2.běh v sec doba výpočtu 3.běh v sec výpočet 2 s měřením času doba výpočtu 1.běh v sec doba výpočtu 2.běh v sec doba výpočtu 3.běh v sec

Pokud není systém Ax=b standardním způsobem řešitelný, pak je oznámena chyba. I v tomto případě lze nalézt přibližné řešení pomocí funkce pinv (Moo‐ re‐Penroseova pseudoinverze). Je nalezeno takové řešení, které minimalizuje kvadratickou chybu řešení e=Ax‐b přesněji kvadratickou normu chyby e tj. min Ax − b 2 (funkce norm) např. x

>> A=magic(5); % matice soustavy původní 5x5 >> rank(A) % hodnost původní matice ans = 5 >> det(A) % hodnota determinantu ans = 5.0700e+006 >> A(:,2)=zeros(5,1); % upravená matice 5x5 >> rank(A) % hodnost upravené matice ans = 4 % je menší než rozměr matice >> det(A) % hodnota determinantu ans = 0 >> b=[1;2;3;4;5]; % vektor pravé strany % standardní řešení % přibližné řešení >> x1=A\b >> x2=pinv(A)*b Warning: Matrix is singular to working precision. x1 = NaN x2 = 0.0273 Inf -0.0000 0.1885 0.1845 0.0114 0.0169 0.0120 0 >> A*x1 >> A*x2 ans = NaN ans = 0.7837 NaN 2.1555 NaN 2.8445 NaN 4.1320 NaN 4.9457

3.1.2 Přeurčená soustava rovnic (žádné řešení) V případě, že máme více rovnic než je neznámých, neexistují hodnoty nezná‐ mých takové, aby všechny rovnice platily. Zdálo by se, že nemá význam se ta‐ kovými soustavami zabývat. Opak je pravdou, mnoho praktických úloh 5 vede na takovouto soustavu rovnic a hledá se „nejlépe vyhovující“ řešení. Nejčastěji se za „nejlépe vyhovující řešení“ volí takové řešení, které minimalizuje kvadra‐ tickou chybu řešení e=Ax‐b. Jde o stejnou formulaci jako v případě soustavy se stejným počtem rovnic a neznámých a singulární maticí soustavy. V MATLABu lze řešení přeurčené soustavy lineárních rovnic spočítat třemi způsoby >> A=[1 2;-2 1;1 -2];

% matice soustavy 3x2

5

Např. proložení naměřených bodů přímkou „metodu nejmenších čtverců“ je nejjednodušší úlohou tzv. lineární regrese – úlohy, která vede na přeurčený systém lineárních rovnic


3‐4 (24) 17.8.11


>> b=[1;1;1]; % vektor pravé strana 3x1 >> x1=A\b >> x2=pinv(A)*b >> x3=inv(A'*A)*A'*b x1 = 0.0400 x2 = 0.0400 x3 = 0.0400 0.1200 0.1200 0.1200 % ověření řešení % chyba % kvadratická norma >> A*x1 >> A*x1-b >> norm(A*x1-b,2) ans = 0.2800 ans = -0.7200 ans = 1.6971 0.0400 -0.9600 -0.2000 -1.2000

3.1.3 Nedourčená soustava rovnic (nekonečně mnoho řešení) V případě, že máme méně rovnic než je neznámých lze nalézt nekonečně mno‐ ho kombinací neznámých, které těmto rovnicím vyhovují. Mezi těmito kombi‐ nacemi je možné vybrat ty, které mají určité speciální vlastnosti. Řešení pomocí levého maticového dělení (asi) vybírá kombinaci s maximálním počtem nulo‐ vých hodnot ve vektoru neznámých, řešení s pseudoinverzí hledá řešení s minimální kvadratickou normou vektoru neznámých stejně jako uvedený maticový výraz >> A=[-1 2 2 1;1 -2 2 1; 1 2 -2 1]; % matice soustavy 3x4 >> b=[1;1;1]; % vektor pravé strany 3x1 >> x1=A\b >> x2=pinv(A)*b >> x3=A'*inv(A*A')*b x1 = 1.0000 x2 = 0.4000 x3 = 0.4000 0.5000 0.2000 0.2000 0.5000 0.2000 0.2000 0 0.6000 0.6000 >> A*x1 >> A*x2 ans = 1.0000 ans = 1.0000 1.0000 1.0000 1.0000 1.0000 >> norm(x1) >> norm(x2) ans = 1.2247 ans = 0.7746

3.2 Základy práce s grafikou V případě, že výsledkem je velké množství čísel, je užitečné grafické zobrazení výsledků – vizualizace výsledků ve formě grafu. MATLAB poskytuje nejen velké množství funkcí pro tvorbu různých typů 2D i 3D grafů ale i prostředky pro ma‐ nipulaci s grafickými objekty, editaci hotových grafů a jejich export do různých formátů. Za zmínku stojí, že původně byl grafický výstup podporován pouze jako tisk na PostScriptové tiskárně a podpora pro grafický výstup na obrazovku byla doplněna až později. Práce s grafikou je plně objektová, ač vznikla v době, kde objektově orientované programování nebylo moc rozšířené. Systém práce s grafikou v MATLABu se nazývá Handle Graphics. Každý grafický objekt je ur‐ čen svém ukazatelem (reference, handler), který je představován reálným čís‐ lem. Grafické funkce vrací handler vytvářeného objektu. Přehled základních funkcí viz help graphics. Všechny grafické objekty jsou odvozeny od objektu Root a jsou umístěny v jednom či více objektech Figure (okno Windows). Vytvoření prázdného nového okna (včetně uložení handleru okna do proměnné hf, může být vyne‐ cháno) je příkazem hf=figure.


3‐5 (24) 17.8.11


Může být vytvořeno více oken ale pouze jedno je „current figure“ – tj. aktivní okno nebo také okno, které má fokus (naposledy vytvořené nebo to, na které se naposledy kliklo myší – okno se zvýrazněným záhlavím). Handler aktivního okna lze získat příkazem hf=gcf. Aktivním se učiní okno příkazem figure(hf). Aktivní okno se zruší příkazem clf nebo close, případně okno s handlerem hf příkazem clf(hf) nebo 6 close(hf). Příkaz close all uzavře všechna okna.

Chceme‐li nalézt nějaký konkrétní objekt (respektive jeho handler) např. v exis‐ tujícím objektu (např. v okně) lze použít příkaz h=findobj(ObjectHandler, PropertyName, Value)

tento příkaz prohledá objekt určený handlerem ObjectHandler a všechny jeho potomky a vrátí handlery (vektor handlerů) všech objektů, jejichž vlast‐ nost PropertyName má hodnotu Value. Blíže viz help findobj. Základem grafů je objekt Axes (obdélníková oblast), do kterého příkazy pro kreslení grafů (např. příkaz plot) umisťují další grafy. Vytvoření nové‐ ho objektu Axes (včetně uložení handleru objektu Axes do proměnné ha, které může být vynecháno) v aktivním okně je příkazem ha=axes(PropertyName, Value)

Pokud nenastavíme vlastnost Position, pak oblast objektu Axes zabírá celou plochu okna 7 . Hodnota vlastnosti NextPlot ovlivňuje, zda následující požada‐ vek na přidání grafu do objektu existující graf přepíše (implicitní hodnota replace) nebo nový graf přidá (hodnota add). Nastavení těchto hodnot je možné explicitně pomocí obecného příkazu pro nastavení vlastností objektů set(handle, PropertyName, Value) nebo speciálním příkazem hold on/off. Tento příkaz mění vlastnost pouze aktivního objektu Axes v aktivním okně 8 tj. následující příkazy jsou ekvivalentní hold on = set(gca,'NextPlot','add') režim přidávání hold off = set(gca,'NextPlot','replace') režim přepisu

Při tvorbě grafů je několik možností kam a jak se graf vytvoří a) žádný objekt Figure (okno) neexistuje – funkce tvořící graf vytvoří objekt Figure (nové okno) včetně objektu Axes a vykreslí do něj graf b) existuje jeden či více objektů Figure – funkce tvořící graf vytvoří objekt Axes v aktivním okně a vykreslí graf c) existuje jeden či více objektů Figure každý obsahující objekt Axes – funkce tvořící graf vykreslí graf do aktivního objektu Axes; je‐li první para‐ metr funkce handler objektu Axes, pak funkce vytvoří graf v zadaném objek‐ tu (podle aktuální hodnoty vlastnosti NextPlot objektu Axes).

6

Příkaz (metoda) close může být přetížený (overloaded) podle typu objektu, na který ukazuje použitý handler

7

Je možné vložit více objektů Axes do jednoho okna

8

Handler objektu Axes v aktivním okně lze získat příkazem ha=gca.


3‐6 (24) 17.8.11


V následujícím příkladě se vytvoří postupně tři okna. Ve druhém okně je vytvo‐ řen prázdný objekt Axes. Příkaz plot vytvoří v aktivním okně č. 3 (posledně vytvořené) graf funkce sinus. Poté je jako aktivní nastaveno okno č. 1 a v něm vytvořen graf funkce cosinus. Následuje vykreslení funkce sinus do objektu Axes v okně č. 2, změna vlastnosti NextPlot tohoto objektu a přidání grafu funkce cosinus. >> >> >> >> >> >> >> >> >> >> >> >> >> >> >>

x=0:0.1:2*pi; ys=sin(x); yc=cos(x); hf1=figure; % okno č.1 hf1 = 1 hf2=figure; % okno č.2 hf2 = 2 ha2=axes; % axes v akt.okně č.2 ha2 = 346.0079 hf3=figure; % okno č.3 hf3 = 3 plot(x,ys) % přepiš obsah aktuálního okna č.3 figure(hf1) % nastav jako aktuální okno č.1 plot(x,yc) % přepiš obsah aktuálního okna č.1 plot(ha2,x,ys) % přepiš obsah okna č.2 set(ha2,'NextPlot','add') % režim přidávání axes č.2 plot(ha2,x,yc) % přidej do okna č.2

3.2.1 2D grafy (dvourozměrné grafy) Nejčastěji používanými jsou 2D grafy tj. grafy v rovině x‐y. Přehled příkazů vztahující se k této problematice viz help graph2d. Základním a nejčastěji používaným příkazem pro tvorbu 2D grafu je příkaz plot s obecnou syntaxí hl=plot(ha,x,y,linespec, …)

kde první parametr ha (může se vynechat) je handler objektu Axes, do které‐ ho se bude kreslit graf, x, y jsou stejně orientované a stejně dlouhé vektory obsahující x‐ové a y‐ové souřadnice bodů, které se mají vykreslit, linespec je maximálně tříznakový řetězec. Tento řetězec (může se vynechat) obsahuje jeden až tři znaky – znak určující barvu čáry, znak určující typ spojovací čáry mezi body a znak určující značku pro vykreslení bodu (libovolný znak se může vynechat, na pořadí nezáleží). V případě, že chceme jedním příkazem vykreslit více průběhů do jednoho grafu, je možné trojice (či dvojice) parametrů opako‐ vat. Návratová hodnota (může se vynechat) je handler na grafický objekt typu Line nebo vektor handlerů v případě vykreslení více průběhů. Blíže viz help plot. Kromě tohoto typu 2D grafu je možné použít příkazy pro vykreslení gra‐ fů s logaritmickým měřítkem na ose x (semilogx), ose y (semilogy) či na obou (loglog). Další možností je použít příkaz hl=polar(ha,theta,rho,linespec)

pro kreslení v polárních souřadnicích (theta je vektor úhlů radiánech a rho vektor poloměrů jednotlivých bodů). Návratová hodnota (může se vynechat) je handler na grafický objekt typu Line. Do jednoho objektu Figure (okna) je možné vložit více objektů Axes s vyme‐ zením obdélníkové oblasti uvnitř okna pro každý objekt. Toto vložení je možné provést dvěma způsoby. Jednoduší způsob využívající rozdělení aktivního okna na matici obdélníkových oblastí o r (≤9) řádcích a s (≤9) sloupcích je příkazem KŘP/IMSW Modelování ve výpočtových software

3‐7 (24) 17.8.11


ha=subplot(r,s,n) nebo

ha=subplot(rsn)

kde n je pořadové číslo oblasti počítáno po řádcích. Příkaz vrací handler vytvo‐ řeného objektu Axes. Druhý způsob s explicitním určením polohy a velikosti objektu Axes uvnitř aktivního okna je příkazem ha=axes('Position',[xL,yB,xW,yH])

kde hodnota parametru je vektor 4 hodnot (hranaté závorky jsou součástí de‐ finice vektoru), xL, yB jsou souřadnice levého dolního rohu obdélníku, do kterého se bude graf kreslit a xW, yH je jeho šířka a výška. Všechna čísla jsou v rozmezí 0‐1 a vyjadřují relativní vzdálenost v rámci velikosti okna. Grafům 9 v objektu Axes s handlerem ha (pokud je vynechán použije se aktivní objekt Axes) lze změnit implicitní automatické měřítko na vlastní hodnoty příkazem

axis(ha,[xmin,xmax,ymin,ymax])

Další možnosti použití příkazu axis viz help axis. Vytvořený graf lze doplnit o mřížku – příkaz grid, nadpis – příkaz title(text), popis osy x– xlabel(text), popis osy y – ylabel(text), legendu – legend(text1, text2, …) atd. Blíže viz help graph2d. Většina výše uvedených příkazů je použita v následujícím příkladu, jehož řešení je uvedeno ve dvou variantách. Každá varianta vytvoří stejné okno. V každém okně tři objekty Axes obsahující různé grafy. V první variantě jsou použity pře‐ devším jednodušší „high level“ funkce, v druhé variantě jsou více využity expli‐ citní nastavení vlastností objektů >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >>

% ------------- první varianta ---------------------------% data spirála polární souřadnice th=0:0.01:4*pi; % úhel ro=th; % poloměr % data spirála kartézské souřadnice xs=th.*cos(th); ys=th.*sin(th); % x, y subplot(2,2,1); % první Axes % první graf - spirála kartézské plot(xs,ys,'b','LineWidth',2) title('Spirála x,y',... 'FontSize',14) % nadpis grid % mřížka legend('kartézské s.',... 'Location','SouthWest') xlabel('osa x'), ylabel('osa y') subplot(2,2,2); % druhý Axes % první graf - spirála polární hp=polar(th,ro,'r'); set(hp,'LineWidth',2) title('Spirála \theta, \rho',... 'FontSize',14); grid % mřížka subplot(2,1,2); % třetí Axes hold on % režim přidávání % první graf - spirála kartézské plot(xs,ys,'b','LineWidth',2) % druhý graf - vývoj souřadnic plot(th,xs,'r-.',th,ys,'g--',...

9

Všechny grafy v objektu Axes sdílejí stejné souřadnice, tj. mají stejné měřítko


3‐8 (24) 17.8.11


'LineWidth',2) >> title('Tři průběhy',... 'FontSize',14) % nadpis >> grid % mřížka >> legend('xs,ys','\theta,xs','\theta,ys') Spirála x,y

10

Spirála θ, ρ 90 20

120

osa y

5

60

10

150

30

0 180

0

-5 210

-10

240

kartézské s. -15 -10

-5

0

330

5

10

270

300

15

osa x

Tři průběhy

15

xs,ys

10

θ,xs θ,ys

5 0 -5 -10 -15 -10

>> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >>

-5

0

5

10

15

% ------------- druhá varianta --------------------% data spirála polární souřadnice th=0:0.01:4*pi; % úhel ro=th; % poloměr % data spirála kartézské souřadnice xs=th.*cos(th); ys=th.*sin(th); % x,y hf=figure; % další okno ha(1)=axes('Position',[0.1,0.55,0.35,0.35]); %subplot(2,2,1); % první Axes ha(2)=axes('Position',[0.55,0.55,0.35,0.35]); %subplot(2,2,2); % druhý Axes ha(3)=axes('Position',[0.1,0.1,0.8,0.35]); %subplot(2,1,2); % třetí Axes % *** grafy % první Axes - spirála kartézské hl(1)=plot(ha(1),xs,ys,'b'); % druhý Axes - spirála polární hl(2)=polar(ha(2),th,ro,'r'); % třetí Axes první graf - spirála kartézské hl(3)=plot(ha(3),xs,ys,'b'); set(ha(3),'NextPlot','add') % hold on % režim přidávání % třetí Axes druhý graf - vývoj souřadnic hl(4:5)=plot(ha(3),th,xs,'r-.',th,ys,'g--'); % hromadné nastavení mřížky set(ha,'XGrid','on'), set(ha,'YGrid','on') % hromadné nastavení tlouštky čar set(hl,'LineWidth',2) % *** nadpisy ht(1)=title(ha(1),'Spirála x,y'); ht(2)=title(ha(2),'Spirála \theta, \rho'); ht(3)=title(ha(3),'Tři průběhy'); % hromadné nastavení velikosti textu


3‐9 (24) 17.8.11


>> >> >> >> >> >> >> >>

set(ht,'FontSize',14) % *** legendy le(1)=legend(ha(1),'kartézské s.'); set(le(1),'Location','SouthWest') le(2)=legend(ha(3),'xs,ys','\theta,xs','\theta,ys'); % *** popisy os hx(1)=xlabel(ha(1),'osa x'); hy(1)=ylabel(ha(1),'osa y');

3.2.2 3D grafy (třírozměrné grafy) Jako 3D grafy se označují grafy zobrazující nějaký útvar v prostoru. Přehled příkazů MATLABu pro práci s 3d grafy viz help graph3d. Podle formy dat potřebných pro zobrazení útvaru jsou v podstatě dva druhy útvarů. První – prostorová křivka („provázek v prostoru“) – je reprezentován uspořádanými trojicemi bodů {x,y,z} (v MATLABu trojice vektorů x, y, z) a příkaz pro vy‐ kreslení h=plot3(x,y,z)

Stejně jako v případě příkazu plot, tak se graf vytvoří v aktuálním objektu Axes. Pokud objekt Figure s objektem Axes neexistuje tak se vytvoří. Příklad použití je uveden na konci kapitoly. Složitější 10 je situace v případě, že útvar, který chceme vykreslit, je plocha v prostoru. V dalším se budeme zabývat pouze jednodušší variantou, kdy pro každý bod roviny x‐y existuje pouze jeden bod plochy. V tomto případě může být vykreslovaná plocha pro potřeby vykreslení popsána maticí Z výšek plochy nad rovinou x‐y v obdélníkové síti bodů a polohou těchto bodů (vektory x, y). Situace je znázorněna na obrázku. Plochu pak můžeme vykreslit jako drátový model (podobně jako na obrázku) příkazem h=mesh(x,y,Z)

x

nebo s vyplněnými plochami příkazem h=surf(x,y,Z)

Na rozdíl od 2D grafů je poněkud složitější příprava matice Z pro všechny kombinace hodnot z vektorů x a y. Máme‐li k dispozici funkční předpis z=f(x, y), můžeme využít funkci MATLABu

z x1

x2

y1

x3

y2

[X,Y]=meshgrid(x,y)

x4

y3

x5

y4

x6

y5

x7

y6

ina rov -y x

která nám z původních vektorů x, y vytvoří matice X, Y takové, že po dosazení (v MATLABu) do původní funkce vytvoří matici Z tj. Z=f(X,Y). Pozor na to, že ač se do funkce dosadí matice nejde o maticové operace tj. je potřeba použít operace prvek po prvku.

y

10

Složitější z pohledu přípravy dat popisujících danou plochu.


3‐10 (24) 17.8.11


V následujícím příkladě ukážeme vykreslení prostorové křivky popsané parametricky zadanou funkcí

y (t ) = t sin(2t ) z (t ) = t 2 pro 0 ≤ t ≤ 4π

x(t ) = t cos(t )

a plochy popsané rovnicí

z ( x, y ) = cos( xy) 2π 2 − x 2 − y 2 >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >>

pro − π ≤ x ≤ π a − π2 ≤ y ≤

π 2

t=0:0.1:4*pi; x=t.*cos(t); y=t.*sin(2*t); z=t.^2; subplot(2,1,1) plot3(x,y,z,'LineWidth',2) title('plot3','FontSize',14) grid x=-pi:0.3:pi; y=-pi/2:0.3:pi/2; [X,Y]=meshgrid(x,y); Z=cos(X.*Y).*(X.^2+Y.^2); subplot(2,2,3) mesh(x,y,Z) % drátový model title('mesh','FontSize',14) xlabel('osa x'), ylabel('osa y'), zlabel('osa z') subplot(2,2,4) surf(x,y,Z) % vyplněné plochy title('surf','FontSize',14) xlabel('osa x'), ylabel('osa y'), zlabel('osa z')

plot3 200

100

0 15

10

5

0

-5

-10

-15

-5

-10

10

10

0

0

-10 -20 2

osa y

5

-10 -20 2

0 -2

-5

15

surf

osa z

osa z

mesh

0

10

5

0

osa y

osa x

5

0

0 -2

-5

osa x


3‐11 (24) 17.8.11


3.2.3 Speciální grafy Kromě základních 2D grafů (plot, polar) a 3D grafů (plot3, mesh, surf) MATLAB obsahuje řadu dalších, specializovaných grafů viz help specgraph. Některé z nich jsou dále uvedeny včetně příkladu použití. h=stairs(x,y) schodovitý průběh křivky

1

Příklad Vykreslete průběh funkce y=sin(x) na intervalu 0 ≤ x ≤ 2π v 120 bodech a pomocí schodové funkce v 12 bodech

0.6 0.4 0.2

Řešení >> x=linspace(0,2*pi,120); >> >> >> >> >> >> >> >>

plot stairs

0.8

x12=linspace(0,2*pi,12); h=plot(x,sin(x)); set(h,'LineWidth',2) hold on h=stairs(x12,sin(x12),'ro-'); set(h,'LineWidth',2) grid legend('plot','stairs')

0

-0.2

-0.4

-0.6

-0.8

-1

0

1

2

3

4

5

6

7

1

h=area(x,y) graf s vyplněnou plochou pod

křivkou

0.6

Příklad Vykreslete průběh funkce y=sin(x) na intervalu 0≤x≤2π s pomocí funkce area v 12 bodech

0.4

0.2

0

-0.2

Řešení >> x12=linspace(0,2*pi,12); >> >> >> >> >> >>

area 0.8

h=area(x12,sin(x12)); set(h,'FaceColor',[.6,.9,1.0]) grid set(h,'LineWidth',2) h=legend('area'); set(h,'FontSize',14)

-0.4 -0.6

-0.8 -1

0

1

2

3

4

5

6

h=bar(x,y) graf s vertikálními sloupci h=barh(x,y) graf s horizontálními sloupci

Příklad Vykreslete průběh funkce y=sin(x) na intervalu 0≤x≤2π s pomocí funkce area v 24 bodech

1 bar 0.5

0

Řešení >> x24=linspace(0,2*pi,24); >> >> >> >> >> >> >> >> >> >>

subplot(2,1,1) h=bar(x24,sin(x24)); set(h,'FaceColor',[.1,.2,.3]) legend('bar') grid subplot(2,1,2) h=barh(x24,sin(x24)); set(h,'FaceColor',[.3,.2,.1]) legend('barh') grid


-0.5

-1 -1

0

1

2

3

4

5

6

7

8 barh 6 4 2 0 -2 -1

3‐12 (24) 17.8.11

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1


h=pie(x,v,p)

koláčový graf 2D

h=pie3(x,v,p)

koláčový graf 3D

a1

a2

a7 a3

A4

Příklad Vykreslete koláčový graf 2D a 3D podle zadaných dat s barevnou paletou „summer“

a6 a5

Řešení >> x=[1,2,3,4,5,6,7]; >> v=[0,0,0,1,0,0,0]; >> p={'a1','a2','a3','A4',... 'a5','a6','a7'}; >> subplot(2,1,1) >> h=pie(x,v,p); >> subplot(2,1,2) >> h=pie3(x,v,p); >> colormap summer

a7 a1 a2 a6

a3

A4

a5

h=errorbar(x,y,e,ls)

graf se sloupci velikosti ±ei u každého bodu {xi, yi}, volitelný řetězec ls určuje barvu, značku a typ čary stejně jako u příkazu plot 5.5

5

Příklad Vykreslete graf podle zadaných dat s chybovými úsečkami, černá čára síla čary 3, označení bodů kolečkem velikosti 12 bodů

4.5

4

3.5

Řešení >> y=[5,4,3,2,3,4,5]; >> >> >> >> >> >>

x=[1,2,3,4,5,6,7]; e=1./y; h=errorbar(x,y,e,'ko-'); grid set(h,'MarkerSize',12) set(h,'LineWidth',3)

3

2.5

2

1.5

0

1

2

3

4

5

6

7

8

N=hist(x,M)

histogram, graf četnosti výskytu čísel ve vektoru x v rozsazích (xmax‐xmin)/M, funkce vrací počet čísel v každém rozsahu. Je‐li vý‐ stupní parametr vynechán vykreslí se histogram

Příklad Nakreslete histogram 10000 čísel s rovnoměrným a normálním rozdě‐ lením N(0,1) ve 12 oblastech Řešení >> x1=rand(1,10000); >> >> >> N1

>> >> >> N2

>>

x2=randn(1,10000); subplot(2,1,1); N1=hist(x1,12) = 835 840 829 826 807 821 864 894 805 839 835 805 hist(x1,12) subplot(2,1,2); N2=hist(x2,12) = 3 32 222 765 1768 2523 2440 1522 560 140 23 2 hist(x2,12)

1000 800 600 400 200 0

0

0.1

0 -4

-3

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

4

5

3000 2500 2000 1500 1000 500 -2

-1

0

1

2

3


3‐13 (24) 17.8.11


[C,h]=contour(x,y,Z)

řezy plochy v prostoru rovnoběžné s rovinou x‐

y, vykresleny jsou obrysy [C,h]=contourf(x,y,Z)

řezy plochy v prostoru rovnoběžné s rovinou x‐y, vykresleny jsou plochy

z ( x, z ) = xe

Vykreslete graf funkce

0 -0.2

-0.5 3

2

− x2 − y 2

a její řezy jako obrysy a jako plochy

1

0

-1

-2

-1

-2

2

1

0

-0.4

3 2

Řešení >> x=-2:0.2:2;

1

-0 .1

0

-0.2 -0.3 -0 .4

-0 .3

-0.2

-0 .1

-1 -2 -2

-1.5

-1

0.4

-0.5

0.1

0.2 0.3

0.3

0.2

0

y=-2:0.2:3; [X,Y]=meshgrid(x,y); Z=X.*exp(-X.^2-Y.^2); subplot(3,1,1) mesh(x,y,Z) colorbar subplot(3,1,2) [C,h]=contour(x,y,Z); set(h,'ShowText','on') grid subplot(3,1,3) [C,h]=contourf(x,y,Z); colorbar

-0 .1

>> >> >> >> >> >> >> >> >> >> >> >> >>

0.2

0

0.1

Příklad

0.4 0.5

0

0.5

0.1

1

1.5

2

3

0.4

2

0.2

1

0

0

-0.2

-1 -2 -2

-1.5

-1

-0.5

0

0.5

1

1.5

2

-0.4

h=fill(x,y,c,…)

jeden nebo více vyplněných 2D polygonů, parametr c určuje barvu výplně

h=fill3(x,y,z,c,…)

jeden nebo více vyplněných 3D polygonů, pa‐

rametr c určuje barvu výplně h=patch(x,y,z,c,…)

přidá jeden nebo více vyplněných 2D nebo 3D polygonů do aktivního objektu Axes, parametr c určuje barvu výplně

Příklad Nakreslete červený a žlutý trojúhelník v rovině (jeden graf) a v prostoru (druhý graf). Potom přidejte třetí zelený trojúhelník do obou grafů

4

3

2

1

Řešení >> x1=[1,2,1]; y1=[1,2,3]; >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >>

z1=[1,2,0]; x2=[3,4,4]; y2=[2,1,3]; z2=[2,0,1]; x3=[1.5,2.5,3.5]; y3=[4.5,2.5,4.5]; z3=[1.5,-1.5,1.5]; subplot(2,1,1) fill(x1,y1,'r',x2,y2,'y') grid axis([0,5,0,4]) subplot(2,1,2) fill3(x1,y1,z1,'r',x2,y2,z2,'y') grid axis([0,5,0,4,-2,2]) subplot(2,1,1), patch(x3,y3,'g') subplot(2,1,2), patch(x3,y3,z3,'g')


0

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

2

0

-2 4

3‐14 (24) 17.8.11

3

2

1

0

0

1

2

3

4

5


3.3 Editace a export grafů Kromě možnosti měnit vlastnosti grafických objektů obecným příkazem set(hObject, PropertyName, PropertyValue). MATLAB poskytuje nástroje k editaci vytvořeného grafu přímo v grafu. Nejprve ale připomeňme základní editaci prostřednictvím příkazu set. Při jeho použití je potřeba znát název vlastnosti (PropertyName) a možné hodnoty (PropertyValue). Název vlastnosti je možné zjistit příkazem a možné hodnoty odhadnout nebo zjistit v nápovědě. Následující příklad ukazuje, jak zjistit vlastnosti objektů Line a rodičovského objektu Axes grafu vytvořeného příkazem plot. >> >> >> hl

x=0:0.1:2*pi; ys=sin(x); yc=cos(x); hl=plot(x,ys,x,yc) % objekty: 1xFigure, 1xAxes 2xLine = 174.0016 % handlery objektů Line 175.0011 >> ha=get(hl(1),'Parent') % handler rodičovského objetu Axes ha = 173.0011 >> ha=get(hl(2),'Parent') % obou objektů Line je stejný ha = 173.0011 % vlastnosti objektu Axes >> get(ha) PlotBoxAspectRatioMode = auto ActivePositionProperty = Projection = orthographic outerposition Position = ALim = [0 1] [0.13 0.11 0.775 0.815] ALimMode = auto TickLength = [0.01 0.025] AmbientLightColor = [1 1 1] TickDir = in Box = on TickDirMode = auto CameraPosition = TightInset = [3.5 0 17.3205] [0.0464286 0.0404762 CameraPositionMode = auto 0.00892857 0.0190476] CameraTarget = [3.5 0 0] Title = [176.001] CameraTargetMode = auto Units = normalized CameraUpVector = [0 1 0] View = [0 90] CameraUpVectorMode = auto XColor = [0 0 0] CameraViewAngle = [6.60861] XDir = normal CameraViewAngleMode = auto XGrid = off CLim = [0 1] XLabel = [177.001] CLimMode = auto XAxisLocation = bottom Color = [1 1 1] XLim = [0 7] CurrentPoint = XLimMode = auto [ (2 by 3) double array] XMinorGrid = off ColorOrder = XMinorTick = off [ (7 by 3) double array] XScale = linear DataAspectRatio = [3.5 1 1] XTick = DataAspectRatioMode = auto [ (1 by 8) double array] DrawMode = normal XTickLabel = 0 1 2 3 4 5 6 7 FontAngle = normal XTickLabelMode = auto FontName = Helvetica XTickMode = auto FontSize = [10] YColor = [0 0 0] FontUnits = points YDir = normal FontWeight = normal YGrid = off GridLineStyle = : YLabel = [178.001] Layer = bottom YAxisLocation = left LineStyleOrder = YLim = [-1 1] LineWidth = [0.5] YLimMode = auto MinorGridLineStyle = : YMinorGrid = off NextPlot = replace YMinorTick = off OuterPosition = [0 0 1 1] YScale = linear PlotBoxAspectRatio = [1 1 1]


3‐15 (24) 17.8.11


YTick = [ (1 by 11) double array] YTickLabel = [ (11 by 4) char array] YTickLabelMode = auto YTickMode = auto ZColor = [0 0 0] ZDir = normal ZGrid = off ZLabel = [179.001] ZLim = [-1 1] ZLimMode = auto ZMinorGrid = off ZMinorTick = off ZScale = linear ZTick = [-1 0 1] ZTickLabel = ZTickLabelMode = auto ZTickMode = auto

BeingDeleted = off ButtonDownFcn = Children = [ (2 by 1) double array] Clipping = on CreateFcn = DeleteFcn = BusyAction = queue HandleVisibility = on HitTest = on Interruptible = on Parent = [1] Selected = off SelectionHighlight = on Tag = Type = axes UIContextMenu = [] UserData = [] Visible = on

% vlastnosti objektů Line >> get(hl(1)) DisplayName: '' Annotation: [1x1 hg.Annotation] Color: [0 0 1] LineStyle: '-' LineWidth: 0.5000 Marker: 'none' MarkerSize: 6 MarkerEdgeColor: 'auto' MarkerFaceColor: 'none' XData: [1x63 double] YData: [1x63 double] ZData: [1x0 double]

>> get(hl(2)) DisplayName: '' Annotation: [1x1 hg.Annotation] Color: [0 0.5000 0] LineStyle: '-' LineWidth: 0.5000 Marker: 'none' MarkerSize: 6 MarkerEdgeColor: 'auto' MarkerFaceColor: 'none' XData: [1x63 double] YData: [1x63 double] ZData: [1x0 double]

BeingDeleted: 'off' ButtonDownFcn: [] Children: [0x1 double] Clipping: 'on' CreateFcn: [] DeleteFcn: [] BusyAction: 'queue' HandleVisibility: 'on' HitTest: 'on' Interruptible: 'on' Selected: 'off' SelectionHighlight: 'on' Tag: '' Type: 'line' UIContextMenu: [] UserData: [] Visible: 'on' Parent: 173.0011 XDataMode: 'manual' XDataSource: '' YDataSource: '' ZDataSource: ''

BeingDeleted: 'off' ButtonDownFcn: [] Children: [0x1 double] Clipping: 'on' CreateFcn: [] DeleteFcn: [] BusyAction: 'queue' HandleVisibility: 'on' HitTest: 'on' Interruptible: 'on' Selected: 'off' SelectionHighlight: 'on' Tag: '' Type: 'line' UIContextMenu: [] UserData: [] Visible: 'on' Parent: 173.0011 XDataMode: 'manual' XDataSource: '' YDataSource: '' ZDataSource: ''

Editace přímo v grafu je možná dvojím způsobem. První způsob, s omezenými možnostmi změn vlastností, je standardní způsob používaný v OS Windows. Po kliknutí na příslušný objekt v grafu pravým tlačítkem myši, výběru vlastnosti KŘP/IMSW Modelování ve výpočtových software

3‐16 (24) 17.8.11


v menu se provede změna 11 její hodnoty. Tento postup je možný až po povole‐ ní režimu editace tlačítkem 1.se šipkou Edit Plot (též Tools->Edit Plot) na nástrojové liště tlačítek (Toolbar) – viz Obrázek 3‐1.

Druhý způsob spočívá v použití nástrojů Property Editor a případně Inspector. První nástroj se spustí volbou položky menu View->Property Editor. Po jeho spuštění se objeví další okno zobrazující (a umožňující editaci) některých vlastností aktuálně vybraného objektu. Všechny vlastnosti vybrané‐ ho objektu jsou pak přístupné v nástroji Inspector, který spustíme pomocí tlačítka 2 More Properties … v pravém horním rohu okna Property Editor. Dále je možné využít nástrojů View->Plot Browser a View>Figure Palette. Záhlaví okna (objekt Figure) grafu obsahuje další funkce a nástroje. Tlačítko se šipkou do kruhu.3 Rotate 3D (též Tools->Rotate 3D)) povoluje režim otáčení 12 grafu v okně pomocí myši (i v případě 2D grafů), tlačítko 4 s ručkou Pan (též Tools->Pan) dovoluje posouvat graf atd. Položka menu Insert obsa‐ huje možnosti vložit či odebrat další objekty grafu (včetně pomocí myši vytvá‐ řených objektů 13 Line, Arrow atd.), položka Tools obsahuje (krom již zmíně‐ ných nástrojů též přístupných z tlačítkové lišty) další nástroje pro editaci grafu jako je zvětšování/zmenšování grafu, rozmisťování a zarovnávání objektů atd.

4

3

Plot Browser

Figure Palette

Object Browser

1

2 Property Editor

Obrázek 3‐1 Editace grafu

11

Výhodou tohoto postupu je, že je možné vybrat pouze přípustné hodnoty dané vlastnosti.

12

Mění interaktivně vlastnosti Camera… objektu Axes.

13

Objekty lze vytvořit také příkazem annotation(ObjectType), kde parametr ObjectType může nabýt hodnot 'rectangle', 'ellipse', 'textbox', 'line', 'arrow', 'doublearrow', nebo'‐ textarrow'.


3‐17 (24) 17.8.11


Položka menu Edit->Copy Figure kopíruje graf do schránky (ClipBoard) a umožňuje tak graf vložit do dalších programů. Vlastnosti kopírování lze nasta‐ vit v položce menu Edit->Copy Options… . Hotový graf lze exportovat do různých formátů v položce menu File->Save As… a následným výběrem požadovaného formátu. Vlastnosti exportu lze nastavit v položce menu File>Export SetUp… . Je možné uložit obrázek v následujících grafických formá‐ tech:

MATLAB Figure

.fig nativní vektorový formát, umožňuje opět obrá‐ zek načíst a dále upravovat Adobe Illustrator file .ai vektorový formát Adobe Illustrator 88 (fonty nejsou vloženy, pouze vektorové obrázky, 72 dpi) Bitmap File .bmp rastrový formát Windows Bitmap, 24‐bit barvy, formát zavedený s Windows 3.0 EPS file .eps vektorový formát Encapsulated PostScript (EPS) jediný podporuje barevný prostor CMYK, nej‐ vhodnější pro DTP Enhanced Metafile .emf vektorový formát Enhanced Windows Metafile, nejvhodnější pro import do Microsoft aplikací (nativní formát Windows) JPEG image .jpg komprimovaný (ztrátově) rastrový formát Joint Photographic Experts Group, 24‐bit barvy, implicitní komprese 75 (100 je minimální komprese), pro jinou kompresy (např. 95) je nutné použít řádkový příkaz print –djpeg95 soubor.jpg Paintbrush 24‐bit file .pcx komprimovaný rastrový formát programu Pa‐ intbrush, poměrně rozšířen Portable Bitmap File .pmb nekomprimovaný rastrový formát, 2‐bit barvy (černobílý),UNIX Portable Dokument Format .pdf formát firmy Adobe používaný pro elektronické publikace Portable Greymap file .pgm nekomprimovaný rastrový formát, 8‐bit barvy (stupně šedi),UNIX Portable Network Graphics file .png komprimovaný (bezztrátově) rastrový formát navržený pro webovou grafiku, 24‐bit barvy Portable Pixmap File .pxm nekomprimovaný rastrový formát, 24‐bit barvy (TrueColor),UNIX Portable inKmap file .pkm rastrový formát, UNIX TIFF image .tif rastrový formát Tagged Image File Format, 24‐ bit, je podporován většinou grafických programů TIFF no compression image .tif bez komprese Kromě uvedené možnosti exportu je možné použít i příkaz print, který jednak umožňuje u některých formátů nastavit parametry a jednak dovoluje další ex‐ portovat do dalších formátů např. print –dhpgl soubor.hpg vektorový formát HPGL (Hewlett Packard Gra‐

phics Language), jazyk pro ovládání ploterů (např. HP7475A) print –dhdf soubor.hdf rastrový formát Hierarchical Data Format


3‐18 (24) 17.8.11


3.4 Příklady na procvičení V kapitole 3.4.1 jsou zadání a v kapitole 3.4.2 řešení jednoduchých příkladů využívající znalosti z prvních tří bloků. Předpokládá se, že při řešení příkladů má čtenář k dispozici MATLAB a aktivně využívá nápovědu k programu.

3.4.1 Zadání příkladů Př01. Řešení soustavy lineárních rovnic. Je dána následující soustava lineárních rovnic 3x1‐2x2+3x4=1‐x1 2x2‐3x3+2x4=2‐x2 ‐x1+3x2+2x3=3‐x3 4x1 ‐x3 +x4=4‐x4 nalezněte řešení a) použitím inverzní matice b) pomocí levého maticového dělení c) pomocí Crammerova pravidla d) využitím LU dekompozice a následnou manuální dvojnásobnou Gaussovou eliminací Př02. Vytvoření grafu aproximace zadaných bodů. Jsou dány body bi={xi,yi},.i=1,..,5 b={[‐2,7], [‐1,7], [0.5,1], [1,3], [3,‐1]} a) určete hodnoty parametrů k, q přímky popsané rovnicí y=k*x+q takové, kte‐ rá dané body prokládá nejlépe ve smyslu minima sumy kvadrátů odchylek b) do grafu vyneste původní body b jako červená kolečka velikosti 12 bodů a 100 bodů odpovídajících nalezené přímce (čárkovaná modrá čára), tloušťka čar 3 body c) na ose x nastavte rozsah od ‐5 do +5 a na ose y od ‐5 do 10 a graf opatřete mřížkou, nadpisem (14 bodů Bold), legendou a popisem os (12 bodů) Př03 Tvorba různých druhů grafů, využití handlerů objektů. Vytvořte řádkový vektor vr naplněný 1000 náhodnými čísly rovnoměrně rozlo‐ ženými v rozsahu ‐1 až 1 a řádkový vektor vn naplněný 1000 náhodnými čísly normální rozložení N(0,1) a vytvořte vektor v tvořený součinem odpovídají‐ cích hodnot ve vektorech vn, vr. a) vytvořte jedno okno obsahující 6 grafů uspořádaných do dvou řad a třech sloupců, handlery jednotlivých objektů uložte do matice H odpovídající struk‐ tury b) s využitím matice handlerů H vytvořte grafy hodnot ve vektorech vr, vn a v (hodnoty ve vektorech jsou y‐nové souřadnice, za x‐ové souřadnice považu‐ jeme index hodnoty ve vektoru) umístěné v horní řadě grafů. Hodnoty jsou znázorněny jako izolované body kolečko. Každý graf označte nadpisem „Rov‐ noměrné“, „Gaussovo“ a „Součin“. Nastavte měřítko v ose x od 0 do 1000, v ose y od ‐4 do +4. c) s využitím matice handlerů H vytvořte histogramy (12 sloupců) hodnot ve vektorech vr, vn a v umístěné v dolní řadě grafů. KŘP/IMSW Modelování ve výpočtových software

3‐19 (24) 17.8.11


Př04 Vytvoření chybového grafu, práce s grafem. Vytvořte vektor x (10 hodnot na intervalu ‐2≤x≤+2), spočítejte vektor y podle vztahu y(x)=x5‐4x4+x+1. Vytvořte vektor e obsahující 10 náhodných čísel s normálním rozložením N(0,1). Do jednoho grafu vyneste pospojované pů‐ vodní body (červené hvězdičky velikosti 16, síla čáry 2) s chybovými úsečkami ± dvojnásobek směrodatné odchylky čísel ve vektoru e a samostatné body získané jako součet vektorů y a e (modrá kolečka velikosti 16 a tloušťka čáry 2). Graf opatřete mřížkou. Př05 Vykreslení průběhů matematických funkcí. Vytvořte jedno okno obsahující čtyři grafy uspořádané do dvou řádků a dvou sloupců. Grafy budou obsahovat průběhy dále uvedených funkcí (síla čary 2) a budou opatřeny nadpisem s názvem grafu a) Descartův list

x(t ) =

3at 1+ t3

y (t ) =

složený ze dvou částí první část (modrá) ‐0.4 ≤ t1 ≤ 999.5 druhá část (červená) ‐1002 ≤ t2 ≤ ‐2.1 b) Bernouliova lemniskáta

x(t ) =

(

3at 2 a=2 1+ t3 1000 bodů logaritmicky 1000 bodů logaritmicky

at 1 + t 2 1+ t4

)

y (t ) =

(

)

at 1 − t 2 1+ t4

a=2

‐ 10 ≤ t ≤ 10 1000 bodů r=1 c) Cykloida r=1 x(t ) = rt − d sin (t ) y (t ) = r − d cos(t ) 1000 bodů 0 ≤ t ≤ 4π tři průběhy pro d=2.5 (červený), d=1 (modrý) a d=0.5 (zelený), doplnit legen‐ du c) Epicykloida

⎛R+r ⎞ t⎟ x(t ) = (R + r ) cos(t ) − d cos⎜ ⎝ r ⎠

⎛R+r ⎞ t⎟ y (t ) = (R + r )sin(t ) − d cos⎜ ⎝ r ⎠

pro R=12, r=4 , d=3 0 ≤ t ≤ 2π

1000 bodů

Př06 Tvorba 3D grafů, práce s handlery grafických objektů. Vytvořte jeden graf obsahující dvě plochy popsané funkcemi

z1 ( x, y ) =

xy 1 + x2 + y2

z2 ( x, y ) =

4 cos( xy 3) 1 + x2 + y 2

První plochu vykreslete jako drátový model a druhou s vyplněnými plochami a doplňte popisy os (řez písma Bold). V aktivním objektu Axes nalezněte handlery všech objektů typu (vlastnost „Type“) plocha (hodnota „Surface“) a ploše odpovídající druhé funkci nastavte Phongovo stínování („FaceLighting“, „phong“) a interpolaci barev („FaceColor“, „interp“). Vytvořte osvětlení v bodu [‐2,2,3] (příkaz light) a nastavte způsob výpočtu osvětlení (příkaz lighting) na Goraudův


3‐20 (24) 17.8.11


3.4.2 Řešení příkladů Př01. Řešení soustavy lineárních rovnic. a)

b)

>> >> >> xa

>> x1 >> x2 >> x3 >> x4

A=[4,-2,0,3; 0,3,-3,2;-1,3,3,0;4,0,1,5]; b=[1;2;3;4]; xa=inv(A)*b b) >> xb=A\b = 0.0390 xb = 0.0390 0.6104 0.6104 0.4026 0.4026 0.6883 0.6883 x1=det([b,A(:,2:4)])/det(A) = 0.0390 x2=det([A(:,1),b,A(:,3:4)])/det(A); = 0.6104 x3=det([A(:,1:2),b,A(:,4)])/det(A); = 0.4026 x4=det([A(:,1:3),b])/det(A) = 0.6883

Ax = b A = LU d)

LUx { = b ⇒ Ly = b a Ux = y y

>> [L,U]=lu(A) % dolni a horni troj.matice L = 1.0000 0 0 0 0 1.0000 0 0 -0.2500 0.8333 1.0000 0 1.0000 0.6667 0.5455 1.0000 U = 4.0000 -2.0000 0 3.0000 0 3.0000 -3.0000 2.0000 0 0 5.5000 -0.9167 0 0 0 1.1667 >> y=zeros(4,1); x=y; % sloupce >> % Gaussova eliminace L*y=b >> y(1)=b(1)/L(1,1); >> y(2)=(b(2)-L(2,1)*y(1))/L(2,2); >> y(3)=(b(3)-L(3,1)*y(1)-L(3,2)*y(2))/L(3,3); >> y(4)=(b(4)-L(4,1)*y(1)-L(4,2)*y(2)-L(4,3)*y(3))/L(4,4); >> y y = 1.0000 2.0000 1.5833 0.8030 >> % Gaussova eliminace U*x=y >> x(4)=y(4)/U(4,4); >> x(3)=(y(3)-U(3,4)*x(4))/U(3,3); >> x(2)=(y(2)-U(2,4)*x(4)-U(2,3)*x(3))/U(2,2); >> x(1)=(y(1)-U(1,4)*x(4)-U(1,3)*x(3)-U(1,2)*x(2))/U(1,1); >> x x = 0.0390 0.6104 0.4026 0.6883 ⎡− 2 1⎤ ⎡7⎤

Př02. Vytvoření grafu aproximace zadaných bodů.

b={[‐2,7], [‐1,7], [0.5,1], [1,3], [3,‐1]} y=kx+q 7=‐2k+q 7=‐k+q 1=0.5k+q


3‐21 (24) 17.8.11

⎢ −1 ⎢ ⎢0.5 ⎢ ⎢1 ⎢⎣ 3

⎢7⎥ 1⎥⎥ ⎡k ⎤ ⎢ ⎥ 1⎥ ⎢ ⎥ = ⎢ 1 ⎥ ⎥ ⎣q ⎦ ⎢ ⎥ 1⎥ ⎢3⎥ ⎢⎣− 1⎥⎦ 1⎥⎦ František Dušek KŘP FEI Univerzita Pardubice

b)

c)

3=k+q ‐1=3k+q

>> x=[-2,-1,0.5,1,3]; 10 >> y=[7,7,1,3,-1]; >> A=[x',ones(5,1)]; >> p=A\y'; >> k=p(1), q=p(2) 5 k = -1.7297 q = 3.9189 >> xx=linspace(min(x),max(x),100); >> yy=k*xx+q; >> hl=plot(x,y,'or',xx,yy,'b--'); 0 >> set(hl(1),'MarkerSize',12) >> set(hl,'LineWidth',3) >> axis([-5,5,-5,10]) >> grid -5 -5 -4 >> legend('body','přímka') >> ht=title('Proložení přímkou'); >> set(ht,'FontSize',14,'FontWeight','Bold') >> xlabel('osa x','FontSize',12) >> ylabel('osa y','FontSize',12)

Proložení přímkou body přímka

osa y

a)

-3

-2

-1

0

1

2

3

4

5

osa x

Př03 Tvorba různých druhů grafů, využití handlerů objektů. a) b) c)

>> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >>

vr=20*rand(1000,1)-10; vn=randn(1000,1); v=vr.*vn; H(1,1)=subplot(2,3,1); H(1,2)=subplot(2,3,2); H(1,3)=subplot(2,3,3); H(2,1)=subplot(2,3,4); H(2,2)=subplot(2,3,5); H(2,3)=subplot(2,3,6); plot(H(1,1),vr,'o') title(H(1,1),'Rovnoměrné') plot(H(1,2),vn,'o') title(H(1,2),'Gaussovo') plot(H(1,3),v,'o') title(H(1,3),'součin') axis(H(1,:),[0,1000,-4,4]) hist(H(2,1),vr,12) hist(H(2,2),vn,12) hist(H(2,3),v,12)

Rovnoměrné

Gaussovo

součin

4

4

4

2

2

2

0

0

0

-2

-2

-2

-4

0

500

1000

120

-4

0

500

-4

1000

200

400

150

300

100

200

50

100

0

500

1000

0

2

100 80 60 40 20 0 -1

0

0 -5

1

0

0 -2

5

Př04 Vytvoření chybového grafu, práce s grafem.

>> x=linspace(-1.2,1.2,10); >> y=x.^5+x.^4-4*x.^2+1; >> e=randn(1,10); >> ep=zeros(size(y))+std(e); >> h=errorbar(x,y,2*ep,'ro-'); >> set(h,'LineWidth',2,… 'MarkerSize',16) >> hold on >> h=plot(x,y+e,'b*'); >> set(h,'LineWidth',2,… 'MarkerSize',16) >> grid

4

2

0

-2

-4

-6

-8 -1.5


3‐22 (24) 17.8.11

-1

-0.5

0

0.5

1

1.5


Př05 Vykreslení průběhů matematických funkcí. >> subplot(2,2,1), a=2; a)

b)

c)

d)

>> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >> >>

t1=logspace(-1,3,1000)-0.5; %-0.49 az 999.5 x1=3*a*t1./(1+t1.^3); y1=3*a*t1.^2./(1+t1.^3); plot(x1,y1,'b','LineWidth',2) hold on t2=-logspace(-1,3,1000)-2; %-2.1 az -1002 x2=3*a*t2./(1+t2.^3); y2=3*a*t2.^2./(1+t2.^3); plot(x2,y2,'r','LineWidth',2) title('Descartuv list'), grid subplot(2,2,2), a=2; t=linspace(-10,10,1000); a=2; x=a*t.*(1+t.^2)./(1+t.^4); y=a*t.*(1-t.^2)./(1+t.^4); plot(x,y,'LineWidth',2) title('Bernouliova lemniskáta') grid subplot(2,2,3) Descartuv list 4 t=linspace(0,4*pi,1000); 2 r=1; d=2*r; x=r*t-d*sin(t); y=r-d*cos(t); 0 plot(x,y,'b','LineWidth',2) -2 hold on d=r; -4 -4 -2 0 x=r*t-d*sin(t); y=r-d*cos(t); plot(x,y,'r','LineWidth',2) Cykloida 3 d=r/2; 2 x=r*t-d*sin(t); y=r-d*cos(t); plot(x,y,'g','LineWidth',2) 1 title('Cykloida'),grid 0 legend('d=2r','d=r','d=r/2') -1 subplot(2,2,4) -5 0 5 R=5; r=1; d=r*2; Rrt=(R+r)/r*t; x=(R+r)*cos(t)-d*cos(Rrt); y=(R+r)*sin(t)-d*sin(Rrt); plot(x,y,'LineWidth',2) title('Cykloida'), grid

Bernouliova lemniskáta 1

0.5

0

-0.5

2

4

-1 -2

-1

0

1

2

5

10

Cykloida 10 d=2r d=r d=r/2

5

0

-5

10

15

-10 -10

-5

0

Př06 Tvorba 3D grafů, práce s handlery grafických objektů.

x=-4:0.2:4; y=-4:0.2:4; [X,Y]=meshgrid(x,y); Z=X.*Y./sqrt(1+X.^2+Y.^2); mesh(x,y,Z), hold on Z=4*cos(X.*Y/3)./sqrt(1+X.^2+Y.^2); surf(x,y,Z) xlabel('osa x','FontWeight','Bold') ylabel('osa y','FontWeight','Bold') zlabel('osa z','FontWeight','Bold') h=findobj(gca,'Type','surface'); set(h(1),'FaceLighting','phong', ... 'FaceColor','interp') >> light('Position',[-2 2 3]); >> lighting gouraud

>> >> >> >> >> >> >> >> >> >> >>


3‐23 (24) 17.8.11


Pojmy k zapamatování Okruhy problémů: soustavy lineárních rovnic, základní 2D a 3D grafy, spe‐ ciální grafy, grafické objekty, editace grafů, export grafů Použité nástroje:

Property Editor, Inspector, Plot Browser, Figure Pallete

Příkazy a funkce: axes, axis, bar, barh, clf, contour, contourf, det, errorbar, figure, fill, fill3, gca, gcf, grid, hist, hold, inv, legend, lu, mesh, meshgrid, norm, patch, pie, pinv, plot, plot3, print, rank, set, stairs, surf, title, xlabel, ylabel

Otázky na procvičení 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.

Co je to přeurčená soustava lineárních rovnic, co se považuje za její řešení a jak se získá v MATLABu? Jak zjistíme v MATLABu hodnost matice a k čemu ji využijeme? Jakou funkcí MATLABu zjistíme, že matice je singulární? Co se označuje pojmem dekompozice matice a jaká dekompozice je využi‐ ta při řešení soustavy n lineárních rovnic o n neznámých? Je dána funkce popisující plochu v prostoru z(x,y)=f(x,y). Jak připravíte matici Z pro příkaz mesh(x,y,Z)? Co označuje pojem Handle Graphics? K čemu slouží funkce get a set? Jak se získá handle aktivního okna? Jaká funkce usnadňuje vytvoření matice s polohou bodů pro vykreslení 3D grafu pomocí funkce mesh či surf? Jak vložíme graf vytvořený v MATLABu např. do dokumentu Wordu?

Odkazy a další studijní prameny •

•

on‐line dokumentace k programu nebo www.mathworks.com/help/ – části „Systems of Linear Equations“, „Basic Plotting Commands“, „Crea‐ ting 3‐D Graphs“, „Creating Specialized Plots“, „Organization of Graphics Objects“, „Accessing Object Handles“, „Editing Plots“, „Overview of Prin‐ ting and Exporting“ elektronická učebnice Learning MATLAB 7– www.mathworks.com/academia/student_version/learnmatlab_sp3.pdf (část „Graphics“)


3‐24 (24) 17.8.11


3 Soustavy lineárních rovnic, tvorba grafů

Recommend Documents