Metody numeryczne semestr zimowy 2018/19

Metody numeryczne (dla informatyków)

• Program bieżącego labu
• Program bieżących ćwiczeń

MN Ćwiczenia gr 5 pt 830-10 i gr. 6 1215-1345 - sale 3170 i co2 tyg. na zmianę gr.6/5 lab 1015-1145 lab 3041 - terminy wg USOSa

Zaliczenie ćwiczeń i labu

Serie zadań domowych - projekty zaliczeniowe z labu

• Pierwsza seria - układy r.lin., LU, QR, LZNK, normy, metody iteracyjne dla r. liniowych, zadanie własne - plik pdf zd1.pdf
- kartkówka 16 XI 2018 pt w czasie ćwiczeń - bez zadania własnego i zadań oznaczonych jako 'dodatkowe' lub 'bardzo trudne'.
• Druga seria - fl, interpolacja wielomianowa, splajny, kwadratury, aproksymacja - plik pdf zd2.pdf
• Projekt z labu - LAB1.pdf - drugi lab w listopadzie tzn. 23 lub 30 listopada 2018
• Drugi projekt z labu - LAB2.pdf - zaliczenie pierwszy lab w 01-2018; (w treści projektu napisano że ortogonalizację należy przeprowadzić algorytmem Gramma-Schmidta ale jeśli ktoś zastosuje inny algorytm (odpowiednią funkcje octave'a) to też będzie OK tylko trzeba tenże algorytm zrozumieć i potrafić wytłumaczyć dlaczego działa etc

Program ćwiczeń

1. Zadania na rozkład LU i Choleskiego.(5X18)
2. Kontynuacja. LZNK.(12X18)
3. LZNK cd(19X18)
4. Normy wektorowe i macierzowe. Iteracyjne metody dla równań liniowych.(26.X.18)
5. Iteracyjne metody dla równań liniowych cd. Numeryczne zadanie własne. (9.XI.18)
6. Interpolacja wielomianowa lagrange'a i NEwtona. Algorytmy i bład interpolacji.
7. Interpolacja splajnowa - liniowa i kubiczna.
8. Aproksymacja średniokwadratowa tzn w normie zadanej iloczynem skalarnym. Rzut ortogonalny. Wielomiany ortogonalne, reguła trójczłonowa. Wielomiany Czebyszewa jako przykład ciągu wielomianów ortogonalnych.
9. Aproksymacja jednostajna: Węzły czebyszewa jako optymalne węzły interpolacji. Oszacowanie błędu interpolacji dla węzłów Czebyszewa. Alternans. Alterans dla f. wypukłej i aproksymacji liniowej. (11.I.2019)
10. Rozwiązywanie równań nieliniowych w 1 wymiarze. Metoda Newtona, bisekcji i iteracji prostych.
11. Kwadratury. Rząd. Kwadratura interpolacyjna. Błąd kwadratury interpolacyjnej. Wszystko na przykładzie kwadratury prostokątów P(a,b;f)=(b-a)f((a+b)/2) \approx \int_a^b f. Złożona kw. prostokątów, wzór, błąd.

• Lab 1 - 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. (pierwszy projekt to rozszerzenie tego zadania)
  mnbasic.m - skrypt z podstawami octave'a
• Lab 2 - bezpośrednie metody rozwiązywania układów równań liniowych; uwarunkowanie macierzy
  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ładów z A1 i A2 . policz macierze odwrotne do A1 i A2 i uwarunkowania tych macierzy.
  2. Funkcja inv(). Przeczytaj help do funkcji: inv() . Przetestuj dla A=[1,1;1,-1] czy macierz B obliczona za pomoca tej funkcji rzeczywiscie jest macierz odwrotna? Policz normy pierwsza i Frobeniusa ||A*B-I|| oraz ||B*A-I|| . Zastosuj funkcje do rozwiazywania ukladu równan liniowych Ax=f dla znanego x (liczymy f=A*x ). Policz y jako iloczyn B z f i policz blad rezydualny ||Ay-f|| i ||x-y|| w normie pierwszej i drugiej. Powtórz to zadanie dla macierzy N\times N losowej (funkcja octave'a rand() zwraca macierz losowa dla n=10,50,250,1250 ze znanym rozwiazaniem - oszacuj czas przy pomocy funkcji tic i toc . Porównaj czas i bledy w normie pierwszej dla rozwiazania uzyskanego przy pomocy operator octave'a.
  3. 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 i czy Q ortogonalna (policz ||Q'Q-I||_1 i ||QQ'-I||_1) itd.
  4. Przetestuj solver octave'a dla macierzy górnotrójkątnej U1 i U2 dla U1 z ze wszystkimi elementami równymi jeden (można ją utworzyć poprzez U1=triu(ones(n,n)) a U2=2*Id-U1 (1 na diagonali, -1 ponaddiagonalą dla N=10,20,30 dla znanego rozwiązania np. losowego, czyli policz Uk*sol=fk - rozwiąż w octavie układ z Uk i fk i policz normy: \|Uk*x-fk\|_\infty i \|xk- sol\|_\infty dla xk=Uk\fk rozwiązania, które wyliczył octave.
     Używając inv() policz macierze odwrotne do U1 i U2. Policz uwarunkowania obu macierzy w normie supremum. Proszę wyswietli macierze U1, U2 i do nich odwrotne dla N=10.
  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 , 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)
  6. Powtórz zadanie z macierza Hilberta ale w arytmetyce pojedynczej precyzji. Musimy to sztucznie wymusic aby zmienne byly zmiennymi pojedynczej precyzji za pomoca funkcji octave'a single() .
  7. 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
  8. 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ą) .
  9. Sprawdzic czy macierz 3diagonalna z 2 na diagonali i -1 na pod i nad diagonali jest symetryczna i dodatnio okreslona z wykorzystaniem funkcji octave'a chol() .
  10. (macierz Schura - nastepny lab?) Zaprogramuj rozwiazywanie ukladu rownan liniowych z macierza B=[A ,b; c^T,a] dla A macierzy trojdiagonalnej (lub ktorej rozklad LU znamy) n X n, b,c wektorami dlugosci n i a skalarem. Funkcja ma sprawdzic czy (numerycznie) macierz B nieosobliwa.
  11. 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 .
  12. 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 .
  13. (pewnie do domu) 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.
  lu18.m - skrypt z rozwiazaniami kilku zad z lu
• Lab 3 - obliczenia w fl.
  1. 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)?
  2. 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)
  3. narysuj wykres funkcji f(x)=(x-a)+a dla a=5^p dla p=10,...,30 - koniecznie liczac zgodnei z nawiasami jak napisano!
  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. Policzyć f(x)=x-\sqrt{1+x*x} raz algorytmem wprost wynikajacym z tego wzoru a raz z wykorzystaniem równoważnego wzoru f(x)=-1/(x+\sqrt{1+x*x}) tzn. Algorytm 1: $a= \sqrt{1+x^2}; w_1=x-a.
     Algorytm 2: a=\sqrt{1+x^2}; w_2=\frac{-1}{x+a} dla x=10^k i k=4,...,10. Czy widać różnice w wyniku? Powtorzyć w arytmetyce pojedynczej precyzji tzn. z wykorzystaniem funkcji single(x).
  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. 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 dla n bliskich 20 jest coraz lepiej - np. można sprawdzić czy przynajmniej obliczone przybliżenia całek spełniają: 1/(6(n+1)) = \int_0^1 x^n/(1+5) dx <= I_n <=\int_0^1 x^n/(0+5) dx = 1/(5*(n+1))
  8. Zastosować bisekcję dla funkcji (x-2)^3 liczonej z wzoru na rozwinięcie dwumianu x^3-...-8 startując z a=2-h i b=2+2h dla h=1e-3, 1e-4,1e-5 z warunkiem stopu że błąd jest mniejszy od 1e-20 (czyli długości odcinka w którym jest rozwiązanie < 2e-20). Narysować wykres tej funkcji liczone z obu wzorów na odcinkach 2+[-h,2*h] dla h=10^{-3},...,10^{-5}. Wsk: metoda bisekcji - szukamy zera funkcji ciaglej f: startujemy z odcinka [a,b] takiego, ze f(a)*f(b)<0 - patrzymy na znak f(x) dla x=(a+b)/2, i zamieniamy x z a jesli znak f(x) taki jak f(a) albo x z b w.p.p (no chyba ze trafimy w zero tzn f(x)=0). Otrzymujemy ciag przedzialów [a_k,b_k] o dlugosci 2^{-k}|b-a| takich ze a_k i b_k zbiegaja do zera f (jakiegos jak jest ich wiecej niz jedno). Zatrzymujemy sie gy |a_k-b_k| < TOL dla zadanej tolerancji TOL. Za przyblizenie zera mozemy wziac a_k,b_k lub (a_k+b_k)/2.
  fl18.m - skrypt z rozwiazaniami zadan z fl z labu nr 2 (fl) - te które byly zadane w labie
