Semestr zimowy 2012/13

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

Metody numeryczne (dla informatyków)

Ćwiczenia/Lab MN

Zaliczenie ćwiczeń

Serie zadań domowych - projekty zaliczeniowe z labu

• Pierwsza seria - metody rozwiązywania równań nieliniowych - plik 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
• Czwarta seria - całkowanie czyli kwadratury złożone i Gaussa - plik pdf
• Projekt z labu - plik pdf - termin ostatni lab w styczniu

Program ćwiczeń

• Ćwicz 1 (8-10-2012) 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 nie zero.
  Zadania nie zrobione na ćwiczeniach będą na następnych ćwiczeniach albo do zrobienia w domu.
• Ćwicz 2 (24-10-2012) 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. Zera wielokrotne: pokaż że dla zera 2-krotnego met. Newtona jest zbieżna lokalnie 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).
  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 .
• Ćwicz 3 (7-11-2012)
  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 note normy dla macierzy zespolonych.
  3. Norma Frobeniusa nie jest normą indukowaną. Z góry szacuje normę drugą. Znajdź stałą równowążności taką, że norma druga szacuje normę Frobeniusa.
  4. Q ortogonalna. 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 .
  6. 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.
  7. 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.
• Ćwicz 4 (21-11-2012)
  1. Eliminacja Gaussa dla macierzy trójdiagonalnej bez permutacji wierszy/kolumn i (do domu) z wyborem elementu głównego w kolumnie (permutacje wierszy).
  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.
  3. Rozpatrzmy macierzy w postaci blokowej G=[A, B;B^T 0] dla A macierzy mxm symetrycznej dodatnio określonej której rozkład LU znamy, B wymiaru mxn o kolumnach tworzących układ liniowow niezależny. Czy ta macierz jest nieosobliwa? jeśli tak to jak układ Gx=b rozwiązać, przy założeniu że potrafimy rozwiązać układ postaci AC=B z C macierzą np. kosztem O(m^3) .
  4. 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. 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. Udowodnij jednoznaczność rozkładu LU o ile istnieje (L dolnotr. z 1 na przekątnej, U górnotr. z elementami dodatnimi na przekątnej)
• Ćwicz 5 (5-12-2012) Kartkówka (30min) 1-2 zadania z zadań domowych (z serii zadań tylko)
  1. Pokaż że metoda Jakobi zbieżna jak macierz silnie diagonalnie dominująca. Podaj oszacowanie w normie maksimum.
  2. 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 ).
  3. 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ą.
  4. (do domu) Powtórz porzednie zadanie 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.
  5. (do domu) 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.
  6. 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?
  7. (do domu) Udowodnij że dla macierzy silnie diagonalnie dominującej metoda Gaussa-Seidla jest zbieżna.
  8. 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 .
  9. (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ą.
• Ćwicz 6 (19-12-2012)
  1. Macierz A nieosobliwa znamy jej rozkład QR. Jak obliczyć możliwie tanio A^{-1} ?
  2. 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) .
  3. 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.
  4. 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.
  5. Czy rozkład QR macierzy nieosobliwej jest jednoznaczny o ile znaki elementów diagonali macierzy R są dodatnie?
  6. 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}?
  7. 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?
  8. 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]?
  9. 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.
• Ćwicz 7(16 I 2013)
  1. Kartkówka 15min - 1 zadanie - zakres: metody iteracyjne, LZNK do zadania wlasnego wlacznie. (czyli seria 2 zad 5-do konca i seria 3 zad 1-6)
  2. 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.
  3. Napisz w pseudokodzie odpowiednią wersję algorytmu Hornera obliczania wartości wielomianu zadanego w bazie Newtona dla znanych węzłów dla danego punktu x.
  4. 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 (do domu)
  5. Wyznacz wymiar przestrzeni S^{2r-2}_{2r-1}(\Pi) dla r=2 oraz jej podprzestrzeni splajnów naturalnych. (do domu)
  6. 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: Interpolacja Hermite'a wpierw na [-2,-1] i analogicznie na [1,2] a potem na [-1,0] i [0,1] .
  7. 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)
  8. Wyznacz wzory na współczynniki w bazie Newtona związanej z lewym końcem każdego pod-odcinka na obcięcie splajnu interpolacyjnego liniowego tzn. dla r=1 dla danej funkcji f .
  9. Wyznacz wzory na bazę splajnów liniowych taką, że każda funkcja z bazy jest zero w jednym węźle a zero w pozostałych. Wyznacz nośnik każdej takiej funkcji. (do domu)
  10. 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) ?
  11. Pokaż, że dla węzłów a=x_0<...< x_N=b zachodzi \|\Pi_k(x-x_k)\|_\infty \leq 0.25 *h^{N+1}N! dla h=\max_k|x_k-x_{k-1}| . Wskazówka: Indukcja. Zastosowanie - oszacowanie błędu interpolacji dla węzłów równomiernych. (do domu)
  12. Oszacuj błąd \|u-L_nu\|_\infty dla L_nu wielomianu interpolacyjnego interpolującego odpowiednio u(x)=\sin(x),\log(1+x) na odcinkach [0,2],[0,10] w węzłach Czebyszewa.
  13. 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 .
  14. 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.
  15. Pokaż, że jeśli s_n f splajn liniowy interpolujący f to \|f-s_nf \|_\infty\leq \|f-s \|_\infty\leq 2 \|f-s_nf \|_\infty \forall s\in S^0_1(\Pi)
      (do domu)
  16. 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
  17. 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 .

