g Zajęcia semestr letni 2012/13

Semestr letni 2014/15

Lab do RRzw (z labem)
Matematyka Obliczeniowa
- Ćwiczenia/Lab
- Program bieżącego labu
- Zadania domowe, projekty na lab, stare zadania z kolokwium
- Wyniki zadań domowych, projektów z labu, kolokwium etc: plik-pdf

Zasady zaliczenia przedmiotu

Wg ustaleń wykładowcy tzn. zaliczenie ćwiczeń/labu od 50% możliwych do uzyskania punktów - na bazie wspólnego kolokwium i ewentualnie kartkówek, projektów z labu etc

Ćwiczenia/Lab MO

Lab piątek 1015-1145, sala 2042 (co 2 tygodnie wg planu USOSie)

Gr 1 ćwiczenia piątek 830-10, sala 5070. Pierwszy lab wg USOSa - drugi wtorek zajęć tzn. 6 marca i potem co 2 tygodnie
Gr 5 ćwiczenia - piątek 1215-1345, sala 5070. Pierwszy lab wg USOSa- pierwszy wtorek zajęć 27 lutego i potem co 2 tygodnie

Uwaga Kolokwium środa 29 kwietnia 2015 (w terminie i miejscu wykładu),
Druga kartkówka ostanie zajęcia (w pt w czasie ćwiczeń) - ok. 20 min - 2 zadania, zakres od zadan z LZNK do zadań z interpolacji wielomianowej włącznie tzn. 3 seria (zad 7,9-11) i 4 seria (zad 1-4) - bez zadań oznaczonych jako trudne.

Wyniki kolokwium/kartkówki/labów: plik-pdf (proszę sprawdzać czy wszystko się zgadza i ewentualnie reklamować)

Zaliczenie ćwiczeń i labu

Ustali wykładowca

Serie zadań domowych i projekty z labu, stare zadania etc

Serie ZD zadawane co ok. kilka zajęć.
Druga kartkówka ostanie zajęcia (w pt w czasie ćwiczeń) - ok. 20 min - 2 zadania, zakres od zadan z LZNK do zadań z interpolacji wielomianowej włącznie tzn. 3 seria (zad 7,9-11) i 4 seria (zad 1-4) - bez zadań oznaczonych jako trudne.

Pierwsza seria ZD równania nieliniowe
Druga seria ZD fl, normy, LU.
Trzecia seria ZD LU cd, Choleski, uwarunkowanie, Householder, LZNK, metoda potegowa (ostatnie 2 zadania - nie wejda w zakres pierwszej kartkówki)
Projekt z labu (prosta najbliższa danym punktom jako LZNK rozwiazane rozkladem QR za pomocą algorytmu ortogonalizacji Gramma-Schmidta) na przedostatni lab czyli koniec maja/poczatek czerwca
Czwarta seria ZD - interpolacja Lagrange'a, interpolacja Hermite'a, splajny liniowe i kubiczne
Piąta seria ZD - aproksymacje średniokwadratowa (w przestrzeniach z iloczynem skalarnym), jednostajna, kwadratury

PROSZĘ NIE przysyłać rozwiązań mailem
Wyniki kolokwium/kartkówki/labów) plik-pdf

Zadania z egzaminów z poprzednich lat

Zadania z kolokwiów (2gr i poprawkowe) i egzaminu w I terminie z roku 2011/12

Zadania z kolokwium(2 grupy) z 2012/13

Program labów

Z ewentualnymi linkami do jakiś prostych skryptów. Będę starał się tu umieszczać krótkie opisy tego co było czy ma być na labie

Lab 1 (gr 5 - 27/02/15; gr 1 - 06/03/15) - 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 @).
mobasic.m - skrypt z przykładami powyższych operacji czyli tym o czym jest ten lab.
1. - Utwórz dowolne macierze 3x4 A i 3x5 B - a następnie macierz 3x9 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 3x3 tzn. 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(x) 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 - równania nieliniowe
1. Przetestuj metodę bisekcji - dla rozwiązania równania cos(x)=0 na odcinku [0,2] - sprawdź jak działa. Mozna samemy zaimplementowac -albo uzyc skryptu testbis.m
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. Znajdź w pomocy funkcję octave'a - rozwiązująca równania nieliniowe - i przetestuj ją na kilku prostych przykładach skalarnych np.
  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.
