Semestr letni 2012/13

• Ćwiczenia/Lab
• Program bieżącego labu
• Zadania domowe, projekty na lab, stare zadania z kolokwium
• Wyniki zadań domowych, projektów z labu, kolokwium etc: plik-pdf

Wykład

Program wykładu

1. Wstęp - czym jest matematyka obliczeniowa (0.5 wykładu).
2. Metody rozwiązywania równania nieliniowego (skalarnego) - bisekcja, metoda newtona i iteracji prostej(1-2 wykłady)
3. Arytmetyka fl (1-2 wykłady)
4. Numeryczna algebra liniowa - rozwiązywanie układów równań liniowych(rozkłady LU, Choleskiego i QR, współczynnik uwarunkowanie macierzy) LZNK (liniowe zadanie najmniejszych kwadratów) - metoda Householdera, numeryczne zadanie własne (metoda potegowa, odwrotna i QR) (5-6 wykładów)
5. Interpolacja wielomianowa Lagrange'a i Hermite'a, interpolacja splajnowa - splajny liniowe i kubiczne (3-4 wykłady)
6. Numeryczne całkowanie: rzad kwadratury, kwadratura interpolacyjna z wzorem na bład, kwadratury złozone, kwadratury Gaussa(2 wykłady)

Zasady zaliczenia przedmiotu

Ćwiczenia/Lab MO

• Gr 1 ćwiczenia piątek 830-10, sala 5070. Pierwszy lab 22 luty i potem co 2 tygodnie
• Gr 2 ćwiczenia - piątek 1215-1345, sala 5070. Pierwszy lab 1 marca i potem co 2 tygodnie

Zaliczenie ćwiczeń i labu

1. Pierwsza seria ZD równania nieliniowe
2. Pierwszy projekt z labu (testy m/ Steffensena) na 1 lab po 21 kwietnia 2013 (26 kwietnia i 10 maja 2013)
3. Druga seria ZD numeryczna algebra liniowa (uklady r. liniowych, eliminacja Gaussa czyli rozklad LU, normy, uwarunkowanie, rozkłady QR, LZNK, zadanie własne)
4. Drugi projekt z labu (prosta najbliższa danym punktom jako LZNK i QR za pomocą ortogonalizacji Gramma-Schmidta) na ostatni lab
5. Zadania z kolokwiów (2gr i poprawkowe) i egzaminu w I terminie z roku 2011/12
6. Zadania z kolokwium(2 grupy) z bieżącego roku tzn. 2012/13
7. Trzecia seria ZD - interpolacja Lagrange'a, interpolacja Hermite'a, splajny liniowe i kubiczne
8. Czwarta seria ZD - kwadratury interpolacyjne, kwadratury Gaussa
Projekt z labu trzeba pokazać w akcji w czasie zajęć labu, za zadania oddane po terminie otrzymują Państwo 50% punktów.

Program labów

