Z Skrypty dla studentów Ekonofizyki UPGOW
Linia 54: | Linia 54: | ||
f = [x; x.^3; x.^5: x.^7]; | f = [x; x.^3; x.^5: x.^7]; | ||
plot(x, f)</source> | 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 < 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 i=1:n | ||
+ | sum =sum+ 1/i + 1/((i+2)*(i+3)) | ||
+ | end | ||
+ | disp( ’ Suma wynosi ’ sum) | ||
+ | </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<sub>0</sub> , | ||
+ | |||
+ | f(x) ≈ f(x<sub>0</sub> ) + f ` (x<sub>0</sub> ) (x − x<sub>0</sub> )(1) | ||
+ | |||
+ | gdzie f `(x<sub>0</sub> ) jest pochodną funkcji dla x = x<sub>0</sub> . 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<sub>0</sub> − f(x<sub>0</sub> ) / 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<sub>i</sub> ) (3) | ||
+ | |||
+ | Jeżeli funkcja f(x) ma w miarę regularny przebieg oraz x<sub>0</sub> jest blisko prawdziwej wartości pierwiastka badanego równania, to sekwencja {x<sub>0</sub> , x<sub>1</sub> , x<sub>2</sub> , x<sub>3</sub> , . . .} szybko zmierza do x<sub>p</sub> . Napiszmy program, który oblicza kolejne wartości x<sub>i+1</sub> 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<sup>4</sup> – 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<sup>4</sup> – 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<sub>m</sub>) – f(x<sub>p</sub>) | < ''epsilon''. Ponieważ x<sub>p</sub> jest pierwiastkiem badanego równania (czyli f(x<sub>p</sub>) = 0), komputer ma przerwać iteracje, gdy |f(x<sub>m</sub>)| < 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<sub>0</sub> leży zbyt daleko od x<sub>p</sub>. | ||
+ | |||
+ | ==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<sub>k</sub> będzie mniejsza od założonej dokładności ''delta'' > 0, tj. |x<sub>k+1</sub> – x<sub>k</sub> | < ''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<sup>2</sup> – 9 oraz kilku innych funkcjach, pamietaj, że należy je definiować jako '''inline'''. |
Wersja z 13:20, 25 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, 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.