Semestr letni 2006/07

Konsultacje:

Obliczenia naukowe

1. czerwiec 2006(obie grupy; plik pdf )
2. wrzesień 2006 (plik pdf)
3. Czerwiec 2007 - Gr 1 i Gr 2
   a tu spakowane zgzipowanym tarem rozwiazania zadań z bibliotek - obie grupy (możliwie proste rozwiązania)

Zaliczenie:

Program wykladu

• wstęp - czym są obliczenia naukowe, co będzie na wykładzie (1 wykład)
  
• wstęp do Octave'a- podstawowe operacje, m-pliki, operacje wektorowe (macierzowe), skrypty, struktury - pętle for,while do until, instrukcje warunkowe if, switch, zmienne globalne, funkcje jako parametry funkcji, kilka użytecznych funkcji, numeryczny elementarz czyli jak przy pomocy octave'a  rozwiązać podstawowe zadania matematyki obliczeniowej (czy inaczej analizy numerycznej), wykresy funkcji w 2D i 3D (część informacji może być na labie tylko!)- 4 wykłady
  
• krótkie przypomnienie/omówienie met. numerycznych rozwiązywania r.r. zwyczajnych - schematy otwarte, liniowe wielokrokowe vs jednokrokowe (Rungego-Kutty), rząd zbieżności, stabilność, ukł. sztywne czyli laczego należy stosować sch. zamknięte (dla u. sztywnych) - literatura jest szeroka; kilka pozycji można znaleźć:tutaj - 2 wykłady
  
