Semestr zimowy 2011/12

1. Metody numeryczne (dla informatyków)
Ćwiczenia/Lab (gr 3)
• Program bieżącego labu
• Program bieżących ćwiczeń
2. Matematyka obliczeniowa 2

Metody numeryczne (dla informatyków)

Ćwiczenia/Lab MN

Zaliczenie ćwiczeń

Serie zadań domowych - projekty zaliczeniowe z labu

• Pierwsza seria - plik pdf termin 3cie ćwiczenia tj. 9 listopada 2010
• Pierwszy projekt - plik pdf termin czwarty lab tj. 16 listopada 2010
• Druga seria - plik pdf termin 5te ćwiczenia tj. 7 grudnia 2010
• Drugi projekt - plik pdf termin ostatni lab tj. 18 stycznia 2010 (o ile nie będzie jakiegoś przesunięcia)
• Trzecia seria - plik pdf termin ostatnie ćwiczenia tj. 11 stycznia 2011

Program ćwiczeń

• Ćwicz 1 (12-10-2011) Skalarne metodu rozwiązywania równań nieliniowych: m. Newtona
  1. Metoda bisekcji - formalny dowód zbieżności.
  2. Metoda Herona (czyli m Newtona dla równania x*x-a=0; a>0), że zbiega dla dowolnego x0; (elementarnie)
  3. Odwracanie bez dzielenia (m. Newtona dla 1/x-a=0 dla jakich x_0 m. Newtona zbiega?).
  4. Zera wielokrotne: pokaż że dla zera 2-krotnego met. Newtona jest zbieżna lokalnie np. dla f(x)=(x-2)^2=0. Uogólnij dla dowolnej funkcji takiej, że f(x)=f'(x)=0, f''(x) różna od zera - wtedy zachodzi lokalna zbieżność Newtona liniowa (o ile w kolejne iteracje dobrze zdefiniowane).
  5. Policz wzór na metodę Newtona dla funkcji g(x)=f(x)/f'(x). Policz pierwszą pochodną tej funkcji w a dla f(a)=f'(a)=0 z f''(a) nie zerowym.
  6. Pokaż, że jeśli f(a)=0=f''(a) a f'(a) nie zero 0 to metoda Newtona zbiega lokalnie kubicznie.
  7. (metoda Halleya) Pokaż że druga pochodna funkcji g(x)=[f(x)]/[\sqrt{f'(x)}] wynosi zero dla x=a zera funkcji f o ile f'(a)>0. Pokaż, że wzór na tę metodę to x_{n+1}=x_n - [2f(x_n)f'(x_n)]/[2(f'(x_n)^2 - f''(x_n)f(x_n)].
  8. Zbadaj zachowanie iteracji metody Newtona dla funkcji f(x)=sgn(x)\sqrt{|x|} dla x_0 nie zero.
  Zrobiliśmy zadania 1-4 na ćwiczeniach. Pozostałe albo na następnych ćwiczeniach albo do domu.
• Ćwicz 2 (26-10-2010) Metodu rozwiązywania równań i układów równań nieliniowych cd. Normy wektorów i macierzy. Dokończenie niektórych zadań z poprzedniej serii.
  1. Jeśli f(x^*)=f''(x^*)=0\not = f'(x^*) metoda Newtona jest zbieżna lokalnie kubicznie.
  2. Metoda iteracji prostych- pokaż że x_n=cos(x_{n-1}) jest zbieżne dla dowolnego x_0 . Jakie musi być n aby otrzymać przybliżenie x_n liczby x^*=cos(x^*) spełniające |x_n-x^*| < 1e-4 .
  3. Metoda Newtona wielowymiarowa. Sprawdź czy dla funkcji f1(x,y)=x*x+y*y-1=0; f2(x,y)=x+2y=0 m. Newtona będzie lokalnie zbieżna w otoczeniu zer tych funkcji. Policz pierwszą iteracje dla (x_0,y_0)=(1,0) . Napisz w pseudokodzie implementacje wiel. met. Newtona, zakładając, że dysponujemy funkcją x=linsolve(A,b) rozwiązującą układ równań liniowych Ax=b .
  4. Optymalne stałe równoważności dla norm p-tych dla p=1,2,\infty .
  5. Pokaż wzór na macierz indukowaną maksimum. Do domu na normę indukowaną pierwszą oraz wzory na note normy dla macierzy zespolonych.
  6. Norma Frobeniusa nie jest normą indukowaną. Z góry szacuje normę drugą. Znajdź stałą równowążności taką, że norma druga szacuje normę Frobeniusa.
  7. Q ortogonalna. Pokaż, że \|QA\|_2=\|AQ\|_2=\|A\|_2 oraz to samo dla normy Frobeniusa.
• Ćwicz 3 (9-11-2010)
  1. Znajdź współczynnik uwarunkowania zadania obliczania iloczynu A*x ( A macierz, x wektor) ze względu na zaburzenie A .
  2. Rozwiąż dwa układy A_kx_k=b z A_1=[1,1;1,0.98] i A_2=[1,1;1,0.99] i b=(1,2)^T . Policz współczynnik uwarunkowania A_1 w normie maksimum.
  3. Eliminacja Gaussa dla macierzy trójdiagonalnej bez permutacji wierszy/kolumn i (do domu) z wyborem elementu głównego w kolumnie (permutacje wierszy).
  4. Rozpatrzmy macierzy w postaci blokowej B=[A, b;c^T a] dla A macierzy mxm nieosobliwej której rozkład LU znamy, b,c wektorów wymiaru m i danego skalaru a. Jak możliwie niskim kosztem rozwiązać url z B (o ile B nieosobliwa)? Podaj warunek na to aby B byla nieosobliwa.
  5. Rozpatrzmy macierzy w postaci blokowej G=[A, B;B^T 0] dla A macierzy mxm symetrycznej dodatnio określonej której rozkład LU znamy, B wymiaru mxn o kolumnach tworzących układ liniowow niezależny. Czy ta macierz jest nieosobliwa? jeśli tak to jak układ Gx=b rozwiązać, przy założeniu że potrafimy rozwiązać układ postaci AC=B z C macierzą np. kosztem O(m^3) .
  6. Udowodnij, że dla m. silnie diagonalnie dominującej wierszowo (|A_{k k}| > \sum_{k \not = l} |A_{k l}| \ \forall k ) LU można wykonać bez wyboru elementu głównego (czyli bez permutacji kolumn czy wierszy).
  7. Treść jak w poprzednim zadaniu ale dla macierzy symetrycznej dodatnio określonej. Wsk: wykorzystaj fakt, że jeśli A=A^T>0 to CAC^T też symetryczna i dodatnio określona o ile C nieosobliwa.
  8. Udowodnij jednoznaczność rozkładu LU o ile istnieje (L dolnotr. z 1 na przekątnej, U górnotr. z elementami dodatnimi na przekątnej)

• Lab 1 (5-10-2011) - 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 Przetestuj metodę bisekcji z poprzedniego zadania do rozwiązania innych problemów: x^2-2=0, exp(x)=2, cos(10*x)=0 startując z odcinka [0,110] - przy dużej ilości zer na danym odcinku znalezione zero jest w sumie przypadkowe. Zad 5 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ładność bezwzględna (błąd |x_n-x^*|=2*|b_n-a_n| 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 maksymalnej ilości iteracji, 2 - wartości w końcach odcinka funkcji mają ten sam znak. Funkcja ma działać również jak podamy tylko trzy czy cztery pierwsze argumenty - i wtedy max. ilość iteracji domyslnie ma być 100 , a epsilon 1e-5 (jak tylko 3 argumenty). Zad 6 znajdź w manualu funkcję octave'a fzero() - rozwiązująca równania skalarne 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 (19-10-2011) - metody rozwiązywania równań nieliniowych
  1. Zad 1 Zaimplementować metodę Newtona w octavie i przetestować jej zbieżnośc dla następujących funkcji: x*x-2 z x0=2, x*x*x-27 z x0=27, exp(x)-2 z x0=10,-10,100, sin(x) dla x0=2, (x-2)^k dla k=2,4,8,16 z x0=3, x*x-2 z x0=1e6, 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.
  2. 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.
  3. Zad 3 Sprawdzić czy metoda iteracji prostych x_n=cos(x_{n-1}) zbiega do x=cos(x). Zbadac czy zbieżność jest liniowa.
  4. Zad 4 Porównać wyniki z zad 1 z funkcjami octave'a fsolve() i fzero() (pomoc: help fsolve/fzero).
  5. 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ładną metodą Newtona (szczególnie ostatnie iteracje) dla funkcji z zadania 1 i z metodą siecznych.
  6. 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?
  7. 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 funkcji odwrotnej. Jak?)
• Lab 3 (2-11-2011) - wielowymiarowe układy równań nieliniowych i obliczenia w fl.
  1. Zaimplementować wielowymiarową metodę Newtona w wersji z dokładnym Jakobianem i (do domu) w wersji gdy Jakobian przybliżany różnicami dzielonymi z parametrem h=1e-8. 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.
  2. Wyznaczyć 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. Porównać z eps! komendą octave'a. Można to zrobić i w C/C++. Czy wyszło to samo co w octavie (dla liczb typu double)?
  3. Powtórzyć to zadanie w pojedyńczej precyzji tzn. wymuszając żeby zmienna była typu single (polecenie x=single(x) konwertuje do poj. precyzji)
  4. 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) \|_\infty czyli \max_k |f1(x(k)-f2(x(k))|.
  5. (do domu) Zastosować bisekcję dla funkcji (x-2)^3 liczonej z wzoru na rozwinięcie dwumianu x^3-...-8 startując z a=2-1e-3 a b=1+2e-3 z warunkiem stopu że błąd jest mniejszy od 1e-20 (czy długości odcinka w którym jest rozwiązanie). Narysować wykres tej funkcji liczone z obu wzorów na odcinkach 2+[-h,2*h] dla h=10^{-3},...,10^{-5}.
  6. Policzyć przybliżenie \exp(x) z rozwinięcia w szereg \sum_{k=0}^\infty 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 liczb różniących się o dziesięć czyli -100, -80,...,-10, 0, 10,...,80, 100. Czy błędy dla liczb ujemnych i dodatnich są tego samego rzędu? Jak to można łatwo poprawić?
  7. (do domu) Policzyć całki I_n=\int_0^1x^n/(5+x)d x dla n=0,..,20 z wzoru I_n+5*I_{n-1}=1/n raz biorąc I_0=\log(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 kiepsko ale okaże się że potem jest coraz lepiej.
• Lab 4 (16-11-2011) - bezpośrednie metody rozwiązywania układów równań liniowych
  1. Przetestuj solver octave'a tzn operator backslash \ dla prostego układu równań liniowych z macierzami nieosobliwymi np A1=[1,1;1,0.98] i b=[2;1] , A2=[1,1;1,0.99] i b=[2;1] oraz osobliwą A=[-2,2;-1,1] z tym samym b . Sprawdź czy otrzymane rozwiązanie x=A\ b! jest rzeczywiście rozwiązaniem tzn policz normę A*x-b . Porównaj rozwiązania ukłądów z A1 i A2 . policz macierze odwrotne do A1 i A2 i uwarunkowania tych macierzy.
  2. Sprawdzić dla konkretnej macierzy 3x3 np A=[1, 2, 3;-1, 0, 1;1, 1, 1] czy odpowiednie macierze górno/dolno trójkątne, permutacji czy ortogonalne uzyskane przy pomocy funkcji octave'a lu() i qr() rzeczywiście zwracają odpowiednie macierze rozkładu LU (PA=LU) lub rozkładu QR (A=QR) np policzyć pierwszą normę różnicy A-QR dla rozkładu QR.
  3. 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 , czyli policz H(n)*sol=f - rozwiąż w octavie układ z H(n) i f i policz normy (1,2 itd) \|H(n)x-f\| i \|x- sol\| dla x rozwiązania które wyliczył octave. Sprawdź czy macierz uzyskana przez inv(hilb(n)) i invhilb(n) są rzeczywiście odwrotne do hilb(n) - policz normę macierzową różnic tzn. np. inhilb(n)*hilb(n)-eye(n)
  4. 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?). Policzyć rozwiązanie układu równań liniowych Ax=b dla różnych znanych x za pomocą y=inv(A)*b i w=A\b (wczesniej obliczymy b=A*x). Porównać normy Ay-b,x-y,w-x (pierwszą, drugą, supremum) dla N=10,100,1e3,1e4 itp
  5. Zaimplementować rozkład LU dla macierzy trójdiagonalnej bez wyboru elementu głównego. Przy czym do funkcji macierz przekazywać 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 ( \ ) Rozwiązać układ dla losowego rozwiązania x=(1,...,1)^T dla N=10,100,1000 czy nawet 10000 . Policzyć błąd rezydualny (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 czy Państwa algorytm czy w najgorszy sposób A^{-1}*f (wcześniej musimy policzyć macierz odwrotną) .
  6. (do domu) Zaimplementować rozkład LU dla macierzy trójdiagonalnej z częściowym wyborem elementu głównego. I dalej testować jak poprzednim zadaniu. Zastosować do rozwiązania układu równań z macierzą z poprzedniego zadania A i z macierzą A-1.5 .
  7. W octavie przetestować eliminacje Gaussa z częściowym wyborem 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 .
  8. Zaprogramować rozkład Choleskiego 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 macierzą z 2 na diagonali i -1 na pod i nad diagonalach.
• Lab 5 (30-11-2011) Macierze rzadkie i iteracyjne rozwiązywanie równań liniowych.
  1. Utwórz macierz trójdiagonalną T_N z a na diagonali i -1 na pod-nadiagonalach z wykorzystaniem funkcji sparse() ale nie tworząc macierzy pełnej tzn podając tylko elementy różne od zera.
  2. Zapoznaj się z funkcją pcg() octave'a. Przetestuj metod cg (funkcja octave'a pcg()) dla macierzy wymiary 25 losowych symetrycznych dodatnio okreslonych: B=rand(N,N), B=B'*B+pI dla różnych wartości p=0,5,10,20,100 tzn na ekranie rysuj wektory zwracane w RESVEC przez pcg() dla różnych wartości p.
  3. Przetestuj metod cg (funkcja octave'a pcg()) dla macierzy wymiary T_N+pI dla różnych wartości p=0,5,10,20,100 i N=10,20,50, tzn na ekranie rysuj wektory zwracane w RESVEC przez pcg() dla różnych wartości p i N.
  4. Zaimplementuj w octave metody Jakobiego i Gaussa Seidla - przetestuj ich zbieżność dla układu równań Ax^*=b dla macierzy A=B+pI gdzie B losowa a p dodatni parametr np. p=[0,5,10,20,30,100] . Tzn proszę dla każdej metody stworzyć wektor (\|b-Ax_k\|_2)_k i następnie narysować wykresy obu wektorów w jednym oknie. Czy widać różnice w szybkości zbieżności obu metod?
  5. Rozpatrzmy macierz 60x60 o wartościach własnych: 1,99,100 - można ją utworzyć np jako A=Q\Lambda Q^T dla Q macierzy ortogonalnej np. macierzy Householdera i \Lambda diagonalnej o wartościach na diagonali równym 1,99,100 Przetestuj zbieżność metody CG. Nastepnie dla znanego rozwiązania Ax^*=b dobierz x^0 tak aby x^*-x^0 ortogonalne do pierwszej kolumny Q np x^0=x^*+Q[0;1;1] i powtórz test cg..
  6. Testy metody Richardsona. Zaimplementuj metodę Richardsona i przetestuj jak w poprzednim zadaniu ale w normie drugiej dla różnych wartości parametru \tau . Oblicz w octavie maksymalną i minimalną wartość własną A i przetestuj czy zbieżność w normie drugiej jest taka jak w teorii. (dla znanego rozwiązania)
• Lab 6 (14-12-2012) - LZNK, metoda Householdera, rozkład QR
  1. Testujemy funkcję qr() octave'a do znalezienia rozkładu macierzy A=[1, 1; 1, 0; 1, -1]. Następnie stosujemy ten rozkład do znalezienia rozwiązania LZNK z A i wektorem prawej strony b=[1;1;1] (powinno wyjść x=[1; 0; 0]), porównujemy z rozwiązaniem operatorem backslash. Stwórz macierz A^TA oraz wektor A^Tb - rozwiąż układ A^TAy=A^Tb - czy y==x? (module błędy zaokrągleń)?
  2. Porównaj 2 metody rozwiązania LZNK dla wektora b i macierzy A m x n takiej, że (A)_{k l}=(x_l)^{k-1} dla k=1,..,m, x_l=l*h l=1,..,n h=1/m. Pierwsza - operator backslash z octave'a, druga tworzymy ukłąd rółnań normalnych A^TA i A^b i rozwiązujemy backslashem układ A^TAy=A^Tb. Przetestować dla m =2*n z n=10,20,30 dla b - pierwsa kolumna A. (Rozwiązanie sol - wersor). Policzyć normę różnicy x_k-sol dla x_k rozwiązania policzonego k-tą metodą. Wsk: Funkcja vander tworzy taką macierz dla m=n (tylko z odpowienio przepermutownymi kolumnami/wierszami) można stworzyć macierz mxm i wyciac odpowiednią podmacierz.
  3. Zastosuj octave'a do rozwiązania zadania znalezienia krzywej zadanej równaniem a x*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 ?
  4. Zaprogramuj funkcję octave'a \begin{lstlisting} function y=H(x,w,nw) \end{lstlisting} która dla danych wektorów \vec{x} i \vec{w} tego samego wymiaru N i skalaru nw=\|w\|_2 obliczy H_w \vec{x} dla H_w=I-2*\frac{1}{nw}\vec{w}\vec{w}^T czyli przekształcenia (macierzy) Householdera. (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 \vec{x} i \vec{w} , sprawdzić czy \|H_w \vec{x} \|_2=\|\vec{x} \|_2. Można napisać ogólniejszą wersje gdzie x=X macierz N\times 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\not=0 wykorzystując tę funkcję.
  5. (testujemy twierdzenie z wykładu) Znajdź wektor Householdera \vec{w} taki że odp przekształcenie Householdera przeprowadza dany wektor \vec{u}\not=0 na dany wektor l*\vec{v}\not=0 dla pewnej dodatniej stałej l . Przetestuj dla dowolnych 2 różnych wektorów czy rzeczywiście tak jest, tzn. czy H_w \vec{u}= l \vec{v} .
  6. Zastosuj metodę Householdera do rozwiązania zadania znalezienia prostej najlepiej przybliżającej N punktów punkty (k,1+2*k+\epsilon_k) dla k=1,...,N gdzie \epsilon(\epsilon_1,...,\epsilon_N) losowy wektor za zakresu [-\epsilon,\epsilon] testować dla \epsilon= 1,1e1,1e2,1e3] . (Funkcja \lstinline+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).
  7. Zaprogramuj metodę Householdera rozwiązywania A\vec{x}=\vec{b} dla A macierzy trójdiagonalnej różnych wymiarów np. 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) ?
• Lab 7 (4-01-2011) - interpolacja, splajny
  1. Interpolacja Lagrange'a: Sprawdź jak działą funkcja polyfit(). Znajdź wielomian interpolacyjny na 4 węzłach -1,0,1,2 dla wielomianów x^3-1 i x^4-2, narysuj wykres wielomianu interpolacyjnego i tej funkcji.
  2. Wykorzystując funkcję polyfit() znajdź wielomiany interpolacyjne 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 dyskretną normę max na siatce 1000 punktowej na tym odcinku, tzn. \max|\sin(x_k)-L_N \sin(x_k)| , gdzie x_k pkt siatki. Następnie narysuj wykresy \sin(x) i tych wielomianów dla różnych N - np. używając funkcji polyval() . Węzły Czebyszewa to zera T_{n+1}(t)=\cos([n+1]\mathrm{arccos}(t)) na [-1,1] odpowiednio przesunięte i przeskalowane.
  3. Powtórz poprzednie zadanie, ale dla \log(1+x) na odcinku [0,10] .
  4. Funkcje octave'a spline() and ppval. Zapoznaj się z tymi funkcjami ( help spline i help ppval). Wykorzystując te funkcje narysuj wykres splajnu kubicznego s_1 na podziale równomiernym odcinka [-3,3] z węzłami \{x_k=k\} dla k=-3,-2,...,3 - ręcznie podając losowe wartości przyjmowane przez te splajny w tych węzłach: np.: s_1(x_k)=(-1)^k , czy odpowiednio s_1(x_k)=1 . Następnie narysuj wykres splajnu s_2 o tych samych wartościach w węzłach ale podając dodatkowo wartości pochodnych w końcowych węzłach równe zero, tzn. wywołaj funkcje spline() podając dwie wartości więcej. Czy otrzymaliśmy te same splajny? Policz przybliżoną normę różnicy s_1-s_2 na odcinku [-3,3] .
  5. Splajn bazowy interpolacyjny - Dla danych węzłów równoodległych \{k\}_{k=-5,-4,...,5} na [-5,5] narysuj wykres splajnu kubicznego takiego, że s(0)=1 i s(k)=0 dla pozostałych węzłów k\not =0 oraz któtry ma pochodne równego zero w węzłach skrajnych tzn. : -5 i 5 . Narysuj jego wykres, jaki jest nośnik tego splajnu? Powtórz zadanie ale dla splajnu typu not-a-knot tzn. nie podając wartości pochodnych w końcach dla funkcji spline().
  6. Splajn kubiczny o minimalnym nośniku. Dla danych węzłów równoodległych \{k\}_{k=-5,-4,...,5} na [-5,5] narysuj wykres splajnu kubicznego takiego, że s(-1)=s(1)=1,s(0)=4 i s(k)=0 dla węzłów k\not \in\{-1,0,1\} oraz ma pochodne równego zero w węzłach skrajnych tzn. : -5 i 5 . Czy poza [-2,2] ten splajn jest równy zero? Policz przybliżone normy maksimum na [-5,-2] i [2,5] dla tego splajnu i narysuj jego wykres.
• Lab 8 (18-01-2011) - interpolacja, splajny - bład, kwadratury złożone - prostokątów, trapezów; zaliczenie projektu nr 2
  1. Interpolacja Lagrange'a przykład Rungego. Interpolacja Lagrange'a: wykorzystując funkcję polyfit() znajdź wielomiany interpolacyjne L_N f dla funkcji f(x)=1/(1+x*x) na [-5,5] dla N węzłów równoodległych i N węzłów Czebyszewa dla N=5,10,20,30 . Oblicz dyskretną normę max na siatce 1000 punktowej na tym odcinku tzn. \max|f(x_k)-L_Nf(x_k)| , gdzie x_k pkt siatki. Następnie narysuj wykresy f(x) i tych wielomianów dla różnych N - np. używając funkcji polyval(). Czy ciąg wielomianów interpolacyjnych zbiega do f ?
  2. Korzystając z funkcji octave'a spline() znajdź współczynniki splajnu interpolacyjnego kubicznego s_N f na N węzłach równoodległych dla f(x)=\sin(x) na odcinku [-pi,2*pi] dla N=5,10,20,40 . Oblicz dyskretną normę maksimum na siatce 1000 punktowej na tym odcinku tzn. E_N=\max|\sin(x_k)-s_N \sin(x_k)| , gdzie x_k pkt siatki, wyprowadź na ekran iloraz E_{2N}/E_N . Następnie narysuj wykresy \sin i tych splajnów dla różnych N . (zadanie pomocnicze - jak obliczyć wartość splajnu otrzymanego z funkcji spline() ?)
  3. Jak zadanie poprzednie ale dla przykładu Rungego czyli odcinka [0,5] i f(x)=1/(1+x*x) . Czy na wykresie widać że splajny zbiegają do funkcji dla dużych N ?
  4. Narysuj wykresy czterech pierwszych wielomianów Czebyszewa: T_k k=0,1,2,3 wygenerowanych przez regułę trójczłonową: T_{n+1}=2 x T_n-T_{n-1} dla n>0 z T_0=1,T_1(x)=x .
  5. Funkcja quad() w octave'ie - policz całkę z \sin(x) na [0,\pi] i [-1,1] .
  6. Policz za pomocą funkcji octave quad() : \int_{-1}^1\frac{1}{\sqrt{1+x*x}}T_n(x)T_m(x)\:d x dla m=2,n=3 .
  7. Zaprogramuj w octave funkcję ze złożoną kwadraturę trapezów: T_n f=h(0.5[f(a)+f(b)]+\sum_{k=1}^n f(a+k*h)) \qquad h=(b-a)/(n+1) - parametry funkcji: wskaźnik do funkcji obliczającej f , a,b,n . Przetestuj dla funkcji f , której całkę znamy dla N_k=2^kN_0 k=0,1,2,... z ustalonym N_0=5 licząc: \frac{E_{2N}}{E_N}, \qquad E_N=|\int_a^b f - T_N f| dla funkcji
     • \sin(x) na [0,\pi] i N_0=5 czyli funkcji analitycznej
     • x^{j+0.5} na [0,1] dla j=0,1,2 czyli funkcji w C^j([0,1])

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

1. I termin - 2010/11
2. I termin - 2009/10
3. II termin - 2009/10

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) pdf.