Z Skrypty dla studentów Ekonofizyki UPGOW

Wersja ZygmuntGburski (dyskusja | edycje) z dnia 17:06, 11 kwi 2010

(różn.) ← poprzednia wersja | przejdź do aktualnej wersji (różn.) | następna wersja → (różn.)

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

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 ?

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

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

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) w przedziale [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.
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.
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.
Rodzinę krzywych Lissajous 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 $\frac{\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 $\frac{\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(4*t);
y = sin(7*t);
comet(x, y)
end

Poczytaj (help comet) o funkcji comet(), eksperymentuj z animacją krzywych Lissajous lub innych.