4. Zaimplementować metodę Newtona w octavie i przetestować jej zbieżność dla następujących funkcji:
  - 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 .
5. 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 .
6. 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.
7. Porównać wyniki z zadania poprzedniego z funkcją octave'a fsolve() (pomoc: help fsolve ).
8. 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ń.
9. 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?
10. 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?)
11. 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
Jak Panstwo nie zdazyli sami napisac to prosze uzyc skryptu testnewton.m w ktorym sa implementacje metody newtona
Lab 3 - 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, 1+2eps 1+eps/2 obliczone w arytmetyce podwójnej precyzji w octavie są 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 ?
Jak Panstwo nie dali rady sami napisac to prosze uzyc skryptu testfl.m z rozwiazaniami niektorych zadan
Lab 4 rozwiązywanie układów liniowych: rozkłady LU, LL', uwarunkowanie macierzy.
1. Normy macierzowe - policz normy 1,2,Frobeniusa, max dla macierzy diagonalnej i konkretnej niesymetrycznej A=[1;2;1;1].
2. 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 ||y-x||_p dla p=1,2,\infty. Powtórz testy dla jakiejś macierzy osobliwej? Co się wtedy dzieje?
3. 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 \ .
4. 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 .
5. Przetestuj solver octave'a dla macierzy górnotrójkątnych U1 i U2=-U1+2*I dla U1 z jedynkami na diagonali i jedynkami na wszystkich pozycjach nad diagonala (polecenie octave'a triu(ones(n,n)) ją tworzy), tzn. stwórz macierze Uk (k=1,2) dla N=10,20,30,50,100 dla znanego rozwiązania np. losowego czyli sol_k=1 policz Uk*sol=fk - rozwiąż w octavie układ z Uk i fk i policz normy (2 itd) ||Uk*xk-fk|| i ||xk- sol|| dla xk rozwiązania które wyliczył octave (k=1,2).
6. 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 .
7. Powtórz zadanie z macierzą Hilberta ale w arytmetyce pojedynczej precyzji. Musimy tu trochę sztucznie wymusić aby zmienne były zmiennymi pojedynczej precyzji za pomocą funkcji octave'a single() .
8. Przetestuj funkcję invhilb() tworzącą macierz odwrotną do macierzy Hilberta (por. zadanie z macierzą Hilberta). 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|| .
9. 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.
10. 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?).
11. Sprawdzić czy macierz 3-diagonalna z 2 na diagonali i -1 na pod i nad diagonali jest symetryczna i dodatnio określona z wykorzystaniem funkcji octave'a chol() .
12. 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 .
13. 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().
14. 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ą).
Jak Panstwo nie dali rady sami napisac to prosze uzyc skryptu uklrlin.m z rozwiazaniami niektorych zadan
Lab 5 - w tym roku mało labów więc zadań b dużo zrobimy tylko niektóre : rozkład QR, regularne LZNK, przekształcenie Householdera, SVD zadanie wlasne - zaliczenie projektu....
1. 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 ||y-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 ||y-x||_2 .
2. Rozkład QR w octave. Funkcje qr(). (o ile nie bylo)
  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 .
3. Dla macierzy z poprzednich zadań rozwišż LZNK - rozwiązując układ r. normalnych - porównaj wyniki.
4. 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ówno-odległ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ą.
5. 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ówno-odległ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ą.
6. 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 ?
7. Zaprogramuj funkcje octave'a
  function y=H(x,w,nw)
  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. Przetestować dla losowych wektorów \vec{x} i \vec{w}, sprawdzić czy ||H_w \vec{x} ||_2=||\vec{x} ||_2. Można napisać ogólniejsza wersje gdzie x=X macierz N x 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 nasza funkcje.
8. (testujemy twierdzenie z wykładu) Znajdź wektor Householdera \vec{w} taki ze 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 tak jest. Tzn czy H_w \vec{u}= l \vec{v}.
9. (SVD i nieregularne LZNK)
  przeczytaj help do funkcji octave'a svd - zastosuj ja do macierzy A1=[1 2 3;1 -1 0], A2=[1 1; 1 2; -3 4] i A3=[1 1 2; 1 0 2;2 1 4] (ta macierz nie ma max rzedu) - znajd dla nich rozkład SVD i za jego pomocš rozwišż LZNK z losowymi wektorami prawej strony (f=randn(3,1)). Znajdz ortonormalne bazy jadra i odpowiednio obrazu tzn R(Ak) tych macierzy. Okresl normy drugie tych macierzy na bazie wyniku svd. Porownaj z wynikami otrzymanymi operatorem backslash \, funkcja norm(Ak,2) itd - sprawdz czy rzeczywiscie otrzymane SVD to dobre rozklady tzn czy \|A-QSV'\|_1==0?
