LM semestr zimowy 2016/17

Semestr zimowy 2016/17

Metody numeryczne (dla informatyków)

Numeryczne równania rózniczkowe

Metody numeryczne (dla informatyków)

Ćwiczenia/Lab z MN pt 830-10 na zmiane - sale 5070 lub lab 2042 - terminy wg USOSa

(co 2 tygodnie na zmianę)
Lista punktów z zadań komputerowych i kartkówek - obecnosic w labie - aktualizowana w miarę na bieżąco - plik pdf

Zaliczenie ćwiczeń i labu

Zaliczenie ćwiczeń tablicowych/labu - wg ustaleń koordynatora (wykładowcy) tzn do egzaminu dopuszcza połowa pktów z ćwi i labu. (O ile wiem ocena: 50% egz, 30% ćwiczenia i 20% lab).
Zaliczenie ćwiczeń - serie zdań domowych - sprawdzane poprzez zapowiedziane 20-30min kartkówki (pewnie 2 labo 3) z 2-3 zadaniami z zadanych serii - trzeba mieć co najmniej 3 obecności na ćwiczeniach by móc poprawi kartkówki w razie niezaliczenia ćwiczeń.
Zaliczenie labu - 2 proste projekty po 7 pkt i po 1 pkt za kazda obecnosc na labie. (jakby komus wyszlo ponad 20p to przeniose nadmiarowe pkt do pkt z cw - a jak ktos i z cw ma max ilosc pkt to juz trudno) Projekt oddany po terminie 50% punktów.

W razie usprawiedliwionej nieobecności na kartkówce proszę o kontakt po inne zadania zastępcze.
Na ostatnim labie bedzie mozliwosc poprawienia jednej (i tylko jednej) kartkówki - 25min - zgloszenie udzialu zeruje pkty uzyskane z danej kartkówki.

Serie zadań domowych - projekty zaliczeniowe z labu

kartkówki - pierwsza kartkówka 18-11-2016 (pt 830) - w terminie ćwiczeń - bez zad na LZNK tzn kartk. z zad 1 do 19 (bez oznaczaczonych jako trudne) - zadania nr 8 i 15.
Druga kartkówka w styczniu pt 830 (drugie zajecia w czasie labu razem z zal II projektu) - zadania od 20 z pierwszej serii i z drugiej.
projekt 25-11-2016 (mozna pokazac wczesniej); drugi projekt - pierwszy lab w styczniu - razem z kartkówka.
Na ostatnim labie 27 stycznia 2017 o 830 bedzie mozna poprawic 1 z kartkówek - przystapienie do poprawy zeruje wynik danej kartkówki.

Pierwsza seria - fl, normy, uwarunkowanie u.r.l., LU, QR, LZNK. - plik pdf ZD1.pdf
Projekt z labu - LAB1.pdf - pierwszy lab w listopadzie 25-11-2016
Druga seria - interpolacja wielomianowa, splajny, kwadratury, aproksymacja jednostajna - plik pdf ZD2.pdf
Trzecia seria - metody rozwiazywania równania nieliniowego - plik pdf ZD3.pdf
Drugi projekt z labu - LAB2.pdf - pierwszy lab w 01-2017

Uwaga!

Druga kartkówka na drugich zajeciach na labie w styczniu 2017 - pt 20-01-2017 830 - druga seria do zadan na kwadratury wlacznie (ale to okresle po cw 13-01-2017) i bez zadan oznaczonych jako trudne oraz zad nr 15,20,21 (z teorii aproksymacji której nie bylo na razie na wykladzie) które te tak oznaczylem.
Pierwszy projekt z labu na 25-11-2016
Drugi projekt - na drugich zajeciach na labie w styczniu 2017 - pt 20-01-2017 830 (po kartkówce)
Na ostatnim labie 27 stycznia 2017 o 830 bedzie mozliwosc poprawienia jednej (i tylko jednej) kartkówki - 25min - zgloszenie udzialu zeruje pkty uzyskane z danej kartkówki.
Punkty za kartkówki - projekty w labie - obecnosci w labie (aktualizowane w miarę na bieżąco)- plik pdf
Prosze sprawdzic punkty i w razie niezgodnosci reklamowac.

Bedę starał się umieszczac krótkie opisy tego co było czy ma być na labie czy ćwiczeniach.

Uwaga!

Terminy zajęć wg USOSa - na zmianę - lab/ćwiczenia tablicowe (choć wolne dni wprowadzają zamieszanie).

Program ćwiczeń

