Semestr zimowy 2005/06

• wstęp - czym jest matematyka obliczeniowa
• podstawowe wiadomości o arytmetyce zmiennopozycyjnej - fl, definicje uwarunkowania względnego zadania i numerycznej poprawności algorytmu
• interpolacja wielomianowa Lagrange'a : postać w bazie Lagrange;'a  i bazie Newtona, różnice dzielone, błąd interpolacji przy max regularności oraz w przypadku ograniczonej regularności - korzystając z tw Jacksona z teorii aproksymacji
• interpolacja Hermite'a - postawienie zadania, uogólnione róznice dzielone, postać wielomianu Hermite'a z wykorzystaniem różnic dzielonych, błąd interpolacji Hermite'a
  
• interpolacja trygonometryczna (zespolona), algorytm szybkiej transformaty Fouriera (FFT) w wersji rekurencyjnej
• interpolacja splajnowa: baza splajnow kubicznych, oszacowania błędu w normie L2
• elementy aproksymacji  w  przestrzeniach unitarnych: def elementu najlepszej aproksymacji, rzut ortogonaln, baza ortogonalna, ukłąd równań normalnych (tzn z macierzą gramma w danej bazie), wielomiany ortogonalne w p. typu L2 - reguła trójczłonowa, w szczególności jako przykład  wielomiany Czebyszewa i ich własności, zastosowanie do interpolacji Lagrange'a
• aproksymacja jednostajna: zadanie, twierdzenie Czebyszewa o alternansie, jednoznaczność, związek z intepolacją Lagrange'a w węzłach Czebyszewa (tj odpowiednio przeskalowanch zerach odp. wiel Czebyszewa)
  
• kwadratury - przybliżone obliczanie calek - rząd kwadratury, kwadratury interpolacyjne, kw. Gaussa, kwadratury złożone m.in. złożona kwadratura trapezów (ewent. kwadratury Romberga)
  
• numeryczna algebra liniowa: bezpośrednie metody rozwiązywania układów równań liniowych m.in. eliminacja Gaussa, rozkład LU, rozkłady QR, uwarunkowanie zadania, liniowe zadanie najmniejszych kwadratów, układ równań normalnych związki z aproksymacją w przestrzeni  unitarnej (tj gdy mamy iloczyn skalarny)
  
• metody rozwiązywania równań nieliniowych (skalarnych): metoda bisekcji, metoda Newtona. Rząd zbieżności metody.  (ewent. metoda siecznych i metoda iteracji prostych (Banachowskich))
  
• zadanie własne (prawie na pewno zabraknie czasu)
  

Literatura:

• Lloyd N.Trefethen, Numerical analysis  (Princeton Companion to Mathematics, 2006, to appear) plik pdf
• Lloyd N.Trefethen, The definition of numerical analysis, (SIAM News, Nov 1992)  zgzipowany poscript.

1. Richard L. Burden,  J. Douglas Faires, Numerical Analysis, Wydanie 7, Brooks/Cole, 2001.
2. K.Moszyński, Metody numeryczne dla informatyków, skrypt, plik ps, 2002
3. L.Plaskota,  Dwanaście wykładów z matematyki obliczeniowej, skrypt, plik pdf, 2002
4. M.Dryja, J.M.Jankowscy, Metody numeryczne, WNT, 1982.

Szkic rozwiązań zadań z kolokwium

