1
Cviˇcení 1a - Úvod od programování v Matlabu/Octave
1.1 1.1.1
Matlab Prostˇredí
ˇ Pro studenty CVUT je dostupná studentská licence ke stažení z https://download.cvut.cz/ (po ˇ pˇrihlášení). Nevýhodou této licence je nutnost pˇripojení z IP adresy CVUT.
Na obrázku vidíte pravdˇepodobné defaultní nastavení Matlabu. Toto hlavní okno mužete ˚ ovládat pomocí funkcí pod prvními tˇremi záložkami (oznaˇceno žlutým rámeˇckem). Aktuálnˇe je zapnuta záložka home, s její nabídkou (modrý rámeˇcek) bychom si mˇeli pro naši potˇrebu vystaˇcit. Zeleným rámeˇckem je oznaˇcena cesta k aktuálnímu pracovnímu adresáˇri, kam se nám bude veškerá práce ukládat. V dolní cˇ ásti okna vidíme: * Current Folder - pruzkumník ˚ souboru˚ * Command Window - obdoba konzole, jednoduché operace lze provádˇet pˇrímo zde * Workspace seznam promˇenných aktuálnˇe uložených v pamˇeti * Editor V editoru je aktuálnˇe otevˇrený skript example01.m (všechny matlabovské funkce a skripty mají koncovku *.m). V editoru mužeme ˚ skripty a funkce upravovat, spouštˇet a debugovat. Nabídka editoru je rozdˇelena do tˇrí založek oznaˇcených cˇ erveným rámeˇckem. Nejduležitˇ ˚ ejší pro nás bude záložka EDITOR.
1
1.1.2
Zpusob ˚ práce
V Matlabu lze pracovat dvˇema zpusoby: ˚ • Práce v Command Window – pˇríkazy zapisujeme pˇrimo do Command Window, Matlab vyhodnocuje operace postupnˇe, tak jak je “odesíláme” – vhodné pouze pro krátké výpoˇcty (jednorázové rˇ ešení soustavy rovnic apod.) • Práce v Editoru – – – – – 1.1.3
umožnuje vytvoˇrit dávku pˇríkazu, ˚ které zapíšeme do skriptu (a funkcí) skript potom spustíme jako celek, výstup se vypíše do Command Window lze pˇridávat komentáˇre skript je možné debugovat lze docílit pˇrehledného zápisu, ke kterému se dá snadno vracet
Dobré vˇedˇet=)
• Klávesa F1 – pomocí F1 spustíte nápovˇedu Matlabu, ve které celkem snadno najdete vše, co budete potˇrebovat • Velikost písma – HOME -> Preferences -> Fonts • Klávesové zkratky – HOME -> Preferences -> Keybord -> Shortcuts • Stˇredníky – stˇredník v editoru potlaˇcuje výstup – pokud pˇríkaz (definici promˇenné, pˇriˇrazení hodnoty promˇenné, volání funkce, . . . ) nezakonˇcíte stˇredníkem, vypíše se celý obsah do Command Window ⇒ nepˇrehlednost, problémy pˇri velkém objemu dat – naopak hodnoty promˇenných, které nás zajímají (výsledky, mezivýsledky, . . . ), mužeme ˚ vypsat tak, že stˇredník neuvedeme • Komentáˇre – komentáˇre v Matlabu zapisujeme za znak % – na zaˇcátek každého skriptu/funkce je vhodné okomentovat, k cˇ emu slouží, jaké jsou vstupy a výstupy ⇒ pokud se budete chtít k práci vrátit po delší dobˇe nebo budete skript sdílet s kolegy, snadneji se zorientujete • Tˇri teˇcky . . . – pro pˇrehlednost kódu mužete ˚ použít odˇrádkování s využitím tˇrí teˇcek . . . 2
– napˇr. pokud zapisujete delší vzorec, mužete ˚ ukonˇcit jeden rˇ ádek pomocí tˇrí teˇcek a pokraˇcovat na dalším rˇ ádku • Jak poznáte, že Matlab stále pracuje? – jakmile spustíte skript, objeví se v levém dolním rohu “Busy” (na úvodním obrázku ružový ˚ rámeˇcek) – po dokonˇcení posledního pˇríkazu ze spuštˇeného skriptu, tento nápis zmizí • Nˇeco je špatnˇe, Matlab je stále “Busy” – vykonání “složitˇejších” programu˚ muže ˚ být samozˇrejmˇe cˇ asovˇe nároˇcné, pokud ale proces trvá déle než byste pˇredpokládali, muže ˚ být chyba nˇekde jinde (nekoneˇcný cyklus, pˇrehnanˇe velká soustava, . . . ) – pokud chcete takový proces ukonˇcit, staˇcí pˇrepnout do Command Window a stisknout Ctrl+C • Chybová hlášení – pˇri spuštˇení skriptu, ve kterém je nˇejaká syntaktická chyba nebo neplatná operace (napˇr. násobení matic nekompatibilních velikostí), vypíše Matlab do Command Window chybová hlášení – v chybovém hlášení je popsáno, v cˇ em je chyba, a jsou uvedeny odkazy na rˇ ádky kódu, kde se chyba nachází • File is not found . . . – pokud jste doposud pracovali v nˇejakém adresáˇri (cesta oznaˇcená na prvním obrázku zelenˇe) a pokusíte se spustit skript, který máte uložený na jiném místˇe, dostanete následující hlášení
– mužete ˚ si vybrat ze dvou variant * volbou Change Folder zmˇeníte aktuální adresáˇr na místo, kde se nachází spouštˇený soubor, cesta k pracovnímu adresáˇri se zmˇení a veškerá další práce bude probíhat v této složce * volbou Add to Path pˇridáte cestu ke spouštˇenému skriptu, ale cesta k pracovnímu adresáˇri zustane ˚ beze zmˇeny * tip: volbu Add to Path používejte pouze pokud spouštˇený skript pˇrímo souvisí s funkcemi, které máte uložené na místˇe, které bylo do ted’ pracovním adresáˇrem, v opaˇcném pˇrípadˇe pˇrepnˇete do složky, ze které skript voláte (všechny novˇe vytvoˇrené funkce se budou ukládat tam a soubory, které spolu souvisí zustanou ˚ pohromadˇe)
3
• Výsledkem je Inf / NaN – Matlab muže ˚ jako výsledek nˇejaké cˇ íselné operace vrátit neˇcíselné hodnoty Inf nebo NaN * Inf (infinity) - zastupuje nekoneˇcno (napˇr. pokud dojde k dˇelení nulou) * NaN (not a number) - zastupuje hodnoty, které nejsou reálnými ani komplexními cˇ ísly (napˇr. pokud dojde k dˇelení nuly nulou)
1.2 1.2.1
Octave Prostˇredí
Octave je volnˇe dostupné ke stažení z https://www.gnu.org/software/octave/. Octave lze ovládat z pˇríkazové rˇ ádky, ale je dostupné i grafické rozhraní, které má obdobné uspoˇrádání jako Matlab. Octave je z velké cˇ ásti kompatibilní s Matlabem, pˇríkazy, které budeme používat, by mˇely fungovat pˇri spuštˇení v Matlabu i Octave stejnˇe.
Na obrázku je ukázáno okno Octave, cˇ erveným rámeˇckem je oznaˇceno pˇrepínání mezi Editorem a Command Window. Práce v tomto prostˇredí je obdobná práci v Matlabu.
4
1.3
Základy programování v Matlabu/Octave
1.3.1
Skaláry, vektory, matice
Název Matlab (“Matrix laboratory”) napovídá, že toto prostˇredí je vhodné pro práci s maticemi, které jsou klíˇcovou datovou strukturou tohoto prostˇredí. Jako matice (nxn) jsou v Matlabu vnímány i skaláry (1x1) a vektory (ˇrádkový 1xn, sloupcový nx1). In [1]: % scalar a = 1 % vector 1x3 v = [1, 2, 3] % matrix 2x3 M = [1, 2, 3; 2, 3, 4] a = v =
1
1
2
3
2 3
3 4
M = 1 2
Matice zapisujeme do hranatých závorek, cˇ árkou (nebo mezerou) oddˇelujeme sloupce, stˇredníkem ukonˇcujeme jednotlivé rˇ ádky matice. S využitím tˇechto znalostí mužeme ˚ vytvoˇrit i rˇ ádkový a sloupcový vektor. In [2]: vrow = [1, 2, 3] vrow = [1 2 3] vcol = [4; 5; 6]
% row vector % row vector % column vector
vrow = 1
2
3
2
3
vrow = 1 vcol = 4 5 5
6
Transpozici získáme pomocí apostrofu. In [3]: vt = vcol' Mt = M'
% vcol = [4; 5; 6] % m = [1, 2, 3; 2, 3, 4]
vt = 4
5
6
Mt = 1 2 3
2 3 4
K jednotlivým prvkum ˚ matic (vektoru) ˚ pˇristupujeme pomocí kulatých závorek a mužeme ˚ tak pˇriˇrazovat prvkum ˚ matice nové hodnoty. K prvkum ˚ matice mužeme ˚ pˇristupovat pomocí indexu˚ rˇ ádku a sloupce (ˇr,s), ménˇe cˇ asté je použití jednoho celkového indexu. Pomocí operátoru : lze pˇristupovat ke všem indexum ˚ daného rˇ ádku/sloupce nebo vybrat urˇcitý rozsah. Pozn: Matlab indexuje od 1 (narozdíl napˇr. od jazyka C, kde je prvním indexem 0). In [4]: a = vt(3) b = Mt(3,2) c = Mt(6) Mcol2 = Mt(:,2) Mrow12 = Mt(1:2,:) Mt(3,2) = 1; Mt Mt(:,1) = [0 0 0]'; Mt a = 6 b = 4 c = 4 Mcol2 = 2 3 4
6
Mrow12 = 1 2
2 3
Mt = 1 2 3
2 3 1
Mt = 0 0 0
2 3 1
Operátor : lze také využít k vytvoˇrení vektoru s ekvidistantnˇe vzdálenými prvky. Defaultním krokem je 1, lze ho zmˇenit uvedením kroku mezi první a poslední požadovaný prvek. In [5]: v1 = 1:5 v2 = 3:-0.25:2 v1 = 1
2
3
4
5
v2 = 3.0000
2.7500
2.5000
2.2500
2.0000
Matice a vektory lze spojovat (v pˇrípadˇe, že mají kompatibilní velikosti). In [6]: M1 = [v1;v2] M2 = [v1 v1] M3 = [M1;M1;v2] M4 = [M3,v1'] M1 = 1.0000 3.0000
2.0000 2.7500
3.0000 2.5000
4.0000 2.2500
M2 = 7
5.0000 2.0000
1
2
3
4
5
1
2
3
4
5
M3 = 1.0000 3.0000 1.0000 3.0000 3.0000
2.0000 2.7500 2.0000 2.7500 2.7500
3.0000 2.5000 3.0000 2.5000 2.5000
4.0000 2.2500 4.0000 2.2500 2.2500
5.0000 2.0000 5.0000 2.0000 2.0000
2.0000 2.7500 2.0000 2.7500 2.7500
3.0000 2.5000 3.0000 2.5000 2.5000
4.0000 2.2500 4.0000 2.2500 2.2500
5.0000 2.0000 5.0000 2.0000 2.0000
M4 = 1.0000 3.0000 1.0000 3.0000 3.0000
1.0000 2.0000 3.0000 4.0000 5.0000
Velikost matice lze zjistit pomocí funkce size, pro vektor ji lze také použít, ale existuje i speciální funkce length. In [7]: sM4 = size(M4) rM4 = size(M4,1) cM4 = size(M4,2) cM2 = size(M2,2) sv1 = length(v1) sM4 = 5 rM4 cM4 cM2 sv1
6 = = = =
5 6 10 5
ˇ Rádek/sloupec matice lze vymazat pˇriˇrazením prázdného vektoru []. In [8]: M4(:,2:4) = [] M4 = 1 3
5 2
1 2 8
1 3 3
5 2 2
3 4 5
Matlab umí automaticky vytvoˇrit nˇekteré speciální matice: nulovou (zeros), jedniˇckovou (ones), jednotkovou (eye). In [9]: Mzeros = zeros(2,3) Mones = ones(2,3) Meye = eye(2,3) Mzeros = 0 0
0 0
0 0
Mones = 1 1
1 1
1 1
Meye = Diagonal Matrix 1 0
0 1
0 0
Matice lze mezi sebou scitat, odcitat, nasobit, delit, nasobit skalarem. Pozor na kompatibilitu velikostí. In [10]: M1 = [1 2; 3 4]; M2 = [3 4; 5 6]; M1 + M2 M1 - M2 M1 * M2 ans = 4 8
6 10
ans =
9
-2 -2
-2 -2
ans = 13 29
16 36
Teˇckový operátor lze použít k násobení (dˇelení, mocnˇení, . . . ) po prvcích. In [11]: M1 .* M2 ans = 3 15
8 24
Matice lze násobit a dˇelit skalárem. In [12]: M3 = ones(2,3) * 6 M4 = M1 / 3 M3 = 6 6
6 6
6 6
M4 = 0.33333 1.00000
0.66667 1.33333
Pro výpis vlastních cˇ ísel (a vlastních vektoru) ˚ lze použít funkci eig().
In [14]: eigvals = eig(M1) % eigenvalues [vecs, vals] = eig(M1) % eigenvectors (in columns) and eigenvalues store eigvals = -0.37228 5.37228 vecs = 10
-0.82456 0.56577
-0.41597 -0.90938
vals = Diagonal Matrix -0.37228 0
0 5.37228
Inverzní matici získáme použitím funkce inv(). In [13]: iM1 = inv(M1) iM1 * M1 % check iM1 = -2.00000 1.50000
1.00000 -0.50000
ans = 1.00000 0.00000
0.00000 1.00000
ˇ 1.3.2 Rešení soustavy Mˇejme soustavu lineárních rovnic Ax = b, kde In [20]: A = [1 2 3; -1 0 2; 1 3 1] b = [1 0 0]' A = 1 -1 1
2 0 3
3 2 1
b = 1 0 0
11
Pro rˇ ešení této soustavy lze využít inverzní matici x = A−1 b. Vhodnˇejší je ale využít operátor , kterým je v Matlabu implementována Gaussova eliminace. In [21]: % inverse x = inv(A)*b % Gauss elimination x = A\b % check A*x-b x = 0.66667 -0.33333 0.33333 x = 0.66667 -0.33333 0.33333 ans = 0.0000e+00 1.1102e-16 5.5511e-17
1.3.3
Cykly a podmínky
Duležitým ˚ nástrojem pro psaní programu jsou cykly a podmínky. Cykly umožnují ˇ provádˇet cˇ ásti kódu opakovanˇe (napˇr. sestavit matici tuhosti pro všechny pruty konstrukce). Nám se bude hodit zejména for cyklus, který je proveden pro pˇredem známý poˇcet opakování. Dalším cyklem je while cyklus, který je vykonáván, dokud je zajištˇena podmínka. Pomocí podmínek mužeme ˚ cˇ lenit program do vˇetví a provádˇet tak cˇ ásti kódu pouze pˇri splnˇení nˇejaké podmínky (napˇr. funkce bude vracet v závislosti na poˇctu Gaussových bodu, ˚ který zadáme, odpovídající souˇradnice a váhy tˇechto bodu). ˚ Pro podmínky používáme klasické relaˇcní operátory >, <, >=, <=, == (rovná se), ∼= (nerovná se), || (nebo), && (a zároven). ˇ ˇ Castá chyba: Nezamˇenujte ˇ operátor pˇriˇrazení = s relaˇcním operátorem ==. Pokud zapíšete a = 1 nejedná se o logický výraz, ale do promˇenné a uložíte hodnotu 1. Naopak výraz a == 1 má návratovou hodnotu true (1), pokud je a rovno jedné, a f alse (0), pokud se a jedné nerovná. for - provede pˇredem známý poˇcet opakování
12
In [1]: % number of cycles n = 5; % allocation of empty vector of size n v = zeros(n, 1); for i = 1:n v(i) = 2*i; end % save transposition of vector v into v_transpose v_transpose = v' j = 1; sum_j = 0; while (j<10) sum_j = sum_j + j; j = j*2; end sum_j v_transpose = 2 sum_j =
4
6
8
10
15
if (elseif) - pokud platí podmínka, provede dané pˇríkazy In [2]: a b c d
= = = =
1; 2; 3; 1;
% Condition 1 disp('C1:'); if (a == 1) disp('a is equal 1'); % function disp('...') prints string in the argum end % Condition 2 disp('C2:'); if (a > d) disp('a is greater than d'); elseif (a < d) disp('a is less than d'); 13
else disp('a is equal to d'); end % Condition 3 disp('C3:'); if (a ~= 1 || c <=3) disp('a is not equal to 1 OR c is less or equal to 3') end % Condition 4 disp('C4:'); if (a == 1 && c <=3) disp('a is equal to 1 AND c is less or equal to 3') end C1: a is C2: a is C3: a is C4: a is
equal 1 equal to d not equal to 1 OR c is less or equal to 3 equal to 1 AND c is less or equal to 3
switch - obdoba if , vhodné pokud chceme v závislosti na jedné promˇenné rozlišit více pˇrípadu˚ (vyhneme se nˇekolikanásobnému použití elseif ) In [17]: type = 'beam'; % string switch type case 'truss' disp('Type of structure: truss') case 'beam' disp('Type of structure: beam') otherwise disp('Unknown type of structure.') end Type of structure: beam
break, continue - tyto pˇríkazy slouží k ukonˇcení aktuálního cyklu (break) nebo k ukonˇcení aktuální iterace a pokraˇcovaní další iterací (continue). In [18]: % break disp ('Use of break in a cycle.'); i = 0; 14
while 1 i = i+1; if i == 4 break; end disp(i); end % continue disp ('Use of continue in a cycle.'); i = 0; while i < 5 i = i+1; if i == 4 continue; end fprintf('%d ',i); % another way how to print the results end Use of break in a cycle. 1 2 3 Use of continue in a cycle. 1 2 3 5
1.3.4
Lokalizace
Pˇri výpoˇctu metodou koneˇcných prvku˚ budeme potˇrebovat uspoˇrádat (lokalizovat) jednotlivé matice tuhosti prvku˚ do celkové matice tuhosti konstrukce, tento proces nazýváme lokalizací. Zpusob, ˚ který budeme používat v našich cviˇceních, využívá lokalizaˇcní matici (loc), kde na každý rˇ ádek zapisujeme kódová cˇ ísla stupnˇ u˚ volnosti elementu. In [3]: loc =
[1 2 3 4; 3 4 5 6; 5 6 7 8]
K = zeros(8,8); for i = 1:3 k = ones(4,4) * i; % element stiffness matrix K(loc(i,:),loc(i,:)) = K(loc(i,:),loc(i,:)) + k; end K loc = 1
2
3
4 15
% global stiffness ma
3 5
4 6
5 7
6 8
1 1 1 1 0 0 0 0
1 1 3 3 2 2 0 0
1 1 3 3 2 2 0 0
K = 1 1 1 1 0 0 0 0
1.3.5
0 0 2 2 5 5 3 3
0 0 2 2 5 5 3 3
0 0 0 0 3 3 3 3
0 0 0 0 3 3 3 3
Funkce
Pro psaní rozsáhlejších programu˚ je vhodné rozˇclenit program na jednotlivé funkce. Kód se tak stává pˇrehlednˇejší a jeho cˇ ásti mužeme ˚ volat opakovanˇe. Funkce jsou ukládány v samostatných souborech (*.m), Matlab požaduje, aby se název souboru shodoval se jménem funkce (jinak hlásí varování). Na prvním rˇ ádku m-souboru je funkce deklarována, následuje nˇekolik komentáˇrových rˇ ádku˚ (zaˇcínají %). Komentované rˇ ádky jsou Matlabem vnímány jako nápovˇeda funkce a lze ji zobrazit voláním help jmenofunkce. Funkce je zakonˇcena end (nepovinné).
In [22]: function [m1, m2] = magicNumbers(num1, num2) % % This function returns two "magic numbers" calculated from two input numb % m1 = floor((num1*10 - num2 * 3) / 6) - 3; m2 = ceil((num1*0.1 + num2 / 2) * 4) + 8; end Funkce je potom volána ze skriptu (nebo z jiné funkce) následovnˇe: In [23]: [magic1, magic2] = magicNumbers(8, 5) a = -8; b = 0; [magic3, magic4] = magicNumbers(a, b) magic1 magic2 magic3 magic4
= 7 = 22 = -17 = 5
16
% function floor() rounds d % function ceil() rounds up
1.3.6
Grafický výstup
Matlab umožnuje ˇ snadnou tvorbu grafu. ˚ Pro vykreslení nˇejaké funkce je vhodná funkce plot. Do jednoho grafu mužeme ˚ nakreslit libovolný poˇcet funkcí, musíme však použít pˇríkaz hold on, který zajistí uchování i pˇredchozí funkce. Volitelnými argumenty funkce plot je barva, vzory bodu˚ a cˇ ar (naleznete v nápovˇedˇe). Pro tisk bodu˚ lze použít funkci scatter. Grafický výstup je v Matlabu kreslen do zvláštních oken (figures), nové okno otevˇreme pˇríkazem figure, pokud chceme všechna otevˇrená okna pozavírat použijeme pˇríkaz close all. Okno s grafem lze rozˇclenit na podokna pomocí funkce subplot. Rozsah hodnot na osách lze nastavit pˇríkazem axis, osy popíšeme pomocí pˇríkazu˚ xlabel, ylabel. Graf lze nadepsat pomocí pˇríkazu title. In [24]: subplot(1, 2, 1) x = [0:0.05:1] * pi; y = cos(x); plot(x,y); hold on; y = sin(x); plot(x,y, '--r*'); subplot(1, 2, 2) scatter(x,y, 'k'); axis([0 pi -1 1]); xlabel('x'); ylabel('sin(x)'); title('Example 2')
17
Existují další funkce, které jsou vhodné pro použití ve 2D (patch), pˇrípadnˇe ve 3D (mesh, surf). Jejich použití ukážeme, pokud na nˇe narazíme v dalších cviˇceních v prubˇ ˚ ehu semestru.
18
