Semestr letni 2010/11

• Ćwiczenia/Lab
• Program bieżącego labu
• Zadania domowe, projekty na lab.

Ćwiczenia/Lab MO

• Gr 3 ćwiczenia piątek 830-10, sala 4070. Pierwszy lab 18 luty i potem co 2 tyg
• Gr 4 ćwiczenia - piątek 1215-1345, sala 4070. Pierwszy lab 25 luty i potem co 2 tyg

Zaliczenie ćwiczeń i labu

1. Pierwsza seria ZD arytmetyka fl - na 18 marca 2011
2. Pierwszy projekt z labu (algorytm różnic dzielonych) na piątek 15 kwietnia gr 3 (jak ktoś odda ten projekt również na kolejnym labie pierwszym po feriach, to nie obniżę punktów z uwagi na Święta/ferie) i piątek 29 kwietnia 2011 gr 4 (przesunięcie z uwagi na ferie)
3. Druga seria ZD interpolacja wielomianowa i splajnowa - na piątek 15 kwietnia 2011
4. Trzecia seria ZD kwadratury i równania nieliniowe - na piątek 20 maja 2011. Zmiana terminu o tydzień.
5. Drugi projekt z labu (met. Steffensena) na piątek 27 maja gr 3 i piątek 3 czerwca 2011 gr 4 (ostatnie zajęcia w labie w danej grupie)
6. Czwarta seria ZD (nie ma zadań pisemnych) - układy równań liniowych, LZNK, normy macierzy
Projekt z labu trzeba pokazać w akcji w czasie zajęć labu, nie gwarantuję, że zadania oddane po terminie będą sprawdzone.

Program labów

• Lab 1 (gr 3 - 18/02/10; gr 4 - 25/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 3 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 3 - 4/03/10; gr 4 - 11/03/10) - arytmetyka fl
  1. Wyznaczyć 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)?
  2. 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
  3. Policzyć f(x)=x-\sqrt{1+x*x} raz algorytmem wprost a drugi raz f(x)=\frac{-1}{x+\sqrt{1+x*x}} 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).
  4. Obliczyć wartości funkcji f1(x)=(x-2)^4 i f2(x)=x^4-.....+16 na siatce równomiernej na [2-a,2+a] w wektorze x dla a=1e-3. Narysować wykresy obu funkcji i policzyć błąd \| f1(x)-f2(x) \|_\infty czyli \max_k |f1(x(k)-f2(x(k))|.
  5. 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ć?
  6. Policzyć całki I_n=\int_0^1x^n/(5+x)d x n=0,..,20 z wzoru I_n+5*I_{n-1}=1/n raz biorąc I_0=\log(6/5) i iterując wprzód tzn dla rosnących n a raz biorąc I_{30}=1/180 i iterując w tył. (Oczywiście dla n=30,29,..) będzie kiepsko ale okaże się że potem jest coraz lepiej.
• Lab 3 (gr 3 - 18/03/10; gr 4 - 25/03/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.
  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 na zbieżność interpolacji) 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 [-pi,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 i tych wielomianów dla różnych N - używając funkcji polyval() i plot().
  4. Powtórz (Zadanie na zbieżność interpolacji) 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)arccos(t)) na [-1,1] odpowiednio przesunięte i przeskalowane. Czy błędy e_N dal 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): [LN,eN]=Lagrangeinterp(FCN,a,b,N,type=0), która za 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ó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 na zbieżność interpolacji) 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. Przykład Rungego. Powtórz poprzednie zadanie 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.
• Lab 4 (gr 3 - 1/04/10; gr 4 - 8/04/10) - interpolacja splajnowa, wielomiany Czebyszewa
  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 kubiczny o minimalnym nośniku czyli B-splajn kubiczny. 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 pozostałych węzłów czyli dla k spoza zbioru {-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.
  3. 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^k*N_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 e_N\e_{2*N}. Czy widać, że e_N\e_{2*N} zbiega do 2^p dla jakiegoś p całkowitego np. p=4,8 lub 16?
  4. (Do domu) 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 e_{N}\e_{2*N} są te same? Tzn. czy szybkość zbieżności \|\sin(x)-S_N\|_\infty jest taka sama?
  5. (Do domu) Testowanie eksperymentalne rzędu zbieżności splajnu interpolacyjnego kubicznego naturalnego (warunek brzegowy - zerowanie się drugich pochodnych w końcach odcinku). Powtórz zadanie na badanie rzędu interpolacji splajnami kubicznymi dla sin(x), 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()
  6. Przykład Rungego czyli f(x)=1/(1+x*x) i odcinek [-5,5] a bież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 e_{N}e_{2*N}. Czy e_N\e_{2*N} przybliża 2^p dla jakiegoś p całkowitego?
     Porównaj z wynikami interpolacji Lagrange'a dla węzłów równoodległych (poprzedni lab).
  7. Funkcja octave'a mkpp(). Zapoznaj się z tą funkcją ( help mkpp() ). Utwórz przy pomocy mkpp() splajn kubiczny s na podziale -1\leq0\leq1 odcinka [-1,1] taki, że s jest wielomianem trzeciego stopnia na całym odcinku [-1,1] np. s(x)=x^2 oraz (x+1)^3.
  8. Funkcja octave'a umkpp(). Zapoznaj się z tą funkcją ( help umkpp() ). Utwórz przy pomocy spline() splajn kubiczny s na podziale -1\leq0\leq1 odcinka [-1,1] taki, że s interpoluje wielomian trzeciego stopnia na całym odcinku [-1,1] np. f(x)=(x+1)^3. Następnie sprawdź współczynniki s w bazie {(x-x_k)^j}_{j=0,1,2,3} za pomocą umkpp() na obu przedziałach tzn. dla x_0=-1 na przedziale [-1,0] i x_1=0 na [0,1].
  9. Znajdź za pomocą umkpp() współczynniki splajny z zadania na B-splajn na wszystkich pododcinkach [k,k+1] w bazie {(x-k)^j}_{j=0,...,3} dla k=-2,...,1 czyli tam gdzie ten B-splajn ma nośnik. Porównaj z wynikami otrzymanymi teoretycznie na ćwiczeniach.
• Lab 5 (gr 3 - 15/04/10; gr 4 - 29/04/10 - przesunięcie związane z Wielkanocą) - zaliczenie pierwszego projektu, wielomiany Czebyszewa, kwadratury.
  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=kwadtrapez(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. 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.
• Lab 6 (gr 3 - 13/05/10; gr 4 - 20/05/10) - równania nieliniowe
  1. Znajdź w pomocy funkcję octave'a - rozwiązująca równania nieliniowe - i przetestuj ją na kilku prostych przykładach skalarnych np. cos(x)=0, x^2-2=0, x^10-2=0, x-sin(x)=0 itp
  2. (MN) 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,100,
     • 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.
  3. 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.
  4. 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.
  5. Porównać wyniki z zadania na metodę Newtona z funkcje octave'a fzero() (pomoc: help fzero) lub fsolve.
  6. 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 zadania~ref{zad:newtone-test}.
  7. 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,ldots,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?
  8. 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 f odwrotnej. Jak?)
• Lab 7 (gr 3 - 27/05/10; gr 4 - 3/06/10) - zaliczenie drugiego projektu, numeryczna algebra liniowa
  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 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()!.
  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\|.

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

Ostatnia aktualizacja: 27 lipca 2011