• krótki kurs C (przy załóżeniu że Państwo znacie Pascala) - struktura programu; pliki wsadowe; biblioteki podstawowe (krótko); podstawowe struktury( pętle while(war){};do{..}while(war) ;for(ini;war;kod){..}; instrukcje warunkowe if(war){..}else{..};  switch(wart){case w1:..break;default:}, zmienne, wskaźniki, alokacja pamięci, jak kompilować pod linuxem itp; - 1 -1.5 wykładu
• biblioteki numeryczne - omówienie ogólne bibliotek w tym b. szczegółowo: blas/lapack; w szczególności jak skompilować program używający funkcji/bibliotek fortranowskich w C; 1-2 wykłady
  
• optymalizacja kodu numerycznego architektura pamięci komputera w szczeg. jaki wpływ na szybkośc obliczeń ma dostęp do pamięci cache i jak należy to uwzględnić w C; kilka wskazówek jak optymalizować kod w C - rozwijanie pętli (loops unrolling), wynoszenie stałych elementów poza pętle (hoisting), rozwijanie funkcji (function inlining) opcje kompilatora gcc 1-2 wykłady
• makefile i pomiar czasu procesora ; niektóre użyteczne opcje kompilatora mające wpływ na szybkość programów numerycznych, make - czyli wygodne kompilowanie dużych projektów; jak zmierzyć czas procesora w C, biblioteki statyczne - 0.5-1 wykład
  
• rozwiązywania nieliniowych układów równań w tym wielowymiarowa met. Newtona i ewentualnie krótkie omówienie met. iteracyjnych dla układów równań liniowych - 2-3 wykłady - (niestety tylko 1 wykład się odbył)

Program labu/ćwiczeń

1. (sala)zapoznanie się z środowiskiem i komendami unixa i octave'a typu kopiowanie, usuwanie, przenoszenie, szukanie plików, tworzenie katalogów.
   
2. wstępne zapoznanie się z octavem http://pl.wikibooks.org/wiki/GNU_Octave - octave jako kalkulator, jak korzystać z helpa, tworzenie wektorów - operator ":", operacje macierzowe np tworzenie wektorów i macierzy poprzez wpisywanie, transponowanie tzn ', funkcje tworzące macierze: ones(),zeros(),hilbert(),vander(), jak z macierzy odzyskać podmacierz? (operator ":"!), mnożenie macierzy również element po elemencie (".*") itd, normy wektorów i macierzy, funkcje od macierzy itp np x=[1 2 3 4] i y=sin(x)+ x.*x, zapisywanie i odczytywanie danych (save i load) itp. 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. 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.
3. Kontynuacja. Skrypty - czyli pliki tekstowe z komendami octave'a - przykładowy skrypt matbasic.m wywoluje sie linii komend octave'a: matbasic. Zad 1: korzystając z funkcji vander() (tworzy macierz Vandermonde'a) oraz polyval() - oblicza wartość wielomianu - rozwiązać zadanie regresji liniowej tzn dla danych x_k,y_k znaleźć a,b takie że \sum_k |ax_k+b -y_k|^2 osiąga minimum. Porównać z wynikami funkcji polyfit(). Zad 2 : znaleźć wiel. interpolacyjny dla węzłów rółnoodległych na odcinku [-5,5] dla funkcji f(x)=1(1+x*x) - zbadać dyskretną normę max tej funkcji i wiel interpolacyjnego (testować dla 5,10,20 i 30 węzłów). Zad 3 treść jak w Zad 2 ale dla węzłów Czebyszewa - przy czym ważne wężły Czebyszewa policzyć wektorowo (bez użycia pętli). Węzły Czebyszewa na [-1,1] to n+1 zer wielmianu Czebyszewa zadanego wzorem cos((n+1)arccos(x)) a na dowolnym [a,b] są odpowiednio przeskalowane i przesunięte. Proste wykresy funkcja: plot.
4. Numeryczna klasyka w Octavie: algebra liniowa - układy r. liniowych, rozkłady LU, QR, obliczanie macierzy odwrotnej, liniowe zadania najmniejszych kwadratów, uwarunkowanie macierzy. Macierze rzadkie. Własne funkcje w octavie - m-pliki. Zmienne globalne. Zad 1: utworzyć swój skrypt tworzący macierz 1-wymiarowego Laplacianu (2 na diagonali - -1 na pod- i naddiagonali) różnymi metodami (pętla, pętla lokowo, funkcja diag()) itp).Porównać czas - funkcje tic i toc. Zapisać do pliku w trybie binarnym i tekstowym. Zad 2 Znaleźć macierz odwrotną do macierzy L z zad 1. Rozwiązać układ z macierzą L raz poprzez operator '\' a potem mnożąc przez macierz odwrotną porównać czas. Policzyć uwarunkowanie L. (dla dużych N np N=100,200, 400). Policzyć wartości i wektory własne macierzy L. Zad 3 utworzyć macierz L w formacie sparse (nietworząc wcześniej macierzy w pełnym formacie). Rozwiązać układ równań z tą macierzą - porównać czas z tym z zad 2. Zad 4 Utworzyć funkcję która dla danego wektora węzłów x i wektora tej samej długości y znajdzie wielomian stopnia N taki jak w funkcji polyfit() - funkcja ma poza tym zwrocic dyskretną norme max na 100 punktach na odcinku [min(x(1),..,x(n)),max(x(1),...,x(n)] - n wymiar x i y, jako parametry ma miec x,y,N ale gdy nie bedzie N ma przyjąc N równe 1 (czyli zad regresji liniowej). jesli funkcja bedzie wywolana z 1 zwracana wartościa - wtedy część kodu licząca normę max ma nie być wykonana (oczywiście należy wykorzystać kod z poprzedniego labu). Zad 5Utworzyć funkcję która liczy normę p-tą dyskretną wektora x - dwa parametry - x i p - wywolana tylko z x ma domyślnie przyjąć p = 2.
5. Dokończenie czyli Zad 3 - (Zad 4 i 5 do domu) z poprzedniego labu.
   Rozwiązywanie równań nieliniowych (funkcja fsolve) . Zmienne globalne - moga w pewnych sytuacjach sluzyc do przekazania parametru (global) Zad 1: rozwiązać kilka prostych r. nieliniowych np x*x-4=0; cos(x)=0, sin(x)=1, z różnymi przybliżeniami startowymi; Zmienić tolerancje względną - funkcja: fsolve_options(). Zad 2: przetestować fsolve na trudniejszych przykładach: x*(1+0.5sin(x))=0; (x-1)^2=0; z i bez podawania pochodnej itp Zad 3: zadanie nieliniowe wielowymiarowe np x1+x2=0; x1^2+x2^2-1=0; Zad 4: przy pomocu fsolve narysować wykres funkcji podanej w sposób uwikłany: np wykres y(x) dla y^2+x^2=1 z y(0)=1 albo y(0)=-1; czy y(x) dla f. zadanej wzorem x^3+y^3-y=0 w otoczeniu y(0)=1 Zad 5: narysować wykres funkcji odwrotnej do danej formułą - na prostych przykładach np sin(x) (porównać z arcsin(x)) itp Czy odwrotna do x^2 na [0,1] - porównać z pierwiastkiem...
6. Całkowanie numeryczne (funkcja quad), całki z osobliwościami, Instrukcja warunkowe if else, switch, pętle while, do until, operatory logiczne (&& and, || - or, !-not, != - nie rowna sie, == - równość itp), w szczególności operatory logiczne element po elemencie Zad 1: zaimplementować wektorowo funkcję daszek (f(x)=f(-x);f(x)=0 poza [-1,1], f(x) =1-x na [0,1]) oraz funkcję 1_[0,1/2](x) czyli funkcję charakterystyczną odcinka [0,1/2].
   Zad 2 Policzyć całki z kilku prostych funkcji np sin(x), wielomianów: x^2, x^p itd, całki po odcinku nieograniczonym np \int_{-\infty}^{+\infty} exp(-x^2)dx czy całki z osobliwościami np \int_{-1}^1 \sqrt{|x|}dx czy ogólniej \int_{-1}^1 |x|^pdx dla -1 < p < 1.
   Zad 3 Narysować wykres na odcinku [-20,20] funkcji F(x)=\int_{-\infty}^x exp(-t*t/2) dt (dystrybuanty rozkładu normalnego). Sprawdzić czy lim_{x->+\infty}F(x) wynosi \sqrt(2\pi) a lim _{x->-\infty}F(x)=0. Rozwiązać równanie F(x)=1/2. Sprawdzić wynik. Jak napisac funkcje ktorej jedynym parametrem jest nazwa funkcji postaci y=f(x) liczącej całkową normę L2? (było na wykładzie). Zad 4 Zaimplementować funkcję i napisać wersję zmodyfikowaną liczącą normę L^p (p-drugi parametr). Zad 5: (do domu zapewne) napisać funkcję liczącą całkę postaci \int_0^1f(x)x^kdx dla funkcji postaci y=f(x) a potem korzystając z tego i macierzy Hilberta znaleźć wielomian najlepszej aproksymacji w L²(0,1) dla funckji daszek czy charakterysycnej odcinka [0,1/2] (zad 1). (czyli rozwiązać układ równań z macierzą Gramma w L²(0,1) dla bazy potęgowej którą w tym przypadku jest m. Hilberta i wektorem prawej strony z odpowiednimi iloczynami skalarnymi \int_0^1 f(x)x^k dx. (28 marca 2007)
7. Rozwiązywanie RRzw w octave'ie czyli funkcja lsode. - przypadki skalarne i wielowymiarowe; Zad 1 Rozwiązać kilka prostych równań różniczkowych np y'=-a*y;y(0)=-y na [0,1] albo [-10,0]. Równanie eksplodujące: y'=y*y; y(0)=1 na odcinku [0,1/2] czy [0,3]. Równanie 2-wymiarowe y''=-y z y(0)=1;y'(0)=0 na [0,2*pi](rozwiązanie cos(x)). Zad 2 Narysować wykres funkcji F(s)=y(1,s) zdefiniowanej jako F(s)=y(1,s) gdzie y(t,s) rozwiązanie równania y'=f(y,t) z y(0)=s dla s w [0,1], [-1,1] dla f=y,10y,-10y 0.1*y*y itp Zad 3 Przy pomocy lsode i fsolve znaleźć s0 takie że rozwiązanie rółnania y'=y*(1-y)*sin(y*y+x) z war pocz. y(0)=s0 spełnia y(1,s0)=0.5. Narysować wykres F(s)=y(1,s) - potok fazowy dla t=1 i y(t,s0) dla t w [0,1]. Zad 4 Przy pomocy lsode i fsolve znaleźć rozwiązanie zadania brzegowego: u''+u*u=f(t) u(0)=u(1)=0 dla różnych f np f=0; f=exp(t) itd Narysować wykres rozwiązania, znaleźć min i max rozwiązania. Przy okazji wykresy funkcji w skali logarytmicznej(loglog) i semilogarytmicznej (semilogy, semilogx), opisy wykresów (3 parametr plot) itp Organizacja swoich skryptów: zmienna LOADPATH. (4 kwietnia 2007)
8. Schematy rowiązywania RRZw Zad 1 Zaimplementować w octavie schemat midpoint : x_n+2=x_n+2h f_n+1. Przetestować dla równania y'=-y; y(0)=1; dla różnych h=1/8;1/16;1/32;1/64;1/128 na różnych odcinkach czasu: [0,1],[0,100],[0,1000].
   Sztywność Zad 2Rozwiązać w octavie (lsode) sztywny układ - schematem adamsa i bdf - porównać czas (lsode_options) np dla układu: x'=998*x+1998*y;y'=-999*x-1999*y z war. pocz. x(0)=y(0)=1 na odcinku [0,10], [0,100] itp Sprawdzić czy ten układ sztywny (jakie są wartości własne maciezy [998 1998;-999 -1999]).
   Wykresy funkcji 2-wymiarowych (mesh,surf,contour czy f. pomocnicza meshgrid), eksport rysunkow na drukarke czy do pliku w różnych formatach np ps czy eps (funkcja print - octave-forge). Zad 3 Narysować wykres funkcji sin(x*x+y*y) na [-2,2]x[-1,3] i poziomnice tego wykresu, poziomnice wyeksportować do pliku eps (kolorowego) cont.eps; wyświetlić plik eps przy pomocy ghostview.
   Zad 4 Chaos: przy pomocy lsode rozwiązać układ zaproponowany przez E.N. Lorentza dx/dt = x - y - x^3; dy/dt = x - x*x*y dla różnych war. początkowych.
   Zad 5 Podsumowanie Octave'a . Z wykorzystaniem funkcji lsode rozwiązać równanie przewodnictwa cieplnego paraboliczne d/dtu=d^2/dx^2u, u(-1,t)=u(1,t)=0, u(x,0)=u0(x) zadane. Dyskretyzujemy po x u_xx ~ h^{-2}(u_{n-1} -2u_n +u_{n+1}) otrzymujemy układ r. zwyczajnych du/dt=h^{-2}Au dla A macierzy 3-diagonalnej z -2 na głównej diagonali i 1 na poddiagonalach. - rozwiązujemy lsode dla różnych u0 np sin(\pi*x) czy funkcją z pikiem u0(x)=0 dla x poza [-eps, eps] i u0(x) =(1/eps)((1-x/eps)(1+x/eps))^3 dla eps=0.5,0.1,0.01,0.001. narysować wykres rozwiązanian na [0,1]x[0,3]. (18 kwietnia 2007)
9. Proste programy w jezyku C. Zastosowanie wszystkich podstawowych struktur: petle while, do while, for, instrukcje warunkowe if, switch, pliki wsadowe, wczytywanie danych z klawiatury czy pliku, zapisywanie do pliku czy na ekran (scanf, fscanf,printf, fprintf). Wskazniki, wektory, macierze, kompilowanie, argumenty programu (czyli funkcji main()). Zad 1 zsumowac wyrazy szeregu harmonicznego 1/x^p dla róznych p (używając różnych pętli: do{...}while(); while(){..},for(){} etc). Zad 2 Zmodyfikować program z Zadanie 1 tak aby podawać p jako parametr programu (czyli jako parametr funkcji main()) Zad 3 Pliki w C. Stworzyć plik z wektorem - pierwsza linia wymiar a potem w 2 linii wartości wektora oddzielone spacjami. Następnie wczytać w kolejnym programie ten wektor do programu. Zad 4 Alokacja pamięci. Alokować dynamicznie (i zwalniać pamięć) na wektor typu float,double etc. Napisać funkcje generującą wektor sin(x) dla x wektora z siatką równoodległą na danym odcinku - parametry: końce odcinka, ilość pktów siatki, oraz adres do alokowanego wektora- pod ten adres powinna być zaalokowana pamięć i wpisany odp. wektor Zad 5 Funkcje prekompilatora - #define #if etc Zmodyfikować program np z zad 2 tak aby używając odp opcji kompilatora kompilował się do pojedyńczej lub podwójnej precyzji.
10. Kolokwium Wyniki kolokwium , treści (obie grupy, plik pdf): kol07.pdf a tu przykładowe rozwiązanie (skrypt octave'a): kol07.m
11. Kontynuacja i jak używać funkcji z fortranu w C na podstawie BLAS Level 1
    Dokończyc Zad 4 i 5i z poprzedniego labu. Zad 1: Stworzyć typ strukturalny z 2 polami - wymiar; i pola wektora. Dynamiczna alokacja pamięci (malloc i free) napisać dwie funkcje tworzące wektor i go wymazujące (tzn alokoujące odp. pamięc i potem usuwające zaalokowaną pamięć) . Zad 2: użyć funkcji z biblioteki fortranowskiej w C na przykładzie możliwie prostej funkcji z blasów dnrm2.f -obliczającej normę euklidesową wektora (to przykład z wykładu) i dscal.f -mnoży wektor przez skalar - napisać swój plik wsadowy (w kodzie programu : #include "mojblas.h") zawierający nagłówki obu funkcji
    
12. Macierze w fortranie; LAPACK - rozwiązywanie układów równań liniowych Zad 1: Macierze w C - w formacie fortranowskim czyli jako wektor w którym ulokowane są elementy macierzy kolumnami. Przemnożyć wektor [1 3] przez macierz A=[ 1 -1; 1 1] używając funkcji z BLASów dgemv.f - napisać prototyp i zastosować do powyższego; do domu zastosować do wymnożenia macierzy [2 2; 1 1; 3 -4] przez wektor [1 1]' (notacja jak z octave'a). Zad 2: Rozwiązać układ równań linowych z macierzą A=[ 1 -1; 1 1] i wektorem prawej strony f=[1; 3] (rozw [ 1;1]) przy wykorzystaniem funkcji z lapacka: dgesv.f Zad 3: Rozwiązać układ równań linowych z macierzą B- trójdiagonalną NxN (dyskretnego Laplaciany 1D) z 2 na diagonali i -1 na pod- i naddiagonalą i wektorem prawej strony [2;-1;0;...;0] dla N=50 (rozwiązanie pierwszy wersor) z wykorzystaniem funkcji z lapacka do rozwiązywania układów r. liniowych z macierzą trójdiagonalną, symetryczną, dodatnio określoną: dptsv.f , sprawdzić rozwiązanie (wersor e1[1;0;0;...;0]). Zad domowe: Znaleźć funckje z Lapacka rozwiązująca układ równań liniowych z macierzą trójdiagonalaną i zastosować do rozwiązania układ równań z macierzą która na diagonali ma 7, na naddiagonali ma -2 a na poddiagonali ma -1 z wektorem prawej strony f=[7;-1;0;0;..;0] dla wymiaru N=100 (rozwiązanie to pierwszy wersor) (23/05/2006)
13. Pomiar czasu w C, proste makefile, tworzenie swoich bibliotek statycznych. Zad 1 Załadować pliki: (biblioteka libtimer.a zawierająca min funkcje tic() i toc() działające jak w Octavie) timer.h - plik nagłówkowy do biblioteki libtimer.a, libtimer.a - biblioteka z min funkcjami tic i toc - timer dla kawałków programów w C. Kompilujemy kod źródłowy z plikiem nagł. timer.h - a potem linkujemy z -L. -ltimer na końcu linii kompilacji o ile libtimer.a jest w lokalnym katalogu np. gcc testtimer.c -L. -ltimer Zmierzyć czas działania funkcji dgesv i dtpsv z lapacka z zadań z poprzedniego labu. Zad 2 Utworzyć swój bardzo prosty Makefile - w którym będzie jedna reguła - kompilujący jakiś Państwa program np ten z Zad 1. Zad 3 Utworzyć bardziej skomplikowany Makefile - z tzw reguła wzorcową tzn działającą np dla wszystkich plików o rozszerzeniu .o skompilowanych przez kompilator C - gcc. Zad 4 Utworzyć swoją biblioteke statyczną (ar rc libnazwa.a plik1.o plik2.o .... , a potem ranlib libnazwa.a) i następnie skompilować program z użyciem tej biblioteki np w bibliotece będą pliki obsługująće macierze (alokujące pamięć, zerujące, tworzące np macierz trójdagonalną itp) a w programie wywołanie funkcji dgesv. Zad 5(do domu) Napisać program który przy pomocy funkcji dtpsv rozwiązuje równanie -u''=f u(0)=0, u(1)=0, przy pomocy różnic dzielonych a tzn na siatce x[k]=k*h h=1/N znamy wartości f[k]=f(x[k]) i szukamy wektora u takiego że u([k]) przybliża u(x[k]) rozwiązując równanie z macierzą która ma 2/(h*h) na diagonali a -1/(h*h) na poddiagonali i wektorem prawej strony f=(f(x[1]),..,f(x[N-1])). Wynik tzn wektor u - zapisać do pliku i wyświetlić używając Octave;a (albo wprost gnuplota). Zmierzyć czas rozwiązywania i kompilować przy użyciu make'a.
14. Ostatni lab. Zad 1 Napisać w Octave'ie implementacje wielowymiarowej metody Newtona (klasycznej) w wersji klasycznej (podaje się Jakobian). Przetestować dla równania Ax+ h*h*x.*x=h*h dla h=1/(N+1) i A macierzy trójdiagonalnej wymiaru NxN z 2 na diagonali i -1 na pod- i naddiagonali. Zmierzyć czas. Zad 2 Treść jak Zad 1 ale w C z wykorzystaniem funkcji z Lapacka dptsv.f - Jakobian to macierz trójdiagonalna A+2*h*h*diag(x) (diag(x) macierz diagonalna z z x[k] na pozycji (k,k))! Porównać czas dla dużych N (100,1000 itd)

Kilka przykładowych skryptów, m-plików (plików funkcyjnych ) octave'a, źródeł prostych pogramów w C, źródłowych funkcji z blasów czy LAPACKa czy użytecznych linków
(dodawanych w miarę postępu wykładu/labu)

Skrypty m-pliki octave'a

Funkcje z blasów (w Fortranie)

Funkcje z LAPACKa (w Fortranie)

Pomiar czasu w C i bibliotyki statyczne

Make