LM semestr zimowy 2017/18

Semestr zimowy 2017/18

Metody numeryczne (dla informatyków)

Numeryczne równania różniczkowe

Metody numeryczne (dla informatyków)

MN Ćwiczenia pt 830-10 - sale 5070 i co2 tyg. lab 1015-1145 lab 3041 (tyg. parzyste) - terminy wg USOSa

Lista punktów z zadań komputerowych i kartkówek - obecnosic w labie - aktualizowana w miarę na bieżąco - plik pdf
Uwaga! ćwiczenia pt 12 stycznia 2018 godz. 8:30-10 - proszę iść do grupy nr 4 w sali 3320.

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, ocena: max(50p egz, 10p ćwiczenia, 20p kolokwium i 20p lab, 2x egzamin - 5p).
Zaliczenie ćwiczeń-labu laczne -min 20p pkt z labu, kolokwium i pozostalych punktów z cwiczen : serie zdań domowych - sprawdzane poprzez zapowiedziane 20-30min kartkówki (pewnie 2) z 2-3 zadaniami z zadanych serii
Zaliczenie labu - kilka (3-4) króciutkich 30min zadanek (z tresci 2 ostatnich labów np zadanie którego nie zrobilismy) sprawdzanych na labie od razu - mniej wiecej co 2 lab 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

Pierwsza seria - fl, normy, uwarunkowanie u.r.l., LU, QR, m/. iteracyjne dla r. liniowych, LZNK, zadanie wlasne - plik pdf zd1.pdf
Druga seria - interpolacja wielomianowa, splajny, kwadratury, aproksymacja jednostajna - plik pdf zd2.pdf

Punkty za kartkówki - zadania z labu (aktualizowane w miarę na bieżąco)- wkrótce 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ęć w labie wg USOSa w tyg. parzystych (wolne dni wprowadzają pewne zamieszanie).

Program ćwiczeń

(plan - na cw moga pojawi sie inne zadania a niektóre zadane do domu za)

