MNwyklad01

Wprowadzenie do metod numerycznych

Metody numeryczne to dziedzina wiedzy zajmująca się problemami obliczeniowymi i konstrukcją algorytmów rozwiązywania zadań matematycznych. Najczęściej, zadania obliczeniowe postawione są w dziedzinie rzeczywistej (lub zespolonej) i dlatego mówimy o zadaniach obliczeniowych matematyki ciągłej (w odróżnieniu od matematyki dyskretnej).

Zadania metod numerycznych

Aby w ogóle mówić w problemie obliczeniowym, musimy najpierw

• określić dane problemu i cel obliczeń, czyli dokładnie

sformułować zadanie w języku matematyki,

• określić środki obliczeniowe dzięki którym chcemy osiągnąć cel,
• dla analizy zadania i sposobów jego rozwiązania wygodnie jest zdefiniować

klasę rozpatrywanych danych oraz model obliczeniowy w obrębie którego będą działać nasze algorytmy.

Wbrew dość powszechnej opinii nie jest prawdą, że głównym celem metod numerycznych jest badanie wpływu błędów zaokrągleń na wynik. Raczej, głównym celem metod numerycznych jest konstrukcja optymalnych (w jasno określonym sensie, np. pod względem wymaganej liczby operacji, lub pod względem ilości niezbędnej informacji, czy też pod względem dokładności uzyskiwanego wyniku) algorytmów rozwiązywania konkretnych zadań matematycznych.

Uwaga

Nasz przedmiot ma różne wcielenia i z tego powodu czasem nosi inne nazwy, w zależności od tego, na jaki aspekt metod obliczeniowych jest położony największy nacisk.

metody numeryczne

--- główny nacisk idzie na aspekty algorytmiczne, ale przy uwzględnieniu pozostałych;

analiza numeryczna

--- przede wszystkim badanie właściwości algorytmów, ich optymalności oraz wpływu arytmetyki zmiennopozycyjnej na jakość uzyskanych wyników;

matematyka obliczeniowa

--- głównie teoretyczna analiza możliwości taniej i dokładnej aproksymacji rozwiązań zadań matematycznych;

obliczenia naukowe

--- nacisk na praktyczne zastosowania metod numerycznych, symulacje, realizacje na komputerach o dużej mocy obliczeniowej.

Oczywiście, granice podziału nie są ostre i najczęściej typowy wykład z tego przedmiotu stara się pokazać pełne spektrum zagadnień z nim związanych. Tak będzie również i w naszym przypadku.

Model obliczeniowy

Tworząc i analizując algorytmy, jakie będą pojawiać w naszym wykładzie, będziemy posługiwać się pewnym uproszczonym modelem obliczeń, dzięki czemu będziemy mogli skoncentrować się na esencji algorytmu, bez detali implementacyjnych --- zostawiając je na inne okazje.

Aby zdefiniować nasz model obliczeniowy, posłużymy się pojęciem programu. Zastosujemy przy tym notację podobną do tej z języka programowania C.

Program składa się z deklaracji, czyli opisu obiektów, których będziemy używać, oraz z poleceń (instrukcji), czyli opisu akcji, które będziemy wykonywać.

Dostępnymi obiektami są stałe i zmienne typu całkowitego (int), rzeczywistego (float i double). Typ logiczny symulujemy tak jak w C wartościami zero-jedynkowymi typu całkowitego.

Zmienne jednego typu mogą być grupowane w wektory albo tablice.

Widzimy więc, że podstawowe algorytmy numeryczne będą bazować na mało skomplikowanych typach danych. Również nieskomplikowane będą instrukcje naszego modelowego języka.

Podstawienie

 
  z = <math>\displaystyle {\cal W}</math>;

gdzie $z$ jest zmienną, a ${𝒲}$ jest wyrażeniem o wartościach tego samego typu co $z$. Jest to polecenie proste.

Wyrażeniem jest pojedyncza stała lub zmienna, albo złożenie skończonej liczby operacji elementarnych na wyrażeniach. Operacje elementarne to:

arytmetyczno--arytmetyczne

$x\mapsto -x$, $(x,y)\mapsto x+y$,

$(x,y)\mapsto x-y$, $(x,y)\mapsto x*y$, $(x,y)\mapsto x//y,y\ne 0$, gdzie $x,y$ są stałymi lub zmiennymi liczbowymi,

arytmetyczno--logiczne

$(x,y)\mapsto x<y$,

$(x,y)\mapsto x\le y$, $(x,y)\mapsto x=y$, $(x,y)\mapsto x\ne y$, gdzie $x,y$ są stałymi lub zmiennymi liczbowymi,

logiczno--logiczne

$p\mapsto \stackrel{}{~}p$,

$(p,q)\mapsto pandq$, $(p,q)\mapsto porq$, gdzie $p,q$ są stałymi lub zmiennymi logicznymi.