• Lab 1 (gr 1 - 22/02/13; gr 2 - 01/03/13) - 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). Rysowanie wykresów funkcji - czyli funkcja plot(). Wskaźniki do funkcji (handle - operator @).
  mobasic.m - skrypt z przykladami powyzszych operacji czyli tym o czym jest ten lab.
  Kilka zadan (aby nabrac wprawy):
  1. Operacje macierzowe
     • Utwórz dowolne macierze 3x4 A i 3x5 B - a następnie macierz 3x9 C, której pierwsze 4 kolumny to A a kolejne to B.
     • Z macierzy C 'wytnij' podmacierz D składającą się z 1 głównego minora 3x3 tzn. od C(1,1) do C(3,3).
     • Zamień kolejność kolumn D.
     • Zamień kolejność wierszy D.
     • Wytnij dolnotrójkątną część macierzy D a potem górnotrójkątną
     • Wstaw D z powrotem do C jako główny minor.
     • Policz \sin(D)=(\sin(D_{ij}) od D.
     • Zapisz D do pliku (binarnego i ASCII) - zamień element D(1,1) na -100 i wczytaj nową macierz do octave'a.
  2. Policz dyskretną normę max od (\sin(x))^2 na [0,1] (wektorowo- czyli bez użycia pętli).
  3. Narysuj wykres sin(x) na odcinku [0,4].
  4. Znajdź przybliżone max i min funkcji x*(3+2*cos(x)) na [-1,5] oraz przybliżenia punktów ekstremalnych. Zrób to samo dla jakiegoś wielomianu stopnia dwa i trzy np x^3+x^2-x-4
  5. Utwórz funkcję w m-pliku np. obliczającą (\sin(x))^2
  6. Utwórz m-plik obliczający wartość funkcji z parametrem np \sin(a*x) - parametr przekaż jako zmienną globalną
  7. Funkcje anonimowe - utwórz funkcję (\sin(x))^2 jako funkcję anonimową - czyli w linii komend
  8. Przy pomocy pętli for i while oblicz sumę 1+...+N dla N=100. Użyj if aby sprawdzić czy równa się tyle ile powinno 0.5*N*(N+1)- wyprowadź na ekran komunikat używając printf().
• Lab 2 (gr 1 - 8/03/13; gr 2 - 15/03/13) - równania nieliniowe
  1. Zaimplementuj metodę bisekcji w skrypcie - dla rozwiązania równania cos(x)=0 na odcinku [0,2] - sprawdź jak działa.
  2. 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 przypadkowo wybrane
  3. Napisz funkcję octave'a: function [y,iter,osz,kod]=bisekcja(nazwa_funkcji,a,b,tola,maxit) w pliku bisekcja.m - której parametrami będą nazwa funkcji (czy funkcja - jak to zrobić odsyłam do pomocy czy stron www), końce przedziału a i b , tola - żądana dokładność bezwzględna (błąd ma być mniejszy od epsilona), maxit - maksymalna ilość iteracji. Funkcja ma zwrócić y - przybliżone rozwiązanie, iter - ilość iteracji, osz - oszacowanie błędu , i kod-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łać również jak podamy tylko trzy czy cztery pierwsze argumenty - i wtedy maksymalna ilość iteracji domyślnie ma być 100, a tolerancja 1e-5 (jak tylko trzy argumenty).
  4. Znajdź w pomocy funkcję octave'a - rozwiązująca równania nieliniowe - i przetestuj ją na kilku prostych przykładach skalarnych np.
     [ f1(x)=]
     1. cos(x)=0 ,
     2. x^2-2=0 ,
     3. x^{10}-2=0 ,
     4. (x-1)*(x-2)*(x-3)=0 ,
     5. x-\sin(x)=0 .
     A dokładniej zastosuj tę funkcję do rozwiązywania powyższych równań dla różnych przybliżeń startowych.
  5. Zaimplementować metodę Newtona w octavie i przetestować jej zbieżność dla następujących funkcji:
     [ f1(x)=]
     • x^2-2 z x_0=2 ,
     • x^3-27 z x_0=27 ,
     • \exp(x)-2 z x_0=10,-10,1e3 ,
     • \sin(x) dla x_0=2 ,
     • cos(x) z x_0=2
     • (x-2)^k x x_0=3 dla k=2,4,8,16 ,
     • x*x-2 z x_0=1e6 , (czy zbiega w ogóle a jak tak to czy kwadratowo)
     • 1/x-a dla danych a=0.5,2,4,100 (oczywiście implementując bez dzielenia) jak dobrać x_0 ?
     Dla wszystkich tych funkcji znamy rozwiązania więc można wyświetlać na ekranie błąd 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 .
  6. Powtórzyć poprzednie zadanie ale zastępując metodę Newtona przez metodę 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 .
  7. Sprawdzić czy metoda iteracji prostych x_n=cos(x_{n-1}) zbiega do x=cos(x) . Zbadać eksperymentalnie czy zbieżność jest liniowa, tzn. czy |x-x_n|/|x-x_{n-1}| zbiega do stałej.
  8. Porównać wyniki z zadania poprzedniego z funkcją octave'a fsolve() (pomoc: help fsolve ).
  9. Zaimplementować przybliżoną metodę 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. h=1e-4,1e-7,1e-10 itp Porównać zbieżność z dokładną metodą Newtona (szczególnie ostatnie iteracje) dla funkcji z poprzednich zadań.
  10. \textsl{Rozwikływanie funkcji:} dla funkcji y(x) zadanej równaniem g(x,y)=2x^2+3y^2-3=0 znaleźć wartości y_k=y(x_k) dla x_k=k*h dla k=-N,...,N i h=1/N (N - 10,20,40,...) Znajdować y_k rozwiązując układ równań g(x_k,y_k)=0 . Jak w kolejnych krokach dobierać przybliżenie startowe w metodzie Newtona?
  11. Odwracanie funkcji: rozpatrzmy daną funkcje np. f(x)=\sin(x)+2*x znaleźć wartości funkcji odwrotnej na odcinku [0,5] na siatce k*h dla k=0,..,100 . Narysować wykres tej funkcji. (Sam wykres można narysować dużo prościej bez wyliczania wartości funkcji odwrotnej. Jak?)
  12. Zaimplementować wielowymiarową metodę Newtona w wersji z dokładnym Jakobianem i w wersji gdy jakobian przybliżany różnicami dzielonymi z parametrem h. 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. Taką funkcje w octavie definiujemy:
      function y=F(x)
      #x- wektor pionowy
      y=[[1 2]*x-1;3.*x(1)*x(1)+x(2)*x(2)-1];
      endfunction
• Lab 3 (gr 1 - 22/03/13; gr 2 - 12/04/13) - dokończenie porzedniego labu i arytmetyka fl
  Celem labu jest wykonanie odpowiednich zadań, które wykażą pewne cechy arytmetyki fl. Domyślnie w octavie wykorzystywane są liczby zmiennopozycyjne o podwójnej precyzji, jakkolwiek w najnowszych wersjach octave'a można również trochę sztucznie wymusić używanie zmiennych w precyzji pojedynczej przy pomocy funkcji single(a)! zwracającej zmienną pojedynczej precyzji z tą samą wartością. Z konstrukcji arytmetyki fl wynika, że odejmowanie dwóch wartości o tym samym znaku o małej różnicy może skutkować dużą utratą dokładności.
  1. Funkcje single(a) i eps. Wywołaj pomoc do tych funkcji w octavie. Sprawdź czy 1+eps obliczone w arytmetyce podwójnej precyzji w octavie jest większe od jeden. Przyjmij, że a=single(eps) i sprawdź czy ponownie 1+a jest większe od 1 .
  2. Wyznacz 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)? A co dla liczb w pojedynczej precyzji? Funkcja octave'a single(x) tworzy takie zmienne. Wykorzystując tę funkcję ponownie wyznacz epsilon maszynowy ale w pojedynczej precyzji.
  3. Narysuj wykres funkcji f(x)=(1+a)-a na [0,1] dla różnych a=1ek dla k=1,2,..,20 . Tutaj ważne aby obliczać wartość f(x) z tego wzoru
  4. Policzyć f(x)=x-\sqrt{1+x*x} raz algorytmem wprost wynikającym z tego wzoru a raz z wykorzystaniem równoważnego wzoru f(x)=\frac{-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óżnicę w wyniku? Powtórzyć w arytmetyce pojedynczej precyzji tzn. z wykorzystaniem funkcji single(x).
  5. Obliczyć wartości wielomianu f_1(x)=(x-2)^4=f_2(x)=x^4-.....+16 dwoma algorytmami:
     Algorytm 1
     a= (x-2), f_1(x)=a^4;
     Algorytm 2
     f_2(x)=x^4-...+16
     na siatce równomiernej 1000 punktowej na [2-a,2+a] dla a=1e-3 . Narysować wykresy obu funkcji i policzyć błąd || f_1(x)-f_2(x) ||_\infty czyli \max_k |f_1(x(k)-f_2(x(k))| . Tu x(k)=2-a+k*h dla h=a/500 . Wektor x można utworzyć w octavie przy pomocy funkcji linspace().
  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 n=0,..,20 dwoma algorytmami ze wzoru I_n+5*I_{n-1}=1/n . Pierwszy algorytm przyjmuje I_0= log(6/5) i oblicza z powyższego wzoru kolejne I_n dla n=1,2,3,... . Drugi algorytm wykorzystuje fakt, że \frac{1}{(n+1)6} = I_n= \frac{1}{(n+1)5} zatem w tym algorytmie przyjmujemy za I_{30} jakąkolwiek wartość z tego przedziału np. I_{30}=1/180 i iterujemy w tył, tzn. dla n=30,29,...,20,...,0 . Porównaj wyniki obu algorytmów dla n=0,...,30 oraz czy wyniki spełniają oszacowanie. Dlaczego jeden z algorytmów działa zdecydowanie lepiej w arytmetyce fl? Jako dodatkowe zadanie pozostawimy uzasadnienie wzoru rekurencyjnego i oszacowania I_n wykorzystywanych w algorytmach.
  8. Zastosować algorytm bisekcji 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 . Jako warunek zakończenia działania algorytmu przyjmujemy, że błąd jest mniejszy od 1e-20 (czy długości odcinka w którym jest rozwiązanie). Czy algorytm zwraca przybliżenie liczby 2 jedynego pierwiastka tego wielomianu? Narysować wykres tej funkcji liczone z obu wzorów. tzn. (x-2)^3 i x^3-...-8 na odcinkach 2+[-h,2*h] dla h=10^{-3},...,10^{-5} . Czy z obu wykresów wynika, że ten wielomian ma tylko jedno miejsce zerowe w otoczeniu 2 ?
• Lab 4 (gr 1 - 5/04/12; gr 2 - 26/04/12) - rozwiązywanie układów liniowych: rozkłady LU, LL', uwarunkowanie macierzy, LZNK
  1. Normy macierzowe - policz normy 1,2,Frobeniusa, max dla macierzy diagonalnej i konkretnej niesymetrycznej A=[1;2;1;1].
  2. Operator octave'a \ służący m.in. do rozwiązywaniu układów równań liniowych w octavie. Przetestuj ten operator dla macierzy A=[1,2;3,4] i x=[1;3] z f=A*x czy y=A\f jest rozwiązaniem tego układu. Policz normy rezydualną ||Ay-f||_p i normę błędu ||y-x||_p dla p=1,2 . Powtórz testy dla jakiejś macierzy osobliwej? Co się wtedy dzieje?
  3. Funkcja inv(). Zapoznaj się z pomocą do funkcji: inv() . Przetestuj dla A=[1,1;1,-1] czy macierz B obliczona za pomocą tej funkcji rzeczywiście jest macierzą odwrotną? Policz normy pierwszą i Frobeniusa ||A*B-I|| oraz B*A-I|| . Zastosuj funkcje do rozwiązywania układu równań liniowych Ax=f dla znanego x (liczymy f=A*x ). Policz y jako iloczyn B z f i policz błędy 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 losową) dla n=10,50,250,1250 ze znanym rozwiązaniem - oszacuj czas przy pomocy funkcji tic i toc . Porównaj czas i błędy w normie pierwszej dla rozwiązania uzyskanego przy pomocy operator octave'a \ .
  4. Funkcje lu() i chol() . Zapoznaj się z pomocą do funkcji: lu() i chol() . Dla macierzy A=[1,1;1,-1] znajdź jej rozkład PA=LU , rozkład Choleskiego A=L_1^TL_1 . Sprawdź te rozkłady licząc normy macierzowe pierwszą i Frobeniusa błędów np. PA-LU czy A-L_1^TL_1 . Zastosuj te rozkłady do znalezienia rozwiązania układu równań liniowych Ax=f ze znanym rozwiązaniem np. x=[1;1] i f=[2;0] . Policz normy pierwszą i drugą wektorowe pomiędzy x a wynikiem algorytmu w polegającym na zastosowaniu odpowiedniego rozkładu oraz takie same normy rezydualne tzn. normy Aw-f .
  5. Przetestuj solver octave'a dla macierzy górnotr. U1 i U2=-U1+2*I dla U1 z 1 na diagonali i 1 na wszystkich pozycjach nad diagonala (polecenie octave'a triu(ones(n,n)) ją tworzy), tzn. stwórz macierze Uk (k=1,2) dla N=10,20,30,50,100 dla znanego rozwiązania np. losowego czyli sol_k=1 policz Uk*sol=fk - rozwiąż w octavie układ z Uk i fk i policz normy (2 itd) ||Uk*xk-fk|| i ||xk- sol|| dla xk rozwiązania które wyliczył octave (k=1,2).
  6. 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. Policz uwarunkowanie macierzy Hilberta dla powyższych N .
  7. 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. Policz uwarunkowanie macierzy Hilberta dla powyższych N .
  8. Powtórz zadanie z macierzą Hilberta ale w arytmetyce pojedynczej precyzji. Musimy tu trochę sztucznie wymusić aby zmienne były zmiennymi pojedynczej precyzji za pomocą funkcji octave'a single() .
  9. W octavie przetestuj eliminacje Gaussa z częściowym wyborem i bez wyboru dla macierzy A=[e, 1; 1, 1] z e=eps/10 (funkcja octave'a eps zwraca epsilon maszynowy) i wektorem prawej strony f=[1;0]^T . Trzeba zaprogramować eliminację Gaussa bez wyboru dla macierzy 2\times 2 samemu.
  10. Przetestuj funkcję invhilb() tworzącą macierz odwrotną do macierzy Hilberta (por. zadanie z macierzą Hilberta). Policz normy macierzowe Frobeniusa i pierwsze dla H(N)*iH(N)-I dla dla N=10,20,30 , gdzie iH(N) macierz odwrotna do macierzy Hilberta obliczona przez invhilb() . Przetestuj wykorzystanie iH(N) do rozwiązania układu równań z macierzą Hilberta H(N) , tzn. stwórz macierz H(N) dla N=10,20,30 dla znanego rozwiązania sol , czyli policz H(N)*sol=f - rozwiąż układ H(N)x=f mnożąc F przez iH(N) tzn. x=iH(N)f i policz normy (np. 1,2 ) ||H(N)x-f|| i ||x- sol|| .
  11. Stworzyć w octavie macierz trójdiagonalną A=(a_{i j})_{i,j=1}^n wymiaru n\times n z 2 na diagonali i -1 na pod- i nad diagonalą tzn.
      a_{i, j}= 2 dla i=j
      -1 dla |i-j|=1
      0 dla|i-j|>1
      dla i,j=1,...,n. 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?).
  12. Sprawdzić czy macierz 3doagonalna z 2 na diagonali i -1 na pod i nad diagonali jest symetryczna i dodatnio określona z wykorzystaniem funkcji octave'a chol() .
  13. Zaimplementować rozkład LU dla macierzy trójdiagonalnej bez wyboru elementu głównego, tzn. stworzyć własną funkcję octave'a
      function [x,d,l]=lu3diag(a,b,c,f,N)
      Parametry funkcji:
      1. a,b,c - przekątna, pod-przekątna i nad-przekątna macierzy A ,
      2. f - wektor prawej strony,
      3. N długość przekątnej a .
      Funkcja zwraca x - rozwiązanie Ax=f oraz czynniki rozkładu macierzy A=LU czyli diagonalę macierzy górnotrójkątnej U w wektorze d oraz pod-diagonalę L - macierzy dolnotrójkątnej, trójdiagonalnej z jedynkami na diagonali czyli wektor l . Pamiętajmy, że nad-diagonala U równa się na-diagonali A tzn. wektorowi c .
  14. Zastosować funkcję z poprzedniego zadania do uzyskania rozkładu LU macierzy trójdiagonalnej dla N=4,10,15 . Porównać z wynikiem uzyskanym przy pomocy funkcji octave'a lu().
  15. Zastosować funkcję z poprzedniego zadania (2 do tyłu) do rozwiązania układu równań z macierzą z z 2 na diagonali i -1 na pod i nad diagonalach Porównać czas rozwiązywania Państwa algorytmem, a algorytmem octave'a czyli operatorem A\ f!. Rozwiązać układ dla dowolnego rozwiązania np. 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ą).