1. Tabelka różnic dzielonych, należalo uzupełnić i zapisać wielomian w bazie Newtona, zadanie interpolacyjne: znaleźć w wiel. st <=5 taki że w(1)=w(2)=w(3)=w(4)=w'(3)=0.5*w''(3)=1 (war. interpolacyjne Hermite'a)
2. Należało zbadać gładkość obu funkcji tzn f1(x)=x₊ⁿ i f2(x)=x_-ⁿsą w C^n-1(R)  (bo n-ta pochodna nieciągła) i następnie oszacować błąd interpolacji - w przypadku gdy p < n-1 z wzoru ||f -Lp f|| <=1/(4(p+1)) ||f^(p+1)|| h^p+1 tzn z wzoru na błąd interp. na węzłach równoodległych, w przypadku gdy p+1>n-1 z wzoru  ||f -Lp f|| <= (1 + 2^p) ||f -Q|| gdzie Q- wiel. stopnia p - a następnie oszacować ||f -Q|| z tw Jacksona. || || to norma supremum na [-1,1].   Błędy interp. w przypadku obu funkcji  są takie same - bo  f1(-x)=(-1)ⁿf2(x)  czyli też Lp f1(-x)=(-1)ⁿLp f2(x) (węzły są symetryczne względem zera, korzystamy z 1-1 zad interp. Lagrange'a)
3. Wystarczyło wyznaczyć tak ai b aby   funkcja byla  klasy C² - badając dla jakich a i b granice lewo i prawostronne w odpowiednich pktach siatki są sobie równe (2 warunki na a i b)
4. Można było skorzystać z tego że O(u * v)=O(u)O(v)  -gdzie O - transformata odwrotna tzn O(F(u))=F(O(u)), a * -oznacza dyskretny splot oraz z tego jak wygląda transformata przesuniętego ciągu (łatwe) - i koniec bo F(u v+)= F( O(F(u)) O(F(v))) =F(O(F(u)*F(v+)))=F(u)*F(v+)  a F(v+)_k=a_k F(v)_k -a_k należało policzyć. Ewentualnie wszystko można było wprost przerachować - to jedyne dość trudne moim zdaniem zadanie.
   

1. Arytmetyka FL i normy wektorów macierzy: Obliczyć wzór na normę pierwszą indukowaną macierzy rzeczywistej A=(a_i,j) wymiaru nxn. (17/10/05) Rozwiązanie: ||Ax||₁=\sum_i |\sum _j a(i,j)x(j) |<= \sum_i |\sum _j |a(i,j) | |x(j)| = \sum_j |x(j)| ( |\sum _i |a(i,j) | )<= \max_j ( |\sum _i |a(i,j) | )  [ \sum_j |x(j)|] =  \max_j ( |\sum _i |a(i,j) | )  ||x||₁, ponieważ max jest osiągane dla jakiegoś i₀ więc dla e_i0 wersora o tym indeksie mamy ||Ae_i0||₁=|\sum _i |a(i₀,j) | = \max_j ( |\sum _i |a(i,j) | ) =||A||₁.
2. Interpolacja: Niech a,b, różne ustalone punkty i w - ustalony dowolny wielomian stopnia n, pokazać że wtedy funkcja  f(x)=w[a,b,x] jest wielomianem (formalnie jest zdefiniowana dla x różnych od a lub b ale jako wielomian da się jednoznacznie przedłużyć do funkcji ciągłej). Określić stopień tego wielomianu  w zależności od n. (7/11/2005). Rozwiązanie: ze wzoru na błąd interpolacji Lagrange'a mamy: dla e(x)= f(x)-L₁f(x)=w[a,b,x](x-a)(x-b) gdzie L₁ wiel. interp. Lagrange'a stopnia <=1 t że f(a)=L₁f(a) i f(b)=L₁f(b), ale prawa strona tzn e(x) to wielomian stopnia n (albo e(x)=0  dla n<=1)  i a,b są miejscami zerowymi e(x) więc  istnieje wielomian st n-2 p(x)  (dla n<=1: p(x)=0) taki że e(x)=p(x)(x-a)(x-b), dzieląc przez (x-a)(x-b)  mamy że dla x różnych od a i b w[a,b,x]=p(x) wiel. st n-2 (w[a,b,x]=0 dla n<=1). 
   
3. Splajny: Niech a=x(0)<...<x(2N)=b, x(k)=a+k*h, dla h=(b-a)/(2N) (węzły równoodległe),    f - funkcja klasy C³taka że |f'''(t)|<1 dla t na [a,b];      s  funkcja ciągła taka że s(x(i))=f(x(i))  i=0...N,  oraz s na [x(2i),x(2i+2)] wielomian kwadratowy dla i=1,2,....N. Podać możliwie dobre oszacowanie błędu w normie sup na [a,b] w zależności od h.  (21/11/2005) Rozwiązanie: Należy zauważyć że s na odcinku [x(i),x(i+2)] jest wielomianem interp Lagrange'a interolującym f w pnktach x(i), x(i+1), x(i+2) - i - parzyste, więc wystarczy skorzystać z wzoru na błąd interp Lagrange'a i dostać że: ||f-s||_\infty <=(1/3)||f'''||_\infty h³ - wzór na błąd interp. L  z wykładu ( podzial równoodległy)
   
4. Aproksymacja w p. z iloczynem skalarnym: Powiemy że krzywa opisana równaniem  p(a,b,x,y)=ax*x+b*y*y-1=0  najlepiej pasuje do pktów (x_k,y_k) k=1,...,n jeśli \sum_k p(a,b,x_k,y_k))² = \min_{c,d} \sum_k p(c,d, x_k,y_k )² . Sformułować to zadanie jako zadanie aproksymacji w przestrzeni Rⁿ ze standardowym iloczynem skalarnym, sprawdzić co trzeba założyć o pktach (x_k,y_k) aby istniało jednoznaczne rozwiązanie, znaleźć  najlepiej pasującą krzywą dla danych: (0,1), (2,0),  (-1,0). (12 grudnia 2005). Rozwiązanie:  Zad. aproksymacyjne: znaleźć element najlepszej aproksymacji w=ax+by dla f=(1,..,1)^Tw V=Lin(x,y) dla  x= (x₁²,...,x_n²)^T,y=(y₁²,..,y_n²)^T, Rⁿ ze standardowym iloczynem skalarnym, wspólczynniki a,b to szukany wynik, jednoznaczność mamy gdy wektory x i y niezależne liniowo.  Dla konkretnych pktów - najprościej chyba rozwiązać układ A^TA(a,b)^T=A^T f  - gdzie A macierz o kolumnach x i y, ewentualnie można zortogonalizować  x i y ale wtedy dostaniemy rozwiązanie w nowym ukłądzie i należy powróić do wyjściowego układu.
5. Kwadratury: Dla funkcji f1(x)=|sin(x)|³ i f2(x)=|x-1|² zaproponować kwadratury Q1, Q2 wykorzystujące możliwie małe  ilości wartości funkcji  takie aby błąd pomiędzy odp. całką z tych funkcji na odcinku [0,10] (waga rho=1) a odp.  kwadraturą był na pewno  mniejszy niż 10^-4.  (uzasadnić) (9 stycznia 2006).  Rozwiązanie:  f1 jest w C² ale już nie w C³ na [0,10] więc stosując złożoną kwadrature trapezów dla której  mamy  że błąd jest jak  c/N² - stała c zależy tylko od f'' i  dł. odcinka, (dokładny wzór w każdej książce i pojawił się  na wykładzie) więc możemy wyliczyć  N - ilośc pktów, (ewent. możemy zastosować złożoną kw prostokątów o analogicznch własnościach), dla f2 -  kwadratura Simpsona (parabol) na 3 punktach daje dokładną wartość (bo jej rząd jest równy 4 a to wielomian stopnia 2, czy jakkakolwiek kwadratura o rzędzie co najmniej 3) - ewent. można zastosować złożoną kw. trapezów/prostokątów ale ilość punktów potrzebnych jest oczywiście spora - i to nie jest optymalne rozwiązanie.
   
