Semestr letni 2011/12

• Ćwiczenia/Lab
• Program bieżącego labu
• Zadania domowe, projekty na lab.
• Wyniki zadań domowych, projektów z labu, kolokwium etc: plik-pdf
• Wyniki egzaminu II termin z % punktów z labu i ćwiczeń: plik-pdf
• Wyniki egzaminu I termin z % punktów z labu i ćwiczeń: plik-pdf

Wykład

• osoby które mają USPRAWIEDLIWIENIE za nieobecność na kolokwium
• osoby które uzyskały mniej niż 15 pktów z kolokwium - przy czym w tym wypadku - otrzymają one max 15 pkt (nawet jak rozwiążą wszystkie zadania) czyli min(x,15) - x ilość uzyskanych pkt z 3 zadań po 10pkt.

Program wykładu

1. Wstęp - czym jest matematyka obliczeniowa (0.5 wykładu).
2. Metody rozwiązywania równania nieliniowego (skalarnego) (1-2 wykłady)
3. Arytmetyka fl (1 wykład)
4. Numeryczna algebra liniowa - rozwiązywanie układów równań liniowych, LZNK (liniowe zadanie najmniejszych kwadratów), numeryczne zadanie własne (5-6 wykładów)
5. Interpolacja wielomianowa Lagrange'a i Hermite'a, interpolacja splajnowa (3-4 wykłady)
6. Numeryczne całkowanie i różniczkowanie (2 wykłady)

Zasady zaliczenia przedmiotu

Ćwiczenia/Lab MO

• Gr 1 ćwiczenia piątek 830-10, sala 5070. Pierwszy lab 17 luty i potem co 2 tyg
• Gr 2 ćwiczenia - piątek 1215-1345, sala 5070. Pierwszy lab 24 luty i potem co 2 tyg

Zaliczenie ćwiczeń i labu

1. Pierwsza seria ZD równania nieliniowe
2. Druga seria ZD arytmetyka fl
3. Pierwszy projekt z labu testami metody Steffensena na piątek 1015: 20 kwietnia/27 kwietnia 2011 (w terminie labu danej grupy)
4. Trzecia seria ZD - układy równań liniowych, LZNK, normy macierzy, zadanie własne
5. Drugi projekt z labu (liniowa regresja czyli LZNK i metoda Householdera) ostatnie zajęcia w labie w danej grupie
6. Czwarta seria ZD interpolacja i kwadratury
7. Zadania z kolokwiów (2gr i poprawkowe) i egzaminu w I terminie
Projekt z labu trzeba pokazać w akcji w czasie zajęć labu, za zadania oddane po terminie otrzymują Państwo 50% punktów.

Program labów

