Matematyka obliczeniowa 2019-20

• Bieżący lab
• Zadania domowe, stare zadania z kolokwium
• Wyniki zadań domowych, kolokwium etc: po 1 kartkówce, kolokwium. plik-pdf
• zadania z cwiczen (czesc. rozwiazane)
  1. cw z 13 marca 2020 Plik pdf zadania wykladowcy na fl
  2. cw z 20 marca 2020 Plik pdf - tutaj sa zadania na te i poprzednie cwiczenia (czesciowo sie pokrywaja) dalej na fl, uwarunkowanie i numeryczna poprawnosc algorytmu
  3. cw z 27 marca 2020 Plik pdf zadania wykładowcy - metody rozw. układów r. liniowych w tym rozkład LU Choleskiego i QR metodą Housholdera - zadania od str 4.1 trzeba przewinać poprzednie zadania. Ja dodatkowo dam Państwu 2 zadania
     1. Napisz w pseudokodzie algorytm eliminacji Gaussa bez wyboru dla url z macierzą 3diagonalną tzn elementy o wsp. (i,j)t. że |i-j|>1 są zero. Policz ilość operacji.
        Wsk:wykonaj pierwszy krok alg. - zobacz ile operacji jest koniecznych -nie musimy zerować elementów już zerowych...
     2. Zakładając że znamy czynniki rozkładu QR macierzy A nxn zaproponuj możliwie tani algorytm rozwiązania url Bx=f z macierza B=[A,2A';0,4A]. Zapis w notacji blokowej octave 0 odpowiada macierzy zerowej nxn a A' macierzy transponowanej. Podaj koszt jako O(n^p) np p=2 lub 1
        Wsk: Czy tanio da się wyznaczyć rozkład QR macierzy B?
  4. cw z 4 i 9 kwietnia 2020 Plik pdf zadania moje - normy i uwarunkowanie zad rozw. ukl. r liniowych.
  5. Cw z 16 kwietnia 2020 Plik pdf zc5.pdf zadania na LZNK (wykladowcy)- strony 5.1-5.6 a koncu. Dodatkowo zadaniacode mnie:
     • Rozwazmy uklad blokowo zapisany jako:

       r-Ax=b
       A^Tr=0

       Dla macierzy A mxn kolumnami regularnej i b wektora mx1. Pokaz ze ten uklad ma jednoznaczne rozwiazanie a skladowa x jest rozwiazaniem LZNK z A i b. Czy ten uklad ma 1-1 rozwiazanie dla dowolnej macierzy A?
     • Niech A macierz 3diagonalna (n+1)xn kolumnami reg. (zera dla elementow o indeksach (i,j) t. Ze |i-j|>1). Zapisz w pseudokodzie wersje alg Householdera znajdujaca rozklad QR macierzy A kosztem O(n). Oczywiscie Q jako odp iloczyn m. Householdera danych odp. wektorami Househ.
       Jaki bedzie koszt rozwiazania LZNK z A i wektorem b z wykorzystaniem tego rozkladu?
       Wsk: popatrz ile elementow A w 1 kroku trzeba wyzerowac, ile bedzie rozne od zera w H1*A i ile niezerowych ma wektor Househ od macierzy Householdera H1.
  6. Cw z 24 kwietnia 2020: zadania na numeryczne zadanie wlasne -wlasnosci teoretyczne zc6pk.pdf w wiekszosci rozwiazane (od wykladowcy) rozdzial 6.
  7. cw 4 maja 2020 - kartkowka 2 zad z 2 pierwszych serii zad dom! Udostepnie nr zadan o 830 (lub ciut wczesniej) na czacie bedzie 45 min od 830 na rozwiazanie 2 zadan zeskanowanie i zaladowanie na moodlea. Poza tym kilka zadan na numeryczne zad wlasne : zc7.pdf czesc rozwiazana (jest kilka literowek poprawie w wolnej chwili chyba rozwiazania powinny byc jasne.)
  8. ćw 15 maja 2020 - zadania na interpolacje Lagrangea zc8.pdf część rozwiązana (to zad na 1.5-2 terminy)
  9. ćw 22 maja 2020 - zadania na splajny liniowe i kubiczne zc9.pdf część rozwiązana
  10. ćw 29 maja 2020 - aproksymacja jednostajna i średniokwadratowa cyli w normie generowanej iloczynem skalarnym - ta druga też na ostanie ćwiczenia bo chyba na wykładzie jeszcze nie będzie jutro (niestety musimy się spieszyć są tylko 2 terminy ćw) zc10.pdf
  11. ćw 5 czerwca 2020 - dokończyć zadania na aproksymację, i zachęcam zrobić zadania z serii na kwadratury: zc11.pdf

  ---

  Ćwiczenia MO

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

  ---

  Zaliczenie ćwiczeń

  Zasady ustalił wykładowca: prof. Przemysław Kiciak Link do strony - można je znaleźć w jego skrypcie: pdf
  
  Ćwiczenia: 20p za 2kartkówki z wcześniej zadanych zadań w czasie ćwiczeń. Możliwość poprawy 1 z kartkówek (z przeskalowaniem 3/4) - na kartkówce trzeba mieć swoje wydruki treści.

  ---

  Serie zadań domowych i projekty z labu

  Serie ZD zadawane co kilka zajęć. Sprawdzane w formie zapowiedzianych 2-3 kartkówek ze znanych zadań.
  LAb - 8 pktów za 2 projekty (po 4p każdy) + 2p za obecności. Minus 1 pkt za każdą 4,5,6 etc nieobecność na labie (max minus 2). Uwaga Druga kartkówka - ostatnie zajęcia pt 5.6.2020 o 830 w czasie zajęć, Zakres : trzecia i czwarta seria zd. (również możliwa porpaw a 1 kartkówki - szceg. wyżej)
  
  1. Pierwsza seria ZD metody rozwiązywania równań nieliniowych.
  2. Druga seria ZD fl, normy, LU.
  3. Trzecia seria ZD uwarunkowanie cd, Householder, LZNK (nie wchodzi w zakres 1 kartkowki) ?
  4. Czwarta seria ZD zadanie własne, interpolacja wielomianowa i splajnowa
  5. Pierwszy projekt z labu testy metody Steffensena - termin do 8 maja 2020 prosze zaladowac na moodlea
  6. Drugi projekt z labu - algorytm różnic dzielonych dla interpolacji Lagrange'a i algorym Hornera w bazie Newtona - ostatni poniedziałek zajęć 8 czerwca 2020 23:59
  PROSZĘ NIE przysyłać rozwiązań mailem
  Wyniki kolokwium/kartkówki: w moodle'u

  ---

  Zadania z egzaminów z poprzednich lat

• Skrypt P. Kiciaka - pierwsze strony - zadania z egzaminów i kolokwiów
• Zadania z kolokwiów (2gr i poprawkowe) i egzaminu w I terminie z roku 2011/12
• Zadania z kolokwium(2 grupy) z 2012/13

---

Program labów

Z ewentualnymi linkami do jakiś prostych skryptów. Będę starał się tu umieszczać krótkie opisy tego co było czy ma być na labie
• Lab 1 - wstępne zapoznanie się z octavem - octave jako kalkulator naukowy, operator dwukropek : tworzenie macierzy, wektorów, zapisywanie/czytanie do/z plików w obu formatach - tekstowym i binarnym, tworzenie macierzy z podmacierzy, wycinanie podmacierzy itp, podstawowe operacje na macierzach - mnożenie, dodawanie, transponowanie, funkcje od macierzy, normy wektorów/macierzy. Skrypty i funkcje (m-pliki) w octavie. Instrukcje warunkowe (if else endif; switch case endswitch), pętle (while( ) do.. endwhile; do.. until( );for.. endfor). Rysowanie wykresów funkcji - czyli funkcja plot(). Wskaźniki do funkcji (handle - operator @).
  mobasic.m - skrypt z przykładami powyższych operacji czyli tym o czym jest ten lab.
  Kilka zadań (aby nabrać wprawy):
  1. Operacje macierzowe
     • Utwórz dowolne macierze 3x4 A i 3x5 B - a następnie macierz 3x9 C, której pierwsze 4 kolumny to A a kolejne to B.
     • Z macierzy C 'wytnij' podmacierz D składającą się z 1 głównego minora 3x3 tzn. od C(1,1) do C(3,3).
     • Zamień kolejność kolumn D.
     • Zamień kolejność wierszy D.
     • Wytnij dolnotrójkątną część macierzy D a potem górnotrójkątną
     • Wstaw D z powrotem do C jako główny minor.
     • Policz sin(D)=(sin(D_{ij}) od D.
     • Zapisz D do pliku (binarnego i ASCII) - zamień element D(1,1) na -100 i wczytaj nową macierz do octave'a.
  2. Policz dyskretną normę max od (sin(x))^2 na [0,1] (wektorowo- czyli bez użycia pętli).
  3. Narysuj wykres sin(x) na odcinku [0,4].
  4. Znajdź przybliżone max i min funkcji x*(3+2*cos(x)) na [-1,5] oraz przybliżenia punktów ekstremalnych. Zrób to samo dla jakiegoś wielomianu stopnia dwa i trzy np x^3+x^2-x-4
  5. Utwórz funkcję w m-pliku np. obliczającą (sin(x))^2
  6. Utwórz m-plik obliczający wartość funkcji z parametrem np sin(a*x) - parametr przekaż jako zmienną globalną
  7. Funkcje anonimowe - utwórz funkcję (sin(x))^2 jako funkcję anonimową - czyli w linii komend
  8. Przy pomocy pętli for i while oblicz sumę 1+...+N dla N=100. Użyj if aby sprawdzić czy równa się tyle ile powinno 0.5*N*(N+1)- wyprowadź na ekran komunikat używając printf().
• Lab 2 - równania nieliniowe
  1. Bisekcja Przetestuj metodę bisekcji - dla rozwiązania równania cos(x)=0 na odcinku [0,2] - sprawdź jak działa. Mozna samemy zaimplementowac -albo uzyc skryptu testbis.m
  2. Przetestuj metodę bisekcji z poprzedniego zadania do rozwiązania innych problemów: x^2-2=0, exp(x)=2, cos(10*x)=0 startując z odcinka [0,110] - przy dużej ilości zer na danym odcinku znalezione zero jest w sumie przypadkowo wybrane
  3. Solver octave'a Znajdź w pomocy funkcję octave'a - rozwiązująca równania nieliniowe - i przetestuj ją na kilku prostych przykładach skalarnych np.
     1. cos(x)=0 ,
     2. x^2-2=0 ,
     3. x^{10}-2=0 ,
     4. (x-1)*(x-2)*(x-3)=0 ,
     5. x-sin(x)=0 .
     A dokładniej zastosuj tę funkcję do rozwiązywania powyższych równań dla różnych przybliżeń startowych.
  4. Metoda Newtona Zaimplementować metodę Newtona w octavie i przetestować jej zbieżność dla następujących funkcji:
     • x^2-2 z x_0=2 ,
     • x^3-27 z x_0=27 ,
     • exp(x)-2 z x_0=10,-10,1e3 ,
     • sin(x) dla x_0=2 ,
     • cos(x) z x_0=2
     • (x-2)^k x x_0=3 dla k=2,4,8,16 ,
     • x*x-2 z x_0=1e6 , (czy zbiega w ogóle a jak tak to czy kwadratowo)
     • 1/x-a dla danych a=0.5,2,4,100 (oczywiście implementując bez dzielenia) jak dobrać x_0 ?
     Dla wszystkich tych funkcji znamy rozwiązania więc można wyświetlać na ekranie błąd e_n=x_n-r ( r - rozwiązanie) i |e_n|/|e_{n-1}|^p dla p=1,2,3 . Proszę dobierać różne wartości startowe x_0
  5. Metoda siecznych Powtórzyć poprzednie zadanie ale zastępując metodę Newtona przez metodę siecznych w szczególności przetestować czy zachodzi e_n/(e_{n-1}e_{n-2}) asymptotycznie zbiega do stałej i czy |e_n|/|e_{n-1}|^p dla p=(1+ \sqrt{5})/2 zbiega do stałej np. dla x*x-2 .
  6. Iteracje proste Sprawdzić czy metoda iteracji prostych x_n=cos(x_{n-1}) zbiega do x=cos(x) . Zbadać eksperymentalnie czy zbieżność jest liniowa, tzn. czy |x-x_n|/|x-x_{n-1}| zbiega do stałej.
  7. Porównać wyniki z zadania poprzedniego z funkcją octave'a fsolve() (pomoc: help fsolve ).
  8. M. Newtona z przybliżeniem pochodnej Zaimplementować przybliżoną metodę Newtona w której pochodną przybliżamy ilorazem różnicowym tzn. x_{n+1}=x_n - f(x_n)*h/(f(x_n+h)-f(x_n)) dla ustalonego h . Przetestować różne h np. h=1e-4,1e-7,1e-10 itp Porównać zbieżność z dokładną metodą Newtona (szczególnie ostatnie iteracje) dla funkcji z poprzednich zadań.
  9. Funkcja odwrotna: Znajdź przybliżenia funkcji odwrotnej do danej f na siatce równomiernej N punktowej na [c,d] dla c=f(a),d=f(b). Wykonaj dla f(x)=x^2 dla a=1,b=2; f(x)=sin(X) na [0,pi/2].
  10. Rozwikływanie funkcji: dla funkcji y(x) zadanej równaniem g(x,y)=2x^2+3y^2-3=0 znaleźć wartości y_k=y(x_k) dla x_k=k*h dla k=-N,...,N i h=1/N (N - 10,20,40,...) Znajdować y_k rozwiązując układ równań g(x_k,y_k)=0 . Jak w kolejnych krokach dobierać przybliżenie startowe w metodzie Newtona?
  Jak Państwo nie dali rade sami napisać to proszę użyć skryptu testnewton.m w ktorym sa implementacje metody newtona
• Lab 3 - arytmetyka fl
  Celem labu jest wykonanie odpowiednich zadań, które wykażą pewne cechy arytmetyki fl. Domyślnie w octavie wykorzystywane są liczby zmiennopozycyjne o podwójnej precyzji, jakkolwiek w najnowszych wersjach octave'a można również trochę sztucznie wymusić używanie zmiennych w precyzji pojedynczej przy pomocy funkcji single(a) zwracającej zmienną pojedynczej precyzji z tą samą wartością. proszę wpisać z linii komand: format long e;d=1.234567890123456
  #a potem
  s=single(d)
  whos
  #czy wszystkie cyfry sie zachowaly?
  d-double(s)
  Z konstrukcji arytmetyki fl wynika, że odejmowanie dwóch wartości o tym samym znaku o małej różnicy może skutkować dużą utratą dokładności i większość zadań jest po to by zauważyć tą własność.
  1. Funkcje single(a) i eps. Wywołaj pomoc do tych funkcji w octavie. Sprawdź czy 1+eps, 1+2eps 1+eps/2 obliczone w arytmetyce podwójnej precyzji w octavie są większe od jeden. Przyjmij, że a=single(eps) i sprawdź czy ponownie 1+a jest większe od 1 .
  2. Wyznacz epsilon maszynowy czyli najmniejszą liczbę fl taką że po dodania jej do jeden dostajemy coś większego od jeden (w fl oczywiście) 2^{-t} dla t - ilości bitów mantysy. Porównać z eps komendą octave'a. Można to zrobić i w C/C++. Czy wyszło to samo co w octavie (dla liczb typu double)? A co dla liczb w pojedynczej precyzji? Funkcja octave'a single(x) tworzy takie zmienne. Wykorzystując tę funkcję ponownie wyznacz największe n takie że 2^(-n) +1>1 ale w pojedynczej precyzji.
  3. Narysuj wykres funkcji f(x)=(1+a)-a na [0,1] dla różnych a=5^k dla k=1,2,..,20,..30 . Tutaj ważne aby obliczać wartość f(x) z tego wzoru.
  4. Policzyć f(x)=x-\sqrt{1+x*x} raz algorytmem wprost wynikającym z tego wzoru a raz z wykorzystaniem równoważnego wzoru f(x)=\frac{-1}{x+\sqrt{1+x*x}} tzn. Algorytm 1 a= \sqrt{1+x^2};
     w_1=x-a Algorytm 2 a=\sqrt{1+x^2};
     w_2=-1/(x+a) dla x=10^k i k=4,...,10 . Czy widać różnicę w wyniku? Policz różnicę bezwględną tzn ab=|x1-x2| i względny wb=1.0./|x2| (dlaczego spodziewamy się że x2 lepiej obliczone?) Powtórzyć w arytmetyce pojedynczej precyzji tzn. z wykorzystaniem funkcji single(x).
  5. Obliczyć wartości wielomianu f_1(x)=(x-2)^4=f_2(x)=x^4-.....+16 dwoma algorytmami:
     Algorytm 1
     a= (x-2), f_1(x)=a^4;
     Algorytm 2
     f_2(x)=x^4-...+16
     na siatce równomiernej 1000 punktowej na [2-a,2+a] dla a=1e-3 . Narysować wykresy obu funkcji i policzyć błąd || f_1(x)-f_2(x) ||_\infty czyli \max_k |f_1(x(k)-f_2(x(k))| . Tu x(k)=2-a+k*h dla h=a/500 . Wektor x można utworzyć w octavie przy pomocy funkcji linspace().
  6. Policzyć przybliżenie exp(x) z rozwinięcia w szereg \sum_{k=0}^\infty x^k/(k!) (biorąc 100 a potem 1000 elementów szeregu) sprawdzić błąd względny |y-exp(x)|/|exp(x)| dla x od -100 do 100 (np. dla liczb różniących się o dziesięć czyli -100, -80,...,-10, 0, 10,...,80, 100 . Czy błędy dla liczb ujemnych i dodatnich są tego samego rzędu? Jak to można łatwo poprawić?
  7. Policzyć całki I_n=\int_0^1x^n/(5+x)d x n=0,..,20 dwoma algorytmami ze wzoru I_n+5*I_{n-1}=1/n . Pierwszy algorytm przyjmuje I_0= log(6/5) i oblicza z powyższego wzoru kolejne I_n dla n=1,2,3,... . Drugi algorytm wykorzystuje fakt, że \frac{1}{(n+1)6} = I_n= \frac{1}{(n+1)5} zatem w tym algorytmie przyjmujemy za I_{30} jakąkolwiek wartość z tego przedziału np. I_{30}=1/180 i iterujemy w tył, tzn. dla n=30,29,...,20,...,0 . Porównaj wyniki obu algorytmów dla n=0,...,30 oraz czy wyniki spełniają oszacowanie. Dlaczego jeden z algorytmów działa zdecydowanie lepiej w arytmetyce fl? Jako dodatkowe zadanie pozostawimy uzasadnienie wzoru rekurencyjnego i oszacowania I_n wykorzystywanych w algorytmach.
  8. Zastosować algorytm bisekcji dla funkcji (x-r)^3 liczonej z wzoru na rozwinięcie dwumianu x^3-...-r^3 startując z a=r-1e-3 a b=r+1e-3 - r losowa liczba z [1,2]: r=1+rand. Jako warunek zakończenia działania algorytmu przyjmujemy, że błąd jest mniejszy od 1e-20 (czy długości odcinka w którym jest rozwiązanie). Czy algorytm zwraca przybliżenie liczby r jedynego pierwiastka tego wielomianu? Narysować wykres tej funkcji liczone z obu wzorów. tzn. (x-r)^3 i x^3-...-r^3 na odcinkach r+[-h,2*h] dla h=10^{-3},...,10^{-5} . Czy z obu wykresów wynika, że ten wielomian ma tylko jedno miejsce zerowe w otoczeniu r ?
     Metoda bisekcji dla f(x^*)=0 z f ciągłą : startujemy z odcinka [a,b] takiego, że f(a) i f(b) mają różne znaki - liczymy x=(a+b)/2 i f(x) - patrzymy na znak jeśli taki jak f(a) zamieniamy x z a wpp x z b a jeśli f(x)=0 kończymy iteracje... i powtarzamy poprzednią procedurę dla nowego odcinka. Otrzymujemy zstępujący ciąg przedziałów dla których ciąg lewych/prawych końców czy środków zbiegają do zera funkcji f. Biorąc n-ty środek jako n-te przybliżenie xn mamy |x^*-xn| < = 2^(-(n+1))|b-a| dla a,b końców wyjściowego odcinka. Oczywiście jak funkcja ma kilka pierwiastków to mamy gwarancję że proces zbiegnie do jakiegoś pierwiastka.
  
  fl20.m - prosty skrypt z niektórymi rozwiązaniami.
• Lab 4 (pt 17 kwietnia 2020) rozwiązywanie układów liniowych: rozkłady LU, LL', uwarunkowanie macierzy.
  1. Normy macierzowe - policz normy 1,2,Frobeniusa, max dla macierzy diagonalnej i konkretnej niesymetrycznej A=[1;2;1;1].
  2. Operator octave'a \ służący m.in. do rozwiązywaniu układów równań liniowych w octavie. Przetestuj ten operator dla macierzy A=[1,2;3,4] i x=[1;3] z f=A*x czy y=A\f jest rozwiązaniem tego układu. Policz normy: rezydualną ||Ay-f||_p i normę błędu ||y-x||_p dla p=1,2,\infty. Powtórz testy dla jakiejś macierzy osobliwej? Co się wtedy dzieje?
  3. Funkcja inv(). Zapoznaj się z pomocą do funkcji: inv() . Przetestuj dla A=[1,1;1,-1] czy macierz B obliczona za pomocą tej funkcji rzeczywiście jest macierzą odwrotną? Policz normy pierwszą i Frobeniusa ||A*B-I|| oraz B*A-I|| . Zastosuj funkcje do rozwiązywania układu równań liniowych Ax=f dla znanego x (liczymy f=A*x ). Policz y jako iloczyn B z f i policz błędy rezydualny ||Ay-f|| i ||x-y|| w normie pierwszej i drugiej. Powtórz to zadanie dla macierzy N\times N losowej (funkcja octave'a rand() zwraca macierz losową) dla n=10,50,250,1250 ze znanym rozwiązaniem - oszacuj czas przy pomocy funkcji tic i toc . Porównaj czas i błędy w normie pierwszej dla rozwiązania uzyskanego przy pomocy operator octave'a \ .
  4. Funkcje lu(), i chol() . Zapoznaj się z pomocą do funkcji: lu() i chol() . Dla macierzy A=[1,1;1,-1] znajdź jej rozkład PA=LU , rozkład Choleskiego A=L_1^TL_1 . Sprawdź te rozkłady licząc normy macierzowe pierwszą i Frobeniusa błędów np. PA-LU czy A-L_1^TL_1 . Zastosuj te rozkłady do znalezienia rozwiązania układu równań liniowych Ax=f ze znanym rozwiązaniem np. x=[1;1] i f=[2;0] . Policz normy pierwszą i drugą wektorowe pomiędzy x a wynikiem algorytmu w polegającym na zastosowaniu odpowiedniego rozkładu oraz takie same normy rezydualne tzn. normy Aw-f9. Porównaj czas rozwiazywania np tic;[l,u,p]=lu(A);x=u\(l\(p*f));toc, porównaj z \ oraz z chol().
  5. Przetestuj solver octave'a dla macierzy górnotrójkątnych U1 i U2=-U1+2*I dla U1 z jedynkami na diagonali i jedynkami na wszystkich pozycjach nad diagonala (polecenie octave'a triu(ones(n,n)) ją tworzy), tzn. stwórz macierze Uk (k=1,2) dla N=10,20,30,50,100 dla znanego rozwiązania np. losowego czyli sol_k=1 policz Uk*sol=fk - rozwiąż w octavie układ z Uk i fk i policz normy (2 itd) ||Uk*xk-fk|| i ||xk- sol|| dla xk rozwiązania które wyliczył octave (k=1,2).
  6. Przetestuj solver octave'a dla macierzy Hilberta H(N) (polecenie octave'a hilb() ją tworzy), tzn. stwórz macierz H(N) dla N=10,20,30 dla znanego rozwiązania np. stałego równego jeden na każdej pozycji czyli sol_k=1 , czyli policz H(N)*sol=f - rozwiąż w octavie układ z H(N) i f i policz normy (1,2 itd) ||H(N)x-f|| i ||x- sol|| dla x rozwiązania które wyliczył octave. Policz uwarunkowanie macierzy Hilberta dla powyższych N .
  7. Powtórz zadanie z macierzą Hilberta ale w arytmetyce pojedynczej precyzji. Musimy tu trochę sztucznie wymusić aby zmienne były zmiennymi pojedynczej precyzji za pomocą funkcji octave'a single() .
  8. Przetestuj funkcję invhilb() tworzącą macierz odwrotną do macierzy Hilberta (por. zadanie z macierzą Hilberta). Policz normy macierzowe Frobeniusa i pierwsze dla H(N)*iH(N)-I dla dla N=10,20,30 , gdzie iH(N) macierz odwrotna do macierzy Hilberta obliczona przez invhilb() . Przetestuj wykorzystanie iH(N) do rozwiązania układu równań z macierzą Hilberta H(N) , tzn. stwórz macierz H(N) dla N=10,20,30 dla znanego rozwiązania sol , czyli policz H(N)*sol=f - rozwiąż układ H(N)x=f mnożąc F przez iH(N) tzn. x=iH(N)f i policz normy (np. 1,2 ) ||H(N)x-f|| i ||x- sol|| .
  9. Dla macierzy A MXM takiej, że ostatnia kolumna i główna przekątna A mają jedynki, wszystkie elementy ściśle pod przekątną są równe minus jeden a pozostałe elementy są równe zero tzn. A(k,k)=A(j,M)=1; A(j,k)=-1 dla j>k; (jak utworzyć taką macierz w 1 linijce kodu bez użycia pętli?) policz dla znanego wektora r: f=A*r rozwiąż backslashem układ Ax=f - policz residuum Ax-f oraz błąd x-r i ich normy dla M=20,25,..,50. Powtórz zadanie w pojedynczej precyzji dla M=10,15,..,40.
  10. W octavie przetestuj eliminacje Gaussa z częściowym wyborem i bez wyboru dla macierzy A=[e, 1; 1, 1] z e=eps/10 (funkcja octave'a eps zwraca epsilon maszynowy) i wektorem prawej strony f=[1;0]^T . Trzeba zaprogramować eliminację Gaussa bez wyboru dla macierzy 2\times 2 samemu.
  11. Stworzyć w octavie macierz trójdiagonalną A=(a_{i j})_{i,j=1}^n wymiaru n\times n z 2 na diagonali i -1 na pod- i nad diagonalą tzn.

      a_{i, j}= 2 dla i=j
      -1 dla |i-j|=1
      0 dla|i-j|>1
      dla i,j=1,...,n.

      przy pomocy funkcji octave'a diag() jak i wprost w pętli. Porównać czas używając tic i toc dla n=10,100,1000 . Policzyć uwarunkowanie macierzy dla różnych N , policzyć macierz odwrotną przy pomocy inv() (czy też jest trójdiagonalna?).
  12. Sprawdzić czy macierz 3-diagonalna z 2 na diagonali i -1 na pod i nad diagonali jest symetryczna i dodatnio określona z wykorzystaniem funkcji octave'a chol() .
  13. Zaimplementować rozkład LU dla macierzy trójdiagonalnej bez wyboru elementu głównego, tzn. stworzyć własną funkcję octave'a
      function [x,d,l]=lu3diag(a,b,c,f,N)
      Parametry funkcji:
      1. a,b,c - przekątna, pod-przekątna i nad-przekątna macierzy A ,
      2. f - wektor prawej strony,
      3. N długość przekątnej a .
      Funkcja zwraca x - rozwiązanie Ax=f oraz czynniki rozkładu macierzy A=LU czyli diagonalę macierzy górnotrójkątnej U w wektorze d oraz pod-diagonalę L - macierzy dolnotrójkątnej, trójdiagonalnej z jedynkami na diagonali czyli wektor l . Pamiętajmy, że nad-diagonala U równa się na-diagonali A tzn. wektorowi c .
  14. Zastosować funkcję z poprzedniego zadania do uzyskania rozkładu LU macierzy trójdiagonalnej dla N=4,10,15 . Porównać z wynikiem uzyskanym przy pomocy funkcji octave'a lu().
  15. Zastosować funkcję z poprzedniego zadania (2 do tyłu) do rozwiązania układu równań z macierzą z z 2 na diagonali i -1 na pod i nad diagonalach Porównać czas rozwiązywania Państwa algorytmem, a algorytmem octave'a czyli operatorem A\ f!. Rozwiązać układ dla dowolnego rozwiązania np. x=(1,...,1)^T dla N=10,100,1000 czy nawet 10000 . Policzyć błąd rezydualny (dla znanego losowego rozwiązania) i błąd rzeczywisty w różnych normach tzn. ||x-\tilde{x}||_2 i ||f-A*\tilde{x}||_2 gdzie f=A*x a \tilde{x} rozwiązanie obliczone przez octave'a czy Państwa algorytm czy w najgorszy sposób A^{-1}*f (wcześniej musimy policzyć macierz odwrotną).
  W skrypcie rozwiązania zadań niektorych: lu20.m
• Lab 5 (4 maja 2020) - rozkład QR, regularne LZNK, przekształcenie Householdera
  1. Operator octave'a \ służący m.in. do rozwiązywaniu układów równań liniowych w octavie al również LZNK.
     • Przetestuj ten operator dla RLZNK dla macierzy A=[1,1;1,-1;1,3] i x=[1;2] z f=A*x czy y=A\f jest rozwiązaniem tego układu? Policz normy rezydualną ||Ay-f||_2 i normę błędu ||y-x||_2 .
     • Przetestuj ten operator dla nieregularnego LZNK dla macierzy A^T z A=[1,1;1,-1;1,3] i x=[1;1;1] z f=A*x czy y=A\f jest rozwiązaniem tego układu? Policz normy rezydualną ||Ay-f||_2 i normę błędu ||y-x||_2 .
  2. Rozkład QR w octave. Funkcje qr(). (o ile nie bylo)
     1. Zapoznaj się z pomocą do funkcji: qr().
     2. Dla macierzy A=[1,1;1,-1;1,3] znajdź jej rozkład A=QR z pomocą funkcji qr().
     3. Sprawdź czy czy uzyskana Q jest ortogonalna - policz normy Frobeniusa QQ^T-I i Q^TQ-I .
     4. Sprawdź ten rozkład licząc normy macierzowe drugą i Frobeniusa błąd A-QR . Zastosuj ten rozkład do znalezienia rozwiązania LZNK Ax=f ze znanym rozwiązaniem np. x=[1;0] i f=[1;1;1].
     5. Policz normę drugą wektorową pomiędzy x a wynikiem algorytmu w polegającym na zastosowaniu odpowiedniego rozkładu oraz takie same normy rezydualne tzn. normy drugie Aw-f oraz Rw-Q^Tf .
  3. Dla macierzy z poprzednich zadań rozwiąż LZNK - rozwiązując układ r. normalnych - porównaj wyniki.
  4. Układ równań normalnych a rozkład QR czy operator \ . Rozpatrzmy A_{2n,k} pod-macierz 2n\times k macierzy A_{2n,2n} Vandermonde'a dla 2n węzłów równo-odległych na [0,1] . Rozwiąż LZNK z A_{2n,k} z wektorem prawej strony równym pierwszej kolumnie tej macierzy (rozwiązanie to pierwszy wersor) - trzema sposobami:
     1. operatorem \ ,
     2. rozkładem QR uzyskanym funkcją qr(),
     3. rozwiązaniem układu równań normalnych: tworzymy macierz układu równań normalnych B=A_{2n,k}^TA_{2n,k} , wektor prawej strony g=A_{2n,k}^T f układu równań normalnych. Rozwiązujemy układ równań normalnych Bx=g operatorem \ .
     Przetestuj dla n=10,20,40,80 i k=2,4,n . Porównaj
     • czas obliczeń
     • błąd - ||x-y||_2
     • błąd rezydualny ||Ax-f||_2
     dla x rozwiązania dokładnego LZNK, f wektora prawej strony LZNK, y przybliżenia rozwiązania uzyskanego daną metodą.
  5. (to bylo na poprzednim labie dla innej macierzy - do domu ) Rozkład QR a operator \ przy rozwiązywaniu układów równań liniowych czyli porównanie QR a LU zastosowanych do rozwiązywania układów równań liniowych Rozpatrzmy macierz A_{n,n} Vandermonde'a dla n węzłów równo-odległych na [0,1] . Rozwiąż układ równań liniowych z wektorem prawej strony równym pierwszej kolumnie tej macierzy (rozwiązanie to pierwszy wersor) - dwoma sposobami:
     1. operatorem \,
     2. rozkładem QR uzyskanym funkcją qr(),
     Przetestuj dla N=10,20,40,80 . Porównaj
     • czas obliczeń
     • błąd - ||x-y||_2
     • błąd rezydualny ||Ax-f||_2
     dla x rozwiązania dokładnego LZNK, f wektora prawej strony LZNK, y przybliżenia rozwiązania uzyskanego daną metodą.
  6. Zastosuj octave'a do rozwiązania zadania znalezienia krzywej zadanej równaniem ax*x+b*y*y=1 najlepiej pasującej do danych N punktów (x_k,y_k) (możemy np. np. wziąć lekko zaburzone punkty z danej elipsy y_k=\sqrt{1-4*x_k*x_k} + zab_k gdzie zab_k zaburzenie wylosowane z [0,1e-2] czyli mamy N punktów x_k np. x_k=1/k czy x_k=-1+h*k dla h=2/N k=1,...,N i z powyższego wzoru wyliczamy y_k +zab_k . Czy obliczone a i b bliskie 4 i 1 ?
  7. Zaprogramuj funkcje octave'a
     function y=H(x,w,nw)
     która dla danych wektorów \vec{x} i \vec{w} tego samego wymiaru N i skalaru nw=||w||_2 obliczy H_w \vec{x} dla H_w=I-2*\frac{1}{nw}\vec{w}\vec{w}^T czyli przekształcenia (macierzy) Householdera. Przetestować dla losowych wektorów \vec{x} i \vec{w}, sprawdzić czy ||H_w \vec{x} ||_2=||\vec{x} ||_2. Można napisać ogólniejsza wersje gdzie x=X macierz N x M i wtedy wynik H*X (jak to zrobić bez użycia pętli?). Sprawdzić czy ||A||_2=||H*A||_2=||A*H||_2 i ||A||_F=||H*A||_2=||A*H||_F dla losowej macierzy A i H przekształcenia Householdera dla losowego wektora w\not=0 wykorzystując nasza funkcje.
  8. (testujemy twierdzenie z wykładu) Znajdź wektor Householdera \vec{w} taki ze odp przekształcenie Householdera przeprowadza dany wektor \vec{u}\not=0 na dany wektor l*\vec{v}\not=0 dla pewnej dodatniej stałej l. Przetestuj dla dowolnych 2 różnych wektorów czy tak jest. Tzn czy H_w \vec{u}= l \vec{v}.
  9. Wykonaj metodz Householdera rozklad QR macierzy tzn.: A=[1,1;1,-1;1,3] i rozwiaz LZNK uzywajac tego rozkladu z f=[1;1;3] - porównaj z rozwiazaniem otrzymanym backslashem.
     Powtórz dla losowej macierzy Nx2.
  10. Porównaj czas obliczania H*A dla macierzy A mxn (m>=n)i macierzy Householdera (z losowym wektorem Householdera) mxm dla m=10,10,..,10^k i m=1,10,n na 2 matematycznie równowazne sposoby - dobry: HA1=A- x*((2/(x'*x))*(x'*A)) i drugi zly: (dlaczego?) HA2=A- (2/(x'*x))*(x*x')*A. Wyjasnienie - policz koszt obliczeniowy poszczegolnych nawiasow i ilosc wymaganej pamieci (koszt pamieciowy).
  W skrypcie rozwiazania zadan ktore udalo sie zrobic na labie: lznk20.m
• Lab 6 - 15 maja 2020
  Interpolacja Lagrange'a
  Funkcja polyval() czyli funkcja obliczająca wartość wielomianu w jednym lub równocześnie wielu punktach czyli wektorowo. Funkcja polyfit() czyli funkcja obliczająca współczynniki wielomianu interpolacyjnego dla zadanych wartości i węzłów. Funkcje ppval() i spline() działające analogicznie ale dla splajnów kubicznych.
  1. Korzystając z funkcji polyval() narysuj wykres wielomianu x^3+x-2 bez wykorzystania pętli czyli wektorowo.
  2. Test funkcji polyfit(). Wykorzystując funkcję polyfit() znajdź wielomian interpolacyjny L_n F dla funkcji F(x)=sin(x) dla węzłów -1,0,1 oraz -1,0,1,10. Policz wartości różnicy F-L_n F w węzłach. Oraz narysuj wykresy F i L_n F na jednym rysunku korzystając funkcji plot(). Wartości wielomianu oblicz bez użycia pętli z wykorzystaniem funkcji polyval().
  3. Test funkcji polyfit(). Wykorzystując funkcję polyfit() znajdź wielomian bazowy lagrange dla 2 wezla od lewej dla węzłów -1,2,3,4. narysuj wykres. Wartości wielomianu oblicz bez użycia pętli z wykorzystaniem funkcji polyval().
  4. (zadanie o zbieżności interpolacji Lagrange'a) Interpolacja Lagrange'a - zbieżność ciągu wielomianów interpolacyjnych dla węzłów równo-odległych i węzłów Czebyszewa:
     • Wykorzystując funkcję polyfit() znajdź wielomiany interpolacyjne L_N f dla funkcji f=sin()i L_N g dla g=sin(4x) dla N węzłów równo-odległych na [0,2*pi] dla N=4,8,12 etc.
     • Oblicz dyskretną normę maksimum różnicy f-L_N f i analogicznie g-L_N g na siatce tysiąc punktowej na tym odcinku tzn. e_N=\max|f(x_k)-L_N f(x_k)|, gdzie x_k pkt siatki. I policz stosunek e_N/e_{2N} dla N=4,..,32. Czy błędy maleją do zera? jak zachowuje się e_N/e_{2N}?
     • Narysuj na ekranie wykresy sin(x) i tych wielomianów dla różnych N - używając funkcji polyval() i plot().
  5. Powtórz (zadanie o zbieżności interpolacji Lagrange'a) dla tej samej funkcji i tego samego odcinka ale dla węzłów Czebyszewa. Węzły Czebyszewa to zera T_{n+1}(t)=\cos((n+1)\mathrm{arccos}(t)) na [-1,1] odpowiednio przesunięte i przeskalowane tzn jesli (t_k) zera T_{n+1} na [-1,1] to (x_k=(a+b)/2+0.5*(b-a)*t_k) to wezly Czebyszewa na [a,b] . Czy błędy e_N dla węzłów Czebyszewa są mniejsze niż dla węzłów równo-odległych?
  6. Interpolacja Hermite'a - wezły 2krotne - szukamy współczynników w w bazie potęgowej: Utwórz macierz 2(n+1)x 2(n+1) że dla danego wektora różnych n+1 węzłów x=[x(1),..,x(n+1)]^T i 2 wektorów wartości y i dy długości n+1 dla układu \sum_(k=0)^(2n+1) a(k) x(j)^k = y(j) j=0,..,n (odpowiada w(x(j))=y(j))
     \sum_(k=1)^(2n+1) a(k)*k x(j)^(k-1) = dy(j) j=0,..,n (odpowiada w'(x(j))=dy(j))
  7. Zastosuj macierz z poprzedniego punktu do znalezienia współczynników wielomianu interpolacyjnego Hermite'a dla 2 2krotnych węzłów interpolujących np sin(x).
  8. sprawdź czy dla f=log(1+x) wielomian Hermite na 2krotnych n+1 wezłach równodległych na [0,10] zbiega do f tzn czy norma supremum f-H_{2n+1}f zmniejsza się? (testuj dla n=2,3,..,20.
  
  interp20.m kilka rozwiazan
• Lab 7 - 29 maja 2020 Splajny kubiczne
  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]. powtórz podając inne wartości pochodnych np 1000 w lewym końcu, czy wykres blisko tego końca się zmienił?
  2. 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<>0 oraz który ma pochodne równe 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().
  3. 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 nie należących do {-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.
  4. 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 \frac{e_{N}}{e_{2*N}} . Czy widać, że \frac{e_{N}}{e_{2*N}}\approx 2^p dla jakiegoś p całkowitego np. p=4, 8 lub 16 ?
     • powtórz ostani punkt ale dla wartości błedu punktowego w danym punkcie - np x=0.1, x=-pi+1e-4 (blisko lewego końca) tzn. zastąp e_N przez ep_N(x)=|sin(x)-S_N sin(x)|.
  5. Testowanie eksperymentalne rzędu zbieżności splajnu interpolacyjnego kubicznego bez warunków brzegowych (splajn typu not-a-knot). Powtórz poprzednie zadanie 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 jestnarysuj wykres splajnu na węzłach (0,1,...,10) równoodległych takiego, który przyjmuje wartość zero w każdym węźle oprócz węzła 6 gdzie przyjmuje wartość jeden. taka sama? W szczególności spraw jak zachowuje się błąd pochodnej w lewym końcu (w przypadku splajnu albo tę pochodną jakoś obliczymy (jak?) albo przybliżymy ilorazem różnicowym (krok ilorazu musi być kilka rzędow wielkości mniejszy niż antycypowany błąd aproksymacji pochodnej inaczej - eksperyment nie ma sensu).
  6. Testowanie eksperymentalne rzędu zbieżności splajnu interpolacyjnego kubicznego naturalnego (warunek brzegowy - zerowanie się drugich pochodnych w końcach odcinku). Powtórz zadanie poprzednie 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()
  7. Przykład Rungego czyli f(x)=1/(1+x*x) i odcinek [-5,5] a zbież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 \frac{e_{N}}{e_{2*N}} . Czy \frac{e_{N}}{e_{2*N}}\approx 2^p dla jakiegoś p całkowitego?
  8. Splajny liniowe : napisz funkcję obliczającą splajn liniowy interpolujący dane watości dane w wektorze y dla pktów danych w x (wymiar x = wymiar y, punkty w x w porządku rosnącym) - splajn interpolacyjnych liniowy reprezentowany w strukturze kawałakmi wielomianowej takiej w jakiej trzymany jest splajn interpolacyjny kubiczny obliczony funkcją spline(), wartość splajnu liniowego interpolacyjnego wtedy może być obliczona za pomocą funkcji ppval(). WSK: użyj mkpp() - funkcję tworzącą strukture kawałkami wielomianową (podajemy podział i współczynniki w odp. bazie na każdym pod-odcinku) (trochę trudniejsze - do domu)



---

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

• Link do stron Octave'a (skd mona cign kolejn dystrybucje - pod linuxa czy windows)
• octave-forge - rozszerzenia octave'a
  
• Manual do octave'a w htmlu (po angielsku)
• Link do skryptu do labu z mo z octave'em plik pdf

Skrypty m-pliki octave'a

• mobasic.m - prosty skrypt octave'a w którym są przedstawione niektóre podstawowe operacje macierzowe, oraz ptle, instrukcje warunkowe, zapisywanie/wczytywanie do/z plików itp
  testnewton.m skrypt w ktorym sa implementacje metody Newtona - mozna testowac rzad metody o ile znamy dokladne rozwiazanie
  fl20.m - prosty skrypt z niektórymi rozwiązaniami z fl
  lu20.m skrypt z niektórymi rozwiązaniami z labu na rozkład lu,qr etc
  lznk20.m kilka rozwiązań z labu na lznk
  interp20.m kilka rozwiązań na interpolacje wielomianową

  Literatura:

  Podręczniki:
  • [KC2006] D.Kincaid, W.Cheney, Analiza numeryczna, WNT, 2006
  • [Kic2014] Przemysław Kiciak, Metody numeryczne dla informatyków /matematyka obliczeniowa plik pdf, skrypt, 2014/15
  • [Mos2002] Krzysztof Moszyński, Metody numeryczne dla informatyków, skrypt, plik ps, 2002
  • [Pla2002] Leszek Plaskota, Trzynaście wykładów z matematyki obliczeniowej, skrypt, plik pdf, 2002
  • [DJ1982] Maksymilian Dryja, Janina i Michał Jankowscy, Metody numeryczne, WNT, 1982.
  • [FMW2005] Z. Fortuna, B. Macukow, J. Wasowski, Metody numeryczne, WNT, 2005. Wydanie 7.
  Pozycje [Kic2014], [Mos2002] i [Pla2002] to skrypty dostępne dla studentów naszego wydziau. A pozycja [FMW2005] to książka skierowana raczej do studentów politechniki ale większość algorytmów jest w niej opisana. [KC2006] jest podstawowym podręcznikiem.
  Wykład jest wg [Kic2014].
  

  Inne użyteczne linki

  Literatura dodatkowa dla osób zainteresowanych metodami numerycznymi, obejmująca materiał częściowo lub często całkowicie poza zakresem wykładu
  Eseje wyjaśniające mam nadzieję czym jest i czym na pewno nie jest Analiza Numeryczna (czy inaczej Metody Numeryczne)
  
  • Lloyd N.Trefethen, Numerical analysis (The Princeton Companion to Mathematics, Princeton University Press, 2007) plik pdf
  • Lloyd N.Trefethen, The definition of numerical analysis, (SIAM News, Nov 1992) plik pdf.
  Inne eseje tego autora o analizie numerycznej i nie tylko people.maths.ox.ac.uk/trefethen/essays.html
  
  

  ---

  Powrót do mojej strony domowej.

  ---

  Ostatnia aktualizacja: 21 maja 2020