Semestr zimowy 2010/11

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

Metody numeryczne (dla informatyków) - wykład czwartki 1215-1345, sala 5440

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

1. Wykład 1 (7-10-2010) Wstęp - czym są metody numeryczne, kilka uwag ogólnych, plan wykładu, kilka pozycji literatury
2. Metody rozwiązywania równań nieliniowych (skalarnych) - zadanie czyli co chcemy obliczy, rzd zbienosi metody (zb. liniowa, kwadratowa) Wykład 2 (14-10-2010) metoda bisekcji - zbieżność z dowodem, metoda Newtona (stycznych) skalarna - wyprowadzenie i tw o zbieżności lokalnej kwadratowej z dowodem,
3. Metody rozwiązywania układów równań nieliniowych (wielowymiarowych) Wykład 3 (21-10-2010) - wielowymiarowa metoda Newtona, (metoda iteracji prostych - Banachowskich - na ćwiczenia) metoda Newtona z przybliżeniem pochodnej czyli przybliżona m. Newtona i m. siecznych (tylko przypadek skalarny), warunki stopu
4. Wykład 4 (28-10-2010) Arytmetyka zmiennopozycyjna - fl - podstawowe własności w tym redukcja cyfr znaczących przy odejmowaniu, definicje uwarunkowania zadania
5. Numeryczna algebra liniowa
   • Wykład 5 (4-11-2010) Bezpośrednie metody rozwiązywania układów równań liniowych (url) eliminacja Gaussa (EG) czyli rozkład LU, algorytm Choleskiego, EG z częściowym i pełnym wyborem elem. głównego Wykładu dnia 11 listopada 2010 nie ma - Święto Państwowe. Wykład 6 (18-11-2010) normy p-te wektorowe, normy macierzowe: indukowane p-te i Frobeniusa, uwarunkowanie zadania rozwiązywania url,
   • Regularne liniowe zadanie najmniejszych kwadratów (RLZNK) - istnienie i jednoznaczność rozwiązania - zastosowanie rozkładów QR do rozwiązywania RLZNK, Wykład 7 (25-11-2010) rozkład QR metodą Householdera zastosowanie rozkładów QR (uzyskanych poprzez m. Householdera) do rozwiązania reg. LZNK i jego szczególnego przypadku tj. układów równań liniowych.
   • Wykład 8 (2-12-2010) Wielkie układy równań liniowych w tym formaty macierzy rzadkich oraz iteracyjne metody rozwiązywania wielkich url: metody klasyczne (Jakobi, Gauss-Seidel, Richardson) z warunkiem dostatecznym zbieżności oraz Wykład 9 (9-12-2010) metody gradientowe na podstawie metody CG (bez wyprowadzenia) Zbieżność metody CG. Prekonditionery - idea prekonditionera dwustronnego i przykład prekonditionera diagonalnego dwustronnego.
6. Aproksymacja
   • Interpolacja wielomianowa Lagrange'a : istnienie i jednoznaczność rozwiązania, postać w bazie Lagrange'a i bazie Newtona w tym metoda Newtona Wykład 10 (16-12-2010) i algorytm różnic dzielonych znajdowania wiel. interp. Lagrange'a, błąd interpolacji przy maksymalnej regularności funkcji przyklad Rungego czyli ciąg wielomianów interpolacyjnych na węzłach równoodległych NIE musi byc zbieżny jednostajnie do funkcji którą interpoluje; węzły Czebyszewa jako optymalne węzły interpolacji
   • Interpolacja splajnowa (funkcjami giętymi) - splany liniowe interpolacyjny w bazie splajnow liniowych; oszacowania błędu Wykład 11 (13-01-2011) splajn interpolacyjny kubiczny w bazach B-splajnow kubicznych, oszacowania błędu w normie max
   • Kwadratury - przybliżone obliczanie całek - rząd kwadratury, kwadratury interpolacyjne w tym trapezów, Wykład 12 (20-01-2011) kwadratury złożone m.in. złożona kwadratura trapezów i Simpsona - wzory na błąd w przypadku gdy funkcja maksymalnej regularności oraz zbieżność w przypadku gdy funkcja tylko ciągła. Zera wielomianów ortogonalnych w przestrzeniach typu L² a konstrukcja kwadratur maksymalnego rzędu czyli kwadratur Gaussa.

Warunki zaliczenia:

ćwiczenia/lab

Egzamin

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.

Ćwiczenia/Lab MN

Zaliczenie ćwiczeń

Serie zadań domowych - projekty zaliczeniowez labu