6. Kwadratury:  znaleźć kwadraturę postaci Q f= Af(a) +Bf(b) taką że dla całki z wagą  1-x^2  na odcinku [-1,1] ma ona maksymalny rząd.  (znaleźć znaczy określić a,b,A,B). Ile ten rząd wynosi.  Podać oszacowanie błędu dla funkcji f takich że |f^(k)(t)| <M dla k=0,1,2,..,10. (10pkt- zadanie dodatkowe - ostanie ćwiczenia w styczniu)
7. Układy równań liniowych : mamy macierz nxn : A=H D gdzie H =I - (2/n)(u u^T) gdzie u= (1,1,..,1)^T  a D macierz diagonalna z diagonalą równa d=(1,2,3,4,..,n)^T.  Pokaż że A nieosobliwa. Oblicz uwarunkowanie A w normie 2.  Zaproponuj sposób rozwiązania układu równań z macierzą A możliwie małym kosztem tzn z możliwie małą ilością operacji arytmetycznych (chodzi o rząd wielkości tzn czy koszt będzie O(n) czy O(n²) czy O(n³) itp). Uzasadnij. (7pkt - zadanie dodatkowe - ostatnie ćwiczenia w styczniu)

1. Arytmetyka FL: Znaleźć takie t  naturalne że  fl(1+ 2^-t )=1 a fl(1+ 2^-t+1 )>1 dla podstawowych typów zmiennych rzeczywistych w jakimś języku programowania (może być Pascal czy C np dla C mamy 2 podstawowe typy double i float). W raporcie - tabelka z nazwą typu i obliczonym t.  (7/11/2005)
2. Interpolacja wielomianowa: obliczyć  wspólczynniki wielomianu Hermite'a  interpolującego funkcję 0.2*x⁵ - 2x  w punktach -2 i 2 z krotnością  2 - w raporcie podać współczynniki tego wielomianu w bazie Newtona (kolejnośc węzłów rosnąca)  oraz  naszkicować wykresy wielomianu interpolacyjnego  i  funkcji na jednym rysunku  na odcinku [-2.2,2.2] oraz wykresy ich pochodnych na drugim rysunku analogicznie - zaznaczyć pkty na wykresach o odciętych równym węzłom interpolacji. Podać jaki algorytm obliczania współczynników Państwo zastosowali (tzn jednym zdaniem napisać np że rozwiązywali państwo układ równań liniowych czy obliczali tabelę różnic dzielonych). Wykresy wielomianów  można  narysować używając gnuplot - standardowy program w dystrubucjach Linuxa, z linii komend wywołuje się go komendą gnuplot. (Pomoc- komenda help na linii komend gnuplota, funkcja gnuplota plot - rysuje wykresy - komenda help plot - podaje opis tej funkcji) (21/11/2005)
3. Aproksymacja średniokwadratowa: obliczyć a,b takie że wielomian w(x)=ax+b jest wielomianem najlepszej aproksymacji stopnia nie większego od jeden dla funkcji f(x) = 4+x -0.1*sin(x) w normie ||g||= (\sum_k=1³ k*g(k)*g(k))^0.5 ( \sum_k=1³  oznacza sumę od k=1 do 3 ) . W raporcie podać obliczone a i b, ||w-f||  oraz |a-1| i |b-4| oraz rysunek z wykresami w i f - zaznaczyć pkty na wykresach  o odciętych k=1,2,3 oraz bardzo krótki opis algorytmu (tzn czy rozwiązywali Państwo układ równań normalnych czy znaleźli wielomiany ortogonalne - jeśli wiel. ortog. to jak je Państwo znaleźli  itp). (5 grudnia 2005).
   