• Lab 1 (gr 1 - 17/02/10; gr 2 - 24/02/10) - 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). Rysowanie wykresów funkcji - czyli funkcja plot(). Wskaźniki do funkcji (handle - operator @).
  1. Operacje macierzowe
     • Utwórz dowolne macierze 3x4 A i 3x5 B - a następnie macierz 3x8 C, której pierwsze 4 kolumny to A a kolejne to B.
     • Z macierzy C '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.
     • Zamień kolejność wierszy D.
     • Wytnij dolnotrójkątną część macierzy D a potem górnotrójkątną
     • Wstaw D z powrotem do C jako główny minor.
     • Policz \sin(D)=(\sin(D_{ij}) od D.
     • Zapisz D do pliku (binarnego i ASCII) - zamień element D(1,1) na -100 i wczytaj nową macierz do octave'a.
  2. Policz dyskretną normę max od (\sin(x))^2 na [0,1] (wektorowo- czyli bez użycia pętli).
  3. Narysuj wykres \sin() na odcinku [0,4].
  4. Znajdź przybliżone max i min funkcji x*(3+2*cos(x)) na [-1,5] oraz przybliżenia punktów ekstremalnych. Zrób to samo dla jakiegoś wielomianu stopnia dwa i trzy np x^3+x^2-x-4
  5. Utwórz funkcję w m-pliku np. obliczającą (\sin(x))^2
  6. Utwórz m-plik obliczający wartość funkcji z parametrem np \sin(a*x) - parametr przekaż jako zmienną globalną
  7. Funkcje anonimowe - utwórz funkcję (\sin(x))^2 jako funkcję anonimową - czyli w linii komend
  8. Przy pomocy pętli for i while oblicz sumę 1+...+N dla N=100. Użyj if aby sprawdzić czy równa się tyle ile powinno 0.5*N*(N+1)- wyprowadź na ekran komunikat używając printf().
• Lab 2 (gr 1 - 29/02/12; gr 2 - 7/03/11) - równania nieliniowe
  1. Zaimplementuj metodę bisekcji w skrypcie - dla rozwiązania równania cos(x)=0 na odcinku [0,2] - sprawdź jak działa.
  2. 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 przypadkowo wybrane
  3. Napisz funkcję octave'a: function [y,iter,osz,kod]=bisekcja(nazwa_funkcji,a,b,tola,maxit) w pliku bisekcja.m - której parametrami będą nazwa funkcji (czy funkcja - jak to zrobić odsyłam do pomocy czy stron www), końce przedziału a i b , tola - żądana dokładność bezwzględna (błąd ma być mniejszy od epsilona), maxit - maksymalna ilość iteracji. Funkcja ma zwrócić y - przybliżone rozwiązanie, iter - ilość iteracji, osz - oszacowanie błędu , i kod-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 maksymalna ilość iteracji domyślnie ma być 100, a tolerancja 1e-5 (jak tylko trzy argumenty).
  4. Znajdź w pomocy funkcję octave'a - rozwiązująca równania nieliniowe - i przetestuj ją na kilku prostych przykładach skalarnych np.
     [ f1(x)=]
     1. cos(x)=0 ,
     2. x^2-2=0 ,
     3. x^{10}-2=0 ,
     4. (x-1)*(x-2)*(x-3)=0 ,
     5. x-\sin(x)=0 .
     A dokładniej zastosuj tę funkcję do rozwiązywania powyższych równań dla różnych przybliżeń startowych.
  5. Zaimplementować metodę Newtona w octavie i przetestować jej zbieżność dla następujących funkcji:
     [ f1(x)=]
     • x^2-2 z x_0=2 ,
     • x^3-27 z x_0=27 ,
     • \exp(x)-2 z x_0=10,-10,1e3 ,
     • \sin(x) dla x_0=2 ,
     • cos(x) z x_0=2
     • (x-2)^k x x_0=3 dla k=2,4,8,16 ,
     • x*x-2 z x_0=1e6 , (czy zbiega w ogóle a jak tak to czy kwadratowo)
     • 1/x-a dla danych a=0.5,2,4,100 (oczywiście implementując bez dzielenia) jak dobrać x_0 ?
     Dla wszystkich tych funkcji znamy rozwiązania więc można wyświetlać na ekranie błąd e_n=x_n-r ( r - rozwiązanie) i |e_n|/|e_{n-1}|^p dla p=1,2,3 . Proszę dobierać różne wartości startowe x_0 .
  6. Powtórzyć poprzednie zadanie ale zastępując metodę Newtona przez metodę 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 .
  7. Sprawdzić czy metoda iteracji prostych x_n=cos(x_{n-1}) zbiega do x=cos(x) . Zbadać eksperymentalnie czy zbieżność jest liniowa, tzn. czy |x-x_n|/|x-x_{n-1}| zbiega do stałej.
  8. Porównać wyniki z zadania poprzedniego z funkcją octave'a fsolve() (pomoc: help fsolve ).
  9. Zaimplementować przybliżoną metodę 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. h=1e-4,1e-7,1e-10 itp Porównać zbieżność z dokładną metodą Newtona (szczególnie ostatnie iteracje) dla funkcji z poprzednich zadań.
  10. \textsl{Rozwikływanie funkcji:} dla funkcji y(x) zadanej równaniem g(x,y)=2x^2+3y^2-3=0 znaleźć wartości y_k=y(x_k) dla x_k=k*h dla k=-N,...,N i h=1/N (N - 10,20,40,...) Znajdować y_k rozwiązując układ równań g(x_k,y_k)=0 . Jak w kolejnych krokach dobierać przybliżenie startowe w metodzie Newtona?
  11. Odwracanie funkcji: rozpatrzmy daną funkcje np. f(x)=\sin(x)+2*x znaleźć wartości funkcji odwrotnej na odcinku [0,5] na siatce k*h dla k=0,..,100 . Narysować wykres tej funkcji. (Sam wykres można narysować dużo prościej bez wyliczania wartości funkcji odwrotnej. Jak?)
  12. Zaimplementować wielowymiarową metodę Newtona w wersji z dokładnym Jakobianem i w wersji gdy jakobian przybliżany różnicami dzielonymi z parametrem h. 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. Taką funkcje w octavie definiujemy:
      function y=F(x)
      #x- wektor pionowy
      y=[[1 2]*x-1;3.*x(1)*x(1)+x(2)*x(2)-1];
      endfunction
• Lab 3 (gr 1 - 16/03/12; gr 2 - 23/03/11) - dokończenie porzedniego labu i arytmetyka fl
  Celem labu jest wykonanie odpowiednich zadań, które wykażą pewne cechy arytmetyki fl. Domyślnie w octavie wykorzystywane są liczby zmiennopozycyjne o podwójnej precyzji, jakkolwiek w najnowszych wersjach octave'a można również trochę sztucznie wymusić używanie zmiennych w precyzji pojedynczej przy pomocy funkcji single(a)! zwracającej zmienną pojedynczej precyzji z tą samą wartością. Z konstrukcji arytmetyki fl wynika, że odejmowanie dwóch wartości o tym samym znaku o małej różnicy może skutkować dużą utratą dokładności.
  1. Funkcje single(a) i eps. Wywołaj pomoc do tych funkcji w octavie. Sprawdź czy 1+eps obliczone w arytmetyce podwójnej precyzji w octavie jest większe od jeden. Przyjmij, że a=single(eps) i sprawdź czy ponownie 1+a jest większe od 1 .
  2. Wyznacz 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)? A co dla liczb w pojedynczej precyzji? Funkcja octave'a single(x) tworzy takie zmienne. Wykorzystując tę funkcję ponownie wyznacz epsilon maszynowy ale w pojedynczej precyzji.
  3. Narysuj wykres funkcji f(x)=(1+a)-a na [0,1] dla różnych a=1ek dla k=1,2,..,20 . Tutaj ważne aby obliczać wartość f(x) z tego wzoru
  4. Policzyć f(x)=x-\sqrt{1+x*x} raz algorytmem wprost wynikającym z tego wzoru a raz z wykorzystaniem równoważnego wzoru f(x)=\frac{-1}{x+\sqrt{1+x*x}} tzn. Algorytm 1 a= \sqrt{1+x^2};
     w_1=x-a Algorytm 2 a=\sqrt{1+x^2};
     w_2=\frac{-1}{x+a} dla x=10^k i k=4,...,10 . Czy widać różnicę w wyniku? Powtórzyć w arytmetyce pojedynczej precyzji tzn. z wykorzystaniem funkcji single(x).
  5. Obliczyć wartości wielomianu f_1(x)=(x-2)^4=f_2(x)=x^4-.....+16 dwoma algorytmami:
     Algorytm 1
     a= (x-2), f_1(x)=a^4;
     Algorytm 2
     f_2(x)=x^4-...+16
     na siatce równomiernej 1000 punktowej na [2-a,2+a] dla a=1e-3 . Narysować wykresy obu funkcji i policzyć błąd || f_1(x)-f_2(x) ||_\infty czyli \max_k |f_1(x(k)-f_2(x(k))| . Tu x(k)=2-a+k*h dla h=a/500 . Wektor x można utworzyć w octavie przy pomocy funkcji linspace().
  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 n=0,..,20 dwoma algorytmami ze wzoru I_n+5*I_{n-1}=1/n . Pierwszy algorytm przyjmuje I_0= log(6/5) i oblicza z powyższego wzoru kolejne I_n dla n=1,2,3,... . Drugi algorytm wykorzystuje fakt, że \frac{1}{(n+1)6} <= I_n<= \frac{1}{(n+1)5} zatem w tym algorytmie przyjmujemy za I_{30} jakąkolwiek wartość z tego przedziału np. I_{30}=1/180 i iterujemy w tył, tzn. dla n=30,29,...,20,...,0 . Porównaj wyniki obu algorytmów dla n=0,...,30 oraz czy wyniki spełniają oszacowanie. Dlaczego jeden z algorytmów działa zdecydowanie lepiej w arytmetyce fl? Jako dodatkowe zadanie pozostawimy uzasadnienie wzoru rekurencyjnego i oszacowania I_n wykorzystywanych w algorytmach.
  8. Zastosować algorytm bisekcji 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 . Jako warunek zakończenia działania algorytmu przyjmujemy, że błąd jest mniejszy od 1e-20 (czy długości odcinka w którym jest rozwiązanie). Czy algorytm zwraca przybliżenie liczby 2 jedynego pierwiastka tego wielomianu? Narysować wykres tej funkcji liczone z obu wzorów. tzn. (x-2)^3 i x^3-...-8 na odcinkach 2+[-h,2*h] dla h=10^{-3},...,10^{-5} . Czy z obu wykresów wynika, że ten wielomian ma tylko jedno miejsce zerowe w otoczeniu 2 ?
• Lab 4 (gr 1 - 30/03/12; gr 2 - 13/04/12 (przesunięcie zajęć z powodu Świąt Wielkanocnych)) - dokończenie fl;
• Lab 5 (gr 1 - 20/04/12; gr 2 - 27/04/12) - rozwiązywanie układów liniowych: rozkłady LU, LL', uwarunkowanie macierzy, LZNK, zaliczenie 1 projektu
  1. Operator octave'a \ służący m.in. do rozwiązywaniu układów równań liniowych w octavie. Przetestuj ten operator dla macierzy A=[1,2;3,4] i x=[1;3] z f=A*x czy y=A\f jest rozwiązaniem tego układu. Policz normy rezydualną ||Ay-f||_p i normę błędu ||x-x||_p dla p=1,2 . Powtórz testy dla jakiejś macierzy osobliwej? Co się wtedy dzieje?
  2. Funkcja inv(). Zapoznaj się z pomocą do funkcji: inv() . Przetestuj dla A=[1,1;1,-1] czy macierz B obliczona za pomocą tej funkcji rzeczywiście jest macierzą odwrotną? Policz normy pierwszą i Frobeniusa ||A*B-I|| oraz B*A-I|| . Zastosuj funkcje do rozwiązywania układu równań liniowych Ax=f dla znanego x (liczymy f=A*x ). Policz y jako iloczyn B z f i policz błędy rezydualny ||Ay-f|| i ||x-y|| w normie pierwszej i drugiej. Powtórz to zadanie dla macierzy N\times N losowej (funkcja octave'a rand() zwraca macierz losową) dla n=10,50,250,1250 ze znanym rozwiązaniem - oszacuj czas przy pomocy funkcji tic i toc . Porównaj czas i błędy w normie pierwszej dla rozwiązania uzyskanego przy pomocy operator octave'a \ .
  3. Funkcje lu() i chol() . Zapoznaj się z pomocą do funkcji: lu() i chol() . Dla macierzy A=[1,1;1,-1] znajdź jej rozkład PA=LU , rozkład Choleskiego A=L_1^TL_1 . Sprawdź te rozkłady licząc normy macierzowe pierwszą i Frobeniusa błędów np. PA-LU czy A-L_1^TL_1 . Zastosuj te rozkłady do znalezienia rozwiązania układu równań liniowych Ax=f ze znanym rozwiązaniem np. x=[1;1] i f=[2;0] . Policz normy pierwszą i drugą wektorowe pomiędzy x a wynikiem algorytmu w polegającym na zastosowaniu odpowiedniego rozkładu oraz takie same normy rezydualne tzn. normy Aw-f .
  4. \label{zad:mathilb} 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. Policz uwarunkowanie macierzy Hilberta dla powyższych N .
  5. W octavie przetestuj eliminacje Gaussa z częściowym wyborem i bez wyboru dla macierzy A=[e, 1; 1, 1] z e=eps/10 (funkcja octave'a eps zwraca epsilon maszynowy) i wektorem prawej strony f=[1;0]^T . Trzeba zaprogramować eliminację Gaussa bez wyboru dla macierzy 2\times 2 samemu.
  6. Powtórz zadanie~\ref{zad:mathilb} ale w arytmetyce pojedynczej precyzji. Musimy tu trochę sztucznie wymusić aby zmienne były zmiennymi pojedynczej precyzji za pomocą funkcji octave'a single() .
  7. Przetestuj funkcję invhilb() tworzącą macierz odwrotną do macierzy Hilberta (por. zadanie~\ref{zad:mathilb}). Policz normy macierzowe Frobeniusa i pierwsze dla H(N)*iH(N)-I dla dla N=10,20,30 , gdzie iH(N) macierz odwrotna do macierzy Hilberta obliczona przez invhilb() . Przetestuj wykorzystanie iH(N) do rozwiązania układu równań z macierzą Hilberta H(N) , tzn. stwórz macierz H(N) dla N=10,20,30 dla znanego rozwiązania sol , czyli policz H(N)*sol=f - rozwiąż układ H(N)x=f mnożąc F przez iH(N) tzn. x=iH(N)f i policz normy (np. 1,2 ) ||H(N)x-f|| i ||x- sol|| .
  8. Stworzyć w octavie macierz trójdiagonalną A=(a_{i j})_{i,j=1}^n wymiaru n\times n z 2 na diagonali i -1 na pod- i nad diagonalą tzn.
     a_{i, j}= 2 dla i=j
     -1 dla |i-j|=1
     0 dla|i-j|>1
     dla i,j=1,...,n. 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?).
  9. Sprawdzić czy macierz z zadania~\ref{zad:mat3diag} jest symetryczna i dodatnio określona z wykorzystaniem funkcji octave'a chol() .
  10. Zaimplementować rozkład LU dla macierzy trójdiagonalnej bez wyboru elementu głównego, tzn. stworzyć własną funkcję octave'a
      function [x,d,l]=lu3diag(a,b,c,f,N)
      Parametry funkcji:
      1. a,b,c - przekątna, pod-przekątna i nad-przekątna macierzy A ,
      2. f - wektor prawej strony,
      3. N długość przekątnej a .
      Funkcja zwraca x - rozwiązanie Ax=f oraz czynniki rozkładu macierzy A=LU czyli diagonalę macierzy górnotrójkątnej U w wektorze d oraz pod-diagonalę L - macierzy dolnotrójkątnej, trójdiagonalnej z jedynkami na diagonali czyli wektor l . Pamiętajmy, że nad-diagonala U równa się na-diagonali A tzn. wektorowi c .
  11. Zastosować funkcję z poprzedniego zadania do uzyskania rozkładu LU macierzy trójdiagonalnej dla N=4,10,15 . Porównać z wynikiem uzyskanym przy pomocy funkcji octave'a lu().
  12. Zastosować funkcję z poprzedniego zadania (2 do tyłu) do rozwiązania układu równań z macierzą z z 2 na diagonali i -1 na pod i nad diagonalach Porównać czas rozwiązywania Państwa algorytmem, a algorytmem octave'a czyli operatorem A\ f!. Rozwiązać układ dla dowolnego rozwiązania np. 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ą).
  13. Operator octave'a \ służący m.in. do rozwiązywaniu układów równań liniowych w octavie al również LZNK.
      • Przetestuj ten operator dla RLZNK dla macierzy A=[1,1;1,-1;1,3] i x=[1;2] z f=A*x! czy y=A\f jest rozwiązaniem tego układu? Policz normy rezydualną \|Ay-f\|_2 i normę błędu \|x-x\|_2 .
      • Przetestuj ten operator dla nieregularnego LZNK dla macierzy A^T z A=[1,1;1,-1;1,3] i x=[1;1;1] z f=A*x czy y=A\f jest rozwiązaniem tego układu? Policz normy rezydualną \|Ay-f\|_2 i normę błędu \|x-x\|_2 .
  14. Rozkład QR w octave. Funkcje qr().
      1. Zapoznaj się z pomocą do funkcji: qr().
      2. Dla macierzy A=[1,1;1,-1;1,3] znajdź jej rozkład A=QR z pomocą funkcji qr().
      3. Sprawdź czy czy uzyskana Q jest ortogonalna - policz normy Frobeniusa QQ^T-I i Q^TQ-I .
      4. Sprawdź ten rozkład licząc normy macierzowe drugą i Frobeniusa błąd A-QR . Zastosuj ten rozkład do znalezienia rozwiązania LZNK Ax=f ze znanym rozwiązaniem np. x=[1;0] i f=[1;1;1] .
      5. Policz normę drugą wektorową pomiędzy x a wynikiem algorytmu w polegającym na zastosowaniu odpowiedniego rozkładu oraz takie same normy rezydualne tzn. normy drugie Aw-f oraz Rw-Q^Tf .
  15. Układ równań normalnych a rozkład QR czy operator \ . Rozpatrzmy A_{2n,k} pod-macierz 2n\times k macierzy A_{2n,2n} Vandermonde'a dla 2n węzłów równoodległych na [0,1] . Rozwiąż LZNK z A_{2n,k} z wektorem prawej strony równym pierwszej kolumnie tej macierzy (rozwiązanie to pierwszy wersor) - trzema sposobami:
      1. operatorem \ ,
      2. rozkładem QR uzyskanym funkcją qr(),
      3. rozwiązaniem układu równań normalnych: tworzymy macierz układu równań normalnych B=A_{2n,k}^TA_{2n,k} , wektor prawej strony g=A_{2n,k}^T f układu równań normalnych. Rozwiązujemy układ równań normalnych Bx=g operatorem \ .
      Przetestuj dla n=10,20,40,80 i k=2,4,n . Porównaj
      • czas obliczeń
      • błąd - \|x-y\|_2
      • błąd rezydualny \|Ax-f\|_2
      dla x rozwiązania dokładnego LZNK, f wektora prawej strony LZNK, y przybliżenia rozwiązania uzyskanego daną metodą.
  16. Rozkład QR a operator \ przy rozwiązywaniu układów równań liniowych czyli porównanie QR a LU zastosowanych do rozwiązywania układów równań liniowych Rozpatrzmy macierz A_{n,n} Vandermonde'a dla n węzłów równoodległych na [0,1] . Rozwiąż układ równań liniowych z wektorem prawej strony równym pierwszej kolumnie tej macierzy (rozwiązanie to pierwszy wersor) - dwoma sposobami:
      1. operatorem \,
      2. rozkładem QR uzyskanym funkcją qr(),
      Przetestuj dla N=10,20,40,80 . Porównaj
      • czas obliczeń
      • błąd - \|x-y\|_2
      • błąd rezydualny \|Ax-f\|_2
      dla x rozwiązania dokładnego LZNK, f wektora prawej strony LZNK, y przybliżenia rozwiązania uzyskanego daną metodą.
  17. Zastosuj octave'a do rozwiązania zadania znalezienia krzywej zadanej równaniem ax*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 ?
• Lab 6 (gr 1 - 4/05/10; gr 2 - 18/05/10 - znów przesunięcie w związku z Juwenaliami) - LZNK w tym met. Householdera dokończenie i zadanie własne - metoda potęgowa i odwrotna potęgowa, QR
  1. Przetestuj funkcje eig() dla prostej macierzy symetrycznej np [1,3;3,1] i losowej symetrycznej np losujemy A i bierzemy B=A^TA lub AA^T czy A+A^T. Znajdź ich wartości i wektory własne i sprawdź czy rzeczywiście takkimi są tzn policz ||A*Q-D*Q||_1 i ||A-Q*D*Q^T||_1 - Q macierz z kolumnami - wektorami własnymi - D macierz diagonalna z wartościami własnymi na diagonali
  2. zaimplementuj metodę potęgową - jako koniec iteracji proszę przyjąć że ilorazy Raileigha spełniają |r_k-r_{k+1}|<10^{-15}.
  3. Przetestuj metode potęgową dla A=[1,3;3,1] na ekranie drukuj iloraz Railegha i x_{2k}. Sprawdź czy obliczone wartości są rzeczywiście przybliżeniami pary własnej tzn oblicz \|A x_{2k}-r_{2k}x_{2k}\|_1. porównaj z wynikiem eig().
  4. Powtórz poprzednie zadanie dla losowej macierzy symetrycznej oraz oraz macierzy o znanych wartościach własnych np 1,2,3,4,..,11 czy 10,20,30,..,110, lub 1,2,..,10,10+b dla b=0,1,10,100.
  5. Zaimplementuj odwrotną iteracje potęgową tzn. metodę potęgową dla (A-a I)^{-1}. przetestuj ją dla A z poprzedniego zadania z różnymi wartościami a np. a=0 czy a dobre przybliżenie wartości własnej (znanej).
  6. zaimplementuj metodę QR (warunek stopu norma macierzy A-D mniejsza od 10^{-12} dla D diagonali A. Przetetuj dla macierzy z zad 1 oraz macierzy o znanych wartościach własnych np 1,2,3,4,..,11 czy 10,20,30,..,110, lub 1,2,..,10,10+b dla b=0,1,10,100.
  7. Zmodyfikuj metodę QR dodając przesunięcie tzn. A_k - m_k I=Q_kR_k i A_{k+1}=R_kQ_k+m_kI. Za m_k=(A_k)_{n n}. I przetestuj jak wyżej.
• Lab 7 (gr 1 - 25/05/12; gr 2 - 01/06/10) - interpolacja Lagrange'a Funkcja polyval() czyli funkcja obliczająca wartość wielomianu w jednym ub równocześnie wielu punktach czyli wektorowo. Funkcja polyfit() czyli funkcja obliczająca współczynniki wielomianu interpolacyjnego dla zadanych wartości i węzłów. Algorytm Hornera w bazie potęgowej i Newtona. Funkcja polyval() czyli funkcja obliczająca wartość wielomianu w jednym ub równocześnie wielu punktach czyli wektorowo. Funkcja polyfit() czyli funkcja obliczająca współczynniki wielomianu interpolacyjnego dla zadanych wartości i węzłów. Algorytm Hornera w bazie potęgowej i Newtona.
  1. Korzystając z funkcji polyval() narysuj wykres wielomianu x^3+x-2 bez wykorzystania pętli czyli wektorowo.
  2. Test funkcji polyfit(). Wykorzystując funkcję polyfit() znajdź wielomian interpolacyjny L_n F dla funkcji F(x)=\sin(x) dla węzłów -1,0,1 oraz -1,0,1,10. Policz wartości różnicy F-L_n F w węzłach. Oraz narysuj wykresy F i L_nF na jednym rysunku korzystając funkcji plot(). Wartości wielomianu oblicz bez użycia pętli z wykorzystaniem funkcji polyval().
  3. (zadanie o zbieżnosci interpolacji Lagrange'a) Interpolacja Lagrange'a - zbieżność ciągu wielomianów interpolacyjnych dla węzłów równoodległych i węzłów Czebyszewa:
     • Wykorzystując funkcję polyfit() znajdź wielomiany interpolacyjne L_N f dla funkcji f=\sin() dla N węzłów równoodległych na [0,2*pi] dla N=4,8,16,32,64.
     • Oblicz dyskretną normę maksimum różnicy f-L_N f na siatce tysiąc punktowej na tym odcinku tzn. e_N=\max|f(x_k)-L_N f(x_k)|, gdzie x_k pkt siatki. I policz stosunek e_N/e_{2N} dla N<64. Czy błędy maleją do zera? jak zachowuje się e_N/e_{2N}?
     • Narysuj na ekranie wykresy \sin(x) i tych wielomianów dla różnych N - używając funkcji polyval() i plot().
  4. Powtórz (zadanie o zbieżnosci interpolacji Lagrange'a) dla tej samej funkcji i tego samego odcinka ale dla węzłów Czebyszewa. Węzły Czebyszewa to zera T_{n+1}(t)=\cos((n+1)\mathrm{arccos}(t)) na [-1,1] odpowiednio przesunięte i przeskalowane. Czy błędy e_N dla węzłów Czebyszewa są mniejsze niż dla węzłów równoodległych?
  5. Napisz funkcję znajdującą współczynniki w bazie potęgowej wielomianu interpolacyjnego zadanego stopnia dla węzłów równoodległych oraz węzłów Czebyszewa dla danej funkcji, odcinka [a,b]. Tzn. bardziej precyzyjnie napisz funkcję octave'a (w m-pliku): function [LN,eN]=Lagrangeinterp(FCN,a,b,N,type=0)!, która dla zadanego wskaźnika funkcyjnego FCN (ang. \emph{function handle}) do funkcji jednego argumentu: function y=f(x)!, a,b - końców odcinka [a,b], N - stopnia wielomianu interpolacyjnego i typu węzłów 0 - równoodległych, 1 - Czebyszewa. Zwróci wektor LN współczynniki L_N f wielomianu interpolującego funkcję w tych węzłach oraz eN przybliżenie dyskretnej normy maksimum różnicy L_Nf - f na tym odcinku - przybliżenie normy liczymy na dyskretnej siatce zawierającej tysiąc punktów.
  6. Powtórz (zadanie o zbieżnosci interpolacji Lagrange'a) dla obu typów węzłów, tzn. znajdowanie wielomianów interpolacyjnych na węzłach równoodległych i węzłach Czebyszewa, ale dla funkcji f(x)=\log(1+x) na odcinkach [0,1] i [0,10]. Czy dla tej funkcji i obu odcinków błędy w normach maksimum maleją wraz ze wzrostem N? Porównaj wyniki otrzymane w tym zadaniu w obliczeniach z oszacowaniami teoretycznymi błędu interpolacji Lagrange'a.
  7. Interpolacja Lagrange'a przykład Rungego. Powtórz (zadanie o zbieżnosci interpolacji Lagrange'a) dla obu typów węzłów, tzn. znajdowanie wielomianów interpolacyjnych na węzłach równoodległych i węzłach Czebyszewa, ale dla funkcji dla f(x)=1/(1+x*x) na [-5,5]. Czy ciąg wielomianów interpolacyjnych zbiega do f jednostajnie?
  8. Napisz funkcję obliczającą wartość wielomianu zadanego w bazie potęgowej tzn. w(x)=\sum_{k=0}^n a_kx^k algorytmem Hornera. Parametrami funkcji ma być wektor współczynników a i macierz wartości x. Na wyjściu funkcja ma zwrócić wartości wielomianu dla wartości w x. Przetestuj dla 1+x^2 oraz 1-2x+x^2 dla węzłów równoodległych na [-1,2] - narysuj wykres funkcji z wykorzystaniem tej funkcji.
  9. Napisz funkcję znajdującą dla danego wielomianu stopnia n: w(x)=\sum_{k=0}^n a_kx^k oraz liczby q współczynniki wielomianu p(x)=\sum_{k=0}^{n-1}b_kx^k i r takie, że w(x)=(x-q)*p(x)+r, obliczone z wykorzystaniem algorytmu Hornera.
  10. Napisz funkcję znajdującą dla danego wielomianu stopnia n w(x)=\sum_{k=0}^n a_kx^k oraz liczby q wartość w(q) i pochodnej w'(q), obliczone z wykorzystaniem algorytmu Hornera.
  11. Algorytm Hornera w bazie Newtona. Różnice dzielone. Zaprogramuj w octavie zmodyfikowany algorytm Hornera obliczający wartość wielomianu zadanego w bazie Newtona dla danych węzłów. Parametrami maja być x pkt w którym obliczamy wielomian (ewentualnie tablica punktów ale wtedy funkcja też musi zwrócić tablice wartości wielomianu w tych punktach), N - stopnień wielomianu, wektor długości N+1 ze współczynnikami wielomianu w bazie Newtona. Przetestuj na kilku b prostych przykładach np. 3 węzły -1,0,1 i wielomian w(x)=x^2 - w bazie Newtona związanej z tymi węzłami: x^2=(x+1)x - (x+1) +1 . Zastosuj do przykładów z poprzednich zadań - dla różnych N, narysuj wykresy funkcji i jej wielomianu interpolacyjnego otrzymanego za pomocą tego algorytmu (wartość wielomianu liczymy algorytmem Hornera).
  12. Napisz funkcję która dla danego wielomianu, którego współczynniki w bazie potęgowej znamy oblicza współczynniki tego wielomianu w bazie Newtona dla zadanych węzłów y, tzn. znamy wektor a=(a_k)_k taki, że w(x)=\sum_{k=0}^n a_kx^k i współczynniki bazy Newtona y_k dla k\geq 0. Chcemy znaleźć współczynniki b_k takie, że w(x)=\sum_{k=0}^n b_k w_k, \qquad w_k(x)=\Pi_{j=0}^{k-1}(x-y_j). Założenie jest takie, że znamy przynajmniej n-1 współczynników y_k, w przeciwnym razie funkcja powinna przerwać działanie sygnalizując błąd na ekranie.
  13. Napisz funkcję która dla danego wielomianu, którego współczynniki w bazie Newtona w_k(x)=\Pi_{j=0}^{k-1}(x-y_j) dla k=0,...,n znamy oblicza współczynniki tego wielomianu w bazie potęgowej 1,x,x^2,...,x^n. Węzły bazy Newtona też muszą być parametrem funkcji. Tzn. znamy wektor b=(b_k)_k taki, że \sum_{k=0}^n b_k w_k i współczynniki bazy Newtona y_k dla k\geq 0. Chcemy znaleźć współczynniki a_k takie, że
      w(x)=\sum_{k=0}^n a_k x^k.
• Lab 8 (gr 1 - 8/06/10) - interpolacja splajnowa, zaliczenie projektu 2
  1. Funkcje octave'a spline() and ppval. Zapoznaj się z tymi funkcjami ( help spline i help ppval). Wykorzystując te funkcje narysuj wykres splajnu kubicznego s_1 na podziale równomiernym odcinka [-3,3] z węzłami \{x_k=k\} dla k=-3,-2,...,3 - ręcznie podając losowe wartości przyjmowane przez te splajny w tych węzłach: np.: s_1(x_k)=(-1)^k , czy odpowiednio s_1(x_k)=1 . Następnie narysuj wykres splajnu s_2 o tych samych wartościach w węzłach ale podając dodatkowo wartości pochodnych w końcowych węzłach równe zero, tzn. wywołaj funkcje spline() podając dwie wartości więcej. Czy otrzymaliśmy te same splajny? Policz przybliżoną normę różnicy s_1-s_2 na odcinku [-3,3] .
  2. Splajn bazowy interpolacyjny - Dla danych węzłów równoodległych \{k\}_{k=-5,-4,...,5} na [-5,5] narysuj wykres splajnu kubicznego takiego, że s(0)=1 i s(k)=0 dla pozostałych węzłów k\not =0 oraz któ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 zadanei ale dla splajnu typu not-a-knot tzn. nie podając wartości pochodnych w końcach dla funkcji spline().
  3. \label{zad:Bsplajnkub} Splajn kubiczny o minimalnym nośniku. Dla danych węzłów równoodległych \{k\}_{k=-5,-4,...,5} na [-5,5] narysuj wykres splajnu kubicznego takiego, że s(-1)=s(1)=1,s(0)=4 i s(k)=0 dla węzłów k\not \in\{-1,0,1\} oraz ma pochodne równego zero w węzłach skrajnych tzn. : -5 i 5 . Czy poza [-2,2] ten splajn jest równy zero? Policz przybliżone normy maksimum na [-5,-2] i [2,5] dla tego splajnu i narysuj jego wykres.
  4. Testowanie eksperymentalne rzędu zbieżności splajnu interpolacyjnego kubicznego z hermitowskimi warunkami brzegowymi. Korzystając z funkcji octave'a spline() znajdź współczynniki splajnu interpolacyjnego kubicznego S_N na N węzłach równoodległych dla funkcji f(x)=\sin(x) na odcinku [-\pi,2*\pi] dla N=2^kN_0 dla N_0=5 i k=1,2,3,4,5 . Następnie
     • narysuj wykresy f(x) i tych splajnów dla różnych N .
     • oblicz dyskretną normę maksimum na siatce równomiernej 1000 punktowej na tym odcinku tzn. e_N=\max|\sin(x_k)-S_N \sin(x_k)| dla x_k punktów siatki.
     • Policz równocześnie współczynnik \frac{e_{N}}{e_{2*N}} . Czy widać, że \frac{e_{N}}{e_{2\,N}}\approx 2^p dla jakiegoś p całkowitego np. p=4,8 lub 16 ?
  5. Testowanie eksperymentalne rzędu zbieżności splajnu interpolacyjnego kubicznego bez warunków brzegowych (splajn typu \emph{ang. not a knot}). Powtórz poprzednie zadanie ~\ref{zad:spline-cub-int} ale dla splajnów interpolacyjnych otrzymanych przez spline() bez podawania wartości pochodnych w skrajnych węzłach. Czy współczynniki \frac{e_{N}}{e_{2*N}} są te same? Tzn. czy szybkość zbieżności \|\sin(x)-S_N\|_\infty jest taka sama?
  6. Testowanie eksperymentalne rzędu zbieżności splajnu interpolacyjnego kubicznego naturalnego (warunek brzegowy - zerowanie się drugich pochodnych w końcach odcinku). Powtórz zadanie poprzednie ale dla splajnów interpolacyjnych naturalnych. Tu trzeba wykorzystać funkcję z octave-forge (czyli rozszerzenia pakietu octave)\\ pp=csape(x,y,'variational')! \\- ostatni argument określa to, że splajn będzie naturalny. Tu link do pomocy do funkcji octave'a csape()
  7. Przykład Rungego czyli f(x)=1/(1+x*x) i odcinek [-5,5] a zbieżność interpolacji splajnami kubicznymi. Przetestuj jak w poprzednich zadaniach czy splajny interpolacyjne kubiczne z podanymi warunkami na pochodne w końcach odcinka zbiegają w normie supremum do f , tzn. korzystając z funkcji octave'a spline()! znajdź współczynniki splajnu interpolacyjnego kubicznego S_N na N węzłach równoodległych dla f na odcinku [-5,5] dla N=2^kN_0 dla N_0=5 i k=1,2,3,4,5 oraz narysuj wykresy f i tych splajnów dla różnych N . Następnie oblicz dyskretną normę max na siatce złożonej z tysiąca punktów na tym odcinku tzn. e_N=\max|\sin(x_k)-S_N \sin(x_k)| dla x_k=-5+k*0.01 z k=0,...,1000. Policz równocześnie współczynnik \frac{e_{N}}{e_{2*N}} . Czy \frac{e_{N}}{e_{2*N}}\approx 2^p dla jakiegoś p całkowitego?

Kilka przykładowych skryptów, m-plików (plików funkcyjnych ) octave'a, czy użytecznych linków
(dodawanych w miarę postępu labu)

Skrypty m-pliki octave'a

Literatura:

• [KC2006] D.Kincaid, W.Cheney, Analiza numeryczna, WNT, 2006
• [Kic2012] Przemysław Kiciak, Metody numeryczne dla informatyków plik pdf, skrypt, 2012
• [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) plik pdf.

Ostatnia aktualizacja: 14 wrzesnia 2012