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

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.

q  =   .7;                          %    inicjalizacja ilorazu szeregu
sum_szer  =  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 wartość zero
 
while sum_szer  >   pom             %    początek pętli while
pom  =  sum_szer;
sum_szer =  sum_szer + q^m;
 m  =  m + 1;
end                                 %     koniec pętli while
 
disp(suma_szer)                     %     pokazuje na ekranie obliczona 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.

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’)

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)

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

Zastosowanie,

a = 1 : 49;                %     tworzenie wektora o składowych od 1 do 49
srednia(a)
ans  =
           25
a = rand(1, 30);           %     tworzy wektor 30 liczb przypadkowych  z   <  0, 1  >
srednia(a)                           
ans  =                     %  wynik zależy od wartości wylosowanych elementów wektora a

Obliczanie sumy

Obliczanie sumy \( \sum_{i=1}^N \frac{1}{i}+\frac{1}{(i+1)(i+2)} \) dla wartości N podanej przez użytkownika

N=input(’Podaj N ')
for i=1:n
sum =sum+ 1/i + 1/((i+2)*(i+3))
end
disp( ’ Suma wynosi ’ sum)

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 xN (przybliżony pierwiastek równania f(x) = 0) .

function x  =  nt(f, pf, x0, N)          %   początek definicji funkcji, którą nazwałem 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 sposobem 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 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) = x⁴ – 16 = 0 jest oczywiste, xp = 2 i -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);      % wykonaj obliczenia również dla epsilon = 0.1, 0.01, 0.00001
disp(xp)

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 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 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 owego przedziału
end
end

Testuj skuteczność metody bisekcji na funkcji f(x) = x² – 9 oraz kilku innych funkcjach, pamietaj, ż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 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 kostka

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 między 0.0 a 6.0 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],

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

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.