Semestr letni 2007/08


Konsultacje:

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

Obliczenia naukowe

(wyklad czw 1215-1345 - sala 4420 i ćwiczenia - 2 grupy - środa 830-10 i 1015-1145 - lab 3041 i sala 4060.


Egzamin II

Termin egzaminu: przy komputerach poniedziałek 1 września 2008 godz 9:00 lab 2042

Wyniki egz II termin ostateczne:
plik pdf
Oceny są zawieszone także i w: USOSie

Proszę wcześniej sprawdzić czy Państwa konta są poprawnie skonfigurowane np czy działa poczta - a przynajmniej jak wysłać z konta maila z załącznikiem, jak pakować tar itp
Zasady podobne jak w I terminie tzn egzamin będzie zawierał zadania z Octave'a i bibliotek (blas/lapack) - być może pomieszane np zrób coś w C wyniki nagraj do pliku zaimportuj do octave'a i coś dalej zrób + kilka pytań teoretycznych (których udzielicie Państwo na kartce). Pierwsza część (max 35min) - egzamin teoretyczny w formie pisemnej (w czasie którego nie wolno korzystać z notatek/ksiązek/komputera) po oddaniu kartki z odpowiedziami - część komputerowa w czasie której można korzystać ze swoich plików oraz z dokumentacji octave'a czy (poprzez np man) funkcji BLAS/LAPACKA ale NIE MOŻNA z programów pocztowych (oczywiście poza końcowym przesłaniem programu do danego pilnującego).
I proszę sprawdzić u mnie po egzaminie czy państwa rozwiązania doszły! (jak nie dojdą to tak jakby Państwo nie brali edziału w egzaminie!
Wszyscy dostaną propozycje ocen. Osoby które uzyskają pewne minimum z części komputerowej będą mogły poprawić zaproponowane oceny na egz ustnym.
Zachęcam do przeczytania i zrobienia zadań z poprzednich lat i I terminu:
Treści zadań z egzaminów z ON z I i II terminu 05/06 i I terminu 07
  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) on07egz1lib.tgz
  4. czerwiec 2008(plik pdf)
    a tu spakowane zgzipowanym tarem rozwiazania zadań z octave'a i bibliotek - obie grupy (możliwie proste rozwiązania) - jeden Makefile : make zad? (?=1 lub 2) skompiluje odpowiednie pliki dla odp grupy on08rozw.tgz

Egzamin I termin

W pliku pliku pdf są wyniki egzaminu I termin - ostateczne!
Oceny są zawieszone także i w: USOSie

Kolokwium przy komputerze:
Możliwe rozwiązania zadań z gr o 10: kol08-10.m i z gr o 830: kol08-08.m

Program ćwiczeń czyli głównie labu

Program najbliższego labu

Kilka skryptów Octave'a czy przykładowych plików C

Program wykladu

(kolejność i zakres punktów; może ulec znaczącej zmianie)

Ćwiczenia/Lab

Większość zajęć (90%) w labie 3041 - rzadko w sali 4060.

Program labu/ćwiczeń

(kolejne punkty dodawane przed odp. zajęciami)
  1. wstępne ćwiczenia w sali 4060 - krótkie omówienie podstaw Unixa - tzn podstwowych narzędzi wywoływanych w terminalu typu: np ls - listowanie zawartości katalogu itd, prostych edytorów tekstowych, pakowaczy, mailerów itp (20/02/2008)
  2. Lab - podstawy octave'a - 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(27/02/2008)
  3. Lab - 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ć). (5/03/2008)
  4. Proste wykresy plot. Zadanie 1: narysuj na 1 wykresie funkcje i wielomiany interpolacyjne z poprzedniego labu (przykład Rungego dla N=5,10,20). 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 (sparse). Własne funkcje w octavie - m-pliki. Zmienne globalne. Zad 2: 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()) itp).Porównać czas - funkcje tic i toc. Zapisać do pliku w trybie binarnym i tekstowym. Zmienić skrypt na m-plik tworzący L - czyli funkcje octave'a z parametrem N - wymiarem L z odp. helpem. Zad 3 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 4 utworzyć macierz L w formacie sparse (nie tworzą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. (12/03/2008)
  5. Pliki funkcyjne (m-pliki, m-files) czyli jak definiować funkcje. Help do funkcji, ilość zwracanych wartości czy ilość argumentów, domyślne wartości argumentów, zmienne określające ilość argumentów i wartośći funkcji (nargin, nargout) itp Rozwiązywanei ukł. r. nieliniowych - funkcje fsolve. Zad 1 utworzyć swój m-plik zawierający funkcje tworzącą macierz 1-wymiarowego Laplacianu (2 na diagonali - -1 na pod- i naddiagonali) w formacie sparse z parametrem N - wymiarem L z odp. helpem (można wykorzystać zad 4 z poprzedniego labu!). Dodać opcjonalny parametr tekstowy - napis full ma zmienić fromat zwracanej macierzy na pełny, a domyślny ma być sparse. Zad 2: 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 3 Utworzyć 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.
  6. 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... Zad 4: narysować wykres funkcji odwrotnej do danej na siatce - na prostych przykładach np sin(x) na 100 pktach na [0,pi/2] (porównać z arcsin(x)) itp Czy odwrotna do x^2 na [0,1] - porównać z pierwiastkiem... Czy funkcji dla której nie znamy explicite wzoru na odwrotną np. (x+1)*exp(x^2) na [1,3] To samo ale 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] - f([0;0])=[0;0] Ważne! Odp dobór punktu startwego pozwoli na szybsze obliczenia w zad 3 i 4!
    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 4: 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 5: (pętla for, while) policzyć na odcinku [0,pi] na siatce z krokiem h=pi/100 elementowej - przybliżoną pochodną sin(x) (df(x) przybliżamy przez (f(x+h)-f(x))/h - przy pomocy pętli for czy while oraz wektorowo. Porównać błąd z cos(x) na tej siatce w normie sup - następnie policzyć pochodną przy pomocy ilorazu centralnego f(x+h)-f(x-h)/2h i policzyć błąd w normie max na tej siatce.
  7. Całkowanie numeryczne (funkcja quad), całki z osobliwościami,
    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? (trudne - trzeba użyć funkcji octave'a: eval() do def funckji w funkcji). 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 L2(0,1) dla funckji daszek czy charakterysycnej odcinka [0,1/2] (zad 1). (czyli rozwiązać układ równań z macierzą Gramma w L2(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.
  8. 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: path,addpath itp. (16 kwietnia 2008)
  9. Proste programy w C: 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. (23 kwietnia 2008)
  10. Proste programy w C i BLASy: Zad 1 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 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.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 (30 kwietnia 2008)
  11. Kolokwium - 7 maja 2008, przy komputerze, można używać helpa i swoje skrypty. Na koniec prześlecie Państwo spakowane skrypty/pliki C mailem do mnie L.Marcinkowski at mimuw.edu.pl. A tu możliwe rozwiązania zadań z gr o 10: kol08-10.m i z gr o 830: kol08-08.m
  12. Biblioteki - macierze w fortranie i lapack. Zad 1 Zad 3 z poprzedniego labu czyli 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). 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.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)
  13. Pomiar czasu w C. - jak zmierzyć czas procesora w C? - funkcja times() i struktura struct tms (opisana w pliku nagłówkowym sys/times.h, nas interesuje pole tms_utime. Ponieważ funkcja times zwraca w polu struktury tms_utime ilość taktów zegara potrzebujemy sysconf(_SC_CLK_TCK) ktora zwraca ilośc taktów zegara na sekundę no i pliku nagłówkowego do tej funkcji unistd.hZad 1: napisać funkcje tic i toc jak w Octavie z wykorzystaniem funkcji times() - warto obejrzeć plik CPUtimer.c - w make.tar.gz. Zad 2: Zmierzyć czas rozwiązywania układu równań z poprzednich ćwiczeń dla macierzy trójdiagonalnej - porównać z Octave'm.Biblioteka statyczna Jak stworzyć swoją biblioteke (statyczną) a potem jak ją używać. Zad 3: stworzyć bibliotekę zawierającą Państwa funkcje np. tic toc czy funkcje dotyczące wektorów/macierzy np funkcje tworzące wektor jako strukturę o 2 polach : wymiar i dane, dodające 2 wektory, mnożącą przez skalar itp nastepnie skompilować i uruchomić jakiś prosty program używający funkcji z bibliotek np tic() i toc(). Zad 4 Przy pomocy odp funkcji lapacka (DSTEV i DGEEV) znaleźć wartości własne macierzy trójdiagonalnej z -1 na pod- i nad-diagonali i 2 na diagonali dla wymiarów: N=10, 100, 1000. Zmierzyć czas - która szybsza? porównać z octave'm dla N=1000.
  14. Make czyli wygodny sposób na kompilacjewhg reguł z Makefile'a - mozliwie prosty Makefile:Makefile-simple. Zad 1 Napisać prosty makefile który kompiluje jakiś Państwa starszy program z labu - najlepiej program w kilku plikach. Zad 2(do domu po jutrzejszym wykładzie) Napisać b skomplikowany Makefile zawierający regułe wzorcową która pasującą do wszystkich plików %.c np z opcją -O0 (brak optymalizacji) Oraz zmodyfikować regułe domyślną (implicite rule) dla plików wykonywalnych modyfikując odpowiednie zmienne make'a w Makefilu. Zad 3(zadanie ze starego egzaminu) Napisać funkcję której parametrami wejściowymi będą macierz A MxN, wektor b dl M, oraz wymiary M,N która zwraca jako parametr wyjściowy rozwiązanie x - wektor długości N problemu Ax=B w sensie LZNK przy czym jeśli M=N i A nieosobliwa te rozwiązanie jest obliczone funckją DGESV() a jeśli M<>N lub A numerycznie osobliwa (tzn DGESV() nie da rady rozwiązać) to przy pomocy funkcji DGELSS() (funkcja ma A i b pozostawić niezmienione). Proszę napisać odp. Makefile. Przetestować dla A=[1 1;-1 1], [1 1; 2 2] dla b=[2;0] i dla A=[1 2; 1 1; 1 0] i b=[1;1;2]. Sprawdzić w octavie czy wyszły dobre rozwiązania! Zmierzyć czas - porównując z octave'm. Zad 4(zadanie ze starego egzaminu) Napisać funkcję w octavie która razwiąże równanie Poissona (-\Laplacian u= f w \Omega i u=0 na brzegu) na kwadracie [0,1]x[0,1]. Wprowadzamy siatke NxN : (i*h,j*h) i,j=1,..,N z h=1/(N+1). I szukamy funkcji u zdefiniowanej na tej siatce spelniającej zdyskretyzowane równanie na siatce:
    h^{-2}[-u((i-1)*h,j*h) - u(i*h,(j-1)*h) +4* u(i*h,j*h) -u((i+1)*h,j*h) - u(ih,(j+1)*h)]=f(i*h,j*h) dla i,j=1,..,N.
    Oczywiście uwzględniamy to że u(0,j*h)=u(i*h,0)=u((N+1)*h,j*h)=u(i*h,(N+1)*h)=0 czyli warunek brzegowy. Na wejściu N i macierz NxN o wartościach f(i*h,j*h) na (i,j)tej pozycji albo nazwa funkcji 2-zmiennych obliczjącej wartość f, funkcja ma zwracać wektor wartości u(i*h,j*h) jako macierz U czyli u(i*h,j*h) ma być na pozycji (i,j)tej w macierzy U czyli U_{i,j}. Do domu: ogólniejsze zadanie w Octavie ale z niezerowymi warunkami brzegowymi (dane wartości funkcji na brzegu przenosimy do wektora prawej strony) oraz analogiczne zadanie w C - układ równań liniowych można rozwiązać przy pomocy odp. funkcji lapacka np dgesv() czy lepiej dposv() (dla symetrycznych dodatnio określonych). Następnie wyeksportować od pliku tekstowego wektor u jako macierz czyli u(i*h,j*h) na pozycji (i,j)-tej w macierzy i wyświetlić tę macierz w octavie.
  15. Wielowymiarowa metoda Newtona w Octavie i C- Zad 1 (pomocnicze) Napisz procedure która dla danej funkcji 1 argumentu (wektorowego) zwróci macierz przybliżonego Jakobianu D_hF(x) przybliżającego DF(x) przez ilorazy różnicowe (tzn A_{ij}=h^{-1}(F_i(x+ h*e_j) - F_i(x)) \approx D_j F_i(x)) - czyli parametry h=1e-7, x- wektor i nazwa funkcji postaci y=F(x) która dla wektora z R^N zwróci wektor z R^N. Przetestować dla F(x)=A*x-b (czyli wtedy DF=A jakaś dowolna macierz 2x2) Zad 2 Zaimplementować metoda Newtona z przybliżonym Jakobianem - tzn funkcję z dwoma parametrami wymaganymi: nazwą funkcji postaci y=F(x) oraz x0 przybliżenie startowym. Pochodna obliczona ma być w sposób przybliżony jak w zad 1 (za parametr h można wziąść 1e-7) a za warunek stopu ||F(x_n)||< tola + rtol*||F(x0)|| - za atol można wziąść domyślnie 1e-8 za rtol - 1e-5. Jeśli przekroczona zostanie ilość max iteracji max_it (powiedzmy równe 30) funkcja ma zwracać ostatnią iterację i info=-1, gdy metoda zbiegnie tj warunek stopu zadziała funkcja ma zwrócić obliczone x_n oraz w info=0. Można dodać opcjonalny trzeci parametr - wektor z tolerancjami (2 elementowy). przetestować dla F(x,y)=(x-y,x*x+y*y-1) startując np z (1,1). Zad 3 (do domu) Zmodyfikować funkcje powyżej - tak aby jakobian był obliczany co ustaloną ilość iteracji (np co 2 iteracje) czyli powiedzmy w pierwszej iteracji obliczamy jakobian i dalej kolejną iteracje Newtonowską a potem w drugiej nie obliczamy jakobianu tylko używamy ten z poprzedniej iteracji. Do zastanowienia się - czy można oszczędzić jakoś na koszcie rozwiązywania układu równań liniowych? (wiemy że octave używa rozkładu LU.... czy rozkład Lu z iteracji pierwszej możmy jakoś użyć w iteracji drugiej? (czy potem rozkład z iteracji 2n-1 w 2n-tej?)) Zad 4 Zaimplementować met z zad w C - używając funkcji DGESV z lapacka do rozwiązywania układu równań liniowych.

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)

