uvod MatLab skripta-2010

[Mart 2009. godine]

PRAKTIKUM

FIZIKA 2

UVOD U PROGRAMSKI PAKET MATLAB

2009 © Katedra za Mikroelektroniku i tehni čku fiziku © M.Krstić i J. Crnjanski

M.Krstić, J.Crnjanski: Uvod u MatLab

2

MatLab je programski paket koji objedinjuje izračunavanja, vizualizaciju i programiranje neophodne za rešavanje raznovrsnih problema prisutnih u prirodno-matemati čkim i tehničkim naukama. Ime je dobio od re či MATrix LABoratory, pošto mu je osnovni element podataka matrica. MatLab se može koristiti za matemati čka izračunavanja, modelovanje i simulacije, analizu i obradu podataka, grafi čko prikazivanje rezultata i razvoj algoritama. Ova skripta namenjena je studentima koji MatLab koriste prvi put, a u programiranju imaju malo ili nimalo iskustva. Skripta je posebno koncipirana da sadrži sve bitne alatke koje će biti potrebne studentima koji prate kurs Praktikum iz fizike 2. Skripta je podeljena u 10 celina: 1. MatLab okruženje 2. Definicija promenljivih 3. Nizovi i operacije sa nizovima 4. Petlje - For, While, If 5. Grafičko predstavljanje rezultata 6. Funkcije 7. Numerička diferencijacija 8. Numerička integracija 9. Rešavanje sistema linearnih jedna čina 10. Rešavanje diferencijalnih jednačina

2009 © Katedra za Mikroelektroniku u tehičku f iziku, Elektrotehnički f akultet, Beograd


3

1. MatLab okruženje Radno okruženje programskog paketa MatLab sastoji se od tri radna prozora: 1. komandni prozor (Command Window), 2. prethodne komande komande (Command History) History) i 3. tekući folder/radni prostor (Current Directory/Workspace).

TEKUĆI FOLDER/ RADNI PROSTOR

KOMANDNI PROZOR

PRETHODNE KOMANDE

U prozoru prethodne komande pamte se sve komande koje su izvršene iz komandnog prozora, slično kao što i bilo koji web pretraživač pamti History svih pose ćenih sajtova. Dvostrukim klikom mišem na bilo koju komandu iz prethodnih komandi, ta komanda se izvršava kao da je izvršena iz komandnog prozora. U prozoru tekući folder/radni prostor nalazi se pregled foldera koji je podešen kao aktivni folder, a ukoliko se izabere i zabere opcija radni prostor, prikazuje se spisak promenljivih koje su izra čunate izvršavanjem komandi u komandnom prozoru ili pokretanjem nekog programa. Vrednost promenljive se može može ispisati u komandnom komandnom prozoru dvostrukim dvostrukim klikom na željenu željenu promenljivu. Komandni prozor omogućava trenutno izvršavanje pojedina čnih komandi, a istovremeno služi i za prikaz rezultata izvršavanja komandi ili kompletnih programa. Izvršavanje komande i prikaz rezultata njenog njenog izvršenja: >> komanda

Jednom izvršena komanda se ne može naknadno modifikovati vra ćanjem kursora u prethodni red komandnog prozora, ve ć u slučaju greške mora biti ponovo uneta i izvršena. Kursorska strelica “na gore” (↑) na tastaturi omogućava izlistavanje prethodno unetih komandi redom i olakšava ponovno izvršenje već unetih komandi. Brisanje sadržaja komandnog prozora postiže se naredbom: >> clc



4

a upisivanjem tačka-zareza (;) nakon komande, rezultat izvršenja se ne prikazuje u komadnom prozoru. Brisanje sadržaja sadržaja radnog prostora, prostora, odnosno radne radne memorije, postiže se komandom >> clear all

dok komanda >> clear promenljiva

briše samo navedenu navedenu promenljivu. promenljivu. Komentari u kodu se ozna čavaju znakom za procenat (%) i služe kao napomene ili podsetnik za onoga ko čita programski programski kod. MatLab ne izvršava komentare i standardno ih obeležava zelenom bojom.

Osim u komandnom prozoru, komande se mogu pisati i izvršavati i iz MatLab Editora. Da bi se pristupilo Editoru potrebno je startovati File/New/M-File. Otvoriće se novi prozor koji se zove Editor i koji pruža mogućnost ispisivanja programskog koda bez toga da se komande neposredno nakon ispisivanja i izvršavaju. Fajl generisan u Editoru može biti sa čuvan i pokrenut u bilo kom kasnijem trenutku. Ekstenzija fajlova generisanih i snimljenih u Editoru je .m. Snimljeni .m fajl se može pokrenuti ili iz samog Editora naredbom Run, ili se može pozvati iz Komandnog prozora pozivanjem imena fajla kao komande. Pokretanjem .m fajla izvršavaju se sukcesivno sve komande navedene u fajlu. Editor se koristi i za generisanje funkcija tj. pravljenje funkcijskih datoteka o čemu će biti reči nešto kasnije.



5

2. Definicija promenljiv promenljivih ih Za razliku od programskog jezika C, u Matlab-u nije neophodno prethodno deklarisati promenljive i njihov tip, već je potrebno samo uneti vrednosti promenljivih pre njihovog pojavljivanja u naredbama. Promenljiva kojoj je dodeljena numerička vrednost može se upotrebljavati u matematičkim izrazima, funkcijama i svim Matlab-ovim iskazima i komandama. Prilikom definisanje imena promenljive treba imati u vidu da: 1. naziv može da sadrži karaktere: slova ( a-z), cifre (0-9) i/ili donju crtu (_) 2. naziv mora po četi slovom 3. MatLab je „case sensitive“ tj. razlikuje velika od malih slova 4. MatLab 6 pamti prvih 31 karaktera imena promenljive; MatLab 7 pamti prvih 63 karaktera imena promenljive. Primeri dobro definisanih promenljivih su: >> brzina, naelektrisanje_elektrona, a1, A1, promenljiva1_pomocna

.

Primeri loše definisanih promenljivih su: >> naelektrisanje-elektrona, 1a, brzina!, _A1 .

Primeri dodele vrednosti nekoj promenljivoj: >> brzina = 100; >> brzina = 100 + brzina*3;

Pojedine često korišćene promenljive automatski su definisane u MatLab-u. Me đu njima su: >> pi % broj

π.

>> eps % najmanja razlika izme đu dva broja koju MatLab može da uo či iznosi 2 -52 >> i % imaginarna jedinica. Alternativno se koristi karakter j . >> NaN % skraćenica od Not a Number. Upotrebljava se kad MatLab ne može da izra čuna numeri čku vrednost. Na primer, rezultat operacije 0/0 je NaN. >> inf % označava beskona čno veliku vrednost.

Imaginarna jedinica ima svoju rezervisanu oznaku i ili j, međutim to ne sprečava korisnika da ova slova koristi i u druge svrhe. Ukoliko je potrebno da se u istom kodu koristi i ili j i kao imaginarna jedinica i kao proizvoljna promenljiva i ukoliko je i ili j već iskorišćeno kao imaginarna jedinica, tada je poželjno pre koriš ćenja ovih oznaka kao proizvoljnih promenljivih ispisati naredbu >> clear i;

(ili clear j;)

koja oslobađa slovo i ili j za unos proizvoljne vrednosti i time sprečava mogućnost greške ukoliko MatLab promenljivu i ili j shvati kao imaginarnu jedinicu. Zadatak 1: Proveriti ta čnost trigonometrijskog izraza, izračunavanjem obe strane jednačine za vrednost x = π/5: co s 2

