Matlab: krótki wstęp
Scilab: krótki wstęp
wer. 1.2

Wojciech Myszka

12 stycznia 2008

Spis treści

1 Uruchomienie programu
2 Okno programu w trybie graficznym
3 Matlab (Scilab) jako kalkulator
4 Operacje matematyczne
5 Specjalne funkcje wbudowane
6 Macierze
7 Sposób wyświetlania wyników
8 Tablice
9 Dane i pliki
10 Wykresy
11 Symulacja liniowych obiektów dynamicznych — rozwiązywanie równań różniczkowych
A Słów kilka o „lewym dzieleniu” czyli rozwiązywanie równań liniowych
B Aproksymacja średniokwadratowa
C Kilka przydatnych funkcji
D Ciut o programowaniu

Wstęp

Ten tekst jest podstawowa dokumentacją do programów Matlab (http://www.mathworks.com/) oraz Scilab (http://www.scilab.org/). Programy nie są identyczne, ale na potrzeby początkujących — magą być traktowane zamiennie.

Scilab to darmowy program, rozpowszechniany na licencji, która pozwala na dosyć szerokie jego wykorzystywanie (również do celów komercyjnych).

System Matlab (http://www.mathworks.com/ to zaawansowane narzędzie do rozwiązywania problemów technicznych. Scilab w wielkim stopniu naśladuje Matlaba, a w zakresie objętym tymi zajęciami — może być uznany za jego zastępnik. . .

Nazwa Matlab pochodzi od słów matrix laboratory, które dobrze oddają to w czym Matlab jest najlepszy: szeroko rozumiana algebra liniowa. . .

Niniejsze opracowanie nie powinno być w żadnym przypadku traktowane jako podstawowe źródło informacji. Jego celem jest szybkie przedstawienie podstawowych możliwości Matlaba jako zaawansowanego kalkulatora i narzędzia do obliczeń numerycznych i symulacyjnych.

Wszystkim zainteresowanym pakietem polecam literaturę zamieszczoną na końcu oraz dokumentację pakietu.

1 Uruchomienie programu

Sposób uruchamiania programu zależy od systemu operacyjnego. W Windows najwygodniej znaleźć odpowiednią ikonę i w nią kliknąć. Podobnie można robić w systemach Uniksowych. Tam można również w trybie poleceń wydać polecenie:

tester@piwo2:~$ matlab  
 
                              < M A T L A B >  
                  Copyright 1984-2005 The MathWorks, Inc.  
                   Version 7.0.4.352 (R14) Service Pack 2  
                              January 29, 2005  
 
 
  To get started, type one of these: helpwin, helpdesk, or demo.  
  For product information, visit www.mathworks.com.  
 
>>

Scilab uruchamiany jest w bardzo podobny sposób (albo należy znaleźć odpowiednią ikonkę, albo wydać odpowiednie polecenie.

        -------------------------------------------  
                         scilab-3.1.1  
 
                  Copyright (c) 1989-2005  
              Consortium Scilab (INRIA, ENPC)  
        -------------------------------------------  
 
 
Startup execution:  
  loading initial environment  
 
-->

Aby zakończyć pracę z systemem albo zamykamy okno korzystając z mechanizmów menedżera okien, albo wydajemy polecenie:

>> exit  
tester@piwo2:~$

2 Okno programu w trybie graficznym

Matlab:

PIC

Scilab:

PIC

3 Matlab (Scilab) jako kalkulator

Najbardziej prymitywne zastosowanie programu.

 
>> 2+2  
ans =  
     4  
>> 3/11  
ans =  
    0.2727  
>> 5^3  
ans =  
   125  
>> log(10)  
ans =  
    2.3026  
>> pi  
ans =  
    3.1416  
>> sin(2*pi)  
ans =  
  -2.4493e-16  
>> sin(pi/2)  
ans =  
    1

-->sin(30)  
 ans  =  
  - 0.9880316  
-->sin(%pi/6)  
 ans  =  
    0.5

Uwaga! Scilab nie rozumie (nie zna) zmiennej pi. Zamiast niej rozpoznaje %pi.

Polecenia wpisujemy po znaku „zachęty” >> (-->). Korzystać możemy z praktycznie wszystkich funkcji matematycznych.

Bardziej „zaawansowane” wykorzystanie Matlaba (Scilaba) pozwala na używanie „pamięci” czyli zmiennych:

>> a=%pi  
a =  
    3.1416  
>> b=%pi/2;

Zwracam uwagę na znak średnika umieszczony na końcu wyrażenia: gdy jego nie ma wynik działania wypisywany jest na terminalu. Po jego podaniu — Matlab (Scilab) jest „cichszy”. Nie ma to zbyt wielkiego zastosowania w przypadku wykonywania bezpośrednich obliczeń, gdyż cóż nam z czegoś takiego:

>> 7/8;

ale zawsze możemy się ratować czymś takim:

>> 3/11;  
>> ans  
ans =  
    0.2727  
>>

zmienna ans zawiera zawsze wynik ostatniego wyliczenia wykonywanego w trybie bezpośrednim.

W przypadku zadania do wykonania dziwnych obliczeń Matlab zachowuje się w sposób następujący:

>> 3/0  
Warning: Divide by zero. (Type ”warning off MATLAB:divideByZero” to  
suppress this warning.)  
 
ans =  
   Inf  
>> ans+1  
ans =  
   Inf

A scilab tak:

-->3/0  
    !--error 27  
division by zero...

No tak. . . „Pamiętaj !@#$^&* nie dziel przez zero” jak głosi stare porzekadło — więc napisał, że jest błąd. Ale wynik, w zasadzie, jest dobry: Inf, czyli nieskończoność. Co więcej można na tym wykonywać dalsze działania! + 1 = ! Ale takie coś nie przechodzi:1

>> ans -3/0  
Warning: Divide by zero. (Type ”warning off MATLAB:divideByZero” to  
suppress this warning.)  
 
ans =  
   NaN

NaN to po naszemu Not a Number — „to nie jest liczba”. Zatem odejmowanie dwu nieskończoności nie daje zera. . . 2

Uwaga! W przypadku Scilaba — takich dziwnych obliczeń wykonać nie można. Ale Scilab poprawnie rozpoznaje zmienne o nazwach %inf (nieskończoność), %nan (Not a Number) czy %eps (najmniejsza rozpoznawalna liczba).

4 Operacje matematyczne

+
suma
-
różnica
*
iloczyn
iloraz
\
„lewe” dzielenie:3
    >> 3/2  
    ans =  
        1.5000  
    >> 3\2  
    ans =  
        0.6667  
  

a\b = b∕a

ˆ
potęgowanie

5 Specjalne funkcje wbudowane

Matlab dostarcza kilku specjalnych funkcji. Są to:

%pi
wartość liczby π
%i
jednostka zespolona (√ ---
  - 1)
j
podobnie jak i Tylko Matlab!
%eps
względna dokładność obliczeń zmiennoprzecinkowych (ɛ = 2-52
realmin
najmniejsza (dodatnia) liczba zmiennoprzecinkowa (2-1022) Tylko Matlab!
realmax
największa liczba zmiennoprzecinkowa ((2 - ɛ)21023) Tylko Matlab!
%inf
nieskończoność
%nan
Not a number

Powyższe nazwy nie są zarezerwowane! Można zatem je podmienić wpisując, na przykład:

>> eps=1.e-4  
 
eps =  
   1.0000e-04

Aby „odzyskać” pierwotne znaczenie funkcji należy wydać polecenie clear:

>> clear eps  
>> eps  
 
ans =  
   2.2204e-16

6 Macierze

Aby wprowadzić stałą można napisać tak:

>> [1 2 3; 4 5 6; 7 8 9]  
ans =  
     1     2     3  
     4     5     6  
     7     8     9

Znak ; (średnik) oddziela kolejne wiersze.

Cechą charakterystyczną Matlaba jest to, że tablic nie trzeba deklarować, nie trzeba też podawać rozmiaru, który zajmują. Zmienna może łatwo zmieniać swój charakter:

>> a=2  
a =  
     2  
>> a=[1 2; 0 1]  
a =  
     1     2  
     0     1

Wszystkie działania arytmetyczne uogólniane są tak, by również dotyczyły tablic:

+
suma; jeżeli składniki są tablicami — muszą mieć jednakowy rozmiar; można dodać skalar do tablicy dowolnego rozmiaru:
  >> [1 2; 0 1]+3  
  ans =  
      4     5  
      3     4  
  

-
różnica; uwagi jak wyżej
*
iloczyn macierzy; obowiązują zasady obowiązujące w rachunku macierzowym; jeżeli chcemy wykonać działanie A * B to liczba kolumn A musi być taka sama jak liczba wierszy B; macierz o dowolnych wymiarach może być mnożona przez skalar.
.*
iloczyn „tablicowy” (aby wykonać C = A. * B obie macierze po prawej stronie znaku równości muszą mieć jednakowe wymiary (chyba, że jedna ze zmiennych jest skalarem); Wynikiem jest macierz, której elementy C(i,j)=A(i,j)*B(i,j)
    >> [1 2; 3 4].*[9 9; 9 9]  
    ans =  
        9    18  
       27    36  
  

\
„lewy” iloraz macierzy: jeżeli A i B są macierzami kwadratowymi to A\B = A-1 *B. Gdy A jest macierzą o wymiarach n × n, a B jest wektorem kolumnowym o wymiarze n wówczas X = B\A jest rozwiązaniem układu równań liniowych A * X = B.

W przypadku gdy A ma wymiary m × n gdzie m = n a B jest wektorem kolumnowym o m składowych (lub macierzą o kilku takich kolumnach) to X = B\A jest uogólnionym rozwiązaniem (w sensie najmniejszych kwadratów) układu równań A * X = B

.\
„lewy” iloraz tablicowy: jeżeli macierze A i B mają jednakowe wymiary to wynikiem działania A.\B jest macierz o elementach B(i,j)∕A(i,j). (A albo B może być skalarem!)
iloraz macierzy A∕B = A*B-1. Można również pokazać, że B∕A = (A\B)(gdzie to symbol transpozycji — patrz dalej).
.∕
ˆ
.


Operacja macierzowa Operacja tablicowa


x
1
2
3
y
4
5
6


x
123
y
456


x + y
5
7
9
x - y
-3
-3
-3


x + 2
3
4
5
x - 2
-3
-3
-3


x * y Errorx. * y
4
10
18


x′* y 32x. * y Error


x * y
4 5 6
8 1012
121518
x. * yError


x * 2
2
4
6
x. * 2
2
4
6


x\y 16/7x.\y
4
5/2
2


2\x
1/2
1
3/2
2.∕x
2
1
2/3


x∕y
001/6
001/3
001/2
x.∕y
1/4
2/5
1/2


x∕2
1/2
1
3/2
x.∕2
1/2
1
3/2


xˆy Errorx.ˆy
1
32
729


xˆ2 Errorx.ˆ2
1
4
9


2ˆx Error2.ˆx
2
4
8


(x + i * y)1 - 4i 2 - 5i 3 - 6i


(x + i * y).1 + 4i 2 + 5i 3 + 6i


7 Sposób wyświetlania wyników

Tylko Matlab! To w jaki sposób Matlab wyświetla wyniki można ustalić wydając polecenie format. Należy go użyć podając dodatkowy parametr definiujący ostateczny efekt:

compact
wyłącza drukowanie pustych linii
short
krótki format wydruków:
>> x = [4/3 1.2345e-6]  
x =  
    1.3333    0.0000  
  

short e
krótki format (5 cyfr znaczących; w postaci wykładniczej): 1.3333e+00 1.2345e-06
short g
krótki format (5 cyfr znaczących; w postaci wykładniczej tylko wtedy gdy taka potrzeba): 1.3333 1.2345e-06
long
format długi — pełna precyzja: 1.33333333333333 0.00000123450000
long e
format długi wykładniczy: 1.333333333333333e+00 1.234500000000000e-06
long g
format długi, postać wykładnicza tylko wtedy gdy potrzebna: 1.33333333333333 1.2345e-06
bank
z dokładnością „do groszy” (dwa miejsca po przecinku): 1.33 0.00
rat
format wymierny: 4/3 1/810045; a π wygląda tak:
>> pi  
ans =  
     355/113  
  

hex
szesnastkowy: 3ff5555555555555 3eb4b6231abfd271

8 Tablice

Jak już wspomniano, tablice wpisuje się tak:

 -->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.  !

W ten sposób powstała macierz o nazwie a i trzech wierszach (wiersze oddzielane są za pomocą średników) i czterech kolumnach (elementy wiersza oddzielane są przecinkami).

Do elementów macierzy odwołujemy się w sposób „konwencjonalny”: a(n,m). W przypadku Scilaba — do macierzy dwuwymiarowej można się dobierać tak jak do jednowymiarowej (pamiętając, że dane pamiętane są kolumnami):

-->a(1)  
 ans  =  
    1.  
-->a(2)  
 ans  =  
    5.  
-->a(4)  
 ans  =  
    2.  
-->a(2,2)  
 ans  =  
    6.

Matlab (i Scilab) mają ciekawą właściwość pozwalającą na łatwe wybieranie kolumn (wierszy) macierzy:

--> -->a(1,:)  
 ans  =  
!   1.    2.    3.    4. !  
-->a(:,1)  
 ans  =  
!   1. !  
!   5. !  
!   9. !

Można, z wydobytych dwu kolumn, stworzyć nową macierz:

-->[a(:,1) a(:,3)]  
 ans  =  
!   1.    3.  !  
!   5.    7.  !  
!   9.    11. !

Analogicznie można stworzyć nową macierz z dwu lub więcej wierszy.

Operator : (dwukropek) służy do generowania liczb z podanego zakresu:

-->1:5  
 ans  =  
!   1.    2.    3.    4.    5. !  
 
-->[1:5]*[1:5]’  
 ans  =  
    55.  
 
-->[1:5]’*[1:5]  
 ans  =  
!   1.    2.     3.     4.     5.  !  
!   2.    4.     6.     8.     10. !  
!   3.    6.     9.     12.    15. !  
!   4.    8.     12.    16.    20. !  
!   5.    10.    15.    20.    25. !

W wersji „rozbudowanej” polecenie ma trzy parametry:

-->1:2:10  
 ans  =  
!   1.    3.    5.    7.    9. !

Pierwszy oznacza początek zakresu, drugi krok, a trzeci — koniec. Parametry nie muszą być liczbami całkowitymi:

-->1:%pi/5:2*%pi  
 ans  =  
         column 1 to 5  
!   1.    1.6283185    2.2566371    2.8849556    3.5132741 !  
         column 6 to 9  
!   4.1415927    4.7699112    5.3982297    6.0265482 !

9 Dane i pliki

Podstawowe polecenia służące do zapisywania (i wczytywania) wartości zmiennych to save i load. Polecenie save zapisuje zawartość wszystkich zdefiniowanych zmiennych w pliku:

  save(filename)

w przypadku, gdy chcemy zachować zawartość tylko niektórych zmiennych — ich listę podajemy po nazwie pliku:

save(filename, x1, x2,..., xn)

Składnia polecenia load jest analogiczna:

  load(filename)  
  load(filename, x1, x2,..., xn)

Polecenia load i save mogą służyć do przechowywania danych między sesjami. Tworzone przez polecenie save pliki są binarne!

Polecenia write służy do zapisu danych w postaci „czytelnej” dla człowieka:

-->b=[1:5]’*[1:5]  
 b  =  
 
!   1.    2.     3.     4.     5.  !  
!   2.    4.     6.     8.     10. !  
!   3.    6.     9.     12.    15. !  
!   4.    8.     12.    16.    20. !  
!   5.    10.    15.    20.    25. !  
-->write(”dane1.txt”,b)

zawartość utworzonego pliku dane.txt jest następująca:

  1.  2.  3.  4.  5.  
  2.  4.  6.  8.  10.  
  3.  6.  9.  12.  15.  
  4.  8.  12.  16.  20.  
  5.  10.  15.  20.  25.

Polecenie read ma nieco więcej parametrów:

-->c=read(”dane1.txt”,2,3)  
 c  =  
 
!   1.    2.    3. !  
!   2.    4.    6. !

Pierwszy z nich oznacza liczbę wierszy pliku, którą chcemy przeczytać, a drugi — liczbę kolumn. Gdy zamiast liczby wierszy podamy -1 — przeczytany zostanie cały plik.

-->c=read(”dane1.txt”,3,5)  
 c  =  
!   1.    2.    3.    4.     5.  !  
!   2.    4.    6.    8.     10. !  
!   3.    6.    9.    12.    15. !  
 
-->c=read(”dane1.txt”,-1,5)  
 c  =  
!   1.    2.     3.     4.     5.  !  
!   2.    4.     6.     8.     10. !  
!   3.    6.     9.     12.    15. !  
!   4.    8.     12.    16.    20. !  
!   5.    10.    15.    20.    25. !

Gdy przeczytać chcemy zbyt wielką liczbę danych — występuje błąd:

-->c=read(”dane1.txt”,3,6)  
                        !--error    62  
end of file at line        2

Pod adresem http://kufel.immt.pwr.wroc.pl/~myszka/matlab/sierpien.txt znajduje się przykładowy plik z danymi. Należy go wczytać do Scilaba...

10 Wykresy

Program Scilab ma kilka poleceń służących do rysowania wykresów. Są to:

Wywołanie każdego z tych poleceń bez parametrów — uruchamia małe „demo”:

-->plot  
Demo of plot  
t=0:%pi/20:2*%pi;  
subplot(211)  
plot(t,sin(t),’ro-.’,t,cos(t),’cya+’,t,abs(sin(t)),’--mo’)  
subplot(212)  
plot([t ;t],[sin(t) ;cos(t)],’xdat’,[1:2])

PIC

-->plot2d  
Demo of plot2d  
x=(0:0.1:2*%pi)’;plot2d(x,[sin(x),sin(2*x),sin(3*x)],style=[-1,-2,3],rect=[0,-2,2*%pi,2]);

PIC

Uwaga! kolejne polecenia plot „rysują” wszystko w jednym okienku. Dodatkowe narzędzia pozwalają zarządzać okienkami graficznymi. Są to:

scf(n)
tworzy okno graficzne o numerze n; wywołanie w postaci: f4=scf(4) pozwala „zapamiętać” parametry okna graficznego w celu późniejszego odwołania się do nich.
clf()
kasuje zawartość aktualnego okna graficznego
f4=scf(4); //tworzy okno o identyfikatorze id==4 i staje się ono aktywne  
f0=scf(0); //tworzy okno o identyfikatorze id==0 i staje się ono aktywne  
plot2d()   // tworzy obrazek w oknie o id=0  
scf(f4);   // okno o identyfikatorze id=4 staje sie aktywne  
plot3d()   // tworzy obrazek w oknie o id=4

Uwaga: Dwa znaki dzielenia // oznaczają początek komentarza.

Parametry polecenia plot mogą być następujące:
plot(y,<LineSpec>,<GlobalProperty>)

<LineSpec> — to stała lub zmienna tekstowa zawierająca informacje o rodzaju linii użytej do rysowania wykresu (można użyć znaków -, --, -., :), kolorze (r g b c m y k w) oraz rodzaju markera (znaczka) stawianego w miejscu punktów (o * . x s lub square d lub diamond ^ v < > pentagram; w przypadku braku parametru — przyjęte będzie none czyli brak znacznika.

<GlobalProperty> to nieobowiązkowy parametr pozwalający zdefiniować domyślne ustawienia dla kolejnych wykresów.

W powyższym przypadku rysowany jest wykres wartości umieszczonych w tablicy y w funkcji numeru.

plot(x,y,<LineSpec>,<GlobalProperty>)

W tym przypadku rysowany jest wykres danych z tablicy y w funkcji x.

W przypadku konieczności wyrysowaniu kilku wykresów jednym poleceniem plot — kolejne pary (trójki lub czwórki parametrów) można powtarzać:

  plot(x,y,x1,y1,x2,y2)

Do rysowania wykresów trójwymiarowych służy polecenie plot3d:

-->plot3d  
Demo of plot3d  
t=-%pi:0.3:%pi;plot3d(t,t,sin(t)’*cos(t),35,45,’X@Y@Z’,[2,2,4]);  
 
-->plot3d1  
Demo of plot3d1  
t=-%pi:0.3:%pi;plot3d1(t,t,sin(t)’*cos(t),35,45,’X@Y@Z’,[2,2,4]);

W wyniku powyższych poleceń otrzymamy następujące dwa wykresy:
PIC PIC

Parametry polecenia plot3d (w wersji podstawowej!) to, kolejno:
plot3d(x,y,z,[theta,alpha,leg,flag,ebox])

11 Symulacja liniowych obiektów dynamicznych — rozwiązywanie równań różniczkowych

Część Scilaba służąca do rozwiązywania tych zagadnień bardzo mocno (co do szczegółów, wyglądu innterfejsu) rózni się od Simulinka dostępnego w Matlabie. Natomiast idea — jest dokładnie taka sama.

Ruch wahadła można opisać następującym równaniem różniczkowycm:

LΘ ′′ + pΘ′ + gΘ =  0
(1)

gdzie: Θ — kąt wychylenia wahadła, L — długość wahadła, g — stała grawitacyjna a p — współczynnik tłumienia.

Układ ten przekształcimy do postaci ciut wygodniejszej do dalszych rozważań:

 ′′    -p  ′  -g
Θ  = - L Θ  - L Θ
(2)

Metoda rozwiązywania tego typu równań różniczkowych wywodzi się z technologii (i czasów) maszyn (komputerów) analogowych. Ponieważ trudno znaleźć układ elektroniczny, który służył będzie do różniczkowania (nawet numerycznie ta operacja jest bardzo trudna!), a stosunkowo łatwo daje się całkować zastosowano następujący trik:

PICT

Co oznacza, mniej więcej tyle, że jeżeli na bloczek całkujący wyślemy „sygnał” x to na wyjściu zostanie on „scałkowany” i dostaniemy x. Jeżeli wrzucimy pochodną z Θ (czyli Θ) to na wyjściu, po scałkowaniu — będzie samo Θ. . .

Podążając tym tropem otrzymamy coś takiego:

PICT

Jak coś takiego „zaprogramować” w Scilabie? Trudne pytanie. (Tym bardziej, że opisać trzeba proces klikania myszą. . . Zaczniemy od czegoś prostszego.

Po pierwsze — należy uruchomić ten fragment Scilaba, który odpowiada za modelowanie. Nazywa się on scicos i takie polecenie wpisujemy. Jego efektem jest otwarcie kolejnego okienka (wszystkie zrzuty ekranu pochodzą z wersji Linuxowej).

PIC

Informacje w tym okienku wprowadzamy za pomocą myszy kopiując gotowe bloczki z tak zwanych palet. Aby się do nich dostać z menu Edit wybieramy polecenie Palettes.

PIC

Z listy możemy wybrać jedną paletę. Proponuję zacząć od Sources czyli źródła.

PIC

Klikamy myszą w interesujący nas obiekt (niech to będzie sinusoid generator) i mysz przenosimy do głównego okienka scicosa — do wskaźnika powinien być podczepiony prostokącik. Umieszczamy go gdzieś na planszy (sugeruję wybranie lewej strony) i klikamy. Powinniśmy uzyskać efekt widoczny na kolejnej ilustracji:

PIC

Kolejny raz wybieramy EditPalettes i wybieramy z listy paletę Sinks:

PIC

Z tej palety przenosimy obiekt opisany MScope uzyskując efekt poniższy:

PIC

Łączymy wyjście generatora (taki czarny „trójkącik” po prawej stronie bloczku) z wejściem oscyloskopu (czarny trójkącik po stronie lewej bloczku):

PIC

Z kolejnej palety (Linear) wybieramy bloczek oznaczony symbolem całki.

PIC

PIC

Dokonujemy kolejnych połączeń uzyskując następujący efekt:

PIC

Uwaga! Nie jest to jeszcze nasze wahadło — to tylko prościutki przykład.

Z paletki Sources będzie nam jeszcze potrzebny taki czerwony zegar,4 żeby uzyskać efekt jak poniżej:

PIC

Wybieramy menu SimulationSettings, a w otwartym okienku, w polu Final integration time zmieniamy liczbę 100000 (sekund) na jakąś mniejszą — 15?

PIC

I teraz możemy już z menu Simulation wybrać Run żeby obejrzeć efekt naszej działalności. (Dla leniwych, którzy się pogubili w tym opisie — gotowy model. Trzeba go „wczytać” do okienka scicosa.)

PIC

To co udało nam się uzyskać to generator przebiegu sinusoidalnego, sygnał wyjściowy generatora trafia na układ całkujący. Na „oscyloskopie” możemy oglądać przebieg oryginalny i po scałkowaniu.

W podobny sposób powinniśmy przygotować nasz model wahadła. Potrzebna nam będzie paletka z układmi liniowymi (z której pobierzemy dwa układy do całkowania (integratory) oraz układ sumacyjny); paletka zatytułowana sink (wyjścia?) — oscyloskop oraz paletka z układami nieliniowymi (dwa bloczki pozwalające na realizację dowolnych operacji) i (do obsłużenia wykresów — „czerwony zegar”) z paletki Sources oraz stałą (taki bloczek z dużą jedynką na tej paletce) — trzy razy, aby zdefiniować L — długość wahadła, p — współczynnik tłumienia oraz przyśpieszenie ziemskie (g).

Wszystko łączymy tak aby uzyskać układ następujący:

PIC

Specjalnej uwagi wymagają bloki „Mathematical Expression”. Po pierwsze mają one tylko dwa wejścia (a my potrzebujemy trzy), po drugie — musimy zdefiniować operację, którą mają realizować (x * α∕β). Po dwukrotnym kliknięciu na bloczek otwiera się okienko:

PIC

W polu number of inputs — zamiast dwójki wstawiamy trójkę. W polu scilab expression powinniśmy wpisać coś takiego: –u1*u2/u3. Wejścia bloku są numerowane jako u1, u2,. . . od góry bloku. Jest jeszcze jeden problem: wolelibyśmy mieć wszystkie wejścia bloku po stronie prawej, a jego wyjście — po lewej. W tym celu z menu Edit wybieramy pozycję flip, a następnie kursorem wskazujemy blok, który powinien być „odwrócony”. Cała reszta powinna pójść łatwo. Jeżeli jednak ktoś nie potrafi sobie poradzić — znajdzie gotowy plik.

Podobnie jak w poprzednim przykładzie w menu Simulation wybieramy Setup i skracamy nieco czas symulacji. Tak do 20–30 sekund. A następnie możemy już wybrać Run.

Jeżeli nic się nie dzieje — to trzeba chwilkę się zastanowić. Zwracam uwagę, że model wahadła jest autonomiczny — nie kontaktuje się ze światem zewnętrznym w żaden sposób. Odpowiada to sytuacji rzeczywistego wahadła, które spokojnie wisi sobie na kulce. Do chwili wytrącenia go z pozycji równowagi nie ma prawa się nic dziać.

W jaki sposób wytrącić wahadło z pozycji równowagi: w świecie rzeczywistym wzielibyśmy wahadło w paluszki i odchylili je o niewielki kąt, a następnie puścili i rozpoczęli obserwację. . .

Czyli w chwili „zero” zmienna Θ ma pewną niezerową wartość. Żeby osiągnąć ten efekt — klikamy dwukrotnie na blok całkujący znajdujący się po prawej stronie naszego schematu i w polu initial condition (warunek początkowy) wpisujemy niewielką liczbę. Teraz powinno działać.

PIC

W efekcie uzyskamy wykresy podobne do poniższych:

PIC

Zakończenie

Najbardziej aktualna wersja tego tekstu dostępna jest pod adresem: http://www.immt.pwr.wroc.pl/~myszka/TI/Matlab — w postaci HTML. Dodatkowo pod adresem http://www.immt.pwr.wroc.pl/~myszka/TI/Matlab/matlab.pdf jest wersja PDF tekstu.

Program Scilab można pobrać albo z adresu http://www.scilab.org/download/index_download.php?page=release.html — wszystkie wersje i kod źródłowy.

A Słów kilka o „lewym dzieleniu” czyli rozwiązywanie równań liniowych

Normalna sytuacja to take, gdy mamy n równań i n niewiadomych. Układ równań (liniowych) zapisujemy wówczas zazwyczaj jakoś tak:

Ax  = b

gdzie: A macierz kwadratowa o wymiarach n × n, b — wektor (kolumnowy) o wymierze n, tak zwane „prawe strony”, x — wektor kolumnowy o wymiarze n czyli niewiadome.

Układ taki rozwiązuje się w następujący sposób:

  1. Mnożymy obie strony równania lewostronnie przez A-1.5
      -1       - 1
A   Ax =  A   b

  2. upraszczamy to co się da uprościć i otrzymujemy:
    x = A -1b

Jeżeli mamy równanie następujące:

3x + 4y = 11
2x - 5y =  - 3

wówczas macierze A i b wyglądają tak:

     [       ]       [    ]
       3   4           11
A =    2  - 5    b =    3

W Matlabie/Scilabie będzie to jakoś tak:

-->A=[3 4;2 -5]  
 A  =  
 
!   3.    4. !  
!   2.  - 5. !  
 
-->b=[11;3]  
 b  =  
 
!   11. !  
!   3.  !  
 
-->x=A^(-1)*b  
 x  =  
 
!   2.9130435 !  
!   0.5652174 !

Zamiast A^(-1) można napisać inv(A):

-->y=inv(A)*b  
 y  =  
 
!   2.9130435 !  
!   0.5652174 !

Gdy mamy więcej równań niż niewiadomych — układ równań, zazwyczaj, nie ma jednoznacznego rozwiązania. Podobnie jest w przypadku gdy mamy więcej niewiadomych niż równań.

W takim przypadku (macierz A nie jest kwadratowa!) powyższa metoda nie ma zastosowania.

W przypadku gdy liczba równań jest większa niż liczba niewiadomych (system nadokreślony) — szukamy rozwiązania metodą najmniejszych kwadratów, czyli takiego, które minimalizuje ∣∣Ax - b∣∣.

Gdy niewiadomych (n) jest więcej niż równań (m) (system niedookreślony) — szukamy rozwiązania o co najmniej m niezerowych współrzędnych.

Polecenie \ („lewe dzielenie”) jest tak skonstruowane, że w każdym przypadku x=Aˉ jest rozwiązaniem równania Ax = b.

Druga metoda postępowania to zastosowanie zamiast funkcji inv funkcji pinv. W=pinv(A) generuje macierz W o wymiarach takich jak A(A transponowane) i taką, że:

AW  A =  A,   W AW   = W

jest to pseudoodwrotność macierzy A.

Zastosowanie pseudoodwrotności daje rozwiązanie układu równań, które minimalizuje

∣∣AX  - b∣∣

(∣∣⋅∣∣ to symbol normy).

B Aproksymacja średniokwadratowa

Kolejne zastosowanie operatora „lewego dzielenia” jest następujące: zmierzyliśmy wielkość y w kilku chwilach czasu t uzyskując następujące rezultaty:

t y
0.00.82
0.30.72
0.80.63
1.10.60
1.60.55
2.30.50

Chcemy zamodelować powyższą zależność za pomocą równania:

y(t) = c1 + c2e-t
(3)

Wprowadzamy dane do Scilaba (Matlaba) za pomocą następujących poleceń:

  t = [0 .3 .8 1.1 1.6 2.3]’;  
  y = [.82 .72 .63 .60 .55 .50]’;

Budujemy układ sześciu równań (o dwu niewiadomych), z których mamy wyznaczyć c1 i c2 tak, aby „jak najlepiej” spełnić zależności:

0.82 =  c1 + c2e-0.0
              -0.3
0.72 =  c1 + c2e-0.8
0.63 =  c1 + c2e
0.60 =  c1 + c2e-1.1
0.55 =  c1 + c2e-1.6
0.50 =  c + c e-2.3
        1   2

Poniższe polecenia utworzą macierz E tego równania:

-->E = [ones(t) exp(-t)]  
 E  =  
 
!   1.    1.        !  
!   1.    0.7408182 !  
!   1.    0.4493290 !  
!   1.    0.3328711 !  
!   1.    0.2018965 !  
!   1.    0.1002588 !

(size(t) zwraca długość (czy ogólnie wymiary) wektora t; ones generuje wektor o zadanej długości złożony z jedynek).

Rozwiązaniem problemu jest bardzo proste:

-->c=E\y  
 c  =  
 
!   0.4759522 !  
!   0.3413195 !

Zobaczmy jak to wygląda:

-->Y=[ones(T) exp(-T)]*c;  
 
-->plot(T,Y,’-’,t,y,’o’)

Ciągła linia to wykres funkcji 0.4759522 + 0.3413195e-t, kółeczka oznaczają dane pomiarowe.

PIC

C Kilka przydatnych funkcji

rand(n,m)
generuje tablicę o wymiarach n × m wypełnioną liczbami losowymi (rozkład jednostajny z zakresu [0, 1)
ones(n,m)
generuje tablicę o wymiarach n × m wypełnioną jedynkami
max(A)
wyznacza maksymalny element w tablicy A
min(A)
wyznacza minimalny element w tablicy A
size(B)
informuje wymiarach tablicy B

D Ciut o programowaniu

Pętla
konstrukcja programistyczna znana pod nazwą pętli realizowana jest tak:
for zmienna=wyrażenie, instrukcja, instrukcja, . . . , end
-->n=5;  
 
-->for i = 1:n, for j = 1:n, a(i,j) = 1/(i+j-1);end;end  
 
-->a  
 a  =  
 
!   1.           0.5          0.3333333    0.25         0.2       !  
!   0.5          0.3333333    0.25         0.2          0.1666667 !  
!   0.3333333    0.25         0.2          0.1666667    0.1428571 !  
!   0.25         0.2          0.1666667    0.1428571    0.125     !  
!   0.2          0.1666667    0.1428571    0.125        0.1111111 !  
  

Literatura

[1]   Marek Czajka. MATLAB: ćwiczenia: opanuj środowisko programistyczne MATLAB-a, napisz programy obliczeniowe, zilustruj wyniki obliczeń wykresami. Helion, Gliwice, 2005.

[2]   Aleksandr I. Jastrebov. Optymalizacja – teoria, algorytmy i ich realizacja w MATLAB-ie. Wydawnictwo Politechniki Świętokrzyskiej, Kielce, 2004.

[3]   Anna Kamińska. Ćwiczenia z Matlab: przykłady i zadania. Mikom, Warszawa, 2002.

[4]   Ryszard Lempka. Modelowanie i symulacja układów dynamicznych: wybrane zagadnienia z przykładami w Matlabie. AGH Uczelniane Wydawnictwa Naukowo-Dydaktyczne, Kraków, 2004.

[5]   Bogumiła Mrozek. MATLAB i Simulink: poradnik użytkownika. Helion, Gliwice, 2004.

[6]   Wiesława Regel. Obliczenia symboliczne i numeryczne w programie Matlab. Mikom, Warszawa, 2004.

[7]   Mirosław Wciślik. Wprowadzenie do systemu MATLAB. Wydawnictwo Politechniki Świętokrzyskiej, Kielce, 2000.

[8]   Mirosław Wciślik. Wprowadzenie do systemu MATLAB. Wydawnictwo Politechniki Świętokrzyskiej, Kielce, 2003.

[9]   Andrzej Zalewski. MATLAB – obliczenia numeryczne i ich zastosowania. Nakom, Poznań, 2002.

[10]   Barbara Łysakowska. Komputerowa symulacja układów automatycznej regulacji w środowisku MATLAB/Simulink. Oficyna Wydawnicza Politechniki Wrocławskiej, Wrocław, 2005.