Semestr letni 2009/10

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

Ćwiczenia/Lab MO

Zaliczenie ćwiczeń

Serie zadań domowych i zadań punktowanych z labu

1. Pierwsza seria ZD - termin 12 marca 2010 Sprawdziłem dokładnie tylko zad 2 i 3 (za zad 1 było 4 punkty za samo oddanie)
2. Druga seria ZD - termin 9 kwietnia 2010 sprawdzone wszystkie zadania
3. Pierwsze zadanie z labu - termin 23 kwietnia gr 3 lub 30 kwietnia gr 2. Trzeba pokazać jak działa w czasie zajęć labu.
4. Trzecia seria ZD - termin 30 kwietnia 2010 sprawdzone zad 1 i 2 (a niesprawdzone zad 3)
5. Czwarta seria ZD - termin 21 maja 2010
6. Drugie zadanie z labu - termin 30 maja gr 3 lub 2 czerwca gr 2 (możliwość zaliczenia: pokój 5010 - 12-13).
7. Piąta seria ZD - termin poniedziałek 31 maja 2010 godzina 16.00 - seria dodatkowa dla chętnych; nie gwarantuję, że zadania oddane po terminie będą sprawdzone
PROSZĘ NIE przysyłać rozwiązań mailem

Program labów

• Lab 1 (gr 3 - 19/02/10; gr 2 - 26/02/10) - 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 @).
  • 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 (plik tekstowy z komendami octave'a w kolejnych liniach - skrypt.m - wywołujemy z octave'a komendą skrypt)- dla rozwiązania równania cos(x)=0 na odcinku [0,2] - sprawdzić jak działa.
  • Zad 4 Przetestuj metodę bisekcji z pprzedniego 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
