Semestr zimowy 2013/14

Matematyka obliczeniowa 2
Metody numeryczne (dla informatyków)

Metody numeryczne (dla informatyków)

Ćwiczenia/Lab MN

(co 2 tygodnie na zmianę)
Lista punktów z zadań komputerowych i kartkówek - aktualizowana w miarę na bieżąco - plik pdf

Zaliczenie ćwiczeń

Zaliczenie ćwiczeń tablicowych - serie zdań domowych - sprawdzane poprzez zapowiedziane 20min (1 do 2) kartkówki z 1-2 zadaniami z zadanych serii lub z zadania zrobionego na ćwiczeniach np. niedokończonego.
W razie usprawiedliwionej nieobecności na kartkówce proszę o kontakt po inne zadania zastępcze. (Dla osób którym mało p-tów zabraknie do zaliczenia ćwiczeń - będą być może dodatkowe serie zadane pod koniec semestru - ale z trudniejszymi zadaniami)
Zaliczenia labu - projekty zaliczeniowe (max. 1-2 projekty) 100%p-tów - projekt oddany po terminie 50% punktów labu)
Uwaga! Do egzaminu dopuszcza wg ustalen wykladowcy 25% z cwiczen i 25% pkt z labu.

Serie zadań domowych - projekty zaliczeniowe z labu

kartkówka - ostatni tydzien stycznia (lab lub cwicz odpowiednio)

Pierwsza seria - metody rozwiązywania równań nieliniowych - plik pdf zd1.pdf
Druga seria - metody rozwiązywania układów równań liniowych (uwarunkowanie, met. bezpośrednie i iteracyjne) - plik pdf
Trzecia seria - metody rozwiązywania LZNK, zadanie własne, interpolacja - plik pdf
Pierwszy projekt z labu - plik pdf - termin 4ty lab
Drugi projekt z labu - plik pdf - pierwszy lab w styczniu

Punkty za zadania domowe - projekty w labie (aktualizowane w miarę na bieżąco)- plik pdf

Bedę starał się umieszczac krótkie opisy tego co było czy ma być na labie czy ćwiczeniach.

Uwaga!

Zajecia beda wg USOSa - na zmiane - lab/cwiczenia tablicowe (choc wolne dni wprowadzaja zamieszanie).

Program ćwiczeń