4. Aproksymacja jednostajna: obliczyć  a,b takie że  ||f - ax-b|| = min_(c,d) ||f-cx-d||  dla normy ||g||=max_k=1,2,3 |g(k)|   i  3  różnych funkcji:   f₁(x) = 4+x -0.1*sin(x), f₂(x)=(x-2)³, f₃=x²-4x.  Podać normy ||f_k-a_kx-b_k||, wyliczone wspólczynniki a_k,b_k k=1,2,3  wraz z krótkim opisem jak Państwo je policzyli. (19 grudnia 2005). Rozwiązanie: tu mamy 3 punkty i wielomian stopnia jeden czyli alternans skłąda się z 3 punktów tzn z (1,2,3) i wystarczyło rozwiązać odp. układ równań liniowych - co powinien zrobić program.
   
5. Kwadratury : obliczyć przybliżenie całek z funkcji:f(x)= x^k, k=0,1,2,3 i sin(x) i cos(x)  na odcinku [0,Pi]. Stosując złożoną kwadraturę prostokątów P_N(f)   z N=20, 40,80,160 - mamy  P_N(f) = \sum_i=0^N f(ih+0.5h)*h dla h=Pi/N. W raporcie drukować: E_N(f) i E_N(f)/E_2N(f) (o ile E_2N(f) różne od zera) dla odpowiednich N i f podanych powyżej. Tutaj E_N(f) = |P_N(f) - I(f)| dla I(f) = całka z f na odcinku [0,Pi].  (17 stycznia 2006)
6. Układy równań liniowych: Zaprogramować metode obliczania kosztem O(n) rozkładu macierzy A=LU dla A=(a_{ij})i,j=1..n macierzy nxn takiej że: a_ij=0 |i-j|>1 a_ii=2, a_{i,i+1}=a_{i+1,i}=-1  i<n która jest symetryczna i dodatnio określona (i trójdiagonalna). Testować  rozwiązując układ równań z tą macierzą    i znanym wylosowanym  rozwiązaniem dla n=10,50,100,500,1000 tzn dla danego n wylosować x wektor w Rⁿ, unormować ten wektor  (tzn obliczyć x= x/||x||) obliczyć b=Ax i rozwiązać układ równań Ax=b - dostajemy x* - obliczone przybliżone rozwiązanie.  W raporcie podać 1 tabelke:w wierszu tabelki podac: n, ||x- x^*||, ||b - Ax*||, dla || x||²=x^Tx - normy euklidesowej (można użyc jakiejś innej normy np normy max lub pierwszej). Skomentować wyniki - tzn czy rząd wielkości błedu (2 kolumna) i błedu residualnego (3 kolumna) są  różne? jeśli tak to dlaczego? Uzasadnić 1-2 zdaniami. (26 stycznia 2006 ) Komentarz: błąd residualny mniejszy bo uwarunkowanie macierzy jest złe dla dużego n (że norma A ok 4 to widać ale norma macierzy odwrotnej duża - Pańśtwo tego nie wiecie ale tylko to może dać tę różnicę).
   
