Matematyka obliczeniowa 2018/19

• Bieżący lab
• Zadania domowe, stare zadania z kolokwium
• Wyniki zadań domowych, kolokwium etc: po 1 kartkówce, kolokwium. plik-pdf

Ćwiczenia MO

• ćwiczenia piątek gr1 830-10 sala 5070 - lab gr1 co 2 tygodnie piątek 1015-1145 sala 2042 (wg USOSa)

Zaliczenie ćwiczeń

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 - 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. 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 ?
     Metoda bisekcji dla f(x^*)=0 z f ciągłą : startujemy z odcinka [a,b] takiego, że f(a) i f(b) mają różne znaki - liczymy x=(a+b)/2 i f(x) - patrzymy na znak jeśli taki jak f(a) zamieniamy x z a wpp x z b a jeśli f(x)=0 kończymy iteracje... i powtarzamy poprzednią procedurę dla nowego odcinka. Otrzymujemy zstępujący ciąg przedziałów dla których ciąg lewych/prawych końców czy środków zbiegają do zera funkcji f. Biorąc n-ty środek jako n-te przybliżenie xn mamy |x^*-xn| < = 2^(-(n+1))|b-a| dla a,b końców wyjściowego odcinka. Oczywiście jak funkcja ma kilka pierwiastków to mamy gwarancję że proces zbiegnie do jakiegoś pierwiastka.
  
  fl19.m - prosty skrypt z niektórymi rozwiązaniami.
• Lab 3 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. Dla macierzy A MXM takiej, że ostatnia kolumna i główna przekątna A mają jedynki, wszystkie elementy ściśle pod przekątną są równe minus jeden a pozostałe elementy są równe zero tzn. A(k,k)=A(j,M)=1; A(j,k)=-1 dla j>k; (jak utworzyć taką macierz w 1 linijce kodu bez użycia pętli?) policz dla znanego wektora r: f=A*r rozwiąż backslashem układ Ax=f - policz residuum Ax-f oraz błąd x-r i ich normy dla M=20,25,..,50. Powtórz zadanie w pojedynczej precyzji dla M=10,15,..,40.
  10. 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.
  11. 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?).
  12. 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() .
  13. 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 .
  14. 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().
  15. 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 rozwiązania zadań które udalo sie zrobić na labie: lu19.m