• Lab 2 (gr 3 - 5/03/10; gr 2 - 12/03/10) Równania nieliniowe cd.
  • Zad 1 Implementacja metody Newtona w octavie, następnie przetestować tę metodę dla następujących zadań z danym x0:
    1. x^2-2=0 x0=2;
    2. x^3-27 z x0=27
    3. exp(x)-2 z x0=10 i x0=1e3 czy x0=-1
    4. sin(x)=0, x0=2;
    5. cos(x)=0 x0=2;
    6. x^2-4x+4=0 x0=3; czy ogólniej (x-2)^k=0 z x0=3 dla k=2,4,5,16
    7. x^2-2=0 x0=1e6 i x0=1e-6.
    Czy szybkość zbieżności powyższych metod się różni? Wyświetlić błędy tzn e(n)=|xn-x^*| oraz e(n-1)/e(n) czy e(n-1)^2/e(n)- czy można określić rząd zbieżności na podstawie wyświetlanych danych?
  • Zad 2 Przetestuj działanie funkcji octave'a fsolve() dla równań z zadania 1, biorąc te same przybliżenia startowe. W szczególności porównaj szybkość działania z Twoją metodą Newtona (zarówno jeśli chodzi o czas jak i ilość iteracji)
  • Zad 3 Przetestuj działanie funkcji octave'a fsolve() dla równań z zadania 1, biorąc te same przybliżenia startowe ale dodatkowo podając pochodną - porównaj czy ilość iteracji uległa zwiększeniu w przypadku naszych zadań?
  • Zad 4 Implementacja metody siecznych w octavie, następnie przetestować tę metodę dla dla równań z zadania 1, biorąc te same przybliżenia startowe a za x_1 brać lekko zaburzone x_0 np. x_1=x_0*(1+1e-8) - porównaj czy ilość iteracji uległa zwiększeniu w przypadku naszych zadań? Wyświetl błąd e(n)=|xn-x^*| oraz e(n-1)e(n-2)/e(n) czy e(n-1)^q/e(n) (dla różnych q).
  • Zad 5 Rozwikływanie funkcji 2 zmiennych. Jak wykorzystująć m. Newtona (czy siecznych) rozwikłąć (w sposób przybliżony funckję tzn np mamy f(x,y)=x*x+y*y-1 chcemy znaleźć wartośći y(x) danej w sposób uwikłany f(x,y(x))=0 na danej siatce x_1,...,x_N i np. narysować wykres tej funkcji.
  • Zad 6 Odwracanie funkcji. Jak wykorzystująć m. Newtona (czy siecznych) narysować wykres funckji odwrotnej x=f(y) do danej funkcji y=f(x) na siatce równomiernej względem osi OY? (na równomiernej względem OX to b łatwe - jak?). Przetestować dla y=sin(x) dla y w [0,0.9]
• Lab 3 Arytmetyka FL i rozkład LU. (gr 3 - 19/03/10; gr 2 - 26/03/10)
  • Zad 1 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 2 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))|. (Oczywiście pomijając błędy fl f1=f2 - dlatego istotne jest aby liczyć funkcje wprost z tych wzorów )
  • Zad 3 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=2+2e-3 z warunkiem stopu że błąd jest mniejszy od 1e-20 (czyli długości odcinka w którym jest rozwiązanie ma być mniejsza od 2e-20). narysuj wykres tej funkcji na odcinku [2-h,2+2*h] dla h=1e-3,1e-4,1e-5.
  • Zad 4 Policzyć przybliżenie exp(x) z rozwinięcia w szereg y= sum_{k=0}^N 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łędy dla liczb ujemnych i dodatnich są tego samego rzędu? Jak to można łatwo poprawic?
  • Zad 5 Policzyć całki I_n= int_0^1 x^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 (mamy 1/(6*(n+1))
  • Zad 6 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ć czynniki rozkładu LU (funkcja lu()) oraz R^T R (R górnotrójkątna - chol()) czy są to macierze diagonalne?
• Lab 4 (gr 3 - 9/04/10; gr 2 - 16/04/10 - przesunięcie terminu Wielkanoc) 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. Do domu?
  • 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. Pewnie do domu.
• Lab 5 (gr 3 - 23/04/10; gr 2 - 30/04/10) LZNK i metoda Householdera. Zaliczenie pierwszego zadania komputerowego.
  • 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||_F=||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 ćwiczeń) 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. Sprawdzić np. dla wektorów u=(1,...,1)^T i v - o kierunku pierwszego wersora
  • 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)?
  < --li>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 (gr 3 - 14/05/10; gr 2 - 21/05/10 - przesunięcie terminu Juwenalia) Interpolacja Lagrange'a, FFT, splajny.
  • Interpolacja Lagrange'a: oblicz wykorzystując funkcję polyfit() znajdź wielomian 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 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 i tych wielomianów dla różnych N - np. używająć funkcji polyval() . Węzły Czebyszewa to zera T_{n+1}(t)=cos((n+1)arccos(t)) na [-1,1] odpowiednio przesunięte i przeskalowane.
  • Powtórz poprzednie zadanie, ale dla log(1+x) na odcinku [0,10] .
  • Interpolacja Lagrange'a 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 ?
  • Korzystając z funkcji octave'a spline() znajdź współczynniki splajnu interpolacyjnego kubicznego S_N na N węzłach równoodległych dla sin na odcinku [-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)| x_k pkt siatki. Następnie narysuj wykresy sin i tych splajnów dla różnych N . (Zadanie pomocnicze - jak obliczyć wartość splajnu otrzymanego z funkcji spline() ?)
  • 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 ?
  • Testy funkcji octave'a z algorytmem FFT fft() i odwrotna FFT ifft()(oblicza wartość transformaty odwrotnej do DFT). Sprawdź czy rzeczywiście operacja wykonywana przez fft() i ifft() są do siebie odwrotne i czy | F_N x |_2= C_N | x |_2 ze stałą C_N . Czy dla danego wektora f=(f_0, ldots,f_N) i wektora g =(f_1, ldots,f_N,f_0) mamy, że DFT od g jest DFT od f pomnożone przez ustalony wektor (po elementach wektora)?
  • Napisz funkcję szybko mnożącą dwa wielomianu stopnia leq N zadane w bazie potęgowej z wykorzystaniem DFT i odwrotnej DFT czyli dyskretnej tranformaty Fouriere'a. (jak na ćwiczeniach) Uwaga Często definicja DFT jest bez 1/N przed macierzą a wtedy macierz odwrotną do DFT należy przemnożyć przez N albo obie przez 1/ sqrt{N} . Sprawdź dla wielomianów p(x)=1+x i q(x)=2+x czyli p*q(x)=2+3x+x^2 a potem dla wielomianów jakiś naprawdę dużych stopni, np. q(x)=x^{50} a p(X)=x^{50}+1 .
• Lab 7 (gr 3 - 28/05/10; gr 2 - niestety nie ma zajęć ) FFT cd, wielomiany Czebyszewa, kwadratury - zaliczenie zadania komputerowego
  • Testy funkcji octave'a z algorytmem FFT:fft() i odwrotna FFT:ifft() (oblicza wartość transformaty odwrotnej do DFT). Sprawdź czy rzeczywiście operacja wykonywana przez te funkcje są do siebie odwrotne i czy ||fft( \vec{x})||_2= C_N||\vec{x}||_2 ze stałą C_N.
  • 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.
  • Policz za pomocą funkcji octave quad(): \int_{-1}^1 1\sqrt{1+x*x} T_n(x)T_m(x)\:d x dla m=2,n=3.
  • 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)), 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: E_{2N}\E_N, E_N=|\int_a^b f - T_N f| dla funkcji: sin(x) na [0,\pi] i N_0=5 czyli funkcji analitycznej, oraz x^{j+0.5} na [0,1] dla j=0,1,2 czyli funkcji w C^j([0,1])
    Porównaj wyniki obliczane kwadraturą trapezów z wynikami funkcji quad(). (Można sprawdzić dla jakiego n błąd obliczony złożoną kwadraturą trapezów jest na poziomie błędu funkcji octave quad() - i porównać wtedy ilość wywołań funkcji f przez obie procedury).

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

Skrypty m-pliki octave'a

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.