Semestr zimowy 2014/15

1. Metody numeryczne (dla informatyków)
Ćwiczenia/Lab
• Program bieżącego labu
• Program bieżących ćwiczeń
2. Numeryczne równania rózniczkowe

Metody numeryczne (dla informatyków)

Ćwiczenia/Lab MN

Zaliczenie ćwiczeń i labu

Serie zadań domowych - projekty zaliczeniowe z labu

• Pierwsza seria - metody rozwiązywania równań nieliniowych - plik pdf zd1.pdf
• Druga seria - normy, metody rozwiązywania równań liniowych, LZNK - plik pdf zd2.pdf
• Trzecia seria - interpolacja i splajny - plik pdf zd3.pdf
• Czwarta seria - aproksymacja, kwadratury, zadanie wlasne, metody iteracyjna dla ukladu rownan liniowych - plik pdf zd4.pdf
• Pierwszy projekt z labu - lab1.pdf - termin 4ty lab po koniec listopada
• Drugi projekt z labu - lab2.pdf - pierwszy lab w styczniu

Program ćwiczeń

• Ćwicz 1 Skalarne metodu rozwiązywania równań nieliniowych: m. Newtona
  1. metoda iteracji prostej dla równania Keplera: x+0.5*cos(x)=10 tzn x_n=10-0.5*cos(x_{n-1}). Dla jakich x_0 zbiega, oszacowanie błędu dla x_0=10 i n=9.
  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?). W ilu iteracjach otrzymamy rozwiązanie dla a z (1,2) biorąc x_0=0.75 takie aby błąd względny był na poziomie podwójnej precyzji?
  4. czy iteracje x_{n+1}=cos(x_n) zbiegną dla dowolnego x_0?
  5. 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).
  6. 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.
  7. (metoda Halleya) Czyli metoda Newtona zastosowana do g(x)=[f(x)]/[\sqrt{f'(x)}]. Pokaż że druga pochodna funkcji g(x)=[f(x)]/[\sqrt{f'(x)}] wynosi zero dla x=a zera funkcji f o ile f'(a) nie =0. Pokaż, że wzór na tę metodę (tzn m Newtona dla g) 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 niezerowego.
  Zrobilismy 4 pierwsze zadania. Zadania nie zrobione na ćwiczeniach będą na następnych ćwiczeniach albo do zrobienia w domu.
• Ćwicz 2 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. Znajd najwieksza, najmniejsza liczbe fl (wsród unormowanych) dla arytmetyki w ktorej mamy 3 bity na ceche, 3 na mantyse i 1 na znak (bias=3). Podaj maksymalna i minimalna roznice pomiedzy kolejnymi liczbami. Powtorz zadani dla arytmetyki podwojnej precyzji.
  2. Uwarunkowanie zadania obliczania prostych funkcji jak np exp(x), (x-2)^2, 1-cos(x), 1-x*x itd
  3. Numeryczna poprawnosc algorytmu sumowania
  4. Numeryczna poprawno obliczania 1-a*a dwoma algorytmami: (1) f=a*a; f=1-f (2) f=1-a; f=f*(1+a) Czy oba algorytmy daja ten sam blad wzgledny o ile a jest reprezentowana dokladnie w fl?
  5. Policz uwarunkowanie oblicznia funkcji \sqrt{x}-\sqrt{x+1} dla x>0, który z alg: (1) wprost ze wzoru (2) z rownowaznego wzoru: 1/(\sqrt{x}+\sqrt{x+1}) jest lepszy z pkt widzenia fl dla x>>0?
