Semestr letni 2015/16

• Obliczenia naukowe
• Matematyka obliczeniowa dla BIOINF (wyklad): wyniki kolokwium, projektu i egzaminu: plik-pdf
• Matematyka Obliczeniowa (przedmiot kursowy; 2 rok matematyki)
  • Ć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

Ćwiczenia/Lab MO

• Gr 1 ćwiczenia piątek 830-10, sala 5070. Lab pt 1015-1145 w 2042. Pierwszy lab wg USOSa i potem co 2 tygodnie
• Gr 2 Lab pt 1015-1145 w 2042. Pierwszy lab wg USOSa i potem co 2 tygodnie wg USOSa
• Gr 3 Lab pt 1215-1345 w 2042. Pierwszy lab wg USOSa i potem co 2 tygodnie wg USOSa

Zaliczenie ćwiczeń i labu

Zadania z egzaminów z poprzednich lat

Program labów

• Lab 1 - 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.
  Kilka zadań (aby nabrać wprawy):
  1. Operacje macierzowe
     • 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 .
     testnewton.m - prosty skrypt umozliwiajacy testy zbieznosci metody Newtona
  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ą. proszę wpisać z linii komand: format long e;d=1.234567890123456
  #a potem
  s=single(d)
  whos
  #czy wszystkie cyfry sie zachowaly?
  d-double(s)
  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 i większość zadań jest po to by zauważyć tą własność.
  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 największe n takie że 2^(-n) +1>1 ale w pojedynczej precyzji.
  3. Narysuj wykres funkcji f(x)=(1+a)-a na [0,1] dla różnych a=5^k dla k=1,2,..,20,..30 . 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=-1/(x+a) dla x=10^k i k=4,...,10 . Czy widać różnicę w wyniku? Policz różnicę bezwględną tzn ab=|x1-x2| i względny wb=1.0./|x2| (dlaczego spodziewamy się że x2 lepiej obliczone?) 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. (coś podobnego jest już zrobione w skrypcie z poprzedniego labu: testbis.m ) Zastosować algorytm bisekcji dla funkcji (x-r)^3 liczonej z wzoru na rozwinięcie dwumianu x^3-...-r^3 startując z a=r-1e-3 a b=r+1e-3 - r losowa liczba z [1,2]: r=1+rand. 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 r jedynego pierwiastka tego wielomianu? Narysować wykres tej funkcji liczone z obu wzorów. tzn. (x-r)^3 i x^3-...-r^3 na odcinkach r+[-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 r ?
  
  Jak Państwo nie dali rady sami napisać to proszę użyć skryptu fl16.m z rozwiązaniami niektórych zadań
• 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 . Porównaj czas rozwiazywania np tic;[l,u,p]=lu(A);x=u\(l\(p*f));toc, porównaj z \ oraz z chol().
  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ą).
  W skrypcie rozwiazania zadan ktore udalo sie zrobic na labie: lu16.m