Dla niektórych zadań wygodnie jest (a czasem koniecznie) uzupełnić ten zbiór o dodatkowe operacje, takie jak obliczanie wartości niektórych standardowych funkcji matematycznych ($\sqrt{},\cos (),\sin (),\exp (),\log (),$ itp.), czy nawet funkcji bardziej skomplikowanych. Na przykład, zastosowanie "szkolnych" wzorów na obliczanie pierwiatków równania kwadratowego byłoby niemożliwe, gdyby pierwiastkowanie było niemożliwe.

Uwaga

Należy pamiętać, że w praktyce funkcje standardowe (o ile są dopuszczalne) są obliczane używając czterech podstawowych operacji arytmetycznych. Dokładniej, jednostka arytmetyczna procesora potrafi wykonywać jedynie operacje $+,-,\times ,÷$, przy czym dzielenie zajmuje kilka razy więcej czasu niż pozostałe operacje arytmetyczne. Niektóre procesory (np. Itanium 2) są w stanie wykonywać niejako jednocześnie operację dodawania i mnożenia (tzw. FMADD, fused multiply and add). Praktycznie wszystkie współczesne procesory mają także wbudowany koprocesor matematyczny, realizujący asemblerowe polecenia wyznaczenia wartości standardowych funkcji matematycznych ($\sqrt{},\cos (),\sin (),\exp (),\log (),$ itp.), jednak wykonanie takiej instrukcji wymaga mniej więcej kilkadziesiąt, a czasem nawet kilkaset razy więcej czasu niż wykonanie operacji dodawania.

Niestety, aby nasz model obliczeniowy wiernie odpowiadał rzeczywistości, musimy w nim uwzględnić fakt, że działania matematyczne (ani tym bardziej obliczanie wartości funkcji matematycznych) nie są wykonywane dokładnie. Czasem uwzględnienie tego faktu wiąże się ze znaczącym wzrostem komplikacji analizy algorytmu i dlatego "w pierwszym przybliżeniu" często pomija się to ograniczenie przyjmując model w którym wszystkie (lub prawie wszystkie) działania arytmetyczne wykonują się dokładnie. Wiedza o tym, kiedy i jak zrobić to tak, by wciąż wyciągać prawidłowe wnioski odnośnie faktycznej realizacji algorytmów w obecności błędów zaokrągleń jest częścią sztuki i wymaga intuicji numerycznej, popartej doświadczeniem.

Mamy trzy podstawowe polecenia złożone: warunkowe, powtarzania i kombinowane.

Warunkowe

 
if(<math>\displaystyle \cal W</math>) 
	<math>\displaystyle {\cal A}_1</math>;
else
	<math>\displaystyle {\cal A}_2</math>;

gdzie ${𝒲}$ jest wyrażeniem o wartościach całkowitych (0 odpowiada logicznemu fałszowi, inne wartości --- logicznej prawdzie), a ${{𝒜}}_1$ i ${{𝒜}}_2$ są poleceniami, przy czym dopuszczamy polecenia puste.

Powtarzane

 
while(<math>\displaystyle {\cal W}</math>)
	<math>\displaystyle {\cal A}</math>;

gdzie $W$ jest wyrażeniem o wartościach logicznych, a ${𝒜}$ jest poleceniem.

Kombinowane

 
{
	<math>\displaystyle {\cal A}_1;\displaystyle {\cal A}_2;\displaystyle \ldots\displaystyle {\cal A}_n;</math>
}

gdzie ${{𝒜}}_j$ są poleceniami.

Na podstawie tych trzech poleceń można tworzyć inne, takie jak pętle for(), czy switch(), itd.

Mamy też dwa szczególne polecenia, które odpowiadają za "wejście" i "wyjście".

Wprowadzanie danych

 
  <math>\displaystyle {\cal IN}</math>(x,t);

gdzie $x$ jest zmienną rzeczywistą, a $t$ "adresem" pewnego funkcjonału Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \displaystyle L:F\toR} należącym to pewnego zbioru $T$. W wyniku wykonania tego polecenia w zmiennej $x$ zostaje umieszczona wartość $L_t(f)$.

Polecenie to pozwala zdobyć informację o danej $f$. Jeśli $F=R^n$ to zwykle mamy $T=\{1,2,\dots ,n\}$ i $L_i(f)=f_i$, co w praktyce odpowiada wczytaniu $i$-tej współrzędnej wektora danych. W szczególności, ciąg poleceń ${ℐ}{𝒩}(x[i],i)$, $i=1,2,\dots ,n$, pozwala uzyskać pełną informację o $f$. Jeśli zaś $F$ jest pewną klasą funkcji Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \displaystyle f:[a,b]\toR} , to możemy mieć np. $T=[a,b]$ i $L_t(f)=f(t)$. W tym przypadku, wykonanie polecenia ${ℐ}{𝒩}(x,t)$ odpowiada w praktyce skorzystaniu ze specjalnej procedury (albo urządzenia zewnętrznego) obliczającej (mierzącego) wartość funkcji $f$ w punkcie $t$.

Wyprowadzanie wyników

 
  <math>\displaystyle {\cal OUT}</math>(<math>\displaystyle {\cal W}</math>);

gdzie ${𝒲}$ jest wyrażeniem o wartościach rzeczywistych. Polecenie to pozwala "wskazać" kolejną współrzędną wyniku.