• Ćwicz 3
  Normy wektorów i macierzy. Rozklad LU, L'L.
  1. Optymalne stałe równoważności dla norm p-tych dla p=1,2,\infty.
  2. Pokaż wzór na macierz indukowaną maksimum. Do domu na normę indukowaną pierwszą oraz wzory na te normy dla macierzy zespolonych.
  3. Norma Frobeniusa nie jest normą indukowaną. Z góry szacuje normę drugą. Znajdź stałą równoważności taką, że norma druga szacuje normę Frobeniusa.
  4. Q ortogonalna (Q'Q=I czyli \|Qx\|_2=\|x\|_2). Pokaż, że \|QA\|_2=\|AQ\|_2=\|A\|_2 oraz to samo dla normy Frobeniusa.
  5. Znajdź współczynnik uwarunkowania zadania obliczania iloczynu A*x ( A macierz, x wektor) ze względu na zaburzenie A (wzgldny)
  6. Pokaz, ze cond(AB)<=cond(A)cond(B).
  7. 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.
  8. Eliminacja Gaussa dla macierzy trójdiagonalnej bez permutacji wierszy/kolumn i (do domu) z wyborem elementu głównego w kolumnie (permutacje wierszy).
  Na cwiczeniach nr 4 planowana jest kartkówka.
• Ćwicz 4
  1. kartkówka (seria 1 i 2 bez QR i LZNK)
  2. (do domu)Znajdź rozkład Choleskiego dla A=[ 1, 0 , 1 , 1; 0 , 1 , -1 , -1; 1 , -1 , 3 , 2; 1 , -1 , 2 , 3]
  3. (do domu) Macierz A nieosobliwa znamy jej rozkład QR. Jak obliczyć możliwie tanio A^{-1} ?
  4. Znajdz wektor Householdera ze macierz Householdera przeprowadzi wektor [-1,1,1,1]' na [4,0,0,0]'. Sprawdz czy rzeczywiscie tak jest.
  5. Dla danych M punktów (x_k,y_k) chcemy znaleźć krzywą o równaniu a x^2+b y^2-1=0 najlepiej pasującą do tych punktów. tzn taką że \sum_k |a x_k^2+ b y_k^2-1|^2 jest minimalne. Sformułuj to zadanie jako LZNK. Co trzeba założyć o punktach aby to zadanie miało jednoznaczne rozwiązanie? Jaki byłby koszt rozwiązania tego LZNK przy pomocy rozkładu QR metodą Householdera jako O(M)? Rozwiąż to zadanie dla 3 punktów (1,0), (1,1) i (-1,0) .
  6. (do domu; lab?) Dla danych punktów (-1,1), (0,2), (2,4), (3,3) znajdź prostą y=ax+b najbliższą tym punktom w sensie: \sum_k|a x_k+b-y_k|^2 jest minimalne.
  7. (do domu) Niech macierz A trójdiagonalna nieosobliwa nxn . Opracować wersję algorytmu rozkładu QR metodą Householdera pozwalającą rozwiązać układ A x=b kosztem O(n) - o ile to możliwe albo uzasadnić dlaczego to niemożliwe.
  8. Czy rozkład QR macierzy nieosobliwej jest jednoznaczny o ile znaki elementów diagonali macierzy R są dodatnie?
  9. Niech macierz 4 X 3 A=[1,1,1;eps*I], (czyli pierwszy wiersz jedynki, podmacierz złożona z ostatnich trzech wierszy = eps*I ) i b=(3,eps,eps,eps)^T, eps - niezerowy parametr. Rozpatrzmy LZNK z A i b . Czy jest ono regularne? Policz układ równań normalnych. Jaki będzie skutek rozwiązania tego LZNK przy pomocy układu równań normalnych w arytmetyce pojedyńczej precyzji eps=1e-4}?