2

=

tan x + sin x . 2 tan x

Rešenje: >> clear all; 2009 © Katedra za Mikroelektroniku u tehičku f iziku, Elektrotehnički f akultet, Beograd

M.Krstić, J.Crnjanski: Uvod u MatLab >> x = pi/5; % unos numericke vrednosti za promenljivu x >> leva_strana = cos(x/2)^2; >> desna_strana = (tan(x) + sin(x))/(2*tan(x)); >> display(leva_strana); % display ispisuje sadržaj promenljive >> display(desna_strana);


6


7

3. Nizovi i operacije sa nizovima Ime programskog paketa MatLab, kao što je ve ć rečeno, potiče od MATrix LABoratory i odražava činjenicu da je MatLab prevashodno orijentisan ka radu sa matricama. Čak i kada se definiše prosta skalarna promenljiva, na primer: >> brzina = 100;

MatLab će tu promenljivu tretirati kao matricu dimenzije 1 × 1. Vektor je specijalna vrsta matrice i sadrži samo jednu kolonu (vektor-kolona) ili jednu vrstu (vektor-vrsta). Primeri vektora su: >> [1 2 3 4 5 6 7 8 9] % vektor-vrsta tj. matrica 1 × 9 >> [1; 2; 3; 4; 5; 6; 7; 8; 9] % vektor-kolona tj. matrica 9 × 1

Reč niz podrazumeva bilo vektor bilo matricu. Ukoliko je niz jednodimenzionalan u pitanju je vektor, pošto ima samo jednu kolonu ili vrstu. Ukoliko je niz dvodimenzionalan dvodimenzionalan onda u pravom smislu reči predstavlja matricu koja ima m kolona i n vrsta. MatLab može raditi i sa matricama koje imaju i više od dve dimenzije. 3.1 Vektori

Kao što je ve ć prikazano, jedan od na čina za formiranje vektora vrste je upisivanje njegovih elemenata, razdvojenih razmakom ili zarezom, unutar uglastih zagrada. Sa druge strane, da bi se napravio vektor kolona, potrebno je upisati levu uglastu zagradu [ i zatim elemente razdvojene tačkom i zarezom, ili Enter posle svakog elementa. Nakon poslednjeg elementa, upisuje se desna uglasta zagrada ]. Me đutim, vektore je moguće generisati i na druge na čine. Jedan od načina zahteva definisanje vrednosti po četnog elementa, koraka i vrednosti poslednjeg elementa. Na primer, vektor 2, 4, 6, 8, 10 sadrži brojeve od 2 do 10, sa korakom 2. U MatLab programskom paketu ovakav vektor generiše se naredbom: naredbom: >> vektor = [2:2:10]; % vektor = [pocetni_el [pocetni_element:korak: ement:korak:krajnji_ele krajnji_element] ment]

Korak ne mora imati pozitivnu vrednost. Ukoliko je korak negativan, generiše se vektor čiji elementi imaju opadajuće vrednosti, na primer: 10, 8, 6, 4, 2: >> vektor = [10:-2:2];

Ukoliko prilikom pozivanja prethodnih komandi izostavi vrednost koraka, MatLab podrazumeva da je korak = 1. Tako Takođe, nije neophodno pisati uglaste zagrade, dozvoljeno je: >> vektor = 2:2:10;

odnosno >> vektor = 10:-2:2;

Vektor se može generisati i pomo ću komande linspace koja zahteva definisanje po četnog elementa, krajnjeg elementa i broja elemenata u vektoru, a MatLab sam odre đuje vrednost ekvidistantnog koraka. Na primer, ukoliko je potrebno definisati vektor koji sadrži 5 elemenata, počinje brojem 2, a završava se brojem 10, potrebno je izvršiti komandu: >> linspace(2,10,5);

% linspace(pocetni_element,krajnji_element,broj_elemenata)

Transponovanje Transponovanje vektora postiže se dodavanjem apostrofa nakon oznake vektora: >> x = 1:5; % definisan je vektor vrsta >> x'; % transponovanje vektora-vrste x u vektor-kolonu sa identi čnim elementima



8

3.2 Matrice

Dvodimenzionalna matrica tj. dvodimenzionalan niz sastoji se od elemenata raspore đenih u vrste i kolone. Kao i prilikom definisanja vektora, matrica se može definisati pozivanjem komande: >> matrica = [1 2 3;4 5 6;7 8 9]; % matrica = [elementi_prve_vrste;elementi_druge_vrste;elementi_trece_vrste]

Ovako definisana matrica predstavlja kvadratnu matricu 3 × 3: matrica =

1 2 3 4 5 6 7 8 9

Definisanje elemenata matrice i vektora može se izvesti direktnim unošenjem vrednosti pojedinačnih elemenata, ali i koriš ćenjem matematičkih izraza ili unapred definisanih promenljivih i funkcija. Tako je, na primer, slede ća matrica pravilno generisana ukoliko su prethodno definisane vrednosti promenljivih a, b i c: >> matrica = [a b^2 c^a; log(a) sqrt(9) c; b-a a-c c*a];

Ugrađene komande zeros(m,n) , ones(m,n) i eye(n) generišu matrice čiji elementi imaju specijalne vrednosti. Komande zeros(m,n) i ones(m,n) generišu matrice sa m vrsta i n kolona, u kojima su svi elementi nule odnosno jedinice, respektivno. Komanda eye(n) generiše jediničnu kvadratnu matricu sa n vrsta i n kolona (elementi ove matrice koji se nalaze na glavnoj dijagonali imaju vrednost jedan, dok ostali elementi imaju vrednost nula). Operator transponovanja kod matrica pretvara njene vrste u kolone i obratno. Određenom elementu matrice pristupa se izvršavanjem komande: >> a = matrica(2,1); matrica(2,1); % a = matrica(redni_broj_vrste, redni_broj_kolone) >> display(a);

3.3. Upotreba dvotačke (:)

Upotreba karaktera dvota čke biće prikazana na primerima. Pozivanjem komande: >> x = niz(:);

promenljivoj x dodeljuju se svi elementi promenljive niz. U slučaju: >> x = niz(m:n);

promenljivoj x dodeljuju se elementi promenljive niz, po čev od elementa sa rednim brojem m, do elementa sa rednim brojem n. Promenljiva x sada predstavlja vektor sa m – n + 1 elemenata. Ukoliko su promenljive višedimenzionalne, dvota čka omogućava jednostavno pristupanje određenim elementima matrice: >> A(:,n)

– označava elemente u svim vrstama kolone n matrice A.

>> A(n,:)

– označava elemente u svim kolonama vrste n matrice A.

>> A(:,m:n)

– označava elemente u svim vrstama izme i zmeđu kolona m i n matrice A.

>> A(m:n,:)

– označava elemente u svim kolonama izme đu vrsta m i n matrice A.

– označava elemente u vrstama od m do n i kolonama od p do q matrice A. Neka je matrica A definisana definisana izvršavanjem izvršavanjem komande: >> A(m:n,p:q)

>> A = [1 3 5 7 9 11; 2 4 6 8 10 12; 3 6 9 12 15 18; 4 8 12 16 20 24;... 2009 © Katedra za Mikroelektroniku u tehičku f iziku, Elektrotehnički f akultet, Beograd


9

5 10 15 20 25 30]

