Semestr zimowy 2009/10

• Wykład
• Ćwiczenia/Lab
  • Program bieżącego labu

Metody numeryczne (dla informatyków)

Zadanie komputerowe:

Literatura:

• [KC2006] D.Kincaid, W.Cheney, Analiza numeryczna, WNT, 2006
• [Mos2002] Krzysztof Moszyński, Metody numeryczne dla informatyków, skrypt, plik ps, 2002
• [Pla2002] Leszek Plaskota, Dwanaście wykładów z matematyki obliczeniowej, skrypt, plik pdf, 2002
• [DJ1982] Maksymilian Dryja, Janina i Michał Jankowscy, Metody numeryczne, WNT, 1982.
• [FMW2005] Z. Fortuna, B. Macukow, J. Wasowski, Metody numeryczne , WNT, 2005. Wydanie 7.

Inne użyteczne linki

• Lloyd N.Trefethen, Numerical analysis (The Princeton Companion to Mathematics, Princeton University Press, 2007) plik pdf
• Lloyd N.Trefethen, The definition of numerical analysis, (SIAM News, Nov 1992) zgzipowany poscript.

Ćwiczenia/Lab MN

• Gr 1 piątek 1015-1145
• Gr 2 piątek 1215-1345

Zaliczenie ćwiczeń

Program ćwiczeń

