MNwyklad01: Różnice pomiędzy wersjami

Z Studia Informatyczne
Przejdź do nawigacjiPrzejdź do wyszukiwania
m
m
 
Linia 1: Linia 1:
  
=Wprowadzenie do metod numerycznych=
 
 
<strong>Metody numeryczne</strong> 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 <strong>matematyki ciągłej</strong> (w odróżnieniu od [[Dodaj WIKIlink|matematyki
 
dyskretnej]]).
 
 
==Zadania metod numerycznych==
 
 
Aby w ogóle mówić w problemie obliczeniowym, musimy najpierw
 
* określić <strong>dane problemu</strong> i <strong>cel obliczeń</strong>, czyli dokładnie
 
sformułować zadanie w języku matematyki,
 
* określić <strong>środki obliczeniowe</strong> dzięki którym chcemy  osiągnąć cel,
 
* dla analizy zadania i sposobów jego rozwiązania wygodnie jest zdefiniować
 
<strong>klasę rozpatrywanych danych</strong> oraz <strong>model obliczeniowy</strong> w obrębie
 
którego będą działać nasze algorytmy.
 
 
 
Wbrew dość powszechnej opinii <strong>nie jest prawdą</strong>, że głównym przedmiotem 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;
 
 
; 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 <strong>esencji</strong> algorytmu, bez detali implementacyjnych ---
 
zostawiając je na inną okazję (dobra implementacja konkretnego algorytmu może być
 
sama w sobie interesującym wyzwaniem programistycznym; często bywa, że dobre
 
implementacje, nawet prostych algorytmów numerycznych, są mało czytelne).
 
 
Aby zdefiniować nasz model obliczeniowy, posłużymy się
 
pojęciem <strong>programu</strong>. Zastosujemy przy tym notację
 
podobną do tej z języka programowania <strong>C</strong>.
 
 
Program składa się z <strong>deklaracji</strong>, czyli opisu obiektów,
 
których będziemy używać, oraz z <strong>poleceń</strong> (<strong>instrukcji</strong>),
 
czyli opisu akcji, które będziemy wykonywać.
 
 
Dostępnymi obiektami są <strong>stałe</strong> i <strong>zmienne</strong> typu
 
całkowitego (<code>int</code>),
 
rzeczywistego (<code>float</code> i <code>double</code>). 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===
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
  z <nowiki>=</nowiki>  <math>\displaystyle {\cal W}</math>;
 
</pre></div>
 
 
gdzie <math>\displaystyle z</math> jest zmienną, a <math>\displaystyle {\cal W}</math> jest <strong>wyrażeniem</strong>
 
o wartościach tego samego typu co <math>\displaystyle z</math>. Jest to polecenie proste.
 
 
Wyrażeniem jest pojedyncza stała lub zmienna, albo złożenie
 
skończonej liczby <strong>operacji elementarnych</strong> na wyrażeniach. 
 
Operacje elementarne to:
 
 
; arytmetyczno--arytmetyczne:
 
:  <math>\displaystyle x\mapsto -x</math>, <math>\displaystyle (x,y)\mapsto x+y</math>,
 
<math>\displaystyle (x,y)\mapsto x-y</math>, <math>\displaystyle (x,y)\mapsto x*y</math>,
 
<math>\displaystyle (x,y)\mapsto x/y, y\ne 0</math>, gdzie <math>\displaystyle x,y</math> są stałymi lub
 
zmiennymi liczbowymi,
 
 
; arytmetyczno--logiczne:
 
:  <math>\displaystyle (x,y)\mapsto x<y</math>,
 
<math>\displaystyle (x,y)\mapsto x\le y</math>, <math>\displaystyle (x,y)\mapsto x=y</math>, <math>\displaystyle (x,y)\mapsto x\ne y</math>,
 
gdzie <math>\displaystyle x,y</math> są stałymi lub zmiennymi liczbowymi,
 
 
; logiczno--logiczne:
 
:  <math>\displaystyle p\mapsto\,\sim \,p</math>,
 
<math>\displaystyle (p,q)\mapsto p\, \wedge \,q</math>, <math>\displaystyle (p,q)\mapsto p\, \vee \,q</math>,
 
gdzie <math>\displaystyle p,q</math> są stałymi lub zmiennymi logicznymi.
 
 
Dla niektórych zadań wygodnie jest (a czasem konieczne)
 
uzupełnić ten zbiór o dodatkowe operacje, takie jak
 
obliczanie wartości niektórych standardowych funkcji matematycznych
 