Matrice B, C, D, E i F koje sadrže odre đeni broj elemenata matrice A, definisane su komandama: >> B = A(:,3) B = 5 6 9 12 15 >> C = A(2,:) C = 2

4

6

8

10

12

>> E = A(2:4,:) E = 2

4

6

8

10

12

3

6

9

12

15

18

4

8

12

16

20

24

>> F = A(1:3,2:4) F = 3

5

7

4

6

8

6

9

12

3.4 Ugrađene funkcije za obradu nizova:

Slede primeri nekih od najznačajnijih ugrađenih funkcija za obradu nizova: •

length(A) – određuje broj elemenata vektora A: >> A = 1:5; >> length(A)

•

size(A) – određuje dimenzije matrice A: >> A = [1 2 3; 4 5 6]; >> size(A)

•

reshape(A,m,n) – preuređuje matricu A sa r vrsta i s kolona tako da ima m Proizvod r i s mora biti jednak proizvodu m i n:

vrsta i n kolona.

>> A = [5 1 6; 8 0 2] A = 5

1

6

8

0

2



10

>> B = reshape(A,3,2) reshape(A,3,2) B =

•

5

0

8

6

1

2

– kada je v vektor, generiše kvadratnu matricu sa elementima vektora v na glavnoj dijagonali: diag(v)

>> v = [7 4 2]; >> A = diag(v) A =

•

7

0

0

0

4

0

0

0

2

diag(A) – kada je A matrica, generiše vektor od elemenata glavne dijagonale matrice A. >> A = [1 2 3; 4 5 6; 7 8 9] A = 1

2

3

4

5

6

7

8

9

>> v = diag(A) v = 1 5 9

3.5 Matematičke operacije sa nizovima

U programskom paketu MatLab ne postoji razlika izme đu skalarnih i vektorskih vrednosti, tj. i najobičnija skalarna promenljiva se u memoriji pamti kao matrica 1 × 1, pa treba razlikovati operacije nad pojedinačnim elementima i operacije sa celim matricama. Sabiranje i oduzimanje

Operacije sabiranja i oduzimanja nizova mogu biti izvedene nad nizovima istih dimenzija. Drugim rečima, nizovi moraju imati jednak broj vrsta i kolona. Zbir tj. razlika dva niza dobija se sabiranjem tj. oduzimanjem odgovaraju ćih elemenata kao što je to predstavljano jedna činom: L M O am1 L a11

A =

a1n

M , B = amn a11 ± b11 11

A ± B =

M a1m ± b1m

L M O bm1 L b11

L O L

b1n

M bmn

a1n ± b1n

M amn ± bmn



11

Ukoliko se neka skalarna vrednost (matrica 1 × 1) sabira ili oduzima sa nekim nizom, ta skalarna vrednost se sabira odnosno oduzima sa svakim elementom niza. Množenje nizova

Operacija množenja nizova (*) izvodi se u skladu sa pravilima linearne algebre. Kod sabiranja tj. oduzimanja bilo je potrebno da nizovi budu istih dimenzija. Prilikom matrič matri čnog množenja potrebno da broj kolona matrice A bude jednak broju vrsta matrice B. Proizvod ovih matrica imać imaće isti broj kolona kao matrica B i isti broj vrsta kao matrica A. Primera radi, ukoliko je matrica A dimenzija m × k , a matrica B dimenzija k × n, množenje je moguć moguće izvesti, jer matrica A ima k kolona, a matrica B ima k vrsta. Proizvod ove dve matrice bić bi će dimenzija m × n. Sledi primer množenja dve matrice dimenzija 4 × 3 i 3 × 2:

A =

A * B =

a11

a12

a13

a21

a22

a23

a31

a32

a33

a41

a42

a43

b11

b12

, B = b21 b22

( a11b11 + a12 b2211 + a13 b3311 ) ( a21b11 + a22 b2211 + a23 b3311 ) ( a31b11 + a32 b2211 + a33b3311 ) ( a41b11 + a42 b2211 + a43 b3311 )

b31

b32

( a11 b12 + a12 b2222 + a13 b3232 ) ( a21 b12 + a22 b2222 + a23 b3322 ) ( a31 b12 + a32 b2222 + a33 b3232 ) ( a41 b12 + a42 b2222 + a43 b3322 )

Treba primetiti da operacija množenja matrica nije komutativna, tj. važi nejednakost A* B ≠ B* A. Takođ Takođe, treba primetiti da je stepenovanje matrice moguć moguće izvesti samo ako je matrica kvadratna, tj. dimenzije n × n, pošto se proizvod A* A može izrač izračunati samo ako je broj kolona prve matrice jednak broju vrsta druge druge matrice. Dva vektora mogu da se pomnože samo ako imaju isti broj elemenata i ako je jedan vektorvrsta, a drugi vektor-kolona. Množenjem vektora-vrste i vektora-kolone dobija se matrica 1 × 1 tj. skalarna vrednost. Ovakvo množenje naziva skalarni proizvod a ugra đena funkcija koju MatLab koristi za izrač izračunavanje skalarnog proizvoda dva vektora je dot(a,b) . Ulazni argumenti a i b mogu biti ili vektori-vrste ili vektori-kolone: >> A =[1 2 3] A = 1

2

3

>> B = [1 2 3] B = 1

2

3

>> dot(A,B) ans = 14 >> A = [1 2 3] A = 1

2

3



12

>> B = [1;2;3] B = 1 2 3 >> dot(A,B) ans = 14

Množenjem vektora kolone i vektora vrste sa po n elemenata dobija se kvadratna matrica n × n: >> A = [1;2;3] A = 1 2 3 >> B = [1 2 3] B = 1

2

3

1

2

3

2

4

6

3

6

9

>> A*B ans =

Deljenje nizova

Kao i kod množenja, tako i kod deljenja važe pravila linearne algebre. U programskom paketu MatLab postoje dve vrste deljenja: deljenje sa leva ( \) i deljenje sa desna ( /). Međ Međutim, pre objašnjenja procedure deljenja nizova, bić bi će uvedeni pojmovi jedinič jedini čne matrice, inverzne operacije i determinante. je kvadratna matrica koja ima jedinice na glavnoj dijagonali i nule na svim ostalim mestima. Ona se može napraviti komandom eye(), kao što je već već pomenuto. Kada se neka matrica ili vektor pomnoži jedinč jedinčnom matricom, dobija se ta ista matrica ili vektor. Specijalno, ako je matrica kvadratna, onda se ona može pomnožiti jedini čnom matricom i sa leva i sa desna, tj. u tom sluč slučaju operacija množenja je komutativna. Matrica B je inverzna matrica matrici A ako je proizvod te dve matrice jedinič jedini čna matrica. U MatLabu inverzna matrica dobija se stepenovanjem na -1 ili funkcijom inv(A) . Inverzna matrica postoji samo ako je matrica kvadratna i ako ako je njena determinanta determinanta različ različita od nule. Jedinična matrica



13

je funkcija pridružena kvadratnim matricama tj. to je funkcija koja svakoj kvadratnoj matrici pridružuje broj koji se zove detereminanta matrice. Bez ulaženja u strogu definiciju determinante, dat je primer određ određivanja determinante matrice dimenzije 2 × 2: Determinanta

A =

a11

a12

a21

a22

det A =

a11

a12

a21

a22

= a11a22 − a12 a21

.

Ugrađ Ugrađena funkcija koja određ odre đuje determinantu matrice A poziva se komandom: >> det(A) . Deljenje s leva (\)