• Lab 5 (gr 1 - 19/04/13; gr 2 - 25 maja 2013) -rozkład QR, regularne LZNK, przekształcenie Householdera
  1. Operator octave'a \ służący m.in. do rozwiązywaniu układów równań liniowych w octavie al również LZNK.
     • Przetestuj ten operator dla RLZNK dla macierzy A=[1,1;1,-1;1,3] i x=[1;2] z f=A*x czy y=A\f jest rozwiązaniem tego układu? Policz normy rezydualną ||Ay-f||_2 i normę błędu ||y-x||_2 .
     • Przetestuj ten operator dla nieregularnego LZNK dla macierzy A^T z A=[1,1;1,-1;1,3] i x=[1;1;1] z f=A*x czy y=A\f jest rozwiązaniem tego układu? Policz normy rezydualną ||Ay-f||_2 i normę błędu ||y-x||_2 .
  2. Rozkład QR w octave. Funkcje qr().
     1. Zapoznaj się z pomocą do funkcji: qr().
     2. Dla macierzy A=[1,1;1,-1;1,3] znajdź jej rozkład A=QR z pomocą funkcji qr().
     3. Sprawdź czy czy uzyskana Q jest ortogonalna - policz normy Frobeniusa QQ^T-I i Q^TQ-I .
     4. Sprawdź ten rozkład licząc normy macierzowe drugą i Frobeniusa błąd A-QR . Zastosuj ten rozkład do znalezienia rozwiązania LZNK Ax=f ze znanym rozwiązaniem np. x=[1;0] i f=[1;1;1] .
     5. Policz normę drugą wektorową pomiędzy x a wynikiem algorytmu w polegającym na zastosowaniu odpowiedniego rozkładu oraz takie same normy rezydualne tzn. normy drugie Aw-f oraz Rw-Q^Tf .
  3. Dla macierzy z poprzednich zadan rozwiaz LZNK - rozwiazujac uklad r. normalnych - porównaj wyniki.
  4. Układ równań normalnych a rozkład QR czy operator \ . Rozpatrzmy A_{2n,k} pod-macierz 2n\times k macierzy A_{2n,2n} Vandermonde'a dla 2n węzłów równoodległych na [0,1] . Rozwiąż LZNK z A_{2n,k} z wektorem prawej strony równym pierwszej kolumnie tej macierzy (rozwiązanie to pierwszy wersor) - trzema sposobami:
     1. operatorem \ ,
     2. rozkładem QR uzyskanym funkcją qr(),
     3. rozwiązaniem układu równań normalnych: tworzymy macierz układu równań normalnych B=A_{2n,k}^TA_{2n,k} , wektor prawej strony g=A_{2n,k}^T f układu równań normalnych. Rozwiązujemy układ równań normalnych Bx=g operatorem \ .
     Przetestuj dla n=10,20,40,80 i k=2,4,n . Porównaj
     • czas obliczeń
     • błąd - ||x-y||_2
     • błąd rezydualny ||Ax-f||_2
     dla x rozwiązania dokładnego LZNK, f wektora prawej strony LZNK, y przybliżenia rozwiązania uzyskanego daną metodą.
  5. Rozkład QR a operator \ przy rozwiązywaniu układów równań liniowych czyli porównanie QR a LU zastosowanych do rozwiązywania układów równań liniowych Rozpatrzmy macierz A_{n,n} Vandermonde'a dla n węzłów równoodległych na [0,1] . Rozwiąż układ równań liniowych z wektorem prawej strony równym pierwszej kolumnie tej macierzy (rozwiązanie to pierwszy wersor) - dwoma sposobami:
     1. operatorem \,
     2. rozkładem QR uzyskanym funkcją qr(),
     Przetestuj dla N=10,20,40,80 . Porównaj
     • czas obliczeń
     • błąd - ||x-y||_2
     • błąd rezydualny ||Ax-f||_2
     dla x rozwiązania dokładnego LZNK, f wektora prawej strony LZNK, y przybliżenia rozwiązania uzyskanego daną metodą.
  6. 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 ?
  7. Zaprogramuj funkcje 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. Przetestować dla losowych wektorów \vec{x} i \vec{w}, sprawdzić czy ||H_w \vec{x} ||_2=||\vec{x} ||_2. Można napisać ogólniejsza 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\not=0 wykorzystując nasza funkcje.
  8. (testujemy twierdzenie z wykładu) Znajdź wektor Householdera \vec{w} taki ze 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 tak jest. Tzn czy H_w \vec{u}= l \vec{v}.
