Obliczenia naukowe 2015/16

wtorek wyklad 1015 w 5870, lab 1215 w 2041(lab)

Konsultacje:

Patrz: Plan. (pokój 5010 - IV piętro - wieża północna przy schodach). w czasie sesji czy ferii tylko po uprzednim kontakcie.

Egzamin i ocena

Ocene planuje na bazie projektów. prosze o kontakt mailowy aby umowic spotkanie. w czasie egzaminu trzeba zaprezentowac dzialanie programu, kilka rozsadnych testow (np. dla znanego rozwiazania zerowego - o ile takie jest) oraz odpowiedziec na kilka prostych pytan dotyczacych materialu z wykladu. (za projekt 80% za znajomosc wykladu 20%).

W polowie wrz bede w czwartek 15/09/2016 od ok 12. Byc moze bede i pn-sr 12-14/09/2016 - jesli tak to informacje umieszcze w weekend.

Program najbliższego labu

Kilka skryptów Octave'a czy przykładowych plików C - beda sie pojawiac stopniowo

Projekt

Literatura

P. Krzyzanowski, Obliczenia inzynierskie i naukowe. Skuteczne, szybkie, efektowne, Wydawnictwo Naukowe PWN, 2011
P. Krzyzanowski, Obliczenia naukowe. Skrypt w html. Wersja w pliku pdf dostepna tamze. Uniwersytet Warszawski. 2010.

Ćwiczenia/Lab

Program labu/ćwiczeń

(kolejne punkty dodawane przed odp. zajęciami)

Lab 1 - podstawy octave'a przypomnienie - 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, ewent. rozwiązywanie układów równań liniowych, LZNK (liniowe zad najmniejszych kwadratów) itd Zad 1 Utwórz dowolne macierze 3x4 A i 3x5 B - a następnie macierz 3x8 C której pierwsze 3 kolumny to A a kolejne to B. Teraz z tej macierzy 'wytnij' podmacierz D składającą się z 1 głównego minora tzn 3x3 od C(1,1) do C(3,3). Zamień kolejność kolumn D. Wstaw D z powrotem do C jako główny minor. Policz sin od D. Zapisz D do pliku (binarnego i ASCII) - zamień element D(1,1) na -100 i wczytaj nową macierz do octave'a jako DD. Policz normy (różne) macierzy DD. Zad 2 Policz dyskretną normę max od (sin(x))^2 na [0,1] (wektorowo)
matbasic.m
Lab 2- podstawy octave'a cd - . Skrypty - czyli pliki tekstowe z komendami octave'a - przykładowy skrypt matbasic.m wywoluje sie linii komend octave'a: matbasic. Zadanie 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. Napisać swój skrypt z rozwiązaniem tego zadania i wywołać go z linii komend octave'a. Zbadać uwarunkowanie macierzy Vandermonde'a. Porównać z wynikami funkcji polyfit(). Zad2: znaleźć wiel. interpolacyjny dla węzłów równoodległych na odcinku [-5,5] dla funkcji f(x)=1(1+x*x) - zbadać dyskretną normę max tej funkcji i wiel interpolacyjnego oraz normę różnicy funkcji i wielomianu interpolacyjnego. Czy norma różnicy maleje wraz ze wzrostem ilości węzłów?(testować dla 5,10,20 i 30 węzłów). Zad 3- treść jak zad 2 ale dla węzłów Czebyszewa (przeskalowanych zer wielomianu T_n(x)=cos (n*arc cos(x)) - można je łatwo wyliczyć).
lab2.m - kilka rozwiazan
Lab 3 i 4 Rozwiązywanie równań nieliniowych (funkcja fsolve) c.d.. Zmienne globalne - moga w pewnych sytuacjach sluzyc do przekazania parametru (global) Zad 1: Przetestować fsolve na trudniejszych przykładach: x*(1+0.5sin(x))=0; (x-1)^2=0; Zad 2: 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 3: narysować wykres funkcji odwrotnej do danej - 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... dla funkcji dla której nie znamy explicite wzoru na odwrotną np. (x+1)*exp(x^2) na [1,3] czy x+cos(x). Zad 4: narysować wykres funkcji odwrotnej do danej na siatce w 2 wymiarach znaleźć np na siatce 10x10 na [0,1]x[0,1] przybliżenie f odwrotnej do f([x;y])= [10x+y*y; 10y+x*y]
Powazniejsza rzecz - rozdzial 3 w skrypcie P.K. O czyli funkcja odwrotna na troche powazniej. Ważne. Odp dobór punktu startowego pozwoli na szybsze obliczenia w zad 2, 3 i 4.
nonlin16.m - kilka prostych przykladow wywolan solverów octave'a dla r. nieliniowych
Proste wykresy plot. Numeryczna klasyka. Calkowanie.
- Zad: narysuj na 1 wykresie funkcje i wielomiany interpolacyjne z poprzedniego labu (przykład Rungego dla N=5,10,20).
- Numeryczna klasyka w Octavie cd: algebra liniowa - układy r. liniowych, rozkłady LU, QR, obliczanie macierzy odwrotnej, liniowe zadania najmniejszych kwadratów, uwarunkowanie macierzy. Macierze rzadkie (sparse). Zmienne globalne, statyczne (persistent a=0), .
- Zad 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 blokowo, funkcja diag(), sparse) itp). Porównać czas - funkcje tic i toc.
- Zad Znaleźć macierz odwrotną do macierzy L z zad 1. Rozwiązać układ z macierzą L raz poprzez operator '\' dla L w formacie sparse i full, a potem mnożąc przez macierz odwrotną porównać czas.
  on16sparse.m - kilka rozwiazan dla macierzy sparse