Deljenjem sa leva može se odrediti rešenje matrič matri čne jednač jednačine AX = B, gde su X i B vektori-kolone. Matematič Matematički gledano, ova jednač jedna čina se rešava množenjem obe strane matricom inverznom matrici A: 1

A− AX

=

1

A− B

IX

ï

=

1

A− B

X

ï

=

1

A− B

Poslednja jednač jednačina definiše rešenje X , a MatLabu se realizuje pomoć pomoću operacije deljenja sa leva: >> X = A \ B

Rešenje je moguć mogu će dobiti i množenjem inverznom matricom. Razlika je u algoritmu koji MatLab koristi u ova dva sluč slučaja. Prilikom množenja inverznom matricom, MatLab prvo traži inverznu matricu, a zatim je množi sa matricom B. U sluč slučaju deljenja s leva, rešenje se dobija numerič numeri čki, primenom Gauss-ove eliminacije. Za rešavanje sistema linearnih jednač jednačina preporuč preporučuje se matrič matrično deljenje s leva, jer je to u sluč slučaju velikih matrica efikasniji metod. Deljenje s desna (/)

Deljenjem sa desna može se rešiti matrič matri čna jednač jednačina XC = D, gde su X i D vektori-vrste. Matematič Matematički, jednač jednačina se može rešiti množenjem obe strane jednač jedna čine matricom koja je inverzna matrici C : XCC −1

=

DC −1

ï

X

=

DC −1

Poslednja jednač jednačina definiše rešenje X , a u MatLabu se realizuje pomoć pomo ću operacije deljenja sa desna: >> X = D / C

Na primerima koji slede bić biće pokazano rešavanje sistema linearnih jednač jedna čina deljenjem sa leva i deljenjem sa desna. Zadatak 2: Primenom matrič matričnih operacija, rešiti sistem jednač jedna čina: 4 x − 2 y + 6 z = 8 2 x + 8y + 2z = 4 6 x + 10 y + 3z = 0. Rešenje: Sistem se može predstaviti u matrič matri čnom obliku i to na dva nač načina, predstavljanjem nepoznatih preko vektora-kolone i predstavljanjem nepoznatih preko vektora-vrste: 4 −2 6 x 2 8 2 y 6 10 3 z

=

8 4 0

ñ

AX

=

B,



4 2 6 x y z −2 8 1 100 6 2 3

14

=

8 4 0

ñ

XC = D.

U prvom sluč slučaju, kada su nepoznate predstavljene preko vektora-kolone, rešenje se dobija pozivanjem komandi: >> A = [4 -2 6; 2 8 2; 6 10 3]; >> B = [8;4;0]; >> X1 = A\B X1 = -1.8049 0.2927 2.6341

Primenom metode množenja sa inverznom matricom, dobija se >> X1prim = inv(A)*B X1prim = -1.8049 0.2927 2.6341

U drugom sluč slučaju, kada su nepoznate predstavljene preko vektora-vrste, rešenje se dobija pozivanjem komandi: >> C = [4 2 6; -2 8 10; 6 2 3]; >> D = [8 4 0]; >> X2 = D/C X2 = -1.8049

0.2927

2.6341

odnosno primenom metode množenja sa inverznom matricom: >> X2prim = D*inv(C) X2prim = -1.8049

0.2927

2.6341

Operacije sa pojedinačnim elementima niza

Ukoliko se ispred operatora množenja (*), deljenja (/) ili stepenovanja (^) stavi tač ta čka (.), MatLab izvršava odgovarajuć odgovarajuću operaciju množenja, deljenja ili stepenovanja element po element, a ne matrič matrično. Primera radi, ako vektor T predstavlja vremensku osu, tj. sadrži vremenske trenutke u kojima se odvijao neki događ događaj, a u nekom matematič matemati čkom izrazu figuriše T 2, potrebno je izvršiti kvadriranje element po element, što se postiže komandom: >> T.^2;

Sledeć Sledeći primer prikazuje matrič matrično kvadridanje kvadratne matrice i kvadriranje element-po-element: >> M = [1 2 3;4 5 6;7 8 9] M = 1

2

3


M.Krstić, J.Crnjanski: Uvod u MatLab 4

5

6

7

8

9

>> M^2 % matricno kvadriranje ans = 30

36

42

66

81

96

102

126

150

>> M.^2 % kvadriranje element-po-element ans = 1

4

9

16

25

36

49

64

81


15


16

4. Petlje - For, While, If Osnovne petlje u programskom paketu MatLab su For, While i If petlje. For-petlja obavlja zadate naredbe više puta u skladu sa broja čem koji je definisan od strane korisnika: >> for brojac >>

telo petlje;

>> end;

Primera radi, formiranje vektora-vrste od 10 elemenata (1, 2, 3, 4, 5, 6, 7, 8, 9, 10) može se posti ći primenom For-petlje: >> for i = 1:1:10 >>

vektor(i) = i;

>> end;

While-petlja obavlja zadate naredbe u okviru tela petlje, sve dok zadati uslov daje logi čki tač tačnu vrednost: >> while uslov >>

telo petlje;

>> end;

Formiranje vektora-vrste iz prethodnog primera može se ostvariti i primenom While-petlje, s tom razlikom što u okviru While-petlje nije automatski definisan brojač broja č, niti postoji moguć mogućnost njegovog automatskog inkrementiranja, već ve ć je ove korake neophodno dodati u kod: >> i=1; >> while i<=10 >>

vektor(i) = i;

>>

i = i+1;

>> end;

If-petlja ima nešto drugač drugačiju namenu i na osnovu ispitivanja da li je zadati uslov ispunjen ili ne, usmerava izvršenje koda u jednom (ili nijednom) od više ponuđ ponu đenih pravaca: >> if uslov >>

ako je uslov ispunjen uraditi ovo;

>> else >>

ako uslov nije ispunjen uraditi ovo;

>> end;

Primena If-petlje prikazana je na sledeć slede ćem primeru. Potrebno prebrojati koliko elemenata niza iz prethodnog primera ima vrednost već veću od 5, a koliko manju ili jednaku 5: >> manji_od5 = 0; >> veci_od5 = 0; >> for i = 1:1:10 >>

if vektor(i)>5

>> >> >>

veci_od5 = veci_od5+1; else manji_od5 = manji_od5 + 1; 2009 © Katedra za Mikroelektroniku u tehičku f iziku, Elektrotehnički f akultet, Beograd

M.Krstić, J.Crnjanski: Uvod u MatLab >>

end

>> end >> manji_od5 manji_od5 = 5 >> veci_od5 veci_od5 = 5


17


18

5. Grafičko predstavljanje rezultata U ovom poglavlju bić biće objašnjeno kako se u programskom paketu MatLab generišu dvodimenzionalni i trodimenzionalni grafikoni. 5.1 Dvodimenzionalni (2D) grafici

Komanda plot omoguć omogućava generisanje 2D grafikona u pravouglom koordinatnom sistemu. U najkrać najkraćem obliku, komanda glasi plot(x,y) , gde su x i y vektori tj. jednodimenzionalni (1D) nizovi sa istim brojem elemenata (istih dimenzija). Pozivanjem ove komande u posebnom prozoru označ označenom kao Figure prikazuje se grafik funkcije y = y( x x), primera radi: >> x = [1 2 3 4 5 6 7 8 9 10]; >> y = [2 4 6 8 10 12 14 16 18 20]; >> plot(x,y);