• Lab 4 - LZNK, metoda Householdera, rozkład QR
  1. zaliczenie projektu nr 1.
  2. 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]), 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ń)?
  3. Porównaj 2 metody rozwiązania LZNK dla wektora b i macierzy A m x n, m>=n - bedacej pierwszymi n kolumnami macierzy hilberta m x m (Funkcja hilb() tworzy taką macierz ze zmienionymi kolejnościami wierszy/kolumn). Pierwsza - operator backslash z octave'a, druga tworzymy układ równań 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 - pierwsza kolumna A. (Rozwiązanie sol - wersor). Policzyć normę różnicy x_k-sol dla x_k rozwiązania policzonego k-tą metodą. Wsk: Funkcja hilb tworzy taką macierz dla m=n (tylko z odpowienio przepermutownymi kolumnami/wierszami) można stworzyć macierz mxm i wyciac odpowiednią podmacierz.
  4. 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 ?
  5. Zaprogramuj funkcję octave'a function y=H(x,w,nw) 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ę.
  6. Weź losowy wektor Householdera - stwórz macierz Hw i sprawdz czy Hw symetryczna (||Hw-Hw'||_1==0?) i Hw*Hw=Id, oraz czy dla losowych wektorów x powiedzmy wymiaru 100 - ||Hw*x||_2=||x||_2 oraz ||H||_2 - tu dla mniejszych rozmiarów.
  7. (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} .
  8. zaprogramuj metode Householdera dla regresji liniowej czy ogolniej dla dowolnej macierzy A 2 kolumnowej i wektora prawej strony y tzn znajdz w1, w2 wektory Householdera takie ze H2H1A=R macierz górnotrójkatna (H1,H2 - macierze householdera dla wektorow Hous. w1,w2 odpowiednio) - zastosuj znaleziony rozklad do rozwiania LZNK z A i y. Porównaj z rozwiazaniem otrzymanym backslashem.
  9. Zastosuj standardowej funkcję octave'a tzn \ oraz metodę korzystającą z funkcji qr(A) która zwraca rozkład QR macierzy ( Q ortogonalna) do rozwiązania zadania znalezienia prostej postaci b+ax najlepiej przybliżającej N punktów punkty (k,1+2*k+\epsilon_k) dla k=1,\ldots,N gdzie \epsilon(\epsilon_1,\ldots,\epsilon_N) losowy wektor za zakresu [-\epsilon,\epsilon] testować dla \epsilon= 1,1e-1,1e-2,1e-3] . (Funkcja rand(n) generuje wektor losowy o rozkładzie jednostajnym na [0,1] w octavie).
  10. 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) ?
  11. Zaprogramuj metode rozwiazywania LZNK z f wektorem prawej strony i A=u*v' dla danych niezerowych wektorów u,v tzn macierzy rzedu jeden z wykorzystaniem 2 macierzy Householdera H1 i H2 takich, ze H1*u=||u||_2*e_1 i H2*v=||v||_2*e_1 (e_1=(1,0,..,0)' wersor w odpowiedniej przestrzeni). Przetestuj dla losowych u,v i f - porównaj z wynikiem uzyskanym backslash. Jak mona rozwiazac zadanie z wykorzystaniem funkcji svd(A)? (tzn jak sie ma wynik svd(A) do H1*A*H2?)
  12. Sprawdx dla losowej macierzy A mxn i losowego b czy rozwiazanie ukladu blokowego (notacja octave'a) [I,A;A^T, 0][r',x']'=[0,b']' daje x bedacy rozwiazaniem LZNK z A i b, tzn stwórz macierz blokowa rozwiaz uklad i porownaj czy ostatnich n skladowych równe rozwiazaniu LZNK otrzymanego operatorem backslash.
  13. Szybko mnoznie macierzy Householdera przez wektor (macierz) policz Hx dla losowego x (macierzy A) i macierzy Householdera H z wektorem np losowym w na 3 sposoby:
      1. stwórz H i mnóz Hx (HA) - policz czas tylko mnzoenie Hx,
      2. y=Hx=I-(2/(w^T*w))*(w*w^T)*x
      3. y=Hx=I-(2/(w^T*w))*w*(w^T*x) (tylko zmienilismy kolejnosc mnozenia)
      Oczywicie mozna policzyc wczesniej a=2/(w^Tw) czy unormowac wektor w - to jakosciowo niczego nie zmieni.
  14. Sprowadzanie macierzy symetrycznej do trójdiagonalnej: wylosuj A - policz B=(A+A') znajdz taki w wektor Householdra ze HBH ma zerowe elementy od 3ciego w pierwszej kolumnie i 1-szym wierszu.
  
  lznk18.m - rozwiązania niektórych zadań z labu w tym na metode Householdera
• Lab 5 - Iteracyjne metody rozwiazywania url. Numeryczne zadanie wlasne.
  1. Zaimplementuj metode Richardsona x_{n+1}=x_n+\tau*r_n dla r_n=b-Ax_n z A=A^T>0 -ze znanymu wartościami własnymi np Q*Lam*Q^T dla Q ortogonalnej i Lam diagonalnej z dodatnią diagonalą. Weź Lam=diag(10,9,..,1) i sprawdź czy metoda zbiega dla losowych \tau\in (0,2/10) - czy szybkosc najwieksza dla \tau^*=0.5(1+10)? Tzn dla znanego rozwiązania tzn b=As liczymy ||x_n-s||_2/||x_0-s||_2 za x_0 można wziąć wektor zerowy. Powtórz zadanie dla 10*I+A. Porównaj z metodą Jakobiego x_{n+1}=x_n+D^{-1}r_n (D diagonala macierzy A) i Gaussa-Seidla x_{n+1}=x_n+(L+D)^{-1}r_n (L dolnotrójkątna część A - bez diagonali czyli D+L - dolnotrójkątna część A) - ten wzó? oczywiście nie jest optymalny z pktu widzenia implementacji.
  2. Napisz funkcje obliczajaca dominujaca w module wartość własna i jej wektor własny za pomoca metody potęgowej x(k+1)=A(x(k))/||Ax(k)||_2. Iloraz Rayleigha r(k)=x(k)'*A*x(k). Warunek stopy gdy różnica między kolejnymi ilorazami Rayleigha mniejsze od zadanej tolerancji np 1e-10. Argumenty: A macierz kwadratowa symetryczna, funkcja ma zwracać przybliżoną wartość własna (iloraz R.) i wektor własny dla niej.
  3. Przetestuj metodę potęgową dla A=[4,1;1,4] z x0=[2;1]/sqrt(3). Policz błędy dla wektora własnego i wartości własnej - czy są takiego samego rzedu?
  4. Napisz funkcje obliczająca wektor własny i wartość wlasna za pomoca odwrotnej m. potegowej - argumenty macierz A i przyblizenie wartosci wlasnej mu - przetestuj dla A symetrycznej z wartosciami wlasnymi 10,9,..,1 (wygeneruj losowa Q 10x10 i A=QDQ' dla D diagonalnej z wartosciami wlasnymi na diagonali) i mu=8+1e-k k=1,..,16 - czy ptzy ej samej zadanej tolerancji ilosc iteracji taka sama? Przetestuj dla mu=9.5 - wtedy metoda zbiega (ma 2 podciagi zbiezne do sumy i odp. róznicy wartosci wlasnych dla 8 i 9.
  5. Przetetuj metode potegowa i odwrotna potegowa dla A diagonalizowalnej ale nie w bazie ortogonalnej tzn A=CDC^{-1} D diagonalna z w.wl. na diagonali. C nieosobliwa. (Wylosuj C i ustal D) Kolumny C to wektory wlasne A. Zobacz czy szybkosc zbieznosci ilorazów Rayleigha i iteracji met. potegowej (odwrotnej potegowej) sa róznego rzedu.
  6. Przetetuj metodę potegową dla A biorąc x0 ortogonalne do v1 wektoru własnego dla dominujacej wartości własnej co do modułu. Czy metoda zbiegnie do wektora własnego drugiego co do wielkości?
  7. Wystartuj metodę potegową dla A (np A symetryczna z w.wl. 50,9,..,1) biorąc za x0 wektor własny dla jakiejś niedominującej wartości własnej np lambda_2=9 - wyświetlaj ilorazy Rayleigha dla kolejnych 50 iteracji (ilość iteracji zależy od stosunku 2 największych w module ww)? Czy stale będą równe lambda_2?
  8. Zaimplemntuj metodę QR: A_0=A; liczymy rozkłąd QR A_k tzn Q_k ortogonalna R_k g-trójkątna : A_k=Q_kR_k i liczymy A_{k+1}=R_kQ_k. Czy dla A symetrycznej iteracja zbiegnie do diagonalnej. Przetestuj dla A z ww : 10,9,8,..,1 a następnie A*A czy A*A*A. Jak znajdując wartości własne A znaleźć też wektory własne? (Wsk: jakie są macierze podobieństwa między A_k a A_0=A)?
     power18.m - funkcja z prostą implementacją metody potegowej i jej testami
     miter18.m - funkcja z prostą implementacją metody Richardsona i kilkoma prostymi testami
• Lab 6 Interpolacja wielomianowa lagrange'a, hermite'a i splajnowa
  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. 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 10000 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 ?
  5. Interpolacja hermite'a - napisz funkcje która dla danych 2 wezlów (x0,x1) i 4 wartsci (f0,df0,f1,df1) znajdzie wielomian interpolacyjny Hermite'a stopnia <=3 taki, ze w(xk)=fk i w'(xk)=dfk dla k=0,1 - narysuj wykres biorac za fk wartosci sin a dfk wartosci sin'=cos w wezlach x0=-1 i x1=3 - narysuj wykresy. (Wskazówka: najprosciej wykreowac macierz w bazie naturalnej - rozwiazac uklad backslashem).
  6. 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] .
  7. 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óry 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().
  8. 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.
  9. 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() ?)
  10. 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 ?
  
  interp18.m - rozwiązania niektórych zadań na interpolacje wielomianową Lagrange'a i splajnami kubicznymi
  
• Lab 7 - dokońćzenie interpolacji wielomianowej, zaliczenie projektu nr 2 i metody rozwiązywania równań nieliniowych
  • Zaliczenie projektu.
  • Powtórz zadania na interpolację Lagrnage'a ale dla węzłów Czebyszewa w szczególności testy dla przykłady Rungego tzn f(x)=1/(1+x^2) na [-5,5]. Czy cia? wiel. interp. Lagrange'a na węzłach Czebyszewa zbiega?
  • Zaimplementować metodę Newtona x_{n+1}=x_n -f(x_n)/df(x_n) 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,
    • (Zerow k-krotne) (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).
    • atan(x)=arctan(x)=0 proszę spróbować dla jakich x0 metoda zbiegnie biorąc kolejno x0=0.5,1,2,4,8,16,32,64,128,..,2^k (i to samo dla x0 ujemnych). jak ustalicie Państwo takie x0 że dla 2k zbiega a dla 2^(k+1) nie - można połowiąc odcinek znaleźć takie a że dla x0 < a mamy zbieżność a dla x0>a nie mamy... Złóżmy 15 Iże rozbieżność zachodzi jak przekroczymy max ilość iteracji (np 60) albo błąd rośnie (znamy rozwiązanie x^*=0)
    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. Proszę dobierać różne wartości startowe x_0 poza zaproponowanymi.
  • 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 (kolejne zadanie).
  • 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. Za x1 proszę wziąć x1=x0-f(x0)*h/(f(x0+h)-f(x0)) dla h małego np h=1e-6 czyli x1 to iteracja m. Newtona z przybliżoną pochodną ilorazem różnicowym.
  • Sprawdzić czy metoda iteracji prostych x_n=cos(x_{n-1}) zbiega do x=cos(x). Zbadac czy zbieżność jest liniowa.
  • Porównać wyniki z zad 1 z funkcjami octave'a fsolve() i fzero() (pomoc: help fsolve/fzero).
  • (wpływ fl na metodę bisekcji) zaimplementuj metodę bisekcji dla f(x)=(x-s)^3=0 ale f obliczanego ze wzoru f(x)=x^3-3*s*x^2+3*s^2*x-s^3. Proszę wziąć odcinek startowy losowy z 0 \in [a,b] (np losując lewy koniec a z s+[-2,-1] a prawy koniec b z s+[3,4]). Zatrzymać iterację gdy długość odcinka w k-tej iteracji [a_k,b_k] jest mniejsza od tol=1e-12,1e-13 etc Proszę testować czy na pewno cały czas a_k < s a b_k > s. (Jak wezmiemy s2 lub 3 to iterracja sie zatrzyma bo mzoemy dostac f(x)==) w fl -stad lepiej s losowe)
  • 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?
  • 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?)
  • Zaimplementować wielowymiarową metodę Newtona w wersji z dokładnym Jakobianem i 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.
  
  

Literatura do matlaba (octave to jego klon.... niektóre funkcje sie nazywaja inaczej itp)
• Desmond J. Higham, Nicholas J. Higham, MATLAB Guide, SIAM, 2005
• Cleve Moler, Numerical Computing with MATLAB SIAM, 2004, wersja on-line - sekcje to pliki pdf - archiwum skryptow uzywanych w ksiazce
• Kermit Sigmon, MATLAB Primer, dostepne za darmo - choc opisuje stara wersje matlaba

Zadania ze starych egzaminów z metod numerycznych:

1. I termin - 2010/11
2. I termin - 2009/10
3. II termin - 2009/10
4. Zadania z kolejnych lat (do 2013/14) w skrypcie o MN i MO prof. P. Kiciaka.
5. Zadania z MN z lat 2014/15 i 2015/16 na stronie prof. Leszka Plaskoty.

Literatura:

• [KC2006] D.Kincaid, W.Cheney, Analiza numeryczna, WNT, 2006
• [Kic2015] Przemyslaw Kiciak, Skrypt do MO i MN, plik pdf, 2015.
• [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.