Semestr letni 2007/08

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) 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

Program wykladu

• wstęp - czym są obliczenia naukowe, co będzie na wykładzie (0.5 wykładu)
  
• 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!)- 3-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 -1.5- 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 wykłady - a tu książka która zawiera teorie wielowymiarowej metody Newtona: C. T. Kelley, Iterative Methods for Linear and Nonlinear Equations

Program labu/ćwiczeń

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 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.
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.h.  Zad 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)

Skrypty m-pliki octave'a

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

Narzędzia potrzebne do pomiaru czasu