Zakładamy, że na początku procesu obliczeniowego wartości wszystkich zmiennych są nieokreślone, oraz że dla dowolnych danych wykonanie programu wymaga wykonania skończonej liczby operacji elementarnych. Wynikiem obliczeń jest skończony ciąg liczb rzeczywistych (albo $puste$), którego kolejne współrzędne pokazywane są poleceniem ${𝒪}{𝒰}{𝒯}$.

Środowisko obliczeniowe

Ze względu na swój utylitarny charakter, w metodach numerycznych niemal równie ważna jak dobór optymalnego algorytmu jest jego efektywna implementacja na konkretnej architekturze.

W praktyce mamy dwie możliwości:

• wykorzystanie standardowych języków programowania (C, Fortran, być może ze wstawkami w asemblerze) oraz specjalizowanych bibliotek
• użycie gotowego środowiska obliczeń numerycznych będącego wygodnym interfejsem do specjalizowanych bibliotek numerycznych

Języki programowania: C i Fortran

Programy numeryczne (a przynajmniej ich jądra obliczeniowe) są zazwyczaj niezbyt wymagające jeśli chodzi o struktury danych, co więcej, prostota struktur danych szybko rewanżuje się efektywniejszym kodem. Dlatego, trawestując Einsteina, w dobrym programie numerycznym należy używać tak prostych struktur danych, jak to możliwe (ale nie prostszych!...)

Językami programowania opartymi na prostych konstrukcjach są: Fortran i C. Dlatego właśnie są to języki dominujące współcześnie pisane programy numeryczne. O ile w przeszłości hegemonia Fortranu była nie do podważenia, o tyle w chwili obecnej coraz więcej oprogramowania numerycznego powstaje w C.

W naszym wykładzie wybieramy C ze względu na jego uniwersalność, doskonałą przenośność i w chwili obecnej całkiem dojrzałe kompilatory. Dodajmy, że funkcje w C można mieszać z np. z gotowymi bibliotekami napisanymi w Fortranie. Fortran, język o bardzo długiej tradycji, wciąż żywy i mający grono wiernych fanów, jest nadal wybierany przez numeryków na całym świecie między innymi ze względu na jego dopasowanie do zadań obliczeniowych (właśnie w tym celu powstał), a także ze względu na doskonałe kompilatory dostępne na superkomputerach, będące efektem wieloletniej ewolucji i coraz lepszego nie tylko dopasowania kompilatora do spotykanych konstrukcji językowych, ale także na odwrót --- coraz lepszego zrozumienia programistów, jak pisać programy, by wycisnąć jak najwięcej z kompilatora.

Inne popularne języki: Java, Pascal (ten język, zdaje się, jest popularny już tylko w obrębie wydziału MIM UW...), VisualBasic nie są zbyt odpowiednie dla obliczeń numerycznych. Mało tego, np. podstawowy typ numeryczny Pascala: real nie jest zgodny z powszechnym standardem IEEE 754. Jednak, ze względu na coraz większą komplikację kodów numerycznych służących np. do prowadzenia zaawansowanych symulacji metodą elementu skończonego, coraz więcej kodów wykorzystuje możliwości obiektowych języków C++ i Fortran90.

W przykładach będziemy najczęściej odnosić się do architektury x86, tzn. 32-bitowej IA-32 procesorów firmy Intel i AMD, najczęściej spotykanej w obecnie używanych komputerach. Należy jednak pamiętać, że obecnie następuje przejście na architekturę 64-bitową. Ze względu jednak na brak pewności co do ostatecznie przyjętych standardów w tym obszarze, ograniczymy się do procesorów 32-bitowych.

Prosty program numeryczny

Napiszemy teraz program obliczający (w niezbyt wyrafinowany sposób) $N$-tą sumę częściową szeregu harmonicznego

$x=\sum _{i=1}^N\frac{1}{i}.$

Przyjmijmy, że parametr $N$ będzie miał wartość równą 2006. W pierwszym odruchu, prawie każdy początkujący student pisze program w rodzaju:

#include <stdio.h>
#define N 2006