Na dobijenom grafiku vrednosti vektora x su na apcisi, a vrednosti y su na ordinati. Automatska boja linije je plava. Opšti oblik komande plot koji omoguć omogućava modifikacije izgleda grafika dat je u formi: plot(x,y,'oznaka plot(x,y,'oznaka linije','ime linije','ime svojstva',vrednost svojstva',vrednost svojstva)

Promena oznake linije i svojstva su opcionalne. Oznaka linije definiše vrstu i boju linije i markera. Markeri su oznake tač tačaka u kojima je definisana vrednost funkcije. Za gornji primer, markeri bi stajali u tač tačkama (1, 2), (2, 4), (3, 6)... U sledeć sledećim tabelama date su vrste linija i njihove oznake, boje linija i njihove oznake oznake i vrste markera i njihove oznake: Vrsta linije

puna (automatska) (automatska) isprekidana tač tačkasta crta-tač crta-tačka

Oznaka

-: -.

Boja linije

crvena zelena plava cijan magenta žuta crna bela

Oznaka r g b c m y k w

Vrsta markera

znak plus kružić kružić zvezdica tač tačka kvadrat romb petokraka zvezda zvezda šestokraka zvezda


Oznaka + ° * . s d p h


19

Ime svojstva i vrednost svojstva omoguć omogu ćavaju definisanje debljine linija, velič veli čine markera, boja ivica markera i boje boje kojom se markeri markeri ispunjavaju. U tabeli tabeli je dat pregled svojstava: svojstava: Ime svojstva LineWidth MarkerSize MarkerEdgeColor MarkerFaceColor

Opis

debljina linije velič veličina markera boja ivice markera popuna markera

Moguće vrednosti svojstva

broj (automatski 0.5) broj oznaka boje – slovo oznaka boje – slovo

Primer korišć korišćenja komande plot u punom obliku, ukoliko se za liniju grafika debljine 2 koristi crvena boja, a markeri su velič veli čine 12 i u obliku kružić kruži ća, ispunjeni žutom i oivič oivičeni zelenom bojom: >> x = [1:1:10]; >> plot(x,sin(x),'-ro','LineWidth',2,... 'MarkerSize',12,'MarkerEdgeColor','g','MarkerFaceColor','y')

Naravno, kako je uzet jako veliki korak (svega deset tač ta čaka), grafič grafički prikaz funkcije sin(x) je veoma neprecizan. Ukoliko postoji potreba da se istovremeno prikaže više grafikona, pre svake naredbe plot, potrebno je otvoriti zaseban prozor za iscrtavanje grafika naredbom figure(n), gde je n redni broj prozora: prozora: >> x = [-2:0.01:4]; >> y = 3.5.^(-0.5*x).*cos(6*x); >> figure(1) >> plot(x,y,':g','LineWidth',2); >> y2 = x.^2+4*sin(2*x)-1; >> figure(2) >> plot(x,y2,'--k','LineWidth',2);



20

Sa druge strane, ukoliko je potrebno da se u istom koordinatnom sistemu grafi čki prikažu dve zavisnosti, umesto otvaranja novog prozora komandom figure(2) koristi se komanda hold all , koja na postojeć postojeću sadržinu Figure prozora dodaje novu grafič grafi čku zavisnost: >> x = [-2:0.01:4]; >> y = 3.5.^(-0.5*x).*cos(6*x); >> figure(1) >> plot(x,y,':g','LineWidth',2); >> hold all; >> y2 = x.^2+4*sin(2*x)-1; >> plot(x,y2,'--k','LineWidth',2);

Konačno, postoji moguć Konač mogućnost da se u okviru istog Figure prozora prikaže više grafika sa nezavisnim koordinatnim osama upotrebom komande subplot. Sintaksa komande >> subplot(m,n,p)

deli prostor Figure prozora u tabelu koja se sastoji od m redova i n kolona. Brojem p definiše se odgovarajuć odgovarajuća pozicija na koju je potrebno da se smesti grafik. Slede primeri koriš ćenja komande subplot:



21

5.2 Formatiranje grafika

Formatiranje grafika podrazumeva opciono dodavanje i formatiranje naslova grafika, imena koordinatnih osa i njihovog opsega, legende i mreže: • • • • • • •

– dodavanje i formatiranje naslova xlabel('tekst') – dodavanje oznake ose na na apcisu ylabel('tekst') – dodavanje oznake oznake ose na na ordinatu legend('tekst1','tekst2',...,pol)– dodavanje legende xlim([x1 x2]) – definisanje opsega vrednosti za apcisu ylim([y1 y2]) – definisanje opsega vrednosti za ordinatu grid on ili grid off – prikazuje pomoć pomoćnu mrežu u skladu sa oznakama na osama title('zeljeni naslov grfikona')

Primer formatiranja grafika: >> x = [-2:0.01:4]; >> y = 3.5.^(-0.5*x).*cos(6*x); >> figure(1) >> plot(x,y,':g','LineWidth',2); >> hold all; >> y2 = x.^2+4*sin(2*x)-1; >> plot(x,y2,'--k','LineWidth',2); >> title('Uporedni prikaz grafika'); >> xlabel('x osa'); >> ylabel('y osa'); 2009 © Katedra za Mikroelektroniku u tehičku f iziku, Elektrotehnički f akultet, Beograd


22

>> legend('y = 3.5 ^-^0^.^5^x cos(6x)','y cos(6x)','y = x^2 + 4sin(2x) - 1',2); % MatLab automatski dodeljuje prvi tekst prvoj funkciji, a drugi tekst drugoj % funkciji >> xlim([-1 3.3]); >> ylim([-5 15]);

Pol

Položaj

-1

Izvan granica osa, desna strana Unutar granica osa, gde najmanje smeta Gornji desni ugao Gornji levi ugao Donji levi ugao Donji desni ugao

0 1 2 3 4

5.3 Logaritamske podele

U mnogim nauč naučnim i tehnič tehničkim primenama, ponekad je potrebno da jedna ili obe ose grafika imaju logaritamsku podelu: >> semilogy(x,y) % logaritamska podela (osnova 10) y ose i linearna podela x ose >> semilogx(x,y) % logaritamska podela (osnova 10) x ose i linearna podela y ose >> loglog(x,y) % logaritamska podela i x i y ose

Sledi primer grafič grafičkog prikazivanja funkcije podelu osa:

y = 2exp(-0.2x + 10)

za linearnu i logaritamsku

>> x = linspace(0.1,60,1000); >> y = 2.^(-0.2*x + 10); >> subplot(2,2,1); >> plot(x,y,'b','LineWidth',2); >> xlabel('Linearna skala'); >> ylabel('Linearna skala'); >> grid on; >> subplot(2,2,2); 2009 © Katedra za Mikroelektroniku u tehičku f iziku, Elektrotehnički f akultet, Beograd


23

>> semilogy(x,y,'b','LineWidth',2); >> xlabel('Linearna skala'); >> ylabel('Log skala'); >> grid on; >> subplot(2,2,3); >> semilogx(x,y,'b','LineWidth',2); >> xlabel('Log skala'); >> ylabel('Linearna skala'); >> grid on; >> subplot(2,2,4); >> loglog(x,y,'b','LineWidth',2); >> xlabel('Log skala'); >> ylabel('Log skala'); >> grid on;

5.4 Grafici u polarnim koordinatama

MatLab omoguć omogućava jednostavan prikaz zavisnosti u polarnim koordinatama primenom komande: polar(ugao,poteg,'oznake linije')