• Ćwicz 6 blad interpolacji L, wielomiany Czebyszewa, splajny, aproksymacja
  1. Pokaż \|u-L_1u\|_\infty\leq C_k(b-a)^k\|u^{(k)}\|_\infty k=1,2 dla interpolanta liniowego L_1u=u(a)+u[a,b](x-a) dla możliwie małych C_k (f \in C^k).
  2. Oszacuj blad w normie max dla f(x)=sin(x) lub sin(10*x) i L_n f wielomianu interpolujacego f na [-\pi,\pi] w 3 wezlach rownoodleglych.
  3. Ile trzeba wezlow Czebyszewa by na [0,4] blad interpolacji Lagrange'a dla \sin(x) byl co najwyzej 2^{-7}? Powtorz dla $\sin(4*x)$.
  4. Oszacuj blad dla interpolacji Lagrange'a na wezlach Czebyszewa dla \log(1+x) na [0,1] i [0,10]. Podaj te wezly.
  5. Znajdź splajn kubiczny B na podziale równomiernym [-3,3] z h=1 taki, że jest tożsamościowo zero poza (-2,2) i przyjmuje wartości B(-1)=B(1)=1 i B(0)=4 . Wsk: Splajn musi by f parzysta. Stad: interpolacja Hermite'a wpierw na [-2,-1] i analogicznie na [1,2] a potem na [-1,0] i [0,1] .
  6. Wyznacz na podziale równomiernym funkcję B_k taką , że jej nośnik to [x_{k-2},x_{k+2}] i przyjmującą w x_k wartość 4 . Wyznacz równania na s(x)=\sum_{k=-1}^{N+1}c_kB_k dla S splajnu naturalnego interpolującego daną f w (x_k)_{k=0}^N . Wykaż, że c_k wyznaczone jednoznacznie.(do domu)
  7. s splajn kubiczny taki, że na odcinku [x_{k-1},x_k] i [x_{k+1},x_{k+2}] jest równy zero. Pokaż, że na [x_k,x_{k+1}] też jest tożsamościowo zero. Czy tak własność jest prawdziwa dla dowolnych splajnów w S^{2r-2}_{2r-1}(\Pi) ?
  8. Pokaż oszacowanie błędu interpolacji splajnami liniowymi. na podziale równomiernym na [a,b] , tzn. jesli \Pi:\{x_k=a+k*h\} , h=(b-a)/N , to \|u- s_N u\|_\infty\leq C_kh^k\|u^{(k)}\|_\infty \quad k=1,2 dla stałych C_k niezależnych od h,a,b,N . Tu s_Nu splajn interpolacyjny liniowy.
  9. Znajdz wielomian najlepszej aproksymacji w L^2_g(-1,1) dla g(x)=1/\sqrt(1-x*x) dla st <= 2 dla f(x)=x^3-2x. Podaj blad aproksymacji. Wsk: znamy wiel. ortogonalne tzn wielomiany Czebyszewa sa ortogonalne z \|T_n\|_{L^2_g(-1,1)}=\pi/2 n>0 i \pi dla n=0
• Ćwicz 7
  1. Kwadratura Gaussa na 1 punkcie tzn kwadratura prostokątów: P_{[a,b]}f=Af(\theta) . Znajdź A,\theta dla którego rząd tej kwadratury dla I(f)=\int_a^bf największy. Oszacuj błąd dla takich A,\theta i f\in C^k z \|f^{(k)} \|_\infty \leq 1 dla k=1,2 , oraz pokaż, że I(f)-P_{[a,b]}f=C(b-a)^3f''(\xi) dla C stałej niezależnej od a,b,f
  2. Dla złożonej kwadratury prostokątów ( x_k=a+k*h , h=(b-a)/N ): P_Nf=h*\sum_{k=0}^{N-1}f(x_k+0.5*h) pokaż, że E_N(f)= \int_a^bf dx -P_N f=C h^2(b-a)f''(\xi) dla pewnego \xi\in(a,b) oraz C stałej niezależnej od h,a,b,f . Policz C oraz oszacuj błąd E_N(f) dla f klasy C^1 .
  3. Znajdz kwadrarture oparta na -1,a,1 o max rzedzie dla calki po [-1,1] z waga jeden.
  4. Znajdz kwadrarture oparta na -1,a o max rzedzie dla calki po [-1,1] z waga jeden.
  5. Pokaz ze jesli kwadratura na n+1 punktch ma rzad >=n+1 to jest to kwadratura interpolacyjna.
  6. Pokaz ze kwadratura Newtona-Cotesa oparta na n+1 punktach dla n+1 nieparzystego ma rzad >=n+2.