• Lab 5 - rozkład QR, regularne LZNK, przekształcenie Householdera
  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. (to bylo na poprzednim labie dla innej macierzy - do domu ) 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. Wykonaj metodz Householdera rozklad QR macierzy tzn.: A=[1,1;1,-1;1,3] i rozwiaz LZNK uzywajac tego rozkladu z f=[1;1;3] - porównaj z rozwiazaniem otrzymanym backslashem.
     Powtórz dla losowej macierzy Nx2.
  10. Porównaj czas obliczania H*A dla macierzy A mxn (m>=n)i macierzy Householdera (z losowym wektorem Householdera) mxm dla m=10,10,..,10^k i m=1,10,n na 2 matematycznie równowazne sposoby - dobry: HA1=A- x*((2/(x'*x))*(x'*A)) i drugi zly: (dlaczego?) HA2=A- (2/(x'*x))*(x*x')*A. Wyjasnienie - policz koszt obliczeniowy poszczegolnych nawiasow i ilosc wymaganej pamieci (koszt pamieciowy).
  W skrypcie rozwiazania zadan ktore udalo sie zrobic na labie: lznk16.m
• Lab 6 SVD, nieregularne LZNK, zadanie własne - metoda potęgowa i odwrotna potęgowa, QR, zaliczenie projektu 2 czerwca 2016
  1. (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) - znajdz dla nich rozkład SVD i za jego pomoca rozwiaz LZNK z losowymi wektorami prawej strony (f=randn(3,1)). Znajdz ortonormalne bazy jadra i odpowiednio obrazu tzn R(A) 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?
  2. (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?
  3. 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
  4. zaimplementuj metodę potęgową - jako koniec iteracji proszę przyjąć że ilorazy Raileigha spełniają |r_k-r_{k+1}|<10^{-15} lub ze iterujemy okreslona ilosc razy np 30 lub 50.
  5. 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().
  6. 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.
  7. 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).
  8. 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.
  9. 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.
  svd16.m - kilka rozwiazan zadan na svd, i zadanie wlasne
Lab 7
Interpolacja Lagrange'a i splajnami kubicznymi. Zaliczenie projektu dla gr z labem 27 V 2016
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. 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_n F na jednym rysunku korzystając funkcji plot(). Wartości wielomianu oblicz bez użycia pętli z wykorzystaniem funkcji polyval().
3. Test funkcji polyfit(). Wykorzystując funkcję polyfit() znajdź wielomian bazowy lagrange dla 2 wezla od lewej dla węzłów -1,2,3,4. narysuj wykres. Wartości wielomianu oblicz bez użycia pętli z wykorzystaniem funkcji polyval().
4. Powtorz zadanie ale dla splajnu kubicznego dla tych samych wezlow i tych samych wartosci interpolacyjnych. Uzyj spline() i ppval().
5. (zadanie o zbieżności interpolacji Lagrange'a) Interpolacja Lagrange'a - zbieżność ciągu wielomianów interpolacyjnych dla węzłów równo-odległ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ówno-odległ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=4,..,32. 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().
6. Powtórz (zadanie o zbieżności 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 tzn jesli (t_k) zera T_{n+1} na [-1,1] to (x_k=(a+b)/2+0.5*(b-a)*t_k) to wezly Czebyszewa na [a,b] . Czy błędy e_N dla węzłów Czebyszewa są mniejsze niż dla węzłów równo-odległych?
7. Napisz funkcję znajdującą współczynniki w bazie potęgowej wielomianu interpolacyjnego zadanego stopnia dla węzłów równo-odległ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. 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ówno-odległ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.
8. Powtórz (zadanie o zbieżności interpolacji Lagrange'a) dla obu typów węzłów, tzn. znajdowanie wielomianów interpolacyjnych na węzłach równo-odległ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.
9. Interpolacja Lagrange'a przykład Rungego. Powtórz (zadanie o zbieżności interpolacji Lagrange'a) dla obu typów węzłów, tzn. znajdowanie wielomianów interpolacyjnych na węzłach równo-odległ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?
10. Znajdz spline kubiczny taki ze we wszystkich 10 wezlach rownooodleglych oprocz 4 rowny zero a w 4tym jest rowny jeden. (no-a-knot lub hermitowski z zerowymi pochodymi w koncach). narysuj wykresy.
11. Splajn o min nosniku. Dla wezlow orownoodleglych x=1,..,11 tj x(k)=k znajdz splajn kubiczny s typu not-a-knot taki ze s(x(5))=1=s(x(7)), s((6))=4, i s(x(k))=0 dla k róznych od 5,6,7. narysuj wykres. Okresl nosnik splajnu. (Czy jest rowny [4,8]?) Powtorz zadanie dla splajnu Hermitowskiego z s'(a)=s'(b)=0.
    rozdziel16.m - prosta implementacja algorytmu róznic dzielonych, Newtona (znajdowania roznic dzielonych) i uogólnionego algorytmu Hornera dla wielomain uw bazie Newtona
    interp16.m - kilka rozwiazan zadan - interpolacja lagrange'a i splajnowa
    testinterp16.m - kilka rozwiazan zadan z labu z 10 czerwca 2016 - interpolacja lagrange'a i splajnowa

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 funkcja implementujaca metode bisekcji i kilkoma testami
  
• testnewton.m skrypt w ktorym sa implementacje metody Newtona - mozna testowac rzad metody o ile znamy dokladne rozwiazanie
  
• fl16.m z rozwiązaniami niektórych zadań na fl
• lu16.m - rozwiazania zadan na uwarunkowanei i rozklad lu, qr dla rozwiazywania ukladow r. liniowych

Literatura:

• [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.

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: 8 kwietnia 2016