Z Skrypty dla studentów Ekonofizyki UPGOW
(→Numeryczne rozwiązywanie równania f(x) = 0 metodą Newtona) |
(→Krzywe Lissajous) |
||
(Nie pokazano 258 wersji pomiędzy niniejszymi.) | |||
Linia 2: | Linia 2: | ||
W tym rozdziale umieszczam przykłady krótkich programów, do szczegółowego omówienia i testowania podczas zajęć w laboratorium komputerowym. Programy można zchowywać w m-plikach lub wpisywać bezpośrednio w linii komend. W podanych przykładach pomijam znak zachęty >> | W tym rozdziale umieszczam przykłady krótkich programów, do szczegółowego omówienia i testowania podczas zajęć w laboratorium komputerowym. Programy można zchowywać w m-plikach lub wpisywać bezpośrednio w linii komend. W podanych przykładach pomijam znak zachęty >> | ||
__TOC__ | __TOC__ | ||
+ | |||
+ | ==Obliczanie średniej arytmetycznej== | ||
+ | Obliczanie średniej arytmetycznej n – składnikowego wektora, z wykorzystaniem do tego celu funkcji zdefiniowanej przez programistę. | ||
+ | <source lang="matlab">function wynik = srednia(a) % początek definicji funkcji o nazwie srednia | ||
+ | % oblicza średnią arytmetyczną elementów wektora a | ||
+ | [k, m] = size(a); | ||
+ | if (~((k = = 1) | (m = = 1)) | ||
+ | error('podaj wektor czyli macierz o jednym wierszu albo jednej kolumnie') | ||
+ | wynik = sum(a)/length(a); % sum() i length() to funkcje wbudowane MATLAB’a | ||
+ | end % koniec definicji funkcji o nazwie srednia | ||
+ | </source> | ||
+ | Przykład obliczania średniej, | ||
+ | <source lang="matlab"> | ||
+ | a = 1 : 49; % tworzenie wektora o składowych od 1 do 49 | ||
+ | srednia(a) | ||
+ | ans = | ||
+ | 25 % wynik uśrednienia (1+2+3+...+49)/49 = 25 | ||
+ | </source> | ||
+ | Drugi przykład, | ||
+ | <source lang="matlab"> | ||
+ | a = rand(1, 30); % tworzy wektor 30 liczb przypadkowych z przedz. <0,1> | ||
+ | srednia(a) | ||
+ | ans = % wynik zależy od wartości wylosowanych elementów wektora a</source> | ||
+ | |||
+ | ==Arytmetyka - samouczek== | ||
+ | Program umożliwiający samodzielne sprawdzanie umiejętności wykonywania prostych operacji arytmetycznych. Przeglądnij kod, od strony programistycznej to jest przykład wykorzystania instrukcji wyboru '''switch'''. | ||
+ | <source lang="matlab"> | ||
+ | disp('ARYTMETYKA - SAMOUCZEK' ) | ||
+ | disp(' ') | ||
+ | a = 't' ; | ||
+ | while a == 't' | a == 'T' | ||
+ | |||
+ | disp(' Wybierz dzialanie: ') | ||
+ | disp(' 1. Dodawanie ') | ||
+ | disp(' 2. Odejmowanie ') | ||
+ | disp(' 3. Mnozenie ') | ||
+ | wybor = input(' '); % tutaj decydujesz, jakie działanie cwiczysz | ||
+ | |||
+ | m = floor(10*rand()) + 1; | ||
+ | n = floor(10*rand()) + 1; | ||
+ | |||
+ | switch wybor | ||
+ | case 1 | ||
+ | suma = m + n; | ||
+ | disp(['Ile wynosi suma liczb: ' num2str(m) ' + ' num2str(n) ' ']) | ||
+ | odpowiedz = input(' ') | ||
+ | if odpowiedz == suma | ||
+ | disp(' Brawo ') | ||
+ | else | ||
+ | disp(' Niestety, niedobrze ') | ||
+ | end | ||
+ | |||
+ | case 2 | ||
+ | roznica = m - n; | ||
+ | disp(['Ile wynosi roznica liczb: ' num2str(m) ' - ' num2str(n) ' ']) | ||
+ | odpowiedz = input(' ') | ||
+ | if odpowiedz == roznica | ||
+ | disp(' OK ') | ||
+ | else | ||
+ | disp(' Niestety, .... ') | ||
+ | end | ||
+ | |||
+ | case 3 | ||
+ | iloczyn = m * n; | ||
+ | disp(['Ile wynosi iloczyn iczb: ' num2str(m) ' * ' num2str(n) ' ']) | ||
+ | odpowiedz = input(' ') | ||
+ | if odpowiedz == iloczyn | ||
+ | disp(' Dobrze ') | ||
+ | else | ||
+ | disp(' Odpowiedz niepoprawna ') | ||
+ | end | ||
+ | |||
+ | end | ||
+ | |||
+ | disp( ' ') | ||
+ | a = input(' Czy chcesz kontynuowac ? T/N ', 's') | ||
+ | if a == 'n' | a == 'N' | ||
+ | disp('Bye, bye, bye, ...') | ||
+ | end | ||
+ | disp(' ') | ||
+ | end | ||
+ | </source> | ||
+ | |||
+ | ==Procent prosty== | ||
+ | Procent prosty w bankowości polega na wypłacaniu odsetek od wkładu pieniężnego z dołu, po zakończeniu okresu rozliczeniowego. Korzysta się z następującego wzoru, | ||
+ | |||
+ | :<math>p = w\cdot(1 + s n)\;</math> | ||
+ | |||
+ | gdzie: | ||
+ | ''p'' - przyszła wartość kapitału, | ||
+ | ''w'' - wejściowa wartość kapitału, | ||
+ | ''n'' - czas, np. w latach, | ||
+ | ''s'' - stopa procentowa, np. roczna. Wygodnie stworzyć odpowiednią funkcję, | ||
+ | <source lang="matlab">function [p,z] = procent_prosty(w, n, s) | ||
+ | % oblicza przyszłą wartość kapitału | ||
+ | % p - przyszła wartość kapitału | ||
+ | % w - wejściowa wartość kapitału | ||
+ | % s - stopa procentowa obowiązująca w jednostce czasu, np. w roku | ||
+ | % n - ilość jednostek czasu, w którym obowiązuje stopa s, np. ilość lat | ||
+ | % z - zysk od kapitału o | ||
+ | |||
+ | p = w*(1 + s*n); | ||
+ | z = p - w; | ||
+ | end</source> | ||
+ | Wykonajmy przykładowe obliczenia dla kwoty 10 000 (zł, $, euro,...) "zamrożonej" na okres ''n'' = 5 lat przy stałym oprocentowaniu 5%, czyli ''s'' = 0.05 . | ||
+ | <source lang="matlab"> | ||
+ | >> [kapital, zysk] = procent_prosty(10000, 5, 0.05) | ||
+ | >> [12500, 2500] % no cóż, student ekonofizyki szuka lepszych ofert | ||
+ | </source> | ||
+ | Uwaga: przy obliczaniu oprocentowania banki w Polsce stosują tzw. regułę bankową. Mówi ona, że czas oblicza się dokładnie w dniach, ale zakłada we wzorach na oprocentowanie, że rok ma 360 dni (a nie 365 lub 366). | ||
+ | |||
+ | ==Procent składany== | ||
+ | Procent składany w bankowości to takie oprocentowanie kapitału (wkładu pieniężnego), że odsetki są doliczane do kapitału, procentują wraz z nim w następnym okresie odsetkowym. Kolejne odsetki naliczane są nie tylko od kapitału początkowego, ale również od zarobionych w poprzednim okresie odsetek. | ||
+ | Oznaczając; ''w'' - wejściowa wartość kapitału, ''p'' - przyszła wartość kapitału, | ||
+ | ''m'' - ilość kapitalizacji w roku, ''n'' - ilość lat, ''s''- stopa procentowa wpisana jako ułamek dziesiętny, dostajemy następujący wzór | ||
+ | |||
+ | : <math>p = w\cdot\left(1+\frac{s}{m}\right)^{mn}\;</math> | ||
+ | |||
+ | W szczególnym przypadku, gdy kapitalizacja następuje tylko raz w roku, czyli ''m'' = 1, | ||
+ | |||
+ | : <math>p = w\cdot\left(1+s\right)^n\;</math> | ||
+ | |||
+ | Kapitalizacja ciągła, formalnie następuje nieskończenie wiele razy w roku: | ||
+ | |||
+ | : <math>p = \lim_{m \to \infty} w\cdot(1+\frac{s}{m})^{mn} = w\cdot \left( \lim_{m \to \infty} (1+\frac{s}{m})^m\right) ^n = w\cdot e^{sn}\;</math> | ||
+ | |||
+ | Warto mieć "pod ręką" kod do wyliczania procentu składanego, | ||
+ | |||
+ | <source lang="matlab">function [p,z] = procent_skladany(w, n, m, s) | ||
+ | % oblicza przyszłą wartość kapitału | ||
+ | % p - przyszła wartość kapitału | ||
+ | % w - wejściowa wartość kapitału | ||
+ | % s - stopa procentowa w ustalonym okresie czasu, np. w roku | ||
+ | % m - ilość kapitalizacji w ustalonym okresie czasu, np. w roku | ||
+ | % n - ilość okresów czasu, w którym obowiązuje stopa s, np. ilość lat | ||
+ | % z - zysk od kapitału wejściowego w | ||
+ | |||
+ | p = (1 + s/m)^(n*m) | ||
+ | p = w * p | ||
+ | z = p - w; | ||
+ | end</source> | ||
+ | |||
+ | Obliczenie dla kapitału wejściowego ''w'' = 10 000, trzymanego na koncie przez ''n'' = 5 lat, kapitalizacja miesięczna czyli ''m'' = 12, stopa procentowa 5% tj. ''s'' = 0.05, | ||
+ | <source lang="matlab"> | ||
+ | [kapital, zysk] = procent_skladany(10000, 5, 12, 0.05) | ||
+ | kapital = 1.2833.6e+004 | ||
+ | zysk = 2833.6 | ||
+ | </source> | ||
+ | |||
+ | ==Obliczanie sumy szeregu== | ||
+ | Obliczanie sumy <math> \sum_{k=1}^N \frac{1}{k}+\frac{1}{(k+1)(k+2)} </math> dla wartości N podanej przez użytkownika | ||
+ | <source lang="matlab"> | ||
+ | suma = 0; | ||
+ | N=input(’Podaj N ') | ||
+ | for k = 1: N | ||
+ | suma =suma + 1/k + 1/((k + 2)*(k + 3)); | ||
+ | end | ||
+ | disp( ’ Suma wynosi ’ suma) | ||
+ | </source> | ||
+ | |||
==Szereg geometryczny== | ==Szereg geometryczny== | ||
- | Szereg geometryczny, jest zbieżny dla | q | < 1. Program poniżej oblicza sukcesywnie sumę wzrastającej liczby składników tego szeregu do momentu, gdy przestaje się ona zmieniać, co naturalnie oznacza osiagnięcie zbieżności. | + | Szereg geometryczny, |
- | <source lang="matlab">q = | + | : <math>\sum_{n=0}^\infty q^n=\frac{1}{1-q}</math> |
- | + | ||
+ | |||
+ | jest zbieżny dla | q | < 1. Program poniżej oblicza sukcesywnie sumę wzrastającej liczby składników tego szeregu do momentu, gdy przestaje się ona zmieniać, co naturalnie oznacza osiagnięcie zbieżności. | ||
+ | <source lang="matlab">q = 0.7; % inicjalizacja ilorazu szeregu | ||
+ | suma_szeregu = 0; % szukana suma szeregu, wartość zero na starcie | ||
m = 0; % liczba iteracji (indeks), zero na starcie | m = 0; % liczba iteracji (indeks), zero na starcie | ||
pom = - 1; % zmienna pomocnicza, aby pętla mogła wystartować | pom = - 1; % zmienna pomocnicza, aby pętla mogła wystartować | ||
- | % na początku wartość | + | % na początku nadaję jej wartość -1 |
- | while | + | while suma_szeregu > pom % początek pętli while |
- | pom = | + | pom = suma_szeregu; |
- | + | suma_szeregu = suma_szeregu + q^m; | |
- | + | m = m + 1; | |
end % koniec pętli while | end % koniec pętli while | ||
- | disp( | + | disp(suma_szeregu) % pokazuje na ekranie obliczoną sumę szeregu |
disp(m) % pokazuje końcową liczbę liczbę iteracji | disp(m) % pokazuje końcową liczbę liczbę iteracji | ||
</source> | </source> | ||
Suma szeregu geometrycznego wynosi 2.0 dla q = ½ , możesz zatem łatwo sprawdzić poprawność kodu. Zobacz także zadanie 1 na końcu tego rozdziału. | Suma szeregu geometrycznego wynosi 2.0 dla q = ½ , możesz zatem łatwo sprawdzić poprawność kodu. Zobacz także zadanie 1 na końcu tego rozdziału. | ||
+ | |||
+ | ==Obliczanie pierwiastków równania kawadratowego== | ||
+ | <source lang="matlab">% A x^2 + B*x + C = = 0 równanie kwadratowe | ||
+ | A = 1; , B = 2; , C = 3; % nadanie wartości stałym równania kwadratowego | ||
+ | delta = B^2 - 4*A*C; % obliczanie delta | ||
+ | if delta > 0 | ||
+ | x1 = ( - B - sqrt(delta))/(2*A); , x2 = (-B + sqrt(delta))/(2*A); | ||
+ | elseif delta = = 0 | ||
+ | x1 =- B/(2*A); , x2 = x1; | ||
+ | else | ||
+ | x1=NaN; , x2=NaN; % skrót od ang. not a number | ||
+ | end | ||
+ | disp(‘pierwszy pierwiastek: ‘), disp(x1), disp(‘drugi pierwiastek: ‘), disp(x2)</source> | ||
==Rysowanie okręgu== | ==Rysowanie okręgu== | ||
Linia 60: | Linia 237: | ||
x = t.*cos(t); y = t.*sin(t); | x = t.*cos(t); y = t.*sin(t); | ||
plot(x, y), axis(‘equal’), axis( [-5 5 -5 5] ), axis(‘square’)</source> | plot(x, y), axis(‘equal’), axis( [-5 5 -5 5] ), axis(‘square’)</source> | ||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
==Numeryczne rozwiązywanie równania f(x) = 0 metodą Newtona== | ==Numeryczne rozwiązywanie równania f(x) = 0 metodą Newtona== | ||
Linia 117: | Linia 252: | ||
x<sub>i+1</sub> = x<sub>i</sub> − f(x<sub>i</sub> ) / f `(x<sub>i</sub> ) (3) | x<sub>i+1</sub> = x<sub>i</sub> − f(x<sub>i</sub> ) / f `(x<sub>i</sub> ) (3) | ||
- | Jeżeli funkcja f(x) ma w miarę regularny przebieg oraz x<sub>0</sub> jest blisko prawdziwej wartości pierwiastka badanego równania, to sekwencja {x<sub>0</sub> , x<sub>1</sub> , x<sub>2</sub> , x<sub>3</sub> , . . .} szybko zmierza do x<sub>p</sub> . Napiszmy program, który oblicza kolejne wartości x<sub>i+1</sub> na podstawie równ. 3, aż do N-tej. Liczbę N podaje użytkownik, program po wykonaniu obliczeń wypisuje x<sub>N</sub> | + | Jeżeli funkcja f(x) ma w miarę regularny przebieg oraz x<sub>0</sub> jest blisko prawdziwej wartości pierwiastka badanego równania, to sekwencja {x<sub>0</sub> , x<sub>1</sub> , x<sub>2</sub> , x<sub>3</sub> , . . .} szybko zmierza do x<sub>p</sub> . Napiszmy program, który oblicza kolejne wartości x<sub>i+1</sub> na podstawie równ. 3, aż do N-tej. Liczbę N podaje użytkownik, program po wykonaniu obliczeń wypisuje x<sub>N</sub> , czyli przybliżony pierwiastek równania f(x) = 0 . |
<source lang="matlab">function x = nt(f, pf, x0, N) % początek definicji funkcji nt | <source lang="matlab">function x = nt(f, pf, x0, N) % początek definicji funkcji nt | ||
Linia 168: | Linia 303: | ||
==Metoda bisekcji== | ==Metoda bisekcji== | ||
- | Metoda Newtona numerycznego znajdywania pierwiastka równania f(x) = 0 nie jest jedyna. Alternatywą jest metoda bisekcji (połowienia). Można ją skutecznie zastosować dla ciągłej funkcji f(x) gdy wiemy, w jakim przedziale znajduje się pierwiastek badanego równania. Załóżmy, że jest to przedział | + | Metoda Newtona numerycznego znajdywania pierwiastka równania f(x) = 0 nie jest jedyna. Alternatywą jest na przykład metoda bisekcji (połowienia). Można ją skutecznie zastosować dla ciągłej funkcji f(x) gdy wiemy, w jakim przedziale znajduje się pierwiastek badanego równania. Załóżmy, że jest to przedział x<sub>1</sub> ≤ x ≤ x<sub>2</sub>. Skoro w tym przedziale znajduje się pierwiastek, to f(x<sub>1</sub>)*f(x<sub>2</sub>) < 0. Obliczamy x<sub>3</sub> = (x<sub>1</sub> + x<sub>2</sub>)/2 dzieląc w ten sposób przedział [x<sub>1</sub>, x<sub>2</sub>] na dwie równe części. Wówczas, albo f(x<sub>1</sub>)*f(x<sub>3</sub>) < 0 albo f(x<sub>3</sub>)*f(x<sub>2</sub>) < 0. Jeżeli f(x<sub>1</sub>)*f(x<sub>3</sub>) < 0, to za następny odcinek do podziału wybieramy [x<sub>1</sub>, x<sub>3</sub>] i dzielimy go na pół w punkcie x<sub>4</sub> = (x<sub>1</sub> + x<sub>3</sub>)/2. Jeżeli natomiast f(x<sub>3</sub>)*f(x<sub>2</sub>) < 0, to x<sub>4</sub> = (x<sub>3</sub> + x<sub>2</sub>)/2. Powtarzamy takie podziały aż do momentu, gdy różnica pomiędzy kolejnymi wartościami x<sub>k</sub> będzie mniejsza od założonej dokładności ''delta'' > 0, tj. |x<sub>k+1</sub> – x<sub>k</sub> | < ''delta''. |
<source lang="matlab">function x = polowienie(f, xp, xk, N) | <source lang="matlab">function x = polowienie(f, xp, xk, N) | ||
- | % | + | % polowienie( ) szuka pierwiastka funkcji f(x) w przedziale <xp, xk> |
- | % | + | % f funcja zdefiniowana przez użytkownika (jako inline) |
% xp punkt początkowy badanego przedziału zmiennej x | % xp punkt początkowy badanego przedziału zmiennej x | ||
% xk punkt końcowy badanego przedziału zmiennej x | % xk punkt końcowy badanego przedziału zmiennej x | ||
% N licznba podziałów odcinka w metodzie połowienia | % N licznba podziałów odcinka w metodzie połowienia | ||
- | % x przybliżony pierwiastek równania f(x) = 0 w przedziale | + | % x przybliżony pierwiastek równania f(x) = 0 w przedziale <xp, xk> |
a = f(xp); b = f(xk); | a = f(xp); b = f(xk); | ||
Linia 185: | Linia 320: | ||
x = (xp + xk)/2; | x = (xp + xk)/2; | ||
y = f(x) | y = f(x) | ||
- | if a*y < 0 | + | if a*y < 0 % porównanie znaków f(xp) i f(x) |
- | xk = x; | + | xk = x; % ustalenie punktu końcowego nowego przedziału |
else | else | ||
- | xp = x; | + | xp = x; % ustalenie punktu początkowego nowego przedziału |
end | end | ||
end | end | ||
</source> | </source> | ||
- | Testuj skuteczność metody bisekcji na funkcji f(x) = x<sup>2</sup> – 9 oraz kilku innych funkcjach, | + | Testuj skuteczność metody bisekcji na funkcji f(x) = x<sup>2</sup> – 9 oraz kilku innych funkcjach, pamiętaj, że należy je definiować jako '''inline'''. |
==Symulacje rzutu monetą== | ==Symulacje rzutu monetą== | ||
Linia 198: | Linia 333: | ||
<source lang="matlab">a = rand | <source lang="matlab">a = rand | ||
if (a < 0.5 ) wynik = 1; % umownie ”reszka” | if (a < 0.5 ) wynik = 1; % umownie ”reszka” | ||
- | else wynik = 0; % | + | else wynik = 0; % umownie ”orzełek” |
end | end | ||
</source> | </source> | ||
- | Przypominam, że w MATLAB nie ma zmiennych typu logicznego, w zamian wyrażeniu zawierającemu znaki logiczne przypisuje się wartość 1 , jeżeli jest prawdziwe, lub 0 gdy nieprawdziwe. Kod przedstawiony powyżej można wobec tego zapisać jeszcze prościej | + | Przypominam, że w MATLAB nie ma zmiennych typu logicznego, w zamian zmiennej lub wyrażeniu zawierającemu znaki logiczne przypisuje się wartość 1 , jeżeli jest prawdziwe, lub 0 gdy nieprawdziwe. Kod przedstawiony powyżej można wobec tego zapisać jeszcze prościej |
<source lang="matlab">a = rand | <source lang="matlab">a = rand | ||
wynik = (a < 0.5) % zmienna wynik przyjmie wartość 1 ( ”reszka” ) albo 0 ( ”orzełek” ) | wynik = (a < 0.5) % zmienna wynik przyjmie wartość 1 ( ”reszka” ) albo 0 ( ”orzełek” ) | ||
Linia 219: | Linia 354: | ||
==Rzut kostką== | ==Rzut kostką== | ||
- | Rzut kostką kubiczną ze ścianami ponumerowanymi od 1 do 6. Symulacja rzutu idealną kostką sześcienną oznacza losowanie jednej liczby ze zbioru {1, 2, 3, 4, 5, 6}. Wektor o nazwie a przypadkowych liczb zmiennoprzecinkowych | + | Rzut kostką kubiczną ze ścianami ponumerowanymi od 1 do 6. Symulacja rzutu idealną kostką sześcienną oznacza losowanie jednej liczby ze zbioru {1, 2, 3, 4, 5, 6}. Wektor o nazwie ''a'' zawierający ''n'' przypadkowych liczb zmiennoprzecinkowych z przedziału < 0, 6 > jest generowany przez kod, |
<source lang="matlab">a = 6 * rand(1, n); </source> | <source lang="matlab">a = 6 * rand(1, n); </source> | ||
Interesują nas liczby całkowite, wbudowana funkcja '''fix'''() obcina tą część liczby, która występuje po przecinku (kropce dziesiętnej), wobec tego instrukcja | Interesują nas liczby całkowite, wbudowana funkcja '''fix'''() obcina tą część liczby, która występuje po przecinku (kropce dziesiętnej), wobec tego instrukcja | ||
<source lang="matlab">a =fix( 6 * rand(1, n))+1; </source> | <source lang="matlab">a =fix( 6 * rand(1, n))+1; </source> | ||
- | tworzy wektor o ''n'' składowych, każda z nich jest losowana ze zbioru liczb naturalnych {1, 2, 3, 4, 5, 6}. Wektor a , obliczony na podstawie instrukcji (13.1), przechowuje zatem symulowane wyniki n rzutów idealną kostką sześcienną. Oto wynik, który ja dostałem symulując dziesięć ( n = 10) rzutów kostką | + | tworzy wektor o ''n'' składowych, każda z nich jest losowana ze zbioru liczb naturalnych {1, 2, 3, 4, 5, 6}. Wektor ''a'', obliczony na podstawie instrukcji (13.1), przechowuje zatem symulowane wyniki ''n'' rzutów idealną kostką sześcienną. Oto wynik, który ja dostałem symulując dziesięć ( ''n'' = 10) rzutów kostką |
a = | a = | ||
- | :4 6 3 5 1 2 1 4 1 5 | + | :4 6 3 5 1 2 1 4 1 5 % jakie liczby wylosował Twój komputer ? |
- | Warto przedstawiony kod nieco uogólnić i zapisać w postaci funkcji, która będzie generowała tablicę/macierz składającą się liczb naturalnych losowanych z przedziału | + | Warto przedstawiony kod nieco uogólnić i zapisać w postaci funkcji, która będzie generowała tablicę/macierz składającą się liczb naturalnych losowanych z przedziału <1, m>, |
<source lang="matlab">function N = losowaN(m, r, s) | <source lang="matlab">function N = losowaN(m, r, s) | ||
% tworzy macierz o wymiarach s * t wypełnioną liczbami naturalnymi losowanymi z | % tworzy macierz o wymiarach s * t wypełnioną liczbami naturalnymi losowanymi z | ||
- | % przedziału | + | % przedziału <1, m> |
N = fix(m*rand(r, s)) + 1; | N = fix(m*rand(r, s)) + 1; | ||
end | end | ||
Linia 245: | Linia 380: | ||
Wywołanie funkcji losowaN() w pętli (np '''for''' lub '''while''') wygeneruje serię macierzy wypełnionych liczbami całkowitymi, zob. zadanie 7. | Wywołanie funkcji losowaN() w pętli (np '''for''' lub '''while''') wygeneruje serię macierzy wypełnionych liczbami całkowitymi, zob. zadanie 7. | ||
- | ==Dynamika zysków/strat== | + | ==Dynamika zysków/strat inwestycji== |
- | Podstawą wielu decyzji inwestycyjnych jest ocena dynamiki przyszłych zysków i strat (niestety). Wskutek złożoności procesów ekonomicznych, w przewidywaniach należy uwzględnić składową losową, reprezentującą brak pełnej informacji (niewiedzę). Zobaczmy jak zmieniają się, zarówno zysk jak i strata, wskutek działania czynnika losowego w bardzo prostym modelu | + | Podstawą wielu decyzji inwestycyjnych jest ocena dynamiki przyszłych zysków i strat (niestety). Wskutek złożoności procesów ekonomicznych, w przewidywaniach należy uwzględnić składową losową, reprezentującą brak pełnej informacji (niewiedzę). Zobaczmy jak zmieniają się, zarówno zysk jak i strata, wskutek działania czynnika losowego w bardzo prostym modelu. Wykonujemu n rzutów idealną monetą, kiedy wypada „reszka” zysk wzrasta o np. 1, gdy „orzeł” to tracimy 1. |
<source lang="matlab">function [zysk, strata] = dynamika(n) | <source lang="matlab">function [zysk, strata] = dynamika(n) | ||
% liczy zysk i stratę w zależności od liczby n losowań (rzutów monetą) | % liczy zysk i stratę w zależności od liczby n losowań (rzutów monetą) | ||
- | a = rand(1, n); | + | a = rand(1, n); % losowanie n liczb zmiennoprzecinkowych z przedziału <0, 1> |
- | a = (a<0.5); | + | a = (a<0.5); % jeżeli wartość składowej wektora a jest mniejsza od 0.5, |
- | + | % to zmień ją na 1, wartości pozostałych składowych zmień na 0 | |
- | a = find(a>0); | + | a = find(a>0); % find() zwraca indeksy/numery tych składowych a, które mają |
- | + | % wartość 1 | |
- | zysk = length(a); | + | zysk = length(a); % zwraca rozmiar wektora a, czyli ilość jego składowych. |
- | + | % Ponieważ aktualnie każda składowa ma wartość 1, to zwraca | |
- | + | % ilość jedynek wylosowanych w n rzutach monetą, tego szukamy | |
strata = n - zysk; | strata = n - zysk; | ||
end</source> | end</source> | ||
Linia 265: | Linia 400: | ||
strata = 5025 % !! | strata = 5025 % !! | ||
</source> | </source> | ||
- | Niestety, bilans = zysk - strata = -50 jest ujemny, no po prostu katastrofa !! Cały dzień rzucania monetą ( | + | Niestety, bilans = zysk - strata = -50 jest ujemny, no po prostu katastrofa !! Cały dzień rzucania monetą (10 000 rzutów) i jestem „na minusie”. Może Ty będziesz miał/a więcej szczęścia w inwestycjach, funkcja ''dynamika''( ) jest do dyspozycji, zob. także zadanie 9. |
Naturalnie, sukces lub porażkę inwestycyjną, w starciu z losowymi siłami rynkowymi, możemy prześledzić na wykresie, | Naturalnie, sukces lub porażkę inwestycyjną, w starciu z losowymi siłami rynkowymi, możemy prześledzić na wykresie, | ||
<source lang="matlab">plot(1: n, zysk) | <source lang="matlab">plot(1: n, zysk) | ||
Linia 271: | Linia 406: | ||
ylabel(‘zysk(n)’) | ylabel(‘zysk(n)’) | ||
</source> | </source> | ||
+ | |||
+ | ==Błądzenie przypadkowe== | ||
+ | Rozpatruję najprostszy wariant błądzenia przypadkowego w jednym wymiarze. Cząstka/osobnik ...itp. startuje z punktu ''x'' = 0, w każdym kroku może się przemieszczać o jednostkową odległość w prawo (δx = +1) lub lewo (δx = -1). Prawdopodobieństwo przemieszczenia w prawo jest takie samo, jak w lewo. W jakim położeniu ''x'' znajdzie się cząstka po ''n'' krokach ? Oto kod, | ||
+ | <source lang="matlab"> | ||
+ | function d = pijak(n) | ||
+ | % oblicza położenie cząstki/osobnika po n losowych krokach, startuje od x = 0 | ||
+ | x = rand(1, n); % losowanie n liczb z przedziału <0, 1> | ||
+ | x = (x<0.5); % jeżeli wartość składowej wektora x jest mniejsza od 0.5, | ||
+ | % to zmień ją na 1, wartości pozostałych składowych zmień na 0 | ||
+ | a = find(x>0); % tworzy wektor a którego składowymi są indeksy/numery tych | ||
+ | % składowych x, które mają wartość 1 | ||
+ | np = length(a); % zwraca rozmiar wektora a. Każda składowa wektora a | ||
+ | % reprezentuje (sygnalizuje) obecnośc jedynki w wektorze x. | ||
+ | % Rozmiar wektora a jest wobec tego równy ilości jedynek w | ||
+ | % wektorze x a zatem ilośc kroków w prawo | ||
+ | nl = n - np; % ilośc kroków w lewo | ||
+ | d = np – nl; % odleglość od położenia początkowego x=0 po n losowych krokach, | ||
+ | % gdy d ma wartość ujemną, cząstka/osobnik jest na lewo od x=0 | ||
+ | end | ||
+ | </source> | ||
+ | Wykonaj obliczenia dla różnych ''n''.Nieco ogólniejszy przypadek występuje, gdy prawdopodobieństwo ''p'' kroku w prawo jest inne, niż w lewo (1 -''p''). Zmodyfikowany program, <source lang="matlab"> | ||
+ | function d = asym() | ||
+ | % oblicza położenie cząstki po n losowych krokach, startuje od x = 0 | ||
+ | % prawdopodobieństwo kroku w prawo jest p, w lewo (1-p) | ||
+ | p = input('Podaj p = ') | ||
+ | n = input('Podaj n = ') | ||
+ | x = rand(1, n); | ||
+ | x = (x<p); | ||
+ | x = find(x>0); | ||
+ | np = length(x); | ||
+ | nl = n - np; | ||
+ | d = np – nl; | ||
+ | disp(['p = ' num2str(p)]) | ||
+ | disp(['n = ' num2str(n)]) | ||
+ | disp( ' ') | ||
+ | disp(['x po n krokach = ' num2str(d)]) | ||
+ | end | ||
+ | </source> | ||
+ | Wykonaj obliczenia dla różnych ''p'' przy ustalonym ''n'', np. ''n'' = 100. Należy dodać, iż koncepcja błądzenia przypadkowego pojawia się także w ekonomii, np. przy modelowaniu fluktuacji cen. | ||
+ | |||
+ | ==Metoda Monte Carlo== | ||
+ | Metoda Monte Carlo (nazwę zawdzięcza wykorzystaniu liczb losowych) jest stosowana w różnych dziedzinach nauki i techniki, np. w matematyce do numerycznego obliczania całek. Chcemy wyliczyć całkę <math>\int\limits_a^b f(x) dx</math>. Jeżeli funkcja <math>f(x)</math> jest rzeczywista i ciągła to zachodzi relacja | ||
+ | ::<math>\int\limits_a^b f(x)\ dx = \bar{f} (b-a)</math> , | ||
+ | |||
+ | gdzie <math>\bar{f}</math> oznacza wartość średnią funkcji <math> f(x) </math> w przedziale < a, b >. Obliczanie całki oznaczonej sprowadza się zatem do wyliczenia wartości średniej <math>\bar{f}</math> funkcji <math> f(x) </math>. Wystarczy wylosować (Monte Carlo !) z przedziału < a, b > ''n'' punktów (liczb), obliczyć sumę wartości funkcji <math> f(x) </math> w tych punktach i podzielić przez ''n''. Oto kod programu, funkcję <math> f(x) </math> definiuję jako '''inline''', aby ją łatwo zmieniać i testować program dla rozmaitych postaci funkcji <math> f(x) </math>. | ||
+ | |||
+ | <source lang="matlab"> | ||
+ | f = inline('x.*x', 'x'); | ||
+ | function r = calka_mc(f, a, b, n) | ||
+ | % metodą Monte Carlo liczy całkę oznaczoną funkcji f w granicach od a do b | ||
+ | % n - liczba punktów (liczb) losowanych z przedziału < a, b > | ||
+ | x = (b - a) * rand(1, n); | ||
+ | srednia = sum(f(x))/n; | ||
+ | r = srednia * (b - a); | ||
+ | end</source> | ||
+ | Jako funkcję podcałkową wpisałem <math> f(x) = x^2 </math>, całka nieoznaczona z takiej funkcji wynosi <math> x^3/3 </math>, całka oznaczona w przedziale < 0, 1 > wynosi 1/3, | ||
+ | |||
+ | <source lang="matlab"> | ||
+ | r = calka_mc(f, 0, 1, 10000) | ||
+ | r = 0.3407 | ||
+ | </source> | ||
+ | |||
+ | Zwiększając liczbę ''n'' ("próbkowań" przedziału < 0, 1 > ) powiniśmy otrzymywać rezultat coraz bliższy dokładnej wartości 0.33333... . | ||
+ | |||
+ | ==Obliczanie liczby <math> \pi </math> metodą MC== | ||
+ | Sztandarowy przykład zastosowania metody Monte Carlo to wyznaczenie wartości liczby <math> \pi </math>. Rozważmy koło o promieniu ''r'' wpisane w kwadrat, bok tego kwadratu musi być wobec tego równy 2 ''r''. Pole tego koła wynosi S<sub>o</sub> = <math>\pi r^2 </math>, natomiast pole kwadratu S<sub>k</sub> = <math> (2r)^2 = 4r^2 </math>. Stąd S<sub>o</sub>/S<sub>k</sub> = <math>\pi/4</math>, a zatem <math>\pi</math> = 4 S<sub>o</sub>/S<sub>k</sub>. Wyobraźmy sobie, iż generujemy losowo bardzo dużą ilość punktów N leżących w kwadracie. Część z nich, oznaczę ją N<sub>o</sub>, będzie również leżeć w kole wpisanym w ten kwadrat. Wówczas stosunek S<sub>o</sub>/S<sub>k</sub> możemy zastąpić przez N<sub>o</sub>/N. Stosunek S<sub>o</sub>/S<sub>k</sub> nie zależy od promienia koła, w granicy dla dużych N to samo dotyczy N<sub>o</sub>/N. Wygodnie wobec tego przyjąć, iż ''r'' = 1. | ||
+ | <source lang="matlab"> | ||
+ | function [p] = PI(n) | ||
+ | % oblicza liczbę pi metodą MC | ||
+ | % n - liczba losowanych punktów czyli par liczb {x, y}, przy czym 0≤x<1, 0≤y<1 | ||
+ | nkolo = 0; | ||
+ | for k = 1 : n | ||
+ | x = rand; | ||
+ | y = rand; | ||
+ | r = sqrt(x^2 + y^2); | ||
+ | if (r <= 1) | ||
+ | nkolo = nkolo +1; | ||
+ | end | ||
+ | end | ||
+ | p = 4*nkolo/n; | ||
+ | end | ||
+ | </source> | ||
+ | |||
+ | Oto przykład wywołania wyżej zdefiniowanej funkcji, | ||
+ | <source lang="matlab"> | ||
+ | >> p = PI(1000) | ||
+ | >> p = | ||
+ | 3.1840 % blisko, ale nie idealnie | ||
+ | </source> | ||
+ | |||
+ | Dla celów dydaktycznych podaję drugą wersję tego programu, | ||
+ | <source lang="matlab"> | ||
+ | function p = PI2() | ||
+ | % oblicza Pi metodą MC | ||
+ | nkolo = 0; | ||
+ | rand('state',sum(100*clock)); % zmienia "seed" przy każdym uruchomieniu | ||
+ | n = input('Podaj n = ') % prosi o podanie liczby losowanych punktów | ||
+ | for k = 1: n | ||
+ | x = rand; | ||
+ | y = rand; | ||
+ | r = sqrt(x^2 + y^2); | ||
+ | if (r <= 1) % jeżeli wylosowany punkt jest w kole o promieniu r =1 | ||
+ | nkolo = nkolo + 1; | ||
+ | end | ||
+ | end | ||
+ | p = 4*nkolo/n; | ||
+ | d = abs(pi - p); % błąd bezwzględny | ||
+ | dw = 100*d/pi; % błąd względny podany w procentach | ||
+ | disp(['Pi = ' num2str(p)]) | ||
+ | disp(['delta = ' num2str(d)]) | ||
+ | disp(['delta % = ' num2str(dw)]) | ||
+ | end | ||
+ | </source> | ||
+ | |||
+ | ==Dopasowanie wielomianu (tzw. fitowanie)== | ||
+ | Procedura fitowania modelowej funkcji do istniejącego zbioru danych występuje niezwykle często w naukach ścisłych, technicznych, ekonomii, ...itp. Poniżej przykład kodu źródłowego programu dopasowującego wielomian 5 stopnia do zbioru danych. | ||
+ | <source lang="matlab"> | ||
+ | % program dopasowuje wielomian podanego stopnia do zbioru danych | ||
+ | |||
+ | x = linspace(0, pi, 70); % generowanie punktów na osi x | ||
+ | |||
+ | dane = cos(x) + 0.2*rand(1,length(x)); % utworzenie zbioru danych, w tym | ||
+ | % przykładzie to "zaszumiona" funkcja cos(x) | ||
+ | |||
+ | wsp = polyfit(x,dane,5); % dopasowanie wielomianu do zbioru danych. | ||
+ | % Funkcja standardowa ployfit() zwraca | ||
+ | % wektor współczynników wielomianu n-tego | ||
+ | % stopnia (tutaj n = 5), dopasowanego do | ||
+ | % wartości elementów wektora dane | ||
+ | |||
+ | fit = polyval(wsp,x); % utworzenie wielomianu fit na podstawie | ||
+ | % dopasowanych przez polyfit()współczynników | ||
+ | % (metoda najmniejszych kwadratów)będących | ||
+ | % skladowymi wektora wsp | ||
+ | |||
+ | plot(x, dane, 'r*', x, fit, 'b+') % na wspólnym wykresie dane wyjściowe (dane) | ||
+ | % oraz dopasowany wielomian (fit) | ||
+ | </source> | ||
+ | |||
+ | Zmodyfikowany program, w którym dane są "produkowane" przez wygodną do modyfikacji funkcję typu '''inline''' | ||
+ | <source lang="matlab"> | ||
+ | % fituje wielomian stopnia n do funkcji f(x) w przedziale xp < x < xk | ||
+ | % N - liczba punktów x w przedziale xp < x < xk | ||
+ | % n - stopień dopasowywanego wielomianu | ||
+ | clear; | ||
+ | close all; % zamknięcie okien wcześniejszej grafiki, jeżeli są otwarte | ||
+ | |||
+ | f = inline('(sin(x) + cos(x))', 'x'); % definicja przykładowej funkcji f(x) | ||
+ | |||
+ | function wsp = fit_bis(f, xp, xk, N, n) | ||
+ | x = linspace(xp, xk, N); | ||
+ | wsp = polyfit(x, f(x), n); | ||
+ | p = polyval(wsp, x); | ||
+ | plot(x, f(x), 'r*', x, p, 'b+') | ||
+ | end | ||
+ | </source> | ||
+ | Wywołanie, | ||
+ | <source lang="matlab"> | ||
+ | f = inline('(sin(x) + cos(x))', 'x'); | ||
+ | fit_bis(f, 0, 2*pi, 50, 3) | ||
+ | </source> | ||
+ | daje porównanie na jednym wykresie funkcji f(x) = sin(x) + cos(x) oraz dopasowanego do niej wielomianu 3 stopnia (metoda najmniejszych kwadratów). | ||
+ | [[Plik:fitbis3.gif]] | ||
+ | Testuj jakość dopasowywania (metoda najmniejszych kwadratów) dla różnych postaci funkcji f(x). | ||
+ | |||
+ | ==Krzywe Lissajous== | ||
+ | Rodzinę krzywych Lissajous (fizyka, elektronika) opisują równania parametryczne<br><br><math>x(t)= A\sin(\omega_1 t + \delta), \quad y(t)= B\sin(\omega_2 t)</math> <br><br>Kształt krzywych zależy od stosunku <math>\omega_1</math>/<math>\omega_2</math>. Oto zbiór kilku instrukcji, dających wykres krzywej typu Lissajous | ||
+ | <source lang="matlab"> close all | ||
+ | t = linspace(0, 2*pi, 3000); | ||
+ | x = cos(4*t); | ||
+ | y = sin(7*t); | ||
+ | plot(x, y) | ||
+ | end</source> | ||
+ | :Animację tej krzywej generuje następujący kod, | ||
+ | <source lang="matlab"> close all | ||
+ | t = linspace(0, 2*pi, 2000); | ||
+ | x = cos(4*t); | ||
+ | y = sin(7*t); | ||
+ | comet(x, y) % zobacz, help comet | ||
+ | end</source> | ||
+ | |||
+ | Wykorzystałem tutaj funkcję '''comet'''() (zob. >> help comet), która podobnie jak '''plot'''() tworzy wykres, z tym, że pokazuje proces rysowania w czasie. Zobacz jeszcze animację rysowania okręgu, uruchamiając poniższy kod | ||
+ | <source lang="matlab"> close all | ||
+ | t = 0:0.1*pi:2*pi; % zmieniaj krok czasowy 0.1 na większy oraz mniejszy, co zobaczysz ? | ||
+ | figure; | ||
+ | axis equal; | ||
+ | axis([-1 1 -1 1]); | ||
+ | hold on | ||
+ | comet(cos(t), sin(t)) | ||
+ | end</source> | ||
=Zadania= | =Zadania= | ||
- | # | + | #Program, przedstawiony w przykładzie [[Programowanie_Przykłady#Szereg_geometryczny|Szereg geometryczny]] nie jest zabezpieczony przed beztroskim używaniem. Uczyń go bardziej dojrzałym. Niech pyta użytkownika o podanie ''q'' tj. ilorazu szeregu. Jeżeli poda on wartość ''q'' ≥ 1.0, w odpowiedzi ma pojawić się napis ” |''q''| musi być mniejszy od 1” i program ponownie czeka na podanie prawidłowej liczby. Po pięciokrotnym wpisaniu nieprawidłowej wartości q, ma pojawić się napis „bye, bye, bye..." i nastąpić zakończenie programu. |
- | + | ||
#Zauważ, iż rozwiązaniem równania f(x) = x<sup>2</sup> – 2 = 0 jest x<sub>p</sub> = <math>\sqrt{2}</math>. Oblicz przybliżoną wartość <math>\sqrt{2}</math> wykorzystując do tego celu metodę Newtona (zob. przykłady [[Programowanie_Przykłady#Obliczanie_sumy|8]] i [[Programowanie_Przykłady#Numeryczne rozwiązywanie równania f(x) = 0 metodą Newtona|9]] ). Uwaga, gdyby program się „zapętlił”, to wiesz co masz zarobić, kombinacja klawiszy Ctr – c lub wprowadzić do kodu odpowiednie zabezpieczenie. | #Zauważ, iż rozwiązaniem równania f(x) = x<sup>2</sup> – 2 = 0 jest x<sub>p</sub> = <math>\sqrt{2}</math>. Oblicz przybliżoną wartość <math>\sqrt{2}</math> wykorzystując do tego celu metodę Newtona (zob. przykłady [[Programowanie_Przykłady#Obliczanie_sumy|8]] i [[Programowanie_Przykłady#Numeryczne rozwiązywanie równania f(x) = 0 metodą Newtona|9]] ). Uwaga, gdyby program się „zapętlił”, to wiesz co masz zarobić, kombinacja klawiszy Ctr – c lub wprowadzić do kodu odpowiednie zabezpieczenie. | ||
- | #Utwórz wykres funkcji f(x) = x<sup>3</sup> − 10 w przedziale | + | #Utwórz wykres funkcji f(x) = x<sup>3</sup> − 10 w przedziale <0, 4>. |
- | #Napisz skrypt (m-plik) programu który na jednym wykresie przedstawi funkcje: x<sup>2</sup>, cos(x), ln(2*x), exp(-x) | + | #Napisz skrypt (m-plik) programu, który na jednym wykresie przedstawi funkcje: x<sup>2</sup>, cos(x), ln(2*x), exp(-x) dla zmiennej x należącej do przedziału <0.5, 2>. |
#Wygeneruj wektor o wymiarze ''n'' , ktorego składowymi są liczby losowe o rozkładzie równomiernym z przedziału < 0, 1>. Przedstaw składowe tego wektora, np. dla n = 100 na wykresie. | #Wygeneruj wektor o wymiarze ''n'' , ktorego składowymi są liczby losowe o rozkładzie równomiernym z przedziału < 0, 1>. Przedstaw składowe tego wektora, np. dla n = 100 na wykresie. | ||
- | #Utwórz dwa wektory o ''n'' składowych każdy. | + | #Utwórz dwa wektory o ''n'' składowych każdy. Którąkolwiek ze składowych pierwszego wektora może być liczba 0 albo 1, reprezentująca wynik rzutu idealną monetą (1 ”reszka”, 0 ”orzełek”). Dokonaj w komputerze symulacji ''n'' rzutów idealną monetą, wpisując wyniki jako składowe pierwszego wektora. Przedstaw ilość zer i jedynek w tym wektorze na wykresie, dla różnych ''n'', od 10 do 500. Składowymi drugiego wektora niech będa zera albo jedynki, losowane dla monety sfałszowanej, pokazującej reszkę średnio 80 razy na 100 rzutów.<br>Na wykresie porównaj ilości zer i jedynek w pierwszym i drugim wektorze dla ''n'' od 10 do 500. Zastosuj różne kolory dla zwiększenia czytelności. |
#Wygeneruj 50 macierzy o wymiarach 4 * 4 każda, wypełnionych liczbami naturalnymi losowanymi z przedziału <0, 10> . Oblicz macierz będącą średnią z sumy tych macierzy, następnie odrzuć mantysy wszystkich elementów tej macierzy „średniej”, aby dostać macierzy „średnią całkowitą” tj wypełnioną elementami całkowitymi. | #Wygeneruj 50 macierzy o wymiarach 4 * 4 każda, wypełnionych liczbami naturalnymi losowanymi z przedziału <0, 10> . Oblicz macierz będącą średnią z sumy tych macierzy, następnie odrzuć mantysy wszystkich elementów tej macierzy „średniej”, aby dostać macierzy „średnią całkowitą” tj wypełnioną elementami całkowitymi. | ||
- | #Napisz program tworzący wektor | + | #Utwórz wektor ''a'' zawierający 20 liczb całkowitych. Utwórz wektor ''b'', którego składowymi są liczby parzyste występujące w ''a''. Następnie utwórz wektor ''c'', zawierający te składowe wektora ''a'', które są liczbami nieparzystymi. |
- | #Utwórz funkcję, która liczy zysk i stratę po n rzutach specjalnie spreparowaną monetą, dla której prawdopodobieństwo pokazania „reszki” jest p, orzełka naturalnie 1 - p ( 0 < p < 1). Wykonaj symulacje rzutów tej monety dla wielu wartości n , np. n = 100: 100: | + | #Utwórz macierz ''A'' o sześciu wierszach i ośmiu kolumnach, której elementami będą liczby całkowite wybrane losowo z przedziału [−12, 12]. |
+ | #Napisz program tworzący wektor o sześciu składowych, który gromadzi liczby naturalne losowane z przedziału < 1, 49 >, przy czy losowane liczby nie mogą się powtarzać. Tak, to symulacja tzw. ciągnienia w loterii typu LOTTO. | ||
+ | #Utwórz funkcję, która liczy zysk i stratę po ''n'' rzutach specjalnie spreparowaną monetą, dla której prawdopodobieństwo pokazania „reszki” jest ''p'', orzełka naturalnie 1 - ''p'' ( 0 < ''p'' < 1). Wykonaj symulacje rzutów tej monety dla wielu wartości ''n'' , np. ''n'' = 100: 100: 1000. Porównaj na jednym wykresie przebiegi bilansów (bilans = zysk – strata) w funkcji ''n''. Obserwuj na wykresie zależność zysku od parametru ''p'' („stopnia nierzetelności”) monety przy ustalonym ''n'', np ''n'' = 1000. Przeczytaj przykład [[Programowanie_Przykłady#Dynamika zysków/strat inwestycji|Dynamika zysków/strat inwestycji]]. | ||
+ | #Narysuj wykres wielomianu f(x) = x<sup>4</sup> − 4x<sup>2</sup> - 3, używając do tego celu 100 punktów z przedziału −3 < x < 3. | ||
+ | #Pokaż na wspólnym wykresie przebiegi następujących funkcji f(t): exp(−2t)cos(ωt), exp(−t)[4cos(ωt) + sin(2ωt)] w przedziale 0 ≤ t ≤ 15, dla trzech wartości ω = 2, 5, 10. | ||
+ | #Testuj skuteczność obliczania całek oznaczonych metodą Monte Carlo. Wykonaj całkowania numeryczne dla następujących funkcji f(x): cos(x), sin(x), x<sup>3</sup>sin(x)cos(x), exp(-x)sin<sup>2</sup>(x)cos<sup>3</sup>(4x<sup>2</sup>) dla ''x'' z przedziału <0, 2<math>\pi</math> >. Pokaż na wykresie, jak zmienia się numerycznie obliczana wartość całki <math>\int\limits_0^\pi sin((x) dx</math> przy rosnącej liczbie ''n'' próbkowań przedziału < 0, <math>\pi</math> >, np dla ''n'' = 10, 100, 1000, 10 000, 100 000, 1 000 000. | ||
+ | #Wykonaj dopasowanie wielomianu ''n''-tego stopnia do funkcji f(x) = 2x<sup>2</sup>cos(x)sin(2x), ''n'' = 2, 3, 4. Zobacz na wykresie, jak zmiana ''n'' wpływa na precyzję dopasowania. | ||
+ | #Dopasuj wielomiany do następujących funkcji f(x): 3x<sup>3</sup>cos<sup>2</sup>(x)sin<sup>2</sup>(4x), x<sup>4</sup>exp(x<sup>2</sup>)sin(x), xexp(x<sup>3</sup>)cos(5x). Przedstaw na jednym wykresie wyniki dopasowania wielomianiu 5-tego stopnia (''n'' = 5)do tych funkcji f(x). | ||
+ | #Rodzinę krzywych Lissajous opisuje układ równań <br><math>x(t)= A\sin(\omega_1 t + \delta)</math> oraz <math> y(t)= B\sin(\omega_2 t).</math> <br> Kształt krzywych zależy od stosunku <math>\omega_1</math>/<math>\omega_2</math>.Wykreśl na jednym wykresie dwie krzywe dla następujących wartości parametrów: a) A = B = 1 oraz <math>\delta = 0 </math>, b) A = B = 1 oraz <math>\delta </math> = <math>\pi</math>/2. Podstaw za <math>\omega_1 </math> i <math>\omega_2</math> kilka różnych liczby naturalnych i pokaż na wspólnym wykresie odpowiednie krzywe Lissajous dla 0 < t ≤ 2 pi. Jak wzrost całkowitych wartości <math>\omega</math> wpływa na kształt tych krzywych. Zobacz co się stanie, gdy zwiększysz zakres ''t'' np. 20 razy. Czy krzywe się zamykają ? Obserwuj krzywe Lissajous w przedziale 0 < t ≤ 2<math>\pi</math>, gdy <math>\omega</math> jest liczbą niewymierną, np. <math>\omega</math> = e, <math>\pi</math>, <math>\sqrt{3} </math> … Zwiększ znacznie zakres ''t'', czy krzywe się zamykają ? Po wykonaniu tych wykresów, być może potrafisz uzupełnić stwierdzenie; krzywe Lissajous są zamknięte, gdy stosunek <math>\omega_1</math>/<math>\omega_2</math> jest ..... , no właśnie ??? | ||
+ | #Jeżeli dość malownicza rodzina krzywych Lissajous wzbudziła Twoje zainteresowanie, uruchom kod, który tworzy wykres krzywej typu Lissajous | ||
+ | <source lang="matlab">close all | ||
+ | t = linspace(0, 2*pi, 3000); | ||
+ | x = cos(4*t); | ||
+ | y = sin(7*t); | ||
+ | plot(x, y) | ||
+ | end</source> | ||
+ | :Animację tej krzywej generuje następujący kod, | ||
+ | <source lang="matlab">close all | ||
+ | t = linspace(0, 2*pi, 2000); | ||
+ | x = cos(6*t); | ||
+ | y = sin(3*t); | ||
+ | comet(x, y) | ||
+ | end</source> | ||
+ | :Poczytaj (help comet) o funkcji '''comet'''(), eksperymentuj z animacją krzywych Lissajous lub innych. |
Aktualna wersja na dzień 11:38, 8 lut 2011
Przykłady
W tym rozdziale umieszczam przykłady krótkich programów, do szczegółowego omówienia i testowania podczas zajęć w laboratorium komputerowym. Programy można zchowywać w m-plikach lub wpisywać bezpośrednio w linii komend. W podanych przykładach pomijam znak zachęty >>
Obliczanie średniej arytmetycznej
Obliczanie średniej arytmetycznej n – składnikowego wektora, z wykorzystaniem do tego celu funkcji zdefiniowanej przez programistę.
function wynik = srednia(a) % początek definicji funkcji o nazwie srednia % oblicza średnią arytmetyczną elementów wektora a [k, m] = size(a); if (~((k = = 1) | (m = = 1)) error('podaj wektor czyli macierz o jednym wierszu albo jednej kolumnie') wynik = sum(a)/length(a); % sum() i length() to funkcje wbudowane MATLAB’a end % koniec definicji funkcji o nazwie srednia
Przykład obliczania średniej,
a = 1 : 49; % tworzenie wektora o składowych od 1 do 49 srednia(a) ans = 25 % wynik uśrednienia (1+2+3+...+49)/49 = 25
Drugi przykład,
a = rand(1, 30); % tworzy wektor 30 liczb przypadkowych z przedz. <0,1> srednia(a) ans = % wynik zależy od wartości wylosowanych elementów wektora a
Arytmetyka - samouczek
Program umożliwiający samodzielne sprawdzanie umiejętności wykonywania prostych operacji arytmetycznych. Przeglądnij kod, od strony programistycznej to jest przykład wykorzystania instrukcji wyboru switch.
disp('ARYTMETYKA - SAMOUCZEK' ) disp(' ') a = 't' ; while a == 't' | a == 'T' disp(' Wybierz dzialanie: ') disp(' 1. Dodawanie ') disp(' 2. Odejmowanie ') disp(' 3. Mnozenie ') wybor = input(' '); % tutaj decydujesz, jakie działanie cwiczysz m = floor(10*rand()) + 1; n = floor(10*rand()) + 1; switch wybor case 1 suma = m + n; disp(['Ile wynosi suma liczb: ' num2str(m) ' + ' num2str(n) ' ']) odpowiedz = input(' ') if odpowiedz == suma disp(' Brawo ') else disp(' Niestety, niedobrze ') end case 2 roznica = m - n; disp(['Ile wynosi roznica liczb: ' num2str(m) ' - ' num2str(n) ' ']) odpowiedz = input(' ') if odpowiedz == roznica disp(' OK ') else disp(' Niestety, .... ') end case 3 iloczyn = m * n; disp(['Ile wynosi iloczyn iczb: ' num2str(m) ' * ' num2str(n) ' ']) odpowiedz = input(' ') if odpowiedz == iloczyn disp(' Dobrze ') else disp(' Odpowiedz niepoprawna ') end end disp( ' ') a = input(' Czy chcesz kontynuowac ? T/N ', 's') if a == 'n' | a == 'N' disp('Bye, bye, bye, ...') end disp(' ') end
Procent prosty
Procent prosty w bankowości polega na wypłacaniu odsetek od wkładu pieniężnego z dołu, po zakończeniu okresu rozliczeniowego. Korzysta się z następującego wzoru,
\[p = w\cdot(1 + s n)\;\]
gdzie: p - przyszła wartość kapitału, w - wejściowa wartość kapitału, n - czas, np. w latach, s - stopa procentowa, np. roczna. Wygodnie stworzyć odpowiednią funkcję,
function [p,z] = procent_prosty(w, n, s) % oblicza przyszłą wartość kapitału % p - przyszła wartość kapitału % w - wejściowa wartość kapitału % s - stopa procentowa obowiązująca w jednostce czasu, np. w roku % n - ilość jednostek czasu, w którym obowiązuje stopa s, np. ilość lat % z - zysk od kapitału o p = w*(1 + s*n); z = p - w; end
Wykonajmy przykładowe obliczenia dla kwoty 10 000 (zł, $, euro,...) "zamrożonej" na okres n = 5 lat przy stałym oprocentowaniu 5%, czyli s = 0.05 .
>> [kapital, zysk] = procent_prosty(10000, 5, 0.05) >> [12500, 2500] % no cóż, student ekonofizyki szuka lepszych ofert
Uwaga: przy obliczaniu oprocentowania banki w Polsce stosują tzw. regułę bankową. Mówi ona, że czas oblicza się dokładnie w dniach, ale zakłada we wzorach na oprocentowanie, że rok ma 360 dni (a nie 365 lub 366).
Procent składany
Procent składany w bankowości to takie oprocentowanie kapitału (wkładu pieniężnego), że odsetki są doliczane do kapitału, procentują wraz z nim w następnym okresie odsetkowym. Kolejne odsetki naliczane są nie tylko od kapitału początkowego, ale również od zarobionych w poprzednim okresie odsetek. Oznaczając; w - wejściowa wartość kapitału, p - przyszła wartość kapitału, m - ilość kapitalizacji w roku, n - ilość lat, s- stopa procentowa wpisana jako ułamek dziesiętny, dostajemy następujący wzór
- \(p = w\cdot\left(1+\frac{s}{m}\right)^{mn}\;\)
W szczególnym przypadku, gdy kapitalizacja następuje tylko raz w roku, czyli m = 1,
- \(p = w\cdot\left(1+s\right)^n\;\)
Kapitalizacja ciągła, formalnie następuje nieskończenie wiele razy w roku:
- \(p = \lim_{m \to \infty} w\cdot(1+\frac{s}{m})^{mn} = w\cdot \left( \lim_{m \to \infty} (1+\frac{s}{m})^m\right) ^n = w\cdot e^{sn}\;\)
Warto mieć "pod ręką" kod do wyliczania procentu składanego,
function [p,z] = procent_skladany(w, n, m, s) % oblicza przyszłą wartość kapitału % p - przyszła wartość kapitału % w - wejściowa wartość kapitału % s - stopa procentowa w ustalonym okresie czasu, np. w roku % m - ilość kapitalizacji w ustalonym okresie czasu, np. w roku % n - ilość okresów czasu, w którym obowiązuje stopa s, np. ilość lat % z - zysk od kapitału wejściowego w p = (1 + s/m)^(n*m) p = w * p z = p - w; end
Obliczenie dla kapitału wejściowego w = 10 000, trzymanego na koncie przez n = 5 lat, kapitalizacja miesięczna czyli m = 12, stopa procentowa 5% tj. s = 0.05,
[kapital, zysk] = procent_skladany(10000, 5, 12, 0.05) kapital = 1.2833.6e+004 zysk = 2833.6
Obliczanie sumy szeregu
Obliczanie sumy \( \sum_{k=1}^N \frac{1}{k}+\frac{1}{(k+1)(k+2)} \) dla wartości N podanej przez użytkownika
suma = 0; N=input(’Podaj N ') for k = 1: N suma =suma + 1/k + 1/((k + 2)*(k + 3)); end disp( ’ Suma wynosi ’ suma)
Szereg geometryczny
Szereg geometryczny,
- \(\sum_{n=0}^\infty q^n=\frac{1}{1-q}\)
jest zbieżny dla | q | < 1. Program poniżej oblicza sukcesywnie sumę wzrastającej liczby składników tego szeregu do momentu, gdy przestaje się ona zmieniać, co naturalnie oznacza osiagnięcie zbieżności.
q = 0.7; % inicjalizacja ilorazu szeregu suma_szeregu = 0; % szukana suma szeregu, wartość zero na starcie m = 0; % liczba iteracji (indeks), zero na starcie pom = - 1; % zmienna pomocnicza, aby pętla mogła wystartować % na początku nadaję jej wartość -1 while suma_szeregu > pom % początek pętli while pom = suma_szeregu; suma_szeregu = suma_szeregu + q^m; m = m + 1; end % koniec pętli while disp(suma_szeregu) % pokazuje na ekranie obliczoną sumę szeregu disp(m) % pokazuje końcową liczbę liczbę iteracji
Suma szeregu geometrycznego wynosi 2.0 dla q = ½ , możesz zatem łatwo sprawdzić poprawność kodu. Zobacz także zadanie 1 na końcu tego rozdziału.
Obliczanie pierwiastków równania kawadratowego
% A x^2 + B*x + C = = 0 równanie kwadratowe A = 1; , B = 2; , C = 3; % nadanie wartości stałym równania kwadratowego delta = B^2 - 4*A*C; % obliczanie delta if delta > 0 x1 = ( - B - sqrt(delta))/(2*A); , x2 = (-B + sqrt(delta))/(2*A); elseif delta = = 0 x1 =- B/(2*A); , x2 = x1; else x1=NaN; , x2=NaN; % skrót od ang. not a number end disp(‘pierwszy pierwiastek: ‘), disp(x1), disp(‘drugi pierwiastek: ‘), disp(x2)
Rysowanie okręgu
Rysowanie okręgu o środku w początku kartezjańskiego układu współrzędnych i jednostkowym promieniu r = 1.
alfa = 0 : pi/60 : 2*pi; y = sin(alfa); x = cos(alfa); plot(x,y), axis(’square’) % bez tej instrukcji powstałaby elipsa wskutek różnego skalowania % osi poziomej i pionowej, zob. help plot oraz help axis control
Wykorzystanie liczb zespolonych do rysowania okręgu
Zobacz także poprzedni przykład
alfa = 0 : pi/60 : 2*pi; z = exp(i*alfa); plot(z), axis(’square’)
Rysowanie kilku funkcji na jednym wykresie
Następujący kod powoduje przedstawienie na jednym wykresie grupy funkcji f(x) = x, x3, x5, x7
x = 0: .2: 2; f1 = x; f2 = x.^3; f3 = x.^5; f4 = x.^7; plot(x, f1, x, f2, x, f3, x, f4)
Uwaga: Jeżeli ostatnią linię zastąpimy przez
plot(x, f1, ‘black’, x, f2, ‘blue’, x, f3, ‘green’, x, f4, ‘red’)
wtedy każda z czterech funkcji będzie rysowana innym kolorem.
Zbiorowy wykres czterech kolejnych nieparzystych funkcji potęgowych otrzymamy na podstawie alternatywnego kodu źródłowego,
x = 0: .2: 2; f = [x; x.^3; x.^5: x.^7]; plot(x, f)
Wykres parametryczny
Wykres parametryczny krzywej opisanej dwoma równaniami: x (t) = t cos(t), y(t) = t sin(t) dla 0 ≤ t ≤ 2 π .
t = 0: .02: 2*pi; x = t.*cos(t); y = t.*sin(t); plot(x, y), axis(‘equal’), axis( [-5 5 -5 5] ), axis(‘square’)
Numeryczne rozwiązywanie równania f(x) = 0 metodą Newtona
Punktem wyjścia jest liniowe przybliżenie rozwinięcia w szereg funkcji f(x) dla wybranego punktu x0 ,
f(x) ≈ f(x0 ) + f ` (x0 ) (x − x0 ) (1)
gdzie f `(x0 ) jest pochodną funkcji dla x = x0 . Ponieważ chcemy znaleźć takie x, dla którego f(x) = 0, więc za lewą stronę powyższego równania podstawiamy zero i rozwiązujemy ze względu na x:
x ≈ x0 − f(x0 ) / f `(x0 ) (2)
Dostaliśmy przepis na znajdywanie przybliżonej wartości x, dla której f(x) = 0. Zaczynamy procedurę od wytypowania wartości startowej x0 , która według dostępnych przesłanek wydaje się być bliska pierwiastkowi xp badanego równania. Następnie generujemy sekwencję wartości {x0 , x1 , x2 , x3 , . . .} według reguły:
xi+1 = xi − f(xi ) / f `(xi ) (3)
Jeżeli funkcja f(x) ma w miarę regularny przebieg oraz x0 jest blisko prawdziwej wartości pierwiastka badanego równania, to sekwencja {x0 , x1 , x2 , x3 , . . .} szybko zmierza do xp . Napiszmy program, który oblicza kolejne wartości xi+1 na podstawie równ. 3, aż do N-tej. Liczbę N podaje użytkownik, program po wykonaniu obliczeń wypisuje xN , czyli przybliżony pierwiastek równania f(x) = 0 .
function x = nt(f, pf, x0, N) % początek definicji funkcji nt % oblicza pierwiastek równania f(x) = 0 , metoda Newtona % f zadana funkcja, podać w formie inline % pf pochodna funkcji f, też podać jako inline % x0 liczba startowa % N ilość kroków (liczba naturalna) x = x0; for m = 1 : N x = x – f(x)/pf(x) % zakodowana formuła (3) end % koniec definicji funkcji nt
Zapisujemy ten kod w pliku nt.m w bieżącym katalogu roboczym. W linii komend definiujemy w formie inline funkcję, np. f(x) = x4 – 16
f = inline(‘ x^4 - 16 ‘)
oraz jej pochodną,
pf = inline(‘ 4 * x^3 ‘)
Ostatni krok to wywolanie fukcji nt(f, pf, x0, N) dla konkretnych wartości parametrów x0 oraz N, na przykład
xp = nt(f, pf, 1.0, 20); disp(xp)
Dobrą praktyką jest wielokrotne wykonanie obliczeń, dla różnych wartości N. Jeżeli od pewnej wartości liczby N, przy dalszym jej zwiększaniu nie obserwuje się już znaczącej zmiany xN, można domniemać, że znaleźliśmy numerycznie rozwiązanie równania f(x) = 0. Naturalnie, w razie potrzeby należy także manipulować wartościami parametru startowego x0. W naszym przykładzie, rozwiązanie równania f(x) = x4 – 16 = 0 jest oczywiste, xp = 2 lub -2. Sprawdź, czy metoda Newtona prowadzi do poprawnego wyniku, jeżeli tak, to ile iteracji wystarcza ? Testuj tą metodę, definiując bardziej skomplikowane funkcje f(x).
Modyfikacja poprzedniego przykładu, nie podajemy z góry liczby iteracji N. Wprowadzamy w zamian małą liczbę epsilon > 0. Chcemy zakończyć iteracje dopiero wtedy, gdy |f(xm) – f(xp) | < epsilon. Ponieważ xp jest pierwiastkiem badanego równania (czyli f(xp) = 0), komputer ma przerwać iteracje, gdy |f(xm)| < epsilon. Oto zmodyfikowany kod
function x = ntmod(f, pf, x0, epsilon) % oblicza pierwiastek równania f(x) = 0 , metoda Newtona % f zadana funkcja, podać w formie inline % pf pochodna funkcji f, też podać jako inline % x0 liczba startowa % epsilon mała liczba rzeczywista dodatnia x = x0; f2 = f(x); while abs(f2) > epsilon x = x – f2/pf(x) f2 = f(x); end
Dalej należy wywołać funkcję ntmod( ) dla wybranych parametrów, np.
xp = nt(f, pf, 1.0, .001); disp(xp)
Wykonaj obliczenia również dla epsilon = 0.1, 0.01, 0.00001. Uwaga, praktyka pokazuje, iż metoda Newtona bywa nieskuteczna, gdy x0 leży zbyt daleko od xp.
Metoda bisekcji
Metoda Newtona numerycznego znajdywania pierwiastka równania f(x) = 0 nie jest jedyna. Alternatywą jest na przykład metoda bisekcji (połowienia). Można ją skutecznie zastosować dla ciągłej funkcji f(x) gdy wiemy, w jakim przedziale znajduje się pierwiastek badanego równania. Załóżmy, że jest to przedział x1 ≤ x ≤ x2. Skoro w tym przedziale znajduje się pierwiastek, to f(x1)*f(x2) < 0. Obliczamy x3 = (x1 + x2)/2 dzieląc w ten sposób przedział [x1, x2] na dwie równe części. Wówczas, albo f(x1)*f(x3) < 0 albo f(x3)*f(x2) < 0. Jeżeli f(x1)*f(x3) < 0, to za następny odcinek do podziału wybieramy [x1, x3] i dzielimy go na pół w punkcie x4 = (x1 + x3)/2. Jeżeli natomiast f(x3)*f(x2) < 0, to x4 = (x3 + x2)/2. Powtarzamy takie podziały aż do momentu, gdy różnica pomiędzy kolejnymi wartościami xk będzie mniejsza od założonej dokładności delta > 0, tj. |xk+1 – xk | < delta.
function x = polowienie(f, xp, xk, N) % polowienie( ) szuka pierwiastka funkcji f(x) w przedziale <xp, xk> % f funcja zdefiniowana przez użytkownika (jako inline) % xp punkt początkowy badanego przedziału zmiennej x % xk punkt końcowy badanego przedziału zmiennej x % N licznba podziałów odcinka w metodzie połowienia % x przybliżony pierwiastek równania f(x) = 0 w przedziale <xp, xk> a = f(xp); b = f(xk); if a*b > 0.0 disp(’f(x) ma taki sam znak na obydwu końcach badanego przedziału’) end for k = 1:N x = (xp + xk)/2; y = f(x) if a*y < 0 % porównanie znaków f(xp) i f(x) xk = x; % ustalenie punktu końcowego nowego przedziału else xp = x; % ustalenie punktu początkowego nowego przedziału end end
Testuj skuteczność metody bisekcji na funkcji f(x) = x2 – 9 oraz kilku innych funkcjach, pamiętaj, że należy je definiować jako inline.
Symulacje rzutu monetą
Funkcja rand wywołana bez podania argumentu zwraca liczbę losową (rozkład równomierny) z przedziału <0, 1>. Wykorzystamy to do symulowania wyniku rzutu idealną monetą, dla której prawdopodobieństwo wypadnięcia reszki albo orzełka jest identyczne, równe 0.5.
a = rand if (a < 0.5 ) wynik = 1; % umownie ”reszka” else wynik = 0; % umownie ”orzełek” end
Przypominam, że w MATLAB nie ma zmiennych typu logicznego, w zamian zmiennej lub wyrażeniu zawierającemu znaki logiczne przypisuje się wartość 1 , jeżeli jest prawdziwe, lub 0 gdy nieprawdziwe. Kod przedstawiony powyżej można wobec tego zapisać jeszcze prościej
a = rand wynik = (a < 0.5) % zmienna wynik przyjmie wartość 1 ( ”reszka” ) albo 0 ( ”orzełek” )
Wyniki n rzutów tą samą monetą (ekwiwalentnie, jednego rzutu n monetami) wygodnie potraktować jako składowe n – wymiarowego wektora. Symulacja wyników takiej serii rzutów nie przedstawia kłopotu
n = 50; a = rand(1, n); wynik = (a < 0.5); % teraz wynik będzie wektorem, każda jego składowa to 0 albo 1
Dla większych n ilość zer i jedynek w wektorze wynik powinna być prawie taka sama, co świadczyłoby o jakości zastosowanego algorytmu obliczania liczb losowych w funkcji rand, zob. zadanie 3.
Uogólnimy nieco rozważania, niech prawdopodobieństwo pojawienia się reszki jest p, orzełka q przy czym oczywiście p + q = 1 czyli q = 1 – p. Mógłby to być przypadek monety nieidealnej (uszkodzona, sfałszowana). Zapiszę ten kod w postaci funkcji,
function wynik = rzut(n, p) a = rand(1, n); wynik = (a < p);
Ilość zer w wektorze wynik będzie się teraz różnić od ilości jedynek, zwłaszcza gdy p znacznie różni się od 0.5, zobacz zadanie nr 6.
Rzut kostką
Rzut kostką kubiczną ze ścianami ponumerowanymi od 1 do 6. Symulacja rzutu idealną kostką sześcienną oznacza losowanie jednej liczby ze zbioru {1, 2, 3, 4, 5, 6}. Wektor o nazwie a zawierający n przypadkowych liczb zmiennoprzecinkowych z przedziału < 0, 6 > jest generowany przez kod,
a = 6 * rand(1, n);
Interesują nas liczby całkowite, wbudowana funkcja fix() obcina tą część liczby, która występuje po przecinku (kropce dziesiętnej), wobec tego instrukcja
a =fix( 6 * rand(1, n))+1;
tworzy wektor o n składowych, każda z nich jest losowana ze zbioru liczb naturalnych {1, 2, 3, 4, 5, 6}. Wektor a, obliczony na podstawie instrukcji (13.1), przechowuje zatem symulowane wyniki n rzutów idealną kostką sześcienną. Oto wynik, który ja dostałem symulując dziesięć ( n = 10) rzutów kostką
a =
- 4 6 3 5 1 2 1 4 1 5 % jakie liczby wylosował Twój komputer ?
Warto przedstawiony kod nieco uogólnić i zapisać w postaci funkcji, która będzie generowała tablicę/macierz składającą się liczb naturalnych losowanych z przedziału <1, m>,
function N = losowaN(m, r, s) % tworzy macierz o wymiarach s * t wypełnioną liczbami naturalnymi losowanymi z % przedziału <1, m> N = fix(m*rand(r, s)) + 1; end
W ramach testowania funkcji, wywołałem losowaN(6, 2, 2), oto wynik
ans =
- 3 2
- 2 4
Wywołanie funkcji losowaN() w pętli (np for lub while) wygeneruje serię macierzy wypełnionych liczbami całkowitymi, zob. zadanie 7.
Dynamika zysków/strat inwestycji
Podstawą wielu decyzji inwestycyjnych jest ocena dynamiki przyszłych zysków i strat (niestety). Wskutek złożoności procesów ekonomicznych, w przewidywaniach należy uwzględnić składową losową, reprezentującą brak pełnej informacji (niewiedzę). Zobaczmy jak zmieniają się, zarówno zysk jak i strata, wskutek działania czynnika losowego w bardzo prostym modelu. Wykonujemu n rzutów idealną monetą, kiedy wypada „reszka” zysk wzrasta o np. 1, gdy „orzeł” to tracimy 1.
function [zysk, strata] = dynamika(n) % liczy zysk i stratę w zależności od liczby n losowań (rzutów monetą) a = rand(1, n); % losowanie n liczb zmiennoprzecinkowych z przedziału <0, 1> a = (a<0.5); % jeżeli wartość składowej wektora a jest mniejsza od 0.5, % to zmień ją na 1, wartości pozostałych składowych zmień na 0 a = find(a>0); % find() zwraca indeksy/numery tych składowych a, które mają % wartość 1 zysk = length(a); % zwraca rozmiar wektora a, czyli ilość jego składowych. % Ponieważ aktualnie każda składowa ma wartość 1, to zwraca % ilość jedynek wylosowanych w n rzutach monetą, tego szukamy strata = n - zysk; end
Wywołując powyższą funkcję dla n = 10 000 otrzymałem
[zysk, strata] = dynamika(10000) zysk = 4975 strata = 5025 % !!
Niestety, bilans = zysk - strata = -50 jest ujemny, no po prostu katastrofa !! Cały dzień rzucania monetą (10 000 rzutów) i jestem „na minusie”. Może Ty będziesz miał/a więcej szczęścia w inwestycjach, funkcja dynamika( ) jest do dyspozycji, zob. także zadanie 9. Naturalnie, sukces lub porażkę inwestycyjną, w starciu z losowymi siłami rynkowymi, możemy prześledzić na wykresie,
plot(1: n, zysk) xlabel(‘ n, liczba rzutów ‘) ylabel(‘zysk(n)’)
Błądzenie przypadkowe
Rozpatruję najprostszy wariant błądzenia przypadkowego w jednym wymiarze. Cząstka/osobnik ...itp. startuje z punktu x = 0, w każdym kroku może się przemieszczać o jednostkową odległość w prawo (δx = +1) lub lewo (δx = -1). Prawdopodobieństwo przemieszczenia w prawo jest takie samo, jak w lewo. W jakim położeniu x znajdzie się cząstka po n krokach ? Oto kod,
function d = pijak(n) % oblicza położenie cząstki/osobnika po n losowych krokach, startuje od x = 0 x = rand(1, n); % losowanie n liczb z przedziału <0, 1> x = (x<0.5); % jeżeli wartość składowej wektora x jest mniejsza od 0.5, % to zmień ją na 1, wartości pozostałych składowych zmień na 0 a = find(x>0); % tworzy wektor a którego składowymi są indeksy/numery tych % składowych x, które mają wartość 1 np = length(a); % zwraca rozmiar wektora a. Każda składowa wektora a % reprezentuje (sygnalizuje) obecnośc jedynki w wektorze x. % Rozmiar wektora a jest wobec tego równy ilości jedynek w % wektorze x a zatem ilośc kroków w prawo nl = n - np; % ilośc kroków w lewo d = np – nl; % odleglość od położenia początkowego x=0 po n losowych krokach, % gdy d ma wartość ujemną, cząstka/osobnik jest na lewo od x=0 end
Wykonaj obliczenia dla różnych n.Nieco ogólniejszy przypadek występuje, gdy prawdopodobieństwo p kroku w prawo jest inne, niż w lewo (1 -p). Zmodyfikowany program,
function d = asym() % oblicza położenie cząstki po n losowych krokach, startuje od x = 0 % prawdopodobieństwo kroku w prawo jest p, w lewo (1-p) p = input('Podaj p = ') n = input('Podaj n = ') x = rand(1, n); x = (x<p); x = find(x>0); np = length(x); nl = n - np; d = np – nl; disp(['p = ' num2str(p)]) disp(['n = ' num2str(n)]) disp( ' ') disp(['x po n krokach = ' num2str(d)]) end
Wykonaj obliczenia dla różnych p przy ustalonym n, np. n = 100. Należy dodać, iż koncepcja błądzenia przypadkowego pojawia się także w ekonomii, np. przy modelowaniu fluktuacji cen.
Metoda Monte Carlo
Metoda Monte Carlo (nazwę zawdzięcza wykorzystaniu liczb losowych) jest stosowana w różnych dziedzinach nauki i techniki, np. w matematyce do numerycznego obliczania całek. Chcemy wyliczyć całkę \(\int\limits_a^b f(x) dx\). Jeżeli funkcja \(f(x)\) jest rzeczywista i ciągła to zachodzi relacja
- \[\int\limits_a^b f(x)\ dx = \bar{f} (b-a)\] ,
gdzie \(\bar{f}\) oznacza wartość średnią funkcji \( f(x) \) w przedziale < a, b >. Obliczanie całki oznaczonej sprowadza się zatem do wyliczenia wartości średniej \(\bar{f}\) funkcji \( f(x) \). Wystarczy wylosować (Monte Carlo !) z przedziału < a, b > n punktów (liczb), obliczyć sumę wartości funkcji \( f(x) \) w tych punktach i podzielić przez n. Oto kod programu, funkcję \( f(x) \) definiuję jako inline, aby ją łatwo zmieniać i testować program dla rozmaitych postaci funkcji \( f(x) \).
f = inline('x.*x', 'x'); function r = calka_mc(f, a, b, n) % metodą Monte Carlo liczy całkę oznaczoną funkcji f w granicach od a do b % n - liczba punktów (liczb) losowanych z przedziału < a, b > x = (b - a) * rand(1, n); srednia = sum(f(x))/n; r = srednia * (b - a); end
Jako funkcję podcałkową wpisałem \( f(x) = x^2 \), całka nieoznaczona z takiej funkcji wynosi \( x^3/3 \), całka oznaczona w przedziale < 0, 1 > wynosi 1/3,
r = calka_mc(f, 0, 1, 10000) r = 0.3407
Zwiększając liczbę n ("próbkowań" przedziału < 0, 1 > ) powiniśmy otrzymywać rezultat coraz bliższy dokładnej wartości 0.33333... .
Obliczanie liczby \( \pi \) metodą MC
Sztandarowy przykład zastosowania metody Monte Carlo to wyznaczenie wartości liczby \( \pi \). Rozważmy koło o promieniu r wpisane w kwadrat, bok tego kwadratu musi być wobec tego równy 2 r. Pole tego koła wynosi So = \(\pi r^2 \), natomiast pole kwadratu Sk = \( (2r)^2 = 4r^2 \). Stąd So/Sk = \(\pi/4\), a zatem \(\pi\) = 4 So/Sk. Wyobraźmy sobie, iż generujemy losowo bardzo dużą ilość punktów N leżących w kwadracie. Część z nich, oznaczę ją No, będzie również leżeć w kole wpisanym w ten kwadrat. Wówczas stosunek So/Sk możemy zastąpić przez No/N. Stosunek So/Sk nie zależy od promienia koła, w granicy dla dużych N to samo dotyczy No/N. Wygodnie wobec tego przyjąć, iż r = 1.
function [p] = PI(n) % oblicza liczbę pi metodą MC % n - liczba losowanych punktów czyli par liczb {x, y}, przy czym 0≤x<1, 0≤y<1 nkolo = 0; for k = 1 : n x = rand; y = rand; r = sqrt(x^2 + y^2); if (r <= 1) nkolo = nkolo +1; end end p = 4*nkolo/n; end
Oto przykład wywołania wyżej zdefiniowanej funkcji,
>> p = PI(1000) >> p = 3.1840 % blisko, ale nie idealnie
Dla celów dydaktycznych podaję drugą wersję tego programu,
function p = PI2() % oblicza Pi metodą MC nkolo = 0; rand('state',sum(100*clock)); % zmienia "seed" przy każdym uruchomieniu n = input('Podaj n = ') % prosi o podanie liczby losowanych punktów for k = 1: n x = rand; y = rand; r = sqrt(x^2 + y^2); if (r <= 1) % jeżeli wylosowany punkt jest w kole o promieniu r =1 nkolo = nkolo + 1; end end p = 4*nkolo/n; d = abs(pi - p); % błąd bezwzględny dw = 100*d/pi; % błąd względny podany w procentach disp(['Pi = ' num2str(p)]) disp(['delta = ' num2str(d)]) disp(['delta % = ' num2str(dw)]) end
Dopasowanie wielomianu (tzw. fitowanie)
Procedura fitowania modelowej funkcji do istniejącego zbioru danych występuje niezwykle często w naukach ścisłych, technicznych, ekonomii, ...itp. Poniżej przykład kodu źródłowego programu dopasowującego wielomian 5 stopnia do zbioru danych.
% program dopasowuje wielomian podanego stopnia do zbioru danych x = linspace(0, pi, 70); % generowanie punktów na osi x dane = cos(x) + 0.2*rand(1,length(x)); % utworzenie zbioru danych, w tym % przykładzie to "zaszumiona" funkcja cos(x) wsp = polyfit(x,dane,5); % dopasowanie wielomianu do zbioru danych. % Funkcja standardowa ployfit() zwraca % wektor współczynników wielomianu n-tego % stopnia (tutaj n = 5), dopasowanego do % wartości elementów wektora dane fit = polyval(wsp,x); % utworzenie wielomianu fit na podstawie % dopasowanych przez polyfit()współczynników % (metoda najmniejszych kwadratów)będących % skladowymi wektora wsp plot(x, dane, 'r*', x, fit, 'b+') % na wspólnym wykresie dane wyjściowe (dane) % oraz dopasowany wielomian (fit)
Zmodyfikowany program, w którym dane są "produkowane" przez wygodną do modyfikacji funkcję typu inline
% fituje wielomian stopnia n do funkcji f(x) w przedziale xp < x < xk % N - liczba punktów x w przedziale xp < x < xk % n - stopień dopasowywanego wielomianu clear; close all; % zamknięcie okien wcześniejszej grafiki, jeżeli są otwarte f = inline('(sin(x) + cos(x))', 'x'); % definicja przykładowej funkcji f(x) function wsp = fit_bis(f, xp, xk, N, n) x = linspace(xp, xk, N); wsp = polyfit(x, f(x), n); p = polyval(wsp, x); plot(x, f(x), 'r*', x, p, 'b+') end
Wywołanie,
f = inline('(sin(x) + cos(x))', 'x'); fit_bis(f, 0, 2*pi, 50, 3)
daje porównanie na jednym wykresie funkcji f(x) = sin(x) + cos(x) oraz dopasowanego do niej wielomianu 3 stopnia (metoda najmniejszych kwadratów). Testuj jakość dopasowywania (metoda najmniejszych kwadratów) dla różnych postaci funkcji f(x).
Krzywe Lissajous
Rodzinę krzywych Lissajous (fizyka, elektronika) opisują równania parametryczne
\(x(t)= A\sin(\omega_1 t + \delta), \quad y(t)= B\sin(\omega_2 t)\)
Kształt krzywych zależy od stosunku \(\omega_1\)/\(\omega_2\). Oto zbiór kilku instrukcji, dających wykres krzywej typu Lissajous
close all t = linspace(0, 2*pi, 3000); x = cos(4*t); y = sin(7*t); plot(x, y) end
- Animację tej krzywej generuje następujący kod,
close all t = linspace(0, 2*pi, 2000); x = cos(4*t); y = sin(7*t); comet(x, y) % zobacz, help comet end
Wykorzystałem tutaj funkcję comet() (zob. >> help comet), która podobnie jak plot() tworzy wykres, z tym, że pokazuje proces rysowania w czasie. Zobacz jeszcze animację rysowania okręgu, uruchamiając poniższy kod
close all t = 0:0.1*pi:2*pi; % zmieniaj krok czasowy 0.1 na większy oraz mniejszy, co zobaczysz ? figure; axis equal; axis([-1 1 -1 1]); hold on comet(cos(t), sin(t)) end
Zadania
- Program, przedstawiony w przykładzie Szereg geometryczny nie jest zabezpieczony przed beztroskim używaniem. Uczyń go bardziej dojrzałym. Niech pyta użytkownika o podanie q tj. ilorazu szeregu. Jeżeli poda on wartość q ≥ 1.0, w odpowiedzi ma pojawić się napis ” |q| musi być mniejszy od 1” i program ponownie czeka na podanie prawidłowej liczby. Po pięciokrotnym wpisaniu nieprawidłowej wartości q, ma pojawić się napis „bye, bye, bye..." i nastąpić zakończenie programu.
- Zauważ, iż rozwiązaniem równania f(x) = x2 – 2 = 0 jest xp = \(\sqrt{2}\). Oblicz przybliżoną wartość \(\sqrt{2}\) wykorzystując do tego celu metodę Newtona (zob. przykłady 8 i 9 ). Uwaga, gdyby program się „zapętlił”, to wiesz co masz zarobić, kombinacja klawiszy Ctr – c lub wprowadzić do kodu odpowiednie zabezpieczenie.
- Utwórz wykres funkcji f(x) = x3 − 10 w przedziale <0, 4>.
- Napisz skrypt (m-plik) programu, który na jednym wykresie przedstawi funkcje: x2, cos(x), ln(2*x), exp(-x) dla zmiennej x należącej do przedziału <0.5, 2>.
- Wygeneruj wektor o wymiarze n , ktorego składowymi są liczby losowe o rozkładzie równomiernym z przedziału < 0, 1>. Przedstaw składowe tego wektora, np. dla n = 100 na wykresie.
- Utwórz dwa wektory o n składowych każdy. Którąkolwiek ze składowych pierwszego wektora może być liczba 0 albo 1, reprezentująca wynik rzutu idealną monetą (1 ”reszka”, 0 ”orzełek”). Dokonaj w komputerze symulacji n rzutów idealną monetą, wpisując wyniki jako składowe pierwszego wektora. Przedstaw ilość zer i jedynek w tym wektorze na wykresie, dla różnych n, od 10 do 500. Składowymi drugiego wektora niech będa zera albo jedynki, losowane dla monety sfałszowanej, pokazującej reszkę średnio 80 razy na 100 rzutów.
Na wykresie porównaj ilości zer i jedynek w pierwszym i drugim wektorze dla n od 10 do 500. Zastosuj różne kolory dla zwiększenia czytelności. - Wygeneruj 50 macierzy o wymiarach 4 * 4 każda, wypełnionych liczbami naturalnymi losowanymi z przedziału <0, 10> . Oblicz macierz będącą średnią z sumy tych macierzy, następnie odrzuć mantysy wszystkich elementów tej macierzy „średniej”, aby dostać macierzy „średnią całkowitą” tj wypełnioną elementami całkowitymi.
- Utwórz wektor a zawierający 20 liczb całkowitych. Utwórz wektor b, którego składowymi są liczby parzyste występujące w a. Następnie utwórz wektor c, zawierający te składowe wektora a, które są liczbami nieparzystymi.
- Utwórz macierz A o sześciu wierszach i ośmiu kolumnach, której elementami będą liczby całkowite wybrane losowo z przedziału [−12, 12].
- Napisz program tworzący wektor o sześciu składowych, który gromadzi liczby naturalne losowane z przedziału < 1, 49 >, przy czy losowane liczby nie mogą się powtarzać. Tak, to symulacja tzw. ciągnienia w loterii typu LOTTO.
- Utwórz funkcję, która liczy zysk i stratę po n rzutach specjalnie spreparowaną monetą, dla której prawdopodobieństwo pokazania „reszki” jest p, orzełka naturalnie 1 - p ( 0 < p < 1). Wykonaj symulacje rzutów tej monety dla wielu wartości n , np. n = 100: 100: 1000. Porównaj na jednym wykresie przebiegi bilansów (bilans = zysk – strata) w funkcji n. Obserwuj na wykresie zależność zysku od parametru p („stopnia nierzetelności”) monety przy ustalonym n, np n = 1000. Przeczytaj przykład Dynamika zysków/strat inwestycji.
- Narysuj wykres wielomianu f(x) = x4 − 4x2 - 3, używając do tego celu 100 punktów z przedziału −3 < x < 3.
- Pokaż na wspólnym wykresie przebiegi następujących funkcji f(t): exp(−2t)cos(ωt), exp(−t)[4cos(ωt) + sin(2ωt)] w przedziale 0 ≤ t ≤ 15, dla trzech wartości ω = 2, 5, 10.
- Testuj skuteczność obliczania całek oznaczonych metodą Monte Carlo. Wykonaj całkowania numeryczne dla następujących funkcji f(x): cos(x), sin(x), x3sin(x)cos(x), exp(-x)sin2(x)cos3(4x2) dla x z przedziału <0, 2\(\pi\) >. Pokaż na wykresie, jak zmienia się numerycznie obliczana wartość całki \(\int\limits_0^\pi sin((x) dx\) przy rosnącej liczbie n próbkowań przedziału < 0, \(\pi\) >, np dla n = 10, 100, 1000, 10 000, 100 000, 1 000 000.
- Wykonaj dopasowanie wielomianu n-tego stopnia do funkcji f(x) = 2x2cos(x)sin(2x), n = 2, 3, 4. Zobacz na wykresie, jak zmiana n wpływa na precyzję dopasowania.
- Dopasuj wielomiany do następujących funkcji f(x): 3x3cos2(x)sin2(4x), x4exp(x2)sin(x), xexp(x3)cos(5x). Przedstaw na jednym wykresie wyniki dopasowania wielomianiu 5-tego stopnia (n = 5)do tych funkcji f(x).
- Rodzinę krzywych Lissajous opisuje układ równań
\(x(t)= A\sin(\omega_1 t + \delta)\) oraz \( y(t)= B\sin(\omega_2 t).\)
Kształt krzywych zależy od stosunku \(\omega_1\)/\(\omega_2\).Wykreśl na jednym wykresie dwie krzywe dla następujących wartości parametrów: a) A = B = 1 oraz \(\delta = 0 \), b) A = B = 1 oraz \(\delta \) = \(\pi\)/2. Podstaw za \(\omega_1 \) i \(\omega_2\) kilka różnych liczby naturalnych i pokaż na wspólnym wykresie odpowiednie krzywe Lissajous dla 0 < t ≤ 2 pi. Jak wzrost całkowitych wartości \(\omega\) wpływa na kształt tych krzywych. Zobacz co się stanie, gdy zwiększysz zakres t np. 20 razy. Czy krzywe się zamykają ? Obserwuj krzywe Lissajous w przedziale 0 < t ≤ 2\(\pi\), gdy \(\omega\) jest liczbą niewymierną, np. \(\omega\) = e, \(\pi\), \(\sqrt{3} \) … Zwiększ znacznie zakres t, czy krzywe się zamykają ? Po wykonaniu tych wykresów, być może potrafisz uzupełnić stwierdzenie; krzywe Lissajous są zamknięte, gdy stosunek \(\omega_1\)/\(\omega_2\) jest ..... , no właśnie ??? - Jeżeli dość malownicza rodzina krzywych Lissajous wzbudziła Twoje zainteresowanie, uruchom kod, który tworzy wykres krzywej typu Lissajous
close all t = linspace(0, 2*pi, 3000); x = cos(4*t); y = sin(7*t); plot(x, y) end
- Animację tej krzywej generuje następujący kod,
close all t = linspace(0, 2*pi, 2000); x = cos(6*t); y = sin(3*t); comet(x, y) end
- Poczytaj (help comet) o funkcji comet(), eksperymentuj z animacją krzywych Lissajous lub innych.