Program labów

• 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.
• Lab 2 - metoda Newtona i obliczenia w fl.
  1. (metoda Newtona) 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. 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. 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}.
  testnewton.m - m-plik z funkcja testnewton() - testujaca metode Newtona
  fl14.m - skrypt z rozwiazaniami niektorych zadan z fl z labu nr 2
• Lab 3 - 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 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)
  5. 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() .
  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ć 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
  7. 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ą) .
  8. (do domu) Używając inv() policz macierze dowrotne dla macierzy górnotrójkątnej U 2diagonalnej z 1 na diagonali i -1 na naddiagonali, potem znajdź macierz odwrotną do (U^{-1}-2I) (na diagonali zamiast 1 mamy -1). Policz uwarunkowania wszystkich macierzy w normie supremum.
  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. (do domu) 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. (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 .
  13. (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.
• Lab 4 - LZNK, metoda Householdera, rozkład QR, metoda potęgowa i QR dla zadania własnego,
  1. zaliczenie pierwszego projektu.
  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; 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 \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ę.
  6. (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} .
  7. 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).
  8. 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 5 - 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. 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 ?
• Lab 6 - interpolacja, splajny - bład, splajny,
  1. 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] .
  2. 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().
  3. 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.
  4. 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() ?)
  5. 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 ?
• Lab 7 Wielomiany Czebyszewa czyli aproksymacja jednostajna wieomianu st. n przez wielomiany st <=n, kwadratury złożone - prostokątów, trapezów. Zadanie wlasn - metoda potegowa i odwrotna potegowa, metoda QR. Macierze rzadkie i iteracyjne rozwiązywanie równań liniowych.
  1. 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 norme sup (dyskretna na 10000 pktów) na [-1,1] wielomianu x^{n+1}+\sum_{k<=n}a_kx^k z a_k losowymi i 2^{-n}T_{n+1}? Porównaj.
  2. Funkcja quad() w octave'ie - policz całkę z \sin(x) na [0,\pi] i [-1,1] .
  3. 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 .
  4. 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])
  5. Sprawdź jak działa eig(A) dla A=[2,1;1,2].
  6. przetestu metodę potęgową dla A=[2,1;1,2] z x0=[1;2]. Można użyć m-pliku: powerr.m z metodą potęgową.
  7. Przetestowac metode potegowa dla A 3x3 z wartociami wlasnymi -4,2,1 i pierwszym wektorem wlasnym [1;1;1] (pozostale jakiekolwiek) startujac z x0=[2;-1;-1] i wykonac ze 30 iteracji metody potegowej (wazne -aby x0 bylo prostopadle do wektora wlasnego dla najwiekszej co do modulu wartosci wlasnej).Czy x_n zbiegaja do wektora wlasnego dla 2?
  8. Zaimplementuj metode odwrotna potegowa (modyfikujac potegowa - ale bez liczenia macierzy odwrotnej explicite). przetestuj na macierzach z poprzednich zadan. Ewentualnie uzyj skryptu: invpow.m
  9. zaprogramuj metode QR i przetestuj ja dla macierzy A=[2,1;1,2]: A_0=A, Z_0=I, obliczamy A_k=Q_kR_k (rozklad QR) i A_{k+1}=R_kQ_k, Z_{k+1}=Z_kQ_k - ciag A_k powinien zbiec do macierzy diagonalej wiec warunek stopu norma ||A-diagonala(A)||_1<=tol dla tol=1e-6, czy sprawdzic Q_0...Q_k zbiega do macierzy ortgonalnej o kolumnach ktore sa wektorami wlasnymi A? Ewentualnei skorzystac ze skryptu: qreig.m

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)

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

Skrypty m-pliki octave'a

1. I termin - 2010/11
2. I termin - 2009/10
3. II termin - 2009/10
4. Zadania z kolejnych lat (do 13/14) w skrypcie o MN prof. P. Kiciaka.

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.