(<math>\displaystyle \sqrt{\;}, \cos(), \sin(), \exp(), \log(),</math> 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 <math>\displaystyle +,\,-,\,\times,\,\div</math>, 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 (<math>\displaystyle \sqrt{\;}, \cos(), \sin(),
 
\exp(), \log(),</math> itp.), jednak wykonanie takiej instrukcji wymaga mniej więcej
 
kilkadziesiąt, a czasem nawet kilkaset  razy więcej czasu
 
niż wykonanie operacji dodawania.
 
}}
 
 
{{uwaga|||
 
Niestety, aby nasz  model obliczeniowy wiernie odpowiadał rzeczywistości, musimy
 
w nim uwzględnić fakt, że działania matematyczne (i, tym bardziej,
 
obliczanie wartości funkcji matematycznych) <strong>nie są</strong> 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, <strong>kiedy</strong> i
 
<strong>jak</strong> 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===
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
if(<math>\displaystyle \cal W</math>)
 
<math>\displaystyle {\cal A}_1</math>;
 
else
 
<math>\displaystyle {\cal A}_2</math>;
 
</pre></div>
 
 
gdzie <math>\displaystyle {\cal W}</math> jest wyrażeniem o wartościach całkowitych (0 odpowiada
 
logicznemu fałszowi, inne wartości --- logicznej prawdzie), a <math>\displaystyle {\cal A}_1</math>
 
i <math>\displaystyle {\cal A}_2</math> są poleceniami, przy czym dopuszczamy polecenia puste.
 
 
===Powtarzane===
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
while(<math>\displaystyle {\cal W}</math>)
 
<math>\displaystyle {\cal A}</math>;
 
</pre></div>
 
 
gdzie <math>\displaystyle W</math> jest wyrażeniem o wartościach logicznych, a <math>\displaystyle \cal A</math>
 
jest poleceniem.
 
 
===Kombinowane===
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
{
 
<math>\displaystyle {\cal A}_1;\displaystyle {\cal A}_2;\displaystyle \ldots\displaystyle {\cal A}_n;</math>
 
}
 
</pre></div>
 
 
gdzie <math>\displaystyle {\cal A}_j</math> są poleceniami.
 
 
Na podstawie tych trzech poleceń można tworzyć inne, takie
 
jak pętle <code>for()</code>, czy <code>switch()</code>, itd.
 
 
Mamy też dwa szczególne polecenia, które odpowiadają
 
za "wejście" i "wyjście".
 
 
===Wprowadzanie danych===
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
  <math>\displaystyle {\cal IN}</math>(x,t);
 
</pre></div>
 
 
gdzie <math>\displaystyle x</math> jest zmienną rzeczywistą, a <math>\displaystyle t</math> "adresem" pewnego
 
funkcjonału <math>\displaystyle L:F\toR</math> należącym to pewnego zbioru <math>\displaystyle T</math>.
 
W wyniku wykonania tego polecenia w zmiennej <math>\displaystyle x</math> zostaje
 
umieszczona wartość <math>\displaystyle L_t(f)</math>.
 
 
Polecenie to pozwala zdobyć <strong>informację</strong> o danej <math>\displaystyle f</math>.
 
Jeśli <math>\displaystyle F=R^n</math> to zwykle mamy <math>\displaystyle T=\{1,2,\ldots,n\}</math> i
 
<math>\displaystyle L_i(f)=f_i</math>, co w praktyce odpowiada wczytaniu <math>\displaystyle i</math>-tej
 
współrzędnej wektora danych. W szczególności, ciąg poleceń
 
<math>\displaystyle {\cal IN}(x[i],i)</math>, <math>\displaystyle i=1,2,\ldots,n</math>, pozwala uzyskać pełną
 
informację o <math>\displaystyle f</math>. Jeśli zaś <math>\displaystyle F</math> jest pewną klasą
 
funkcji <math>\displaystyle f:[a,b]\toR</math>, to możemy mieć np. <math>\displaystyle T=[a,b]</math> i <math>\displaystyle L_t(f)=f(t)</math>.
 
W tym przypadku, wykonanie polecenia <math>\displaystyle {\cal IN}(x,t)</math> odpowiada
 
w praktyce skorzystaniu ze specjalnej procedury (albo urządzenia
 
zewnętrznego) obliczającej (mierzącego) wartość funkcji <math>\displaystyle f</math>
 
w punkcie <math>\displaystyle t</math>.
 
 
===Wyprowadzanie wyników===
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
  <math>\displaystyle {\cal OUT}</math>(<math>\displaystyle {\cal W}</math>);
 
</pre></div>
 
 
gdzie <math>\displaystyle {\cal W}</math> 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 <math>\displaystyle puste</math>), którego kolejne współrzędne
 
pokazywane są poleceniem <math>\displaystyle {\cal OUT}</math>.
 
 
==Ś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 <strong>implementacja</strong> 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 wyspecjalizowanych bibliotek
 
* użycie gotowego środowiska obliczeń numerycznych będącego wygodnym interfejsem do  specjalizowanych bibliotek numerycznych
 
 
Zaleta pierwszego podejścia to (zazwyczaj) szybko działający kod wynikowy, ale kosztem długotrwałego i żmudnego programowania. W drugim przypadku jest na odwrót: developerka i testowanie --- wyjątkowo ważne w przypadku programu numerycznego --- postępują bardzo szybko, ale czasem kosztem ogólnej efektywności uzyskanego produktu.
 
 
====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
 
 
<blockquote  style="background-color:#fefeee"> 
 
używać tak prostych struktur danych, jak to
 
możliwe (ale nie prostszych!...)
 
</blockquote>
 
 
Językami opartymi na prostych konstrukcjach programistycznych 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: [[Dodaj WIKIlink|Java]], Pascal (ten język, zdaje się, jest popularny już
 
tylko w obrębie [http://www.mimuw.edu.pl  wydziału MIM UW]...), VisualBasic i inne,  nie są zbyt odpowiednie dla
 
obliczeń numerycznych. Mało tego, np. podstawowy typ numeryczny Pascala:
 
<code>real</code> nie jest zgodny z powszechnym standardem [[Dodaj WIKIlink|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 [http://developer.intel.com  Intel] i [http://developer.amd.com  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. W przykładach będziemy korzystać z kompilatora [http://gcc.gnu.org  GCC], który jest omówiony w wykładzie [[Dodaj WIKIlink|Środowisko programistyczne]].
 
 
====Prosty program numeryczny w C====
 
 
Napiszemy teraz program obliczający (w niezbyt wyrafinowany sposób) <math>\displaystyle N</math>-tą sumę
 
częściową szeregu harmonicznego
 
<center><math>\displaystyle
 
x = \sum_{i=1}^N \frac{1}{i}.
 
</math></center>
 
 
Przyjmijmy, że parametr <math>\displaystyle N</math> będzie miał wartość równą 2006. W pierwszym odruchu,
 
prawie każdy początkujący student pisze program w rodzaju:
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
#include <stdio.h>
 
#define N 2006
 
 
int main(void)
 
{
 
float x;
 
unsigned int i;
 
 
x <nowiki>=</nowiki>  0.0;
 
for(i <nowiki>=</nowiki>  1; i <<nowiki>=</nowiki>  N; i++)
 
x <nowiki>=</nowiki>  x + 1/i;
 
printf("Wartość sumy x <nowiki>=</nowiki>  1 + 1/2 + ... + 1/&#37;d jest równa &#37;g\n",  N, x);
 
 
return(0);
 
}
 
</pre></div>
 
 
Sęk w tym, że ten program <strong>nie działa!</strong> To znaczy: owszem, kompiluje się i
 
uruchamia, ale twierdzi uparcie, że nasza suma wynosi... 1.
 
 
Winę za to ponosi linijka
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
x <nowiki>=</nowiki>  x + 1/i;
 
</pre></div>
 
 
w której wykonujemy dzielenie <code>1/i</code>. Obie liczby są typu <code>int</code> 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 <math>\displaystyle i >
 
1</math>, 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
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
x <nowiki>=</nowiki>  x + 1.0/i;
 
</pre></div>
 
 
albo bardziej uniwersalnie, rzutując choć jedną z liczb na odpowiedni typ:
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
x <nowiki>=</nowiki>  x + 1/((float) i);
 
</pre></div>
 
 
Poprawny kod miałby więc postać
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
#include <stdio.h>
 
#define N 2006
 
 
int main(void)
 
{
 
float x;
 
unsigned int i;
 
 
x <nowiki>=</nowiki>  0.0;
 
for(i <nowiki>=</nowiki>  1; i <<nowiki>=</nowiki>  N; i++)
 
x <nowiki>=</nowiki>  x + 1.0/i;
 
printf("Wartość sumy x <nowiki>=</nowiki>  1 + 1/2 + ... + 1/&#37;d jest równa
 
return(0);
 
}
 
</pre></div>
 
 
====Typy numeryczne w C====
 
 
W języku C mamy dostępnych sześć typów numerycznych:
 
* stałoprzecinkowe, dla reprezentacji liczb całkowitych
 
** <code>int</code> oraz <code>long int</code>.  W realizacji [http://gcc.gnu.org  GCC] na komputery klasy PC, oba typy: <code>int</code> i <code>long int</code> są identyczne (32-bitowe) i ich zakres wynosi około <math>\displaystyle -2.1\cdot 10^9\ldots +2.1\cdot 10^9</math>. Typ <code>int</code> i jemu pokrewne  odnoszą się do liczb całkowitych ze znakiem (dodatnich lub ujemnych). Ich warianty bez znaku: <code>unsigned int</code>, itp. odnoszą się do liczb bez znaku (nieujemnych), dlatego np. zakresem <code>unsigned int</code> będzie w przybliżeniu
 
<math>\displaystyle 0\ldots +4.2\cdot 10^9</math>.
 
** <code>long long int</code> (64-bitowy) o zakresie w przybliżeniu <math>\displaystyle -9.2\cdot 10^{18}\ldots +9.2\cdot 10^{18}</math>.
 
* zmiennopprzecinkowe, dla reprezentacji liczb rzeczywistych
 
** <code>float</code>, pojedynczej precyzji, 32-bitowy, gwarantuje precyzję około <math>\displaystyle 10^{-7}</math> i zakres liczb reprezentowalnych w przybliżeniu <math>\displaystyle 10^{-38}\ldots 10^{38}</math>;
 
** <code>double</code>, podwójnej precyzji, 64-bitowy, ma precyzję na poziomie <math>\displaystyle 10^{-16}</math> przy orientacyjnym zakresie <math>\displaystyle 10^{-308}\ldots 10^{308}</math>;
 
** <code>long double</code>, rozszerzonej podwójnej precyzji, na pecetach 80-bitowy, ale w pamięci zajmuje on 12 bajtów) o precyzji rzędu <math>\displaystyle 10^{-20}</math> i odpowiadający standardowi double extended IEEE 754.
 
 
Powyższe typy zmiennoprzecinkowe w realizacji na PC odpowiadają
 
[http://www.cs.berkeley.edu/&nbsp;wkahan/ieee754status/754story.html  standardowi IEEE 754].
 
Standardowo, operacje arytmetyczne
 
na obu typach <code>float</code> i <code>double</code> 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
 
<code>long double</code>), a następnie dopiero wynik zapisywany do zmiennej
 
reprezentowany jest w stosownym typie . Jednakże typ
 
pojedynczej precyzji <code>float</code> oferuje znacznie większe możliwości
 
optymalizacji uzyskanego kodu a ponadto, obiekty typu <code>float</code> zajmują
 
dwukrotnie mniej miejsca w pamięci niż <code>double</code>, 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 (<math>\displaystyle \sin, \cos, \exp,
 
\ldots</math>, itp.) nie są składnikami języka C, lecz w zamian są zaimplementowane w
 
tzw. standardowej bibliotece matematycznej <code>libm.a</code>; prototypy tych
 
funkcji oraz definicje rozmaitych stałych matematycznych: <math>\displaystyle \pi, e, \ldots</math> itp.
 
znajdują się w pliku nagłówkowym <code>math.h</code>. 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ę <code>#include <math.h></code>
 
* przy linkowaniu dołączyć bibliotekę matematyczną za pomocą opcji <code>-lm</code>
 
 
<div style="margin-top:1em; padding-top,padding-bottom:1em;">
 
<span  style="font-variant:small-caps;">Przykład</span>
 
<div class="solution">
 
 
Oto przykładowy prosty program numeryczny w C; drukuje on tablicę wartości
 
sinusów losowo wybranych liczb z przedziału <math>\displaystyle [0,2\pi]</math>.
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
[{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 <nowiki>=</nowiki>  0; i < N; i++)
 
{
 
x <nowiki>=</nowiki>  rand()/(double)RAND_MAX;
 
x *<nowiki>=</nowiki>  2.0*M_PI;
 
/* oczywiście, wystarczyłoby x <nowiki>=</nowiki> (2.0*M_PI*rand())/RAND_MAX; */
 
y <nowiki>=</nowiki>  sin(x);
 
fprintf(stderr, "(&#37;3d) x <nowiki>=</nowiki>  &#37;10.5e sin(x) <nowiki>=</nowiki>  &#37;10.5e\n", i, x, y);
 
}
 
return(0);
 
}
 
</pre></div>
 
 
Zwróćmy uwagę na linię <code>x <nowiki>=</nowiki>  rand()/(double)RAND_MAX;</code> Funkcja <code>rand()</code>
 
zwraca losową liczbę całkowitą z przedziału [0,<code>RAND_MAX</code>], więc iloraz z
 
predefiniowaną stałą <code>RAND_MAX</code> będzie liczbą z przedziału <math>\displaystyle [0,1]</math>. Dla
 
prawidłowego uzyskania losowej liczby z przedziału <math>\displaystyle [0,1]</math> kluczowe jest jednak
 
zrzutowanie jednej z dzielonych liczb na typ <code>double</code>! Gdyby tego nie zrobić,
 
uzyskalibyśmy zawsze <math>\displaystyle x=0</math> lub sporadycznie <math>\displaystyle x=1</math>, zgodnie z regułą C <strong>typ
 
wyniku jest zgodny z typem argumentów</strong>. Rzeczywiście, w naszym wypadku,
 
(błędna) linia <code>x <nowiki>=</nowiki>  rand()/RAND_MAX;</code> zostałaby wykonana tak: ponieważ
 
wynikiem funkcji <code>rand()</code> jest <code>int</code> i stała <code>RAND_MAX</code> jest także typu
 
<code>int</code>, to wynik ma być również typu <code>int</code> -- zostanie więc wykonane dzielenie
 
całkowite; ponieważ mamy <code>rand() <math>\displaystyle \leq</math> RAND_MAX</code>, to wynikiem będzie
 
albo 0, albo 1, i taki rezultat, po zamianie na typ <code>double</code>, zostanie przypisany
 
zmiennej <math>\displaystyle x</math>, co oczywiście nie jest naszym zamiarem. Natomiast, gdy
 
przynajmniej jedna z dzielonych liczb jest typu <code>double</code>, to oczywiście wynikiem
 
też musi być liczba typu <code>double</code>, zostanie więc wykonane zwyczajne dzielenie
 
dwóch liczb zmiennoprzecinkowych (wynik <code>rand()</code> automatycznie zostanie
 
zrzutowany na typ <code>double</code>).
 
 
Kompilujemy ten program, zgodnie z uwagami uczynionymi na początku, poleceniem
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
gcc -o sinusy sinusy.c -lm
 
</pre></div>
 
 
</div></div>
 
 
====Środowisko obliczeń numerycznych: MATLAB i jego klony====
 
 
[[Image:MNmatlab-screenshot.png|frame|400px|center|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, pod nazwą [http://www.mathworks.com  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.
 
 
[[Image:MNoctave-screenshot.png|frame|400px|center|Screenshot Octave. W terminalu otwarta sesja w
 
ascetycznym
 
trybie tekstowym, grafika wyświetlana z wykorzystaniem [http://www.gnuplot.info  Gnuplota]]]
 
 
Kierując się podobnymi przesłankami co C.&nbsp;Moler, oraz bazując na wielkim
 
sukcesie MATLABa, John W.&nbsp;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: [http://www.octave.org  Octave]. Wersja 1.0 pakietu Octave ukazała się w 1996 roku i jest
 
intensywnie rozwijana do dziś.
 
 
[[Image:MNscilab-screenshot.png|frame|400px|center|Screenshot Scilaba.]]
 
 
Drugim udanym klonem MATLABa jest francuski [http://www.scilab.org  Scilab] ,
 
opracowany w laboratoriach INRIA i wciąż doskonalony. W
 
subiektywnej ocenie autorów 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ę.
 
 
Korzystając z dowolnego 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żynierskich, 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 <strong>na bieżąco</strong> 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.
 
 
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.
 
 
Wybór pomiędzy darmowymi klonami MATLABa: Octave i Scilabem, jest właściwie kwestią osobistych preferencji.
 
W ninijeszym kursie będziemy posługiwać 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.
 
 
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 <code>octave</code>,
 
 
<div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
 
 
[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>
 
</pre></div>
 
 
 
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:
 
 
<div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
 
 
octave:1> 1900/2000
 
ans <nowiki>=</nowiki>  0.95000
 
octave:2> sin(pi)
 
ans <nowiki>=</nowiki>  1.2246e-16
 
octave:3> e^(i*pi)
 
ans <nowiki>=</nowiki>  -1.0000e+00 + 1.2246e-16i
 
</pre></div>
 
 
Ostatnie dwa wyniki dają nam namacalny dowód, że obliczenia wykonywane są
 
<strong>numerycznie</strong>, ze skończoną precyzją: w Octave niektóre tożsamości matematyczne
 
są spełnione <strong>jedynie w przybliżeniu</strong>, np. <math>\displaystyle \sin(\pi) \approx 0</math> oraz
 
<math>\displaystyle e^{i\pi} \approx -1</math>. Przy okazji widzimy, że Octave dysponuje podstawowymi
 
stałymi matematycznymi (oczywiście, są to także wartości przybliżone!):
 
<code>e <math>\displaystyle \approx</math> 2.71828182845905</code>,  <code>pi <math>\displaystyle \approx</math> 3.14159265358979</code>
 
oraz jednostką urojoną <code>i</code> <math>\displaystyle = \sqrt{-1}</math>.
 
 
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 <code>help</code>.
 
 
Podstawowym obiektem, z jakim mamy do czynienia w Octave, jest macierz
 
dwuwymiarowa <math>\displaystyle n\times m</math>, o elementach rzeczywistych lub zespolonych,
 
<center><math>\displaystyle
 
\begin{pmatrix}
 
a_{11} & \cdots & a_{1m}\\
 
\vdots &        & \vdots\\
 
a_{n1} & \cdots & a_{nm}
 
\end{pmatrix} .
 
</math></center>
 
 
W szczególności, biorąc <math>\displaystyle m=1</math> albo <math>\displaystyle n=1</math>, dostajemy wektor (kolumnowy lub
 
wierszowy
 
<center><math>\displaystyle
 
\begin{pmatrix}
 
a_{1} \\
 
\vdots\\
 
a_{n}
 
\end{pmatrix} , \qquad\text{lub}\qquad
 
\begin{pmatrix}
 
a_{1} & \cdots & a_{m}
 
\end{pmatrix} .
 
</math></center>
 
 
Z kolei dla <math>\displaystyle m=n=1</math> 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:
 
 
<div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
 
 
octave:6> A <nowiki>=</nowiki>  [2 -1 0; -1 3 -2; 2 2.71828 3.14] 
 
A <nowiki>=</nowiki>
 
 
  2.00000  -1.00000  0.00000
 
  -1.00000  3.00000  -2.00000
 
  2.00000  2.71828  3.14000
 
 
</pre></div>
 
 
definiuje macierz kwadratową <math>\displaystyle 3\times 3</math> o elementach równych
 
<center><math>\displaystyle
 
\begin{pmatrix}
 
2 & -1 & 0\\ -1 & 3 & -2\\ 2 & 2.71828 & 3.14
 
\end{pmatrix} .
 
</math></center>
 
 
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:
 
 
<div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
 
 
octave:3> B(1,1) <nowiki>=</nowiki>  4
 
B <nowiki>=</nowiki>  4
 
 
octave:4> B(2,1) <nowiki>=</nowiki>  3 + B(1,1)
 
B <nowiki>=</nowiki>
 
 
  4
 
  7
 
 
octave:5> B(3,2) <nowiki>=</nowiki>  28
 
B <nowiki>=</nowiki>
 
 
  4  0
 
  7  0
 
  0  28
 
</pre></div>
 
 
Pierwsza komenda określa <code>B</code> jako macierz <math>\displaystyle 1\times 1</math>, czyli zwykłą
 
liczbę, druga -- rozszerza nasz obiekt do wektora, natomiast trzecia --
 
powoduje, że <code>B</code> zostaje rozszerzony do macierzy <math>\displaystyle 3\times 2</math>. Zauważmy, że
 
elementom <code>B</code> 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 <strong>docelowego wymiaru</strong> (u nas to była
 
macierz <math>\displaystyle 3\times 2</math>) na zero; potem możemy już działać jak poprzednio:
 
-->
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
B <nowiki>=</nowiki>  zeros(3,2)
 
B(1,1) <nowiki>=</nowiki>  3 + B(1,1)
 
B(2,1) <nowiki>=</nowiki>  pi
 
B(3,2) <nowiki>=</nowiki>  28
 
</pre></div>
 
 
Jeśli nie chcemy za każdym razem obserwować echa prowadzonych przez nas działań
 
na macierzy <code>B</code>, komendy powinniśmy kończyć średnikiem: średnik po
 
komendzie Octave blokuje echo wyniku jej wykonania. Tak więc, znacznie
 
wygodniej będzie nam napisać
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
B <nowiki>=</nowiki>  zeros(3,2);
 
B(1,1) <nowiki>=</nowiki>  3 + B(1,1);
 
B(2,1) <nowiki>=</nowiki>  pi;
 
B(3,2) <nowiki>=</nowiki>  28;
 
B
 
</pre></div>
 
 
Ostatnie polecenie -- <code>B</code> <strong>bez</strong> średnika na końcu -- spowoduje
 
wypisanie zawartości <code>B</code> na ekran.
 
 
Tak utworzoną macierz możemy zapisać do pliku, wykorzystując jeden z kilku
 
formatów:
 
* tekstowy
 
* binarny
 
* HDF5
 
 
Zapis w formacie tekstowym to po prostu zapis do formatu ASCII. Ze względu na późniejszą [[Dodaj WIKIlink|konwersję z zapisu w
 
postaci dziesiętnej do dwójkowej]], wiąże się z drobnymi
 
niedokładnościami w odczycie.
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
save macierzb.dat B
 
</pre></div>
 
 
W wyniku dostajemy
 
<div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
 
# 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
 
</pre></div>
 
 
Zapis binarny to prostu zrzut stanu pamięci
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
save -binary  macierzb.dat B
 
</pre></div>
 
 
Format [http://hdf.ncsa.uiuc.edu/HDF5  HDF5] to format gwarantujący pełną przenośność danych pomiędzy różnymi
 
architekturami, akceptowany przez wiele zewnętrznych aplikacji, m.in.
 
[http://www.opendx.org  OpenDX].
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
save -hdf5  macierzb.dat B
 
</pre></div>
 
 
Dane z pliku wczytujemy z powrotem poleceniem
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
load macierzb.dat
 
</pre></div>
 
 
Sesję Octave kończymy poleceniem <code>quit</code>.
 
 
====Notacja dwukropkowa Molera====
 
 
Do istniejącej macierzy możemy odwoływać się tradycyjnie, czyli do pojedynczych
 
elementów, np. <code>alfa <nowiki>=</nowiki>  D(2,2)</code>, lub do jej większych bloków, używając popularnej tzw. <strong>notacji dwukropkowej</strong> Molera.
 
 
[[Image:MNcolon-matrix-1-color.png|frame|400px|center|Wybór bloku w macierzy <code>D</code> przy użyciu notacji dwukropkowej]]
 
 
Jest ona
 
szalenie intuicyjna, a mianowicie: pisząc na przykład <code>v <nowiki>=</nowiki>  D(2:5,2)</code>
 
definiujemy (wektor) <math>\displaystyle v</math>, który zawiera wiersze macierzy <math>\displaystyle D</math> od 2 do 5 wybrane z
 
drugiej kolumny. Podobnie, pisząc <code>W <nowiki>=</nowiki>  D(2:3,5:7)</code> definiujemy macierz
 
<math>\displaystyle W</math> wymiaru
 
<math>\displaystyle 2\times 3</math>, o elementach, które zostały wybrane z przecięcia wierszy od 2 do 3 z
 
kolumnami od 5 do 7 macierzy <math>\displaystyle D</math>.
 
 
[[Image:MNcolon-matrix-2-color.png|frame|400px|center|Wybór zestawu kolumn w macierzy <code>D</code>]]
 
 
Dodatkowo, aby odwołać się do całego 4. wiersza (odpowiednio: 5. kolumny)
 
macierzy <math>\displaystyle D</math>, można użyć skrótowej notacji <code>D(4,:)</code> (odpowiednio:
 
<code>D(:,5)</code>).
 
 
Odwołanie się do całej macierzy jest także możliwe (przez użycie jej nazwy, np.
 
<code>A <nowiki>=</nowiki>  2*D</code> lub <code>G <nowiki>=</nowiki>  log(abs(D))</code>).
 
 
Notację dwukropkową można także wykorzystać do wygenerowania samodzielnych zbiorów indeksów:
 
 
<div class="output" style="background-color:#e0e8e8; padding:1em"><pre>
 
 
octave:1> N <nowiki>=</nowiki>  10;
 
octave:2> idx <nowiki>=</nowiki>  1:N
 
idx <nowiki>=</nowiki>
 
  1  2  3  4  5  6  7  8  9  10
 
 
octave:3> idx2 <nowiki>=</nowiki>  1:2:N
 
idx2 <nowiki>=</nowiki>
 
  1  3  5  7  9
 
 
octave:4> nidx <nowiki>=</nowiki>  N:-1:1
 
nidx <nowiki>=</nowiki>
 
  10  9  8  7  6  5  4  3  2  1
 
</pre></div>
 
 
i dlatego pętle, które w C zapisywalibyśmy
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
for(i<nowiki>=</nowiki> 0; i<<nowiki>=</nowiki> N; i++)
 
{
 
...instrukcje...
 
}
 
for(i<nowiki>=</nowiki> N; i><nowiki>=</nowiki> 1; i--)
 
{
 
...instrukcje...
 
}
 
</pre></div>
 
 
w Octave zapiszemy
 
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
for i<nowiki>=</nowiki> 0:N
 
...instrukcje...
 
end
 
for i<nowiki>=</nowiki> N:-1:1
 
...instrukcje...
 
end
 
</pre></div>
 
 
Za chwilę przekonamy się, jak można w ogóle pozbyć się potrzeby stosowania większości pętli w kodach Octave i MATLABa.
 
 
====Wektoryzacja w Octave====
 
 
Ponieważ podstawowymi obiektami w Octave są wektory i macierze, predefiniowane
 
operacje matematyczne wykonują się od razu na całej macierzy. Bez żadnej
 
przesady możena stwierdzić, że umiejętność  <strong>wektoryzacji</strong> i
 
<strong>blokowania</strong> 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.
 
 
{| border=1
 
|+ <span style="font-variant:small-caps"> </span>
 
|-
 
|
 
  ||  Tradycyjny kod w Octave, używający pętli  ||  Efektywny kod wektorowy (macierzowy) w Octave
 
|-
 
|
 
 
Alg 1  ||
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
s <nowiki>=</nowiki>  0;
 
for i <nowiki>=</nowiki>  1:size(x,1)
 
s <nowiki>=</nowiki>  s + abs(x(i));
 
end
 
</pre></div>
 
 
||
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
s <nowiki>=</nowiki>  sum(abs(x));
 
</pre></div>
 
 
|-
 
|
 
 
Alg 2  ||
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
N <nowiki>=</nowiki>  500; h <nowiki>=</nowiki>  (b-a)/(N-1);
 
for i <nowiki>=</nowiki>  1:N
 
x(i) <nowiki>=</nowiki>  a + (i-1)*h;
 
y(i) <nowiki>=</nowiki>  sin(x(i));
 
end
 
plot(x,y);
 
</pre></div>
 
 
||
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
N <nowiki>=</nowiki>  500;
 
x <nowiki>=</nowiki>  linspace(a,b,N);
 
y <nowiki>=</nowiki>  sin(x);
 
plot(x,y);
 
</pre></div>
 
 
|-
 
|
 
Alg 3a  ||
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
for i <nowiki>=</nowiki>  1:size(C,1),
 
for j <nowiki>=</nowiki>  1:size(C,2),
 
for k <nowiki>=</nowiki>  1:size(A,2),
 
C(i,j) <nowiki>=</nowiki>  C(i,j) + A(i,k)*B(k,j);
 
end
 
end
 
end
 
</pre></div>
 
 
||
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
C <nowiki>=</nowiki>  C + A*B;
 
</pre></div>
 
 
|-
 
|
 
 
Alg 3b  ||
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
for i <nowiki>=</nowiki>  1:size(C,1),
 
for j <nowiki>=</nowiki>  1:size(C,2),
 
C(i,j) <nowiki>=</nowiki>  C(i,j) + A(i,:)*B(:,j);
 
end
 
end
 
</pre></div>
 
 
||
 
<div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
 
 
C <nowiki>=</nowiki>  C + A*B;
 
</pre></div>
 
 
|-
 
|
 
 
|}
 
 
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 <code>x</code> nie tylko gdy <code>x</code>
 
jest wektorem (co musieliśmy milcząco założyć w kodzie w lewej kolumnie), ale
 
tak samo dobrze zadziała, gdy <code>x</code> jest macierzą!
 
 
Szczególnie wiele na czytelności i wygodzie zyskujemy w drugim przykładzie,
 
gdzie najpierw funkcja <code>linspace</code> pozwala nam uniknąć żmudnego
 
wyznaczania <math>\displaystyle N</math> równoodległych węzłów <math>\displaystyle x_i</math> w odcinku <math>\displaystyle [a,b]</math>, a następnie
 
funkcja <code>sin</code> zaaplikowana do całego wektora <math>\displaystyle x</math> 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: (3a) i (3b), 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 (3a) w wersji z potrójną pętlą
 
<code>for</code> naśladuje sposób programowania znany nam z C lub Pascala,
 
natomiast przykład (3b) zdaje się być napisany odrobinę w duchu wektorowym
 
(brak trzeciej, wewnętrznej pętli, zastąpionej operacją wektorową:  iloczynem
 
skalarnym <math>\displaystyle i</math>-tego wiersza <math>\displaystyle A</math> i <math>\displaystyle j</math>-tej kolumny <math>\displaystyle B</math>). Poniżej porównanie
 
czasów działania tych trzech implementacji w przypadku macierzy <math>\displaystyle 64\times 64</math>
 
(czasy dla PC z procesorem Celeron 1GHz):
 
* Dla pętli postaci <code>C(i,j)<nowiki>=</nowiki> C(i,j)+A(i,k)*B(k,j)</code> uzyskano czas
 
21.6s,
 
* Dla pętli postaci  <code>C(i,j)<nowiki>=</nowiki> C(i,j)+A(i,:)*B(:,j)</code> --- 0.371s
 
* Dla pętli postaci <code>C<nowiki>=</nowiki> C+A*B</code> 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
 
<code>C <nowiki>=</nowiki>  C + A*B</code>. Po wektoryzacji wewnętrznej pętli, program doznaje
 
kilkudziesięciokrotnego przyspieszenia, lecz nadal jest ponadstukrotnie
 
wolniejszy od kodu macierzowego!
 

Aktualna wersja na dzień 11:28, 14 maj 2009