Elementy 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.
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 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 programu. Zastosujemy przy tym notację podobną do tej z języka programowania C.
Program składa się z deklaracji, czyli opisu obiektów, których będziemy używać, oraz z poleceń (instrukcji), czyli opisu akcji, które będziemy wykonywać.
Dostępnymi obiektami są stałe i zmienne typu całkowitego (int), rzeczywistego (float i double). Typ logiczny symulujemy tak jak w C wartościami zero-jedynkowymi typu całkowitego. Zmienne jednego typu mogą być grupowane w wektory albo tablice. Widzimy więc, że podstawowe algorytmy numeryczne będą bazować na mało skomplikowanych typach danych. Również nieskomplikowane będą instrukcje naszego modelowego języka. Dążenie do prostoty (ale nie za wszelką cenę) jest charakterystyczne dla numeryki. Typowe jest zastępowanie skomplikowanych tablic rekordów prostszymi macierzami, eleganckich rekurencji --- zwykłymi pętlami działającymi na danych w miejscu. Jedną z myśli przewodnich numeryki jest bowiem szybkość, a rezygnacja z barokowych konstrukcji językowych jest zgodna z filozofią architektury procesora RISC: efektywność przez prostotę.
(Na pewno zastanawiasz się teraz, jaka jest druga myśl przewodnia numeryki. To dokładność.)
Podstawienie
Najprostsza z instrukcji, bez której nie można się obejść:
=
Tutaj
jest zmienną, a
jest wyrażeniem o wartościach tego samego typu co
. 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
,
,
,
,
, gdzie
są stałymi lub zmiennymi liczbowymi,
arytmetyczno--logiczne
,
,
,
, gdzie
są stałymi lub zmiennymi liczbowymi,
logiczno--logiczne
,
,
, gdzie
są stałymi lub zmiennymi logicznymi.
Dla niektórych zadań wygodnie jest uzupełnić ten zbiór o dodatkowe operacje, takie jak obliczanie wartości niektórych standardowych funkcji matematycznych (
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 przy użyciu czterech podstawowych operacji arytmetycznych. Dokładniej, jednostka arytmetyczna procesora potrafi wykonywać jedynie operacje
, 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 (
itp.), jednak wykonanie takiej instrukcji wymaga mniej więcej kilkadziesiąt 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 nawet podstawowe działania arytmentyczne (i, tym bardziej, obliczanie wartości funkcji matematycznych) nie są wykonywane dokładnie. Czasem uwzględnienie tego faktu wiąże się ze znaczącym wzrostem komplikacji analizy algorytmu i dlatego "w pierwszym przybliżeniu" często pomija się to ograniczenie przyjmując model, w którym wszystkie (lub prawie wszystkie) działania arytmetyczne wykonują się dokładnie. Wiedza o tym, kiedy i jak zrobić to tak, by wciąż wyciągać prawidłowe wnioski odnośnie faktycznej realizacji algorytmów w obecności błędów zaokrągleń jest częścią sztuki i wymaga intuicji numerycznej, popartej doświadczeniem.
Mamy trzy podstawowe polecenia złożone: warunkowe, powtarzania i kombinowane.
Warunkowe
if
;
else
;
gdzie
jest wyrażeniem o wartościach całkowitych (0 odpowiada logicznemu fałszowi, inne wartości --- logicznej prawdzie), a
i
są poleceniami, przy czym dopuszczamy polecenia puste.
Powtarzane
while
;
gdzie
jest wyrażeniem o wartościach logicznych, a
jest poleceniem.
Kombinowane
{
;
;An
}
gdzie
są poleceniami.
Na podstawie tych trzech poleceń można tworzyć inne, takie jak pętle for(), czy instrukcje wariantowe switch(), itd.
Mamy też dwa szczególne polecenia, które odpowiadają za "wejście" i "wyjście".
Wprowadzanie danych
;
gdzie
jest zmienną rzeczywistą, a
"adresem" pewnego funkcjonału
należącym to pewnego zbioru T. W wyniku wykonania tego polecenia w zmiennej x zostaje umieszczona wartość Lt(f).
Polecenie to pozwala zdobyć informację o danej
. Jeśli
to zwykle mamy
i
, co w praktyce odpowiada wczytaniu
-tej współrzędnej wektora danych. W szczególności, ciąg poleceń
,
, pozwala uzyskać pełną informację o
. Jeśli zaś
jest pewną klasą funkcji
, to możemy mieć np.
i
. W tym przypadku, wykonanie polecenia
odpowiada w praktyce skorzystaniu ze specjalnej procedury (albo urządzenia zewnętrznego) obliczającej (mierzącego) wartość funkcji
w punkcie
.
Wyprowadzanie wyników
W;
gdzie
jest wyrażeniem o wartościach rzeczywistych. Polecenie to pozwala "wskazać" kolejną współrzędną wyniku.
Zakładamy, że na początku procesu obliczeniowego wartości wszystkich zmiennych są nieokreślone, oraz że dla dowolnych danych wykonanie programu wymaga wykonania skończonej liczby operacji elementarnych. Wynikiem obliczeń jest skończony ciąg liczb rzeczywistych (albo
), którego kolejne współrzędne pokazywane są poleceniem
.
Środowisko obliczeniowe
Ze względu na swój utylitarny charakter, w metodach numerycznych niemal równie ważna jak dobór optymalnego algorytmu jest jego efektywna implementacja na konkretnej architekturze.
W praktyce mamy dwie możliwości:
wykorzystanie standardowych języków programowania (C, Fortran, być może ze wstawkami w asemblerze) oraz 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: deweloperka 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
używać tak prostych struktur danych, jak to możliwe (ale nie prostszych!...)
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 teraz coraz więcej oprogramowania numerycznego powstaje w C.
W naszym wykładzie wybieramy C ze względu na jego uniwersalność, doskonałą przenośność i całkiem dojrzałe kompilatory. Dodajmy, że funkcje w C można mieszać 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 i inne, nie są zbyt odpowiednie dla obliczeń numerycznych. Mało tego, np. podstawowy typ numeryczny Pascala: real nie jest zgodny z powszechnym standardem IEEE 754. Jednak, ze względu na coraz większą komplikację kodów numerycznych służących np. do prowadzenia zaawansowanych symulacji metodą elementu skończonego, coraz więcej aplikacji numerycznych 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. W przykładach będziemy korzystać z kompilatora GCC, który jest omówiony w wykładzie Środowisko programisty.
Prosty program numeryczny w C
Napiszemy teraz program obliczający (w niezbyt wyrafinowany sposób)
-tą sumę częściową szeregu harmonicznego
Przyjmijmy, że parametr
będzie miał wartość równą 2006. W pierwszym odruchu, prawie każdy pisze program w rodzaju:
#include <stdio.h>
#define N 2006
int main(void)
{
float x;
unsigned int i;
x = 0.0;
for(i = 1; i <= N; i++)
x = x + 1/i;
printf("Wartość sumy x = 1 + 1/2 + ... + 1/%d jest równa %g\n", N, x);
return(0);
}
Sęk w tym, że ten program nie działa! To znaczy: owszem, kompiluje się i uruchamia, ale twierdzi uparcie, że nasza suma wynosi... 1.
Winę za to ponosi linijka
x = x + 1/i;
w której wykonujemy dzielenie 1/i. Obie liczby są typu int i dlatego, zgodnie z regułami C, wynik ich dzielenia także jest całkowity: dostajemy część całkowitą z dzielenia tych dwóch liczb, czyli, gdy tylko
, po prostu zero.
Prawidłowy program musi więc wymusić potraktowanie choć jednej z tych liczb jako liczby zmiennoprzecinkowej, co najprościej uzyskać albo przez
x = x + 1.0/i;
albo bardziej uniwersalnie, rzutując choć jedną z liczb na odpowiedni typ:
x = x + 1/((float) i);
Poprawny kod miałby więc postać
#include <stdio.h>
#define N 2006
int main(void)
{
float x;
unsigned int i;
x = 0.0;
for(i = 1; i <= N; i++)
x = x + 1.0/i;
printf("Wartość sumy x = 1 + 1/2 + ... + 1/%d jest równa %g\n", N, x);
return(0);
}
W tym krótkim przykładzie zetknęliśmy się z największą trudnością pisania programów wykonujących obliczenia numeryczne:
Bardzo łatwo jest napisać jakiś program, który działa i generuje ładnie wyglądające, choć całkowicie fałszywe wyniki.
Pisanie programów numerycznych (podobnie jak innych), wymaga żelaznej dyscypliny i precyzji. Wielokrotne testy, sprawdzające nie tylko kilka konkretnych danych ("czy działa dla
?"), ale także spodziewane jakościowe własności produkowanych rezultatów, są smutną koniecznością w realnym życiu numeryka-programisty.
Typy numeryczne w C
W języku C mamy dostępnych sześć typów numerycznych:
stałoprzecinkowe, dla reprezentacji liczb całkowitych
int oraz long int. W realizacji GCC na komputery x86, oba typy: int i long int są identyczne (32-bitowe) i ich zakres wynosi około
. Typ int i jemu pokrewne odnoszą się do liczb całkowitych ze znakiem (dodatnich lub ujemnych). Ich warianty bez znaku: unsigned int, itp. odnoszą się do liczb bez znaku (nieujemnych), dlatego np. zakresem unsigned int będzie w przybliżeniu
.long long int (64-bitowy) o zakresie w przybliżeniu
.zmiennoprzecinkowe, dla reprezentacji liczb rzeczywistych
float, pojedynczej precyzji, 32-bitowy, gwarantuje precyzję około
i zakres liczb reprezentowalnych w przybliżeniu
;double, podwójnej precyzji, 64-bitowy, ma precyzję na poziomie
przy orientacyjnym zakresie
;long double, rozszerzonej podwójnej precyzji, 80-bitowy, o precyzji rzędu
i odpowiadający standardowi double extended IEEE 754.
Powyższe typy zmiennoprzecinkowe w realizacji na PC odpowiadają standardowi IEEE 754, o którym będzie mowa w wykładzie o arytmetyce zmiennoprzecinkowej. Standardowo, operacje arytmetyczne na obu typach float i double są tak samo pracochłonne, gdyż wszystkie obliczenia w C wykonywane są z maksymalną dostępną precyzją (czyli, na procesorach architektury IA-32 Intela i AMD: w precyzji oferowanej przez typ long double), a następnie dopiero wynik zapisywany do zmiennej rzutowany jest na stosowny typ.
Mimo to, czasem można skorzystać z typu pojedynczej precyzji float, który oferuje znacznie większe możliwości optymalizacji uzyskanego kodu. Ponadto, obiekty typu float zajmują dwukrotnie mniej miejsca w pamięci niż double, dając możliwość lepszego wykorzystania pamięci podręcznej (cache) i przetwarzania wektorowego.
Stałe matematyczne i podstawowa biblioteka matematyczna
Język C jest językiem "małym" i, jak wiadomo, nawet proste operacje wejścia-wyjścia są w istocie nie częścią języka, ale funkcjami bibliotecznymi lub makrami. 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 (
, itp.) nie są składnikami języka C, lecz w zamian są zaimplementowane w tzw. standardowej bibliotece matematycznej libm.a; prototypy tych funkcji oraz definicje rozmaitych stałych matematycznych:
itp. znajdują się w pliku nagłówkowym math.h. Aby więc skorzystać z tych funkcji w programie, należy
w nagłówku pliku, w którym korzystamy z funkcji lub stałych matematycznych, umieścić linię #include <math.h>
przy linkowaniu dołączyć bibliotekę matematyczną za pomocą opcji -lm
Przykład
Oto prosty program numeryczny w C; drukuje on tablicę wartości sinusów losowo wybranych liczb z przedziału
.
#include <math.h>
#include <stdio.h>
#include <stdlib.h> /* zawiera definicję funkcji rand() i stałej RAND_MAX */
#define N 15 /* ile liczb wydrukować */
int main(void)
{
int i;
double x,y;
for( i = 0; i < N; i++)
{
x = rand()/(double)RAND_MAX;
x *= 2.0*M_PI;
/* oczywiście, wystarczyłoby x =(2.0*M_PI*rand())/RAND_MAX; */
y = sin(x);
fprintf(stderr, "(%3d) x = %10.5e sin(x) = %10.5e\n", i, x, y);
}
return(0);
}
Zwróćmy uwagę na linię x = rand()/(double)RAND_MAX; Funkcja rand() zwraca losową liczbę całkowitą z przedziału [0,RAND_MAX], więc iloraz z predefiniowaną stałą RAND_MAX będzie liczbą z przedziału
. Dla prawidłowego uzyskania losowej liczby z przedziału
kluczowe jest jednak zrzutowanie jednej z dzielonych liczb na typ double! Gdyby tego nie zrobić, uzyskalibyśmy zawsze
lub sporadycznie
, zgodnie z regułą C: "typ wyniku jest zgodny z typem argumentów".
Kompilujemy ten program, zgodnie z uwagami uczynionymi na początku, poleceniem
gcc -o sinusy sinusy.c -lm
Środowisko obliczeń numerycznych:
Stworzone przez Cleve Moler: 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ą MATLAB (od: "MATrix LABoratory").
Typowa sesja MATLABa. Zwróć uwagę na edytor kodu źródłowego na bieżąco interpretujący go i wychwytujący potencjalne błędy
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).
Kierując się podobnymi przesłankami co Moler, oraz bazując na wielkim sukcesie MATLABa, John W. Eaton z Wydziału Inżynierii Chemicznej Uniwersytetu Wisconsin w USA, zaczął w 1994 roku opracowywać darmowe (udostępniane na tzw. licencji GPL) oprogramowanie o funkcjonalności maksymalnie zbliżonej do MATLABa: Octave. Wersja 1.0 pakietu Octave ukazała się w 1996 roku i jest intensywnie rozwijana do dziś.
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.
Literatura
D. Kincaid, W. Cheney, Analiza numeryczna, Wydawnictwa Naukowo-Techniczne, Warszawa 2006. ISBN 83-204-3078-X
A. Björck, G. Dahlquist, Metody numeryczne, Państwowe Wydawnictwo Naukowe, Warszawa 1987. ISBN 83-01-04276-1
Wyszukiwarka