- operatory logiczne (&& and, || - or, !-not, != - nie rowna sie, == - równość itp), w szczególności operatory logiczne element po elemencie
- Zad: 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].
- Całkowanie numeryczne (funkcja quad(), quadv() etc trapz(), dblquad(), triplequad()), całki z osobliwościami
- Zad 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 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 dazy do -\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ę L^2?
- Zad Zaimplementować funkcję i napisać wersję zmodyfikowaną liczącą normę L^p (p-drugi parametr).
- Zad : 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.
- Zad policz calki 2wymiarowe: \int_0^1\int_0^1 sin(x*y) dx dy lub \int_([0,1]^2) sin(x+y^2)*(x-3y) dx dy lub pole elipsy: 2x^2+3y^2-1=0 (y(x) jest dana jako f uwiklana... )
Poniewaz nie bylo ostaniego labu - rozwiazemy z 1-2 zadania z poprzedniego na calkowanie. 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 < span style="font-weight: bold;">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: path,addpath itp.
- Zad 5 Znajdz s takie ze f(s)=\int_c^b|y(t;s)-h(t)|^2 dt minimalne dla y rozwiazania y''=F(t,y,y'); y(a)=0;y'(a)=s; h- zadana funkcja na [c,b] (a < c < b) i F dane pole wektorowe - dana funkcja 3 argumentów e.g F=-sin(x) (penduulum). Mozna zastapic zad poczatkowe przez brzegowe : np niech y rozwiazanie y''=F(t,y,y'); y(0)=0;y'(b)=s; etc
Kontynuacja rozwiazywania RRZw - model trzesienia ziemim transformator (rozdzial 6 w skrypcie PK), minimalizacja zadania zadanego jakozad poczatkowe lub brzegowe (w 1 wymiarze) etc szczególy na stronie po angielsku.
Zadanie optymalnej kontroli - kontrolowanego przez r. eliptyczne w 1 wymiarze i zadanie paraboliczne w 1 wymiarze - kontrolujemy lewy warunek brzegowy (np war Dirichelta) i szukamy takiej wartości by zadane rozwiazanei było możliwie bliskie zadanej funkcji na kawałku odcinka w sensie normy L^2 - szczegóły na stronie po angielsku.
Proste programy w C i BLASy: Zad 1 Napisać program który policzy sumę szeregu harmonicznego \sum 1/n^p dla p>1 przy wykorzystaniu różnych pętli for(){ }, while(){ }, do{ }while(), ewentualnie w wersji z p podanym z linii komend. Zad 2 Pliki tekstowe i binarne w C. Napisać program, który zapisze wektor (tablica double) do pliku tekstowego lub do pliku binarnego a następnie drugi program który wczyta ten wektor i wyświetli na ekranie (to czy plik tekstowy czybinarny podajecie państwo z klawiatury). Przetestować dla wektora [1 2 3 4 -9.45]. Do domu: program który wczytuje macierze w formacie octave'a do programu w C - macierz np jako podwójna tablica double odp alokując dynamicznie pamięć i wyprowadza wektor na ekran po czym zwalnia pamięć. Przetestować eksportując jakiś wektor z octave'a i wczytując do programu. Napisać program który exportuje do pliku tekstowego w formacie Octave'a macierz z C. Przestestować na macierzy Hilberta A=(1/(i+j+1))_{i,j=0,..,N} - wyeksportować i wczytać do sesji octave'a. Zad 1 Alokacja pamięci. Alokować dynamicznie (i zwalniać pamięć) na wektor typu float,doubleetc. 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 a następnie (ewent dw domu) wyeksportować wektor sin(x) do pliku binarnego (fwrite()) i potem wczytać do programu ponownie pod inny adres. Zad 2 Funkcje prekompilatora - #define #if etc Zmodyfikować program np z zad 1 tak aby używając odpowiednich opcji kompilatora kompilował się do pojedyńczej lub podwójnej precyzji. Albo zmieniając pierwszą linie kodu albo podając odpowiednią opcje prekompilatora z linii komend? Zad 3 napisać makro (jak na wykładzie) do indeksowania macierzy trzymanej kolumnami w w jednym wektorze, oraz funkcje które z wykorzystaniem tego makra wypisuje macierz na ekran, mnoży wektor przez macierz, dodaje 2 macierze, oblicza normę 2 wektora itp. Zad 4(pewnie do domu) stworzyć nowy typ strukturalny - pierwsze pole wymiar wektora drugie - dane - pola wektora jako wskaźnik do double - oraz funkcje tworzące wektor (alokujące pamięć do drugiego pola zmiennej strukturalnej) span style="font-weight: bold;">. Zad 5: użyć funkcji z biblioteki fortranowskiej w C na przykładzie możliwie prostej funkcji z blasów DNRM2(): 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; kompilujemy: gcc naszprogram.c -lgfortran -lblas
tblas.c - prosty program w C z funkcja BLASów dnrm2() i dynamiczna alokacja pamieci na wektor (tzn wykorzystujemy funkcje z libstdlib.a: malloc(); free())
Biblioteki - macierze w fortranie i lapack.
1. Zad 1 makro IJ(i,j,M) zwracające nr indeksu elementu macierzy trzymanej kolumnami o długości M w wektorze w C. Przy pomocy odpowiedniej funkcji a BLASów dgemv.f przemnóż w programie C macierz [1 2;-1 1] przez wektory [0;1], [1;0], [1;1]. Wynik wyprowadź na ekran. Do domu zastosować do wymnożenia macierzy [2 2; 1 1; 3 -4] przez wektor [1;-1] (notacja jak z octave'a).
2. Zad 2: Rozwiązać układ równań liniowych z macierzą A=[ 2 -2; 1 1] i wektorem prawej strony f=[0; 2] (rozw [ 1;1]) przy wykorzystaniem funkcji z lapacka DGESV(): dgesv.f
3. Zad 3: Rozwiązać układ równań liniowych 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(): dptsv.f , sprawdzić rozwiązanie (wersor e1[1;0;0;...;0]).
4. Zad: 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)
5. Zad: regularne LZNK: Uzywajac DGELS(): dgels.f rozwiaz RLZNK: Ax=f. Testuj dla A=[1,1;1,2;1,3] i f=[2;3;4] (x=[1;1]) Rozwiaz Aa=y zadanie regresji liniowej, i.e. x,y dane wektory, szukamy a=arg\min_a \sum_k|a(1)+a(2)x(k)-y(k)|^2. Napisz funkcje C/C++ rozwiazujaca liniowa regresje z uzyciem DGELS(). Testuj dla x=[1;...;10] i y=x+1 dla M=2,3,5,0,120. (a=[1;1]). Funkcja za parametry ma miec: (INPUT) integer M ilosc rownan,(wierzy macierzy, wymiar x,y), wskazniki do double: *x,*y, (OUTPUT) wskaznik do double a - rozwiazanie LZNK, wskaznik do integer info - kod wyjscia funkcji - 0 OK, 1,2,etc bledy