Ćwicz 1 Skalarne metodu rozwiązywania równań nieliniowych: m. Newtona
1. Metoda Herona (czyli m Newtona dla równania x*x-a=0; a>0), że zbiega dla dowolnego x0; (elementarnie)
2. Odwracanie bez dzielenia (m. Newtona dla 1/x-a=0 dla jakich x_0 m. Newtona zbiega?).
3. 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).
4. 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.
5. (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)].
6. Zbadaj zachowanie iteracji metody Newtona dla funkcji f(x)=sgn(x)\sqrt{|x|} dla x_0 niezerowego.
Zadania nie zrobione na ćwiczeniach będą na następnych ćwiczeniach albo do zrobienia w domu.
Ćwicz 2 (18/25 października 2013) 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. dla metody x_k=F(x_{k+1}) mamy F(X%*)=x^* i F^{(k})(x^*)=0 k=1,..,s czy lokalnie można pokazać zbieżność z wykładnikiem s+1? (F gładka)
2. Jeśli f(x^*)=f''(x^*)=0\not = f'(x^*) metoda Newtona jest zbieżna lokalnie kubicznie.
3. Metoda iteracji prostych- pokaż że x_n=cos(x_{n-1}) jest zbieżne dla dowolnego x_0 . Jakie musi być n aby otrzymać przybliżenie x_n liczby x^*=cos(x^*) spełniające |x_n-x^*| < 1e-4 .
4. Metoda Newtona wielowymiarowa. Sprawdź czy dla funkcji f1(x,y)=x*x+y*y-1=0; f2(x,y)=x+y=0 m. Newtona będzie lokalnie zbieżna w otoczeniu zer tych funkcji. Policz pierwszą iteracje dla (x_0,y_0)=(1,0) . Napisz w pseudokodzie implementacje wiel. met. Newtona, zakładając, że dysponujemy funkcją x=linsolve(A,b) rozwiązującą układ równań liniowych Ax=b .
5. Jak rozwiązać układ 2 r. nieliniowych w którym w równaniu 2gim nie ma zależności od pierwszej zmiennej. Metoda Gaussa-Seidla i Jakobiego dla układu 2 równań
6. dla wyrażenie a^2-b^2 wypisz wynik obu agorytmów w arytmetyce fl: (1:x=a*a;2:y=b*b;3:w=x-y) i (1:x=a+b;2:y=a-b;3:w=a*b). oraz błąd bezwzględny.
Ćwicz 3
1. Zadanie na fl policz wzór na błąd bezwzględny wyniku algorytmu obliczania a^2-b^2: x=a*a;y=b*b;w=x-y wykonanego w arytmetyce zmiennopozycyjnej.
2. Optymalne stałe równoważności dla norm p-tych dla p=1,2,\infty.
3. Pokaż wzór na macierz indukowaną maksimum. Do domu na normę indukowaną pierwszą oraz wzory na note normy dla macierzy zespolonych.
4. 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.
5. Q ortogonalna. Pokaż, że \|QA\|_2=\|AQ\|_2=\|A\|_2 oraz to samo dla normy Frobeniusa.
6. Znajdź współczynnik uwarunkowania zadania obliczania iloczynu A*x ( A macierz, x wektor) ze względu na zaburzenie A (wzgldny)
7. Pokaz, ze cond(AB)<=cond(A)cond(B).
8. 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.
9. Eliminacja Gaussa dla macierzy trójdiagonalnej bez permutacji wierszy/kolumn i (do domu) z wyborem elementu głównego w kolumnie (permutacje wierszy).
Ćwicz 4
1. Znajdź rozkład Choleskiego dla A=[ 1, 0 , 1 , 1; 0 , 1 , -1 , -1; 1 , -1 , 3 , 2; 1 , -1 , 2 , 3]
2. Rozpatrzmy macierzy w postaci blokowej B=[A, b;c^T a] dla A macierzy mxm nieosobliwej której rozkład LU znamy, b,c wektorów wymiaru m i danego skalaru a. Jak możliwie niskim kosztem rozwiązać url z B (o ile B nieosobliwa)? Podaj warunek na to aby B byla nieosobliwa. (zad ze skryptu P.K.)
3. Rozpatrzmy macierzy w postaci blokowej G=[A, B;B^T, 0] dla A macierzy m\times m symetrycznej dodatnio określonej której rozkład LU znamy, B wymiaru m \times n o kolumnach tworzących układ liniowo niezależny. Czy ta macierz jest nieosobliwa? jeśli tak to jak układ Gx=b rozwiązać, przy założeniu że wpierw policzymy L czynnik rozkładu Choleskiego A i znajdziemy C takie, że LC=B . Policz koszt w zaleznosci od m/n. (zad ze skryptu P.K.)
4. (do domu zapewne) Udowodnij, że dla m. silnie diagonalnie dominującej wierszowo (|A_{k k}| > \sum_{k \not = l} |A_{k l}| \ \forall k ) LU można wykonać bez wyboru elementu głównego (czyli bez permutacji kolumn czy wierszy).
5. (do domu zapewne) Treść jak w poprzednim zadaniu ale dla macierzy symetrycznej dodatnio określonej. Wsk: wykorzystaj fakt, że jeśli A=A^T>0 to CAC^T też symetryczna i dodatnio określona o ile C nieosobliwa.
6. Pokaż że istnieje rozkład A=LDL^T dla D diagonalnej a L dolnotrójkątnej z 1 na diagonali (jednoznaczny o ile elementy D dodatnie na diagonali)
7. Udowodnij jednoznaczność rozkładu LU o ile istnieje (L dolnotr. z 1 na przekątnej, U górnotr. z elementami dodatnimi na przekątnej)
8. Pokaż że metoda Jakobi zbieżna jak macierz silnie diagonalnie dominująca. Podaj oszacowanie w normie maksimum.
9. Niech A=A^T>0 taka że jej wartości własne są w przedziale [c,d] . Znajdź \tau takie aby metoda Richardsona: x^{k+1}=x^k-\tau(Ax^k-b) była zbieżna do x^* rozwiązania Ax^*=b . Oszacuj szybkość zbieżności w normie \|\cdot\|_2 . Wyznacz wartości parametru dla których oszacowanie szybkości zbieżności w normie drugiej będzie najlepsze (przyjmując \lambda_{min}(A)=c <=\lambda_{max}(A)=d ).
10. Pokaż że jeśli A=A^T>0 to istnieje dokładnie jedna macierz symetryczna dodatnio określona B taka, że B^2=A oznaczana jako A^{1/2} . Wsk: istnieje Q ortogonalna taka że A=QDQ^T z D diagonalną.
11. Powtórz porzednie zadanie na oszacowanie błędu m. Richardsona ale w normie energetycznej \|x\|_A=\sqrt{x^TAx} . Wsk: Pokaż że \|I-\tau A\|_A=\|I-\tau A\|_2 korzystając z istnienia A^{1/2} takiej, że A^{1/2}A^{1/2}=A - poprzednie zadanie.
Ćwicz 5
1. Mamy uk. r. liniowych Ax=b dla A symetrycznej i >0 i C nieosobliwa taka ze cond_2(C^{-T}AC^{-1}) << cond_2(A). Stosujemy m. Richardsona z optymalnym parametrem do C^{-T}AC^{-1}y=C^{-T}b otrzymujac ciag y_k zbiezny do y. Pokaz ze przyjmujac Cx_k=y_k dostaniemy wzory na x_k wymagajace tylko mnozenia przez S^{-1} dla S=C^TC (zakladamy ze mamy procedure szybko mnozaca przez macierz S^{-1}) i A wektorów i dodawania wektorów. Pokaz ze norma bledu ||x_k-x||_A =||y_k-y||_{C^{-T}AC^{-1}}, zastosuj do oszacowania szybkosci zbieznosci bledu w tej normie zakladajac ze znamy cond_2(C^{-T}AC^{-1}) lub jego dobre oszacowanie. (zad ze skryptu P.K.)
2. Niech A=A^T>0 taka że jej wartości własne to \{1,2,...,10\} . Chcemy rozwiązać układ Ax^*=b przy pomocy metody Richardsona i znamy x^0 takie że x^0-x^* jest ortogonalne do wektorów własnych dla wartości własnych (10,9,8,7,6,5) . Czy dla \tau=1/3 metoda Richardsona: x^{k+1}=x^k-\tau(Ax^k-b) jest zbieżna do x^* ? Oszacuj szybkość zbieżności w normie \|\cdot\|_2 . Wyznacz wartości parametru dla których szybkość zbieżności w normie drugiej będzie największa przy taki założeniu.
3. Niech A nieosobliwa. Do rozwiązania układu A^TAx=A^Tb zastosowano następującą metodę iteracyjną: x^{k+1}=x^k+\alpha_k r^k \qquad r^k=A^T(b-Ax^k) dla \alpha_k minimalizującego funkcję g(\tau)=\|x^k+\tau r^k -x^*\|_2 . Oblicz wzór na \alpha_k ,sprawdź czy metoda popranie określona tzn. czy \alpha_k można policzyć bez znajomości x^* . Czy metoda jest zbieżna dla dowolnego x^0 i czy można oszacować \|x^k-x^*\|_2/\|x^0-x^*\|_2?
4. (do domu) Udowodnij że dla macierzy silnie diagonalnie dominującej metoda Gaussa-Seidla jest zbieżna.
5. Znajdź dla konkretnego wektora x wymiaru 3, taki wektor Householdera w, że odp. macierz Householdera przeprowadza x na wektor o kierunku pierwszego wersora. Sprawdź czy rzeczywiście Hx=(a,0,..,0)^T gdzie a>0 .
6. (do domu) Pokaż że jeśli A=A^T>0 to istnieje dokładnie jedna macierz symetryczna dodatnio określona B taka, że B^2=A oznaczana jako A^{1/2} . Wsk: istnieje Q ortogonalna taka że A=QDQ^T z D diagonalną.
7. Macierz A nieosobliwa znamy jej rozkład QR. Jak obliczyć możliwie tanio A^{-1} ?
8. 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) .
9. 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.
10. 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.
11. Czy rozkład QR macierzy nieosobliwej jest jednoznaczny o ile znaki elementów diagonali macierzy R są dodatnie?
12. 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}?
13. Macierz A kolumnami regularna mxn dla m>n . Znamy jej rozkład QR. Czy dualne LZNK tzn. LZNK dla macierzy A^T jest regularne? Czy ma rozwiązanie? Czy jeśli istnieje rozwiązanie to czy jest ono jednoznaczne? Jak wykorzystując ten rozkład policzyć rozwiązanie LZNK z A^T, a jeśli niejednoznaczne to rozwiązanie o najmniejszej normie?
Ćwicz 6
1. Dla macierzy A=[2,1;1,2] znajdź wektory i wartości własne. (można zgadnąć wektory...) Czy ciąg iteracji metoda potęgowa zbiegnie dla x_0=[2;1]? Czy jakiś jego podciąg zbiegnie? Co się stanie jak za x_0 weźmiemy [1;1-1e-6]? Czy odwrotna metoda potęgowa ((A-aI)x_k'=x_{k-1};x_k=x_k'/||x_k'||) zbiegnie dla parametru a=0.5 dla x_0=[2;1]?
2. Dla danej konkretnej tabelki np x_k :(-1,0,1,2) a y_k : (1, 0, 3,10) znajdź współczynniki wielomianu interpolacyjnego: w bazie jednomianów rozwiązując wprost układ równań liniowych: \sum_{k=0}^3a_kx_l^k=y_l l=0,1,2,3, w bazie Newtona związanej z \{x_k\}_{k=0}^2 przy pomocy algorytmu różnic dzielonych i metodą Newtona. Sprawdź wynik.
3. dla danych 3 węzłów (0,1,3) o krotności (1,3,1) i zadanych odpowiednich wartości funkcji i pochodnych (np dla f(x)=x^4) w tych węzłach znajdź współczynniki wielomianu interpolacyjnego Hermite'a : w bazie Newtona związanej z \{x_k\}_{k=0}^2 przy pomocy algorytmu różnic dzielonych. Sprawdź wynik.
4. Napisz w pseudokodzie odpowiednią wersję algorytmu Hornera obliczania wartości wielomianu zadanego w bazie Newtona dla znanych węzłów dla danego punktu x.
5. Niech S^{2r-2}_{2r-1}(\Pi) dla danego podziału \Pi:a=x_0\leq...\leq x_N=b będzie przestrzenią funkcji klasy C^{2r-2}(\mathbb{R}) i taką, że s_{|[x_k,x_{k+1}]}\in P_{2r-1} oraz s_{|(-\infty,a]},s_{|[b,\infty)}\in P_{2r-1} Wykaż, wzór na splajny w S^{2r-2}_{2r-1}(\Pi) dla r=2 : tzn. s\in S^{2}_{3}(\Pi) wtedy i tylko wtedy gdy s(x)=p(x)+\sum_{k=0}^na_k(x-x_k)_+^3 gdzie (x)_+= 0 dla x <= 0 i x dla x>= 0 a a_k współczynniki a p\in P_3 . Jeśli s splajn naturalny to dodatkowo p \in P_1 i \sum_{k=0}^n a_k x_k^j=0 dla j=0,1
6. Wyznacz wymiar przestrzeni S^{2r-2}_{2r-1}(\Pi) dla r=2 oraz jej podprzestrzeni splajnów naturalnych.
Ćwicz 7 (ostatnie)

Zadania, których nie zdążymy zrobić są do domu.

Program labów

Z ewentualnymi linkami do jakiś prostych skryptó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 (1015 : 18-10-2013 lub 1215: 25-10-2013) - metody rozwiązywania równań nieliniowych
- Zad 1 Zaimplementować metodę Newtona w octavie i przetestować jej zbieżnośc dla następujących funkcji: x*x-2 z x0=2, x*x*x-27 z x0=27, exp(x)-2 z x0=10,-10,100, sin(x) dla x0=2, (x-2)^k dla k=2,4,8,16 z x0=3, x*x-2 z x0=1e6, 1/x-a dla danych a=0.5,2,4,100 (oczywiście implementując bez dzielenia). Dla wszystkich tych funkcji znamy rozwiązania więc można wyświetlać na ekranie bład e_n=x_n-r (r - rozwiązanie) i |e_n|/|e_{n-1}|^p dla p=1,2,3. Dobierać różne wartości startowe x_0.
- Zad 2 Jak w zad 1 ale dla metody siecznych w szczególności przetestować czy zachodzi e_n/(e_{n-1}e_{n-2}) asymptotycznie zbiega do stałej i czy | e_n|/|e_{n-1}|^p dla p=(1+ \sqrt(5))/2 zbiega do stałej np dla x*x-2.
- Zad 3 Sprawdzić czy metoda iteracji prostych x_n=cos(x_{n-1}) zbiega do x=cos(x). Zbadac czy zbieżność jest liniowa.
- Zad 4 Porównać wyniki z zad 1 z funkcjami octave'a fsolve() i fzero() (pomoc: help fsolve/fzero).
- Zad 5 Zaimplementować przybliżoną metode Newtona w której pochodną przybliżamy ilorazem różnicowym tzn x_{n+1}=x_n - f(x_n)*h/(f(x_n+h)-f(x_n)) dla ustalonego h. Przetestować różne h np 1e-4,1e-7,1e-10 itp Porównać zbieżność z dokładną metodą Newtona (szczególnie ostatnie iteracje) dla funkcji z zadania 1 i z metodą siecznych.
- Zad 6 Rozwikływanie funkcji: dla funkcji y(x) zadanej równaniem g(x,y)=2x^2+3y^2-3=0 znaleźć wartości yk=y(x(k)) dla xk=k*h dla k=-N,...,N i h=1/N (N - 10,20,40,...) Znajdować yk rozwiązując układ równań g(xk,yk)=0. Jak w kolejnych krokach dobierać przybliżenie startowe w metodzie Newtona?
- Zad 7 Odwracanie funkcji: mamy daną funkcje np sin(x)+2*x znależć wartości funkcji odwrotnej na odcinku [0,5] na siatce k*h dla k=0,..,100. Narysować wykres funkcji. (Sam wykres można narysować dużo prościej bez wyliczania wartości funkcji odwrotnej. Jak?)
Lab 3 - wielowymiarowe układy równań nieliniowych i obliczenia w fl.
1. Zaimplementować wielowymiarową metodę Newtona w wersji z dokładnym Jakobianem i (do domu) w wersji gdy Jakobian przybliżany różnicami dzielonymi z parametrem h=1e-8. Zastosować do rozwiązania układu f1(x,y)=x+2y=1; f2(x,y)=3*x^2+y^2=1 dla różnych przybliżeń początkowych.
2. Wyznaczyć epsilon maszynowy czyli najmniejszą liczbę fl taką że po dodania jej do jeden dostajemy coś większego od jeden (w fl oczywiście) 2^{-t} dla t - ilości bitów mantysy. Porównać z eps! komendą octave'a. Można to zrobić i w C/C++. Czy wyszło to samo co w octavie (dla liczb typu double)?
3. Powtórzyć to zadanie w pojedyńczej precyzji tzn. wymuszając żeby zmienna była typu single (polecenie x=single(x) konwertuje do poj. precyzji)
4. Obliczyć wartości funkcji f1(x)=(x-2)^4 i f2(x)=x^4-.....+16 na siatce równomiernej na [2-a,2+a] w wektorze x dla a=1e-3. Narysować wykresy obu funkcji i policzyć błąd \| f1(x)-f2(x) \|_\infty czyli \max_k |f1(x(k)-f2(x(k))|.
5. (do domu) Zastosować bisekcję dla funkcji (x-2)^3 liczonej z wzoru na rozwinięcie dwumianu x^3-...-8 startując z a=2-1e-3 a b=1+2e-3 z warunkiem stopu że błąd jest mniejszy od 1e-20 (czy długości odcinka w którym jest rozwiązanie). Narysować wykres tej funkcji liczone z obu wzorów na odcinkach 2+[-h,2*h] dla h=10^{-3},...,10^{-5}.
6. Policzyć przybliżenie \exp(x) z rozwinięcia w szereg \sum_{k=0}^\infty x^k/k! (biorąc 100 a potem 1000 elementów szeregu) sprawdzić błąd względny |y-\exp(x)|/|\exp(x)| dla x od -100 do 100 (np. dla liczb różniących się o dziesięć czyli -100, -80,...,-10, 0, 10,...,80, 100. Czy błędy dla liczb ujemnych i dodatnich są tego samego rzędu? Jak to można łatwo poprawić?
7. Policzyć całki I_n=\int_0^1x^n/(5+x)d x dla n=0,..,20 z wzoru I_n+5*I_{n-1}=1/n raz biorąc I_0=\log(6/5) i iterując wprzód tzn dla rosnących n a raz biorąc I_{30}=1/180 i iterując w tył. (Oczywiście dla n=30,29,..) będzie kiepsko ale okaże się że potem jest coraz lepiej.
Lab 4 - zaliczenie projektu nr 1, 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. Sprawdzić dla konkretnej macierzy 3x3 np A=[1, 2, 3;-1, 0, 1;1, 1, 1] czy odpowiednie macierze górno/dolno trójkątne, permutacji czy ortogonalne uzyskane przy pomocy funkcji octave'a lu() i qr() rzeczywiście zwracają odpowiednie macierze rozkładu LU (PA=LU) lub rozkładu QR (A=QR) np. policzyć pierwszą normę różnicy A-QR dla rozkładu QR.
3. Przetestuj solver octave'a dla macierzy Hilberta H(n) (polecenie octave'a hilb() ją tworzy), tzn. stwórz macierz H(N) dla N=10,20,30 dla znanego rozwiązania np. stałego równego jeden na każdej pozycji czyli sol_k=1 , czyli policz H(n)*sol=f - rozwiąż w octavie układ z H(n) i f i policz normy (1,2 itd) \|H(n)x-f\| i \|x- sol\| dla x rozwiązania które wyliczył octave. Sprawdź czy macierz uzyskana przez inv(hilb(n)) i invhilb(n) są rzeczywiście odwrotne do hilb(n) - policz normę macierzową różnic tzn. np. inhilb(n)*hilb(n)-eye(n)
4. Stworzyć w octavie macierz trójdiagonalną A z 2 na diagonali i -1 na pod- i nad diagonalą tzn A_{i, i}=2 dla i=1,...,n i A_{i, i+1}=A_{i+1, i}=-1 dla i=1,..,n-1 a pozostałe elementy zero, zarówno przy pomocy funkcji octave'a diag() jak i wprost w pętli. Porównać czas używając tic i toc dla n=10,100,1000 . Policzyć uwarunkowanie macierzy dla różnych N , policzyć macierz odwrotną przy pomocy inv() (czy też jest trójdiagonalna?). Policzyć rozwiązanie układu równań liniowych Ax=b dla różnych znanych x za pomocą y=inv(A)*b i w=A\b (wczesniej obliczymy b=A*x). Porównać normy Ay-b,x-y,w-x (pierwszą, drugą, supremum) dla N=10,100,1e3,1e4 itp
5. Zaimplementować rozkład LU dla macierzy trójdiagonalnej bez wyboru elementu głównego. Przy czym do funkcji macierz przekazywać jako trzy wektory z przekątnymi. Zastosować do rozwiązania układu równań z macierzą z poprzedniego zadania. Porównać czas rozwiązywania Państwa algorytmem a algorytmem octave'a ( \ ) Rozwiązać układ dla losowego rozwiązania x=(1,...,1)^T dla N=10,100,1000 czy nawet 10000 . Policzyć błąd rezydualny (dla znanego losowego rozwiązania) i błąd rzeczywisty w różnych normach tzn \|x-\tilde{x}\|_2 i \|f-A*\tilde{x}\|_2 gdzie f=A*x a \tilde{x} rozwiązanie obliczone przez octave'a czy Państwa algorytm czy w najgorszy sposób A^{-1}*f (wcześniej musimy policzyć macierz odwrotną) .
6. 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.
7. W octavie przetestować eliminacje Gaussa z częściowym wyborem i bez wyboru dla macierzy A=[e, \; 1; 1, \;1] z e=eps/10 (tu eps funkcja octave'a zwracająca epsilon maszynowy) i wektorem prawej strony f=[1;0]^T .
8. (pewnie 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 .
9. (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.
Lab 5 Macierze rzadkie i iteracyjne rozwiązywanie równań liniowych.
1. Utwórz macierz trójdiagonalną T_N z a na diagonali i -1 na pod-nadiagonalach z wykorzystaniem funkcji sparse() ale nie tworząc macierzy pełnej tzn podając tylko elementy różne od zera.
2. Zapoznaj się z funkcją pcg() octave'a. Przetestuj metod cg (funkcja octave'a pcg()) dla macierzy wymiaru 25 losowych symetrycznych dodatnio okreslonych: B=rand(N,N), B=B'*B+pI dla różnych wartości p=0,5,10,20,100 tzn na ekranie rysuj wektory zwracane w RESVEC przez pcg() dla różnych wartości p.
3. Przetestuj metod cg (funkcja octave'a pcg()) dla macierzy wymiary T_N+pI dla różnych wartości p=0,5,10,20,100 i N=10,20,50, tzn na ekranie rysuj wektory zwracane w RESVEC przez pcg() dla różnych wartości p i N.
4. Zaimplementuj w octave metody Jakobiego i Gaussa Seidla - przetestuj ich zbieżność dla układu równań Ax^*=b dla macierzy A=B+pI gdzie B losowa a p dodatni parametr np. p=[0,5,10,20,30,100] . Tzn proszę dla każdej metody stworzyć wektor (\|b-Ax_k\|_2)_k i następnie narysować wykresy obu wektorów w jednym oknie. Czy widać różnice w szybkości zbieżności obu metod?
5. Rozpatrzmy macierz 61x61 o wartościach własnych: 1 (1 krotnej) i 99,100(30 krotnych) - można ją utworzyć np jako A=Q\Lambda Q^T dla Q macierzy ortogonalnej np. macierzy Householdera i \Lambda diagonalnej o wartościach na diagonali równym 1,99,100 Przetestuj zbieżność metody CG.
6. Wygeneruj jak w poprzednim zadaniu macierz wymiaru N=10 o wartociach wsnych 100,1,..,9. Nastepnie dla znanego rozwiązania Ax^*=b dobierz x^0 tak aby x^*-x^0 ortogonalne do pierwszej kolumny Q np x^0=x^*+Q[0;1;1] i powtórz test cg..
7. Testy metody Richardsona. Zaimplementuj metodę Richardsona i przetestuj jak w poprzednim zadaniu ale w normie drugiej dla różnych wartości parametru \tau. Oblicz w octavie maksymalną i minimalną wartość własną A i przetestuj czy zbieżność w normie drugiej jest taka jak w teorii. (dla znanego rozwiązania)
Lab 6 - LZNK, metoda Householdera, rozkład QR, metoda potęgowa i QR dla zadania własnego,
1. Testujemy funkcję qr() octave'a do znalezienia rozkładu macierzy A=[1, 1; 1, 0; 1, -1]. Następnie stosujemy ten rozkład do znalezienia rozwiązania LZNK z A i wektorem prawej strony b=[1;1;1] (powinno wyjść x=[1; 0; 0]), porównujemy z rozwiązaniem operatorem backslash. Stwórz macierz A^TA oraz wektor A^Tb - rozwiąż układ A^TAy=A^Tb - czy y==x? (module błędy zaokrągleń)?
2. Porównaj 2 metody rozwiązania LZNK dla wektora b i macierzy A m x n, 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 vander tworzy taką macierz dla m=n (tylko z odpowienio przepermutownymi kolumnami/wierszami) można stworzyć macierz mxm i wyciac odpowiednią podmacierz.
3. Zastosuj octave'a do rozwiązania zadania znalezienia krzywej zadanej równaniem a x*x+b*y*y=1 najlepiej pasującej do danych N punktów (x_k,y_k) (możemy np np wziąć lekko zaburzone punkty z danej elipsy y_k=\sqrt{1-4*x_k*x_k} + zab_k gdzie zab_k zaburzenie wylosowane z [0,1e-2] czyli mamy N punktów x_k np. x_k=1/k czy x_k=-1+h*k dla h=2/N k=1,...,N i z powyższego wzoru wyliczamy y_k +zab_k . Czy obliczone a i b bliskie 4 i 1 ?
4. Zaprogramuj funkcję octave'a \begin{lstlisting} function y=H(x,w,nw) \end{lstlisting} która dla danych wektorów \vec{x} i \vec{w} tego samego wymiaru N i skalaru nw=\|w\|_2 obliczy H_w \vec{x} dla H_w=I-2*\frac{1}{nw}\vec{w}\vec{w}^T czyli przekształcenia (macierzy) Householdera. (skalar może być parametrem opcjonalnym - jeśli funkcja będzie wywołana z dwoma tylko parametrami co można sprawdzić nargin<3 to tę normę można obliczyć w tej funkcji). Przetestować dla losowych wektorów \vec{x} i \vec{w} , sprawdzić czy \|H_w \vec{x} \|_2=\|\vec{x} \|_2. Można napisać ogólniejszą wersje gdzie x=X macierz N\times M i wtedy wynik H*X (jak to zrobić bez użycia pętli?). Sprawdzić czy \|A\|_2=\|H*A\|_2=\|A*H\|_2 i \|A\|_F=\|H*A\|_2=\|A*H\|_F dla losowej macierzy A i H przekształcenia Householdera dla losowego wektora w\not=0 wykorzystując tę funkcję.
5. (testujemy twierdzenie z wykładu) Znajdź wektor Householdera \vec{w} taki że odp przekształcenie Householdera przeprowadza dany wektor \vec{u}\not=0 na dany wektor l*\vec{v}\not=0 dla pewnej dodatniej stałej l . Przetestuj dla dowolnych 2 różnych wektorów czy rzeczywiście tak jest, tzn. czy H_w \vec{u}= l \vec{v} .
6. Zastosuj 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).
7. Zaprogramuj metodę Householdera rozwiązywania A\vec{x}=\vec{b} dla A macierzy trójdiagonalnej różnych wymiarów np. z 2 na diagonali a -1 na pod- i naddiagonalą i jakimś losowym wektorem rozwiązania tak aby koszt znalezienia rozkładu i rozwiązania układu równań liniowych był jak O(N) ?
8. Sprawdź jak działa eig(A) dla A=[2,1;1,2].
9. 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ą.
10. Przetestowc 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?
11. 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

Lab 7 - interpolacja, splajny; zaliczenie projektu nr 2 (gr o 1215) lub kartkówka (gr 1015)

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.
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.
Powtórz poprzednie zadanie, ale dla log(1+x) na odcinku [0,10] .
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] .
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().
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.
Interpolacja Lagrange'a przykład Rungego. Interpolacja Lagrange'a: wykorzystując funkcję polyfit() znajdź wielomiany interpolacyjne L_N f dla funkcji f(x)=1/(1+x*x) na [-5,5] dla N węzłów równoodległych i N węzłów Czebyszewa dla N=5,10,20,30 . Oblicz dyskretną normę max na siatce 1000 punktowej na tym odcinku tzn. \max|f(x_k)-L_Nf(x_k)| , gdzie x_k pkt siatki. Następnie narysuj wykresy f(x) i tych wielomianów dla różnych N - np. używając funkcji polyval(). Czy ciąg wielomianów interpolacyjnych zbiega do f ?
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() ?)
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 ?

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)

Tutaj link do stron Octave'a (skąd można ściągnąć kolejną dystrybucje - pod linuxa czy windows)

octave-forge

manual do octave'a w htmlu

Skrypty m-pliki octave'a

I termin - 2010/11
I termin - 2009/10
II termin - 2009/10
Zadania z kolejnych lat w skrypcie wykladowcy (prof. P. Kiciak)

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

Literatura dodatkowa

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.

http://people.maths.ox.ac.uk/trefethen/essays.html

wykłady i ćwiczenia

Powrót do mojej strony domowej.

Ostatnia aktualizacja: 14 stycznia 2014