Ćwicz 1
1. LU dla macierzy trójdiagonalnej (do domu)
2. 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.
3. rozszerzamy uklad Ax=vb o dodatkowa kolumne i wiersz - jak znajac rozklad LU macierzy A szybko rozwiazac nowy uklad (o ile istnieje rozwiazanie)?
4. rozkład PA=LU kossh-copy-id nkretnej macierzy 3x3 lub 4x4 -
5. Znajdź rozkład Choleskiego dla konkretnej A=[ 1, 0, 1, 1; 0, 1, -1, -1; 1, -1, 3, 2; 1, -1, 2, 3]
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.)
Ćwicz 2 Kontynuacja. Arytmetyka fl.
1. Wzór S-M: B=A+u*v' wzór na macierz odwrotn, jak tanio rozwizac uklad r . z B przy zalozeniu z tanio rozwiazujemy uklady Ay=f
2. 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). (do domu)
3. Z którego wzoru lepiej w fl liczyc: 1-cos(x) czy 2*sin^2(x/2) dla malych x?
4. Numeryczna poprawnosc algorytmu obliczania iloczynu skalarnego
5. Podaj maksymalna i minimalna roznice pomiedzy kolejnymi liczbami. Powtorz zadanie dla arytmetyk podwojnej/pojedynczej precyzji.(do domu)
6. 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? (do domu)
7. Optymalne stałe równoważności dla norm p-tych dla p=1,2,\infty.
Ćwicz 3 Normy wektorów i macierzy.
1. Optymalne stałe równoważności dla norm p-tych dla p=1,2,\infty.
2. policz stale równowaznosci miedzy norma supremum a ||x||=|x_1|+a*|x_2| dla 0< a < 1.
3. Normy p-te macierzy u*v^T p=1,2,\infty.
4. Dla ukladu Ax^*=b mamy dla pewnego x (obliczonego jakims algorytmem) r=b-Ax z ||r||_2=1e-12. Pokaz ze istnieje macierz E postaci u*v^T taka, ze (A+E)x=b z ||E||_2=1e-12 tzn x jest dokladnym rozwiazaniem zadanai z zaburzonymi danymi (a konkretnie zaburzana macierza A) czyli tak jak powinno byc dla algorytmu numerycznie poprawnego.
5. dla macierzy U g-trójkatnej dwudiagonalnej z 1 na glównej diagonali i -1 na superdiagonali znajd macierz odwrotn i policz wsp. uwarunkowania macierzy w normie maksimum.
6. Dla macierzy diagonalnej D pokaz ze ||D||_2=max_k |d_k} dla d_k wartosci na diagonali
7. Pokaz ze dla Q ortogonalnej (QQ^T=I) zachodzi $||QA||_2=||AQ||_2=||A||_2 dla dowolnej macierzy A.
8. Pokazac ||A||_2 <= ||A||_F czyli 2-ga norma macierzy jest mniejsza równa od normy Frobeniusa
9. Pokaż wzór na macierz indukowaną maksimum ||A||_\infty=\max_k \sum_l |a_{kl}|. (do domu) WSK: pokaz najpierw ze ||A||_\infty<=\max_k \sum_l |a_{kl}| a potem znajdz wektor x<>0 taki ze ||Ax||_\infty=||A||_\infty||x||_\infty - co pokazuje nierónosc w 2ga strone.
10. 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)
Ćwicz 4 Uwarunkowanie - metody iteracyjne
1. Uwarunkowanie wzgledne zadania obliczania prostych funkcji jak np exp(x), (x-2)^2, 1-cos(x), 1-x*x itd
2. 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?
3. Pokaz, ze cond(AB)<=cond(A)cond(B).(do domu)
4. Metoda Richardsona - analiza - dla jakich wartosci parametru mamy zbieznosc dla A=A^T>0. Jaki jest parametr optymalny i jaka jest w tesdy szybkosc zbieznosci w || ||_2.
5. A symetryczna z wartosciami wlasnymi w zbiorze {-2,-2,0,1,2,3}. Okresl dla jakich wartosci parametru metoda x^{(n+1)=x^{(n)}+tau*(x^{(n)}+A*A*x^{(n)}-b) zbiezna do x^*. Jaki uklad r. liniowych rozwiazuje wtedy ta granica x^*? Dla jakiej wartosci tau zbieznosc w normie euklidesowej najszybsza? Dla tego optymalnego parametru oszacuj ||x_9-x^*||_2/ ||x_0-x^*||_2.
6. A symetryczna nieosobliwa majaca wartosci wlasne -1 i 3. Znajdz x_0 dla którego metoda Richardsona rozbiega (nie jest zbiezna) w tej sytuacji dla kazdej wartosci parametru tau.
7. Rozwiaż 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. (do domu/nastepne cwiczenia)
Ćwicz 5 Met iteracyjne cd - LZNK
1. Metoda najszybszego spakku. Dla ukladu Ax^*=b z A=A'>0 i danego x_n szukamy x_{n+1}=x_n+tau*(b-Ax_n) takie ze ||x_{n+1}-x^*\|_A=min_tau || x_n+tau*(b-Ax_n)||_A. Pokaz ze x_{n+1} mozna wyzanczyc bez znajomosci x^* - w kazdej iteracji liczac co najwyzej r az iloczyc A razy wektor. Pokaz ze x_n zbiega do x^*. (Mozna i pokazac oszacowanie w || ||_A ale to trudniejsze)
2. Pokaz ze rozwiazanie x^* RLZNK z A mxn i b spelnia uklad równan normalnych A'Ax^*=A'b.
3. Pokaz ze rozwiazanie x^* RLZNK z A mxn i b spelnia uklad (m+n)x(m+n) zapisany blokowo: r+Ax^*=b; A'r=0 dla pewnego wektora r. Czy ten uklad ma 1-1 rozwiazanie?
4. Dla macierzy A=[0ssh-copy-id ,0,2,0,0,0,0; 0,1,0,1,0,1,-1] rozwiaz metoda Householdera LZNK z A i b=[0,0,0,0,0,2,1] tzn znajdz rozklad QR i zastosuje go. (nie trzeba tworzyc Q - wystarczy podac wektory H.)
5. dla danych punktów (x_k,y_k) szukamy takich a,b by te pkty pasowaly do krzywej zadanej równaniem ax^2+b^y^2-1=0 tzn a,b minimalizujace: \sum_k| ax^2+b^y^2-1|^2. Sformuluj zadanie jako LZNK i podaj warunek na pkty by to bylo RLZNK.
Cwicz 6 LZNK cd. Zadanie wlasne.
1. Kartkówka z zadan z serii pierwszej (bez LZNK i zadania wlasnego)
2. dla danych punktów (x_k,y_k) szukamy takich a,b by te pkty pasowaly do krzywej zadanej równaniem ax^2+b^y^2-1=0 tzn a,b minimalizujace: \sum_k| ax^2+b^y^2-1|^2. Sformuluj zadanie jako LZNK i podaj warunek na pkty by to bylo RLZNK.
3. Jak znalezc rozklad QR macierzy A wersja metody Householdera macierzy 3diagonalnej kosztem liniowym O(n)
4. LZNK dla A takiej ze A' jest kolumnami regularna (dualne LZNK) - do domu
5. 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) (do domu)
6. 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]'
7. 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ącej do tych punktów. tzn takie ż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 bylby 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)
8. Niech macierz A trójdiagonalna nieosobliwa nxn . Opracować wersję algorytmu rozkładu QR metod? Householdera pozwalające rozwiązać układ A x=b kosztem O(n) - o ile to możliwe albo uzasadnic dlaczego to niemozliwe.
Cwicz 7 nieregularne LZNK i SVD
1. A=uv^T jak znaleźć SVD A za pomocą co najwyżej 2 mnożeń przez macierze householdera. Zastosuj do rozwiązania LZNK z A i f. Znajdź bazę ortonormalną jądra A.
2. Jak dla A dowolnej macierzy m x n za pomocą m. Householdera znaleźć Q1,Q2 ortogonalne takie, że Q1*A*Q2 jest macierzą g-trójkątną dwudiagonalną.
3. Jak dla A=A^T dowolnej macierzy symetrycznej n x n za pomocą m. Householdera znaleźć Q ortogonalne takie, że Q*A*Q^T jest macierzą trzydiagonalną.
4. (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.
5. (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}?
6. Dla A=[4 1;1 4] symetrycznej znajdź pary własne i określ czy metoda dla x0=[1;2] jest zbieżna - ma podciągi zbieżne.
Cwicz 8 metoda potęgowa - poprawa kolowkium
1. Zadania z kolokwium
2. A symetryczna taka że jeje wartości własne spełniają a_1=-a_2 >|a_3>=..>=|a_n|. Pokaż że x_k iteracje z met. potęgowej ma 2 podciągi zbieżne v1 i v2. Określ jak znając v1 i v2 policzyć wektory własne dla a_1 i a_2. (zakładamy, że x0 nie jest ortogonalne do tych 2 wektorów)
Cwicz 9 i 10 Interpolacja wielomianowa, splajnowa - zadanie najlepszej aproksymacji w normie zadanej iloczynem skalarnym
1. Pokaż że baza Lagrange'a to rzeczywiście baza...
2. 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.
3. 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.
4. (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.
5. 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)
6. Oszacuj blad interpolacji wielomianem liniowym w zaleznosci od regularnosci f.
7. 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.
8. 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)$.
9. Oszacuj blad dla interpolacji Lagrange'a na wezlach Czebyszewa dla \log(1+x) na [0,1] i [0,10]. Podaj te wezly.
10. 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).
11. 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].
12. 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.
13. 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)$.
14. Oszacuj blad dla interpolacji Lagrange'a na wezlach Czebyszewa dla \log(1+x) na [0,1] i [0,10]. Podaj te wezly.
15. 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] .
16. 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)
17. 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) ?
18. 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.
19. Dla f=x^3-x znajdź wielomian najlepszej aproksymacji kwadratowy w normie \int_{-1}^1 |f|^2 dx: na 2sposoby - rozwiązując ukąd z macierzą Gramma - i generująć wielomianu ortgogonalne (kolejnych stopni) - albo procesem ortgonalizacji G_S alb regułą trójczłonową.
20. 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 11 aproksymacja jednostajna wielomianowa tzn w normie sup, wielomiany Czebyszewa, alternans. Optymalne węzły interpolacji tzw wzły Czebyszewa. Cwiczenia 12 stycznia - proszę iść do grupy nr 4 w sali 3320.

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}. Wsk: metoda bisekcji - szukamy zera funkcji ciaglej f: startujemy z odcinka [a,b] takiego, ze f(a)*f(b)<0 - patrzymy na znak f(x) dla x=(a+b)/2, i zamieniamy x z a jesli znak f(x) taki jak f(a) albo x z b w.p.p (no chyba ze trafimy w zero tzn f(x)=0). Otrzymujemy ciag przedzialów [a_k,b_k] o dlugosci 2^{-k}|b-a| takich ze a_k i b_k zbiegaja do zera f (jakiegos jak jest ich wiecej niz jedno). Zatrzymujemy sie gy |a_k-b_k| < TOL dla zadanej tolerancji TOL. Za przyblizenie zera mozemy wziac a_k,b_k lub (a_k+b_k)/2.
fl17.m - skrypt z rozwiazaniami zadan nr 1-6 z fl z labu nr 2 (fl) - wszystkie oprócz pracy domowej czyli zad 7
Lab 3 - bezpośrednie metody rozwiązywania układów równań liniowych; uwarunkowanie macierzy
1. Zadanie przy komputerze - oceniane.
2. 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.
3. 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.
4. 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.
5. Przetestuj solver octave'a dla macierzy górnotrójkątnej U1 i U2 dla U1 z ze wszystkimi elementami równymi jeden (można ją utworzyć poprzez U1=triu(ones(n,n)) a U2=2*Id-U1 (1 na diagonali, -1 ponaddiagonalą 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: \|Uk*x-fk\|_\infty i \|xk- sol\|_\infty dla xk=Uk\fk rozwiązania, które wyliczył octave.
  Używając inv() policz macierze odwrotne do U1 i U2. Policz uwarunkowania obu macierzy w normie supremum. Proszę wyswietli macierze U1, U2 i do nich odwrotne dla N=10.
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. 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)
7. 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() .
8. 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
9. 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ą) .
10. 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() .
11. (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.
12. (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 .
13. (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 .
14. (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.
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?)
Lab 5 - Numeryczne zadanie wlasne.
1. Zadanie komputerowe : zad 5 a labu nr 3
2. Napisz funkcje obliczajaca dominujaca w module wartosc wlasna i jej wektor wlasny za pomoca metody potegowej x(k+1)=A(x(k))/||Ax(k)||_2. Iloraz Rayleigha r(k)=x(k)'*A*x(k). Warunek stopy gdy róznica miedzy kolejnymi ilorazami Rayleigha mniejsze od zadanej tolerancji np 1e-10. Argumenty A macierz kwadratowa, funkcja ma zwraca przyblizona wartosc wlasna (iloraz R.) i wektor wlasny dla niej.
3. Przetestuj metode potegowa dla A=[4,1;1,4] z x0=[2;1]/sqrt(3). Policz bledy dla wektora wlasnego i wartosci wlasnej - czy sa takiego samego rzedu?
4. Napisz funkcje obliczajaca wektor wlasny i wartosc wlasna za pomoca odwrotnej m. potegowej - argumenty macierz A i przyblizenie wartosci wlasnej mu - przetestuj dla A symetrycznej z wartosciami wlasnymi 10,9,..,1 (wygeneruj losowa Q 10x10 i A=QDQ' dla D diagonalnej z wartosciami wlasnymi na diagonali) i mu=8+1e-k k=1,..,16 - czy ptzy ej samej zadanej tolerancji ilosc iteracji taka sama? Przetestuj dla mu=9.5 - wtedy metoda zbiega (ma 2 podciagi zbiezne do sumy i odp. róznicy wartosci wlasnych dla 8 i 9.
5. Przetetuj metode potegowa i odwrotna potegowa dla A diagonalizowalnej ale nie w bazie ortogonalnej tzn A=CDC^{-1} D diagonalna z w.wl. na diagonali. C nieosobliwa. (Wylosuj C i ustal D) Kolumny C to wektory wlasne A. Zobacz czy szybkosc zbieznosci ilorazów Rayleigha i iteracji met. potegowej (odwrotnej potegowej) sa róznego rzedu.
6. Przetetuj metode potegowa dla A biorac x0 ortogonalne do v1 wektoru wlasnego dla dominujacej wartosci wlasnej co do modulu. Czy metoda zbiegnie do wektora wlasnego drugiego co do wielkosci?
7. Wystartuj metode potegowa dla A (np A symetryczna z w.wl. 50,9,..,1) biorac za x0 wektor wlasny dla jakiejs niedominujacej wartosci wlasnej np lambda2=9 - wyswietlaj ilorazy Rayleigha dla kolejnych 50 iteracji (ilosc iteracj izalezy do stosunku 2 najwiekszych w module ww)? Czy stale beda rowne lambda2?
  
  power17.m - funkcja z prosta implementacja metody potegowej
  ipower17.m - funkcja z prosta implementacja odwrotnej metody potegowej
  testpower17.m - skrypt z testami m potegowej i odwrotnej potegowej
Lab 6
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 ?
intL17.m - rozwiazania zadan na interpolacje wielomianowa
spline17.m - rozwiazania zadan na interpolacje splajnowa (kubiczna z war. brzegowymi not-a knot i hermitowskimi)
Lab 7 - metody rozwiązywania równań nieliniowych
- 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,
  - (Zerow k-krotne) (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).
  - atan(x)=arctan(x)=0 proszę spróbować dla jakich x0 metoda zbiegnie biorąc kolejno x0=0.5,1,2,4,8,16,32,64,128,..,2^k (i to samo dla x0 ujemnych). jak ustalicie Państwo takie x0 że dla 2k zbiega a dla 2^(k+1) nie - można połowiąc odcinek znaleźć takie a że dla x0 < a mamy zbieżność a dla x0>a nie mamy... Złóżmy że rozbieżność zachodzi jak przekroczymy max ilość iteracji (np 60) albo błąd rośnie (znamy rozwiązanie x^*=0)
  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. Proszę dobierać różne wartości startowe x_0 poza zaproponowanymi.
- 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 (kolejne zadanie).
- 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. Za x1 proszę wziąć x1=x0-f(x0)*h/(f(x0+h)-f(x0)) dla h małego np h=1e-6 czyli x1 to iteracja m. Newtona z przybliżoną pochodną ilorazem różnicowym.
- 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).
- (wpływ fl na metodę bisekcji) zaimplementuj metodę bisekcji dla f(x)=(x-s)^3=0 ale f obliczanego ze wzoru f(x)=x^3-3*s*x^2+3*s^2*x-s^3. Proszę wziąć odcinek startowy losowy z 0 \in [a,b] (np losując lewy koniec a z s+[-2,-1] a prawy koniec b z s+[3,4]). Zatrzymać iterację gdy długość odcinka w k-tej iteracji [a_k,b_k] jest mniejsza od tol=1e-12,1e-13 etc Proszę testować czy na pewno cały czas a_k < s a b_k > s. (Jak wezmiemy s2 lub 3 to iterracja sie zatrzyma bo mzoemy dostac f(x)==) w fl -stad lepiej s losowe)
- 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?
- 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?)
- Zaimplementować wielowymiarową metodę Newtona w wersji z dokładnym Jakobianem i 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.
testnewton17.m - m-plik z funkcja testnewton17() - testujaca metode Newtona, siecznych i m Newtona z przyblizona pochodna (potrzebuje f newton17.m sieczne17.m newtonp17.m)
newton17.m - m-plik z funkcja newton17() - metode Newtona z testami
newtonp17.m - m-plik z funkcja newtonp17() - metode Newtona z przyblizona pochodna z testami
sieczne17.m - m-plik z funkcja sieczne17() - metode siecznych z testami

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: 15 I 2018