Z Skrypty dla studentów Ekonofizyki UPGOW

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, x³, x⁵, x⁷

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 x₀ ,

f(x) ≈ f(x₀ ) + f ` (x₀ ) (x − x₀ ) (1)

gdzie f `(x₀ ) jest pochodną funkcji dla x = x₀ . 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 ≈ x₀ − f(x₀ ) / f `(x₀ ) (2)

Dostaliśmy przepis na znajdywanie przybliżonej wartości x, dla której f(x) = 0. Zaczynamy procedurę od wytypowania wartości startowej x₀ , która według dostępnych przesłanek wydaje się być bliska pierwiastkowi x_p badanego równania. Następnie generujemy sekwencję wartości {x₀ , x₁ , x₂ , x₃ , . . .} według reguły:

x_i+1 = x_i − f(x_i ) / f `(x_i ) (3)

Jeżeli funkcja f(x) ma w miarę regularny przebieg oraz x₀ jest blisko prawdziwej wartości pierwiastka badanego równania, to sekwencja {x₀ , x₁ , x₂ , x₃ , . . .} szybko zmierza do x_p . Napiszmy program, który oblicza kolejne wartości x_i+1 na podstawie równ. 3, aż do N-tej. Liczbę N podaje użytkownik, program po wykonaniu obliczeń wypisuje x_N , 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) = x⁴ – 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 x₀ 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 x_N, 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 x₀. W naszym przykładzie, rozwiązanie równania f(x) = x⁴ – 16 = 0 jest oczywiste, x_p = 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(x_m) – f(x_p) | < epsilon. Ponieważ x_p jest pierwiastkiem badanego równania (czyli f(x_p) = 0), komputer ma przerwać iteracje, gdy |f(x_m)| < 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 x₀ leży zbyt daleko od x_p.

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ł x₁ ≤ x ≤ x₂. Skoro w tym przedziale znajduje się pierwiastek, to f(x₁)*f(x₂) < 0. Obliczamy x₃ = (x₁ + x₂)/2 dzieląc w ten sposób przedział [x₁, x₂] na dwie równe części. Wówczas, albo f(x₁)*f(x₃) < 0 albo f(x₃)*f(x₂) < 0. Jeżeli f(x₁)*f(x₃) < 0, to za następny odcinek do podziału wybieramy [x₁, x₃] i dzielimy go na pół w punkcie x₄ = (x₁ + x₃)/2. Jeżeli natomiast f(x₃)*f(x₂) < 0, to x₄ = (x₃ + x₂)/2. Powtarzamy takie podziały aż do momentu, gdy różnica pomiędzy kolejnymi wartościami x_k będzie mniejsza od założonej dokładności delta > 0, tj. |x_k+1 – x_k | < 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) = x² – 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 S_o = $\pi r^2 $, natomiast pole kwadratu S_k = $ (2r)^2 = 4r^2 $. Stąd S_o/S_k = $\pi/4$, a zatem $\pi$ = 4 S_o/S_k. Wyobraźmy sobie, iż generujemy losowo bardzo dużą ilość punktów N leżących w kwadracie. Część z nich, oznaczę ją N_o, będzie również leżeć w kole wpisanym w ten kwadrat. Wówczas stosunek S_o/S_k możemy zastąpić przez N_o/N. Stosunek S_o/S_k nie zależy od promienia koła, w granicy dla dużych N to samo dotyczy N_o/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) = x² – 2 = 0 jest x_p = $\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) = x³ − 10 w przedziale <0, 4>.
Napisz skrypt (m-plik) programu, który na jednym wykresie przedstawi funkcje: x², 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) = x⁴ − 4x² - 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³sin(x)cos(x), exp(-x)sin²(x)cos³(4x²) 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) = 2x²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³cos²(x)sin²(4x), x⁴exp(x²)sin(x), xexp(x³)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