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, x3, x5, x7
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 x0 ,
f(x) ≈ f(x0 ) + f ` (x0 ) (x − x0 )(1)
gdzie f `(x0 ) jest pochodną funkcji dla x = x0 . 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 ≈ x0 − f(x0 ) / f `(x0 ) (2)
Dostaliśmy przepis na znajdywanie przybliżonej wartości x, dla której f(x) = 0. Zaczynamy procedurę od wytypowania wartości startowej x0 , która według dostępnych przesłanek wydaje się być bliska pierwiastkowi xp badanego równania. Następnie generujemy sekwencję wartości {x0 , x1 , x2 , x3 , . . .} według reguły:
xi+1 = xi − f(xi ) / f `(xi ) (3)
Jeżeli funkcja f(x) ma w miarę regularny przebieg oraz x0 jest blisko prawdziwej wartości pierwiastka badanego równania, to sekwencja {x0 , x1 , x2 , x3 , . . .} szybko zmierza do xp . Napiszmy program, który oblicza kolejne wartości xi+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) = x4 – 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) = x4 – 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(xm) – f(xp) | < epsilon. Ponieważ xp jest pierwiastkiem badanego równania (czyli f(xp) = 0), komputer ma przerwać iteracje, gdy |f(xm)| < 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 x0 leży zbyt daleko od xp.
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 xk będzie mniejsza od założonej dokładności delta > 0, tj. |xk+1 – xk | < 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) = x2 – 9 oraz kilku innych funkcjach, pamietaj, że należy je definiować jako inline.