Sledi primer korišć korišćenja ove komande: >> ugao = linspace(0,2*pi,200); >> poteg = 3*cos(0.5*ugao).^2 + ugao; >> figure(1) >> polar(ugao,poteg,'b');



24

Pored navedenih moguć mogu ćnosti, programski paket MatLab omoguć omogu ćava generisanje i nekih specijalnih vrsta grafika: vertikalni trakasti grafikon, horizontalni trakasti grafikon, stepenasti grafikon, grafikon diskretnih podataka, kružni dijagram, histogram... Detaljne informacije o ovim tipovima grafika se mogu nać naći u MatLab Help-u. 5.5 Trodimenzionalni (3D) grafici

Generisanje 3D grafikona postiže se komandama mesh i surf. Komanda mesh generiše trodimenzionalnu zavisnost u formi mreže, dok naredba surf prikazuje trodimenzionalnu zavisnost kao ispunjenu površinu. Razlika izmeđ izme đu ove dve komande se najlakše može uo čiti na konkretnom primeru: >> [x,y] = meshgrid(-8:0.5:8); meshgrid(-8:0.5:8); >> r = sqrt(x.^2+y.^2)+eps; >> z = sin(r)./r; >> figure(1);mesh(z); >> figure(2);surf(z);

Naredba meshgrid se koristi za generisanje podele x i y ose u formi matrica x i y odgovarajuć odgovarajućih diemenzija. Promenljiva z definisana je kao sinc funkcija (podseć (podseća na meksič meksički šešir – sombrero). >> mesh(x,y,z) >> surf(x,y,z)



25

6. Funkcije Programski paket MatLab raspolaže velikim brojem ugrađ ugra đenih funkcija kao što su exp(x) , sqrt(x), sin(x)... Međ Međutim, MatLab omoguć omogućava i da korisnik samostalno definiše funkciju, generisanjem .m fajla u MatLab Editoru i snimanjem ovog fajla u obliku funkcijske datoteke. Komandna linija koja definiše funkciju ima oblik: function [izlaz1, izlaz2, ...] = imefunkcije(ulaz1,ulaz2, ...)

Promenljive izlaz1, izlaz2 ,... su promenljive u koje se smešta rezultat funkcije, dok su ulazne promenljive vrednosti za koje funkcija treba da obavi traženu operaciju. Imefunckije je proizvoljno korisnič korisničko ime pod kojim funkcija mora da se zapamti. Kasnije, pod tim imenom funkcija se poziva iz glavnog programa ili Komandnog prozora. Obratiti pažnju da posle definisanja funkcije ne stavlja ; !!! Zarad preglednosti, funkcije će biti objašnjene na konkretnom primeru. Zadatak 3: Napraviti funkciju koja rešava kvadratnu jednač jednačinu oblika: ax

2

+ bx + c =

0,

gde su a, b, c koeficijenti, a x promenljiva. Rešenje kvadratne jednač jednačine dato je izrazom: x1/ 2

=

−b ±

b

2

2a

−

4ac

.

Sledi kod funkcije koja rešava zadati problem, a zatim i glavni program iz koga se poziva generisana generisana funkcija za određ određene koeficijente a, b, c. Funkcija kvadratna_jednacina: kvadratna_jednacina: function [x1 x2] = kvadratna_jednacina(a,b,c) %izlazne promenljive su x1 i x2. Ime funckije je kvadratna_jednacina i pod tim imenom se snima funkcija kao funkcijska datoteka tj. kao .m fajl. Isto ime se kasnije koristi za pozivanje funkcije iz glavnog programa. Ulazne promenljive su a, b, c koje predstavljaju koeficijente koeficijent e jednacine x1 = (-b + sqrt(b^2 - 4*a*c))/(2*a); x2 = (-b - sqrt(b^2 - 4*a*c))/(2*a);

Nakon generisanja funkcije, potrebno pozivanjem File/Save As snimiti funkciju pod imenom kvadratna_jednacina. kvadratna_jednacina. Glavni program: clear all; close all a = 5;b = 8;c = 12; [resenje1 resenje2] = kvadratna_jednacina(a,b,c); % Promenljiva resenje1 prihvata vrednost iz promenljive x1, a promenljiva resenje2 vrednost iz promenljive x2. Nakon pokretanja glavnog programa u Komandnom prozoru prikazuju se resenja: >>resenje1 resenje1 = -0.8000 + 1.3266i >> resenje2 resenje2 = -0.8000 - 1.3266i 2009 © Katedra za Mikroelektroniku u tehičku f iziku, Elektrotehnički f akultet, Beograd


26

7. Numerička diferencijaci diferencijacija ja Prema definiciji, izvod funkcije f određ određen je izrazom: df dx

=

lim h →0

f ( x + h ) − f ( x ) h

.

Ukoliko je x vektor: x = [x(1) x(2) ... x(n)]

onda komanda diff(x) daje vektor razlika dva susedna elementa vektora x: diff(x) = [x(2)-x(1) x(3)-x(2) ... x(n)-n(n-1)]. x(n)-n(n-1)].

Treba napomenuti da dobijeni vektor ima jedan element manje nego vektor x. Neka su dati vektori x i y koji predstavljaju argumente i vrednosti funkcije y(x), respektivno. Izvod dy/dx se može odrediti pozivanjem komande: >> diff(y)./diff(x);

Komanda diff(y,n)./diff(x,n) generiše n-ti izvod funkcije y. Zadatak 4: Ispitati funkciju 3x3 – x2 + 6, čiji je izvod po x jednak 9x2 – 2x. Nacrtati u istom prozoru izvod kao grafik funkcije 9x2 – 2x i izvod posmatrane funkcije koristeć koristeći komandu diff i time pokazati da su ova dva grafika ista. x = linspace(0,50,5000); linspace(0,50,5000); y = 3*x.^3 - x.^2 + 6; figure(1);subplot(2,1,1); plot(x,9*x.^2 - x,'LineWidth',2); % Prethodna komanda iscrtava funkciju y(x). Da bi se grafi čki prikazao i vektor diff(y)./diff(x) potrebno je definisati x osu sa jednim elementom manje: for i=1:length(x)-1 xosa(i) = x(i); end subplot(2,1,2); plot(xosa,diff(y)./diff(x),'LineWidth',2);



27

8. Numerička integracija Određ Određeni integral funkcije f(x) od a do b, predstavlja površinu koju kriva zahvata sa apcisom (x osom) u granicama [a,b]. U numerič numeričkoj matematici postoji više algoritama tj. metoda za izrač izračunavanje određ određenog integrala. MatLab nudi tri ugrađ ugra đene funkcije za izrač izračunavanje određ određenog integrala: quad, quadl, trapz. Komanda quad koristi Simpson-ovu metodu integraljenja i ima sledeć sledeći oblik: quad('funkcija',a,b)

ili quad(@funkcija,a,b) .

Komanda quadl koristi Lobatto-ovu metodu integraljenja. Ima identič identičan oblik kao i prethodna funkcija: quadl('funkcija',a,b)

ili quadl(@funkcija,a,b)

Komanda trapz se koristi za integraljenje funkcije zadate u obliku skupa tač ta čaka. Koristi metodu trapeza. Oblik komande je: trapz(x,y)

gde su x i y vektori x koordinata i y koordinata. Vektori moraju imati isti broj elementa. Zadatak 5: Primenom ugrađ ugra đenih funkcija odrediti numerič numeričku vrednost integrala: 8