• Pierwsza seria - plik pdf termin 3cie ćwiczenia tj. 5 listopada 2010
• Pierwszy projekt - plik pdf termin trzeci lab tj. 12 listopada 2010 (a jeśli 12/11/2010 dniem wolnym od zajęć to 19/11/2010)
• Druga seria - plik pdf zmiana terminu na szóste ćwiczenia tj. 17 grudnia 2010
• Drugi projekt - plik pdf termin piąty lab tj. 10 grudnia 2010
• Trzeci projekt - plik pdf termin ostatni lab tj. 21 stycznia 2011
• Trzecia seria - plik pdf termin ostatnie ćwiczenia tj. 14 stycznia 2011

Program ćwiczeń

• Ćwicz 1 (8-10-2010) Skalarne metodu rozwiązywania równań nieliniowych: m. Newtona Z1: Metoda Herona (czyli m Newtona dla równania x*x-a=0), że zbiega dla dowolnego x0>0 (na palcach) z2: Odwracanie bez dzielenia (m. Newtona dla 1/x-a=0 dla jakich x_0 m. Newtona zbiega?). Z3: 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). (Część dla dowolnej funkcji do domu - ewent. zrobimy na następnych ćwiczeniach)
• Ćwicz 2 (22-10-2010) Metodu rozwiązywania równań i układów równań nieliniowych cd. Metody bezpośrednie rozwiązywania uk. r. liniowych: Dokończenie Z3 z poprzednich ćwiczeń. Z1: Ile potrzeba iteracji m. bisekcji aby dla przedziału [3,4] metoda bisekcji zastosowana do sin(x) zwróciła przybliżenie x_n liczby \pi spełniające |x_n-\pi| < 1e-4. Z2: 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. Z3: Metoda Newtona wielowymiarowa. Sprawdź czy dla funkcji f1(x,y)=x*x+y*y-1=0; f2(x,y)=x+2y=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. Zadania domowe: ZD1 Pokaż że jeśli f(x)=f''(x)=0 z f'(x) różnym od zera to metoda Newtona dla równania f(x)=0 zbiega kubicznie.
• Ćwicz 3 (5-11-2010) Z1: Metoda eliminacji Gaussa jako rozkład LU (przypomnienie) na konkretnej macierzy. Z2: LU dla macierzy trójdiagonalnej. Z3: 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). Z4: Treść jak w poprzednim zadaniu ale dla macierzy symetrycznej dodatnio określonej. Wsk: pokaż, że jeśli A=A^T>0 to CAC^T też symetryczna i dodatnio określona o ile C nieosobliwa. Z5: Udowodnij jednoznaczność rozkładu LU o ile istnieje (L dolnotr. z 1 na przekątnej, U górnotr. z elementami dodatnimi na przekątnej) Zadanie Z3 do domu (w gr na 1015 również zadanie Z5)
• Ćwicz 4 (19-11-2010) Z1: 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. Z2:Udowodnij wzór na indukowaną normę pierwszą macierzy. Z3: Mamy dwa równoważne układu równań lin. A_k x= f_k k=1,2 z A_1=[1 1;-1 1] i A_2=[1e-3 1e-3;-1 1]. Wiemy że ||A_kx_k-f_k||_1=1e-6 dla pewnych wektorów x_k. Dla którego k prawdziwe jest oszacowanie ||x_k-x||_1 < 1e-4? Z4: Pokaż, że norma druga macierzy Q ortogonalnej (Q^TQ=I) wynosi jeden a norma druga QA równa się tej normie dla A. Z5: Rozpartrz w fl (przyjmując że fl(1+\eps)=1 i fl(-1/eps+1)=-1/eps) algorytm eliminacji Gaussa z i bez częściowego wyboru elementu głównego dla układu \eps*x+y=1 i x+y=0. Z6: Mamy macierz A=[a1,...,an] nieosobliwa. Znamy rozkład PA=LU. Jak możliwie niskim kosztem rozwiązać układ A1x=b z macierzą: A1=[b1,a2,...,an]? Przyjmujemy że b1 takie że A1 nieosobliwa. Z7: Pokaż że jeśli macierz kwadratowa górnotrójkątna R jest ortogonalna to R diagonalna z 1 lub -1 na diagonali. Wywnioskuj z tego że rozkład QR macierzy nieosobliwej (Q ortogonalna, R górnotrójkątna z elementami dodatnimi na diagonali) jest jednoznaczny.
• Ćwicz 5 (3-12-2010) Z1: Pokaż, że dla LZNK dla A kolumnami regularnej i wektora b - rozwiązanie tego LZNK spełnia: A^TAx=A^Tb (czyli tzw układ równań normalnych). Wsk: skorzystaj z istnienia rozkładu QR macierzy A. Z2: 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 a liczba niezerowa. Z4: Załóżmy że mamy macierz Q wymiaru mxn i macierz górnotrójkątną R nxn takie, że A=Q R. Pokaż, że A kolumnami regularna oraz że LZNK z macierzą A można rozwiązać rozwiązując odpowiedni układ równań liniowych z macierzą R. Z5: 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). Z6: 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. Z7 A kolumnami regularna mxn dla m>n. Znamy jej rozkład QR. Czy LZNK dla macierzy A^T jest regularne? Czy ma rozwiązanie? Czy jeśli istnieje rozwiązanie to jest ono jednoznaczne? Jak wykorzystując ten rozkład policzyć rozwiązanie LZNK z A^T? Z8: Niech 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. Z9: A nieosobliwa znamy jej rozkład QR. Jak obliczyć możliwie tanio A^{-1}? Z10 Czy rozkład QR jest jednoznaczny o ile znaki elementów diagonali macierzy R są dodatnie? Prawdopodobnie zdążymy zrobić max. 4-5 zadań - więc pozostałe do domu.
• Ćwicz 6
  1. 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ą.
  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$ oraz w normie energetycznej $\|x\|_A=\sqrt{x^TAx}$. Wyznacz wartości parametru dla których szybkość zbieżności w odpowiednich normach największa.
  3. Niech $A=A^T>0$ taka że jej wartości własne to $\{1,2,\ldots,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 pierwszych czterech wartości własnych. Znajdź $\tau$ takie aby metoda Richardsona: $x^{k+1}=x^k-\tau(Ax^k-b)$ była 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.
  4. Udowodnij że dla macierzy silnie diagonalnie dominującej metod Gaussa-Seidla jest zbieżna.
  5. Niech $A$ nieosobliwa. Do rozwiązania układu $A^TAx=b$ zastosowano następującą metodę iteracyjną: $$ x^{k+1}=x^k+\alpha_k r^k \qquad r^k=b-A^TAx^k $$ dla $\alpha_k$ minimalizującego funkcję $g(\tau)=\|A^TA(x^k+\tau r^k)-b\|_2$. Oblicz wzór na $\alpha_k$ oraz oszacuj $\|A^TAx^k-b\|_2/\|A^TAx^0-b\|_2$. Czy metoda jest zbieżna?
  6. 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.
  7. Napisz w pseudokodzie odpowiednią wersję algorytmu Hornera obliczania wartości wielomianu zadanego w bazie Newtona dla znanych węzłów dla danego punktu $x$.
• Ćwicz 7
  1. Pokaż, że dla węzłów $a=x_0<\ldots< 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.
  2. 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 i węzłach równoodległych.
  3. 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$.
  4. 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.
  5. Pokaż, że dla $s_N u$ splajnu liniowego interpolacyjnego dla $u\in C([a,b])$ widzimy, że $$ \inf_s\|s-u\|_\infty\leq \|s-s_N u\|\leq 2\inf_s\|s-u\|_\infty. $$
  6. 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$
  7. 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$.
  8. Kwadratura: $Qf=Af(\theta)$ dla $\int_0^1f(t)(t+1)\,dt$. Znajdź $A,\theta$ dla którego rząd tej kwadratury jest największy. Oszacuj błąd dla takich $A,\theta$ i $f\in C^k$ z $\|f^{(k)} \|_\infty \leq 1$ dla $k=1,2$.
Zadania, których nie zdążymy zrobić są do domu.

• Lab 1 (15-10-2010) - 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. Zad 4 Przetestuj metodę bisekcji z poprzedniego zadania do rozwiązania innych problemów: x^2-2=0, exp(x)=2, cos(10*x)=0 startując z odcinka [0,110] - przy dużej ilości zer na danym odcinku znalezione zero jest w sumie przypadkowe. Zad 5 Napisz funkcję octave'a w pliku bisekcja.m - której parametrami będą nazwa funkcji (czy funkcja - jak to zrobić odsyłam do manuala czy stron www), końce przedziału a i b, epsilon - żądana dokłądnośc bezwzględna (błąd ma być mniejszy od epsilona), max. ilośc iteracji. Funkcja ma zwrócić przybliżone rozwiązanie, ilość iteracji, oszacowanie błędu, i kod wyniku - 0 metoda zbiegła, 1 - metoda zatrzymała się z powodu przekroczenia max ilości iteracji, 2 - wartości w końcach odcinka funkcji mają ten sam znak. Funkcja ma działać również jak podamy tylko trzy czy cztery pierwsze argumenty - i wtedy max. ilość iteracji domyslnie ma być 100 , a epsilon 1e-5 (jak tylko 3 argumenty). Zad 6 znajdź w manualu funkcję octave'a -solve() - rozwiązująca równania nieliniowe - i przetestuj ją na kilku prostych przykładach skalarnych np cos(x)=0, x^2-2=0, x^10-2=0, x-sin(x)=0 itp
• Lab 2 (29-10-2010) - metody rozwiązywania równań nieliniowych cd.
  • Zad 1 Zaimplementować metodę Newtona w octavie i przetestować jej zbieżnośc dla następujących funkcji: x*x-4, x*x*x-27, exp(x)-2, sin(x) dla x=\pi, (x-2)^k dla k=2,4,8,16, 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ść liniowa.
  • Zad 4 Porównać wyniki z zad 1 z funkcją octave'a fsolve() (pomoc: help fsolve).
  • 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ą met 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 f odwrotnej. Jak?)
