|
|
| (Nie pokazano 12 pośrednich wersji utworzonych przez tego samego użytkownika) |
| Linia 1: |
Linia 1: |
|
| |
|
| ==Wprowadzenie do metod numerycznych==
| |
|
| |
| ''Metody numeryczne'' to dziedzina wiedzy zajmująca się
| |
| problemami obliczeniowymi i konstrukcją algorytmów rozwiązywania
| |
| zadań matematycznych. Najczęściej, zadania obliczeniowe postawione są w
| |
| dziedzinie rzeczywistej (lub zespolonej) i dlatego mówimy o zadaniach
| |
| obliczeniowych ''matematyki ciągłej'' (w odróżnieniu od matematyki
| |
| dyskretnej).
| |
|
| |
| ===Zadania metod numerycznych===
| |
|
| |
| Aby w ogóle mówić w problemie obliczeniowym, musimy najpierw
| |
| * określić ''dane problemu'' i ''cel obliczeń'', czyli dokładnie
| |
| sformułować zadanie w języku matematyki,
| |
| * określić ''środki obliczeniowe'' dzięki którym chcemy osiągnąć cel,
| |
| * dla analizy zadania i sposobów jego rozwiązania wygodnie jest zdefiniować
| |
| ''klasę rozpatrywanych danych'' oraz ''model obliczeniowy'' w obrębie
| |
| którego będą działać nasze algorytmy.
| |
|
| |
| Wbrew dość powszechnej opinii ''nie jest prawdą'', że głównym celem metod
| |
| numerycznych jest badanie wpływu błędów zaokrągleń na wynik. Raczej, głównym
| |
| celem metod numerycznych jest konstrukcja optymalnych (w jasno określonym
| |
| sensie, np. pod względem wymaganej liczby operacji, lub pod względem ilości
| |
| niezbędnej informacji, czy też pod względem dokładności uzyskiwanego wyniku)
| |
| algorytmów rozwiązywania konkretnych zadań matematycznych.
| |
|
| |
| {{uwaga|||
| |
| Nasz przedmiot ma różne wcielenia i z tego powodu czasem nosi inne nazwy, w
| |
| zależności od tego, na jaki aspekt metod obliczeniowych jest położony największy
| |
| nacisk.
| |
|
| |
| ; metody numeryczne
| |
| : --- główny nacisk idzie na aspekty algorytmiczne, ale przy uwzględnieniu pozostałych;
| |
|
| |
| ; analiza numeryczna
| |
| : --- przede wszystkim badanie właściwości algorytmów, ich optymalności oraz wpływu arytmetyki zmiennopozycyjnej na jakość uzyskanych wyników;
| |
|
| |
| ; matematyka obliczeniowa
| |
| : --- głównie teoretyczna analiza możliwości taniej i dokładnej aproksymacji rozwiązań zadań matematycznych;
| |
|
| |
| ; obliczenia naukowe
| |
| : --- nacisk na praktyczne zastosowania metod numerycznych, symulacje, realizacje na komputerach o dużej mocy obliczeniowej.
| |
|
| |
| Oczywiście, granice podziału nie są ostre i najczęściej typowy wykład z tego
| |
| przedmiotu stara się pokazać pełne spektrum zagadnień z nim związanych. Tak
| |
| będzie również i w naszym przypadku.
| |
|
| |
| }}
| |
|
| |
| ===Model obliczeniowy===
| |
|
| |
| Tworząc i analizując algorytmy, jakie będą pojawiać w naszym wykładzie, będziemy
| |
| posługiwać się pewnym uproszczonym modelem obliczeń, dzięki czemu będziemy mogli
| |
| skoncentrować się na ''esencji'' algorytmu, bez detali implementacyjnych ---
| |
| zostawiając je na inne okazje.
| |
|
| |
| Aby zdefiniować nasz model obliczeniowy, posłużymy się
| |
| pojęciem ''programu''. Zastosujemy przy tym notację
| |
| podobną do tej z języka programowania ''C''.
| |
|
| |
| Program składa się z ''deklaracji'', czyli opisu obiektów,
| |
| których będziemy używać, oraz z ''poleceń'' (''instrukcji''),
| |
| czyli opisu akcji, które będziemy wykonywać.
| |
|
| |
| Dostępnymi obiektami są ''stałe'' i ''zmienne'' typu
| |
| całkowitego (<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 ''wyrażeniem''
| |
| 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 ''operacji elementarnych'' 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\,\tilde{}\,p</math>,
| |
| <math>\displaystyle (p,q)\mapsto p\,and\,q</math>, <math>\displaystyle (p,q)\mapsto p\,or\,q</math>,
| |
| gdzie <math>\displaystyle p,q</math> są stałymi lub zmiennymi logicznymi.
| |
|
| |
| Dla niektórych zadań wygodnie jest (a czasem koniecznie)
| |
| uzupełnić ten zbiór o dodatkowe operacje, takie jak
| |
| obliczanie wartości niektórych standardowych funkcji matematycznych
| |
| (<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.
| |
| }}
| |
|
| |
| Niestety, aby nasz model obliczeniowy wiernie odpowiadał rzeczywistości, musimy
| |
| w nim uwzględnić fakt, że działania matematyczne (ani tym bardziej
| |
| obliczanie wartości funkcji matematycznych) ''nie są'' wykonywane dokładnie.
| |
| Czasem
| |
| uwzględnienie tego faktu wiąże się ze znaczącym wzrostem komplikacji analizy
| |
| algorytmu i dlatego "w pierwszym przybliżeniu" często pomija się to
| |
| ograniczenie przyjmując model w którym wszystkie (lub prawie wszystkie)
| |
| działania arytmetyczne wykonują się dokładnie. Wiedza o tym, ''kiedy'' i
| |
| ''jak'' zrobić to tak, by wciąż wyciągać prawidłowe wnioski odnośnie faktycznej
| |
| realizacji algorytmów w obecności błędów zaokrągleń jest częścią sztuki i wymaga
| |
| intuicji numerycznej, popartej doświadczeniem.
| |
|
| |
| Mamy trzy podstawowe polecenia złożone: warunkowe, powtarzania i kombinowane.
| |
|
| |
| ====Warunkowe====
| |
|
| |
| <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ć ''informację'' 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 ''implementacja'' na
| |
| konkretnej architekturze.
| |
|
| |
| W praktyce mamy dwie możliwości:
| |
| * wykorzystanie standardowych języków programowania (C, Fortran, być może ze wstawkami w asemblerze) oraz specjalizowanych bibliotek
| |
| * użycie gotowego środowiska obliczeń numerycznych będącego wygodnym interfejsem do specjalizowanych bibliotek numerycznych
| |
|
| |
| ====Języki programowania: C i Fortran====
| |
|
| |
| Programy numeryczne (a przynajmniej ich jądra obliczeniowe) są zazwyczaj niezbyt
| |
| wymagające jeśli chodzi o struktury danych, co więcej, prostota struktur danych
| |
| szybko rewanżuje się efektywniejszym kodem. Dlatego, trawestując Einsteina, w
| |
| dobrym programie numerycznym należy ''używać tak prostych struktur danych, jak to
| |
| możliwe (ale nie prostszych!''...)
| |
|
| |
| Językami programowania opartymi na prostych konstrukcjach są: Fortran i C. Dlatego właśnie są to języki
| |
| dominujące współcześnie pisane programy numeryczne. O ile w przeszłości
| |
| hegemonia Fortranu była nie do podważenia, o tyle w chwili obecnej coraz więcej
| |
| oprogramowania numerycznego powstaje w C.
| |
|
| |
| W naszym wykładzie wybieramy C ze względu na jego uniwersalność,
| |
| doskonałą przenośność i w chwili obecnej całkiem dojrzałe kompilatory. Dodajmy,
| |
| że funkcje w C można mieszać z np. z gotowymi bibliotekami napisanymi w
| |
| Fortranie. Fortran,
| |
| język o bardzo długiej tradycji, wciąż żywy i mający grono wiernych fanów, jest
| |
| nadal wybierany przez numeryków na całym świecie między innymi ze względu na
| |
| jego dopasowanie do zadań obliczeniowych (właśnie w tym celu powstał), a także
| |
| ze względu na doskonałe kompilatory dostępne na superkomputerach, będące efektem
| |
| wieloletniej ewolucji i coraz lepszego nie tylko dopasowania kompilatora do
| |
| spotykanych konstrukcji językowych, ale także na odwrót --- coraz lepszego
| |
| zrozumienia programistów, jak pisać programy, by wycisnąć jak najwięcej z
| |
| kompilatora.
| |
|
| |
| Inne popularne języki: Java, Pascal (ten język, zdaje się, jest popularny już
| |
| tylko w obrębie wydziału MIM UW...), VisualBasic nie są zbyt odpowiednie dla
| |
| obliczeń numerycznych. Mało tego, np. podstawowy typ numeryczny Pascala:
| |
| <code>real</code> nie jest zgodny z powszechnym standardem IEEE 754. Jednak, ze względu na coraz większą
| |
| komplikację kodów numerycznych służących np. do prowadzenia zaawansowanych
| |
| symulacji metodą elementu skończonego, coraz więcej kodów wykorzystuje
| |
| możliwości obiektowych języków C++ i Fortran90.
| |
|
| |
| W przykładach będziemy najczęściej odnosić się do architektury x86, tzn. 32-bitowej IA-32
| |
| procesorów firmy Intel i AMD, najczęściej spotykanej w obecnie używanych
| |
| komputerach. Należy jednak pamiętać, że obecnie następuje przejście na
| |
| architekturę 64-bitową. Ze względu jednak na brak pewności co do ostatecznie
| |
| przyjętych standardów w tym obszarze, ograniczymy się do procesorów
| |
| 32-bitowych.
| |
|
| |
| ====Prosty program numeryczny====
| |
|
| |
| Napiszemy teraz program obliczający (w niezbyt wyrafinowany sposób) <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
| |
| return(0);
| |
| }
| |
| </pre></div>
| |
|
| |
| Sęk w tym, że ten program ''nie działa!'' To znaczy, kompiluje się i
| |
| uruchamia, ale twierdzi uparcie, że nasza suma wynosi... 1.
| |
|
| |
| Winę za to ponosi linijka
| |
| <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>
| |
|
| |
| 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 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 ''typ
| |
| wyniku jest zgodny z typem argumentów''. 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|thumb|400px||Typowa sesja MATLABa. Zwróć uwagę na edytor kodu źródłowego na
| |
| bieżąco interpretujący go i wychwytujący potencjalne błędy]]
| |
|
| |
| W końcu lat siedemdziesiątych ubiegłego wieku, Cleve Moler wpadł na pomysł
| |
| stworzenia prostego interfejsu do ówcześnie istniejących bibliotek numerycznych
| |
| algebry liniowej: pakietów LINPACK i EISPACK
| |
|
| |
| (których był współautorem). Stworzone przez niego:
| |
| język skryptów i łatwe w użyciu narzędzia manipulacji macierzami, zdobyły
| |
| wielką popularność w środowisku naukowym. W 1984 roku Moler skomercjalizował swe
| |
| dzieło, zakładając wspólnie z J. Little'm firmę MathWorks
| |
| .
| |
| Jej główny produkt nazwano [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|thumb|400px||Screenshot Octave. W terminalu otwarta sesja w
| |
| ascetycznym
| |
| trybie tekstowym, grafika wyświetlana z wykorzystaniem Gnuplota]]
| |
|
| |
| Kierując się podobnymi przesłankami co C. Moler, oraz bazując na wielkim
| |
| sukcesie MATLABa, John W. Eaton z Wydziału Inżynierii Chemicznej Uniwersytetu
| |
| Wisconsin w USA, zaczął w 1994 roku opracowywać darmowe (udostępniane na tzw.
| |
| licencji GPL) oprogramowanie o funkcjonalności maksymalnie
| |
| zbliżonej do MATLABa: [http://www.octave.org Octave]. Jak wyjaśnia twórca Octave,
| |
| , nazwa nie ma nic wspólnego z muzyką, ale z
| |
| pewną anegdotą. Wersja 1.0 pakietu Octave ukazała się w 1996 roku i jest
| |
| intensywnie rozwijana do dziś.
| |
|
| |
| [[Image:MNscilab-screenshot.png|thumb|400px||Screenshot Scilaba]]
| |
|
| |
| Drugim udanym klonem MATLABa jest francuski [http://www.scilab.org Scilab] ,
| |
| opracowany w laboratoriach INRIA i wciąż doskonalony. W
| |
| subiektywnej ocenie autora niniejszych notatek, nie dorównuje on elegancji i
| |
| szlachetnej prostocie Octave. Na plus tego pakietu należy zaliczyć m.in.
| |
| znacznie bardziej rozbudowany zestaw funkcji podstawowych, na minus -- przede
| |
| wszystkim znacznie mniejszy stopień zgodności z MATLABem, a poza tym: niewygodny
| |
| system pomocy oraz "toporną" (choć o dużym potencjale) grafikę.
| |
|
| |
| ====Porównanie MATLABa, Octave i Scilaba====
| |
|
| |
| Korzystając z omówionych powyżej pakietów, otrzymujemy:
| |
| * możliwość obliczania funkcji matematycznych (nawet dość egzotycznych),
| |
| finansowych, analizy sygnałów, itp.;
| |
| * bardzo szeroki zakres nowoczesnych narzędzi umożliwiających
| |
| wykonywanie podstawowych zadań numerycznych, takich jak: rozwiązywanie równań,
| |
| obliczanie całek, itd.;
| |
| * efektowną wizualizację wyników w postaci wykresów dwu- i trójwymiarowych,
| |
| a także funkcje do manipulacji obrazem i dźwiękiem;
| |
| * możliwość nieograniczonych rozszerzeń przy użyciu funkcji i skryptów
| |
| tworzonych przez osoby trzecie (lub własnych).
| |
|
| |
| Jeśli uwzględnić dodatki wielce ułatwiające życie, w tym rozbudowane funkcje
| |
| wejścia-wyjścia, to narzędzie robi się rzeczywiście intrygujące, jako środowisko
| |
| obliczeń numerycznych "dla każdego" (nie tylko dla profesjonalnego numeryka).
| |
| Stąd wielka popularność MATLABa w środowiskach inżynierkich, gdyż umożliwia
| |
| szybką implementację rozmaitych --- zwłaszcza testowych --- wersji algorytmu np.
| |
| przeprowadzania skomplikowanej symulacji. Po zweryfikowaniu wyników, można
| |
| ewentualnie rozpocząć implementację w docelowym środowisku (np. masywnie
| |
| równoległym komputerze, w Fortranie 95 z wykorzystaniem komercyjnych bibliotek i
| |
| specyficznymi optymalizacjami kodu), choć bardzo często wyniki uzyskane w
| |
| MATLABie w zupełności wystarczą.
| |
|
| |
| Zarówno MATLAB, jak Octave i Scilab, opierają się w dużym zakresie na darmowych
| |
| --- acz profesjonalnych --- zewnętrznych bibliotekach numerycznych (LAPACK,
| |
| ATLAS, ARPACK, itp.)
| |
|
| |
| MATLAB jest wciąż niedościgniony jeśli chodzi o narzędzia programistyczne (np.
| |
| ma wbudowany --- działający ''na bieżąco'' w czasie edycji kodu źródłowego! ---
| |
| debugger) oraz możliwości wizualizacji danych i wyników obliczeń. Jeśli nie dbać
| |
| o koszta zakupu dodatkowych modułów, ma także najbogatszy zestaw funkcji
| |
| numerycznych. W chwili pisania niniejszego tekstu, (stosunkowo tanie) wersje
| |
| studenckie MATLABa, nie są oferowane w Polsce przez producenta, choć z drugiej
| |
| strony są powszechnie dostępne (nawet w wersji profesjonalnej) w komputerowych
| |
| laboratoriach akademickich w Polsce.
| |
|
| |
| Wybór pomiędzy darmowymi klonami MATLABa: Octave i Scilabem, jest właściwie kwestią osobistych preferencji.
| |
| W ninijeszym kursie będziemy posłygiwać się
| |
| Octave, który jest bardzo wiernym klonem MATLABa, przez co daje zaskakująco wysoki
| |
| stopień kompatybilności z programami napisanymi w MATLABie. To jest bardzo
| |
| ważne w praktyce, ponieważ jest dostępnych wiele (nawet darmowych) pakietów
| |
| numerycznych opracowanych właśnie w języku poleceń MATLABa. Umiejąc posługiwać
| |
| się Octave, bez trudu można "przesiąść się" na MATLABa i korzystać z jego
| |
| bogatszych możliwości, droga w przeciwną stronę też jest nietrudna.
| |
|
| |
| Inną zaletą Octave jest swobodniejsza niż w MATLABie składnia: przykładowo, w
| |
| Octave zarówno <code>'Newton'</code> jak i <code>"Newton"</code> oznacza łańcuch
| |
| znakowy (string) --- w MATLABie możliwa jest tylko pierwsza wersja; w
| |
| Octave możemy pisać pętle <code>for ... endfor</code> jak i <code>for ... end</code>
| |
| --- w MATLABie możliwa jest tylko druga postać; w Octave możemy pisać zarówno
| |
| <code>z <nowiki>=</nowiki> z+5</code> jak i, wzorem języka C, <code>z +<nowiki>=</nowiki> 5</code> --- MATLAB nie
| |
| dopuszcza drugiego z wariantów; znakiem komentarza w Octave może być <code>#</code>
| |
| oraz używany w tym celu MATLABie znak <code>definiowanie funkcji w trybie interaktywnym, natomiast funkcje MATLABa muszą
| |
| być zapisane w zewnętrznym pliku tekstowym.
| |
|
| |
| Zarówno w MATLABie, jak i w Octave, jest możliwe pisanie funkcji, które będą
| |
| prekompilowane, dając znaczące przyspieszenie działania.
| |
|
| |
| Octave jest dołączany do większości popularnych dystrybucji Linuxa, najczęściej
| |
| jednak użytkownik musi samodzielnie doinstalować go z płytki instalacyjnej.
| |
| Ponadto, kody źródłowe najświeższej wersji Octave są dostępne na stronie
| |
| http://www.octave.org. Dodatkowo, pod http://octave.sourceforge.net znajdziemy
| |
| pakiet rozszerzeń do Octave, pod nazwą Octave-forge. Wymaga on od użytkowników
| |
| Linuxa samodzielnej (przebiegającej bezproblemowo) instalacji. Octave można
| |
| także zainstalować pod Windows, korzystając z programu instalacyjnego
| |
| dostępnego pod adresem internetowym http://octave.sourceforge.net. W
| |
| Windowsowej wersji Octave, pakiet Octave-forge jest standardowo dołączony.
| |
|
| |
| ====Pierwsze kroki w Octave====
| |
|
| |
| Octave uruchamiamy poleceniem <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ą
| |
| ''numerycznie'', ze skończoną precyzją: w Octave niektóre tożsamości matematyczne
| |
| są spełnione ''jedynie w przybliżeniu'', 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>.
| |
|
| |
| ====Dokumentacja Octave====
| |
|
| |
| Octave ma dobrą dokumentację, dostępną w trakcie sesji
| |
| Octave; dla każdej funkcji, stałej itp. można uzyskać jej dobry opis przy
| |
| użyciu polecenia <code>help</code>.
| |
|
| |
| ====Macierz: postawowy obiekt w Octave====
| |
|
| |
| 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 ''docelowego wymiaru'' (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> ''bez'' ś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
| |
| : Zapis do formatu ASCII. Ze względu na późniejszą konwersję z zapisu w
| |
| postaci dzisiętnej do dwójkowej, zob. \link{xxx}, wiąże się z drobnymi
| |
| niedokładnościami w odczycie.
| |
| <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>
| |
|
| |
| ; binarny
| |
| : Po prostu zrzut stanu pamięci
| |
| <div class="code" style="background-color:#e8e8e8; padding:1em"><pre>
| |
|
| |
| save -binary macierzb.dat B
| |
| </pre></div>
| |
|
| |
| ; HDF5
| |
| : Format gwarantujący pełną przenośność danych pomiędzy różnymi
| |
| architekturami, akceptowany przez wiele zewnętrznych aplikacji, m.in. 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>.
| |
|
| |
| ====Wektoryzacja====
| |
|
| |
| Ponieważ podstawowym obiektem w Octave jest macierz, predefiniowane operacje matematyczne
| |
| wykonują się od razu na całej macierzy.
| |
| Bez żadnej przesady możena stwierdzić, że umiejętność ''wektoryzacji'' i ''blokowania'' algorytmów jest podstawą pisania efektywnych
| |
| implementacji algorytmów w Octave.
| |
|
| |
| Zobaczmy kilka prostych przykładów zawartych w tabeli poniżej. W
| |
| pierwszej kolumnie tabeli, dane zadanie zaimplementujemy w Octave przy użyciu
| |
| pętli (zwróćmy przy okazji uwagę na to, jak programuje się pętle w Octave, one
| |
| nie są zupełnie bezużyteczne...). W drugiej kolumnie zobaczymy, jak to samo
| |
| zadanie można wykonać korzystając z operatorów lub funkcji macierzowych.
| |
|
| |
| {| border=1
| |
| |+ <span style="font-variant:small-caps"> </span>
| |
| |-
| |
| |
| |
| || Tradycyjny kod w Octave, używający pętli || Efektywny kod wektorowy (macierzowy) w Octave
| |
| |-
| |
| |
| |
|
| |
| 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>
| |
|
| |
| |-
| |
| |
| |
|
| |
| 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>
| |
|
| |
| |-
| |
| |
| |
| 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>
| |
|
| |
| |-
| |
| |
| |
|
| |
| 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: trzeci i czwarty, które pokażą nam prawdziwą moc funkcji
| |
| macierzowych, unikających wielokrotnie zagnieżdżonych pętli. Oba dotyczą
| |
| operacji mnożenia dwóch macierzy. Przykład trzeci w wersji z potrójną pętlą
| |
| <code>for</code> naśladuje sposób programowania znany nam z C lub Pascala,
| |
| natomiast przykład czwarty zdaje się być napisany odrobinę w duchu wektorowym
| |
| (brak trzeciej, wewnętrznej pętli, zastąpionej operacją wektorową: iloczynem
| |
| skalarnym <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!
| |