1. 9-10-2009. Skalarne metodu rozwiązywania równań nieliniowych: metoda Herona (czyli m Newtona dla równania x*x-a=0), odwracanie bez dzielenia (m. Newtona dla 1/x-a=0 dla jakich x_0 m. Newtona zbiega?), metoda iteracji prostych- pokaż że x_n=cos(x_{n-1}) zbieżne dla dowolnego x_0. Zera wielokrotne: pokaż że dla zera 2-krotnego met. Newtona jest zbieżna lokalnie.
2. 23-10-2009 Równania nieliniowe cd. Jak badając pochodne funkcji iteracji \Phi dla m iteracyjnej x_n=\Phi(x_{n-1}) w p-cie stałym iteracji x=\Phi(x) zbadać rząd zbieżności metody? Zastosować do m. Newtona przy standardowych założeniach. Jak zmodyfikować m Newtona aby zbiegała kwadratowo dla zera 2-krotnego. Metoda wielowymiarowa Newtona (na przykładzie prostej funkcji F(x,y)=(x+y;2x*x+y*y-4)^T- gdzie kolejna iteracja dobrze określona, czy zbieżna kwaratowow lokalnie w otoczeniu rozwiązania.
3. 06-11-2009 Równania liniowe. Rozkład LU czyli eliminacja Gaussa. Zad 1: LU dla macierzy silnie diagonalnie dominującej można wykonać bez permutacji (czyli wyboru elementu głównego). Zad 2: Rozkład LU dla macierzy trójdiagonalnej bez permutacji (do domu z permutacjami wierszy czyli częściowym wyborem elementu głównego). Do domu: LU dla m. trójdagonalnej z niezerowymi rogami (czyli niezerowe elementy mogą sie pojawić na diagonali, pod- i naddiagonali i w rogach). Zad 3: Mając szybki solver dla konkretnej macierzy NxN A, jak rozwiązać szybko z wykorzystaniem tego solvera układ Ax+by=f, cx+ay=g gdzie c - wektor poziomy, b - wektor pionowy a,g- skalary - [x;y] rozwiązanie - tu oczywiście x- wektor a y - skalar. Podać warunek konieczny i dostateczny na to aby ten układ miał rozwiązanie dla dowolnych f i g. Zad 4: (tylko grupa 2) Dla macierzy Q ortogonalnej (Q*Q^T=Q^TQ=I) i A dowolnej kwadratowej pokaż że QA i AQ mają takie same normy drugie jak A (to samo dla normy Frobeniusa \|A\|_F=\sqrt{\sum_{i,j}|A_{ij}|^2}).
4. 20-11-2009 Kartkówka z prac domowych. LU c.d. i LZNK. Zad 1: Metoda Choleskiego tzn rozkład LU bez wyboru elementu czyli tak naprawdę A=LDL^T to L dolnotrójkątna z 1 na diagonali a D - m diagonalna dla A=A^T>0 (do domu A=L^T L - L dolnotr. - policzyć wprost wzory ne elementy L z tego rozkładu i że dla A 3-diagonalnej macierz dolnotr. L jest 2-diagonalna). Zad 2: Zadanie znajdowania krzywej opisanej równaniem (a,b, parametry) ax^2+by^2=1 najlepiej pasującej do zadanych M punktów (xk,yk) jako LZNK (na labieteż to w octavie zrobimy).
5. 04-12-2009 Macierze rzadkie. Metody iteracyjne rozwiązywania układów równań liniowych. Zad 1 Zbieżność met. G-S dla m. silnie diag. dominujących w dyskretnej normie supremum. Zad 2 Dla jakich wartości parametru \tau met. Richardsona zbiega dla A=A^T>0 takiej że widmo A w przedziale [c,d]. Dla jakiego \tau szybkość zbieżnośic jest największaw normie drugiej i energetycznej \|x\|_A? Mamy met. iteracyjną x_{n+1}=x_n+2*a*A*Ax_n + a*g przy czym A=A^T o wartościach własnych {-1,2,5}. Dla jakich wartości a metoda zbiega? Określić w terminach A,a,g do czego zbiegnie x_n (o ile zbiegnie). Podać a takie że zbieżność w normie dwa metody najszybsza. Czy jeśli x_0=c*x_1+d*x_2 dla pewnych stałych c,d i wektorów własnych x_1,x_2 dla wart. własnych -1,2 to można dla pewnego a przyspieszyć zbieżność metody? Zad 3 Napisać w pseudokodzie mnożenie przez macierz w formacie spakowanych wierszy (do domu kolumn).
6. 11-12-2009Zadanie własne i interpolacja
   1. zaproponuj algorytm który przy pomocy p. Householdera sprowadza A=A^T do macierzy podobnej 3diagonalnej. Oszacuj koszt.
   2. pokaż że wielomian charakterystyczny dla macierzy [0 1;-a -b] to x^2+bx+a. Czy dla danego wielomianu (-1)^n+a_{n-1} x^{n-1}+...+a_0 istnieje macierz A taka że jej wiel. charakterystyczny to ten wielomian? (Macierz ma ostatni wiersz: (-a_0,..,-a_{n-1}) a tak 1 na pierwszej nadprzekątnej (czy 1 na pozycjach (i,i+1) i=1,..,n) i poza tym elementy zerowe).
   3. dla funkcji która w węzłach -1,0,1,2 ma wartości 0, 1, 2, 9 znajdź wiel interpolacyjny lagrange'a na 3 sposoby - w bazie Lagrange'a, w bazie Newtona korzystając z wzoru Newtona i przy pomocy tablic różnic dzielonych - porównaj wyniki. Oszacuj błąd interpolacji w normie max przy założeniu że \|f^(4)\|_{max}<1. (Oszacowanie do domu)
   4. udowodnij wzór na różnice dzielone F[x_0..x_N]= \sum_{k=1}^N F(x_k)/((x_k-x_0)...(x_k-x_{k-1})(x_k-x_{k+1})...(x_k-x_N) Wsk: Skorzystaj z wzoru Newtona wstawiająć postać wiel. interpol. L w bazie Lagrange'a z prawej strony wzoru Newtona. (do domu dla chętnych)
   5. niech L_N f wiel. interp dla f(x)=sin(x) na siatce równoodległej na [0,2pi]. Oszacuj szybkość zbieżności L_N f i S_N f do f w normie supremum. (tzn podaj błędy) (do domu to samo ale dla sin(2x))
7. 15-01-2009 Kartkówka z prac domowych. Splajny i kwadratury.
   1. niech S_N f splajn kubiczny z warunkami hermitowskimi dla f(x)=sin(x) na siatce równoodległej na [0,2pi]. Oszacuj szybkość zbieżności S_N f do f i pochodnej S_Nf do f' w normie supremum. (tzn podaj błędy) (do domu to samo ale dla dla sin(2x) i x^{3/2} - ta druga funkcja nie jest C^2!)
   2. niech dla równomiernego podziału \PI: x_k=k*h h =2\pi/N odcinka [0,2\pi] S_N f to splajn z S^0_2 interpolujący f(x)=sin(x) w węzłąch i środkach podocinków. Czy jest wyznaczony jednoznacznie? Oszacuj szybkość zbieżności S_N f do f w normie supremum. (tzn podaj błędy) (do domu to samo ale dla sin(2x))
   3. rząd i zbieżność złożónej kwadratury prostokątów dla funkcji klasy C^1 i C^2.

Program labów

• Lab 1 (2-10-09) - wstępne zapoznanie się z octavem - octave jako kalkulator naukowy, operator dwukropek : tworzenie macierzy, wektorów, zapisywanie/czytanie do/z plików w obu formatach - tekstowym i binarnym, tworzenie macierzy z podmacierzy, wycinanie podmacierzy itp, podstawowe operacje na macierzach - mnożenie, dodawanie, transponowanie, funkcje od macierzy, normy wektorów/macierzy. Skrypty i funkcje (m-pliki) w octavie. Instrukcje warunkowe (if else endif; switch case endswitch), pętle (while( ) do .. endwhile; do .. until( );for .. endfor). Wskaźniki do funkcji (handle - operator @). Bisekcja. Zad 1 Utwórz dowolne macierze 3x4 A i 3x5 B - a następnie macierz 3x8 C której pierwsze 3 kolumny to A a kolejne to B. Teraz z tej macierzy 'wytnij' podmacierz D składającą się z 1 głównego minora tzn 3x3 od C(1,1) do C(3,3). Zamień kolejność kolumn D. Wstaw D z powrotem do C jako główny minor. Policz sin od D. Zapisz D do pliku (binarnego i ASCII) - zamień element D(1,1) na -100 i wczytaj nową macierz do octave'a jako DD. Policz normy (różne) macierzy DD. Zad 2 Policz dyskretną normę max od (sin(x))^2 na [0,1] (wektorowo- czyli bez użycia pętli). Zad 3 Zaimplementuj metodę bisekcji w skrypcie - dla rozwiązania równania cos(x)=0 na odcinku [0,2] - sprawdź jak działa. Zad 4 Napisz funkcję octave'a w pliku bisekcja.m - której parametrami będą nazwa funkcji (czy funkcja - jak to zrobić odsyłam do manuala czy stron www), końce przedziału a i b, epsilon - żądana dokłądnośc bezwzględna (błąd ma być mniejszy od epsilona), max. ilośc iteracji. Funkcja ma zwrócić przybliżone rozwiązanie, ilość iteracji, oszacowanie błędu, i kod wyniku - 0 metoda zbiegła, 1 - metoda zatrzymała się z powodu przekroczenia max ilości iteracji, 2 - wartości w końcach odcinka funkcji mają ten sam znak. Funkcja ma działąć również jak podamy tylko trzy czy cztery pierwsze argumenty - i wtedy max. ilośc iteracji domyslnie ma być 100 , a epsilon 1e-5 (jak tylko 3 argumenty). Zad 5 (to i tak będzie) - znajdź w manualu funkcję octave'a - rozwiązująca równania nieliniowe - i przetestuj ją na kilku prostych przykładach skalarnych np cos(x)=0, x^2-2=0, x^10-2=0, x-sin(x)=0 itp
• Lab 2 (16-10-2009)- metody rozwiązywania równań nieliniowych cd.
  • Zad 1 Zaimplementować metodę Newtona w octavie i przetestować jej zbieżnośc dla następujących funkcji: x*x-4, x*x*x-27, exp(x)-2, sin(x) dla x=\pi, (x-2)^k dla k=2,4,8,16, 1/x-a dla danych a=0.5,2,4,100 (oczywiście implementując bez dzielenia). Dla wszystkich tych funkcji znamy rozwiązania więc można wyświetlać na ekranie bład e_n=x_n-r (r - rozwiązanie) i |e_n|/|e_{n-1}|^p dla p=1,2,3. Dobierać różne wartości startowe x_0.
  • Zad 2 Jak w zad 1 ale dla metody siecznych w szczególności przetestować czy zachodzi e_n/(e_{n-1}e_{n-2}) asymptotycznie zbiega do stałej i czy | e_n|/|e_{n-1}|^p dla p=(1+ \sqrt(5))/2 zbiega do stałej np dla x*x-2.
  • Zad 3 Sprawdzić czy metoda iteracji prostych x_n=cos(x_{n-1}) zbiega do x=cos(x). Zbadac czy zbieżność liniowa.
  • Zad 4 Porównać wyniki z zad 1 z funkcją octave'a fsolve() (pomoc: help fsolve).
  • Zad 5 Zaimplementować przybliżoną metode Newtona w której pochodną przybliżamy ilorazem różnicowym tzn x_{n+1}=x_n - f(x_n)*h/(f(x_n+h)-f(x_n)) dla ustalonego h. Przetestować różne h np 1e-4,1e-7,1e-10 itp Porównać zbieżność z dokłądną met Newtona (szczególnie ostatnie iteracje) dla funkcji z zadania 1.
  • Zad 6 Rozwikływanie funkcji: dla funkcji y(x) zadanej równaniem g(x,y)=2x^2+3y^2-3=0 znaleźć wartości yk=y(x(k)) dla xk=k*h dla k=-N,...,N i h=1/N (N - 10,20,40,...) Znajdować yk rozwiązując układ równań g(xk,yk)=0. Jak w kolejnych krokach dobierać przybliżenie startowe w metodzie Newtona?
  • Zad 7 Odwracanie funkcji: mamy daną funkcje np sin(x)+2*x znależć wartości funkcji odwrotnej na odcinku [0,5] na siatce k*h dla k=0,..,100. Narysować wykres funkcji. (Sam wykres można narysować dużo prościej bez wyliczania wartości f odwrotnej. Jak?)
  newton.m - prosta implementacja skalarnej metody Newtona
• Lab 3 (30-10-2009)- metody rozwiązywania układów równań nieliniowych i obliczenia w fl.
  • Zad 1 Zaimplementować wielowymiarową metode Newtona w wersji z dokładnym jakobianem i w wersji gdy jakobian przybliżany różnicami dzielonymi z parametrem h (jak w zad 5 z poprzednego labu). Zastosować do rozwiązania układu f1(x,y)=x+2y=1; f2(x,y)=3*x^2+y^2=1 dla różnych przybliżeń początkowych.
  • Zad 2 Zastosuj metode iteracji prostych x_n=x_n-aF(x) do układu z zad 1 tzn F=[f1;f2] a a parametr!=0 np a=1,-1,1e-1 itp Czy zawsze zbiega i jeśli tak to jak szybko?
  • Zad 3 Do układu równań z zad 1 zastosować fsolve() - porównać z Pańśtwa implementacją metody Newtona!
  • Zad 4 Wyznaczyc epsilon maszynowy czyli najmniejszą liczbę fl taką że po dodania jej do jeden dostajemy coś większego od jeden (w fl oczywiście) 2^{-t} dla t - ilości bitów mantysy. Można to zrobić i w C/C++. Czy wyszło to samo co w octavie (dla liczb typu double)?
  • Zad 5 Obliczyć wartości funkcji f1(x)=(x-2)^4 i f2(x)=x^4-.....+16 na siatce równomiernej na [2-a,2+a] w wektorze x dla a=1e-3. narysować wykresy obu funkcji i policzyć błąd ||f1(x)-f2(x)||_\inf czyli max_k |f1(x(k)-f2(x(k))|. (Oczyiwście pomijając błędy fl f1=f2 - dlatego istotne jest aby liczyć funkcje wprost z tych wzorów!)
  • Zad 6 Zastosować bisekcję dla funckji (x-2)^3 liczonej z wzoru na rozwinięcie dwumianu x^3-...-8 startując z a=2-1e-3 a b=1+1e-3 z warunkiem stopu że błąd jest mniejszy od 1e-20 (czy długości odcinka w którym jest rozwiązanie ma być mniejsza od 2e-20).
  • Zad 7 Policzyć przybliżenie exp(x) z rozwinięcia w szereg \sum_{k=0}x^k/k! (biorąc 100 a potem 1000 elementów szeregu) sprawdzić błąd względny |y-exp(x)|/|exp(x)| dla x od -100 do 100 (np dla -100 -80 ..-10 0 10 .. 80 100). Czy błedy dla liczb ujemnych i dodatnich są tego samego rzędu? Jak to można łatwo poprawic?
  • Zad 8 Policzyć całki I_n=\int_0^1x^n/(5+x)dx n=0,..,20 z wzoru I_n+5*I_{n-1}=1/n raz biorąc I_0=ln(6/5) i iterując wprzód tzn dla rosnących n a raz biorąc I_30=1/180 i iterując w tył. (Oczywiście dla n=30,29,..) będzie kiesko ale okaże się że potem jest coraz lepiej.
• Lab 4 (13-11-2009) Rozkład LU. Uwarunkowanie zadania rozwiązywania url.
  • Zad 1 Stworzyć w octavie macierz trójdiagonalną A z 2 na diagonali i -1 na pod- i nad diagonalą tzn A_{i, i}=2 dla i=1,...,n i A_{i, i+1}=A_{i+1, i}=-1 dla i=1,..,n-1 a pozostałe elementy zero, zarówno przy pomocy funkcji octave'a diag() jak i wprost w pętli. Porównać czas używając tic i toc dla n=10,100,1000. Policzyć uwarunkowanie macierzy dla różnych N , policzyć macierz odwrotną przy pomocy inv() (czy też jest trójdiagonalna?).
  • Zad 2 Zaimplementować rozkład LU dla macierzy trójdiagonalnej bez wyboru elementu głównego. Przy czym do funkcji macierz przekazywac jako trzy wektory z przekątnymi. Zastosować do rozwiązania układu równań z macierzą z poprzedniego zadania. Porównać czas rozwiązywania Państwa algorytmem a algorytmem octave'a (operator '\'). Rozwiązać układ dla losowego rozwiązania x=(1,\ldots,1)^T dla N=10,100,1000 czy nawet 10000 . Policzyć błąd residualny (dla znanego losowego rozwiązania) i błąd rzeczywisty w różnych normach tzn \|x-\tilde{x}\|_2 i \|f-A*\tilde{x}\|_2 gdzie f=A*x a \tilde{x} rozwiązanie obliczone przez octave'a (A\f) czy Państwa algorytm czy w najgorszy sposób A^{-1}*f (wczesniej musimy policzyć macierz odwrotną) .
  • Zad 3 Zaimplementować rozkład LU dla macierzy trójdiagonalnej z częściowym wyborem elementu głównego. I dalej testować jak w zadaniu 2. Zastosować do rozwiązania układu równań z macierzą z poprzedniego zadania A i z macierzą A-1.5*I .
  • Zad 4 W octavie przetestowac eliminacje Gaussa z częściowym wyborem elementu głównego (permutacja wierszy) i bez wyboru dla macierzy A=[e, 1; 1, 1] z e=eps/10 (tu eps funkcja octave'a zwracająca epsilon maszynowy) i wektorem prawej strony f=[1;0]^T.
  • Zad 5 Przetestuj solver octave'a dla macierzy Hilberta H(n) (polecenie octave'a hilb() ją tworzy), tzn stwórz macierz H(N) dla N=10,20,30 dla znanego rozwiązania np. stałego równego jeden na każdej pozycji czyli sol(k)=1, policz H(n)*sol=f - rozwiąż w octavie uklad z H(n) i f i policz normy (1,2 itd) błędu residualnego ||H(n)x-f|| i błędu ||x- sol|| - (w praktyce jesteśmy w stanie tylko obliczać błąd residualny bo zazwyczaj nie znamy sol!).
  • Zad 6 Zaprogramować rozkład Cholskiego dla A=A^T>0 trójdiagonalnej tzn policzyć L dolnotrójkątną z jedną poddiagonalą (czyli reprezentowaną przez dwa wektory z przekątną i podprzekątną) taka że A=L^TL . Zastosować do rozwiązania układu równań z zad 1 i 2. Testować jak w zadaniu 2.
• Lab 5 (27-11-2009) macierze rzadkie (sparse) w octavie czyli polecenie sparse()i inne (np spdiag()). LZNK i metoda Householdera.
  • Zad 1 Zaprogramuj funkcję octave'a " function y=H(x,w,inw)" która dla danych wektorów x i w tego samego wymiaru N i skalaru inw=1/||w||_2 obliczy H_w *x dla H_w=I-2*inw*w *w^T czyli działanie przekształcenia (macierzy) Householdera na wektorze x. (skalar może być parametrem opcjonalnym - jeśli funkcja będzie wywołana z dwoma tylko parametrami co można sprawdzić (nargin<3) to tę normę można obliczyć w tej funkcji). Przetestować dla losowych wektorów x i w , sprawdzić czy ||H_w *x||_2=||x ||_2. Można napisać ogólniejszą wersje gdzie x=X macierz N X M i wtedy wynik H*X (jak to zrobić bez użycia pętli?). Sprawdzić czy ||A||_2=||H*A||_2=||A*H||_2 i ||A||_F=||H*A||_2=||A*H||_F dla losowej macierzy A i H przekształcenia Householdera dla losowego wektora w<>0 wykorzystując tę funkcję.
  • Zad 2(testujemy twierdzenie z wykładu) Znajdź wektor Householdera w taki że odp przekształcenie Householdera przeprowadza dany wektor u<>0 na dany wektor l* v <>0 dla pewnej dodatniej stałej l. Przetestuj dla dowolnych 2 różnych wektorów czy rzeczywiście tak jest! Tzn czy H_w u= l v.
  • Zad 3 Zastosuj octave'a do rozwiązania zadania znalezienia krzywej zadanej równaniem ax*x+b*y*y=1 najlepiej pasującej do danych N punktów (x_k,y_k) (możemy np np wziąć lekko zaburzone punkty z danej elipsy y_k=\sqrt{1-4*x_k*x_k} + zab_k gdzie zab_k zaburzenie wylosowane z [0,1e-2] czyli mamy N punktów x_k np. x_k=1/k czy x_k=-1+h*k dla h=2/N k=1,...,N i z powyższego wzoru wyliczamy y_k +zab_k. Czy obliczone a i b bliskie 4 i 1?
  • Zad 4 Zastosuj metodę Householdera do rozwiązania zadania znalezienia prostej najlepiej przybliżającej N punktów punkty (k,1+2*k+ los_k) dla k=1,...,N gdzie los=( los_1,..., los_N) losowy wektor za zakresu [- los, los] testować dla los= 1,1e1,1e2,1e3] . (Funkcja rand(n) generuje wektor losowy o rozkładzie jednostajnym na [0,1] w octavie). Porównaj z wynikami otrzymanymi za pomocą standardowej funkcji octave'a tzn \ oraz korzystając z funkcji qr(A).która zwraca rozkład QR macierzy (Q ortogonalna!).
  • Zad 5 Zaprogramuj wersję metody Houlsholdera rozwiązywania Ax=b dla A macierzy trójdiagonalnej różnych wymiarów z 2 na diagonali a -1 na pod- i naddiagonalą i jakimś losowym wektorem rozwiązania tak aby koszt znalezienia rozkładu i rozwiązania układu równań liniowych był jak O(N)?
  • Zad 6 Stwórz macierz rzadką L_N o 2 na diagonali i -1 na pod- i nad diagonali w formacie rzadkim (nie tworząc wpierw macierzy pełnej!) i rozwiąż ją dla N=1000,1e4,1e5 z wektorem z prawej strony f=(2,-1,0...,0)^T . Porównaj czas z solverem octave'a czy solverem Państwa zastosowanym do takiej macierzy pełnej (poprzedni lab!).
• Lab 6 (11-12-2009) metody iteracyjne rozwiązywania układów równań liniowych i zadania własnego czyli metoda potęgowa i odwrotna potęgowa.
  • Zad 1 Zaimplementuj metody iteracyjne Jakobiego i Gaussa-Seidla i zastosuj je do macierzy L_{N,a} trójdigonalnej 2+a na diagonali -1 na nad- poddiagonali (W przypadku m. Jakobiego można porównać dwa podejścia - z macierzami rzadkimi i funkcjami które wprost w tym przypadku zaimplementują te metody dla tej danej macierzy) Wektor prawej strony jak w poprzednim zadaniu. Testować zbieżność dla różnych wartości a<>-2 np a=0,-1,0,1,10 w normie max i drugiej dla x_0=0. Można a zastąpić przez wektor a_k np [0,1,...,N] tzn na pozycji (k,k) mamy 2+k czy przeskalować macierz L_{N,0} przez macierz diagonalną z diagonalą [1,...,N] czy [1, 2^2,...,N^2].
  • Zad 2 Metoda Richardsona dla macierzy trójdiagonalnej (jak w zad 1 dla a=0 i a=10) i dla macierzy 2x2 symetrycznej dodatnio określonej - zbadać zbieżność dla różnych wartości parametru. Dla N=100 znaleźć przy pomocy eig() wartości własne i sprawdzić czy szybkość zbieżności zgadza się z teorią dla x_0=0. (Wylosować x^* i wziąść b=A*x^*)
  • Zad 3 Metoda potęgowa dla macierzy z poprzedniego zadania - warto przy okazji sprawdzić jak działa solver octave'a dla zadania własnego (eig() i eigs()). Można też przetestować metodę dla jakiejś symetrycznej macierzy losowej (losujemy A a potem A=(A+A')/2). Jako warunek stopu |r_k-r_{k+1}|
  • Zad 4 Metoda odwrotna potęgowa dla macierzy z zad 1 i ewnt. losowych symetrycznych. Czyli znajdujemy najmniejszą wartość własną. Można też odwrotną potęgową z parametrem.
  • Zad 5 Metoda QR. Zaprogramuj nastepującą metode dla A_0=A=A^T: A_k=QR (Q - ortogonalna, R górnotrójkątna) A_{k+1}=RQ i warto wymnażać macierze Q. Teoria mówi że dla A symetrycznej ciąg A_k zbiegnie do macierzy diagonalnej (i A_k są podobne do A czyli mają te same wartości własne). Przy okazji możemy wyliczyć macierz wektorów własnych (jak?). Jako warunek stopu weźmy to że suma kwadratów elementów poza diagonalą jest mniejsza od zadanej tolerancji np. \sum_{k,l;k<>l} A_{k l}^2
  • Zad 6 Wykorzystując qr() octave'a lub swoją procedury tworzące odp. macierz Householdera napisz procedure znajdującą dla danej macierzy kwadratowej A taką Q ortogonalą i B 3-diagonalną że QBQ^T=A.