7. Układy równań liniowych: Dla zadanego wektora u=(u₁,...,u₃₆₁) ^T  t. że u _i =(-1) ⁱ  dla i=1,..,361, napisać program który wykorzystuje możliwie mało pamięci i możliwie małym kosztem  umożliwia dla dowolnego wektora z  R³⁶¹ znalezienie jego obrazu przy przeksztalceniu Householdera H (przypomnienie: przekształcenie Housholdera H=Id - 2w*w ^T dla wektora Householdera w takiego że  ||w || ₂ =1) takiego  że Hu = a*e ₁  gdzie a =||u || ₂  a  e ₁ =(1,0,..,0)^T- wersor. W raporcie podać cztery   liczby:  normy drugą i pierwszą wektora  z=H x  dla x =(1,1,..,1)^T, normę pierwszą wektora Householdera w  oraz  z₁₂₇ (czyli z - wektor który jest obrazem wektora x przy przekształceniu H). Wszystko wykonać w arytmetyce double. Napisać jak  Państwo zaimplementowali  mnożenie przez H w szczeg. jak Państwo poradzili sobie z umieszczeniem w pamięci  H. (7pktów zadanie dodatkowe -  ostatnie ćwiczenia)
8. Układy równań liniowych: zaimplementuj rozkład QR macierzy A nxn metodą ortogonalizacji Gramma Schmidta. Przetestuj na macierzach Vandermonde'a dla n=5,10,20,40 dla xk równoodległych punktów na odcinku [0,1] tzn x k=(k-1)*h, k=1,..,nh=1/(n-1).  W raporcie(1) Czy obliczona Q rzeczywiście ortogonalna tzn podać normy Frobeniusa:  ||Q Q^T-I ||_F, ||Q^T Q-I||_F.  (2) Co dostajemy gdy zastosujemy rozkład do rozwiązywania Ax=b  czyli dla danego rozwiązania x - (można wylosować) - obliczamy b=Ax, nastepnie przy pomocy uzyskanego QR obliczyć przybliżone rozwiązanie y, i  policzyć ||x-y||₂/||x|||₂ i ||b-Ay||₂/(||A||_F.||x|||₂ ). ||B ||_F- norma Frobeniusa. Uwaga:ocenione będzie również sposób doboru testów w szczeg. x (czy równoważnie b)!  (10pktów zadanie dodatkowe na ostatnie ćwiczenia)