Tutaj link do stron Octave'a (skąd można ściągnąć kolejną dystrybucje - pod linuxa czy windows)
octave-forge - rozszerzenia octave'a

Skrypty m-pliki octave'a

matbasic.m - przykładowy skrypt z tworzeniem macierzy, rozwiązywaniem ukł .r. liniowych itd Wywołujemy komendą matbasic (z octave'a)- modyfikujemy dowolnym edytorem tekstowym np wywołując gedit matbasic.m (z lini komend unixa npz terminala)
sesja3.m - przyklad jak rozwiązać ukl r liniowych w Octavie min zastosowanie octave'a do znalezienia wiel. interpolacyjnego czy wielomianu pasującego do zadanych p-tów (f. polyfit albo przy pomocy narzędzi algebry liniowej)
sesja4.m - przyklady na proste całkowanie i rozwiązywanie równania nieliniowych w Octavie, operatory logiczne itp
sesja5.m - rozwiązywanie RRzw - f. lsode() - dwa proste przykłady
approx.m - proste zastosowanie - znalezienie wiel. najlepszej aproksymacji w normie L^2(0,1) (calkowej); m.in. przykład jak definiować funkcje w funkcji. Zadanie: zmodyfikować przykład dla L^2(a,b) czy dla innej miary np z waga Gaussowska exp(-x^2/2) na R - odrobinę trudniejsze
zagbrzeg.m - jak można rozwiązać zagadnienie brzegowe w octavie

Przykłady z C (C++) w tym funkcje BLASów i LAPACKa

plik.c - możliwie prosty przykład na zapis i odczyt do/z pliku w trybie tekstowym i bianrnym, program zapisuje macierz liczb o podwójnej precyzji (double) trzymaną wierszami w jednej tablicy tekstowo do pliku a następnie ją wczytuje do innego fragmentu pamięci (innej tablicy) a potem to samo do/z pliku binarnego!
dnrm2.f - możliwie prosta - oblicza normę drugą wektora
ddot.f - możliwie prosta - oblicza iloczyn skalarny (standardowy) dwu wektorów

blas-simplest.c - bardzo prosty program w C z wykorzystaniem dwu prostych funkcji fortranowskich z BLASow, kompilujemy: gcc blas-simplest.c -lg2c -lblas
blasbasic.cc - przyklad w C++ na użycie funkcji z blasów - DNRM2 - kompilujemy np g++ blasbasic.cc -lblas
macwek.c - przykładowy prosty program na użycie funkcji dgemv z BLAS level 2 czyli jednym z parametrów jest macierz - obliczającej y=alpha*A*x+beta*y dla alpha,beta -skalarów, A macierzy MxN, X wektora długości N i y - wektora długości M
macwek.cc - a tu wersja w C++
dgesv.f - funkcja z Lapacka - rozwiązywanie układów równań liniowych z macierzą ogólną (pełną wg formatu Fortranu tzn macierz jako wektor - umieszczony kolumnami)
dptsv.f - funkcja z Lapacka - rozwiązywanie układów równań liniowych z macierzą trójdiagonalną, symetryczną dodatnio określoną, format macierzy - dwa wektory - diagonala i poddiagonala (równa naddiagonali)

trisv.c - przykładowy prosty program na użycie funkcji dptsv z LAPACKa rozwiązującej układ równań liniowych A*X=B z macierzą A trójdiagonalną dodatnio określoną (zadaną przez diagonale i poddiagonale)

Narzędzia potrzebne do pomiaru czasu

times() funkcja opisana w pliku nagłówkowym sys/times.h - zwraca w struct tms - czas użytkownika w taktach (czy milionach taktów) zegaru
sysconf(_SC_CLK_TCK) - zwraca ilośc taktów zegara na sekundę (czy raczej 1000 czy milionów taktów wbrew opisom funkcji ale ważne że zgodnie z tym co zwraca funkcja times().
unistd.h. - pliku nagłówkowy do funkcji sysconf

Wyklad - Środowko programisty gdzie można znależć materiały o języku C jak i Makefile'ach

Powrót do mojej strony domowej


Ostatnia aktualizacja: 1 września 2008

Dziś jest