• Lab 1 (3 X 2012) - 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 (17-10-2012) - metody rozwiązywania równań nieliniowych
  1. 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.
  2. 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.
  3. Zad 3 Sprawdzić czy metoda iteracji prostych x_n=cos(x_{n-1}) zbiega do x=cos(x). Zbadac czy zbieżność jest liniowa.
  4. Zad 4 Porównać wyniki z zad 1 z funkcjami octave'a fsolve() i fzero() (pomoc: help fsolve/fzero).
  5. 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.
  6. 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?
  7. 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 (31-10-2012) - 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. (do domu) 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. (do domu) 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 (14-11-2012) - 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. (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 .
  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. 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 (28-11-2012) 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 (12-12-2012) - LZNK, metoda Householdera, rozkład QR
  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 takiej, że (A)_{k l}=(x_l)^{k-1} dla k=1,..,m, x_l=l*h l=1,..,n h=1/m. (Funkcja vander() 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 - pierwsa 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) ?
• Lab 7 - metoda potęgowa i QR dla zadania własnego, interpolacja, splajny
  1. Sprawdź jak działa eig(A) dla A=[2,1;1,2].
  2. 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ą.
  3. 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.
  4. 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.
  5. Powtórz poprzednie zadanie, ale dla \log(1+x) na odcinku [0,10] .
  6. Funkcje octave'a spline() and ppval. Zapoznaj się z tymi funkcjami ( help spline i help ppval). Wykorzystując te funkcje narysuj wykres splajnu kubicznego s_1 na podziale równomiernym odcinka [-3,3] z węzłami \{x_k=k\} dla k=-3,-2,...,3 - ręcznie podając losowe wartości przyjmowane przez te splajny w tych węzłach: np.: s_1(x_k)=(-1)^k , czy odpowiednio s_1(x_k)=1 . Następnie narysuj wykres splajnu s_2 o tych samych wartościach w węzłach ale podając dodatkowo wartości pochodnych w końcowych węzłach równe zero, tzn. wywołaj funkcje spline() podając dwie wartości więcej. Czy otrzymaliśmy te same splajny? Policz przybliżoną normę różnicy s_1-s_2 na odcinku [-3,3] .
  7. Splajn bazowy interpolacyjny - Dla danych węzłów równoodległych \{k\}_{k=-5,-4,...,5} na [-5,5] narysuj wykres splajnu kubicznego takiego, że s(0)=1 i s(k)=0 dla pozostałych węzłów k\not =0 oraz któ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 zadanie ale dla splajnu typu not-a-knot tzn. nie podając wartości pochodnych w końcach dla funkcji spline().
  8. Splajn kubiczny o minimalnym nośniku. Dla danych węzłów równoodległych \{k\}_{k=-5,-4,...,5} na [-5,5] narysuj wykres splajnu kubicznego takiego, że s(-1)=s(1)=1,s(0)=4 i s(k)=0 dla węzłów k\not \in\{-1,0,1\} oraz ma pochodne równego zero w węzłach skrajnych tzn. : -5 i 5 . Czy poza [-2,2] ten splajn jest równy zero? Policz przybliżone normy maksimum na [-5,-2] i [2,5] dla tego splajnu i narysuj jego wykres.
• Lab 8 - interpolacja, splajny - bład, kwadratury złożone - prostokątów, trapezów; zaliczenie projektu nr 1
  1. 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 ?
  2. 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() ?)
  3. 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 ?
  4. 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 .
  5. Funkcja quad() w octave'ie - policz całkę z \sin(x) na [0,\pi] i [-1,1] .
  6. 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 .
  7. 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])

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)

Skrypty m-pliki octave'a

1. I termin - 2010/11
2. I termin - 2009/10
3. II termin - 2009/10

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.