• Lab 7 (8-01-2009)BLAS/LAPACK czyli kilka przykłądów funckji z bibliotek numerycznych i interpolacja wielomianowa Lagrange'a (ewent. splajnowa)
  • Zad 1 BLAS Level 1. Napisz możliwie prosty program w C (ewnt. C++) obliczający wpierw iloczyn skalarny dwu wektorów długości trzy x=[1 1 1]^T i y=[1 0 1]^T oraz kwadrat normy drugiej y (y^T*y)(czy jakichkolwiek innych) z wykorzystaniem funkcji z BLAS Level 1 DDOT(). A następnie ortogonalizyjący te wektory tzn obliczający x=x-(x^Ty)*y/(y^T*y) z wykorzystaniem DAXPY() A tu przykładowy program wykorzystujący funkcje DNRM2() i DSCAL() do unormowania wektora: blas-simplest.c
  • Zad 2 BLAS Level 2. Napisz program mnożący macierz A=[1 2;-2 3] i potem A^T przez wektor x=[1;1] z wykorzystaniem funkcji z BLASów DGEMV() .
  • Zad 3 LAPACK. Napisz program znajdujący rozwiązanie AX=F dla macierzy A=[1 2;-2 3] i F=[3;1] z wykorzystaniem funkcji z lapacka DGESV(). A następnie niech ten program obliczy A^{-1} z(z wykorzystaniem tej samej funkcji!)
  • Zad 4 LAPACK. Napisz program znajdujący rozwiązanie AX=F dla macierzy trójdiagonalnej NxN symetrycznej dodatnio określonej która ma 2 na diagonali i -1 na pod- i naddiagonali z wykorzystaniem funkcji z lapacka DPTSV(). Przetestuje dla F=[2 -1 0...0]^T (roz. pierwszy wersor) i N=10,100,1000 etc
  • Zad 5 Interpolacja lagrange'a- oblicz wykorzystując funkcję polyfit() wiel. interpolacyjny dla funkcji sin() dla N węzłów równoodległych i N węzłów Czebyszewa na [-pi,2*pi] dla N=5,10,20,30. Oblicz dyksretną normę max na siatce 1000 punktowej na tym odcinku tzn \max|sin(x_k)-L_N sin(x_k)| x_k pkt siatki. Następnie narysuj wykresy sin i tych wielomianów dla różnych N.
  • Interpolacja L przykład Rungego - Wykonaj poprzednie zadanie ale dla f(x)=1/(1+x*x) na [-5,5]. Czy ciąg wielomianów interpolacyjnych zbiega do f?