0 ( x exp ( −3x + 10) + 2x + 10) dx . >> quad('x.*exp(-3*x+10) + 2*x + 10',0,8) ans = 2.5914e+003 >> quadl('x.*exp(-3*x+10) + 2*x + 10',0,8) ans = 2.5914e+003

Postoji još jedan nač način na koji se mogu izvršiti prethodne komande i on podrazumeva prethodno definisanje funkcijske datoteke koja sadrži definiciju funkcije f unkcije f = x*exp(3*x+10) + 2*x + 10 , a zatim pozivanje komande u quad obliku: >> quad(@f,0,8)



28

9. Rešavanje sistema linearnih jedna čina Sistem jednač jednačina definisan izrazima: a11 x1 + a12 x2

+ a13 x3 = b1

a21 x1 + a22 x2

+ a23 x3 = b2

a31 x1 + a32 x2

+ a33 x3 = b3

Matrič Matrični zapis ovog sistema glasi: A ⋅ X

a11

a12

a13

A = a21

a22

a23 ,

a31

a32

a33

=

B,

b1 B = b2 , b3

1

X

=

x2 , x3

pa je moguć moguće primenom operacije deljenja sa leva (\) rešiti jednač jedna činu AX = B: >> X = A\B;

što je ekvivalentno matematič matemati čkom rešenju X = A-1 B. Dobijeni vektor X je vektor kolona koji u prvom redu ima rešenje rešenje po x1, u drugom po x2, a u treć trećem po x3.



29

10. Rešavanje diferencijalnih jedna čina Diferencijalne jednač jednačine igraju jako važnu ulogu u nauci i tehnici jer predstavljaju osnovu gotovo svakog fizič fizičkog procesa. I ne samo fizič fizi čki procesi, već već i razni ekonomski, biološki, sociološki, klimatološki i drugi procesi, modeluju se određ odre đenim diferencijalnim jednač jednačinama ili sistemima diferencijalnih jednač jednačina. Problem je, međ međutim, što se jako mali broj tih jednač jednačina može rešiti analitič analitički. Programski paket MatLab pruža širok opseg funkcija tj. "solvera" koji služe za numerič numeričko rešavanje složenih problema. U narednom razmatranju bić bi će dat prikaz rešavanja obič obi čnih diferencijalnih jednač jednačina (Ordinary Differential Equations, ODE) sa poč početnim uslovima. Obič Obična diferencijalna jednač jednačina data je u obliku: dx dt

= f

( x, t ) ,