int main(void)
{
float x;
unsigned int i;

x = 0.0;
for(i = 1; i <= N; i++) 
	x = x + 1/i;
printf("Wartość sumy x = 1 + 1/2 + ... + 1/\%d jest równa 
return(0);
}

Sęk w tym, że ten program nie działa! To znaczy, kompiluje się i uruchamia, ale twierdzi uparcie, że nasza suma wynosi... 1.

Winę za to ponosi linijka

 
	x = x + 1/i;

w której wykonujemy dzielenie 1/i. Obie liczby są typu int i dlatego, zgodnie z regułami C, wynik ich dzielenia także jest całkowity: dostajemy część całkowitą z dzielenia tych dwóch liczb, czyli, gdy tylko $i>1$, po prostu zero.

Prawidłowy program musi więc wymusić potraktowanie choć jednej z tych liczb jako liczby zmiennoprzecinkowej, co najprościej uzyskać albo przez

 
	x = x + 1.0/i;

albo bardziej uniwersalnie, rzutując choć jedną z liczb na odpowiedni typ:

 
	x = x + 1/((float) i);

Poprawny kod miałby więc postać

#include <stdio.h>
#define N 2006

int main(void)
{
float x;
unsigned int i;

x = 0.0;
for(i = 1; i <= N; i++) 
	x = x + 1.0/i;
printf("Wartość sumy x = 1 + 1/2 + ... + 1/\%d jest równa 
return(0);
}

W języku C mamy dostępnych sześć typów numerycznych:

• stałoprzecinkowe, dla reprezentacji liczb całkowitych
  • int oraz long int. W realizacji GCC na komputery klasy PC, oba typy: int i long int są identyczne (32-bitowe) i ich zakres wynosi około $-2.1\cdot 10^9\dots +2.1\cdot 10^9$. Typ int i jemu pokrewne odnoszą się do liczb całkowitych ze znakiem (dodatnich lub ujemnych). Ich warianty bez znaku: unsigned int, itp. odnoszą się do liczb bez znaku (nieujemnych), dlatego np. zakresem unsigned int będzie w przybliżeniu

$0\dots +4.2\cdot 10^9$.

• • long long int (64-bitowy) o zakresie w przybliżeniu $-9.2\cdot 10^{18}\dots +9.2\cdot 10^{18}$.
• zmiennopprzecinkowe, dla reprezentacji liczb rzeczywistych
  • float, pojedynczej precyzji, 32-bitowy, gwarantuje precyzję około $10^{-7}$ i zakres liczb reprezentowalnych w przybliżeniu $10^{-38}\dots 10^{38}$;
  • double, podwójnej precyzji, 64-bitowy, ma precyzję na poziomie $10^{-16}$ przy orientacyjnym zakresie $10^{-308}\dots 10^{308}$;
  • long double, rozszerzonej podwójnej precyzji, na pecetach 80-bitowy, ale w pamięci zajmuje on 12 bajtów) o precyzji rzędu $10^{-20}$ i odpowiadający standardowi double extended IEEE 754.

Powyższe typy zmiennoprzecinkowe w realizacji na PC odpowiadają standardowi IEEE 754. Standardowo, operacje arytmetyczne na obu typach float i double są tak samo pracochłonne, gdyż wszystkie obliczenia w C wykonywane są z maksymalną dostępną precyzją (czyli, na procesorach architektury IA-32 Intela i AMD: w precyzji oferowanej przez typ long double), a następnie dopiero wynik zapisywany do zmiennej reprezentowany jest w stosownym typie . Jednakże typ pojedynczej precyzji float oferuje znacznie większe możliwości optymalizacji uzyskanego kodu a ponadto, obiekty typu float zajmują dwukrotnie mniej miejsca w pamięci niż double, dając możliwość lepszego wykorzystania pamięci podręcznej (cache ) i przetwarzania wektorowego.

Stałe matematyczne i podstawowa biblioteka matematyczna

Język C jest językiem małym i, jak wiadomo, nawet proste operacje wejścia-wyjścia są w istocie nie częścią języka, ale funkcjami (makrami?) bibliotecznymi. Z drugiej strony jednak, jak zorientowaliśmy się, nie stwarza to programiście żadnych specjalnych niedogodności. Podobnie rzecz ma się z funkcjami matematycznymi. Podstawowe funkcje matematyczne ($\sin ,\cos ,\exp ,\dots$, itp.) nie są składnikami języka C, lecz w zamian są zaimplementowane w tzw. standardowej bibliotece matematycznej libm.a; prototypy tych funkcji oraz definicje rozmaitych stałych matematycznych: $\pi ,e,\dots$ itp. znajdują się w pliku nagłówkowym math.h. Aby więc skorzystać z tych funkcji w programie, należy

• w nagłówku pliku, w którym korzystamy z funkcji lub stałych matematycznych,

umieścić linię #include <math.h>

• przy linkowaniu dołączyć bibliotekę matematyczną za pomocą opcji -lm

Przykład

Oto przykładowy prosty program numeryczny w C; drukuje on tablicę wartości sinusów losowo wybranych liczb z przedziału $[0,2\pi ]$.

 [{Tablica losowych sinusów}]
#include <math.h>
#include <stdio.h> 
#include <stdlib.h> /* zawiera definicję funkcji rand() i stałej RAND_MAX */
#define N 15 /* ile liczb wydrukować */

int main(void)
{
int i;
double x,y;

	for( i = 0; i < N; i++)
	{
		x = rand()/(double)RAND_MAX;
		x *= 2.0*M_PI;
		/* oczywiście, wystarczyłoby x =(2.0*M_PI*rand())/RAND_MAX; */
		y = sin(x);
		fprintf(stderr, "(\%3d) x = \%10.5e sin(x) = \%10.5e\n", i, x, y);
	}
	return(0);
}

Zwróćmy uwagę na linię x = rand()/(double)RAND_MAX; Funkcja rand() zwraca losową liczbę całkowitą z przedziału [0,RAND_MAX], więc iloraz z predefiniowaną stałą RAND_MAX będzie liczbą z przedziału $[0,1]$. Dla prawidłowego uzyskania losowej liczby z przedziału $[0,1]$ kluczowe jest jednak zrzutowanie jednej z dzielonych liczb na typ double! Gdyby tego nie zrobić, uzyskalibyśmy zawsze $x=0$ lub sporadycznie $x=1$, zgodnie z regułą C typ wyniku jest zgodny z typem argumentów. Rzeczywiście, w naszym wypadku, (błędna) linia x = rand()/RAND_MAX; zostałaby wykonana tak: ponieważ wynikiem funkcji rand() jest int i stała RAND_MAX jest także typu int, to wynik ma być również typu int -- zostanie więc wykonane dzielenie całkowite; ponieważ mamy rand() $\le$ RAND_MAX, to wynikiem będzie albo 0, albo 1, i taki rezultat, po zamianie na typ double, zostanie przypisany zmiennej $x$, co oczywiście nie jest naszym zamiarem. Natomiast, gdy przynajmniej jedna z dzielonych liczb jest typu double, to oczywiście wynikiem też musi być liczba typu double, zostanie więc wykonane zwyczajne dzielenie dwóch liczb zmiennoprzecinkowych (wynik rand() automatycznie zostanie zrzutowany na typ double).

Kompilujemy ten program, zgodnie z uwagami uczynionymi na początku, poleceniem

 
gcc -o sinusy sinusy.c -lm

Środowisko obliczeń numerycznych: MATLAB i jego klony

Typowa sesja MATLABa. Zwróć uwagę na edytor kodu źródłowego na bieżąco interpretujący go i wychwytujący potencjalne błędy

W końcu lat siedemdziesiątych ubiegłego wieku, Cleve Moler wpadł na pomysł stworzenia prostego interfejsu do ówcześnie istniejących bibliotek numerycznych algebry liniowej: pakietów LINPACK i EISPACK

(których był współautorem). Stworzone przez niego: język skryptów i łatwe w użyciu narzędzia manipulacji macierzami, zdobyły wielką popularność w środowisku naukowym. W 1984 roku Moler skomercjalizował swe dzieło, zakładając wspólnie z J. Little'm firmę MathWorks . Jej główny produkt nazwano MATLAB (od: "MATrix LABoratory").

MATLAB wkrótce rozrósł się potężnie, implementując (lub wykorzystując) wiele najlepszych z istniejących algorytmów numerycznych, a także oferując bogate możliwości wizualizacji wyników. Dzięki swemu interfejsowi, składającemu się z prostych, niemal intuicyjnych funkcji, oraz ogromnym możliwościom jest jednym z powszechniej używanych pakietów do prowadzenia symulacji komputerowych w naukach przyrodniczych (i nie tylko).

Interpreter MATLABa jest dosyć wolny, dlatego podstawową regułą pisania efektywnych kodów w MATLABie jest maksymalne wykorzystanie gotowych (prekompilowanych) funkcji MATLABa oraz --- zredukowanie do minimum obecności takich struktur programistycznych jak pętle.

Screenshot Octave. W terminalu otwarta sesja w ascetycznym trybie tekstowym, grafika wyświetlana z wykorzystaniem Gnuplota

Kierując się podobnymi przesłankami co C. Moler, oraz bazując na wielkim sukcesie MATLABa, John W. Eaton z Wydziału Inżynierii Chemicznej Uniwersytetu Wisconsin w USA, zaczął w 1994 roku opracowywać darmowe (udostępniane na tzw. licencji GPL) oprogramowanie o funkcjonalności maksymalnie zbliżonej do MATLABa: Octave. Jak wyjaśnia twórca Octave, , nazwa nie ma nic wspólnego z muzyką, ale z pewną anegdotą. Wersja 1.0 pakietu Octave ukazała się w 1996 roku i jest intensywnie rozwijana do dziś.

Screenshot Scilaba

Drugim udanym klonem MATLABa jest francuski Scilab , opracowany w laboratoriach INRIA i wciąż doskonalony. W subiektywnej ocenie autora niniejszych notatek, nie dorównuje on elegancji i szlachetnej prostocie Octave. Na plus tego pakietu należy zaliczyć m.in. znacznie bardziej rozbudowany zestaw funkcji podstawowych, na minus -- przede wszystkim znacznie mniejszy stopień zgodności z MATLABem, a poza tym: niewygodny system pomocy oraz "toporną" (choć o dużym potencjale) grafikę.

Porównanie MATLABa, Octave i Scilaba

Korzystając z omówionych powyżej pakietów, otrzymujemy:

• możliwość obliczania funkcji matematycznych (nawet dość egzotycznych),

finansowych, analizy sygnałów, itp.;

• bardzo szeroki zakres nowoczesnych narzędzi umożliwiających

wykonywanie podstawowych zadań numerycznych, takich jak: rozwiązywanie równań, obliczanie całek, itd.;

• efektowną wizualizację wyników w postaci wykresów dwu- i trójwymiarowych,

a także funkcje do manipulacji obrazem i dźwiękiem;

• możliwość nieograniczonych rozszerzeń przy użyciu funkcji i skryptów

tworzonych przez osoby trzecie (lub własnych).

Jeśli uwzględnić dodatki wielce ułatwiające życie, w tym rozbudowane funkcje wejścia-wyjścia, to narzędzie robi się rzeczywiście intrygujące, jako środowisko obliczeń numerycznych "dla każdego" (nie tylko dla profesjonalnego numeryka). Stąd wielka popularność MATLABa w środowiskach inżynierkich, gdyż umożliwia szybką implementację rozmaitych --- zwłaszcza testowych --- wersji algorytmu np. przeprowadzania skomplikowanej symulacji. Po zweryfikowaniu wyników, można ewentualnie rozpocząć implementację w docelowym środowisku (np. masywnie równoległym komputerze, w Fortranie 95 z wykorzystaniem komercyjnych bibliotek i specyficznymi optymalizacjami kodu), choć bardzo często wyniki uzyskane w MATLABie w zupełności wystarczą.

Zarówno MATLAB, jak Octave i Scilab, opierają się w dużym zakresie na darmowych --- acz profesjonalnych --- zewnętrznych bibliotekach numerycznych (LAPACK, ATLAS, ARPACK, itp.)

MATLAB jest wciąż niedościgniony jeśli chodzi o narzędzia programistyczne (np. ma wbudowany --- działający na bieżąco w czasie edycji kodu źródłowego! --- debugger) oraz możliwości wizualizacji danych i wyników obliczeń. Jeśli nie dbać o koszta zakupu dodatkowych modułów, ma także najbogatszy zestaw funkcji numerycznych. W chwili pisania niniejszego tekstu, (stosunkowo tanie) wersje studenckie MATLABa, nie są oferowane w Polsce przez producenta, choć z drugiej strony są powszechnie dostępne (nawet w wersji profesjonalnej) w komputerowych laboratoriach akademickich w Polsce.

Wybór pomiędzy darmowymi klonami MATLABa: Octave i Scilabem, jest właściwie kwestią osobistych preferencji. W ninijeszym kursie będziemy posłygiwać się Octave, który jest bardzo wiernym klonem MATLABa, przez co daje zaskakująco wysoki stopień kompatybilności z programami napisanymi w MATLABie. To jest bardzo ważne w praktyce, ponieważ jest dostępnych wiele (nawet darmowych) pakietów numerycznych opracowanych właśnie w języku poleceń MATLABa. Umiejąc posługiwać się Octave, bez trudu można "przesiąść się" na MATLABa i korzystać z jego bogatszych możliwości, droga w przeciwną stronę też jest nietrudna.

Inną zaletą Octave jest swobodniejsza niż w MATLABie składnia: przykładowo, w Octave zarówno 'Newton' jak i "Newton" oznacza łańcuch znakowy (string) --- w MATLABie możliwa jest tylko pierwsza wersja; w Octave możemy pisać pętle for ... endfor jak i for ... end --- w MATLABie możliwa jest tylko druga postać; w Octave możemy pisać zarówno z = z+5 jak i, wzorem języka C, z += 5 --- MATLAB nie dopuszcza drugiego z wariantów; znakiem komentarza w Octave może być # oraz używany w tym celu MATLABie znak definiowanie funkcji w trybie interaktywnym, natomiast funkcje MATLABa muszą być zapisane w zewnętrznym pliku tekstowym.

Zarówno w MATLABie, jak i w Octave, jest możliwe pisanie funkcji, które będą prekompilowane, dając znaczące przyspieszenie działania.

Octave jest dołączany do większości popularnych dystrybucji Linuxa, najczęściej jednak użytkownik musi samodzielnie doinstalować go z płytki instalacyjnej. Ponadto, kody źródłowe najświeższej wersji Octave są dostępne na stronie http://www.octave.org. Dodatkowo, pod http://octave.sourceforge.net znajdziemy pakiet rozszerzeń do Octave, pod nazwą Octave-forge. Wymaga on od użytkowników Linuxa samodzielnej (przebiegającej bezproblemowo) instalacji. Octave można także zainstalować pod Windows, korzystając z programu instalacyjnego dostępnego pod adresem internetowym http://octave.sourceforge.net. W Windowsowej wersji Octave, pakiet Octave-forge jest standardowo dołączony.

Pierwsze kroki w Octave

Octave uruchamiamy poleceniem octave,

 
[przykry@bit MN] octave
GNU Octave, version 2.9.5 (i686-redhat-linux-gnu).
Copyright (C) 2006 John W. Eaton.
This is free software; see the source code for copying conditions.
There is ABSOLUTELY NO WARRANTY; not even for MERCHANTIBILITY or
FITNESS FOR A PARTICULAR PURPOSE.  For details, type `warranty'.

Additional information about Octave is available at http://www.octave.org.

Please contribute if you find this software useful.
For more information, visit http://www.octave.org/help-wanted.html

Report bugs to <bug@octave.org> (but first, please read
http://www.octave.org/bugs.html to learn how to write a helpful report).

octave:1>

Naturalnie, Octave zachowuje wszystkie możliwości kalkulatora naukowego, wykonując podstawowe operacje arytmetyczne i obliczając wartości wielu funkcji matematycznych; oto zapis takiej sesji Octave:

 
octave:1> 1900/2000
ans = 0.95000
octave:2> sin(pi)
ans =  1.2246e-16
octave:3> e^(i*pi)
ans = -1.0000e+00 + 1.2246e-16i

Ostatnie dwa wyniki dają nam namacalny dowód, że obliczenia wykonywane są numerycznie, ze skończoną precyzją: w Octave niektóre tożsamości matematyczne są spełnione jedynie w przybliżeniu, np. $\sin (\pi )\approx 0$ oraz $e^{i\pi }\approx -1$. Przy okazji widzimy, że Octave dysponuje podstawowymi stałymi matematycznymi (Oczywiście, są to także wartości przybliżone!): e $\approx$ 2.71828182845905, pi $\approx$ 3.14159265358979 oraz jednostką urojoną i $=\sqrt{-1}$.

Dokumentacja Octave

Octave ma dobrą dokumentację, dostępną w trakcie sesji Octave; dla każdej funkcji, stałej itp. można uzyskać jej dobry opis przy użyciu polecenia help.

Macierz: postawowy obiekt w Octave

Podstawowym obiektem, z jakim mamy do czynienia w Octave, jest macierz dwuwymiarowa $n\times m$, o elementach rzeczywistych lub zespolonych,

$\left(\begin{array}{ccc}{a}_{11}& \cdots & a_{1m}\\ ⋮& & ⋮\\ a_{n1}& \cdots & a_{nm}\end{array}\right).$

W szczególności, biorąc $m=1$ albo $n=1$, dostajemy wektor (kolumnowy lub wierszowy

$\left(\begin{array}{c}{a}_1\\ ⋮\\ a_n\end{array}\right),\text{lub}\left(\begin{array}{ccc}{a}_1& \cdots & a_m\end{array}\right).$

Z kolei dla $m=n=1$ mamy do czynienia ze "zwykłymi" liczbami rzeczywistymi lub zespolonymi.

Octave umożliwia, za swoim pierwowzorem --- MATLABem, bardzo wygodną pracę z macierzami.

Wprowadzanie macierzy do Octave jest bardzo intuicyjne i możemy zrobić to na kilka sposobów, w zależności od potrzeb. Można robić to interakcyjnie, podając elementy macierzy z linii komend, na przykład:

 
octave:6> A = [2 -1 0; -1 3 -2; 2 2.71828 3.14]  
A =

   2.00000  -1.00000   0.00000
  -1.00000   3.00000  -2.00000
   2.00000   2.71828   3.14000

definiuje macierz kwadratową $3\times 3$ o elementach równych

$\left(\begin{array}{ccc}2& -1& 0\\ -1& 3& -2\\ 2& 2.71828& 3.14\end{array}\right).$

Tak więc, macierz w Octave definiujemy przez podanie jej elementów. Elementy wprowadzamy kolejno wierszami, elementy w wierszu oddzielamy spacjami (lub ewentualnie przecinkami), natomiast kolejne wiersze oddzielamy średnikami.

Macierz możemy także tworzyć stopniowo, wprowadzając jeden za drugim jej kolejne elementy:

 
octave:3> B(1,1) = 4
B = 4

octave:4> B(2,1) = 3 + B(1,1)
B =

  4
  7

octave:5> B(3,2) = 28
B =

   4   0
   7   0
   0  28

Pierwsza komenda określa B jako macierz $1\times 1$, czyli zwykłą liczbę, druga -- rozszerza nasz obiekt do wektora, natomiast trzecia -- powoduje, że B zostaje rozszerzony do macierzy $3\times 2$. Zauważmy, że elementom B o nieokreślonych przez nas wartościach zostaje przypisana domyślna wartość zero.

W przypadku, gdy wypełniana macierz jest duża, znacznie korzystniej jest prealokować macierz, ustawiając na samym początku wszystkie elementy macierzy docelowego wymiaru (u nas to była macierz $3\times 2$) na zero; potem możemy już działać jak poprzednio:

 
B = zeros(3,2)
B(1,1) = 3 + B(1,1)
B(2,1) = pi
B(3,2) = 28

Jeśli nie chcemy za każdym razem obserwować echa prowadzonych przez nas działań na macierzy B, komendy powinniśmy kończyć średnikiem: średnik po komendzie Octave blokuje echo wyniku jej wykonania. Tak więc, znacznie wygodniej będzie nam napisać

 
B = zeros(3,2);
B(1,1) = 3 + B(1,1);
B(2,1) = pi;
B(3,2) = 28;
B

Ostatnie polecenie -- B bez średnika na końcu -- spowoduje wypisanie zawartości B na ekran.

Tak utworzoną macierz możemy zapisać do pliku, wykorzystując jeden z kilku formatów:

tekstowy

Zapis do formatu ASCII. Ze względu na późniejszą konwersję z zapisu w

postaci dzisiętnej do dwójkowej, zob. \link{xxx}, wiąże się z drobnymi niedokładnościami w odczycie.

 
save macierzb.dat B

W wyniku dostajemy

# Created by Octave 2.9.5, Mon Jul 31 17:15:25 2006 CEST <przykry@bit>
# name: B
# type: matrix
# rows: 3
# columns: 2
 3 0
 3.14159265358979 0
 0 28

binarny

Po prostu zrzut stanu pamięci

 
save -binary  macierzb.dat B

HDF5

Format gwarantujący pełną przenośność danych pomiędzy różnymi

architekturami, akceptowany przez wiele zewnętrznych aplikacji, m.in. OpenDX.

 
save -hdf5  macierzb.dat B

Dane z pliku wczytujemy z powrotem poleceniem

 
load macierzb.dat

Sesję Octave kończymy poleceniem quit.

Wektoryzacja

Ponieważ podstawowym obiektem w Octave jest macierz, predefiniowane operacje matematyczne wykonują się od razu na całej macierzy. Bez żadnej przesady możena stwierdzić, że umiejętność wektoryzacji i blokowania algorytmów jest podstawą pisania efektywnych implementacji algorytmów w Octave.

Zobaczmy kilka prostych przykładów zawartych w tabeli poniżej. W pierwszej kolumnie tabeli, dane zadanie zaimplementujemy w Octave przy użyciu pętli (zwróćmy przy okazji uwagę na to, jak programuje się pętle w Octave, one nie są zupełnie bezużyteczne...). W drugiej kolumnie zobaczymy, jak to samo zadanie można wykonać korzystając z operatorów lub funkcji macierzowych.

	Tradycyjny kod w Octave, używający pętli	Efektywny kod wektorowy (macierzowy) w Octave
1 || s = 0; for i = 1:size(x,1) s = s + abs(x(i)); end	s = sum(abs(x));
2 || N = 500; h = (b-a)/(N-1); for i = 1:N x(i) = a + (i-1)*h; y(i) = sin(x(i)); end plot(x,y);	N = 500; x = linspace(a,b,N); y = sin(x); plot(x,y);
3a || for i = 1:size(C,1), for j = 1:size(C,2), for k = 1:size(A,2), C(i,j) = C(i,j) + A(i,k)*B(k,j); end end end	C = C + A*B;
3b || for i = 1:size(C,1), for j = 1:size(C,2), C(i,j) = C(i,j) + A(i,:)*B(:,j); end end	C = C + A*B;

Zwróćmy uwagę na to, że kod wektorowy Octave oraz, w jeszcze większym stopniu, kod macierzowy jest znacznie bardziej elegancki i czytelniejszy od tradycyjnego. Widać to dobrze na wszystkich powyższych przykładach.

Pierwszy przykład pokazuje też, że kod macierzowy jest elastyczniejszy: albowiem obliczy sumę modułów wszystkich elementów x nie tylko gdy x jest wektorem (co musieliśmy milcząco założyć w kodzie w lewej kolumnie), ale tak samo dobrze zadziała, gdy x jest macierzą!

Szczególnie wiele na czytelności i wygodzie zyskujemy w drugim przykładzie, gdzie najpierw funkcja linspace pozwala nam uniknąć żmudnego wyznaczania $N$ równoodległych węzłów $x_i$ w odcinku $[a,b]$, a następnie funkcja sin zaaplikowana do całego wektora $x$ daje wartość sinusa w zadanych przez nas węzłach.

Kod wektorowy lub (zwłaszcza) macierzowy jest też znacznie szybszy. Spójrzmy teraz na przykłady: trzeci i czwarty, które pokażą nam prawdziwą moc funkcji macierzowych, unikających wielokrotnie zagnieżdżonych pętli. Oba dotyczą operacji mnożenia dwóch macierzy. Przykład trzeci w wersji z potrójną pętlą for naśladuje sposób programowania znany nam z C lub Pascala, natomiast przykład czwarty zdaje się być napisany odrobinę w duchu wektorowym (brak trzeciej, wewnętrznej pętli, zastąpionej operacją wektorową: iloczynem skalarnym $i$-tego wiersza $A$ i $j$-tej kolumny $B$). Poniżej porównanie czasów działania tych trzech implementacji w przypadku macierzy $64\times 64$ (czasy dla PC z procesorem Celeron 1GHz):

• Dla pętli postaci C(i,j)=C(i,j)+A(i,k)*B(k,j) uzyskano czas

21.6s,

• Dla pętli postaci C(i,j)=C(i,j)+A(i,:)*B(:,j) --- 0.371s
• Dla pętli postaci C=C+A*B kod działał jedynie 0.00288s!

Widzimy, jak beznadziejnie wolny jest kod oparty na trzech zagnieżdżonych pętlach: jest on kilka tysięcy razy wolniejszy od implementacji macierzowej C = C + A*B. Po wektoryzacji wewnętrznej pętli, program doznaje kilkudziesięciokrotnego przyspieszenia, lecz nadal jest ponadstukrotnie wolniejszy od kodu macierzowego!