1 Matlab. Wiadomoci wst

,

epne.

1

´

Cwiczenia z symulacji komputerowych. Zaj

,

ecia 4.

Dla dalszych informacji o Matlabie polecam stron

,

e t

,

a.

1

Matlab. Wiadomoci wst

,

epne.

Kod mo˙ze by´c wprowadzany z wiersza polece´

n linia po linii lub uruchamiany z

pliku m, kt´ory jest przez Matlaba interpretowany (nie ma kompilacji czy kodu
po´sredniego).

Linie kodu nie musz

,

a ko´

nczy´c si

,

e ´srednikiem. Je´sli ko´

ncz

,

a si

,

e ´srednikiem,

to wynik operacji nie jest wypisywany. Kilka linii mo˙zna umie´sci´c w jednej,
w´owczas nale˙zy je rozdzieli´c przecinkiem. Przecinek sÃlu˙zy r´ownie˙z do separacji
wsp´oÃlrz

,

ednych w macierzach oraz argument´ow w funkcjach. Z kolei wielokropek

(...) pozwala na przeniesienie cz

,

e´sci instrukcji do kolejnej linii.

Nie ma deklaracji zmiennych liczbowych, od razu si

,

e je definiuje, a domy´slny

typ to double lub matrix. Zmienne typu symbolicznego deklaruje si

,

e poleceniem

syms, np.

syms x f;
f=3*x;
diff(f,x)

Nawiasy zwykÃle oznaczaj

,

a r´ownocze´snie odwoÃlanie do tabeli (wektora), ma-

cierzy jak i wywoÃlanie funkcji, co mo˙ze by´c powodem konfuzji. SÃlu˙z

,

a te˙z do

okre´slania kolejno´sci wykonywania dziaÃla´

n.

W Matlabie nie ma wska´znik´ow ani referencji co praktycznie uniemo˙zliwia

stosowania bardziej skomplikowanych struktur danych.

´

Cwiczenie: Oblicz (symbolicznie) transformat

,

e Laplace’a dla funkcji: 1,

sin(t), t, e

−t

. Skorzystaj z funkcji laplace (informacje o niej: help laplace).

Uwaga: Matlab (etymologia: Matrix laboratory) z zaÃlo˙zenia jest ´srodowiskiem

do oblicze´

n numerycznych raczej ni˙z symbolicznych (programy bardziej zorien-

towane ku obliczeniom symbolicznym to Maple czy Mathematica).

2

Macierze.

Polecenie

a=[3 2; 3 1; 4 5]

stworzy macierz o dw´och wierszach i trzech kolumnach. ´

Srednik oddziela od

siebie kolejne wiersze. Spacja (lub przecinek) rozdziela kolejne pola w wierszu.

Do p´ol macierzy mamy dost

,

ep przez przecinek i nawiasy zwykÃle. Np

a(3,2)

to pole w trzecim wierszu i drugiej kolumnie. W szczeg´olno´sci podstawienie

a(4,5)=6

2 Macierze.

2

automatycznie poszerzy macierz a do macierzy o czterech wierszach i pi

,

eciu

kolumnach i luki uzupeÃlni zerami. Natomiast polecenie

b=a(4,5)

automatycznie zadeklaruje zmienn

,

a b ale tylko je´sli

a(4,5)

ma sens.

Uwaga: wszystkie macierze numerowane s

,

a od 1.

Wektor poziomy zawieraj

,

acy pierwszy wiersz macierzy a mo´zna uzyska´c

przez

b=a(1,:)

natomiast wektor pionowy zawieraj

,

acy jej drug

,

a kolumn

,

e przez

b=a(:,2)

Wektor poziomy z ∈ R

11

taki, ˙ze z

i

= 1 + 0.5(i − 1), i = 1, . . . , 11 mo˙zna

utworzy´c przez

z=1:0.5:6

natomiast i = (1, 2, 3, 4, 5) przez

i=1:5

Operacja prawego dzielenia

x=A\b

rozwi

,

azuje ukÃlad r´owna´

n Ax = b gdzie A jest macierz

,

a kwadratow

,

a, a b wekto-

rem pionowym (czyli macierz

,

a skÃladaj

,

ac

,

a si

,

e z jednej kolumny). Z kolei

x=b/A

rozwi

,

a˙ze nam r´ownwnie xA = b gdzie b musi by´c wektorem wierszowym.

Inne operacje macierzowe to transpozycja (’), mno˙zenie macierzy (*), pot

,

egowanie

(^), dodawanie (+), odejmowanie (-), mno˙zenie wyraz po wyrazie (.*), pot

,

egowanie

wyraz po wyrazie (.^).

Zwr´o´cmy uwag

,

e, ˙ze instrukcja

A+2

zwi

,

ekszy nam o 2 wszystkie pola macierzy lub wektora A.

Inne przydatne funkcje to norm, eig, rank, det (opis ka˙zdej funkcji mo˙zna

przeczyta´c wpisuj

,

ac help nazwaFunkcji).

3 Instrukcje steruj

,

ace.

3

3

Instrukcje steruj

,

ace.

P

,

etla for (przykÃlad)

for i = 1:20
a(i)=i*i;
end

W przykÃladzie tym i przyjmuje kolejno warto´sci kwadrat´ow liczb naturalnych
do 20. W p

,

etli i mo˙ze zmienia´c si

,

e po dowolnym wektorze (wierszowym).

P

,

etla while (przykÃlad)

i=5;
while i>1
a(i)=i-1;
i=i-1;
end

W warunkach logicznych mog

,

a wyst

,

api´c operatory: mniejsze (<), mniejsze lub

r´owne (<=), wi

,

eksze (>), wi