• Lab 3 (12-11-2010) - nieliniowe cd i obliczenia w fl.
  1. Zaimplementować wielowymiarową metodę Newtona w wersji z dokładnym Jakobianem i 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. 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))|.
  4. 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},\ldots,10^{-5}.
  5. 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,\ldots,-10, 0, 10,\ldots,80, 100. Czy błędy dla liczb ujemnych i dodatnich są tego samego rzędu? Jak to można łatwo poprawić?
  6. 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.
  7. 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,\ldots,n i A_{i, i+1}=A_{i+1, i}=-1 dla i
• Lab 4 (26-11-2010) - bezpośrednie metody rozwiązywania układów równań liniowych
  1. 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.
  2. 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)
  3. 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,\ldots,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
  4. 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,\ldots,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ą) .
  5. (do domu) Zaimplementować rozkład LU dla macierzy trójdiagonalnej z częściowym wyborem elementu głównego. I dalej testować jak w zadaniu 2. Zastosować do rozwiązania układu równań z macierzą z poprzedniego zadania A i z macierzą A-1.5 .
  6. 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 .
  7. (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 zad 1 i 2. Testować jak w zadaniu 2.
• Lab 5 (10-12-2010) - LZNK, metoda Householdera, iteracyjne metody rozwiązywania układów równań liniowych
  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. Pierwsza - operator backslash z octave'a, druga tworzymy 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 metodę Householdera do rozwiązania zadania znalezienia prostej 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,1e1,1e2,1e3] . (Funkcja \lstinline+rand(n)+ generuje wektor losowy o rozkładzie jednostajnym na [0,1] w octavie). Porównaj z wynikami otrzymanymi za pomocą standardowej funkcji octave'a tzn "\" oraz korzystając z funkcji qr(A).która zwraca rozkład QR macierzy ( Q ortogonalna).
  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 6 (7-01-2011) - macierze rzadkie, metody iteracyjne,
  1. Utwórz macierz trójdiagonalną 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. Przetestuj metod cg (funkcja octave'a pcg()) dla macierzy wymiary 2,4,8 (losowych symetrycznych dodatnio okreslonych -np. B=rand(N,N), B=B'*B) i dla macierzy z poprzedniego punktu z a=2,10 i dla N=2,8,16,64,256.
  3. Zaimplementuj w octave metody Jakobiego i Gaussa Seidla - przetestuj ich zbieżność w normie $\|\cdot\|_\infty$ dla układu równań $Ax^*=b$ dla macierzy $A$ trójdiagonalnej z $a$ na diagonali a $-1$ na pod i nad-diagonali dla różnych $a$ np $a=-10,-5,-1,0,1,2,10,100$. Tzn proszę drukować na ekranie $\frac{\| x^k-x^*\|_\infty}{\| x^{k-1}-x^*\|_\infty}$.
  4. Testy metody Richardsona. Zaimplementuj metodę Richardsona i przetestuj jak w poprzednim zadaniu ale w normie drugiej dla $a=2,20,100$ 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 i normie energetycznej jest taka jak w teorii. (dla znanego rozwiązania)
  5. Rozpatrzmy macierz $3x3$ o wartościach własnych: $1,99,100$ - można ją utworzyć np jako $A=Qdiag([1,99,100])Q^T$ dla $Q$ macierzy ortogonalnej np. macierzy Householdera. 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]$. Przetestuj zbieżność metody Richardsona dla różnych $\tau$.
• Lab 7 (21-01-2011) - interpolacja, splajny, kwadratury złożone - prostokątów, trapezów i Simpsona; zaliczenie projektu nr 3
  1. Interpolacja Lagrange'a: 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.
  2. Powtórz poprzednie zadanie, ale dla $\log(1+x)$ na odcinku $[0,10]$.
  3. Interpolacja Lagrange'a przykład Rungego. Wykonaj poprzednie zadanie, ale dla $f(x)=1/(1+x*x)$ na $[-5,5]$. Czy ciąg wielomianów interpolacyjnych zbiega do $f$?
  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$?
  6. 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$.
  7. Funkcja $quad()$ w octave'ie - policz całkę z $\sin(x)$ na $[0,\pi]$ i $[-1,1]$.
  8. 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$.
  9. 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,\ldots$ 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