• Lab 6 (gr 1 - 17/05/10; gr 2 - NIE (1 lab mniej) - -met. Householdera dokończenie z poprzedniego labu (o ile tego nie było) i zadanie własne - metoda potęgowa i odwrotna potęgowa, QR
  1. Przetestuj funkcje eig() dla prostej macierzy symetrycznej np [1,3;3,1] i losowej symetrycznej np losujemy A i bierzemy B=A^TA lub AA^T czy A+A^T. Znajdź ich wartości i wektory własne i sprawdź czy rzeczywiście takimi są tzn policz ||A*Q-D*Q||_1 i ||A-Q*D*Q^T||_1 - Q macierz z kolumnami - wektorami własnymi - D macierz diagonalna z wartościami własnymi na diagonali
  2. zaimplementuj metodę potęgową - jako koniec iteracji proszę przyjąć że ilorazy Raileigha spełniają |r_k-r_{k+1}|<10^{-15}.
  3. Przetestuj metodę potęgową dla A=[1,3;3,1] na ekranie drukuj iloraz Railegha i x_{2k}. Sprawdź czy obliczone wartości są rzeczywiście przybliżeniami pary własnej tzn oblicz ||A x_{2k}-r_{2k}x_{2k}||_1. porównaj z wynikiem eig().
  4. Powtórz poprzednie zadanie dla losowej macierzy symetrycznej oraz oraz macierzy o znanych wartościach własnych np 1,2,3,4,..,11 czy 10,20,30,..,110, lub 1,2,..,10,10+b dla b=0,1,10,100.
  5. Zaimplementuj odwrotną iteracje potęgową tzn. metodę potęgową dla (A-a I)^{-1}. przetestuj ją dla A z poprzedniego zadania z różnymi wartościami a np. a=0 czy a dobre przybliżenie wartości własnej (znanej).
  6. zaimplementuj metodę QR (warunek stopu norma macierzy A-D mniejsza od 10^{-12} dla D diagonali A. Przetestuj dla macierzy z zad 1 oraz macierzy o znanych wartościach własnych np 1,2,3,4,..,11 czy 10,20,30,..,110, lub 1,2,..,10,10+b dla b=0,1,10,100.
  7. zaimplementuj metodę równoczesnych iteracji tzn dla danej Z_k ortogonalnej (startujmy z Z_0=I) liczymy Y_{k+1}=AZ_k i znajdujemy rozkład QR macierzy Y_{k+1}=Z_{k+1}R_{k+1} tzn Z_{k+1} ortogonalna R_{k+1} górnotrójkątna (warunek stopu norma macierzy Z_k^TAZ_k - diagonala(Z_k^TAZ_k) mniejsza od 10^{-12} dla D diagonali Z_k^TAZ_k.) Przetestuj dla macierzy z zad 1 oraz macierzy o znanych wartościach własnych np 1,2,3,4,..,11 czy 10,20,30,..,110, lub 1,2,..,10,10+b dla b=0,1,10,100.
• Lab 7 (gr 1 - 31/05/12; gr 2 - 07/06/10) - ocena projektu nr 2; interpolacja Lagrange'a i splajnami kubicznymi
  Funkcja polyval() czyli funkcja obliczająca wartość wielomianu w jednym ub równocześnie wielu punktach czyli wektorowo. Funkcja polyfit() czyli funkcja obliczająca współczynniki wielomianu interpolacyjnego dla zadanych wartości i węzłów. Funkcje ppval() i spline() działajace analogicznie ale dla splajnów kubicznych.
  1. (tylko gr 2) Wykonaj 2-3 pierwsze zadania z poprzedniego labu (eig() i metoda potęgowa i ewent. odwrotna potęgowa)
  2. Korzystając z funkcji polyval() narysuj wykres wielomianu x^3+x-2 bez wykorzystania pętli czyli wektorowo.
  3. Test funkcji polyfit(). Wykorzystując funkcję polyfit() znajdź wielomian interpolacyjny L_n F dla funkcji F(x)=\sin(x) dla węzłów -1,0,1 oraz -1,0,1,10. Policz wartości różnicy F-L_n F w węzłach. Oraz narysuj wykresy F i L_n F na jednym rysunku korzystając funkcji plot(). Wartości wielomianu oblicz bez użycia pętli z wykorzystaniem funkcji polyval().
  4. (zadanie o zbieżności interpolacji Lagrange'a) Interpolacja Lagrange'a - zbieżność ciągu wielomianów interpolacyjnych dla węzłów równoodległych i węzłów Czebyszewa:
     • Wykorzystując funkcję polyfit() znajdź wielomiany interpolacyjne L_N f dla funkcji f=\sin() dla N węzłów równoodległych na [0,2*pi] dla N=4,8,16,32,64.
     • Oblicz dyskretną normę maksimum różnicy f-L_N f na siatce tysiąc punktowej na tym odcinku tzn. e_N=\max|f(x_k)-L_N f(x_k)|, gdzie x_k pkt siatki. I policz stosunek e_N/e_{2N} dla N=4,..,32. Czy błędy maleją do zera? jak zachowuje się e_N/e_{2N}?
     • Narysuj na ekranie wykresy \sin(x) i tych wielomianów dla różnych N - używając funkcji polyval() i plot().
  5. Powtórz (zadanie o zbieżności interpolacji Lagrange'a) dla tej samej funkcji i tego samego odcinka ale dla węzłów Czebyszewa. Węzły Czebyszewa to zera T_{n+1}(t)=\cos((n+1)\mathrm{arccos}(t)) na [-1,1] odpowiednio przesunięte i przeskalowane. Czy błędy e_N dla węzłów Czebyszewa są mniejsze niż dla węzłów równoodległych?
  6. Napisz funkcję znajdującą współczynniki w bazie potęgowej wielomianu interpolacyjnego zadanego stopnia dla węzłów równoodległych oraz węzłów Czebyszewa dla danej funkcji, odcinka [a,b]. Tzn. bardziej precyzyjnie napisz funkcję octave'a (w m-pliku): function [LN,eN]=Lagrangeinterp(FCN,a,b,N,type=0)!, która dla zadanego wskaźnika funkcyjnego FCN (ang. function handle) do funkcji jednego argumentu: function y=f(x)!, a,b - końców odcinka [a,b], N - stopnia wielomianu interpolacyjnego i typu węzłów 0 - równoodległych, 1 - Czebyszewa. Zwróci wektor LN współczynniki L_N f wielomianu interpolującego funkcję w tych węzłach oraz eN przybliżenie dyskretnej normy maksimum różnicy L_Nf - f na tym odcinku - przybliżenie normy liczymy na dyskretnej siatce zawierającej tysiąc punktów.
  7. Powtórz (zadanie o zbieżności interpolacji Lagrange'a) dla obu typów węzłów, tzn. znajdowanie wielomianów interpolacyjnych na węzłach równoodległych i węzłach Czebyszewa, ale dla funkcji f(x)=log(1+x) na odcinkach [0,1] i [0,10]. Czy dla tej funkcji i obu odcinków błędy w normach maksimum maleją wraz ze wzrostem N? Porównaj wyniki otrzymane w tym zadaniu w obliczeniach z oszacowaniami teoretycznymi błędu interpolacji Lagrange'a.
  8. Interpolacja Lagrange'a przykład Rungego. Powtórz (zadanie o zbieżności interpolacji Lagrange'a) dla obu typów węzłów, tzn. znajdowanie wielomianów interpolacyjnych na węzłach równoodległych i węzłach Czebyszewa, ale dla funkcji dla f(x)=1/(1+x*x) na [-5,5]. Czy ciąg wielomianów interpolacyjnych zbiega do f jednostajnie?
  9. 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] .
  10. 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 zadanei ale dla splajnu typu not-a-knot tzn. nie podając wartości pochodnych w końcach dla funkcji spline().
  11. 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 w {-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.
  12. Testowanie eksperymentalne rzędu zbieżności splajnu interpolacyjnego kubicznego z hermitowskimi warunkami brzegowymi. Korzystając z funkcji octave'a spline() znajdź współczynniki splajnu interpolacyjnego kubicznego S_N na N węzłach równoodległych dla funkcji f(x)=sin(x) na odcinku [-pi,2*pi] dla N=2^k*N_0 dla N_0=5 i k=1,2,3,4,5 . Następnie
      • narysuj wykresy f(x) i tych splajnów dla różnych N .
      • oblicz dyskretną normę maksimum na siatce równomiernej 1000 punktowej na tym odcinku tzn. e_N=max|sin(x_k)-S_N sin(x_k)| dla x_k punktów siatki.
      • Policz równocześnie współczynnik \frac{e_{N}}{e_{2*N}} . Czy widać, że \frac{e_{N}}{e_{2*N}}\approx 2^p dla jakiegoś p całkowitego np. p=4, 8 lub 16 ?
  13. Testowanie eksperymentalne rzędu zbieżności splajnu interpolacyjnego kubicznego bez warunków brzegowych (splajn typu not-a-knot). Powtórz poprzednie zadanie ale dla splajnów interpolacyjnych otrzymanych przez spline() bez podawania wartości pochodnych w skrajnych węzłach. Czy współczynniki {e_{N}}\{e_{2*N}} są te same? Tzn. czy szybkość zbieżności ||sin(x)-S_N||_\infty jest taka sama?
  14. Testowanie eksperymentalne rzędu zbieżności splajnu interpolacyjnego kubicznego naturalnego (warunek brzegowy - zerowanie się drugich pochodnych w końcach odcinku). Powtórz zadanie poprzednie ale dla splajnów interpolacyjnych naturalnych. Tu trzeba wykorzystać funkcję z octave-forge (czyli rozszerzenia pakietu octave)\\ pp=csape(x,y,'variational') \\- ostatni argument określa to, że splajn będzie naturalny. Tu link do pomocy do funkcji octave'a csape()
  15. Przykład Rungego czyli f(x)=1/(1+x*x) i odcinek [-5,5] a zbieżność interpolacji splajnami kubicznymi. Przetestuj jak w poprzednich zadaniach czy splajny interpolacyjne kubiczne z podanymi warunkami na pochodne w końcach odcinka zbiegają w normie supremum do f , tzn. korzystając z funkcji octave'a spline()! znajdź współczynniki splajnu interpolacyjnego kubicznego S_N na N węzłach równoodległych dla f na odcinku [-5,5] dla N=2^kN_0 dla N_0=5 i k=1,2,3,4,5 oraz narysuj wykresy f i tych splajnów dla różnych N . Następnie oblicz dyskretną normę max na siatce złożonej z tysiąca punktów na tym odcinku tzn. e_N=\max|\sin(x_k)-S_N \sin(x_k)| dla x_k=-5+k*0.01 z k=0,...,1000. Policz równocześnie współczynnik \frac{e_{N}}{e_{2*N}} . Czy \frac{e_{N}}{e_{2*N}}\approx 2^p dla jakiegoś p całkowitego?

Kilka przykładowych skryptów, m-plików (plików funkcyjnych ) octave'a, czy użytecznych linków
(dodawanych w miarę postępu labu)

• Link do stron Octave'a (skąd można ściągnąć kolejną dystrybucje - pod linuxa czy windows)
• octave-forge - rozszerzenia octave'a
  
• Manual do octave'a w htmlu (po angielsku)
• Link do skryptu do labu z mo z octave'em plik pdf

Skrypty m-pliki octave'a

• mobasic.m - prosty skrypt octave'a w którym są przedstawione niektóre podstawowe operacje macierzowe, oraz pętle, instrukcje warunkowe, zapisywanie/wczytywanie do/z plików itp
• bisekcja.m - m-plik z funkcją z implementacją metody bisekcji
• testnewton.m - m-plik z funkcja z implementacją m. newtona i testami rzędu zbieżności (do testów trzeba znać rozwiązanie równania)
• testsieczne.m - m-plik z funkcja z implementacją m. siecznych i testami rzędu zbieżności (trzeba znać rozwiązanie równania)
• fl.m - skrypt z rozwiązaniami niektórych zadań z labu o fl
• lablu.m - skrypt z rozwiązaniami niektórych zadań z labu o LU i uwarunkowaniu
• lznk.m - skrypt z rozwiązaniami niektórych zadań z labu o QR i LZNK
• simit.m - skrypt z metodą równoczesnych iteracji dla zadania własnego (symetrycznego)

Literatura:

• [KC2006] D.Kincaid, W.Cheney, Analiza numeryczna, WNT, 2006
• [Kic2012] Przemysław Kiciak, Metody numeryczne dla informatyków plik pdf, skrypt, 2012
• [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) plik pdf.

Ostatnia aktualizacja: 11 lipca 2013