Z Skrypty dla studentów Ekonofizyki UPGOW

Wersja z 14:14, 26 mar 2010

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 przedz. <0,1>
srednia(a)                           
ans  =                     %  wynik zależy od wartości wylosowanych elementów wektora a

Obliczanie sumy

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

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

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 (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 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 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, 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 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 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],

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

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 a 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ą (10000 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)’)

Zadania

1. . Program, przedstawiony w przykładzie Szereg geometryczny nie jest zabezpieczony przed bezmyślnym używanie. 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.
2. 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.
3. Utwórz wykres funkcji f(x) = x³ − 10 w przedziale [0, 4].
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].
5. 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.
6. Utwórz dwa wektory o n składowych każdy. Jakąkolwiek składową 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.
7. 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.
8. Napisz program tworzący wektor 1*6 , 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.
9. 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: 500. 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 13