10. (SVD dla macierzy rzędu jeden) Dla macierzy A=u*v' u wektor dlugosci m np losowy m=5 i v wektor losowy dlugosci np n=3. Policz svd za pomoca funkcji svd(). Znajdz wektory Householdera w_k takie że H_1AH_2=S - macierz z warosciami szczegolnymi z rozkładu svd. Sprawdz czy \|S-H_1AH_2\|_1==0?
11. 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 takimi 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
12. zaimplementuj metodę potęgową - jako koniec iteracji proszę przyjąć że ilorazy Raileigha spełniają |r_k-r_{k+1}|<10^{-15}.
13. Przetestuj metodę 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().
14. 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.
15. 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).
16. zaimplementuj metodę QR (warunek stopu norma macierzy A-D mniejsza od 10^{-12} dla D diagonali A. Przetestuj 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.
17. zaimplementuj metodę równoczesnych iteracji tzn dla danej Z_k ortogonalnej (startujmy z Z_0=I) liczymy Y_{k+1}=AZ_k i znajdujemy rozkład QR macierzy Y_{k+1}=Z_{k+1}R_{k+1} tzn Z_{k+1} ortogonalna R_{k+1} górnotrójkątna (warunek stopu norma macierzy Z_k^TAZ_k - diagonala(Z_k^TAZ_k) mniejsza od 10^{-12} dla D diagonali Z_k^TAZ_k.) Przetestuj 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.
Lab 6
Interpolacja Lagrange'a i splajnami kubicznymi.
Funkcja polyval() czyli funkcja obliczająca wartość wielomianu w jednym lub 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. Funkcje ppval() i spline() działające analogicznie ale dla splajnów kubicznych.
1. 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()
2. 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?

Zadania, których nie zdążą Państwo zrobić na labie, są do domu

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

Link do stron Octave'a (skąd można ściągnąć kolejną dystrybucje - pod linuxa czy windows)
octave-forge - rozszerzenia octave'a
Manual do octave'a w htmlu (po angielsku)
Link do skryptu do labu z mo z octave'em plik pdf

Skrypty m-pliki octave'a

mobasic.m - prosty skrypt octave'a w którym są przedstawione niektóre podstawowe operacje macierzowe, oraz pętle, instrukcje warunkowe, zapisywanie/wczytywanie do/z plików itp
testbis.m - skrypt octave'a z funckja implementujaca metode bisekcji i kilkoma testami
testnewton.m skrypt w ktorym sa implementacje metody Newtona (z dokladna i przyblizona pochodna) i metody siecznych - mozna testowac rzad metod o ile znamy dokladne rozwiazanie

Skrypty będą się pojawiały stopniowo.
W razie znalezienia błędów proszę o kontakt (częśc błędów może się brać ze zmian w octave)

Literatura:

Podręczniki:

[KC2006] D.Kincaid, W.Cheney, Analiza numeryczna, WNT, 2006
[Kic2014] Przemysław Kiciak, Metody numeryczne dla informatyków /matematyka obliczeniowa plik pdf, skrypt, 2014/15
[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.

Pozycje [Kic2014], [Mos2002] i [Pla2002] to skrypty dostępne dla studentów naszego wydziału. A pozycja [FMW2005] to książka skierowana raczej do studentów politechniki ale większość algorytmów jest w niej opisana. [KC2006] jest podstawowym podręcznikiem.
Wykład jest wg [Kic2014].

Inne użyteczne linki

Literatura dodatkowa dla osób zainteresowanych metodami numerycznymi, obejmująca materiał częściowo lub często całkowicie poza zakresem wykładu
Eseje wyjaśniające mam nadzieję czym jest i czym na pewno nie jest Analiza Numeryczna (czy inaczej Metody Numeryczne)

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.

Inne eseje tegoż autora o analizie numerycznej i nie tylko http://people.maths.ox.ac.uk/trefethen/essays.html

A tu link do wykładu z Metod Numerycznych on-line na ważniaku: wykłady i ćwiczenia

Powrót do mojej strony domowej.

Ostatnia aktualizacja: 30 kwietnia 2015