Z Skrypty dla studentów Ekonofizyki UPGOW

Nie masz uprawnień do edytowania tej strony, z następującego powodu:

Tekst źródłowy strony można podejrzeć i skopiować.

=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 >> __TOC__ ==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. <source lang="matlab">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 </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. ==Rysowanie okręgu== Rysowanie okręgu o środku w początku kartezjańskiego układu współrzędnych i jednostkowym promieniu r = 1. <source lang="matlab">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 </source> ==Wykorzystanie liczb zespolonych do rysowania okręgu== Zobacz także poprzedni przykład <source lang="matlab">alfa = 0 : pi/60 : 2*pi; z = exp(i*alfa); plot(z), axis(’square’)</source> ==Rysowanie kilku funkcji na jednym wykresie== Następujący kod powoduje przedstawienie na jednym wykresie grupy funkcji f(x) = x, x³, x⁵, x⁷ <source lang="matlab">x = 0: .2: 2; f1 = x; f2 = x.^3; f3 = x.^5; f4 = x.^7; plot(x, f1, x, f2, x, f3, x, f4)</source> Uwaga: Jeżeli ostatnią linię zastąpimy przez <source lang="matlab">plot(x, f1, ‘black’, x, f2, ‘blue’, x, f3, ‘green’, x, f4, ‘red’)</source> 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, <source lang="matlab">x = 0: .2: 2; f = [x; x.^3; x.^5: x.^7]; plot(x, f)</source> ==Wykres parametryczny== Wykres parametryczny krzywej opisanej dwoma równaniami: x (t) = t cos(t), y(t) = t sin(t) dla 0 ≤ t ≤ 2 π . <source lang="matlab">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’)</source> ==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> ==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> Zastosowanie, <source lang="matlab"> 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</source> ==Obliczanie sumy== Obliczanie sumy <math> \sum_{i=1}^N \frac{1}{i}+\frac{1}{(i+1)(i+2)} </math> dla wartości N podanej przez użytkownika <source lang="matlab">N=input(’Podaj N ') for k = 1: n suma =suma + 1/k + 1/((k + 2)*(k + 3)) end disp( ’ Suma wynosi ’ suma) </source> ==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<sub>0</sub> ) (x − x<sub>0</sub> )(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<sub>0</sub> ) (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<sub>0</sub> , która według dostępnych przesłanek wydaje się być bliska pierwiastkowi x<sub>p</sub> badanego równania. Następnie generujemy sekwencję wartości {x<sub>0</sub> , x<sub>1</sub> , x<sub>2</sub> , x<sub>3</sub> , . . .} według reguły: x<sub>i+1</sub> = x<sub>i</sub> − f(x<sub>i</sub> ) / 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) . <source lang="matlab">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 </source> Zapisujemy ten kod w pliku nt.m w bieżącym katalogu roboczym. W linii komend definiujemy sposobem '''inline''' funkcję, np. f(x) = x⁴ – 16 <source lang="matlab">f = inline(‘ x^4 - 16 ‘)</source> oraz jej pochodną, <source lang="matlab">pf = inline(‘ 4 * x^3 ‘)</source> Ostatni krok to wywolanie fukcji nt(f, pf, x0, N) dla konkretnych wartości parametrów x0 oraz N, na przykład <source lang="matlab">xp = nt(f, pf, 1.0, 20); disp(xp)</source> 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). <br> <br> 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 <source lang="matlab">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</source> Dalej należy wywołać funkcję ntmod( ) dla wybranych parametrów, np. <source lang="matlab">xp = nt(f, pf, 1.0, .001); % wykonaj obliczenia również dla epsilon = 0.1, 0.01, 0.00001 disp(xp)</source> 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 &le; x &le; 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''. <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 % 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 </source> 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. <source lang="matlab">a = rand if (a < 0.5 ) wynik = 1; % umownie ”reszka” else wynik = 0; % umownie ”orzełek” end </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 <source lang="matlab">a = rand wynik = (a < 0.5) % zmienna wynik przyjmie wartość 1 ( ”reszka” ) albo 0 ( ”orzełek” ) </source> 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 <source lang="matlab">n = 50; a = rand(1, n); wynik = (a < 0.5); % teraz wynik będzie wektorem, każda jego składowa to 0 albo 1 </source> 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, <source lang="matlab">function wynik = rzut(n, p) a = rand(1, n); wynik = (a < p); </source> 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, <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 <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ą a = :4 &nbsp; 6 &nbsp; 3 &nbsp; 5 &nbsp; 1 &nbsp; 2 &nbsp; 1 &nbsp; 4 &nbsp; 1 &nbsp; 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], <source lang="matlab">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 </source> W ramach testowania funkcji, wywołałem losowaN(6, 2, 2), oto wynik ans = :&nbsp;3 &nbsp; 2 :&nbsp;2 &nbsp; 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. <source lang="matlab">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</source> Wywołując powyższą funkcję dla n = 10 000 otrzymałem <source lang="matlab">[zysk, strata] = dynamika(10000) zysk = 4975 strata = 5025 % !! </source> 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, <source lang="matlab">plot(1: n, zysk) xlabel(‘ n, liczba rzutów ‘) ylabel(‘zysk(n)’) </source> =Zadania= #. Program, przedstawiony w przykładzie [[Programowanie_Przykłady#Szereg_geometryczny|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 &ge; 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 = <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³ − 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. 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.<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. #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. #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 [[Programowanie_Przykłady#Dynamika_zysków/strat|13]]

Wróć do strony Programowanie Przykłady.