,

eksze lub r´owne (>=), r´owne (==), r´o˙zne (~=), i (&),

lub (|), nie (~). Pojedynczy znak r´owno´sci (=) to instrukcja podstawienia.

Instrukcja if (przykÃlad)

i=2;
if i==3
i=i+1;
elseif i==4
i=i-1;
else
i=0
end

´

Cwiczenie: Napisa´c funkcj

,

e, kt´ora jako argument we´zmie wektor i zwr´oci

wektor zawieraj

,

acy te same elementy w kolejno´sci posortowanej. Skorzystaj z

funkcji length.

4

Funkcje.

SkÃladnia jest nast

,

epuj

,

aca:

function [val1, val2, ...] = nazwa(arg1,arg2,...)

tresc

Je´sli zwracana jest tylko jedna warto´s´c, to nawias kwadratowy nie jest konieczny.
PrzykÃlad:

function b = oJedenWiecej(a)
b=a+1;

5 Solvery MATLABowe dla zada´n pocz

,

atkowych.

4

5

Solvery MATLABowe dla zada´

n pocz

,

atkowych.

Funkcje rozwi

,

azuj

,

ace numerycznie r´ownania zwyczajne to (mi

,

edzy innymi) ode23,

ode45. Funkcja ode23 implementuje metode Rungego i Kutty rz

,

edu 2 z adap-

tacj

,

a czasow

,

a, natomiast ode45 Rungego i Kutty rz

,

edu 4 z adaptacj

,

a czasow

,

a.

Dzialanie jest takie (na przykÃladzie ode23):

ode23(@f,[t1 t2],x0)

f jest tu funkcj

,

a prawej strony r´ownania i powinna mie´c skÃladni

,

e y = f(t,x) -

koniecznie czas jako pierwszy argument. Z kolei [t1 t2] to przedziaÃl czasowy na
kt´orym szukamy rozwi

,

azania u a x0 to warto´s´c warunku pocz

,

atkowego (u(t1)).

Funkcja ode23 zwraca dwie warto´sci: wektor (kolumnowy) punkt´ow czasowych
i wektor kolumnowy warto´sci rozwi

,

azania w tych punktach. PrzykÃlad pliku m

rozwi

,

azuj

,

acego problem u

0

= −2u + t, u(0) = 1:

function [t, y] = test
[t y] = ode23(@f,[0 1],0);

function y=f(t,x)
y=-2*x+t;

Skrypt wywoÃlujemy przez polecenie [a, b] = test.

Solvery ODE dziaÃlaj

,

a dla przypadku wielowymiarowego (ukÃladu r´owna´

n)

analogicznie jak dla jednowymiarowego przy czym x0 powinien by´c teraz wek-
torem wierszowym warunk´ow pocz

,

atkowych, a funkcja f powinna pobiera´c jako

drugi argument wektor (o wymiarze takim samym jak x0) i zwraca´c wektor (te˙z
o takim samym wymiarze) kolumnowy (a nie wierszowy taki jak jest domy´slnie).

PrzykÃlad:

function [t, y] = test

[t y] = ode23(@f,[0 1],[0,0]);

function y=f(t,x)

y(1,1)=-2*x(1)+t+x(2);
y(2,1)=-3*x(2)+t+x(1);

Zwracana warto´s´c y jest macierz

,

a, kt´ora w poszczeg´olnych kolumnach za-

wiera warto´sci rozwi

,

azania dla poszczeg´olnych czas´ow.

Funkcj

,

e wywoÃlujemy przez

[t,y]=test;
plot(t,y(:,1),’-bo’,t,y(:,2),’-ro’)

6

Zadania.

1. Napisa´c funkcj

,

e realizuj

,

ac

,

a caÃlkowanie numeryczne metod

,

a trapez´ow

Z

b

a

f (x) dx ≈

n−1

X

i=0

b − a

2n

(f (a + i

b − a

n

) + f (a + (i + 1)

b − a

n

)).

6 Zadania.

5

Parametrami powinny by´c: funkcja, przedziaÃl caÃlkowania, oraz ilo´s´c krok´ow
siatki. Przetestowa´c dziaÃlanie tej funkcji dla f (x) = sin(x) na przedziale (0,

π

2

).

2. Ci

,

ag logistyczny dany jest rekurencyjnym wzorem a

n

= ka

n−1

(1 − a

n−1

).

Zbadaj eksperymentalnie (za pomoc

,

a wykresu) zachowanie tego ci

,

agu dla a

1

=

0.5 i warto´sci k = 2.8, 3.1, 3.7. Znajd´z maksymaln

,

a warto´s´c k dla kt´orej ci

,

ag

wydaje si

,

e zbiega´c. Znajd´z mimimaln

,

a warto´s´c k dla kt´orej ci

,

ag wydaje si

,

e

zachowywa´c nieprzewidywalnie.

3. Znajd´z warto´sci wÃlasne i odpowiadaj

,

ace im wektory wÃlasne dla macierzy

B=[2 1 0; -1 2 3; -1 3 1]

4. Napisz funkcj

,

e rozwi

,

azuj

,

ac

,

a numerycznie ukÃlad ODE (jest to problem

typu drapie˙znik-ofiara).

u

0

= au − buv − cu

2

,

v

0

= duv − ev

gdzie a = 0.05, b = 0.0002, c = 0.00001, d = 0.0003, e = 0.06. Rozwi

,

a˙z ukÃlad dla

t ∈ (0, 300) oraz warunk´ow pocz

,

atkowych u = 300, v = 100 oraz u = 0, v = 100.

Narysuj w obu przypadkach wykresy obrazuj

,

ace zmiany u i v w czasie.