Ćwicz 1 Arytmetyka fl. Normy wektorów i macierzy.
1. Znajdz najwieksza, najmniejsza liczbe fl (wsród unormowanych) dla arytmetyki w ktorej mamy 3 bity na ceche, 3 na mantyse i 1 na znak (bias=3). Podaj maksymalna i minimalna roznice pomiedzy kolejnymi liczbami. Powtorz zadanie dla arytmetyk podwojnej/pojedynczej precyzji.
2. Uwarunkowanie wzgledne zadania obliczania prostych funkcji jak np exp(x), (x-2)^2, 1-cos(x), 1-x*x itd
3. Numeryczna poprawnosc algorytmu sumowania
4. Numeryczna poprawno obliczania 1-a*a dwoma algorytmami: (1) f=a*a; f=1-f (2) f=1-a; f=f*(1+a) Czy oba algorytmy daja ten sam blad wzgledny o ile a jest reprezentowana dokladnie w fl?
5. Policz uwarunkowanie oblicznia funkcji \sqrt{x}-\sqrt{x+1} dla x>0, który z alg: (1) wprost ze wzoru (2) z rownowaznego wzoru: 1/(\sqrt{x}+\sqrt{x+1}) jest lepszy z pkt widzenia fl dla x>>0?
6. Optymalne stałe równoważności dla norm p-tych dla p=1,2,\infty.
Ćwicz 2
1. Numeryczna poprawnosc róznicy a^2-b^2 w normie max 2 algorytmami(do domu w 2giej i pierwszej do domu)
2. Pokaż wzór na macierz indukowaną maksimum. Do domu
3. Dla A=A'> 0 istnieje rozkład Choleskiego LDL'(1 na diagonali L D diagonalna z dodatnia), który można otrzymać poprzez klasyczny algorytm eliminacji Gaussa (bez permutacji wierszy/kolumn) i LL' (klasyczny rozkład Choleskiego). Pokaz ze kazdy z nich wynika z drugiego.
4. Znajdz A^{-1} i wspolczynnik uwarunkowania macierz yL g-trójkatnej z 1 na diagonali i -1 na pierwszej naddiagonali (macierz 2 diagonalna). Porównaj koszt 2 algorytmów rozwiazywania Ax=b: a/ x=A^{-1}b b/ Rozwiazujemy uklad z A podstawiajac w tyl tzn x_n=b_n i potem x_k=b_k+x_{k+1} kMamy uklady Ak*x=f_k dla 2 konkretnych macierzy A1=[1 -1;1 1] i A2=[1 -1; 1e-6 1e-6]. Znamy 2 wektory yk takie ze ||fk-A*yk||_\infty=1e-7. Oszacuj ek=||x-y_k||_\infty. W szczególnosci czy bedziemy
5. rozkład PA=LU konkretnej macierzy 3x3 lub 4x4
6. Znajdź rozkład Choleskiego dla A=[ 1, 0, 1, 1; 0, 1, -1, -1; 1, -1, 3, 2; 1, -1, 2, 3] (do domu)
Ćwicz 3 Kontynuacja norm, rozkładów. metody dla macierzy blokowych
1. Norma Frobeniusa nie jest normą indukowaną. Z góry szacuje normę drugą. Znajdź stałą równoważności taką, że norma druga szacuje normę Frobeniusa. (do domu)
2. ( Q ortogonalna (Q'Q=I czyli \|Qx\|_2=\|x\|_2). Pokaż, że \|QA\|_2=\|AQ\|_2=\|A\|_2 oraz to samo dla normy Frobeniusa.
3. Pokaz, ze cond(AB)<=cond(A)cond(B).
4. Rozwiąż dwa układy A_kx_k=b z A_1=[1,1;1,0.98] i A_2=[1,1;1,0.99] i b=(1,2)^T . Policz współczynnik uwarunkowania A_1 w normie maksimum.
5. Rozpatrzmy macierzy w postaci blokowej B=[A, b;c^T a] dla A macierzy mxm nieosobliwej której rozkład LU znamy, b,c wektorów wymiaru m i danego skalaru a. Jak możliwie niskim kosztem rozwiązać url z B (o ile B nieosobliwa)? Podaj warunek na to aby B byla nieosobliwa.
6. Rozpatrzmy macierzy w postaci blokowej G=[A, B;B^T, 0] dla A macierzy m\times m symetrycznej dodatnio określonej której rozkład LU znamy, B wymiaru m \times n o kolumnach tworzących układ liniowo niezależny. Czy ta macierz jest nieosobliwa? jeśli tak to jak układ Gx=b rozwiązać, przy założeniu że wpierw policzymy L czynnik rozkładu Choleskiego A i znajdziemy C takie, że LC=B . Policz koszt w zaleznosci od m/n. (zad ze skryptu P.K.)
7. (do domu raczej) Udowodnij jednoznaczność rozkładu LU o ile istnieje (L dolnotr. z 1 na przekątnej, U górnotr. z elementami dodatnimi na przekątnej)
8. Czy rozkład QR macierzy nieosobliwej jest jednoznaczny o ile znaki elementów diagonali macierzy R są dodatnie? (Do domu) Czy rozklad QR macierzy mxn kolumnami regularnej jest 1-1 (Q mxm ortogonalna)? Czy rozklad QR dla Q kolumnami ortonormlanje jest 1-1?
9. (do domu) Macierz A nieosobliwa znamy jej rozkład QR, LU lub LL^T. Jak obliczyć możliwie tanio A^{-1}?
10. Znajdz wektor Householdera taki ze Hx=y dla x i y 2 wektorow o tej samej normie drugiej.
11. Czy istnieje baza w której macierz Householdera dla danego wektora H. jest diagonalna? Jakie wartosci sa na diagonali? Czy z dokladnoscia do znaku wektory bazy sa wyznaczone jednoznacznie?
Ćwicz 4 LZNK
1. Kartkówka -z 8 i 15 z pierszej serii
2. Znajdz wektor Householdera ze macierz Householdera przeprowadzi wektor a=[1,1,1,-1]' na [4,0,0,0]'. Sprawdz czy rzeczywiscie tak jest. Znajdz metoda Householdra rozklad QR macierzy A=(a) - 1 kolumnowej - i rozwiaz LZNK z b=[1,0,2,0]'. Q wystarczy przedstawic jako iloczyn macierzy Householdera (a kazda m Householdera jest reprezentowana przez wektor Householdera)
3. Znajdz metoda Householdra rozklad QR macierzy A=[1,1,1,-1;0,1,-1,0]' - 2 kolumnowej - i rozwiaz LZNK z b=[1,0,2,0]'
4. Dla danych M punktów (x_k,y_k) chcemy znaleźć krzywą o równaniu a x^2+b y^2-1=0 najlepiej pasującą do tych punktów. tzn taką że \sum_k |a x_k^2+ b y_k^2-1|^2 jest minimalne. Sformułuj to zadanie jako LZNK. Co trzeba założyć o punktach aby to zadanie miało jednoznaczne rozwiązanie? Jaki byłby koszt rozwiązania tego LZNK przy pomocy rozkładu QR metodą Householdera jako O(M)? Rozwiąż to zadanie dla 3 punktów (1,0), (1,1) i (-1,0) . (Pierwsza cz. byla - tylko koszt)
5. (do domu; lab?) Dla danych punktów (-1,1), (0,2), (2,4), (3,3) znajdź prostą y=ax+b najbliższą tym punktom w sensie: \sum_k|a x_k+b-y_k|^2 jest minimalne.
6. (do domu) Niech macierz A trójdiagonalna nieosobliwa nxn . Opracować wersję algorytmu rozkładu QR metodą Householdera pozwalającą rozwiązać układ A x=b kosztem O(n) - o ile to możliwe albo uzasadnić dlaczego to niemożliwe.
7. (do domu) Niech macierz 4 X 3 A=[1,1,1;eps*I], (czyli pierwszy wiersz jedynki, podmacierz złożona z ostatnich trzech wierszy = eps*I ) i b=(3,eps,eps,eps)^T, eps - niezerowy parametr. Rozpatrzmy LZNK z A i b . Czy jest ono regularne? Policz układ równań normalnych. Jaki będzie skutek rozwiązania tego LZNK przy pomocy układu równań normalnych w arytmetyce pojedyńczej precyzji eps=1e-4}?

Ćwicz 5 i 6 nieregularne LZNK, interpolacja Lagrange'a, blad interpolacji L, wielomiany Czebyszewa, splajny, aproksymacja, kwadratury - rzad

Znajdz x^* o najmniejszej normie drugiej rozwiazanie LZNK z macierza mxn A=(a_{i,j}) taka, ze a_{1,1}<>0, a a_{i,j}=0 w innym przypadku (macierz zerowa poza pozycja w lewycm górnym rogu)i wektorem prawej strony b. Jak wyznaczyc ortonormalna baze jadra A?
Jak znalezc x^* o najmniejszej normie drugiej rozwiazanie LZNK z wektorem prawej strony f i macierza mxn znajac B=Q_1,Q_2 macierze ortogonalne takie, ze B=Q_1AQ_2 jest macierza diagonalna o 1 niezerowym elemencie na diagonali tj. B_{1,1}<>0 - jak znalezc baze ortonormalna jadra A (N(A))?
Dla macierzy A m x n rzedu jeden o co najmniej 2 kolumnach (n>2) pokaz ze istnieja niezerowe wektory u,v takie ze A=u*v^T
jak za pomoca 2 macierzy Householdera (z znanymi wektorami H.) sprawdzic LZNK z macierza A=u*v^T i wektorem prawej strony f do równowaznego LZNK z macierza diagonalna (o 1 niezerwoym elemencie na diagonali)
Dla danej konkretnej tabelki np x_k :(-1,0,1,2) a y_k : (1, 0, 3,10) znajdź współczynniki wielomianu interpolacyjnego: w bazie jednomianów rozwiązując wprost układ równań liniowych: \sum_{k=0}^3a_kx_l^k=y_l l=0,1,2,3, w bazie Newtona związanej z \{x_k\}_{k=0}^2 przy pomocy algorytmu różnic dzielonych. Sprawdź wynik.
dla danych 3 węzłów (0,1,3) o krotności (1,3,1) i zadanych odpowiednich wartości funkcji i pochodnych (np dla f(x)=x^4) w tych węzłach znajdź współczynniki wielomianu interpolacyjnego Hermite'a : w bazie Newtona związanej z \{x_k\}_{k=0}^2 przy pomocy algorytmu różnic dzielonych. Sprawdź wynik.
(do domu) Napisz w pseudokodzie odpowiednią wersję algorytmu Hornera obliczania wartości wielomianu zadanego w bazie Newtona dla znanych węzłów dla danego punktu x.
pokaz ze roznica dzielona nie zalezy od kolejnosci wezlów wprost ze wzoru L_nf=L_{n-1}f + f[x_0,...,x_n]\Pi(x-x_j)
Oszacuj blad interpolacji wielomianem liniowym w zaleznosci od regularnosci f.
Oszacuj blad w normie max dla f(x)=sin(x) lub sin(10*x) i L_n f wielomianu interpolujacego f na [-\pi,\pi] w 3 wezlach rownoodleglych.
Ile trzeba wezlow Czebyszewa by na [0,4] blad interpolacji Lagrange'a dla \sin(x) byl co najwyzej 2^{-7}? Powtorz dla $\sin(4*x)$.
Oszacuj blad dla interpolacji Lagrange'a na wezlach Czebyszewa dla \log(1+x) na [0,1] i [0,10]. Podaj te wezly.
Pokaż \|u-L_1u\|_\infty\leq C_k(b-a)^k\|u^{(k)}\|_\infty k=1,2 dla interpolanta liniowego L_1u=u(a)+u[a,b](x-a) dla możliwie małych C_k (f \in C^k).
Oszacuj blad interpolacji Hermite'a dla sin(x) i wezlow -\pi,0,2*\pi - pierwszy ostatni o krotnosci dwa, srodkowy o krotnosci cztery w normie max na [-pi,2*pi].
Oszacuj blad w normie max dla f(x)=sin(x) lub sin(10*x) i L_n f wielomianu interpolujacego f na [-\pi,\pi] w 3 wezlach rownoodleglych.
Ile trzeba wezlow Czebyszewa by na [0,4] blad interpolacji Lagrange'a dla \sin(x) byl co najwyzej 2^{-7}? Powtorz dla $\sin(4*x)$.
Oszacuj blad dla interpolacji Lagrange'a na wezlach Czebyszewa dla \log(1+x) na [0,1] i [0,10]. Podaj te wezly.
Znajdź splajn kubiczny B na podziale równomiernym [-3,3] z h=1 taki, że jest tożsamościowo zero poza (-2,2) i przyjmuje wartości B(-1)=B(1)=1 i B(0)=4 . Wsk: Splajn musi by f parzysta. Stad: interpolacja Hermite'a wpierw na [-2,-1] i analogicznie na [1,2] a potem na [-1,0] i [0,1] .
Wyznacz na podziale równomiernym funkcję B_k taką , że jej nośnik to [x_{k-2},x_{k+2}] i przyjmującą w x_k wartość 4 . Wyznacz równania na s(x)=\sum_{k=-1}^{N+1}c_kB_k dla S splajnu naturalnego interpolującego daną f w (x_k)_{k=0}^N . Wykaż, że c_k wyznaczone jednoznacznie.(do domu)
s splajn kubiczny taki, że na odcinku [x_{k-1},x_k] i [x_{k+1},x_{k+2}] jest równy zero. Pokaż, że na [x_k,x_{k+1}] też jest tożsamościowo zero. Czy tak własność jest prawdziwa dla dowolnych splajnów w S^{2r-2}_{2r-1}(\Pi) ?
Pokaż oszacowanie błędu interpolacji splajnami liniowymi. na podziale równomiernym na [a,b] , tzn. jesli \Pi:\{x_k=a+k*h\} , h=(b-a)/N , to \|u- s_N u\|_\infty\leq C_kh^k\|u^{(k)}\|_\infty \quad k=1,2 dla stałych C_k niezależnych od h,a,b,N . Tu s_Nu splajn interpolacyjny liniowy.
Znajdz wielomian najlepszej aproksymacji w L^2_g(-1,1) dla g(x)=1/\sqrt(1-x*x) dla st <= 2 dla f(x)=x^3-2x. Podaj blad aproksymacji. Wsk: znamy wiel. ortogonalne tzn wielomiany Czebyszewa sa ortogonalne z \|T_n\|_{L^2_g(-1,1)}=\pi/2 n>0 i \pi dla n=0

Ćwicz 7

Kwadratura Gaussa na 1 punkcie tzn kwadratura prostokątów: P_{[a,b]}f=Af(\theta) . Znajdź A,\theta dla którego rząd tej kwadratury dla I(f)=\int_a^bf największy. Bylo ale prosze zrobic zgodnie z teria kwadratur Gaussa... Oszacuj błąd dla takich A,\theta i f\in C^k z \|f^{(k)} \|_\infty \leq 1 dla k=1,2 , oraz pokaż, że I(f)-P_{[a,b]}f=C(b-a)^3f''(\xi) dla C stałej niezależnej od a,b,f
Dla złożonej kwadratury prostokątów ( x_k=a+k*h , h=(b-a)/N ): P_Nf=h*\sum_{k=0}^{N-1}f(x_k+0.5*h) pokaż, że E_N(f)= \int_a^bf dx -P_N f=C h^2(b-a)f''(\xi) dla pewnego \xi\in(a,b) oraz C stałej niezależnej od h,a,b,f . Policz C oraz oszacuj błąd E_N(f) dla f klasy C^1 .
Znajdz kwadrarture oparta na -1,a,1 o max rzedzie dla calki po [-1,1] z waga jeden.
Znajdz kwadrarture oparta na -1,a o max rzedzie dla calki po [-1,1] z waga jeden. (do domu)
Pokaz ze jesli kwadratura na n+1 punktch ma rzad >=n+1 to jest to kwadratura interpolacyjna. (do domu)
Pokaz ze kwadratura Newtona-Cotesa oparta na n+1 punktach dla n+1 nieparzystego ma rzad >=n+2. (do domu)
Ćwicz Skalarne metodu rozwiązywania równań nieliniowych: m. Newtona

metoda iteracji prostej dla równania Keplera: x+0.5*cos(x)=10 tzn x_n=10-0.5*cos(x_{n-1}). Dla jakich x_0 zbiega, oszacowanie błędu dla x_0=10 i n=9.
Metoda Herona (czyli m Newtona dla równania x*x-a=0; a>0), że zbiega dla dowolnego x0; (elementarnie)

Odwracanie bez dzielenia (m. Newtona dla 1/x-a=0 dla jakich x_0 m. Newtona zbiega?). W ilu iteracjach otrzymamy rozwiązanie dla a z (1,2) biorąc x_0=0.75 takie aby błąd względny był na poziomie podwójnej precyzji?

czy iteracje x_{n+1}=cos(x_n) zbiegną dla dowolnego x_0?
Zera wielokrotne: pokaż że dla zera 2-krotnego met. Newtona jest zbieżna lokalnie np. dla f(x)=(x-2)^2=0. Uogólnij dla dowolnej funkcji takiej, że f(x)=f'(x)=0, f''(x) różna od zera - wtedy zachodzi lokalna zbieżność Newtona liniowa (o ile w kolejne iteracje dobrze zdefiniowane).

Policz wzór na metodę Newtona dla funkcji g(x)=f(x)/f'(x). Policz pierwszą pochodną tej funkcji w a dla f(a)=f'(a)=0 z f''(a) nie zerowym.

Pokaz dla cos(x)=0 dla x0 bliskiego \pi/2 metoda newton zbiezna kubicznie.
(metoda Halleya) Czyli metoda Newtona zastosowana do g(x)=[f(x)]/[\sqrt{f'(x)}]. Pokaż że druga pochodna funkcji g(x)=[f(x)]/[\sqrt{f'(x)}] wynosi zero dla x=a zera funkcji f o ile f'(a) nie =0. Pokaż, że wzór na tę metodę (tzn m Newtona dla g) to x_{n+1}=x_n - [2f(x_n)f'(x_n)]/[2(f'(x_n)^2 - f''(x_n)f(x_n)].

Zbadaj zachowanie iteracji metody Newtona dla funkcji f(x)=sgn(x)\sqrt{|x|} dla x_0 niezerowego.

Zadania nie zrobione na ćwiczeniach będą na następnych ćwiczeniach albo do zrobienia w domu.

Zadania, których nie zdążymy zrobić są do domu.

Program labów

Z ewentualnymi linkami do jakiś prostych skryptó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). Wskaźniki do funkcji (handle - operator @). Bisekcja. Zad 1 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. Teraz z tej macierzy '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. Wstaw D z powrotem do C jako główny minor. Policz sin od D. Zapisz D do pliku (binarnego i ASCII) - zamień element D(1,1) na -100 i wczytaj nową macierz do octave'a jako DD. Policz normy (różne) macierzy DD. Zad 2 Policz dyskretną normę max od (sin(x))^2 na [0,1] (wektorowo- czyli bez użycia pętli). Zad 3 Zaimplementuj metodę bisekcji w skrypcie - dla rozwiązania równania cos(x)=0 na odcinku [0,2] - sprawdź jak działa. (pierwszy projekt to rozszerzenie tego zadania)
mnbasic.m - skrypt z podstawami octave'a
Lab 2 - obliczenia w 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. Powtórzyć to zadanie w pojedyńczej precyzji tzn. wymuszając żeby zmienna była typu single (polecenie x=single(x) konwertuje do poj. precyzji)
3. narysuj wykres funkcji f(x)=(x-a)+a dla a=5^p dla p=10,...,30 - koniecznie liczac zgodnei z nawiasami jak napisano!
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ć f(x)=x-\sqrt{1+x*x} raz algorytmem wprost wynikajacym z tego wzoru a raz z wykorzystaniem równoważnego wzoru f(x)=-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óżnice w wyniku? Powtorzyć w arytmetyce pojedynczej precyzji tzn. z wykorzystaniem funkcji single(x).
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 dla 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 dla n bliskich 20 jest coraz lepiej - np. można sprawdzić czy przynajmniej obliczone przybliżenia całek spełniają: 1/(6(n+1)) = \int_0^1 x^n/(1+5) dx <= I_n <=\int_0^1 x^n/(0+5) dx = 1/(5*(n+1))
8. Zastosować bisekcję dla funkcji (x-2)^3 liczonej z wzoru na rozwinięcie dwumianu x^3-...-8 startując z a=2-h i b=2+2h dla h=1e-3, 1e-4,1e-5 z warunkiem stopu że błąd jest mniejszy od 1e-20 (czyli długości odcinka w którym jest rozwiązanie < 2e-20). Narysować wykres tej funkcji liczone z obu wzorów na odcinkach 2+[-h,2*h] dla h=10^{-3},...,10^{-5}.
fl16.m - skrypt z rozwiazaniami niektorych zadan z fl z labu nr 2 (fl)
Lab 3 - bezpośrednie metody rozwiązywania układów równań liniowych; uwarunkowanie macierzy
1. Przetestuj solver octave'a tzn operator backslash \ dla prostego układu równań liniowych z macierzami nieosobliwymi np A1=[1,1;1,0.98] i b=[2;1] , A2=[1,1;1,0.99] i b=[2;1] oraz osobliwą A=[-2,2;-1,1] z tym samym b . Sprawdź czy otrzymane rozwiązanie x=A\ b! jest rzeczywiście rozwiązaniem tzn policz normę A*x-b . Porównaj rozwiązania układów z A1 i A2 . policz macierze odwrotne do A1 i A2 i uwarunkowania tych macierzy.
2. Funkcja inv(). Przeczytaj help do funkcji: inv() . Przetestuj dla A=[1,1;1,-1] czy macierz B obliczona za pomoca tej funkcji rzeczywiscie jest macierz odwrotna? Policz normy pierwsza i Frobeniusa ||A*B-I|| oraz ||B*A-I|| . Zastosuj funkcje do rozwiazywania ukladu równan liniowych Ax=f dla znanego x (liczymy f=A*x ). Policz y jako iloczyn B z f i policz blad 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 losowa dla n=10,50,250,1250 ze znanym rozwiazaniem - oszacuj czas przy pomocy funkcji tic i toc . Porównaj czas i bledy w normie pierwszej dla rozwiazania uzyskanego przy pomocy operator octave'a.
3. Sprawdzić dla konkretnej macierzy 3x3 np A=[1, 2, 3;-1, 0, 1;1, 1, 1] czy odpowiednie macierze górno/dolno trójkątne, permutacji czy ortogonalne uzyskane przy pomocy funkcji octave'a lu() i qr() rzeczywiście zwracają odpowiednie macierze rozkładu LU (PA=LU) lub rozkładu QR (A=QR) np. policzyć pierwszą normę różnicy A-QR dla rozkładu QR i czy Q ortogonalna (policz ||Q'Q-I||_1 i ||QQ'-I||_1) itd.
4. Przetestuj solver octave'a dla macierzy dla macierzy górnotrójkątnej U1 i U2 dla U1 z ze wszystkimi elementami ró?nymi jeden (mozna ją utworzyć U1=triu(ones(n,n)) a U2=2I-U1 dla N=10,20,30 dla znanego rozwiązania np. losowego, czyli policz Uk*sol=fk - rozwiąż w octavie układ z Uk i fk i policz normy (1,2 itd) \|Uk*x-fk\| i \|xk- sol\| dla xk=Uk\fk rozwiązania które wyliczył octave.
  Używając inv() policz macierze odwroten do U1 i U2. Policz uwarunkowania obu macierzy w normie supremum.
5. 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. Sprawdź czy macierz uzyskana przez inv(hilb(n)) i invhilb(n) są rzeczywiście odwrotne do hilb(n) - policz normę macierzową różnic tzn. np. inhilb(n)*hilb(n)-eye(n)
6. Powtórz zadanie z macierza Hilberta ale w arytmetyce pojedynczej precyzji. Musimy to sztucznie wymusic aby zmienne byly zmiennymi pojedynczej precyzji za pomoca funkcji octave'a single() .
7. Stworzyć w octavie macierz trójdiagonalną A z 2 na diagonali i -1 na pod- i nad diagonalą tzn A_{i, i}=2 dla i=1,...,n i A_{i, i+1}=A_{i+1, i}=-1 dla i=1,..,n-1 a pozostałe elementy zero, zarówno 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?). Policzyć rozwiązanie układu równań liniowych Ax=b dla różnych znanych x za pomocą y=inv(A)*b i w=A\b (wczesniej obliczymy b=A*x). Porównać normy Ay-b,x-y,w-x (pierwszą, drugą, supremum) dla N=10,100,1e3,1e4 itp
8. Zaimplementować rozkład LU dla macierzy trójdiagonalnej bez wyboru elementu głównego. Przy czym do funkcji macierz przekazywać jako trzy wektory z przekątnymi. Zastosować do rozwiązania układu równań z macierzą z poprzedniego zadania. Porównać czas rozwiązywania Państwa algorytmem a algorytmem octave'a ( \ ) Rozwiązać układ dla losowego rozwiązania 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ą) .
9. Sprawdzic czy macierz 3diagonalna z 2 na diagonali i -1 na pod i nad diagonali jest symetryczna i dodatnio okreslona z wykorzystaniem funkcji octave'a chol() .
10. (macierz Schura - nastepny lab?) Zaprogramuj rozwiazywanie ukladu rownan liniowych z macierza B=[A ,b; c^T,a] dla A macierzy trojdiagonalnej (lub ktorej rozklad LU znamy) n X n, b,c wektorami dlugosci n i a skalarem. Funkcja ma sprawdzic czy (numerycznie) macierz B nieosobliwa.
11. (do domu) W octavie przetestować eliminacje Gaussa z częściowym wyborem i bez wyboru dla macierzy A=[e, \; 1; 1, \;1] z e=eps/10 (tu eps funkcja octave'a zwracająca epsilon maszynowy) i wektorem prawej strony f=[1;0]^T .
12. (do domu) Zaimplementować rozkład LU dla macierzy trójdiagonalnej z częściowym wyborem elementu głównego. I dalej testować jak poprzednim zadaniu. Zastosować do rozwiązania układu równań z macierzą z poprzedniego zadania A i z macierzą A-1.5 .
13. (do domu) Zaprogramować rozkład Choleskiego dla A=A^T>0 trójdiagonalnej tzn policzyć L dolnotrójkątną z jedną poddiagonalą (czyli reprezentowaną przez dwa wektory z przekątną i podprzekątną) taka że A=L^TL . Zastosować do rozwiązania układu równań z macierzą z 2 na diagonali i -1 na pod i nad diagonalach.
lu.m - skrypt z rozwiazaniami kilku zad z lu i uwarunkowania macierzy
Lab 4 - LZNK, metoda Householdera, rozkład QR
1. Testujemy funkcję qr() octave'a do znalezienia rozkładu macierzy A=[1, 1; 1, 0; 1, -1]. Następnie stosujemy ten rozkład do znalezienia rozwiązania LZNK z A i wektorem prawej strony b=[1;1;1] (powinno wyjść x=[1; 0; 0]), porównujemy z rozwiązaniem operatorem backslash. Stwórz macierz A^TA oraz wektor A^Tb - rozwiąż układ A^TAy=A^Tb - czy y==x? (module błędy zaokrągleń)?
2. Porównaj 2 metody rozwiązania LZNK dla wektora b i macierzy A m x n, m>=n - bedacej pierwszymi n kolumnami macierzy hilberta m x m (Funkcja hilb() tworzy taką macierz ze zmienionymi kolejnościami wierszy/kolumn). Pierwsza - operator backslash z octave'a, druga tworzymy układ równań normalnych A^TA i A^b i rozwiązujemy backslashem układ A^TAy=A^Tb. Przetestować dla m =2*n z n=10,20,30 dla b - pierwsza kolumna A. (Rozwiązanie sol - wersor). Policzyć normę różnicy x_k-sol dla x_k rozwiązania policzonego k-tą metodą. Wsk: Funkcja hilb tworzy taką macierz dla m=n (tylko z odpowienio przepermutownymi kolumnami/wierszami) można stworzyć macierz mxm i wyciac odpowiednią podmacierz.
3. Zastosuj octave'a do rozwiązania zadania znalezienia krzywej zadanej równaniem a x*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 ?
4. Zaprogramuj funkcję 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. (skalar może być parametrem opcjonalnym - jeśli funkcja będzie wywołana z dwoma tylko parametrami co można sprawdzić nargin<3 to tę normę można obliczyć w tej funkcji). Przetestować dla losowych wektorów \vec{x} i \vec{w} , sprawdzić czy \|H_w \vec{x} \|_2=\|\vec{x} \|_2. Można napisać ogólniejszą wersje gdzie x=X macierz N\times 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 tę funkcję.
5. Wez losowy wektor Householdera - stwórz macierz Hw i sprawdz czy Hw symetryczna (||Hw-Hw'||_1==0?) i Hw*Hw=Id, oraz czy dla losowych wektorów x powiedzmy wymiaru 100 - ||Hw*x||_2=||x||_2.
6. (testujemy twierdzenie z wykładu) Znajdź wektor Householdera \vec{w} taki że 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 rzeczywiście tak jest, tzn. czy H_w \vec{u}= l \vec{v} .
7. zaprogramuj metode Householdera dla regresji liniowej czy ogolniej dla dowolnej macierzy A 2 kolumnowej i wektora prawej strony y tzn znajdz w1, w2 wektory Householdera takie ze H2H1A=R macierz górnotrójkatna (H1,H2 - macierze householdera dla wektorow Hous. w1,w2 odpowiednio) - zastosuj znaleziony rozklad do rozwiania LZNK z A i y. Porównaj z rozwiazaniem otrzymanym backslashem.
8. Zastosuj standardowej funkcję octave'a tzn \ oraz metodę korzystającą z funkcji qr(A) która zwraca rozkład QR macierzy ( Q ortogonalna) do rozwiązania zadania znalezienia prostej postaci b+ax najlepiej przybliżającej N punktów punkty (k,1+2*k+\epsilon_k) dla k=1,\ldots,N gdzie \epsilon(\epsilon_1,\ldots,\epsilon_N) losowy wektor za zakresu [-\epsilon,\epsilon] testować dla \epsilon= 1,1e-1,1e-2,1e-3] . (Funkcja rand(n) generuje wektor losowy o rozkładzie jednostajnym na [0,1] w octavie).
9. Zaprogramuj metodę Householdera rozwiązywania A\vec{x}=\vec{b} dla A macierzy trójdiagonalnej różnych wymiarów np. z 2 na diagonali a -1 na pod- i naddiagonalą i jakimś losowym wektorem rozwiązania tak aby koszt znalezienia rozkładu i rozwiązania układu równań liniowych był jak O(N) ?
10. Zaprogramuj metode rozwiazywania LZNK z f wektorem prawej strony i A=u*v' dla danych niezerowych wektorów u,v tzn macierzy rzedu jeden z wykorzystaniem 2 macierzy Householdera H1i H2 takich, ze H1*u=||u||_2*e_1 i H2*v=||v||_2*e_1 (e_1=(1,0,..,0)' wersor w odpowiedniej przestrzeni). Przetestuj dla losowych u,v i f - porównaj z wynikiem uzyskanym backslash. Jak mona rozwiazac zadanie z wykorzystaniem funkcji svd(A)? (tzn jak sie ma wynik svd(A) do H1*A*H2?)
lznk16.m - rozwiazania zadan z labu w tym na metode Householdera - w tym regresja liniowa (LZNK z macierza 2 kolumnowa rozwiazujaca zadanie minimalizacji po a,b \sum_k |a xk+bk-yk|^2 czyli LZNK z A=[\vec{x},\vec{1}] i wektorem prawej strony \vec{y}
Lab 5 - interpolacja wielomianowa lagrange'a, hermite'a i splajnowa
1. Interpolacja Lagrange'a: Sprawdź jak działą funkcja polyfit(). Znajdź wielomian interpolacyjny na 4 węzłach -1,0,1,2 dla wielomianów x^3-1 i x^4-2, narysuj wykres wielomianu interpolacyjnego i tej funkcji.
2. Wykorzystując funkcję polyfit() znajdź wielomiany interpolacyjne dla funkcji sin() dla N węzłów równoodległych i N węzłów Czebyszewa na [-pi,2*pi] dla N=5,10,20,30 . Oblicz dyskretną normę max na siatce 1000 punktowej na tym odcinku, tzn. \max|\sin(x_k)-L_N \sin(x_k)| , gdzie x_k pkt siatki. Następnie narysuj wykresy sin(x) i tych wielomianów dla różnych N - np. używając funkcji polyval() . Węzły Czebyszewa to zera T_{n+1}(t)=\cos([n+1]\mathrm{arccos}(t)) na [-1,1] odpowiednio przesunięte i przeskalowane.
3. Powtórz poprzednie zadanie, ale dla log(1+x) na odcinku [0,10] .
4. Interpolacja Lagrange'a przykład Rungego. Interpolacja Lagrange'a: wykorzystując funkcję polyfit() znajdź wielomiany interpolacyjne L_N f dla funkcji f(x)=1/(1+x*x) na [-5,5] dla N węzłów równoodległych i N węzłów Czebyszewa dla N=5,10,20,30. Oblicz dyskretną normę max na siatce 10000 punktowej na tym odcinku tzn. \max|f(x_k)-L_Nf(x_k)| , gdzie x_k pkt siatki. Następnie narysuj wykresy f(x) i tych wielomianów dla różnych N - np. używając funkcji polyval(). Czy ciąg wielomianów interpolacyjnych zbiega do f ?
5. Interpolacja hermite'a - napisz funkcje która dla danych 2 wezlów (x0,x1) i 4 wartsci (f0,df0,f1,df1) znajdzie wielomian interpolacyjny Hermite'a stopnia <=3 taki, ze w(xk)=fk i w'(xk)=dfk dla k=0,1 - narysuj wykres biorac za fk wartosci sin a dfk wartosci sin'=cos w wezlach x0=-1 i x1=3 - narysuj wykresy. (Wskazówka: najprosciej wykreowac macierz w bazie naturalnej - rozwiazac uklad backslashem).
6. 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] .
7. Splajn bazowy interpolacyjny - Dla danych węzłów równoodległych \{k\}_{k=-5,-4,...,5} na [-5,5] narysuj wykres splajnu kubicznego takiego, że s(0)=1 i s(k)=0 dla pozostałych węzłów k\not =0 oraz który ma pochodne równego zero w węzłach skrajnych tzn. : -5 i 5 . Narysuj jego wykres, jaki jest nośnik tego splajnu? Powtórz zadanie ale dla splajnu typu not-a-knot tzn. nie podając wartości pochodnych w końcach dla funkcji spline().
8. Splajn kubiczny o minimalnym nośniku. Dla danych węzłów równoodległych \{k\}_{k=-5,-4,...,5} na [-5,5] narysuj wykres splajnu kubicznego takiego, że s(-1)=s(1)=1,s(0)=4 i s(k)=0 dla węzłów k\not \in\{-1,0,1\} oraz ma pochodne równego zero w węzłach skrajnych tzn. : -5 i 5 . Czy poza [-2,2] ten splajn jest równy zero? Policz przybliżone normy maksimum na [-5,-2] i [2,5] dla tego splajnu i narysuj jego wykres.
9. Korzystając z funkcji octave'a spline() znajdź współczynniki splajnu interpolacyjnego kubicznego s_N f na N węzłach równoodległych dla f(x)=\sin(x) na odcinku [-pi,2*pi] dla N=5,10,20,40 . Oblicz dyskretną normę maksimum na siatce 1000 punktowej na tym odcinku tzn. E_N=\max|\sin(x_k)-s_N \sin(x_k)|, gdzie x_k pkt siatki, wyprowadź na ekran iloraz E_{2N}/E_N . Następnie narysuj wykresy sin i tych splajnów dla różnych N . (zadanie pomocnicze - jak obliczyć wartość splajnu otrzymanego z funkcji spline() ?)
10. Jak zadanie poprzednie ale dla przykładu Rungego czyli odcinka [0,5] i f(x)=1/(1+x*x) . Czy na wykresie widać że splajny zbiegają do funkcji dla dużych N ?
interp16.m - rozwiazania zadan na interpolacje wielomianowa i splajnowa (kubiczna z war. brzegowymi not-a knot)
Lab 6 kartkówka - zaliczenie projektu nr 2. Wielomiany Czebyszewa czyli aproksymacja jednostajna wieomianu st. n przez wielomiany st <=n, alternans dla aproksymacji jednostajenj wiel. liniowym, kwadratury złożone - prostokątów, trapezów.
1. Narysuj wykresy czterech pierwszych wielomianów Czebyszewa: T_k k=0,1,2,3 wygenerowanych przez regułę trójczłonową: T_{n+1}=2 x T_n-T_{n-1} dla n>0 z T_0=1,T_1(x)=x . Policz norme sup (dyskretna na 10000 pktów) na [-1,1] wielomianu x^{n+1}+\sum_{k<=n}a_kx^k z a_k losowymi i 2^{-n}T_{n+1}? Porównaj.
2. Szukanie p=c*x+d wiel. liniowego najlepszej aproksymacji jednostajnej na [a,b] dla f scisle wypuklej: wtedy alternans [a,x,b] gdzie x ekstremum wewnatrz (a,b) dla f-p tzn taki ze df(x)=(f(b)-f(a))/(b-a)=c (rownanie nieliniowe - prawa strona c wynika z r. alternansu patrz dalej - rozwiazujemy je np funkcja fzero() lub bisekcja (projekt nr 1)) - znajac x mozemy policzyc wspolczynniki wiel. najl. ap. jednostajnej p(x)=c*x+d z wzorów na alternans tzn f(a)-c*a-d=er;f(x)-c*x-d=-er;f(b)-c*b-d=er -? er to blad aproksymacji.
3. zastosuj swoja funkcje z poprzedniego zadania lub skrypt unipoly.m by znalec alternans i wiel. najl. ap. jednostajnej dla f(x)=cos(x) + x^2 na [-1,2.5]. Sprawdz czy alternans nim rzeczywiscie jest.
4. Funkcja quad() w octave'ie - policz całkę z \sin(x) na [0,\pi] i [-1,1] .
5. Policz za pomocą funkcji octave quad() : \int_{-1}^1\frac{1}{\sqrt{1+x*x}}T_n(x)T_m(x)\:d x dla m=2,n=3 .
6. Zaprogramuj w octave funkcję ze złożoną kwadraturę trapezów: T_n f=h(0.5[f(a)+f(b)]+\sum_{k=1}^n f(a+k*h)) \qquad h=(b-a)/(n+1) - parametry funkcji: wskaźnik do funkcji obliczającej f , a,b,n . Przetestuj dla funkcji f , której całkę znamy dla N_k=2^kN_0 k=0,1,2,... z ustalonym N_0=5 licząc: \frac{E_{2N}}{E_N}, \qquad E_N=|\int_a^b f - T_N f| dla funkcji
  - \sin(x) na [0,\pi] i N_0=5 czyli funkcji analitycznej
  - x^{j+0.5} na [0,1] dla j=0,1,2 czyli funkcji w C^j([0,1])
7. Powtorz zadanie dla kwadratury Simpsona
8. Zaimplementuj kwadrature Czebyszewa na [-1,1] (wezly to zera wiel Czebyszewa ale waga w(x)=1/sqrt(1+x*x) - co wazne wspolczynniki sa sobie równe tzn KC_n=A \sum_{k=1}^N f(x_k) dla A=\int_{-1}^1 w(x) dx /n (x_k zera T_n - ntego wiel Czebyszewa)
unipoly.m - skrypt z funkcja znajdujaca wielomian liniowy najlepszej aproskymacji jednostajnej na [a,b] dla f scisle wypuklej np f(x)=cos(x) +2*x^2 lub x^4 + 3*x^2-10
Lab 7 - metody rozwiązywania równań nieliniowych Zadanie własne - metoda potegowa i odwrotna potegowa, metoda QR.
- Zaimplementować metodę Newtona x_{n+1}=x_n -f(x_n)/df(x_n) w octavie i przetestować jej zbieżnośc dla następujących funkcji: x*x-2 z x0=2, x*x*x-27 z x0=27, exp(x)-2 z x0=10,-10,100, sin(x) dla x0=2, (x-2)^k dla k=2,4,8,16 z x0=3, x*x-2 z x0=1e6, 1/x-a dla danych a=0.5,2,4,100 (oczywiście implementując bez dzielenia). Dla wszystkich tych funkcji znamy rozwiązania więc można wyświetlać na ekranie bład e_n=x_n-r (r - rozwiązanie) i |e_n|/|e_{n-1}|^p dla p=1,2,3. Dobierać różne wartości startowe x_0.
- Jak w zad 1 ale dla metody 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.
- Sprawdzić czy metoda iteracji prostych x_n=cos(x_{n-1}) zbiega do x=cos(x). Zbadac czy zbieżność jest liniowa.
- Porównać wyniki z zad 1 z funkcjami octave'a fsolve() i fzero() (pomoc: help fsolve/fzero).
- (do domu) Zaimplementować przybliżoną metode 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 1e-4,1e-7,1e-10 itp Porównać zbieżność z dokładną metodą Newtona (szczególnie ostatnie iteracje) dla funkcji z zadania 1 i z metodą siecznych.
- (do domu) Rozwikływanie funkcji: dla funkcji y(x) zadanej równaniem g(x,y)=2x^2+3y^2-3=0 znaleźć wartości yk=y(x(k)) dla xk=k*h dla k=-N,...,N i h=1/N (N - 10,20,40,...) Znajdować yk rozwiązując układ równań g(xk,yk)=0. Jak w kolejnych krokach dobierać przybliżenie startowe w metodzie Newtona?
- (do domu) Odwracanie funkcji: mamy daną funkcje np sin(x)+2*x znależć wartości funkcji odwrotnej na odcinku [0,5] na siatce k*h dla k=0,..,100. Narysować wykres funkcji. (Sam wykres można narysować dużo prościej bez wyliczania wartości funkcji odwrotnej. Jak?)
- (do domu) Zaimplementować wielowymiarową metodę Newtona w wersji z dokładnym Jakobianem i (do domu) w wersji gdy Jakobian przybliżany różnicami dzielonymi z parametrem h=1e-8. 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.
testnewton.m - m-plik z funkcja testnewton() - testujaca metode Newtona
testNmultzero.m - m-plik z funkcja testNmultzero() - testuje zbieznosc m Newtona dla (x-2)^p(exp(x)+1) czyli funkcji z p-krotnym zerem równym dwa. p-parametr do ustawienia (p=1 zbieznosc kwadratowa; p calkowite <>1 zbieznosc liniowa ze stala = 1- 1/p

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

Tutaj link do stron Octave'a (skąd można ściągnąć kolejną dystrybucje - pod linuxa czy windows)

octave-forge

manual do octave'a w htmlu

Literatura do matlaba

Desmond J. Higham, Nicholas J. Higham, MATLAB Guide, SIAM, 2005
Cleve Moler, Numerical Computing with MATLAB SIAM, 2004, wersja on-line - sekcje to pliki pdf - archiwum skryptow uzywanych w ksiazce
Kermit Sigmon, MATLAB Primer, dostepne za darmo - choc opisuje stara wersje matlaba

Skrypty m-pliki octave'a

I termin - 2010/11
I termin - 2009/10
II termin - 2009/10
Zadania z kolejnych lat (do 2013/14) w skrypcie o MN i MO prof. P. Kiciaka.
Zadania z MN z lat 2014/15 i 2015/16 na stronie prof. Leszka Plaskoty.

Literatura:

[KC2006] D.Kincaid, W.Cheney, Analiza numeryczna, WNT, 2006
[Kic2015] Przemyslaw Kiciak, Skrypt do MO i MN, plik pdf, 2015.
[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

Literatura dodatkowa

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.

http://people.maths.ox.ac.uk/trefethen/essays.html

wykłady i ćwiczenia

Powrót do mojej strony domowej.

Ostatnia modyfikacja: 17 I 2017