• Lab 4 - dokończenie zadań na układy równań liniowych, rozkład QR, regularne LZNK, przekształcenie Householdera
  1. Zaliczenie projektu.
  2. Zrobić zadanie z poprzedniego labu: 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 .
  3. 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 .
  4. 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 .
  5. Dla macierzy z poprzednich zadań rozwiąż LZNK - rozwiązując układ r. normalnych - porównaj wyniki.
  6. 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ą.
  7. (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ą.
  8. 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 ?
  9. 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.
  10. (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}.
  11. 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.
  12. 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:
• Lab 5
  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. 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. (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()i L_N g dla g=sin(4x) dla N węzłów równo-odległych na [0,2*pi] dla N=4,8,12 etc.
     • Oblicz dyskretną normę maksimum różnicy f-L_N f i analogicznie g-L_N g 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().
  5. 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?
  6. Zapoznaj się z funkcja spline() i ppval(). Znajdź splajn kubiczny typu not-a-knot przyjmujący w węzłach: (-1,0,1,2,3) następujące wartości (1,-1,4,6,-7). narysuj wykres splajnu. Powtórz zadanie ale dodatkow dla splajnu mającego zerowe pochdne w końćach odcinku.
  7. narysuj wykres splajnu na węzłach (0,1,...,10) równoodległych takiego, który przyjmuje wartość zero w każdym węźle oprócz węzła 6 gdzie przyjmuje wartość jeden.
  8. narysuj wykres splajnu na węzłach (0,1,...,10) równoodległych takiego, który przyjmuje wartość zero w każdym węźle oprócz 3 kolejnych węzłów (3,4,5) w których przyjmuje wartości odpowiednio (1,4,1). Jaki jest nośnik tego splajnu.
  9. Znajdź s(N) splajn interpolacyjny dla N węzłów równoodległych na [-5,5] dla f=1/(1+x^2) policz błąd e(N)=||f-s(N)|| w dyskretnej normie supremum (na 10^4 pktów). Przetestuj dla N=10,20,40,80. Pokaż na ekranie e(N) oraz e(N)/e(2N). Przetestuja dla splajnu typu not-a-know.
     Powtórz badanie biorąc N0-2 losowych węzłów z (-5,5] i dodając końce odcinków tzn -5< x(1)< .. < x(N-2)< 5. następnie zagęszczaj węzły dokładając pkty (x(k)+x(k+1))/2 otrzymując kolejno 2N0, 4N0 etc węzłów.
     interp19.m kilka rozwiazan
• Lab 6 kwadratury, aproksymacja średniokwdratowa.
  1. Zapoznaj się z pomocą do funkcji octave'a quad(). Policz przy pomocy quad() całkę \int_a^b f(t) dt dla następujących funkcji i odcinków:
     • x^2 na [-1,1]
     • x^9 na [0,1]
     • sin(x) na [0,\pi]
     • cos(100*x) na [0,\pi]
     • cos(100*x)*cos(1000*x) na [0,\pi]
     Czy wyniki są zgodne z teorią?
  2. Złożona kwadratura trapezów.
     • Zaprogramuj w octavie funkcję function c=trap(FCN,a,b,n) obliczającą złożoną kwadraturą trapezów:
       T_n f=h\left(0.5[f(a)+f(b)]+\sum_{k=1}^{n-1} f(a+k*h)\right) h=(b-a)/n
       Parametry funkcji:
       • FCN - wskaźnik do funkcji octave'a function y=f(x) obliczającej wartość funkcji podcałkowej
       • a - lewy koniec odcinka
       • b - prawy koniec odcinka
       • n - ilość obliczeń wartości funkcji podcałkowej w kwadraturze trapezów minus jeden (tak jak we wzorze powyżej)
       Funkcja zwraca przybliżoną wartość całki obliczoną za pomocą powyższego wzoru tzn. złożonej kwadratury trapezów. Funkcja powinna działać również jeśli ją wywołamy tylko z trzema parametrami. Wtedy n powinno domyślnie przyjąć wartość dwieście.
     • Przetestuj dla funkcji f, której całkę znamy dla N_k=2^k*N_0 k=0,1,2,.... z ustalonym N_0=5 licząc:
       E_N/E_{2N}
       dla błędu
       E_N=|\int_a^b f dt - T_N f|
       dla następujących funkcji
       • x^2 na [0,1]
       • sin(x) na [0,pi] i N_0=5 czyli funkcji analitycznej
       • sin(100*x) na [0,pi] i N_0=5 czyli funkcji analitycznej ale silnie oscylującej (duże wartości drugiej pochodnej)
       • x^{j+0.5} na [0,1] dla j=0,1,2 czyli funkcji w C^j([0,1]) i o eksplodującej w zerze j+1 pochodnej.
       Porównaj wyniki obliczane kwadraturą trapezów z wynikami funkcji quad(). (Można sprawdzić dla jakiego n błąd obliczony złożoną kwadraturą trapezów jest na poziomie błędu funkcji octave quad() - i porównać wtedy ilość wywołań funkcji f przez obie procedury).
  3. Złożona kwadratura prostokątów.
     • Zaprogramuj w octave funkcję function c=kwadprost(FCN,,a,b,n) ze złożoną kwadraturę prostokątów
       P_n f=h \sum_{k=1}^n f(a+ (k-0.5)*h)) h=(b-a)/n.
       Parametry funkcji: wskaźnik do funkcji obliczającej FCN, końce odcinków a,b i ilość punktów n. Funkcja zwraca przybliżoną wartość całki obliczoną za pomocą powyższego wzoru tzn. złożonej kwadratury prostokątów. Funkcja powinna działać również jeśli ją wywołamy bez ostatniego parametru. Wtedy n powinno domyślnie przyjąć wartość np. sto lub dwieście.
     • Przetestuj tak samo jak kwadraturę trapezów tzn. dla tych samych funkcji f, której całkę znamy dla N_k=2^kN_0 k=0,1,2,.... z ustalonym N_0=5 licząc: E_N\E_{2N}, E_N=|\int_a^b f - P_N f|
       Porównaj wyniki obliczane kwadraturą prostokątów z wynikami funkcji quad().
  4. Policz dla danej f wielomian najlepszej aproksymacji co najwyżej n w L^2(-1,1) - iloczyn skalarny \int_{-1}^1 fg dx przybliżymy kwadraturą trapezów na 1000 pktach na [-1,1] zastosowaną do f(x)*g(x) - lub stosując funkcję quad(). Stwórz macierz Gramma w bazie jednomianów (w zasadzie znamy tę macierz - jaka funkcja ją tworzy?) - wektor prawej strony (f,x^j) j=0,...,n. Rozwiąż układ. narysuj wykres f i tego wielomianu dla f=sin,sin(4x), funkcja charakterystyczna [-0.5,0.4], |x-0.3|, sqrt(|x-0.5|) etc
  5. Porównaj wyniki obliczane kwadraturą prostokątów z wynikami funkcji obliczającej całkę za pomocą złożonej kwadratury trapezów. Porównaj błędy dla f(x) z poprzedniego zadania czyli całek, których wartość teoretycznie znamy z wynikami dla obu kwadratur:
     • dla tej samej ilości wywołań funkcji f, tzn dla P_N f z T_{N-1} f,
     • porównaj obie kwadratury dla tego samego h tzn. T_Nf z P_Nf aby ocenić błąd w terminach parametru h=(b-a)/N.
  
  calkow19.m testy rzędu zbieżności kwadratury trapezów

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

  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, Trzynaś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 people.maths.ox.ac.uk/trefethen/essays.html
  
  

  ---

  Powrót do mojej strony domowej.

  ---

  Ostatnia aktualizacja: 14 czerwca 2019