|
|
| 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/%d jest równa %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/%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/ 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, "(%3d) x <nowiki>=</nowiki> %10.5e sin(x) <nowiki>=</nowiki> %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. 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: [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!
| |