• Lab 8 (22-01-2009) Zaliczenie projektu zaliczeniowego. Kwadratury i mnożenie wielomianów z użyciem DTF i odwrotnej DTF w octavie.
  • Zaprogramuj w octave złożoną kwadraturę prostokątów $P_N f=\sum_{k=1}^N f((k-0.5)*h)*h$ dla $h=(b-a)/N$. Porównaj tę kwadraturę z wynikami funkcji quad() np dla $exp(-x*x/2) czy sin(k*x)$ na różnych odcinkach i różnych k np dla odcinka $[0,10*pi] i k=2,8,64$ itp Zbadaj rząd eksperymentalnie dla całki $sin(x)$ na $[0,\pi] $równej dwa tzn Oblicz $P_N \sin$ dla N=10,20,40, itd i potem dla $E_N=|2-P_N|$ obliczaj $E_N/E_{2N}$ - czy wychodzi w przybliżeniu cztery? Przetestuj tak samo tę kwadraturę dla funkcji \sqrt{x} na [0,1] (tu rzeczywiście brak regularności jako że pochodna blisko zera zbiega do nieskończoności!) Jaki rząd tej kwadratury dostajemy z eksperymentu?
  • Treść, jak powyżej ale dla złożonej kwadratury trapezów i Simpsona. (ewent. do domu)
  • Korzystając z funkcji fft() i ifft() napisać procedure wymnażającą dwa wielomiany $w(z)=\sum_k a_k z^k i p(z)=\sum_k b_k z^k$ st<=N przez siebie tzn. jako argumenty dwa wektory $(a_k) i (b_k)$ ze współczynnikami wielomianów (N może być podane albo wzięte jako max z długości obu wektorów) a na wyjściu wektor $(c_k)$ ze współczynnikami wielomianu $r(z)=w(z)*p(z)=\sum_k c_k z^k$ stopnia <=2N.
  • (do domu - wymaga przeczytania rozdziału z książki) Zaprogramuj kwadraturę Romberga (jak opisano w książce Kincaid,Cheney Analiza numeryczna)

Kilka przykładowych skryptów, m-plików (plików funkcyjnych ) octave'a, źródeł prostych pogramów w C, źródłowych funkcji z blasów czy LAPACKa czy użytecznych linków
(dodawanych w miarę postępu labu)

Skrypty m-pliki octave'a

ZAchęcam żeby poza projektem do napisania funckji z tą metodą w C++ i potem skompilowania jako plik oct octave'a czyli dołączenia do dystrybucji octave'a (octave-forge).

---

Powrót do mojej strony domowej.

---

Ostatnia aktualizacja: 8 marca 2010