Matematyka obliczeniowa 2017/18

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

Ćwiczenia MO

• ćwiczenia piątek gr1 830-10 sala 5070 i 1015-1145 gr 2 12:15-13:45 sala 3250 - lab gr2/1 na zmianę piątek 1015-1145 sala 2042

Zaliczenie ćwiczeń

1. Pierwsza seria ZD fl, normy, LU.
2. Druga seria ZD uwarunkowanie cd, Householder, LZNK
3. Trzecia seria ZD - interpolacja Lagrange'a, interpolacja Hermite'a,
4. Czwarta seria ZD - aproksymacje średniokwadratowa (w przestrzeniach z iloczynem skalarnym), jednostajna,
5. Piata seria ZD splajny liniowe i kubiczne i kwadratury
6. Szósta seria ZD metody rozwiazywania równan nieliniowych.

Zadania z egzaminów z poprzednich lat

Program labów

• Lab 1 - wstępne zapoznanie się z octavem - octave jako kalkulator naukowy, operator dwukropek : tworzenie macierzy, wektorów, zapisywanie/czytanie do/z plików w obu formatach - tekstowym i binarnym, tworzenie macierzy z podmacierzy, wycinanie podmacierzy itp, podstawowe operacje na macierzach - mnożenie, dodawanie, transponowanie, funkcje od macierzy, normy wektorów/macierzy. Skrypty i funkcje (m-pliki) w octavie. Instrukcje warunkowe (if else endif; switch case endswitch), pętle (while( ) do.. endwhile; do.. until( );for.. endfor). Rysowanie wykresów funkcji - czyli funkcja plot(). Wskaźniki do funkcji (handle - operator @).
  mobasic.m - skrypt z przykładami powyższych operacji czyli tym o czym jest ten lab.
  Kilka zadań (aby nabrać wprawy):
  1. Operacje macierzowe
     • Utwórz dowolne macierze 3x4 A i 3x5 B - a następnie macierz 3x9 C, której pierwsze 4 kolumny to A a kolejne to B.
     • Z macierzy C 'wytnij' podmacierz D składającą się z 1 głównego minora 3x3 tzn. od C(1,1) do C(3,3).
     • Zamień kolejność kolumn D.
     • Zamień kolejność wierszy D.
     • Wytnij dolnotrójkątną część macierzy D a potem górnotrójkątną
     • Wstaw D z powrotem do C jako główny minor.
     • Policz sin(D)=(sin(D_{ij}) od D.
     • Zapisz D do pliku (binarnego i ASCII) - zamień element D(1,1) na -100 i wczytaj nową macierz do octave'a.
  2. Policz dyskretną normę max od (sin(x))^2 na [0,1] (wektorowo- czyli bez użycia pętli).
  3. Narysuj wykres sin(x) na odcinku [0,4].
  4. Znajdź przybliżone max i min funkcji x*(3+2*cos(x)) na [-1,5] oraz przybliżenia punktów ekstremalnych. Zrób to samo dla jakiegoś wielomianu stopnia dwa i trzy np x^3+x^2-x-4
  5. Utwórz funkcję w m-pliku np. obliczającą (sin(x))^2
  6. Utwórz m-plik obliczający wartość funkcji z parametrem np sin(a*x) - parametr przekaż jako zmienną globalną
  7. Funkcje anonimowe - utwórz funkcję (sin(x))^2 jako funkcję anonimową - czyli w linii komend
  8. Przy pomocy pętli for i while oblicz sumę 1+...+N dla N=100. Użyj if aby sprawdzić czy równa się tyle ile powinno 0.5*N*(N+1)- wyprowadź na ekran komunikat używając printf().
• Lab 2 - arytmetyka fl
  Celem labu jest wykonanie odpowiednich zadań, które wykażą pewne cechy arytmetyki fl. Domyślnie w octavie wykorzystywane są liczby zmiennopozycyjne o podwójnej precyzji, jakkolwiek w najnowszych wersjach octave'a można również trochę sztucznie wymusić używanie zmiennych w precyzji pojedynczej przy pomocy funkcji single(a) zwracającej zmienną pojedynczej precyzji z tą samą wartością. proszę wpisać z linii komand: format long e;d=1.234567890123456
  #a potem
  s=single(d)
  whos
  #czy wszystkie cyfry sie zachowaly?
  d-double(s)
  Z konstrukcji arytmetyki fl wynika, że odejmowanie dwóch wartości o tym samym znaku o małej różnicy może skutkować dużą utratą dokładności i większość zadań jest po to by zauważyć tą własność.
  1. Funkcje single(a) i eps. Wywołaj pomoc do tych funkcji w octavie. Sprawdź czy 1+eps, 1+2eps 1+eps/2 obliczone w arytmetyce podwójnej precyzji w octavie są większe od jeden. Przyjmij, że a=single(eps) i sprawdź czy ponownie 1+a jest większe od 1 .
  2. Wyznacz epsilon maszynowy czyli najmniejszą liczbę fl taką że po dodania jej do jeden dostajemy coś większego od jeden (w fl oczywiście) 2^{-t} dla t - ilości bitów mantysy. Porównać z eps komendą octave'a. Można to zrobić i w C/C++. Czy wyszło to samo co w octavie (dla liczb typu double)? A co dla liczb w pojedynczej precyzji? Funkcja octave'a single(x) tworzy takie zmienne. Wykorzystując tę funkcję ponownie wyznacz największe n takie że 2^(-n) +1>1 ale w pojedynczej precyzji.
  3. Narysuj wykres funkcji f(x)=(1+a)-a na [0,1] dla różnych a=5^k dla k=1,2,..,20,..30 . Tutaj ważne aby obliczać wartość f(x) z tego wzoru.
  4. Policzyć f(x)=x-\sqrt{1+x*x} raz algorytmem wprost wynikającym z tego wzoru a raz z wykorzystaniem równoważnego wzoru f(x)=\frac{-1}{x+\sqrt{1+x*x}} tzn. Algorytm 1 a= \sqrt{1+x^2};
     w_1=x-a Algorytm 2 a=\sqrt{1+x^2};
     w_2=-1/(x+a) dla x=10^k i k=4,...,10 . Czy widać różnicę w wyniku? Policz różnicę bezwględną tzn ab=|x1-x2| i względny wb=1.0./|x2| (dlaczego spodziewamy się że x2 lepiej obliczone?) Powtórzyć w arytmetyce pojedynczej precyzji tzn. z wykorzystaniem funkcji single(x).
  5. Obliczyć wartości wielomianu f_1(x)=(x-2)^4=f_2(x)=x^4-.....+16 dwoma algorytmami:
     Algorytm 1
     a= (x-2), f_1(x)=a^4;
     Algorytm 2
     f_2(x)=x^4-...+16
     na siatce równomiernej 1000 punktowej na [2-a,2+a] dla a=1e-3 . Narysować wykresy obu funkcji i policzyć błąd || f_1(x)-f_2(x) ||_\infty czyli \max_k |f_1(x(k)-f_2(x(k))| . Tu x(k)=2-a+k*h dla h=a/500 . Wektor x można utworzyć w octavie przy pomocy funkcji linspace().
  6. Policzyć przybliżenie exp(x) z rozwinięcia w szereg \sum_{k=0}^\infty x^k/(k!) (biorąc 100 a potem 1000 elementów szeregu) sprawdzić błąd względny |y-exp(x)|/|exp(x)| dla x od -100 do 100 (np. dla liczb różniących się o dziesięć czyli -100, -80,...,-10, 0, 10,...,80, 100 . Czy błędy dla liczb ujemnych i dodatnich są tego samego rzędu? Jak to można łatwo poprawić?
  7. Policzyć całki I_n=\int_0^1x^n/(5+x)d x n=0,..,20 dwoma algorytmami ze wzoru I_n+5*I_{n-1}=1/n . Pierwszy algorytm przyjmuje I_0= log(6/5) i oblicza z powyższego wzoru kolejne I_n dla n=1,2,3,... . Drugi algorytm wykorzystuje fakt, że \frac{1}{(n+1)6} = I_n= \frac{1}{(n+1)5} zatem w tym algorytmie przyjmujemy za I_{30} jakąkolwiek wartość z tego przedziału np. I_{30}=1/180 i iterujemy w tył, tzn. dla n=30,29,...,20,...,0 . Porównaj wyniki obu algorytmów dla n=0,...,30 oraz czy wyniki spełniają oszacowanie. Dlaczego jeden z algorytmów działa zdecydowanie lepiej w arytmetyce fl? Jako dodatkowe zadanie pozostawimy uzasadnienie wzoru rekurencyjnego i oszacowania I_n wykorzystywanych w algorytmach.
  8. Zastosować algorytm bisekcji dla funkcji (x-r)^3 liczonej z wzoru na rozwinięcie dwumianu x^3-...-r^3 startując z a=r-1e-3 a b=r+1e-3 - r losowa liczba z [1,2]: r=1+rand. Jako warunek zakończenia działania algorytmu przyjmujemy, że błąd jest mniejszy od 1e-20 (czy długości odcinka w którym jest rozwiązanie). Czy algorytm zwraca przybliżenie liczby r jedynego pierwiastka tego wielomianu? Narysować wykres tej funkcji liczone z obu wzorów. tzn. (x-r)^3 i x^3-...-r^3 na odcinkach r+[-h,2*h] dla h=10^{-3},...,10^{-5} . Czy z obu wykresów wynika, że ten wielomian ma tylko jedno miejsce zerowe w otoczeniu r ?
     Metoda bisekcji dla f(x^*)=0 z f ciągłą : startujemy z odcinka [a,b] takiego, że f(a) i f(b) mają różne znaki - liczymy x=(a+b)/2 i f(x) - patrzymy na znak jeśli taki jak f(a) zamieniamy x z a wpp x z b a jeśli f(x)=0 kończymy iteracje... i powtarzamy poprzednią procedurę dla nowego odcinka. Otrzymujemy zstępujący ciąg przedziałów dla których ciąg lewych/prawych końców czy środków zbiegają do zera funkcji f. Biorąc n-ty środek jako n-te przybliżenie xn mamy |x^*-xn| < = 2^(-(n+1))|b-a| dla a,b końców wyjściowego odcinka. Oczywiście jak funkcja ma kilka pierwiastków to mamy gwarancję że proces zbiegnie do jakiegoś pierwiastka.
  
  fl18.m - prosty skrypt z niektórymi rozwiazaniami.
• Lab 3 (kartkówka - nastepny lab)
  rozwiązywanie układów liniowych: rozkłady LU, LL', uwarunkowanie macierzy.
  1. Normy macierzowe - policz normy 1,2,Frobeniusa, max dla macierzy diagonalnej i konkretnej niesymetrycznej A=[1;2;1;1].
  2. Operator octave'a \ służący m.in. do rozwiązywaniu układów równań liniowych w octavie. Przetestuj ten operator dla macierzy A=[1,2;3,4] i x=[1;3] z f=A*x czy y=A\f jest rozwiązaniem tego układu. Policz normy: rezydualną ||Ay-f||_p i normę błędu ||y-x||_p dla p=1,2,\infty. Powtórz testy dla jakiejś macierzy osobliwej? Co się wtedy dzieje?
  3. Funkcja inv(). Zapoznaj się z pomocą do funkcji: inv() . Przetestuj dla A=[1,1;1,-1] czy macierz B obliczona za pomocą tej funkcji rzeczywiście jest macierzą odwrotną? Policz normy pierwszą i Frobeniusa ||A*B-I|| oraz B*A-I|| . Zastosuj funkcje do rozwiązywania układu równań liniowych Ax=f dla znanego x (liczymy f=A*x ). Policz y jako iloczyn B z f i policz błędy rezydualny ||Ay-f|| i ||x-y|| w normie pierwszej i drugiej. Powtórz to zadanie dla macierzy N\times N losowej (funkcja octave'a rand() zwraca macierz losową) dla n=10,50,250,1250 ze znanym rozwiązaniem - oszacuj czas przy pomocy funkcji tic i toc . Porównaj czas i błędy w normie pierwszej dla rozwiązania uzyskanego przy pomocy operator octave'a \ .
  4. Funkcje lu(), i chol() . Zapoznaj się z pomocą do funkcji: lu() i chol() . Dla macierzy A=[1,1;1,-1] znajdź jej rozkład PA=LU , rozkład Choleskiego A=L_1^TL_1 . Sprawdź te rozkłady licząc normy macierzowe pierwszą i Frobeniusa błędów np. PA-LU czy A-L_1^TL_1 . Zastosuj te rozkłady do znalezienia rozwiązania układu równań liniowych Ax=f ze znanym rozwiązaniem np. x=[1;1] i f=[2;0] . Policz normy pierwszą i drugą wektorowe pomiędzy x a wynikiem algorytmu w polegającym na zastosowaniu odpowiedniego rozkładu oraz takie same normy rezydualne tzn. normy Aw-f . Porównaj czas rozwiazywania np tic;[l,u,p]=lu(A);x=u\(l\(p*f));toc, porównaj z \ oraz z chol().
  5. Przetestuj solver octave'a dla macierzy górnotrójkątnych U1 i U2=-U1+2*I dla U1 z jedynkami na diagonali i jedynkami na wszystkich pozycjach nad diagonala (polecenie octave'a triu(ones(n,n)) ją tworzy), tzn. stwórz macierze Uk (k=1,2) dla N=10,20,30,50,100 dla znanego rozwiązania np. losowego czyli sol_k=1 policz Uk*sol=fk - rozwiąż w octavie układ z Uk i fk i policz normy (2 itd) ||Uk*xk-fk|| i ||xk- sol|| dla xk rozwiązania które wyliczył octave (k=1,2).
  6. Przetestuj solver octave'a dla macierzy Hilberta H(N) (polecenie octave'a hilb() ją tworzy), tzn. stwórz macierz H(N) dla N=10,20,30 dla znanego rozwiązania np. stałego równego jeden na każdej pozycji czyli sol_k=1 , czyli policz H(N)*sol=f - rozwiąż w octavie układ z H(N) i f i policz normy (1,2 itd) ||H(N)x-f|| i ||x- sol|| dla x rozwiązania które wyliczył octave. Policz uwarunkowanie macierzy Hilberta dla powyższych N .
  7. Powtórz zadanie z macierzą Hilberta ale w arytmetyce pojedynczej precyzji. Musimy tu trochę sztucznie wymusić aby zmienne były zmiennymi pojedynczej precyzji za pomocą funkcji octave'a single() .
  8. Przetestuj funkcję invhilb() tworzącą macierz odwrotną do macierzy Hilberta (por. zadanie z macierzą Hilberta). Policz normy macierzowe Frobeniusa i pierwsze dla H(N)*iH(N)-I dla dla N=10,20,30 , gdzie iH(N) macierz odwrotna do macierzy Hilberta obliczona przez invhilb() . Przetestuj wykorzystanie iH(N) do rozwiązania układu równań z macierzą Hilberta H(N) , tzn. stwórz macierz H(N) dla N=10,20,30 dla znanego rozwiązania sol , czyli policz H(N)*sol=f - rozwiąż układ H(N)x=f mnożąc F przez iH(N) tzn. x=iH(N)f i policz normy (np. 1,2 ) ||H(N)x-f|| i ||x- sol|| .
  9. Dla macierzy A MXM takiej, że ostatnia kolumna i główna przekątna A mają jedynki, wszystkie elementy ściśle pod przekątną są równe minus jeden a pozostałe elementy są równe zero tzn. A(k,k)=A(j,M)=1; A(j,k)=-1 dla j>k; (jak utworzyć taką macierz w 1 linijce kodu bez użycia pętli?) policz dla znanego wektora r: f=A*r rozwiąż backslashem układ Ax=f - policz residuum Ax-f oraz błąd x-r i ich normy dla M=20,25,..,50. Powtórz zadanie w pojedynczej precyzji dla M=10,15,..,40.
  10. W octavie przetestuj eliminacje Gaussa z częściowym wyborem i bez wyboru dla macierzy A=[e, 1; 1, 1] z e=eps/10 (funkcja octave'a eps zwraca epsilon maszynowy) i wektorem prawej strony f=[1;0]^T . Trzeba zaprogramować eliminację Gaussa bez wyboru dla macierzy 2\times 2 samemu.
  11. Stworzyć w octavie macierz trójdiagonalną A=(a_{i j})_{i,j=1}^n wymiaru n\times n z 2 na diagonali i -1 na pod- i nad diagonalą tzn.
      a_{i, j}= 2 dla i=j
      -1 dla |i-j|=1
      0 dla|i-j|>1
      dla i,j=1,...,n. przy pomocy funkcji octave'a diag() jak i wprost w pętli. Porównać czas używając tic i toc dla n=10,100,1000 . Policzyć uwarunkowanie macierzy dla różnych N , policzyć macierz odwrotną przy pomocy inv() (czy też jest trójdiagonalna?).
  12. Sprawdzić czy macierz 3-diagonalna z 2 na diagonali i -1 na pod i nad diagonali jest symetryczna i dodatnio określona z wykorzystaniem funkcji octave'a chol() .
  13. Zaimplementować rozkład LU dla macierzy trójdiagonalnej bez wyboru elementu głównego, tzn. stworzyć własną funkcję octave'a
      function [x,d,l]=lu3diag(a,b,c,f,N)
      Parametry funkcji:
      1. a,b,c - przekątna, pod-przekątna i nad-przekątna macierzy A ,
      2. f - wektor prawej strony,
      3. N długość przekątnej a .
      Funkcja zwraca x - rozwiązanie Ax=f oraz czynniki rozkładu macierzy A=LU czyli diagonalę macierzy górnotrójkątnej U w wektorze d oraz pod-diagonalę L - macierzy dolnotrójkątnej, trójdiagonalnej z jedynkami na diagonali czyli wektor l . Pamiętajmy, że nad-diagonala U równa się na-diagonali A tzn. wektorowi c .
  14. Zastosować funkcję z poprzedniego zadania do uzyskania rozkładu LU macierzy trójdiagonalnej dla N=4,10,15 . Porównać z wynikiem uzyskanym przy pomocy funkcji octave'a lu().
  15. Zastosować funkcję z poprzedniego zadania (2 do tyłu) do rozwiązania układu równań z macierzą z z 2 na diagonali i -1 na pod i nad diagonalach Porównać czas rozwiązywania Państwa algorytmem, a algorytmem octave'a czyli operatorem A\ f!. Rozwiązać układ dla dowolnego rozwiązania np. x=(1,...,1)^T dla N=10,100,1000 czy nawet 10000 . Policzyć błąd rezydualny (dla znanego losowego rozwiązania) i błąd rzeczywisty w różnych normach tzn. ||x-\tilde{x}||_2 i ||f-A*\tilde{x}||_2 gdzie f=A*x a \tilde{x} rozwiązanie obliczone przez octave'a czy Państwa algorytm czy w najgorszy sposób A^{-1}*f (wcześniej musimy policzyć macierz odwrotną).
  W skrypcie rozwiązania zadań które udalo sie zrobić na labie: lu18.m
• Lab 4 -kartkówka przy komputerze (lekko zmienione zadanie z poprzednich labów) 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: lznk18.m
• Lab 5
  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. Napisz funkcję liczącą błąd pomiędzy wielomianem interpolacyjnym dla kolejnych stopni dla węzłów równo-odległych oraz węzłów Czebyszewa dla danej funkcji, odcinka [a,b]. Tzn. bardziej precyzyjnie napisz funkcję octave'a (w m-pliku): function err=Lagrangeinterp(FCN,a,b,type=0,M=20), która dla zadanego wskaźnika funkcyjnego FCN (ang. function handle) do funkcji jednego argumentu: function y=f(x)!, a,b - końców odcinka [a,b], typu węzłów 0 - równo-odległych, 1 - Czebyszewa i M max stopnia wielomianu Zwróci macierz która będzie miala w 1szej kolumnie kolejne stopnie wielomianu, w drugiej w j-tym wierszu max|f(x(k)-L_jf(x(k))| dla x(k) j+1 węzłów interpolacji, w trzeciej kolumnie w j-tym wierszu przybliżenie dyskretnej normy maksimum różnicy L_j f - f na tym odcinku - przybliżenie normy liczymy na dyskretnej siatce zawierającej tysiąc punktów
• Lab 6 Kontynuacja interpolacji Lagrange'a i aproksymacja jednostajna i średnokwadratowa
  • Powtórz (zadanie o zbieżności interpolacji Lagrange'a) dla obu typów węzłów, tzn. znajdowanie wielomianów interpolacyjnych na węzłach równo-odległych i węzłach Czebyszewa, ale dla funkcji f(x)=log(1+x) na odcinkach [0,1] i [0,10]. Czy dla tej funkcji i obu odcinków błędy w normach maksimum maleją wraz ze wzrostem N? Porównaj wyniki otrzymane w tym zadaniu w obliczeniach z oszacowaniami teoretycznymi błędu interpolacji Lagrange'a.
  • Interpolacja Lagrange'a przykład Rungego. Powtórz (zadanie o zbieżności interpolacji Lagrange'a) dla obu typów węzłów, tzn. znajdowanie wielomianów interpolacyjnych na węzłach równo-odległych i węzłach Czebyszewa, ale dla funkcji dla f(x)=1/(1+x*x) na [-5,5]. Czy ciąg wielomianów interpolacyjnych zbiega do f jednostajnie?
  • Napisz funkcję znajdującą dla danego wielomianu stopnia n: w(x)=\sum_{k=0}^n a_kx^k oraz liczby q współczynniki wielomianu p(x)=\sum_{k=0}^{n-1}b_kx^k i r takie, że w(x)=(x-q)*p(x)+r, obliczone z wykorzystaniem algorytmu Hornera.
  • Napisz funkcję znajdującą dla danego wielomianu stopnia n w(x)=\sum_{k=0}^n a_kx^k oraz liczby q wartość w(q) i pochodnej w'(q), obliczone z wykorzystaniem algorytmu Hornera.
  • Algorytm Hornera w bazie Newtona. Różnice dzielone. Zaprogramuj w octavie zmodyfikowany algorytm Hornera obliczający wartość wielomianu zadanego w bazie Newtona dla danych węzłów. Parametrami maja być x pkt w którym obliczamy wielomian (ewentualnie tablica punktów ale wtedy funkcja też musi zwrócić tablice wartości wielomianu w tych punktach), N - stopnień wielomianu, wektor długości N+1 ze współczynnikami wielomianu w bazie Newtona. Przetestuj na kilku b prostych przykładach np. 3 węzły -1,0,1 i wielomian w(x)=x^2 - w bazie Newtona związanej z tymi węzłami: x^2=(x+1)x - (x+1) +1 . Zastosuj do przykładów z poprzednich zadań - dla różnych N, narysuj wykresy funkcji i jej wielomianu interpolacyjnego otrzymanego za pomocą tego algorytmu (wartość wielomianu liczymy algorytmem Hornera).
  
  interp18.m - kilka rozwiazan zadan - interpolacja Lagrange'a i Hermite'a (podwójne węzły)
• Lab 7 Ostatni test na labie. Interpolacja Hermite'a Splajny kubiczne. Aproksymacja średniokwadratowa.
  1. 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))
  2. 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).
  3. 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.
  4. Wygeneruj pierwszych wielomianów Czebyszewa na [-1,1] z 1 przy najwyższej potędze jednomianu tzm x^n+... tzn współczynniki w bazie potęgowej korzystając z z reguły trójczłonowej. Narysuj ich wykresy - oblicz przybliżone maksimum i minimum na [-1,1]. Czy teoria zgadza się z teorią.
  5. wylosuj 100 wsp. wielomianów stopnia =5. postaci x^5+\sum_{k<5} a(k)x^k (a(k) losowe) i policz ich przybliżone normy sup na [-1,1] czy są większe niż ||T_5/2^{-4}||_\infty=2^{-4}?
  6. przeczytaj help do funkcji quad() - za jej pomocą sprawdź czy wielomiany Czebyszewa T_n są ortogonalne w iloczynie skalarnym \int_{-1}^1 f*g*1/sqrt(1+x*x) dx.
  7. znajdź wielomian najlepszej aproksymacji stopnia <= 5 w normie generowanej przez iloczyn skalarny \int_{-1}^1 f*g*1/sqrt(1+x*x) dx dla: sin(x), sin(10x), sqrt(|x|), funkcja charakterystyczna odcinka [0,1]
  8. 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ł?
  9. 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().
  10. 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.
  11. 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)|.
  12. 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 jest 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).
  13. 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()
  14. 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?
  
  spline18.m - kilka rozwiazan zadan - interpolacja splajnami kubicznymi