gde je f ( x x,t ) poznata funkcija od x i t . Rešavanje ovakve diferencijalne jednač jedna čine prvog reda (jednač (jednačina je prvog reda pošto u njoj figuriše prvi izvod promenljive x) zahteva poznavanje jednog poč početnog uslova, odnosno vrednost promenljive x0 u nekom trenutku t = t 0. Diferencijalne jednač jednačine se mogu podeliti na linearne i nelinearne jednač jedna čine. Linearna diferencijalna jednač jednačina je ona jednač jedna čina u kojoj nepoznata, odnosno zavisna promenljiva x i njeni izvodi figurišu u linearnom smislu. Prema tome, jednač jednačina: dx dt

= c1 x + c2

u kojoj su c1 i c2 konstante, je linearna, jer su x i dx/dt prvog stepena. Štaviše, prethodna jednač jedna čina bila bi linearna čak i onda kada bi c1 i c2 bili komplikovane funkcije vremena t , jer je od znač zna čaja samo na koji nač način u jednač jednačini figurišu zavisna promenljiva x i njeni izvodi. Primeri linearnih diferencijalnih jednač jednačina su: dx

=

dt dx

(1 − x ) x.

=

dt x

dx

=

dt

sin x, x + t ...

Nelinearne diferencijalne jednač jednačine se, generalno gledano, mogu podeliti na dve klase: autonomne i neautonomne jednač jedna čine. Ukoliko se ovom trenutku razmatranje ogranič ograniči na nelinearne diferencijalne jednač jednačine prvog reda, autonomne jednač jednačine su one kod kojih je brzine promene x prosto funkcija samog samog rešenja x i ne zavisi eksplicitno od vremena: dx dt

= f

( x) .

Kod neautonomnih jednač jednačina, brzina promene x zavisi i od vremena t . Analitič Analitička rešenja su moguć moguća samo za specijalne sluč slu čajeve. 10.1 Postupak rešavanja diferencijalne jednačine

Prilikom rešavanja diferencijalnih jednač jednačina poželjno je pratiti proceduru koja će biti izložena korak po korak. 2009 © Katedra za Mikroelektroniku u tehičku f iziku, Elektrotehnički f akultet, Beograd


30

Korak1 - Postavka problema

Za konkretan fizič fizički problem generisati matematič matematičko-fizič ko-fizički model, odnosno formulisati odgovarajuć odgovarajuće diferencijalne (i ukoliko je potrebno, alegabarske) jednač jednačine. Korak 2 - Formiranje funkcijske datoteke Svaku diferencijalnu jednač jednačinu predstaviti u skladu sa sintaksom MatLab-a u odgovaraju ćoj funkcijskoj datoteci, što omoguć omogu ćava jednostavno pozivanje jednač jedna čina iz bilo kog drugog MatLab programa, bilo iz Editora ili komandnog prozora. Slič Slično kao i kod formiranja bilo koje druge funkcije, funkcijska datoteka se definiše zaglavljem koje ima oblik: function lokalna_promenljiva = ime_funkcije(promenljive)

gde je

lokalna_promenljiva lokalna promenljiva u koju se smešta rezultat izvršenja funkcije, a ime_funkcije je proizvoljno korisnič korisničko ime po kome će se funkcija pozivati iz glavnog programa.

Formiranje jednostavne funkcijske datoteke bić bi će prikazano na primeru autonomne nelinearne jednač jednačine prvog reda kojom se mogu modelovati epidemski procesi. Zadatak 6: Jedna kolonija od 1000 bakterija razmnožava se brzinom od r = 0.8 jedinki na sat. Koliko će bakterija biti u koloniji nakon 10 sati? Rešenje: Proces razmnožavanja bakterija se može modelovati diferencijalnom jednač jednačinom oblika (korak 1): dN dt

=

rN ,

gde je N broj individua tj. populacija bakterija, a poč po četni uslov je definisan sa N (0) (0) = 1000. Ova jednač jednačina se može rešiti analitič analitički i ima dobro poznato rešenje u obliku ekponencijalne funkcije: N ( t ) = N ( 0 ) exp ( rt ) .

Rešenje ove jednač jednačine se može dobiti i numerič numeri čki, koristeć koristeći integrisan MatLab ODE solver. Nakon postavke problema, pristupa se pisanju funkcijske datoteke (korak 2) u MatLab editoru: function dNdt = bakterije(t,N) % zaglavlje funkcijske datoteke % u lokalnu promenljivu dNdt smesta se rezltat izvrsenja funkcije; % ime funkcije je bakterije i pod tim imenom ce funkcija biti pozivana iz nekog drugog programa ili direktno iz komandnog prozora % promenljive su t i N - vreme i populacija bakterija. % Obratiti paznju da se ne kraju komandne linije u kojoj se definise funkcija NE stavlja ta čka-zarez (;) ! dNdt = 0.8*N; % telo funkcije

Ovako napisanu funkcijsku datoteku potrebno je snimiti kao .m fajl pod imenom koje odgovara imenu funkcije, u konkretnom sluč slu čaju bakterije.m. Osim toga, funkciju je uvek potrebno snimiti u isti folder gde se nalazi i glavni program iz kog će ta funkcija biti pozvana. Korak 3 - Odabir solvera

MatLab nudi više ugrađ ugrađenih solvera za rešavanje diferencijalnih jednač jedna čina. Svaki od njih koristi neki drugi algoritam tj. neku drugu numerič numeri čku proceduru i pogodan je za određ odre đeni tip 2009 © Katedra za Mikroelektroniku u tehičku f iziku, Elektrotehnički f akultet, Beograd


31

jednač jednačina. U sledeć sledećoj tabeli dat je pregled ODE solvera i date su napomene kada se dati solver primenjuje: Solver

Tip problema

Preciznost

Kada koristiti

ode45

Linearan

Srednja

Prvo probati ovaj solver (Runge-Kuta metoda - single step metoda)

ode23

Linearan

Niska

Kada se ne traži velika preciznost (Euler-ova metoda - single step)

ode113

Linearan

Visoka

Multistep metoda – kada se traži visoka preciznost

ode15s

Nelinearan

Srednja

Ako ode45 ne radi

ode23s

Nelinearan

Niska

Kada se ne traži velika preciznost

Nis Niska


Niska


ode23t Osre Osred dnje nje neli nelin neara earan n ode23tb

Nelinearan

Korak 4 - Pokretanje ODE solvera iz

glavnog programa i podešavanje parametara ODE solvera

Komanda koja pokreć pokreće izabrani solver i poziva funkciju definisanu u okviru funkcijske datoteke ima oblik: [t,x] = odgovarajuci_solver('ime funkcije',[interval_t],x0);

Odnosno, za sluč slučaj konkretnog primera iz zadatka 6: >> [t,N] = ode23('bakterije',[0 10],1000); % promenljiva t prihvata vremensku osu, promenljiva N prihvata odgovaraju vremensku zavisnost rešenja, odnosno broja bakterija. % upotrebljen ‘bakterije’

je

solver

ode23,

koji

poziva

prethodno

definisanu

ću

funkciju

% vremenski interval u toku kog se posmatra dinamika sistema je definisan sa [0 10] s obzirom na to da se u zadatku trazi populacija bakterija nakon 10 sati, a da je zadat pocetni uslov da je populacija u trenutku t = 0 bila 1000 jedinki

Ceo programski kod glavnog programa koji pokre pokrećće ODE solver i prikazuje rešenje jednač jedna čine izgleda ovako: >> clear all, close all % definisanje pocetnog i krajnjeg trenutka vremena >> tpoc = 0;tkraj = 10; % definisanje vremenske ose sa rezolucijom od 100 tacaka >> vreme = linspace(tpoc,tkraj,100); >> N0 = 1000; % pocetni uslov N(0)= 1000 % pokretanje ode solvera. >> [t,N]=ode23('bakterije',[tpoc tkraj],N0); % Upotrebljen je ode23, zato sto je ovo problem koji numericki nije zahtevan (sto potvrdjuje cinjenica da za ovaj problem postoji analiticko resenje). Solver ode23 nema veliku preciznost ali daje veoma dobre rezultate za resavanje jednacina koje cije resenje je eksponencijalna funkcija. Osim toga ode23 nije vremenski zahtevan solver pa se resenje dobija relativno brzo. % prikaz rezultata 2009 © Katedra za Mikroelektroniku u tehičku f iziku, Elektrotehnički f akultet, Beograd


32

>> figure(1); plot(t,N,'LineWidth',2); >> title('Populacija bakterija nakon 10h - Numericko resenje'); >> xlabel('Vreme'); >> ylabel('Populacija'); % analiticko resenje jednacine >> Nanaliticko = N0*exp(0.8*vreme); >> figure(2);plot(vreme,Nanaliticko,'LineWidth',2); >> xlabel('Vreme'); >> ylabel('Populacija'); >> title('Analiticko resenje'); 6

3

6

Populacija bakterija nakon 10h

x 10

3

2.5

2.5

2

2

a j i c a l u 1.5 p o P

a j i c a l u 1.5 p o P

1

1

0.5

0.5

0

Analiticko resenje resenje

x 10

0

1

2

3

4

5 Vreme

6

7

8

9

10 10

0

0

1

2

3

4

5 Vreme

6

7

8

9

10 10

Sledi još jedan primer u okviru kog će biti prikazano rešavanje neautonomne nelinearne diferencijalne jednač jednačine prvog reda. Zadatak 7: Na jednom kraju trkač trka čke staze postavljen je sigurnosni odbojnik koji zaustavlja vozilo nad kojim je vozač voza č izgubio kontrolu. Sila kojojm odbojnik deluje na vozilo funkcija je i brzine i pomaka prednje ploč ploče odbojnika: F

=

3

Kv 3 ( x + 1) ,

gde je K konstanta proporcionalnosti iznosi 30 s ÿkg/m5. Vozilo čija je masa m = 1500 kg, udara u odbojnik brzinom od 90 km/h. Nacrtati krivu brzine vozila kao funkciju položaja na intervalu intervalu od 0 do 3 m.

v x Korak 1 – Postavka problema: U trenutku kada automobil udari u odbojnik, koji je na slici prikazan kao opruga, na automobil deluje sila u smeru koji je suprotan od smera kretanja automobila i usporava automobil do njegovog zaustavljanja. Na osnovu II Newton-ovog zakona: 2009 © Katedra za Mikroelektroniku u tehičku f iziku, Elektrotehnički f akultet, Beograd


33

ma = − F ma = − Kv 3 (1 + x )

3

K 3 a = − v3 (1 + x ) , m

Kako se u zadatku traži grafik brzine u zavisnosti od pređ pre đenog puta (koordinate x), a ne vremena t , potrebno je izvršiti smenu promenljivih u diferencijalnoj jednač jednačini: a=

dv

=

dt

dv

1

⋅ =

dt

dv dx

=

dt dx

dv dx dx dt

=

dv dx

v.

pa diferencijalna jedna jednačina postaje: v

dv

=−

dx

dv dx

=−

K m

K m

v 3 (1 + x )

2 v (1 + x )

3

3

Korak 2 – Formiranje funkcijske datotetke: function dvdx = odbojnik(x,v) global k m % komanda global definise da su k i m globalne promenljive. Ovim se omogucava da se vrednost promenljivih definise u glavnom programu a da se koriste i u ovoj funkcijskoj datoteci. Moguce je i definisati ove promenljive diretno u samoj funkcijskoj datoteci tj. umesto komande global k m definisati k = 30; m = 1500; ali se time smanjuje fleksibilnost programa dvdx = -k/m*v^2*(x+1)^3; -k/m*v^2*(x+1)^3; % telo funkcije

Korak 3 i 4 – Odabir solvera i glavni program: >> global k m >> k = 30; m = 1500; v0 = 90*1000/3600; 90*1000/3600; >> [x v] = ode45('odbojnik',[0:0.1:3],v0); >> plot(x,v,'LineWidth',2); >> title('Brzina u funkciji od pomeraja x'); >> xlabel('x [m]'); ylabel('Brzina ylabel('Brzina [m/s]');



34 Brzina u funkciji od pomeraja x

25

20

] 15 s / m [ a n i z r B 10

5

0

0

0.5

1

1.5 x [m]

2

2. 5

3


uvod MatLab skripta-2010

Recommend Documents