plik


Obliczenia z u|yciem programu MATLAB Skrypt powstaB w ramach prac Centrum Korepetycji KoBa Naukowego Studentw Politechniki  Gambrinus . www.gambrinus.pwr.wroc.pl Autorzy: Karol BiaBows  rozdziaBy 1,2,4,7 Anna Borowska  rozdziaB 3 PaweB Gra  rozdziaB 5 Natalia Gemza - rozdziaB 6 Spis tre[ci: 1. Podstawy, proste obliczenia ....................................................................... 1 1.1 Interfejs ........................................................................................................ 1 1.2 Proste obliczenia .......................................................................................... 2 1.3 Wbudowane funkcje .................................................................................... 2 1.4 Zmienne ....................................................................................................... 3 1.5 Przydatne polecenia ..................................................................................... 4 1.6 Przenoszenie danych z arkusza kalkulacyjnego do Matlaba ....................... 5 2. Macierze i wektory ..................................................................................... 6 2.1 Sposoby definiowania macierzy i wektorw ............................................... 6 2.2 Dostp do elementw macierzy/wektorw .................................................. 7 2.3 Operacje na wybranym wierszu/kolumnie ................................................... 9 2.4 Podstawowe operacje rachunku macierzy ................................................... 10 2.5 Rozwizywanie ukBadw rwnaD liniowych metod macierzy odwrotnej .. 11 3. Wielomiany ................................................................................................ 12 3.1 Obliczanie warto[ci wielomianu (polyval) .................................................. 12 3.2 Obliczanie pierwiastkw wielomianu (roots) .............................................. 12 3.3 Mno|enie i dzielenie wielomianw (conv i deconv) ................................... 13 3.4 Pochodna wielomianu (polyder) .................................................................. 13 3.5 Aproksymacja wielomianowa (polyfit) ....................................................... 14 4. Funkcje i skrypty ........................................................................................ 15 4.1 Skrypty ........................................................................................................ 15 4.2 Funkcje ........................................................................................................ 17 4.3 Zmienne globalne ........................................................................................ 21 4.4 Funkcja jako argument innej funkcji ........................................................... 22 4.5 Znajdowanie miejsc zerowych dowolnej funkcji (fzero) ............................ 23 5. Ptle i instrukcje warunkowe ..................................................................... 25 5.1 Ptla for ....................................................................................................... 25 5.2 Ptla while ................................................................................................... 26 5.3 Instrukcja warunkowa if .............................................................................. 28 6. Grafika w Matlabie .................................................................................... 31 6.1 Funkcja plot ................................................................................................. 31 6.2 Edycja wykresu ............................................................................................ 31 6.3 Hold on ........................................................................................................ 35 6.4 Subplot ......................................................................................................... 36 6.5 Interpolacja .................................................................................................. 37 6.6 Funkcja fplot ................................................................................................ 38 7. CaBkowanie ................................................................................................ 40 7.1 Funkcja trapz ............................................................................................... 40 7.2 Funkcja ode45 ............................................................................................. 42 7.3 Funkcja ode45  caBkowanie ukBadw rwnaD ........................................... 45 7.4 Zdarzenia (Events) ...................................................................................... 46 1. Podstawy, proste obliczenia Matlab to zaawansowane [rodowisko do obliczeD matematycznych. Udostpnia ogromn ilo[ specjalistycznych funkcji oraz jzyk skryptowy pozwalajcy nawet tworzy programy z GUI (graficznym interfejsem u|ytkownika). W kursie tym postaramy si przybli|y niewielk cz[ najpotrzebniejszych in|ynierowi  chemikowi funkcji programu. 1.1 Interfejs Zacznijmy od przyjrzenia si interfejsowi. Okno po uruchomieniu programu Matlab. Elementy, ktre bd nas interesowa wyr|niono czerwonym kolorem. 1  przycisk uruchamiajcy edytor M-plikw 2  przycisk uruchamiajcy Simulink 3  Command Window 4  okno Workspace O dwch pierwszych pozycjach wicej powiemy w dalszych cz[ciach kursu. Na razie korzysta bdziemy jedynie z: % Command Window  okno, w ktrym mo|emy bezpo[rednio wpisywa polecenia i w ktrym zobaczymy wyniki ich wykonania % Workspace  w tym oknie znajduj si zmienne utworzone w czasie aktualnej sesji 1 1.2 Proste obliczenia Pierwsza operacja jak wykonamy to zwykBe dodawanie. Wpiszmy w oknie Command Window: >> 3+2 i wci[nijmy enter. Powinni[my zobaczy nastpujcy wynik: >> 3+2 ans = 5 >> Zwr uwag, |e w oknie Workspace pojawiBa si zmienna o nazwie ans i warto[ci 5. Do zmiennej tej bdziemy mogli si pzniej odwoBywa. Analogicznie mo|emy oblicza warto[ci wszystkich innych wyra|eD algebraicznych. Sprbuj teraz samodzielnie wykona nastpujce operacje i sprawdz ich wyniki: 3-2 3*2 3/2 3^2 Ostatnia z tych operacji to podnoszenie do potgi. Sprbujmy teraz obliczy warto[ nieco bardziej skomplikowanego wyra|enia: 2.3"e3.46"10-3 1.23-3 4.2 >> (2.3*exp(3.4) + 6e-3)/(1.23 - 3*sqrt(4.2)) ans = -14.0140 Zauwa|, |e w Matlabie mozesz stosowa nawiasy. Kolejno[ wykonywania dziaBaD jest standardowa, tj. potgowanie, nastpnie mno|enie/dzielenie itd., z uwgldnieniem nawiasw. W przykBadzie tym zwr uwag na to jak zapisano wyra|enie 6*10-3 - 6e-3. Jest to tak zwany zapis wykBadniczy. Warto pamita |e liczby mo|na w matlabie zapisywa wBa[nie w ten sposb. 1.3 Wbudowane funkcje Kolejna istotna rzecz pokazana w tym przykBadzie to funkcje. exp(3.4) oraz sqrt(4.2) Aatwo domy[li si co one oznaczaj  exp to funkcja eksponent czyli  e do x a sqrt to pierwiastek kwadratowy z liczby. Matlab udostpnia jeszcze wiele funkcji. Kilkana[cie czsto u|ywanych zebrano w tabeli poni|ej. Funkcja Co oblicza: exp(x) e^x sqrt(x) Pierwiastek kwadratowy log(x) Logarytm naturalny 2 log10(x) Logarytm dziesitny abs(x) Warto[ bezwzgldna sin(x), cos(x), tan(x) Warto[ci funkcji trygonometrycznych (x  w radianach) asin(x), acos(x, atan(x) Funkcje cyklometryczne, wynik w radianach deg2rad(x), rad2deg(x) Przeliczanie odpowiednio stopni na radiany (deg2rad) i radianw na stopnie(rad2deg) sinh(x),cosh(x),tanh(x) Funkcje hiperboliczne mod(x,y) Reszta z dzielenia x przez y round(x) Zaokrglanie do liczby caBkowitej floor(x) Zaokrglanie do najbli|szej liczby caBkowitej w gr ceil(x) Zaokrglanie do najbli|szej liczby caBkowitej w dB sum(X) Oblicza sum wszystkich elementw macierzy/wektora min(X), max(X) Warto[ najmniejszego/najwikszego elementu macierzy/wektora Sprbuj wykona kilka dziaBaD u|ywajc ka|dej z tych funkcji. W jednym z kolejnych rozdziaBw nauczysz si tak|e pisa wBasne funkcje. 1.4 Zmienne Podczas obliczeD w Matlabie mo|emy rwnie| wykorzystywa zmienne. Oto przykBad ich u|ycia: >> a=3 a = 3 >> b=2 b = 2 >> a+b ans = 5 >> c=log(a) c = 1.0986 >> d=b^c d = 2.1415 Zmienne pozwalaj zapamitywa warto[ci pod postaci r|nych symboli. Nie musz to by pojedyncze litery, mog to by praktycznie dowolne cigi znakw  wa|ne |eby znaki w nazwie nie byBy rozdzielone spacj. Po wykonaniu powy|szego cigu poleceD przyjrzyj si zawarto[ci okna Workspace. Powinny si w nim pojawi nazwy i warto[ci naszych zmiennych. Ze zmiennych bdziemy korzysta do samego koDca kursu zatem jeszcze zd|ysz si z nimi oswoi. Na razie przydadz Ci si one jedynie podczas wykonywania dBu|szych cigw obliczeD, gdy z danej warto[ci musisz skorzysta wielokrotnie. Zmienne przydadz Ci si rwnie| podczas obliczania skomplikowanych wyra|eD. Czasem Batwiej rozbi takie wyra|enie na kilka cz[ci i ka|d z nich obliczy oddzielnie  w krotkich wzorach Batwiej unikn bBdu. Oto prosty przykBad wykorzystania zmiennych w tym wBa[nie celu: 3 Sprbujmy obliczy warto[ wyra|enia: 3.2"ln2.2 2 3 2"sin 3 2 2"sin 2 3e2.1 2 2" " "e 2-ln2.8 2 >> a=(3+exp(2.1))/(2-log(2.8)) a = 11.5070 >> b=(2*sin(3/2*pi))/sqrt(2) b = -1.4142 >> c=exp(3.2*log(2.2)/b) c = 0.1680 >> wynik=2*a*b^2*c wynik = 7.7305 >> Powy|sze wyra|enie bez u|ycia zmiennych do przechowywania wynikw obliczeD po[rednich wyglda tak: >> wynik=2*(3+exp(2.1))/(2-log(2.8))*((2*sin(3/2*pi))/sqrt(2))^2* (exp(3.2*log(2.2)/((2*sin(3/2*pi))/sqrt(2)))) wynik = 7.7305 >> Przy takim zapisie znacznie Batwiej o bBd ktry mo|e by trudny do znalezienia, np.: >> wynik=2*(3+exp(2.1)/(2-log(2.8)))*((2*sin(3/2*pi))/sqrt(2))^2* (exp(3.2*log(2.2)/((2*sin(3/2*pi))/sqrt(2)))) wynik = 7.6690 1.5 Przydatne polecenia Na koniec jeszcze kilka przydatnych poleceD/cech Matlab-a o ktrych warto wiedzie. Historia poleceD: wciskajc klawisz strzaBki w gr w oknie Command Window cofamy si do wykonywanych wcze[niej operacji i mo|emy je powtrzy lub lekko zmodyfikowa i wykona t zmodyfikowan wersj. clc  komenda, przy pomocy ktrej czy[cimy zawarto[ okna Command Window who  wypisanie w oknie Command Window wszystkich utworzonych zmiennych ; - wyniki poleceD gdy lini zakoDczymy [rednikiem nie bd wy[wietlane, np.: >> a=3+2; >> a=3+2 a = 5 4 bdzie to przydatne pzniej  przy pisaniu wBasnych funkcji. format long/short  format wy[wietlania liczb rzeczywistych, np.: >> a=sqrt(2); >> a a = 1.4142 >> format long >> a a = 1.414213562373095 >> format short >> a a = 1.4142 clear  usuwa wszystkie zmienne utworzone w danej sesji clear nazwa_zmiennej  usuwa wybran zmienn 1.6 Przenoszenie danych z arkusza kalkulacyjnego do Matlab-a Aby przenie[ dane z wikszo[ci popularnych arkuszy kalkulacyjnych do programu Matlab, wystarczy te dane zaznaczy i po klikniciu prawego klawisza myszki wybra opcj kopiuj, nastpnie w matlabie nale|y utworzy tablic, do ktrej skopiujemy dane: >> tab1=[] tab1 = [] Nazwa tablicy pojawi si w okienku Workspace - tam nale|y klikn na ni dwukrotnie. Obok powinna pojawi si wtedy tabela w ktrej mo|emy w prosty sposb edytowa zawarto[ tablicy. Nale|y wtedy klikn prawym przyciskiem na komrk tabeli i z rozwijanego meu wybra paste aby wklei dane skopiowane wcze[niej z arkusza. Kopiowanie danych z arkusza kalkulacyjnego do Matlab-a 5 2. Macierze i wektory Warto wiedzie, |e nazwa programu Matlab pochodzi od Matrix laboratory. Zgodnie z tym co nazwa sugeruje, Matlab umo|liwia bardzo proste wykonanie wielu operacji na macierzach i wektorach. W zasadzie macierze i wektory stanowi typ danych, ktrym najcz[ciej bdziesz si posBugiwaB. 2.1 Sposoby definiowania macierzy i wektorw Definiowanie macierzy w Matlabie jest bardzo proste: >> A=[1 2 3 4;5 6 7 8;9 10 11 12] A = 1 2 3 4 5 6 7 8 9 10 11 12 >> Spacja oddziela kolejne wyrazy w wierszu, natomiast przej[cie do nowej linii lub [rednik oznaczaj rozpoczcie kolejnego wiersza. Jeszcze jeden przykBad: >> B=[1 2 3 4 5 6 7 8] B = 1 2 3 4 5 6 7 8 >> Istniej jeszcze trzy bardzo przydatne sposoby tworzenia wektorw (macierzy posiadajcych tylko jeden wiersz lub tylko jedn kolumn) w programie Matlab. Tymi sposobami otrzymamy wektory zawierajce warto[ci z okre[lonego przedziaBu w odstpach rwnych, lub w postaci cigu arytmetycznego. Sposb pierwszy prezentuje poni|szy przykBad: >> C=[0:10:50] C = 0 10 20 30 40 50 >> D=[0:5:23] D = 0 5 10 15 20 >> Oglnie mo|emy to zapisa jako [pocztek : odstp : warto[ maksymalna] Zauwa|, |e do wektora trafiaj tylko wyrazy mniejsze lub rwne warto[ci maksymalnej. 6 Dwa kolejne sposoby tworzenia wektora to funkcje linspace i logspace. >> linspace(1,50,5) ans = 1.0000 13.2500 25.5000 37.7500 50.0000 >> pierwszy i drugi argument funkcji linspace to warto[ci pierwszego i ostatniego elementu wektora, ostatni  liczba elementw wektora. Elementy wektora tworzonego funkcj linspace s od siebie rwno oddalone. Funkcja logspace generuje wektor w ktrym kolejne elementy to warto[ci 10^x. Najlepiej zilustrowa dziaBanie tej funkcji przykBadem: >> logspace(-2,2,5) ans = 0.0100 0.1000 1.0000 10.0000 100.0000 >> Otrzymali[my wektor zawierajcy 5 warto[ci, od 10^-2 do 10^2. Warto[ci to kolejno 10^-2, 10^-1, 10^0 itd. Przydatne polecenia do tworzenia macierzy to ones(m), eye(m), zeros(m) oraz te same funkcje wywoBywane z dwoma argumentami, tj.: ones(m,n) itd.... Funkcje te tworz macierze kwadratowe o wymiarze m gdy wywolamy je z jednym argumentem lub macierze posiadajace m wierszy i n kolumn je[li wywoBamy je z dwoma argumentami. Sprbuj u|y ka|dej z tych funkcji w obu formach aby zobaczy jakie macierze powstaj w wyniku ich dziaBania. 2.2 Podstawowe operacje rachunku macierzy Najprostsze operacje rachunku macierzy to dodawanie i odejmowanie macierzy. Macierze musz mie ten sam wymiar. Dodawanie i odejmowanie macierzy wyglda dokBadnie tak samo jak dodawanie i odejmowanie od siebie dwu zmiennych. Oto przykBad dodawania macierzy: >> A=[1 2 3;4 5 6] A = 1 2 3 4 5 6 >> B=[1 2 3;7 8 9] B = 1 2 3 7 8 9 >> C=A+B C = 2 4 6 11 13 15 >> Analogicznie wyglda odejmowanie macierzy (sprbuj odj od siebie macierze A i B). Mno|enie macierzy przez liczb jest rwnie proste: >> A=[1 2 3;4 5 6] A = 1 2 3 4 5 6 7 >> A*2 ans = 2 4 6 8 10 12 >> Mno|enie macierzy i podnoszenie jej do potgi jest ju| nieco bardziej skomplikowane. Mamy tu do wyboru dwie opcje. Albo mno|ymy macierze zgodnie z zasadami rachunku macierzy albo mno|ymy/dzielimy przez siebie odpowiadajce sobie elementy macierzy. W poni|szej tabeli podano odpowiednie przykBady. Najpierw zdefiniujemy dwie macierze: >> A=[1 2 3;4 5 6;7 8 9] A = 1 2 3 4 5 6 7 8 9 >> B=[11 12 13;14 15 16;17 18 19] B = 11 12 13 14 15 16 17 18 19 >> Wyniki wykonania komend z lewej kolumny znajduj si w prawej kolumnie tabeli. >> A*B 90 96 102 216 231 246 342 366 390 >> A.*B 11 24 39 56 75 96 119 144 171 >> A./B 0.0909 0.1667 0.2308 0.2857 0.3333 0.3750 0.4118 0.4444 0.4737 Zatem zapis bez kropki poprzedzajcej symbol mno|enia/dzielenia oznacza zwykBe mno|enie/dzielenie macierzy, natomiast w przypadku u|ycia operatora .* lub ./ wymna|ane/dzielone przez siebie s odpowiadajce sobie elementy macierzy. Analogicznie wyglda potgowanie. W wypadku u|ycia operatora ^, Matlab wykonuje zwykle mno|enie macierzy, natomiast .^ oznacza podniesienie ka|dego elementu macierzy do danej potgi. Aby lepiej zrozumie r|nic, przyj|yj si poni|szemu przykBadowi: >> A=[2 2 2; 2 2 2; 2 2 2] A = 2 2 2 2 2 2 2 2 2 >> B=A^3 B = 72 72 72 72 72 72 72 72 72 >> C=A.^3 C = 8 8 8 8 8 8 8 8 8 8 Matlab umo|liwia tak|e proste wykonanie operacji takich jak transpozycja macierzy, wyliczenie wyznacznika, warto[ci wBasnych oraz wyznaczenie macierzy odwrotnej. Aby otrzyma macierz transponowan, nale|y u|y operatora ' : >> A=[1 1 1;2 2 2;3 3 3] A = 1 1 1 2 2 2 3 3 3 >> B=A' B = 1 2 3 1 2 3 1 2 3 >> PozostaBe operacje: % inv(X)  wyznaczenie macierzy odwrotnej do macierzy X % det(X)  obliczenie wyznacznika macierzy X % eig(X)  wyznaczenie warto[ci wBasnch macierzy X 2.3 Dostp do elementw macierzy/wektorw Na pocztek utwrzmy przykBadow macierz wielko[ci 7x7: >> A=magic(7) A = 30 39 48 1 10 19 28 38 47 7 9 18 27 29 46 6 8 17 26 35 37 5 14 16 25 34 36 45 13 15 24 33 42 44 4 21 23 32 41 43 3 12 22 31 40 49 2 11 20 >> W przykBadzie poszli[my troch na skrty  komenda magic tworzy macierz kwadratow o zadanym wymiarze, w ktrej sumy wszystkich wierszy i kolumn s rwne i w ktrej nie powtarza si |aden element. Odczytanie pojedynczego elementu: >> A(1,1) ans = 30 >> A(4,5) ans = 34 >> Jako pierwszy argument podajemy wiersz, jako drugi kolumn, wiersze i kolumny macierzy numerowane s od 1. 9 Odczytanie wybranego fragmentu macierzy wykonujemy tak: >> A(2,2:5) ans = 47 7 9 18 Operator : pozwala zdefiniowa zakres interesujcych nas elementw. Tu odczytali[my elementy od 2 do 5 z wiersza 2. Operatora dostpu : mo|emy u|y rwnie| do odczytania caBych wierszy/kolumn  wtedy nie podajemy zakresu jaki chcemy odczyta a jedynie wstawiamy sam operator :, tak jak w poni|szych przykBadach: >> A(:,2) ans = 39 47 6 14 15 23 31 >> A(2,:) ans = 38 47 7 9 18 27 29 >> Jak si Batwo domy[li A(:,:) oznacza caB zawarto[ macierzy A. 2.4 Operacje na wybranym wierszu/kolumnie Zajmiemy si teraz operacjami w ktrych chcemy wykona dziaBania tylko na wybranym wierszu lub wybranej kolumnie macierzy, takimi jak: mno|enie wiersza/kolumny przez liczb, dodawanie wierszy/kolumn, usuwanie wierszy/kolumn >> A=magic(7); >> A(3,:)=2*A(3,:) A = 30 39 48 1 10 19 28 38 47 7 9 18 27 29 92 12 16 34 52 70 74 5 14 16 25 34 36 45 13 15 24 33 42 44 4 21 23 32 41 43 3 12 22 31 40 49 2 11 20 >> A(4,:)=A(4,:)-A(3,:) A = 30 39 48 1 10 19 28 38 47 7 9 18 27 29 92 12 16 34 52 70 74 -87 2 0 -9 -18 -34 -29 13 15 24 33 42 44 4 21 23 32 41 43 3 12 22 31 40 49 2 11 20 >> W powy|szym przykBadzie najpierw pomno|ono wiersz 3 macierzy przez 2 a nastpnie od wiersza 4 odjto wiersz 3. Wybrany wiersz/kolumn macierzy usuniemy w nastpujcy sposb: 10 >> A(:,3)=[] A = 30 39 1 10 19 28 38 47 9 18 27 29 92 12 34 52 70 74 -87 2 -9 -18 -34 -29 13 15 33 42 44 4 21 23 41 43 3 12 22 31 49 2 11 20 >> W przykBadzie tym usunli[my kolumn 3 macierzy bdcej wynikiem poprrzednich dwu operacji. 2.5 Rozwizywanie ukBadw rwnaD liniowych metod macierzy odwrotnej - przykBad Na koniec rozwi|emy jeszcze prosty ukBad rwnaD wykorzystujc metod macierzy odwrotnej: x-2y3z=-7 3x y4z=5 2x5yz=18 UkBad zapiszemy w postaci macierzowej: A"X=B 1 -2 3 x -7 3 1 4"y= 5 2 5 1 z 18 Rozwizanie ukBdu znajdziemy dziki wBasno[ciom macierzy odwrotnej: A"X=B! X =A-1"B Nastpnie zdefiniujmy w Matlabie odpowiednie macierze A i B oraz wyznaczymy rozwizanie ukBadu: >> A=[1 -2 3;3 1 4;2 5 1] A = 1 -2 3 3 1 4 2 5 1 >> B=[-7;5;18] B = -7 5 18 >> X=inv(A)*B X = 2.0000 3.0000 -1.0000 Zatem rozwizanie ukBadu rwnaD to: x=2, y=3 i z=-1. 11 3. Wielomiany Wielomiany w matlabie przechowywane s w postaci wierszowego wektora wspBczynnikw. Pierwszy wyraz wektora zawiera wspBczynnik przy najwy|szej potdze zmiennej niezale|nej, kolejne wyrazy to wspBczynniki przy kolejnych, coraz ni|szych potgach zmiennej, ostatni to wyraz wolny wielomianu. Dla przykBadu wielomian: Wprowadzamy do matlaba w postaci wektora: >> p=[3,2,-1] 3.1 Obliczanie warto[ci wielomianu (polyval): Funkcja polyval(p,x) oblicza warto[ wielomianu o wspBczynnikach zawartych w wektorze p w punktach wyspecyfikowanych w zmiennej x. Dla przykBadu posBu|ono si powy|szym wielomianem p(x): >> p=[3,2,-1] p = 3 2 -1 >> polyval(p,1) ans = 4 Oczywi[cie mo|na od razu wprowadzi kilka zmiennych x w postaci wektora: >> polyval(p,[1,3,0,7]) ans = 4 32 -1 160 3.2 Obliczanie pierwiastkw wielomianu (roots): Funkcja roots pozwala obliczy pierwiastki wielomianu. PrzykBad: >> roots(p) ans = -1.0000 0.3333 Oczywi[cie Matlab oblicza tak|e pierwiastki zespolone: >> q=[3,1,0,-1]; >> roots(q) ans = -0.4658 + 0.5834i -0.4658 - 0.5834i 0.5982 12 3.3 Mno|enie i dzielenie wielomianw (conv i deconv): Funkcja conv pozwala obliczy wektor wspBczynnikw wielomianu bdcego iloczynem wielomianw reprezentowanych przez jej argumenty, tj: >> p=[1,1]; >> q=[1,2,2]; >> c=conv(p,q) c = 1 3 4 2 Funkcja deconv pozwala obliczyc iloraz oraz reszt z dzielenia jednego wielomianu przez drugi: >> p=[2,2,2,1]; >> q=[1,2]; >> [w,r]=deconv(p,q) w = 2 -2 6 r = 0 0 0 -11 3.4 Pochodna wielomianu (polyder): Funkcja polder(p) oblicza wektor wspBczynnikw wielomianu bdcego pochodn wielomianu reprezentowanego przez p. >> p=[3,2,-1]; >> polyder(p) ans = 6 2 Mo|na tak|e obliczy pochodn iloczynu wielomianw  po prawej stronie znaku rwno[ci znajduje si pojedyncza zmienna: >> p=[3,2,-1]; >> q=[4,0,-1,1]; >> k=polyder(p,q) k = 60 32 -21 2 3 Powy|szy wynik mo|na uzyska stosujc znan ju| funkcj conv: >> p=[3,2,-1]; >> q=[4,0,-1,1]; >> polyder(conv(p,q)) ans = 60 32 -21 2 3 13 Oczywi[cie w Batwy sposb mo|na obliczy pochodn ilorazu wektorw. SkBadnia tej funkcji jest taka sama jak w przypadku funkcji deconv  po prawej stronie znaku rwno[ci znajduje si tablica z dwiema zmiennymi: >> p=[3,2,-1]; >> q=[4,0,-1,1]; >> [w,r]=polyder(q,p) w = 12 16 -9 -6 -1 r = 9 12 -2 -4 1 gdzie: w  wynik dzielenia wielomianw, r  reszta z dzielenia. 3.5 Aproksymacja wielomianowa (polyfit): GBwn funkcj wykorzystywan przy aproksymacji wielomianowej jest polyfit. Oblicza ona wspBczynniki wielomianu zadanego stopnia aproksymujcego dane wej[ciowe w sensie minimum sumy kwadratw. SkBadnia funkcji: >> p=polyfit(x,y,n) x,y  wektory zawierajce dane n  stopieD wielomianu jaki chcemy otrzyma >> x=[1 3 5 7 9 11]; >> y=[1 13 17 23 31 39]; >> p=polyfit(x,y,3) p = 0.0463 -0.8512 7.8505 -5.5853 >> x2=1:0.5:11; >> y2=polyval(p, x2); >> plot(x,y,'x',x2,y2) 14 15 4. Funkcje i skrypty Matlab posiada wBasny, rozbudowany jzyk skryptowy pozwalajcy Batwo tworzy zaawansowane programy. W tym rozdziale omwimy krtko tworzenie prostych skryptw i funkcji. Skrypt to cig poleceD zapisany w pliku, ktre po uruchomieniu skryptu zostaj kolejno przez Matlab-a wykonane. Z funkcji korzystaBe[ ju| wielokrotnie zatem pewnie wiesz ju| dobrze co to jest i do czego sBu|y. Skrypty i funkcje zapisywane s w postaci tzw. m-plikw, czyli plikw z rozszerzeniem .m. W przypadku funkcji, ka|da funkcja musi by zdefiniowana w oddzielnym m- pliku o nazwie takiej, jak nazwa danej funkcji. Aby otworzy edytor plikw m wybierz z menu: File ! New ! M-File Powinien zosta uruchomiony edytor (u Ciebie mo|e wyglda troch inaczej): Okno edytora m-plikw 4.1 Skrypty Napiszmy skrypt obliczajcy czas wypBywu cieczy z cylindrycznego zbiornika. Jako dane potrzebne bd (w nawiasach podano przyjte w skrypcie warto[ci): D  [rednica zbiornika (3m) h  wysoko[ do ktrej zbiornik jest napeBniony (4m) d  [rednica otworu, przez ktry wypBywa ciecz (100mm = 0.1m) - wspBczynnik wypBywu (0.8) g  przyspieszenie ziemskie (9.81 m/s2) CaBkowity czas wypBywu ze zbiornika w sekundach obliczymy ze wzoru: 2"D2 twyp= " h u"d2 2g 16 W edytorze m-plikw wprowadz nastpujc tre[: D=3 h=4 d=0.1 u=0.8 g=9.81 t_wyp=(2*D^2)/(u*d^2*sqrt(2*g))*sqrt(h) Teraz z menu edytora m-plikw wybieramy Debug ! Safe file and Run (lub podobn opcj  nazwy mog si nieznacznie r|ni w zale|no[ci od wersji programu Matlab). Zostaniesz zapytany pod jak nazw zapisa skrypt  nazwij skrypt t_wyplywu. Po zapisaniu skrypt zostanie uchomiony. Przejdz do okna Command Window. Tam znajdziesz wyniki obliczeD: D = 3 h = 4 d = 0.1000 u = 0.8000 g = 9.8100 t_wyp = 1.0159e+03 >> Inna metoda wywoBania skryptu to wpisanie jego nazwy w wierszu poleceD. Teraz wr do edytora i zmieD wysoko[ cieczy w zbiorniku na 2 m. Ponownie wybierz z menu Debug opcj Safe file and Run. Otrzymasz czas opr|niania zbiornika dla nowych danych. Tym razem program nie powinien pyta o nazw skryptu (skrypt zostanie ponownie zapisany pod podan wcze[niej nazw). Wynik powinien by nastpujcy: D = 3 h = 2 d = 0.1000 u = 0.8000 g = 9.8100 t_wyp = 718.3697 >> Sprbuj jeszcze doda znak [rednika (;) na koDcu pierwszych piciu linii: D=3; h=4; d=0.1; u=0.8; g=9.81; t_wyp=(2*D^2)/(u*d^2*sqrt(2*g))*sqrt(h) 17 Sprawdz jak wtedy bdzie wygldaB wynik uruchomienia skryptu. Je[li nie wiesz dlaczego tak jest, wr do cz[ci 1.5 pierwszego rozdziaBu kursu  tam dowiesz si czym skutkuje zakoDczenie linii [rednikiem. 4.2 Funkcje Funkcjom po[wicimy troch wicej uwagi. Bdziesz je wykorzystywaB i pisaB bardzo czsto podczas pracy z programem Matlab. WBasne funkcje bdziesz pisaB na przykBad podczas numerycznego obliczania caBek. Z menu (edytora m-plikw lub Matlaba) wybierz: File ! New ! M-File Funkcja, ktr stworzymy na pocztek bdzie oblicza objto[ molow gazu doskonaBego: R"T p"v=R"T ! v p , T = p function wynik=v(p,T) %funkcja oblicza objto[ molow [m^3/mol] %p - ci[nienie [Pa] %T - temperatura [K] R=8.314; %[J/(mol*K)] wynik=R*T/p; end Teraz zapisujemy nasz plik. Plik musi nazywa si dokBadnie tak samo jak funkcja  czyli w naszym przypadku bdzie to v.m i tak wBa[nie nazw zaproponuje nam Matlab. Zapamitaj, w jakim dokBadnie folderze zapisujesz plik. Teraz przejdz do okna Command Window i wpisz: >> v(101300,273.15) ans = 0.0224 >> Obliczyli[my w ten sposb objto[, zajmowan przez mol gazu doskonaBego w warunkach normalnych  22.4 dm3 czyli 0.0224m3. Je[li otrzymaBe[ komunikat: >> v(101300,273.15) ??? Undefined function or method 'v' for input arguments of type 'double'. >> Nie przejmuj si. Prawdopodobnie folder w ktrym zapisaBe[ swoj funkcj nie znajduje si na li[cie folderw w ktrych Matlab szuka funkcji (tzw. Path). Musisz ten folder do [cie|ki doda. Z menu w gBwnym oknie programu wybierz: File ! Set Path... W oknie, ktre si pojawi wybierz Add Folder... a nastpnie wybierz folder, w ktrym zapisaBe[ swoj funkcj. Teraz wszystko powinno dziaBa. Przyjrzyjmy si teraz bardziej szczegBowo napisanej przez nas funkcji. 18 function wynik=v(p,T) Plik z definicj funkcji powinien si zaczyna sBowem kluczowym function  informuje ono Matlab-a, |e zaraz zdefiniujemy funkcj. Po sBowie kluczowym function podajemy nazwy zmiennych, ktrych warto[ci funkcja zwrci. W naszym przypadku jest to pojedyncza zmienna o nazwie wynik ale mo|e to by rwnie| wektor lub macierz (przykBady w dalszej cz[ci rozdziaBu). Nastpnie znak rwno[ci i nazwa funkcji, po ktrej w nawiasie znajduje si lista argumentw przyjmowanych przez funkcj (u nas s to p i T). %funkcja oblicza objto[ molow [m^3/mol] %p - ci[nienie [Pa] %T - temperatura [K] W kolejnych trzech linijkach zaczynajcych si od znaku % umieszczono komentarze mwice co funkcja liczy oraz w jakich jednostkach podawa agumenty wej[ciowe. Linie zaczynajce si od znaku % s ignorowane przez Matlab-a, mo|emy w nich wpisa co chcemy  zwykle wBa[nie jednostki lub jaki[ krtki opis dziaBania znajdujcego si w kolejnych linijkach kodu. Komentarze takie nie s konieczne do dziaBania funkcji ale warto je pisa, |eby potem nie mie wtpliwo[ci co dokBadnie robi nasza funkcja i jakich parametrw od nas oczekuje. R=8.314; %[J/(mol*K)] wynik=R*T/p; end Dalej podobnie jak w skrypcie piszemy cig poleceD, ktre Matlab ma wykona, w tej cz[ci obliczamy warto[ci ktre funkcja ma zwrci  my w tej cz[ci definiujemy zmienn R  uniwersaln staB gazow oraz obliczamy wynik  molow objto[ gazu. W obliczeniach wykorzystujemy przekazane do funkcji argumenty (w naszym przykBadzie argumenty to p i T). Definiowanie funkcji powinni[my koDczy sBowem kluczowym end. Poni|ej przykBad funkcji, zwracajcej kilka warto[ci: function [f g]=fg(x) %funkcja oblicza warto[ci dwu funkcji : %f(x)=2x+2 %g(x)=x^2+2x+1 f=2*x+2; g=x^2+2*x+1; end >> [j k]=fg(2) j = 6 k = 9 >> a=fg(2); >> Przyj|yj si ostatniemu wydanemu poleceniu  jak my[lisz, co powinno si znalez w zmiennej a? Pewnie spodziewasz si, |e a bdzie wektorem zawierajcym warto[ci 6 i 9, czyli: [6 9] Nie. Zmienna a bdzie miaBa warto[ 6. Zostanie pod ni podstawiona pierwsza warto[ z wektora zwracanego przez funkcj. Pamitaj o tym  dziki temu unikniesz wielu bBdw. Kolejnym czsto 19 popeBnianym bBdem jest u|ywanie w funkcjach operatorw mno|enia i potgowania macierzowego zamiast operatorw wykonujcych te operacje  element po elemencie . Aby zilostrowa ten problem, stwrzmy kolejn funkcj  podobn do fg(x): function [f g]=fg1(x) %funkcja oblicza warto[ci dwu funkcji : %f(x)=2x+2 %g(x)=x^2+2x+1 f=2.*x+2; g=x.^2_2.*x+1; end R|nica midzy fg i fg1 polega wyBcznie na tym, |e operatory mno|enia i potgowania macierzowego (* i ^) zastpiono w funkcji fg1 operatorami mno|enia i potgowania  element po elemencie (.^, .*). Wykonanie poni|szych poleceD pozwoli nam zobaczy r|nice midzy obiema funkcjami. >> fg1(2) ans = 6 >> [j k]=fg1(2) j = 6 k = 9 >> x=[1:1:5] x = 1 2 3 4 5 >> [j k]=fg(x) ??? Error using ==> mpower Matrix must be square. Error in ==> fg at 6 g=x^2+2*x+1; >> [j k]=fg1(x) j = 4 6 8 10 12 k = 4 9 16 25 36 >> Jak wida dla pojedynczego argumentu funkcje fg i fg1 daj dokBadnie takie same rezultaty. W kolejnej operacji tworzymy wektor x. Chcemy obliczy warto[ci funkcji f(x) i g(x) dla ka|dego elementu wektora. W przypadku u|ycia funkcji fg otrzymujemy bBd: ??? Error using ==> mpower Matrix must be square. Error in ==> fg at 6 g=x^2+2*x+1; BBd wynika z tego, |e chcemy podnosi do potgi macierz, ktra nie jest kwadratow, co zgodnie z definicj potgowania macierzy jest operacj niedozwolon. Funkcja fg1(x) daje po|dany rezultat  w wektorach k i j znajduj si odpowiednio warto[ci funkcji f i g dla kolejnych liczb w wektorze x. Podobnie w przypadku operatorw mno|enia i dzielenia: 20 function z=f(x,y) z=y*x; end function z=f1(x,y) z=y.*x; end >> fun1(1,2) ans = 2 >>fun(1,2) ans = 2 >> fun1(x,y) ans = 0 0.5000 2.0000 4.5000 8.0000 >> fun(x,y) ??? Error using ==> mtimes Inner matrix dimensions must agree. Error in ==> fun at 2 z=y*x; >> W funkcjach zwykle powinni[my u|ywa operatorw wykonujcych operacje  element po elemencie , zatem zapamitaj aby operatorw macierzowych (czyli tych normalnych,  bez kropki , ktrych poza Matlabem u|ywamy praktycznie zawsze) u|ywa tylko wtedy kiedy naprawd chodzi nam o wykonanie operacji potgowania/mno|enia/dzielenia macierzy. Poni|ej znajduje si prosty przykBad pokazujcy dziaBanie [rednika na koDcu linii wewntrz funkcji. Zmodfikuj funkcj fg1(x), usuwajc [rednik z koDca linii, w ktrej liczona jest warto[ g(x): function [f g]=fg1(x) %funkcja oblicza warto[ci dwu funkcji : %f(x)=2x+2 %g(x)=x^2+2x+1 f=2.*x+2; g=x.^2_2.*x+1 end >> x=[1:1:5] x = 1 2 3 4 5 >> k=fg1(x) g = 4 9 16 25 36 k = 4 6 8 10 12 >> k=fg1(x); g = 4 9 16 25 36 >> Zwr uwag, |e gdy linijka w funkcji nie jest zakoDczona [rednikiem, wynik przeprowadzonej w niej operacji zawsze zostanie wy[wietlony (zwr uwag na ostanie wywoBanie). Sprbuj usun [rednik rwnie| z linii, w ktrej liczona jest warto[ funkcji f(x). 21 4.3 Zmienne globalne Oprcz argumentw, r|ne warto[ci mo|emy przekazywa do funkcji tak|e jako zmienne globalne. Nie jest to  elegancki sposb i raczej powinno si go unika, jednak czasami sposb ten pozwala znacznie upro[ci skrypty i funkcje. Na pocztek utwrz funkcj: function fn1(k) k=k+1 m=m+2 end Nastpnie wykonaj nastpujce polecenia: >> j=4 j = 4 >> m=3 m = 3 >> fn1(j) k = 5 ??? Undefined function or variable 'm'. Error in ==> fn1 at 3 m >> j j = 4 >> Tworzymy dwie zmienne  j oraz m. Nastpnie wywoBujemy funkcj, przekazujc jej jako argument m. W funkcji do argumentu dodawana jest liczba jeden ale zauwa|, |e po zakoDczeniu dziaBania funkcji warto[ j pozostaBa niezmieniona. Funkcja zgBasza rwnie| bBd  nie zdefiniowano zmiennej m. Teraz zmodyfikuj funkcj fn1: function fn1(k) global m k=k+1 m=m+1 end A nastpnie wydaj nastpujce polecenia w oknie Command Window: >> clear >> global m >> j=3 j = 3 >> m=4 m = 4 >> fn1(j) k = 4 m = 5 >> m 22 m = 5 >> j j = 3 Teraz zmienn m definiujemy jako globaln. Pamitaj, |e aby u|y jakiej[ zmiennej jako globalnej musisz zdeklarowa j jako globaln zarwno w swoim skrypcie/w Command Window jak te| w ka|dej funkcji, ktra z niej korzysta. Nastpnie tworzymy zmienn j o warto[ci 3 oraz zmiennej globalnej m nadajemy warto[ 4. WywoBujemy funkcj fn1  podobnie jak poprzednio. Zauwa|, |e w wyniku dziaBania funkcji warto[ zmiennej m zmieniBa si. Zatem kolejn rzecz o ktrej musisz pamita je[li u|ywasz zmiennych jest to, |e po zakoDczeniu dziaBania funkcji operujcej na takich zmiennych ich warto[ci nie powracaj automatycznie do stanu sprzed wywoBania funkcji. 4.4 Funkcja jako argument innej funkcji Przekazanie funkcji jako argumentu do innej funkcji jest bardzo przydatn mo|liwo[ci. Skorzystasz z niej np. podczas caBkowania numerycznego funkcji. Stworz nastpujce pliki z funkcjami: function y=f1(x) y=2*x.^2+3; end function y=f2(x) y=2*exp(x+1); end function y=f3(x) y=2*x+1; end Teraz stworzymy przykBadow funkcj, do ktrej bdziemy przekazywa nasze funkcje jako argumenty. Funkcja ta zwrci r|nic warto[ci funkcji prekazanej jako fn_1 i funkcji przekazanej jako fn_2 dla argumentu x. function y=diff_fun(fn_1,fn_2,x) y=fn_1(x)-fn_2(x); end Nastpnie przejdz do okna Command Window i sprbuj wykona nastpujce operacje: >>format long >> x=[1:1:10] x = 1 2 3 4 5 6 7 8 9 10 >> diff_fun(@f1,@f2,x) ans = 1.0e+05 * Columns 1 through 3 -0.000097781121979 -0.000291710738464 -0.000881963000663 Columns 4 through 6 -0.002618263182052 -0.007538575869855 -0.021182663168569 Columns 7 through 9 -0.058609159740835 -0.160751678551508 -0.438879315896134 23 Column 10 -1.195452834303956 >> f1(1)-f2(1) ans = -9.778112197861301 >> -0.000097788e5 ans = -9.778800000000000 >> diff_fun(f2,f3,x) ??? Input argument "x" is undefined. Error in ==> f2 at 2 y=2*exp(x+1); >> diff_fun(@f2,@f3,x) ans = 1.0e+05 * Columns 1 through 3 0.000117781121979 0.000351710738464 0.001021963000663 Columns 4 through 6 0.002878263182052 0.007958575869855 0.021802663168569 Columns 7 through 9 0.059469159740835 0.161891678551508 0.440339315896134 Column 10 1.197272834303956 >> Zwr uwag na sposb wywoBania funkcji, do ktrej jako parametry przekazujemy inne funkcje. diff_fun(@f1,@f2,x) Przed nazw przekazywanej funkcji musimy umie[ci znak @. Jest to operator zwracajcy tzw. uchwyt funkcji, wystarczy jednak, |eby[ zapamitaB |e gdy przekazujesz jako argument funkcj, przed jej nazw powinienne[ umie[ci znak @ - inaczej mo|esz otrzyma bBdy, takie jak np. ten: >> diff_fun(f2,f3,x) ??? Input argument "x" is undefined. Error in ==> f2 at 2 y=2*exp(x+1); Unikniesz ich stawiajc znak @ przed nazw funkcji przekazywanej jako argument do innej funkcji. 4.5 Znajdowanie miejsc zerowych dowolnej funkcji (fzero) Jednym z przykBadw funkcji do ktrych przekazujesz jako argument wBasn funkcj, jest wbuowana w Matlab-a fzero. Pozwala ona znalez punkt, w ktrym Twoja funkcja zwraca warto[ zero. Funkcj t mo|na zastosowa do funkcji jednej zmiennej. Jej wywoBanie wyglda tak: x_szukane = fzero(@fn,x0) pod zmienn x_szukane podstawiony zostanie wynik  x, przy ktrym warto[ zwracana przez funkcj fn jest rwna zero. x0 to punkt, od ktrego nale|y rozpocz poszukiwania. Miejsce zerowe powinno znalez si w okolicy tego punktu. Je[li wiemy na przykBad do jakiego przedziaBu ma nale|e poszukiwana warto[  mo|emy jako x0 poda punkt, bdcy [rodkiem tego przedziaBu. DziaBanie funkcji polega na znalezieniu w okolicy punktu x0 przedziaBu, na ktrego koDcach 24 warto[ci zwracane przez funkcj maj r|ne znaki a nastpnie dokBadne zlokalizowanie miejsca zerowego w tym przedziale. Rozwi|emy nastpujcy problem: Majc dane rwnanie: 2 11.17"1-x x = " 1-x 1-1.17"x Nale|y znalez warto[ x z przedziaBu (0;1), dla ktrej =1.76 . Aby mc u|y funkcji fzero, musimy przeksztaBci rwnanie do postaci: 2 11.17"1-x x 0= " - 1-x 1-1.17"x Tworzymy funkcj: function y=f4_5(x) a=x./(1-x); b=(1+1.17*(1-x))./(1-1.17*x); y=a.^2*b  1.76; end Zapisujemy j i przechodzimy do okna CommandWindow. Tam wydajemy komend: >> x=fzero(@f4_5,0.5) x = 0.4217 >> Otrzymujemy punkt, w ktrym warto[ funkcji f4_5 jest rwna zero. 25 5. Ptle i instrukcje warunkowe 5.1 Ptla for ZaB|my, |e chcemy wypeBni jednowymiarow tablic liczbami od 1 do 10. Wystarczy wic, |e zastosujemy nastpujcy fragment kodu: >> A(1) = 1; >> A(2) = 2; >> A(3) = 3; i tak dalej a| do 10. Nieco prostszym i efektywniejszym sposobem jest np. zastosowanie nastpujcego kodu: >> A = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]; Wci| jednak nie jest to jednak metoda najefektywniejsza. Gdyby[my musieli wypeBni tablic, ktra posiada 1000 pozycji, jej wypeBnianie byBoby dosy mczce. Tutaj pojawia si mo|liwo[ zastosowania ptli. DziaBanie ptli mo|na okre[li w bardzo prosty sposb: powtarza pewien fragment programu okre[lon ilo[ razy. DziaBanie ptli przybli| teraz przykBadem: function f = fun(x) f = x^2 +12*x -3; Tworzymy skrypt o nazwie run.m: k=5; for a=1:k fun(a) end Nastpnie w oknie Command Window wydajemy polecenie: >> run Naciskamy enter i na ekranie powinien pojawi si nastpujcy kod: ans = 10 ans = 25 ans = 42 ans = 61 ans = 82 Teraz wyja[nienia: Mamy funkcj dan nastpujcym wzorem: 26 f(x) = x2 + 12x  3 Chcemy obliczy warto[ tej funkcji dla x od 1 do 5 i u|y do tego celu programu Matlab. W tym celu: % tworzymy plik  fun.m , w tym pliku umieszczamy funkcj, ktra nastpnie bdzie wywoBywana dla kolejnych warto[ci  x % tworzymy plik  run.m , ten plik jest plikiem wykonawczym i w nim umieszczamy ptl % w  command window wpisujemy  run i otrzymujemy wyniki dziaBania kodu Poniewa| zakBadam, |e czytelnik posiadB ju| umiejtno[ pisania wBasnych funkcji, przejd do omwienia punktu, w ktrym tworzymy plik  run.m . W 1. linijce kodu definiujemy zmienn  k i przypisujemy jej warto[  5 . W 3. linijce kodu rozpoczyna si ptla  for . Budowa ptli for jest nastpujca: n= 1; k= 5; for a = n:k  code end Ptla zaczyna si sBowem kluczowym  for . Nastpnie podawany jest warunek dziaBania ptli. Ptla bdzie si wykonywa k-n+1 razy. W wy|ej wymienionym przypadku bdzie to 5 razy, ale za k i n mo|na podstawi dowolne liczby naturalne oraz 0. W miejscu gdzie znajduje si wyraz  code , umieszcza si wszystkie instrukcje, ktre maj by powtarzane. W rozpatrywanym wy|ej przykBadzie w miejscu  code umieszczone jest wywoBanie funkcji od zmiennej  a . Zmienna  a , nazywana te| licznikiem ptli, w trakcie dziaBania ptli zmienia swoj warto[ w zakresie od n do k z krokiem rwnym  1 . Tak wic za ka|dym razem gdy zostanie wywoBana nasza funkcja, warto[ zmiennej  a bdzie o  1 wiksza. Na koDcu znajduje si sBowo kluczowe  end , ktre koDczy ptl for. Wyja[niBem ju| na przykBadzie podstawy dziaBania ptli for. Teraz pora na bardziej szczegBowe omwienie. Ptla zawsze wykonuje si okre[lon ilo[ razy. Ilo[ ta zawsze znana jest przed definicj. Kolejne iteracje ptli s zale|ne od speBnienia warunku ilo[ci wykonanych iteracji a nie warunku logicznego, co ma miejsce w przypadku ptli  while". Zmienna, ktra jest licznikiem ptli (w naszym przykBadzie zmienna  a ), w ka|dej iteracji zmienia si o  1 . 5.2 Ptla while Ptla while podobnie jak ptla for sBu|y do wielokrotnego powtarzania kodu. W odr|nieniu od ptli for ilo[ iteracji nie jest znana w momencie pisania kodu. Dla zobrazowania dziaBania ptli posBu|ymy si poprzednim przykBadem. Za pomoc ptli while bdziemy oblicza warto[ funkcji: f(x) = x2 + 12x  3 w tym celu piszemy kod: function f = fun(x) f = x^2 +12*x -3; 27 oraz tworzymy skrypt run2.m n=1; k=5; a=1; while a<=k fun(a) a=a+1; end I uruchamiamy go, wpisujc w oknie Command Window: >> run2 wyniki: ans = 10 ans = 25 ans = 42 ans = 61 ans = 82 Jak wida zastosowanie ptli while w tej postaci daje ten sam efekt jak ptla for. W tym jednak przypadku powinni[my zastosowa ptl while do innych zadaD. Dobrym przykBadem bdzie zastosowanie ptli while do obliczenia po ilu iteracjach warto[ funkcji osignie lub przekroczy oczekiwan warto[. Przyjmijmy, |e warto[ci graniczn w naszym przykBadzie bdzie z=1000. Teraz piszemy kod: function f = fun(x) f = x^2 +12*x -3; tworzymy skrypt run3.m z=1000; a=1; y=0; while y>=z y = fun(a) a=a+1; end a I uruchamiamy go wpisujc w oknie Command Window: >> run3 wyniki: a = 28 28 Otrzymana warto[  28 oznacza, |e funkcja osignBa lub przekroczyBa warto[  1000 po 28 iteracjach, czyli nasza funkcja f(x) dla x=28 przyjmuje warto[: f(28) >= 1000. 5.3 Instrukcja warunkowa if Instrukcje warunkowe, podobnie jak ptle, spotyka si we wszystkich jzykach programowania. Stosuje si je w miejscach, w ktrych chcemy w zale|no[ci od warto[ci jakie[ zmiennej podj jedno z kilku dziaBaD jakie opisali[my w programie. Dobr analogi do instrukcji warunkowej jest skrzy|owanie, na ktrym mo|emy skrci w prawo lub lewo. Najprostszym przykBadem instrukcji warunkowej jest taka konstrukcja: Skrypt przyklad.m if w_ktora_strone == 'w prawo' ulica = 'ul. Zielona' end Nastpnie przechodzimy do Command Window: >> w_ktora_strona = 'w prawo'; >> przyklad ulica = ul. Zielona W przykBadzie stworzyli[my plik przyklad.m , w ktrym umieszczona zostaBa instrukcja warunkowa. Nastpnie po zapisaniu pliku definiujemy zmienn  w_ktora_strona w command window. Po nadaniu jej warto[ci wywoBujemy plik z instrukcj warunkow. W tym przypadku warunek zostaB speBniony i wykonany zostaB kod znajdujcy si wewntrz instrukcji. Oprcz prostych konstrukcji, jak ta wy|ej, mo|na tworzy bardzo rozbudowane wyra|enia z wieloma opcjami. Aby to lepiej zobrazowa zmodyfikujemy nasz przykBad ze skrzy|owaniem. ZaB|my, |e chcemy jako[ uregulowa ruch na skrzy|owaniu czterech ulic. Poniewa| nie chcemy, |eby w ka|d ulic mgB wjecha dowolny samochd musimy postawi jakie[ znaki, ktre pewnym pojazdom zabraniaByby wjazdu na dan ulic. Okre[lili[my, |e bdziemy si kierowa 3 kryteriami - mas pojazdu - wysoko[ci pojazdu - czy wBa[ciciel jest mieszkaDcem osiedla W tym celu tworzymy 3 zmienne: masa [ton] wysokosc [m] wlasciciel [tak/nie] Teraz modyfikujemy nasz plik przyklad.m do takiej postaci: if (masa < 3.5)&(wysokosc < 5) wynik = 'moze jechac w prawo' 29 elseif (masa < 3.5)&(wlasciciel == 'tak') wynik = 'moze jechac w lewo' elseif wysokosc < 3.5 wynik = 'moze jechac prosto' else wynik = 'nie moze jechac w zadnym kierunku' end Nastpnie przechodzimy do okna Command Window: >> masa = 3; >> wysokosc = 3.5; >> wlasciciel = 'tak'; >> przyklad wynik = moze jechac w prawo >> masa = 2.5; >> wysokosc = 1.8; >> wlasciciel = 'nie'; >> przyklad wynik = moze jechac w prawo >> masa = 5; >> wysokosc = 4; >> wlasciciel = 'nie'; >> przyklad wynik = nie moze jechac w zadnym kierunku Jak wida za pomoc bardziej zBo|onych instrukcji warunkowych mo|emy znacznie precyzyjniej okre[la jaki kod ma zosta wykonany w zale|no[ci od grupy zmiennych. Jest to bardzo przydatne przy pisaniu bardziej zBo|onych algorytmw. Dobrym przykBadem jest np. dokBadno[ rwnaD stanu gazw w zale|no[ci od ci[nienia i temperatury panujcych w ukBadzie. Gdyby[my chcieli obliczy zale|no[ ci[nienia od temperatury i objto[ci i narysowa wykres dla bardzo du|ego zakresu ci[nienia, wwczas musieliby[my u|y kilku rwnaD stanu gazu, po to |eby zminimalizowa bBd obliczeD. Tworzymy skrypt gaz.m % deklarujemy zmienne n = 1; % licznosc = 1 mol R = 8.314; % staBa gazowa V = 0.0224; % objeto[ 1 mola w cum warunkach normalnych % wspolczynniki rownania sztywnych kul i van der Waals'a dla metanu a = 0.228; b = 42.8e-6; % tworzymy petle for T = 1:1500 %instrukcja warunkowa - wybieranie rownania stanu gazu if T<500 % rownanie stanu gazu doskonalego p = (n*R*T)/V; elseif (T>=500)&(T<1000) % rownanie sztywnych kul p = (n*R*T)/(V-n*b); elseif T>=1000 % rownanie van der Waals'a p = ((n*R*T)/(V-n*b)) + (a*n^2)/(V^2); 30 end; temperatura(T) = T; cisnienie(T) = p; end; plot(temperatura, cisnienie); W Command Window uruchamiamy nasz skrypt. >> gaz Jako efekt jego dziaBania powinni[my otrzyma nastpujcy wykres: 5 x 10 6 5 4 3 2 1 0 0 500 1000 1500 W kodzie uzale|nili[my wybr rwnania stanu gazu od temperatury panujcej w ukBadzie. W rezultacie otrzymali[my tablic punktw, z ktrej wydrukowali[my wykres. 31 6. Grafika w Matlabie KoDcowym etapem wikszo[ci obliczeD in|ynierskich jest graficzna prezentacja otrzymanych wynikw. Interpretacja niektrych, bardzo skomplikowanych czasami, obliczeD nie byBaby mo|liwa gdyby nie istnienie funkcji graficznych umo|liwiajcych w prosty sposb tworzenie wykresw czy diagramw. Jest to rwnie| jeden z prostszych sposobw na sprawdzenie poprawno[ci obliczeD. Operujc macierzami o kilkunastu czy kilkudziesiciu wymiarach nie jeste[my w stanie stwierdzi, czy otrzymane warto[ci maj sens fizyczny a ich zmiany w czasie nie przecz prawom natury. 6.1 Funkcja plot Funkcja  plot jest podstawowym narzdziem sBu|cym do tworzenia dwuwymiarowych wykresw. Argumentami tej funkcji s cigi liczb, reprezentowanych przez nazw lub wpisywanych wprost. Wpisujc w okno poleceD: >> a=[1.23 3.45 5.98 7.43 9.56]; >> b=[1 3 4 6 8]; >> plot(a,b) w osobnym oknie otworzy si wykres automatycznie nazwany  Figure1 : 8 7 6 5 4 3 2 1 1 2 3 4 5 6 7 8 9 10 Wektory a i b s odpowiednio warto[ciami rzdnych i odcitych. Gdy argumentem jest natomiast tylko jeden wektor, rzdn staje si jego indeks lub, w przypadku liczb zespolonych, OX staje si osi rzeczywist, a OY osi urojon. 6.2 Edycja wykresu Jak wida na powy|szym przykBadzie MATLAB domy[lnie formatuje wykresy w pewien okre[lony 32 sposb, ich wygld mo|na zmieni na dwa sposoby. Pierwszy z nich to edycja otrzymanego ju| wykresu poprzez polecenia w pasku narzdzi  INSERT oraz  TOOLS' EDIT PLOT : wstawianie nazwy osi OX wstawianie nazwy osi OY wstawienie tytuBu wstawienie legendy Jest to sposb na szybkie zatytuBowanie wykresu bez potrzeby wchodzenia w tryb edycji. Tryb edycji umo|liwia midzy innymi zmian: -tytuBu wykresu i osi -grubo[ci, koloru, typu linii i znacznikw -skala i zakres warto[ci na osi -rodzaj i wielko[ czcionki -koloru tBa -siatki -a nawet danych (z dostpnych w pamici-o ile zgadza si wymiar obu wektorw) Zmiany te mog by wprowadzane rwnie| odpowiednimi instrukcjami z okna poleceD, u|ywanie trybu edycji wykresu jest jednak znacznie prostsze dla u|ytkownikw poruszajcych si na co dzieD w [wiecie Windowsa. Drugim sposobem formatowania wykresw jest wprowadzenie dodatkowej zmiennej do funkcji PLOT. Decyduje ona o rodzaju, wielko[ci, kolorze znacznikw oraz linii interpolacyjnej. Zmiany te 33 specyfikuje si u|ywajc odpowiednich oznaczeD wpisywanych w pojedynczych nawiasach: Symbol Kolor Symbol Styl linii r Czerwony - Linia cigBa b Niebieski -- Linia przerywana m R|owy : Linia kropkowana k czarny -. Linia punktowa g Zielony y {Bty W domy[lnych wykresach nie s stosowane znaczniki, co jest korzystne gdy operujemy na du|ej ilo[ci punktw jednak zwykle warto dokBadnie zaznaczy punkty, gdy na jednym wykresie umieszczamy kilka serii danych. Znacznikami mog by:  + ,  " ,  * ,  x ,  . ,  d -romb,  s - kwadrat. Dziki tym oznaczeniom mo|na polepszy wygld powy|szego wykresu, wpisujc w okno poleceD: >> plot(a,b,'*k--')%kolejno[ wpisywania symboli nie ma znaczenia otrzymamy: 8 7 6 5 4 3 2 1 1 2 3 4 5 6 7 8 9 10 Jak wida grubo[ linii jest zawsze taka sama, kolor znacznikw natomiast zawsze odpowiada kolorowi linii. Mo|na to zmieni stosujc nastpujce polecenia: polecenie DziaBanie  markersize Zmiana wielko[ znacznikw  linewidth Zmiana grubo[ci linii  markerfacecolor Zmiana koloru wypeBnienia znacznikw  markeredgecolor Zmiana koloru krawdzi znacznikw 34 Polecenia te wpisuje si w pojedynczych nawiasach, oddzielajc je od innych oznaczeD przecinkiem (warto[ci grubo[ci linii czy wielko[ci znacznikw nale|y jednak wpisa bez nawiasw), np. >> plot(a,b,'*k--', markersize ,20) Dodatkowe komendy pomocne w tworzeniu wykresw: >> title('tytul wykresu') >> xlabel('nazwa osi OX') >> ylabel('nazwa osi OY') >> xlim([2 67])%ustawienie zakresu osi OX% >> ylim('auto')%przywrcenie automatycznego ustawienia zakresu osi OY% Aby jednocze[nie ustawi zakres obu osi mo|na posBu|y si poleceniem  axis , ktrego argumentem jest wektor z zakresem warto[ci dla obu osi: >> axis([2 4 6 9]) Dziki komendom  axis mo|na sterowa skalami osi: polecenie funkcja axis manual obecne zakresy osi s stosowane do nastpnych wykresw tworzonych na tym samym wykresie (z zastosowaniem polecenia  hold on axis auto przywracane jest automatyczne dobieranie skali osi przez program axis tight zakres skali jest [ci[le dostosowany do warto[ci danych na wykresie axis square skala osi zmienia si by pole wykresu byBo kwadratem 6 5 4 3 2 1 0 -1 -2 -3 -10 -5 0 5 10 asix equal skale obu osi staj si takie same 35 8 6 4 2 0 -2 -4 -10 -8 -6 -4 -2 0 2 4 6 8 10 axis image analogicznie do  axis equal skale obu osi s identyczne natomiast pole wykresu jest dostosowywane [ci[le do osiganych 6 warto[ci 4 2 0 -2 -10 -8 -6 -4 -2 0 2 4 6 8 10 W praktyce in|ynierskiej jak wiadomo czsto stosuje si skale logarytmiczne. I w tym przypadku MATLAB przychodzi nam z pomoc. Aby jedna lub obie osie przestawi na skale logarytmiczn wystarczy zastpi polecenie  plot jednym z nastpujcych: polecenie funkcja semilogx skala logarytmiczna dla osi OX semilogy skala logarytmiczna dla osi OY loglog obie osie w skali logarytmicznej 6.3.Hold on Czsto rwnie| na jednym wykresie chcemy zamie[ci kilka krzywych. Mo|na to osign na dwa sposoby. Jeden z nich wymaga u|ycia nowego polecenia:  hold on lub  hold all . Dziki niemu ka|de kolejne u|ycie  plot (lub innego analogicznego polecenia) spowoduje dodanie nowej krzywej do istniejcego ju| wykresu do momentu gdy zastosujemy  hold off , przy czym komenda  hold all zapewnia przydzielanie ka|dej kolejnej krzywej innego koloru. U|ycie  hold on natomiast oznacza identyczne formatowanie ka|dej kolejnej krzywej. >> plot(y,x) >> plot(x,y) >> hold all >> hold on >> plot(x,y) >> plot(y,x) 36 10 10 9 9 8 8 7 7 6 6 5 5 4 4 3 3 2 2 1 1 0 0 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10 6.4 Subplot W momencie gdy serie danych ktre chcemy umie[ci na wykresie r|ni si warto[ciami na tyle, |e jedne z nich s praktycznie niewidoczne warto u|y komendy  subplot dziki ktrej w jednym oknie mo|emy umie[ci kilka wykresw. Argumentami tej funkcji s trzy warto[ci: liczba kolumn, liczba wierszy oraz numer pozycji w ktrej ma si znalez dany wykres: >> subplot(2,2,3)%tworzone s cztery pozycje, a dany wykres znajduje si w trzeciej% Pozycje numerowane s od strony lewej do prawej i od wiersza grnego w dB. Je|eli chcemy dany wykres umie[ci w zajtym ju| miejscu wystarczy do polecenia wpisa  replace : >> subplot(2,2,3,'replace') Po poleceniu  subplot musi wystpi kolejne precyzujce co i jak chcemy wykre[li: >> subplot(2,1,1);fplot(@(x)x^2+x-12,[1 3.5]) >> subplot(2,1,2);fplot(@(x)exp(x)-3.5*x,[0 2.4]) Po powy|szych poleceniach otrzymamy nastepujce dwa wykresy: 5 0 -5 -10 1 1.5 2 2.5 3 3.5 3 2 1 0 -1 0 0.5 1 1.5 2 37 Podobny efekt mo|na uzyska edytujc powstaBy ju| wykres: Dodawanie nowych wykresw do otwartego ju| okna 6.5 Interpolacja Przy niezwykBej mnogo[ci poleceD umo|liwiajcych formatowanie wykresw jednej ich cechy nie mo|na zmieni wprost-charakteru linii interpolacyjnej. Przy du|ej liczbie punktw nie ma to znaczenia gdy| bdzie praktycznie niewidoczna jednak dla kilkunastu punktw wykres mo|e okaza si zbyt  kanciasty . Aby to zmieni nale|y posBu|y si funkcj  interp1 , ktra sBu|y do interpolacji funkcji jednej zmiennej. Przy pracy z funkcjami dwch zmiennych nale|y posBugiwa si  interp2 a dla ambitnych prbujcych dziaBa w przestrzeni trzech zmiennych powstaBo  interp3 . Algorytm postpowania w przypadku ka|dej z powy|szych funkcji jest taki sam. Tworzc wykres, na ktry nanosimy punkty do[wiadczalne automatycznie s one Bczone prost, dlatego przed dodaniem krzywej interpolacyjnej nale|y te prosta usun. >> x=0:10; >> y=atan(x); >> plot(x,y,'o') Kolejnym krokiem jest stworzenie punktw dla danego zakresu argumentw: >> x1=linspace(0,10,200);%aby krzywa interpolacyjna byBa pBynna nale|y w danym przedziale wygenerowa do[ du| liczb argumentw % >> y1=interp1(x,y,x1,'spline'); W efekcie otrzymamy nastpujcy wykres (dla porwnania zamieszczono rwnie| wykres tej samej funkcji bez u|ycia interpolacji): 1.5 1.5 1 1 0.5 0.5 0 0 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10 Argumentami funkcji  interp1 s jak wida w powy|szym przykBadzie: punkty midzy ktrymi 38 ma przebiega krzywa interpolacyjna, wygenerowane argumenty oraz element decydujcy o rodzaju interpolacji. Oprcz zastosowanego w przykBadzie polecenia  spline ( jedna z interpolacji kwadratowych) mo|na wybra rwnie|: polecenie rodzaj krzywej linear interpolacja liniowa-powstaBe odcinki s wielomianem pierwszego stopnia cubic interpolacja kwadratowa-powstaBe odcinki s wielomianami drugiego stopnia nearest metoda  najbli|szego ssiada - wygenerowane argumenty przyjmuj warto[ dla najbli|szego punktu Je|eli podczas interpolacji nie u[ci[limy jej rodzaju wykonane obliczenia bd dotyczy interpolacji liniowej, czyli uzyskany efekt bdzie prawie identyczny jak w przypadku braku jakiejkolwiek naszej ingerencji. 6.6 Funkcja fplot Korzystajc ze znanych ju| narzdzi , aby stworzy wykres zdefiniowanej funkcji nale|y obliczy jej warto[ci dla zadanego zbioru argumentw: >> x=linspace(2,10,5); >> y=exp(x); >> plot(x,y) Proces ten znacznie uBatwia funkcja FPLOT, ktra wykonuje praktycznie caBa prac za nas. Argumentami tej funkcji s: zdefiniowana wcze[niej funkcja oraz zakres dla ktrego chcemy stworzy wykres. Zakres ten mo|e dotyczy zarwno argumentw jak i warto[ci funkcji: >> fplot(@sin,[2 5 -0.6 0.6]) >> fplot(@sin,[2 5]) %wykres funkcji f(x)=sin(x) dla %wykres funkcji f(x)=sin(x) dla argumentw z przedziaBu (2,5) z argumentw z przedziaBu (2,5) ograniczeniem warto[ci funkcji od -0.6 do 0.6% 1 0.5 0.8 0.4 0.6 0.3 0.4 0.2 0.2 0.1 0 0 -0.1 -0.2 -0.2 -0.4 -0.3 -0.6 -0.4 -0.5 -0.8 -1 2 2.5 3 3.5 4 4.5 5 2 2.5 3 3.5 4 4.5 5 Jak wida a przykBadzie program automatycznie dobiera przedziaB warto[ci funkcji (od minimalnej 39 do maksymalnej osiganej przez dana funkcj w wybranym przedziale). Je[li jednak chcemy osign efekt  zoomu nale|y u[ci[li oba przedziaBy. Funkcja u|yta w poleceniu powinna by wcze[niej zdefiniowana jednak mo|liwe jest wypisanie matematycznej formuBy (gdy przykBadowo wiemy, |e nie bdziemy w przyszBo[ci korzysta z danej funkcji): >> fplot(@(x)sin(x)+cos(x)-1,[2 10]) 0.5 0 -0.5 -1 -1.5 -2 -2.5 2 3 4 5 6 7 8 9 10 40 7. CaBkowanie Matlab oferuje wiele funkcji slu|cych do caBkowania numerycznego. Omwimy dwie z nich  ode45  najpopularniejsz i najbardziej uniwersaln funkcj caBkujc dostpn w Matlab-ie oraz funkcj trapz  imlpementujc obliczanie caBek metod trapezw. Na pocztek rozwa|my prost funkcj: y=x2 CaBka nieoznaczona z tej funkcji to oczywi[cie 1 z= x3 3 7.1 Funkcja trapz Zaczniemy od trapz. Funkcj, tak jak wikszo[ funkcji programu Matlab mo|emy wywoBa na wiele sposobw. SzczegBowe informacje i opis funkcji otrzymasz, wpisujc help trapz w oknie Command Window. My bdziemy korzysta z wywoBania, w ktrym jako argumenty przekazujemy wektory warto[ci X i Y. Obliczanie caBki metod trapezw zilustrowano na poni|szym rysunku: Ilustracja caBkowania metod trapezw CaBkowanie t metod polega na przybli|aniu przebiegu funkcji pomidzy kolejnymi znanymi punktami liniami prostymi i obliczeniu sumarycznego pola powstaBych tym sposobem trapezw. Funkcja ta przydaje si gdy mamy np. punkty pomiarowe i musimy obliczy przybli|one pole pod krzyw utworzon z tych punktw. PrzykBad u|ycia funkcji trapz poni|ej: >> x1=[0:1:10]; >> x2=[0:0.1:10]; >> x3=[0:0.01:10]; >> y1=x1.^2; >> y2=x2.^2; >> y3=x3.^2; >> z1=trapz(x1,y1) z1 = 335 >> z2=trapz(x2,y2) z2 = 333.3500 >> z3=trapz(x3,y3) z3 = 333.3335 >> 1/3*(10^3) ans = 333.3333 41 >> Pierwsze trzy operacje to utworzenie wektorw dla zmiennej niezale|nej (x) o rznym zagszczeniu punktw w tym samym przedziale. Nastpnie dla ka|dego z tych wektorw liczymy warto[ci funkcji  y, odpowiadajace wszystkim jego wyrazom i obliczamy caBk przy pomocy funkcji trapz  jej warto[ zapisuejmy w zmiennej z. W ostatniej linii obliczono dokBadn warto[ caBki z funkcji y=x2 w przedziale 0  10. Jak mo|na si byBo spodziewa  dla coraz wikszego zagszczenia punktw otrzymujemy za pomoc funkcji trapz coraz dokBadniejszy wynik. Obliczmy, wykorzystujc funkcj trapz, minimaln ilo[ ciepBa potrzebn do ogrzania 1 kg naftalenu od temperatury 25 st. C (298.15 K)do 80 st. C (353.15 K) w procesie izobarycznym. Szukana ilo[ ciepBa bdzie rwna zmianie entalpii naftalenu, zatem: T=353.15 K Q= H = C T dT +" p T=298.15 K Poka|emy dwa sposoby rozwizania tego zadania  gdy dysponujemy tablic z warto[ciami ciepBa wBa[ciwego naftalenu pod stalym ci[nieniem dla r|nych temperatur oraz gdy dysponujemy korelacj pozwalajac te warto[ci oblicza. Zaczniemy od przypadku w ktrym mamy warto[ci ciepBa wBa[ciwego w r|nych temperaturach zebrane w tabeli. W pierwszm kroku musimy posiadane informacje na temat ciepBa i T [K] C p [J/(kg*K] temperatury wprowadzi do programu Matlab. W tym celu utworzymy tablic i skopiujemy do niej warto[ci temperatury i ciepBa wBa[ciwego 298,15 1570,10 sposobem, omwionym w pierwszym rozdziale. Ilustracja na rysunku 303,15 1586,02 poni|ej: 308,15 1601,93 313,15 1617,85 318,15 1633,76 323,15 1649,68 328,15 1665,59 333,15 1681,51 338,15 1697,42 343,15 1713,34 348,15 1729,25 353,15 1745,17 Wklejanie tablicy temperatura  ciepBo wBa[ciwe Nastpnie wykonujemy caBkowanie przy pomocy funkcji trapz  w ten sposb obliczymy ilo[ energii, ktr nale|y dostarczy. >> trapz(t_cp(:,1),t_cp(:,2)) ans = 9.1170e+04 >> Zatem wymagana energia to okoBo 91.17kJ. Gdy dysponujemy zale|no[ci funkcyjn rwnie| mo|emy u|y funkcji trapz  podobnie jak wtedy, gdy liczyli[my przybli|on warto[ caBki z funkcji y=x2. Korelacja na ciepBo wBa[ciwe pod 42 staBym ci[nieniem dla naftalenu ma posta: J c =621.093.183"T p [ ] kg"K A obliczenia mo|emy wykona w nastpujcy sposb: >> T=[298.15:1:353.15]; >> cp=621.09 + 3.183 * T; >> trapz(T,cp) ans = 9.1170e+04 >> Najpierw tworzymy tablic z warto[ciami temperatury, dla ktrych w kolejnej linii liczymy warto[ci ciepBa wBa[ciwego. Nastpnie korzystamy z funcji trapz aby obliczy warto[ szukanej caBki. 7.2 Funkcja ode45 Najcz[ciej u|ywan funkcj sBu|c do caBkowania w Matlab-ie jest funkcja ode45. Ode to rodzina funkcji caBkujchych, w skBad ktrej wchodzi wiele funkcji nadajcych si do rozwizywania r|nych klas problemw (sztywne ukBady rwnaD r|niczkowych, r|na wymagana dokBadno[ obliczeD). ode45 jest najbardziej uniwersaln funkcj spo[rd nich. Wszystkie funkcje z rodziny ode wywoBuje si w praktycznie identyczny sposb. Informacje na temat innych funkcji z rodziny ode znajdziesz wpisujc w Command Window help ode45 i pod|ajc za odno[nikami do innych funkcji z rodziny ode, znajdujcymi si pod koniec opisu funkcji ode45. Funkcje z rodziny ode pozwalaj nam uzyska przebieg caBki funkcji, mo|emy dziki nim caBkowa ukBady rownaD r|niczkowych a tak|e wszukiwa punkty, dla ktrych caBka osiga okre[lon warto[. Podstawowy sposb wywoBania funkcji ode45 (i innych funkcji z rodziny ode) wyglda nastpujco: [x,y] = ode45(@funkcja_do_scaBkowania,[xp xk],yp) Funkcja zwraca dwa wektory  x i y. Wektor x zawiera warto[ci zmienej niezale|nej a wektor y  odpowiadajce im warto[ci caBki z funkcji funkcja_do_scaBkowania w przedziale od xp do xk. Jako parametr yp przekazujemy warto[ funkcji ktr caBkujemy dla punktu xp. Obliczenie warto[ci caBki z funkcji y=x2 w przedziale 0  10 za pomoc funkcji ode45 obejmuje dwa kroki. W pierwszym - tworzymy plik z nasz funkcj: function dxdy=f1(x,y) dxdy=x.^2; Nastpnie przechodzimy do Command Window i tam obliczamy warto[ caBki za pomoc funkcji 43 ode45: >> [x,y]=ode45(@f1,[0 10],0); >> y(length(y)) ans = 333.3333 >> x(length(x)) ans = 10 >> x(1) ans = 0 Je[li otrzymaBe[ bBd, to prawdopodobnie dlatego, ze plik z funkcj f1 nie znajduje si w jednym z katalogw, w ktrych Matlab poszukuje m-plikw z funkcjami. Problem ten i jego rozwizanie opisano w rozdziale 4.2. y(length(y)) oznacza odwoBanie si do ostatniego elementu wektora y. Analogicznie dla wektora x. Zauwa|, |e jako wynik dostali[my od razu do[ dokBadn warto[ caBki (przynajmniej do 4-go miejsca po przecinku), ale to samo i to nawet pro[ciej, bo bez poleceD typu y(length(y)) mogli[my uzyska przy pomocy funkcji trapz. Aby pozna mo|liwo[ci jakie zyskujemy dziki funkcji ode45, zacznijmy od prostego przykBadu, w ktrym rozwa|ymy model cylindrycznego zbiornika przepBywowego, przedstawiony na poni|szym rysunku: Zbiornik przepBywowy Rwnanie opisujce zmiany ilo[ci cieczy w zbiorniku, przy zaBo|eniu |e ciecz wpBywajca do zbiornika ma tak sam gsto[ jak ciecz wypBywajca ze zbiornika, ma posta: dV =Fin-Fout dt ZakBadamy, |e nat|enie strumienia wpBywajcego do zbiornika jest staBe. Nat|enie strumienia, wypBywajcego ze zbiornika jest proporcjonalne do pierwiastka kwadratowego wysoko[ci cieczy w zbiorniku, zatem: A"h d =Fin-" h dt dh A =Fin-" h dt Fin " dh = - h dt A A 44 Gdzie A to powierzchnia dna zbiornika dana wzorem: "D2 A= 4 Aby uzyska zale|no[ wysoko[ci cieczy w zbiorniku od czasu musimy scaBkowa rwnanie: Fin " dh = - h dt A A Rwnanie to jest na tyle proste, |e bez problemu daje si scaBkowa analitycznie ale my u|yjemy do jego rozwizania funkcji ode45. Dziki niej otrzymamy zale|no[ wysoko[ci cieczy w zbiorniku od czasu w przedziale od t = 0 [s] do t = 400 [s]. ZakBadamy, |e w chwili pocztkowej w p k zbiorniku nie byBo cieczy, czyli: ht =0[s]=0[m] p Tak jak poprzednio najpierw tworzymy m-plik z nasz funkcj obliczajc dh/dt: function dhdt=wysokosc(t,h) Fi=0.07; %m^3/s A=2.5; %m^2 beta=0.08; %m^(5/2)/s dhdt=(Fi-beta.*sqrt(h))./A; i zapisujemy go. Nastpnie u|ywajc funkcji ode45 obliczamy szukan zale|no[ wysoko[ci od czasu a uzyskane warto[ci nanosimy na wykres (funkcj plot): >> [t,h]=ode45(@wysokosc,[0 400],0); >> plot(t,h) Oto uzyskany wykres zmian wysoko[ci cieczy w zbiorniku w czasie. Wykres zmian poziomu cieczy w zbiorniku Gdyby[my chcieli teraz powtrzy obliczenia dla zbiornikw o innych wymiarach, za ka|dym razem musieliby[my edytowa plik z funkcj wysokosc. Znacznie wygodniej byBoby gdyby[my zwyczajnie przekazywali wymiary zbiornika jako argument dla tej funkcji. Funkcja wygldaBa by wtedy tak: 45 function dhdt=wysokosc2(t,h,Fi,A,beta) dhdt=(Fi-beta.*sqrt(h))./A; Aby funkcja ode45 przekazaBa do funkcji ktr ma scaBkowa dodatkowe parametry nale|y j wywoBa w nastpujcy sposb: [x,y] = ode45(@funkcja_do_scaBkowania,[xp xk],yp,[],p1,p2,p3,itd...) Jak wida na koDcu naszego wywoBania dodajemy pusty wektor a za nim kolejne parametry przekazywane do funkcji. W miejsce pustego wektora mo|emy poda dodatkowe opcje dla funkcji rozwizujcej nasze rwnanie  przykBad takiej opcji poka|emy pzniej. CaBkowanie naszej funkcji wysokosc2 wyglda tak: >> [t,h]=ode45(@wysokosc2,[0 400],0,[],0.07,2.5,0.08); >> plot(t,h) >> Powinni[my otrzyma dokBadnie taki sam wykres jak poprzednio. Zapamitaj, |e jako dodatkowy argument mo|esz tak|e przekaza uchwyt funkcji. Mo|e si to czasem przyda  szczeglnie przy pisaniu uniwersalnych skryptw w ktrych bdziesz chciaB mie mo|liwo[ Batwej zmiany na przykBad funkcji obliczajcych parametry mediw biorcych udziaB w procesie. 7.3 Funkcja ode45  caBkowanie ukBadw rwnaD CaBkowanie ukBadw rwnaD r|niczkowych rwnie| mo|na wykona przy pomocy funkcji ode45. Rozwa|my nastpujcy przykBad: W izotermicznym reaktorze zbiornikowym pracujcym w sposb okresowy biegn nastpujce reakcje elementarne: 1: A B 2 :B C 3 :C D StaBe szybko[ci poszczeglnych reakcji w temperaturze procesu wynosz: k = 3.46 . 10-3 [1/s] 1 k = 6.51 . 10-2 [1/s] 2 k = 1.38 . 10-4 [1/s] 2 St|enie A w reaktorze w chwili t=0 wynosi 20 kmol/m3. Nale|y obliczy st|enia poszczeglnych skBadnikw w zale|no[ci od czasu prowadzenia procesu. Zmiany st|enia skBadnikw opisuje ukBad rwnaD r|niczkowych: d CA =-k1"CA dt d CB =k1"C -k2"CB A dt d CC =k2"CB-k3"CC dt d CD =k3"CC dt 46 St|enie skBadnika A w chwili 0 wynosiBo 20 kmol/m3, natomiast pozostaBe skBadniki nie byly obecne w mieszaninie reakcyjnej. Tworzymy plik z funkcj z naszym ukBadem rwnaD: function dcdt=stezenia(t,c) k1=3.46e-3; %1/s k2=6.51e-2; %1/s k3=1.38e-4; %1/s dcdt=[-k1*c(1), k1*c(1)-k2*c(2), k2*c(2)-k3*c(3), k3*c(3)]; Nale|y zwrci uwag, |e w przypadku caBkowania ukBadu rownaD, jako argument bdcy zmienn zale|n (w naszym przypadku c, czyli wektor zawierajacy st|enia skBadnikw) funkcja otrzymuje nie pojedyncza liczb jak to si dziaBo dotychczas ale tablic liczb, zawierajac st|enia poszczeglnych skBadnikw w chwili t. W naszym przypadku pierwszy element tej tablicy  c(1)  odpowiada skBadnikowi A, drugi (c(2))  skBadnikowi B itd. Zauwa|, |e jako warto[ funkcji tak|e zwracany jest wektor, ktrego kolejne pola odnosz si do szybko[ci zmian st|enia poszczeglnych skBadnikw, przy czym kolejno[ skBdnikw jest taka, jak w wektorze c. WywoBanie funkcji ode45 w celu rozwizania ukBadu rownaD r|ni si tym, |e musimy poda wektor warto[ci pocztkowych dla wszystkich zmiennych zale|nych (w naszym przypadku  warto[ci st|eD wszystkich skBadnikw w chwili t=0s, czyli [20 0 0 0]). >> [t,c]=ode45(@stezenia,[0 500],[20 0 0 0]); Nastpnie na wykres nanosimy przebiegi zale|no[ci stzenia od czasu dla poszczeglnych reagentw  po ka|dej komendzie plot spjrz na wykres, |eby zobaczy ktra linia odpowiada ktremu skBadnikowi. >> hold on >> plot(t,c(:,1)) >> plot(t,c(:,2)) >> plot(t,c(:,3)) >> plot(t,c(:,4)) >> 7.4 Zdarzenia (Events) Niekiedy majc dany r|niczkowy model jakiego[ zjawiska (np. model zbiornika z ciecz, jak w jednym z poprzednich podpunktw), chcemy wiedzie dla jakiej warto[ci zmiennej niezale|nej, zmienna zale|na osignie okre[lon warto[ (np. po jakim czasie uzyskamy okre[lony poziom cieczy w zbiorniku). Mo|emy w takim wypadku narysowa wykres i odczyta z niego interesujac nas wartosc z pewn dokBadno[ci, jednak lepiej skorzysta w takiej sytuacji z obecnego w matlabie mechanizmu zdarzeD. Nale|y w tym celu najpierw napisa funkcj, ktra wykryje dla nas wystpienie zdarzenia a nastpnie  poinformowa o jej istnitniu funckj caBkujc  w naszym przypadku ode45. Rozwa|ymy dwa przykBady  znajdziemy czas, po ktrym poziom cieczy w zbiorniku z punktu 7.2 47 wyniesie 0.5 metra oraz czasy, po ktrych stopieD przereagowania skBadnika A w przykBadzie z punktu 7.3 wyniesie 0.5, 0.8 i 0.9. Zacznijmy zatem od zbiornika. Funkcja opisujca szybko[ zmian wysoko[ci cieczy w zbiorniku ma posta: function dhdt=wysokosc(t,h) Fi=0.07; %m^3/s A=2.5; %m^2 beta=0.08; %m^(5/2)/s dhdt=(Fi-beta.*sqrt(h))./A; Musimy teraz stworzy dodatkow funkcj, ktra bdzie odpowiedzialna za wykrycie wystpienia zdarzenia  w naszym przypadku bdzie to moment, w ktrym poziom cieczy w zbiorniku wyniesie 0.5 metra. function [wartosc, koniec, kierunek]=ev_05m(t,h) wartosc = h-0.5; koniec=0; kierunek=0; Funkcja musi zwraca trzy argumenty  pierwszy to warto[. Matlab przyjmuje, ze zdarzenie nastpiBo wtedy, kiedy argument wartosc jest rwny zero  dlatego wBa[nie od wysoko[ci odejmujemy wysoko[ ktra nas interesuje. Kolejny argument okre[la, czy po wystpieniu zdarzenia matlab ma przerwac liczenie caBki (wtedy zmiennej koniec nadajemy warto[ 1) czy ma kontynuowa obliczanie caBki- wtedy tak jak w tym przykBadzie zmiennej nadajemy warto[ 0. Zmienna kierunek pozwala dodatkowo okre[li czy interesuje nas zdarzenie wystpujce tylko wtedy, gdy funkcja jest rosnca (wtedy jej pochodna jest dodatnia a my musimy nada zmiennej kierunek warto[ 1), malejca (-1) lub czy jest to nam obojtne (wtedy tak, jak w naszym przypadku zmiennej kierunek przypisujemy warto[ 0). Ostatnim krokiem przed wywoBaniem funkcji caBkujcej z wykorzystaniem zdarzenia jest utworzenie specjalnego zbioru dodatkowych opcji dla funkcji ode45, przy pomocy ktrego przeka|emy informacj o interesujcym nas zdarzeniu: >> opt=odeset('Events',@ev_05m); odeset tworzy taki wBa[nie zbir opcji. Wszystkie poza tymi, ktre wymienimy w wywoBaniu funkcji (w tym przypadku jest to opcja Events), otrzymuj warto[ci domy[lne. Jako kolejny argument przekazujemy uchwyt naszej funkcji ktra ma sBu|y do wykrycia zdarzenia. Wpisujc w Command Window help odeset dowiesz si jakie inne opcje mo|esz przekazywa do funkcji caBkujcej. Nastpnie wykonujemy caBkowanie: >> [t,h,t_05m,h_05m]=ode45(@wysokosc,[0 400],0,opt); t_05m i h_05m to zmienne, w ktrych znajd si odpowiednio czas, w ktrym wystpiBo zdarzenie i wysoko[ dla ktrej wykryto wystpienie zdarzenia (w naszym przypadku h_05m powinno by rowne 0.5). Sprbuj teraz wykona polecenie: >> plot(t,h,'-',t_05m,h_05m,'o') 48 Aby zobaczy gdzie na wykresie znajduje si znaleziony przy pomocy mechanizmu zdarzeD punkt. Teraz sprbuj samodzielnie zmodyfikowa zdarzenie tak, |eby otrzyma czas, po ktrym poziom cieczy w zbiorniku wyniesie 0.65 m (powinienne[ otrzyma wynik 88.7118 s). Aby lepiej oswoi si z mechanizmem zdarzeD rozwa|ymy jeszcze jeden przykBad  znajdziemy czasy, dla ktrych stopieD przereagowania skladnika A w przykBadzie z rozdziaBu 7.3 wyniesie 0.5, 0.8 i 0.9. Oczwi[cie skorzystamy z napisanej wcze[niej funkcji opisujcej zmiany st|eD skBadnikw w czasie: function dcdt=stezenia(t,c) k1=3.46e-3; %1/s k2=6.51e-2; %1/s k3=1.38e-4; %1/s dcdt=[-k1*c(1), k1*c(1)-k2*c(2), k2*c(2)-k3*c(3), k3*c(3)]; Piszc funkcj wykrywajc zdarzenia, musimy pamita |e w wektorze c mamy st|enia a nie stopieD przereagowania. Musimy zatem przeliczy stopeiD przereagowania na odpowiadajce mu st|enie, zatem: c0-c = c0 std: c=c0"1- St|enie pocztkowe substratu a wynosiBo 20 kmol/m3, zatem : kmol"1-0.5=10 kmol c05=20 m3 m3 kmol"1-0.8=4 kmol c08=20 m3 m3 kmol"1-0.9=2 kmol c09=20 m3 m3 Funkcja wykrywajca zdarzenia bdzie wic miaBa posta: function [wartosc, koniec, kierunek]=evnt(t,c) wartosc = [c(1)-10; c(1)-4; c(1)-2]; koniec=[0;0;0]; kierunek=[0;0;0]; Nastpnie tworzymy struktur z opcjami dla funkcji ode45, do ktrej przekazujemy informacje o naszej funkcji: >> opt=odeset('Events',@evnt); Oraz caBkujemy: >> [t,c,te,ce]=ode45(@stezenia,[0 1000],[20 0 0 0],opt); >> te 49 te = 200.3315 465.1555 665.4871 >> ce ce = 10.0000 0.5613 9.3058 0.1329 4.0000 0.2244 15.1773 0.5983 2.0000 0.1123 16.8434 1.0444 W wektorze ce znajduj si st|enia wszystkich zwizkw w reaktorze w chwili wystpienia zdarzenia (osignicia zadanego st|enia skBadnika A) a w odpowiadajcych im wierszach wektora te czasy, po ktrych osigane s szukane stopnie przereagowania. 50

Wyszukiwarka

Podobne podstrony:
skrypty matlab
MATLAB cw Skrypty
matlab skrypty
20121113 MATLABJ skrypty
Matlab operacje na macierzach, skrypty
Skrypt PODSTAWY PROGRAMOWANIA W JĘZYKU MATLAB
przykład MATLAB skrypt
Matlab Podstawy programowania skrypt
Matlab Podstawy programowania skrypt
8 37 Skrypty w Visual Studio (2)
syst oper skrypty 2
Skrypt Latex
skrypt rozdz 2 4
SIMULINK MATLAB to VHDL Route
IMiR NM2 Introduction to MATLAB
Biochemia zwierząt skrypt UR

więcej podobnych podstron