• Lab aproksymacja sredniokwadratowa, wielomiany Czebyszewa, kwadratury.
  1. Napisz funkcję liczącą dyskretny iloczyn skalarny na [a,b] postaci (f,g)=\sum_{k=0}^{N-1} hf(x(k)+0.5*h)g(x(k)+0.5*h) dla h=(b-a)/N (przeskalowana przez h=(b-a)/N iloczyn wartości f i g w środkach odcinków ró?nomiernego podziału [a,b] na N odcinków) - domyslne N=200.
  2. przetestuj rząd zbieżności kwadratury prostokątów P_Nf =\sum_{k=0}^{N-1} hf(x(k)+0.5*h) (przybliżjącej całkę metodą połowienia kroku H (lub podwajania N). Dla f której całkę znamy np sin, cos, x^p etc policz
     Następnie stwórz macierz Gramma G dla bazy jednomianów (1,x,..,x^M) M < < N np. M=20 i policz F(k)=(f,x^k) dla k=0,..,M -dla f - wielomianu stopnia k<=M, sin(x), funkcji charakterystycznej podocinka [a,0.5(a+b)] etp rozwiąż układy Ga=F dla różnych f otrzymując wektory współczynników wielomianu optymalnego w normie generowanej w tym iloczynie skalarnym. narysuj wykres f i odpowiedniego wielomianu (coraz wyższych stopni).
  3. Dla iloczynu skalarnego (f,g)=\int_{-1}^1 (\sqrt(a-x*x))^(-1)f(x)g(x) dx - wygeneruj pierwsze 20 wielomianów ortogonalnych tzn wielomianów Czebyszewa Tn z reguły rekurencyjnej T_{n+1}=2xT_n-T_{n-1} T_1=x T_0=1. Policz (T_n,T_m) dla n,m=0,...,20 z wykorzystaniem funkcji quad() - trzeba przeczytać help i np stworzyć funkcję która z wykorzystaniem polyval liczy T_n(x)T_m(x)/\sqrt(1+x*x) czyli funkcję podcałkową w iloczynie skalarnym.
     Dla funkcji takich jak wielomian st 4, sin(x), f. charakterystyczna [1,0.5(a+b)] itp policz a(n)=(f,T_n) - korzystając z tego znajdź wiel. najlepszej aproksymacji stopnia m - narysuje wykresy f i WNA odpowiedniego stopnia dla kolejnych stopni.
  4. Zapoznaj się z pomocą do funkcji octave'a quad(). Policz przy pomocy quad() całkę \int_a^b f(t) dt dla następujących funkcji i odcinków:
     • x^2 na [-1,1]
     • x^9 na [0,1]
     • sin(x) na [0,\pi]
     • cos(100*x) na [0,\pi]
     • cos(100*x)*cos(1000*x) na [0,\pi]
     Czy wyniki są zgodne z teorią?
  5. Złożona kwadratura trapezów.
     • Zaprogramuj w octavie funkcję function c=kwadtrapez(FCN,a,b,n) obliczającą złożoną kwadraturą trapezów:
       T_n f=h\left(0.5[f(a)+f(b)]+\sum_{k=1}^{n-1} f(a+k*h)\right) h=(b-a)/n
       Parametry funkcji:
       • FCN - wskaźnik do funkcji octave'a function y=f(x) obliczającej wartość funkcji podcałkowej
       • a - lewy koniec odcinka
       • b - prawy koniec odcinka
       • n - ilość obliczeń wartości funkcji podcałkowej w kwadraturze trapezów minus jeden (tak jak we wzorze powyżej)
       Funkcja zwraca przybliżoną wartość całki obliczoną za pomocą powyższego wzoru tzn. złożonej kwadratury trapezów. Funkcja powinna działać również jeśli ją wywołamy tylko z trzema parametrami. Wtedy n powinno domyślnie przyjąć wartość dwieście.
     • Przetestuj dla funkcji f, której całkę znamy dla N_k=2^k*N_0 k=0,1,2,.... z ustalonym N_0=5 licząc:
       E_N/E_{2N}
       dla błędu
       E_N=|\int_a^b f dt - T_N f|
       dla następujących funkcji
       • x^2 na [0,1]
       • sin(x) na [0,pi] i N_0=5 czyli funkcji analitycznej
       • sin(100*x) na [0,pi] i N_0=5 czyli funkcji analitycznej ale silnie oscylującej (duże wartości drugiej pochodnej)
       • x^{j+0.5} na [0,1] dla j=0,1,2 czyli funkcji w C^j([0,1]) i o eksplodującej w zerze j+1 pochodnej.
       Porównaj wyniki obliczane kwadraturą trapezów z wynikami funkcji quad(). (Można sprawdzić dla jakiego n błąd obliczony złożoną kwadraturą trapezów jest na poziomie błędu funkcji octave quad() - i porównać wtedy ilość wywołań funkcji f przez obie procedury).
  6. Złożona kwadratura prostokątów.
     • Zaprogramuj w octave funkcję function c=kwadprost(FCN,,a,b,n) ze złożoną kwadraturę prostokątów
       P_n f=h \sum_{k=1}^n f(a+ (k-0.5)*h)) h=(b-a)/n.
       Parametry funkcji: wskaźnik do funkcji obliczającej FCN, końce odcinków a,b i ilość punktów n. Funkcja zwraca przybliżoną wartość całki obliczoną za pomocą powyższego wzoru tzn. złożonej kwadratury prostokątów. Funkcja powinna działać również jeśli ją wywołamy bez ostatniego parametru. Wtedy n powinno domyślnie przyjąć wartość np. sto lub dwieście.
     • Przetestuj tak samo jak kwadraturę trapezów tzn. dla tych samych funkcji f, której całkę znamy dla N_k=2^kN_0 k=0,1,2,.... z ustalonym N_0=5 licząc: E_N\E_{2N}, E_N=|\int_a^b f - P_N f|
       Porównaj wyniki obliczane kwadraturą prostokątów z wynikami funkcji quad().
  7. Porównaj wyniki obliczane kwadraturą prostokątów z wynikami funkcji obliczającej całkę za pomocą złożonej kwadratury trapezów. Porównaj błędy dla f(x) z poprzedniego zadania czyli całek, których wartość teoretycznie znamy z wynikami dla obu kwadratur:
     • dla tej samej ilości wywołań funkcji f, tzn dla P_N f z T_{N-1} f,
     • porównaj obie kwadratury dla tego samego h tzn. T_Nf z P_Nf aby ocenić błąd w terminach parametru h=(b-a)/N.
  
  integ18.m - kwadratura prostokatow zlozona - testy rzedu zbieznosci

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

• Link do stron Octave'a (skąd można ściągnąć kolejną dystrybucje - pod linuxa czy windows)
• octave-forge - rozszerzenia octave'a
  
• Manual do octave'a w htmlu (po angielsku)
• Link do skryptu do labu z mo z octave'em plik pdf

Skrypty m-pliki octave'a

• mobasic.m - prosty skrypt octave'a w którym są przedstawione niektóre podstawowe operacje macierzowe, oraz pętle, instrukcje warunkowe, zapisywanie/wczytywanie do/z plików itp

  Literatura:

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

  Inne użyteczne linki

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

  ---

  Powrót do mojej strony domowej.

  ---

  Ostatnia aktualizacja: 11 czerwca 2018