Metody Numeryczne w Fizyce 1 Błędy obliczeń numerycznych 1
Wprowadzenie, błędy obliczeń numerycznych
Błędy
Sposób zapisu liczb
Na początek zajmiemy się sposobem zapisu liczb w pamięci komputera. Wszyscy oczywiście wiedzą, że komputer
posługuje się zapisem w systemie dwójkowym, jednak aby bliżej zająć się tym zagadnieniem potrzebujemy więcej
szczegółów. Zajmiemy się wymyślonym komputerem prowadzącym obliczenia na liczbach rzeczywistych i przeznaczającym
na każdą liczbę 1 bajt (8 bitów)  jest więc to komputer ośmio bitowy. W rzeczywistości np. w Turbo Pascalu liczba typu
real zajmuje 6 bajtów (48 bitów). Bity słowa naszego wymyślonego komputera są interpretowane następująco:
1  bit znaku liczby  wartość 0 oznacza plus, wartość 1 oznacza minus.
2,3,4,5  bity mantysy (określenie to nie ma właściwie nic wspólnego ze znaną nam mantysą logarytmu).
Interpretuje się je jako cztery cyfry po kropce  dziesiętnej (w układzie dwójkowym należałoby ją może
nazwać  kropką dwójkową ).
6  bit znaku wykładnika, interpretacja taka sama jak bitu znaku liczby.
7,8  bity wykładnika (zwanego też czasem  cechą , podobnie jak w przypadku logarytmów). Cyfry te
interpretuje się jako liczbę całkowitą.
W ten sposób bajt postaci 11111111 (zazwyczaj będziemy to zapisywać biorąc bity znaków w nawiasy: (1)1111(1)11
 jest to ten sam bajt, tylko Å‚atwiejszy do zinterpretowania  na oko ) zostanie zinterpretowany jako: 0.11112 2 112 (tu
i dalej dolne indeksy 2 będą oznaczać, że mamy do czynienia z liczbą zapisaną w układzie dwójkowym)
1 1 1 1 15
0.11112 2 1 2 2 2 3 2 4 0.9375
2 4 8 16 16
1
112 21 20 3 ; 2 3
8
1
0.9375 0.1171875
8
Jest to więc  0.1171875.
Rozważmy teraz dwie bliskie sobie liczby: (1)1101 (0)10 oraz (1)1100 (0)10 (jak widać nie mogą mniej się od siebie
różnić). Obliczmy ich wartości:
1 1 1
(1)1101 (0)10 0.11012 2 102 4 3.25
2 4 16
1 1
(1)1100 (0)10 0.11002 2 102 4 3
2 4
Jak widać w żaden sposób nie uda nam się zapisać liczby  3.1 (!)  maszyna musi zapamiętać ją jako  3.0 (jest bliższe
prawidłowej wartości niż  3.25).
Ogólnie wprowadzane liczby mogą różnić się od siebie o nie więcej niż 2 liczba bitów mantysy 2cecha
Możemy więc przyjąć, że każda wprowadzona do komputera liczba musi być obarczona jakimś błędem  jego wartość
bezwzględna zależy od wielkości liczby, ale jego wartość odniesiona do wartości wpisywanej liczby jest zawsze podobna
 w naszym przykładzie jest to około 2 4 (1/16 czyli 0,0625). Inaczej można to zinterpretować jako ograniczenie liczby
cyfr znaczących wprowadzanych liczb. I tu przechodzimy do pierwszego rodzaju błędów:
Błędy wejściowe
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Błędy obliczeń numerycznych 2
Nazwa od razu kojarzy nam się z sytuacją gdy dane wejściowe są np. wynikami pomiarów fizycznych, które z natury
rzeczy obarczone są błędami. Ale chodzi tu przede wszystkim o błędy wynikające z wstępnego zaokrąglenia
wprowadzanych (nawet zupełnie dokładnych!) wartości (dotyczy to zwłaszcza liczb niewymiernych)  a wynika z powyżej
opisanego sposobu zapisu liczb.
Przykład: Rozważamy nadal 5-cio bitową mantysę (wraz z bitem znaku) i 3-y bitowy wykładnik (też wraz z bitem
znaku).
Chcemy zsumować 3 liczby (w zapisie dziesiętnym):
0.48
0.24
0.12
(wynik powinien oczywiście wynosić 0.84).
W zapisie dwójkowym mają one nieskończone rozwinięcia:
0.(01111010111000010100)
0.0(01111010111000010100)
0.00(01111010111000010100)
Ponieważ jednak nasz wymyślony komputer dysponuje tylko czterema bitami dla przedstawienia wartości liczby,
będzie musiał wprowadzane liczby zapamiętać jako:
(0)1111 (1)01
(0)1111 (1)10
(0)1111 (1)11
czyli tak naprawdę będziemy sumować: 0.46875, 0.234375 i 0.1171875 !
Gdybyśmy liczby te mogli dokładnie dodać w systemie dwójkowym to wynik powinien wynosić 0.8203125. A zatem
tylko z powodu występowania błędów wejściowych popełnimy błąd bezwzględny wynoszący około 0.02 ! (a
odpowiadający mu błąd względny wynosi 2,3%).
Błędy zaokrągleń
Kontynuujmy poprzedni przykład: w wyniku zsumowania liczb dwójkowych:
0.011112 + 0.0011112 + 0.00011112
powinniśmy otrzymać wynik:
0.11010012 (czyli 0.820312510)
Jednak ponieważ mamy do dyspozycji tylko 5 cyfr mantysy (a tak na prawdę tylko 4, bo jedna jest przeznaczona
na znak) będziemy sumować liczby:
0.011112 + 0.001112 + 0.000112
(to, że zapisane liczby mają po pięć cyfr po przecinku jest najzupełniej zgodne z założeniem o przeznaczeniu 4 cyfr
na mantysę. Liczby te różnią się wykładnikami. Faktem jest, że drugą i trzecią liczbę można by zapisać dokładniej
(z innymi wykładnikami), ale w celu ich dodania musimy najpierw uzgodnić wykładniki  w przeciwnym razie
procesor nie zdoła dodać wszystkich liczb.
W wyniku tej operacji otrzymalibyśmy
0.110012 (czyli 0.7812510 !)
ale liczba ta ma 5 cyfr (dwójkowych)  tyle nasz komputer nie jest w stanie zapamiętać. Wynik zostanie więc
zapisany jako:
(0)1100(0)00 czyli 0.112 (w zapisie dziesiętnym: 0.75).
Błąd bezwzględny wynosi więc  0.09 (błąd względny  10.7% !).
Kolejną ilustracją może być następujący przykład: chcemy dodać liczby:
1.87
0.93
0.47
(w zapisie maszynowym (przez zapis maszynowy będziemy rozumieć zapis dwójkowy z uwzględnieniem
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Błędy obliczeń numerycznych 3
ograniczenia liczby cyfr):
(0)1111 (0)01
(0)1111 (0)00
(0)1111 (1)01
co odpowiada faktycznie wartościom:
1.875
0.9375
0.46875
Mamy więc faktycznie dodawanie takich liczb dwójkowych:
1.1112
0.11112
0.011112
Po ich dodaniu powinniśmy dostać:
11.010012 czyli: 3.28125 (nawiasem mówiąc zamiast 3.27).
Tymczasem ze względu na obecność ograniczonej liczby cyfr znaczących wykładnika maszyna będzie musiała
przedstawić wszystkie te liczby z jednakowym wykładnikiem. W przypadku dodawania jakichkolwiek liczb
uzgodnienie wykładników zawsze będzie polegać na zrównaniu ich do najwyższego. Gdy zmieniamy wykładnik
liczby na wyższy (tzn. gdy mnożymy ją przez dodatnią potęgę podstawy  2) jej cyfry przesuwamy w prawo, a
przed nimi dopisujemy zera. Z powodu ograniczonej liczby cyfr ostatnie z nich  ginÄ… powodujÄ…c powstanie
pewnego błędu. Gdybyśmy chcieli uzgodnić wykładniki  w dół  do któregoś z niższych (tzn. pomnożyć liczbę
przez ujemną potęgę dwójki), musielibyśmy przesuwać cyfry w lewo, przy czym, znów z powodu ograniczenia
liczby cyfr musielibyśmy  gubić cyfry początkowe, a to spowodowałoby znacznie większe błędy. Jeśli więc
zaczniemy dodawanie od liczby największej to będzie to:
(0)1111 (0)01 + (0)1111 (0)00 + (0)1111 (1)01
po dopasowaniu wykładników:
(0)1111 (0)01 + (0)0111 (0)01 + (0)0011 (0)01 = (0)1100 (0)10
1.1112 + 0.1112 + 0.0112 = 11.002 (byłoby to 11.0012 ale ostatnia  1 została obcięta)
1.875 + 0.875 + 0.375 = 3 (na prawdę 3.125, ale  patrz wyżej ).
Błąd bezwzgledny wynosi 0.27 (a względny ok. 8%).
Jeśli zaś dodawanie zaczniemy od najmniejszego składnika dostaniemy:
(0)1111 (1)01 + (0)1111 (0)00
po dopasowaniu wykładników:
(0)0111 (0)00 + (0)1111 (0)00 = (0)1011 (0)01
0.01112 + 0.11112 = 1.0112 (faktycznie: 1.011012, ale obcięte)
0.4375 + 0.9375 = 1.375
i w drugim kroku (nie trzeba dopasowywać pierwszej z liczb bo wykładniki się zgadzają):
(0)1011 (0)01 + (0)1111 (0)01 = (0)1101 (0)10
1.0112 + 1.1112 = 11.012
1.375 + 1.875 = 3.25
dostajemy więc wynik wyraz
nie lepszy niż dodając w niewłaściwej kolejności  błąd bezwzględny wynosi 0.02,
a względny 0.6%.
Z powyższych rozważań wynika praktyczna reguła:
Podczas sumowania wielu liczb należy zadbać o to aby dodawać liczby w kolejności rosnącej
(najpierw małe potem duże).
Jeszcze większe błędy względne możemy uzyskać w przypadku odejmowania liczb o zbliżonych wartościach. Np.
chcemy wykonać działanie:
0.5  0.46875
0.12  0.011112 = 0.000012 (dokładnie: 1/32)
(0)1000 (0)00  (0)1111 (1)01 =
(0)1000 (0)00  (0)0111 (0)00 = (0)0001 (0)00 (czyli 1/16!)
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Błędy obliczeń numerycznych 4
w ten sposób otrzymaliśmy błąd względny wynoszący 100%!. Na takie sytuacje właściwie nie ma rady  trzeba
o tym pamiętać pisząc programy i tak zapisywać algorytmy, aby nie zachodziła konieczność odejmowania liczb o
zbliżonych wartościach.
Zróbmy jeszcze jeden  eksperyment : sprawdzimy czy 3.2  3.0 = 0.3  0.1 .
Pierwsze odejmowanie. Najpierw sprawdzimy jak obie liczby zostaną zapisane w pamięci naszego komputera
(podkreślona została wartość obarczona mniejszym błędem):
3.210 = 11.001100112 = (0)1100(0)10 = 3.0
= (0)1101(0)10 = 3.25
Liczba 3.2 zostanie zapisana jako (0)1101(0)10  czyli jako 3.25.
3.010 = 11 = (0)1100(0)10 = 3.0 (zapis dokładny)
(0)1101(0)10  (0)1100(0)10 = (0)0001(0)10 = 0.25
Wynikiem pierwszgo odejmowania będzie 0.25 (zamiast 0.2 !).
Drugie odejmowanie. Najpierw zapis obu liczb:
0.3 = 0.01001100110 = (0)1001(1)01 = 0.28125
= (0)1010(1)01 = 0.3125
0.1 = 0.00011001100 = (0)1100(1)11 = 0.09375
= (0)1101(1)11 = 0.1015625
Obie liczby zostaną zapisane z błędami.
(0)1010(1)01  (0)1101(1)11 = (0)1010(1)01  (0)0011(1)01 = (0)0111(1)01 = 0.21875
Jak widać wynik tym razem wynosi 0.21875  też z błędem, ale tym razem innym.
Z  eksperymentu tego wynika pewnien  morał : W przypadku liczb zmiennoprzecinkowych należy unikać sprawdzania
ich równości  jest bowiem bardzo prawdopodobne, że liczby o których wiemy na pewno, że są sobie równe (bo dwa różne
ciągi działań arytmetycznych przeprowadzone analitycznie (=na karteczce) dały identyczny wynik), w pamięci komputera
wcale równe być nie muszą.
W opisie Turbo Pascala można znalez
ć informację, że liczba bajtów przeznaczanych na zapis liczb rzeczywistych i ich
zakresy wynoszÄ… odpowiednio:
real 6 od 2.9×10 39 do 1.7×1038
single 4 od 1.5×10 45 do 3.4×1038
double 8 od 5.0×10 324 do 1.7×10308
extended 10 od 3.4×10 4932 do 1.1×104932)
Jednak oczywiście nie znaczy to, że dodanie najmniejszej możliwej liczby do największej cokolwiek zmieni! Liczby
z dopuszczalnego zakresu pamiętane są tylko z określoną liczbą cyfr znaczących (real  11, single  7, double  15
i extended  19).
Przykładem skutków błędów wynikających z
ograniczonego zakresu komputerowej reprezentacji
liczb rzeczywistych może być wykonanie następującej
operacji: do jednostki (danego typu) dodajemy małą
liczbÄ™ h (tego samego typu) po czym od wyniku
odejmujemy jednostkę. Oczekujemy, że w wyniku
otrzymamy właśnie małą liczbę h. Analitycznie jest to
twierdzenie bezwzględnie prawdziwe, jednak kiedy
przeprowadzimy tÄ™ operacjÄ™ na komputerze,
stwierdzimy, że gdy tylko h stanie się odpowiednio
małe zaczną występować błędy. Wykres obok przed-
stawia błędy względne uzyskiwane podczas działania
programu w języku Turbo Pascal wykonującego
właśnie taką operację. Błędy o wartości 1 odpowiadają
całkowicie błędnym wynikom (tzn. gdy dodanie h do
jedynki daje w wyniku jedynkÄ™). WykonywanÄ… przez
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Błędy obliczeń numerycznych 5
program procedurę można przedstawić tak:
h 1 h 1
´
h
Analogiczną  operację można przeprowadzić w Scilabie. Posłuży nam do tego następujący program:
// Ilustracja bł dów wynikaj cych z ograniczonej dokładno ci zapisu liczb.
//
clear
xdel(winsid());
//
h = logspace(-17,-13,100);
delta = (h - abs((h+1)-1) )./h;
xinit()
xset('thickness',2);
plot2d(h,zeros(h),logflag='ln',style=5);
xset('thickness',3);
plot2d(h,delta,logflag='ln',style=13)
xset('thickness',1);
a = get("current_axes");
// równie dobrze a = gca;
a.font_size = 4;
xtitle(sprintf('Wykres bł dów wynikaj cych z dokładno ci zapisu liczb\nwarto %%eps =
%g',%eps));
W wyniku jego działania otrzymamy następujący wykres:
Jak widać wprawdzie zmienna %eps Scilaba, która w zasadzie informuje nas o dokładności wynosi około 10 16 to jednak
bÅ‚Ä™dy wyrażenia ´ przekraczajÄ… 10% już o rzÄ…d wielkoÅ›ci wczeÅ›niej. OczywiÅ›cie wartość parametru h, przy której bÅ‚Ä…d
całego wyrażenia przestanie być akceptowalnie mały zależy od rozmiaru słowa maszynowego. Przy okazji warto zwrócić
uwagę, że %eps można zinterpretować jako najmniejszą liczbę jaką uda się dodać do jedności i zauważyć skutek tego
dodawania. Jeżeli chcielibyśmy dodawać ją do liczby większej (np. do 1000) to oczywiście taka najmniejsza  zauważalna
liczba będzie odpowiednio większa (o trzy rzędy wielkości).
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Błędy obliczeń numerycznych 6
Błędy obcięcia
Tak określay błędy występujące przy obliczaniu sum szeregów nieskończonych za pomocą skończonej liczby wyrazów
(wÅ‚aÅ›ciwie wszystkie funkcje sÄ… obliczane jako sumy szeregów, np. ex = 1 + x + ½x2 + ... + 1/N! xN). Wielkość takiego bÅ‚Ä™du
zależy od tego jak szybko zbieżny jest wykorzystywany szereg. Szczególnie łatwe jest oszacowanie błędu obcięcia w
przypadku szeregu, którego wyrazy są monotonicznie malejące i mają przemienne znaki  w takiej sytuacji błąd jest
dokładnie oszacowany z góry przez wartość pierwszego odrzuconego składnika.
Przykład: W przykładzie tym zrezygnujemy z przeprowadzania obliczeń w układzie dwójkowym i w reprezentacji
maszynowej. Założymy po prostu, że każda liczba jest zaokrąglana do trzech cyfr znaczących (nie do zadanego
miejsca po przecinku, a właśnie do zadanej liczby cyfr znaczących  tzn. 5.67 + 7.89 = 13.6 (chociaż dokładnie
powinno być 13.56)). Chcemy obliczyć ex za pomocą rozwinięcia ex = 1 + x + x2/2! + x3/3! + ... Rozważmy
obliczenia dla x = 0.5.
n wartości suma suma
wyrazów w dół w górę
dokł. zaokr. dokł. zaokr. dokł. zaokr.
0 1 1 1 1 1.64872 1.648
1 0.5 0.5 1.5 1.5 0.64872 0.648
2 0.125 0.125 1.625 1.62 0.14872 0.148
3 0.02083 0.0208 1.64583 1.64 0.02372 0.0236
4 0.0026 0.0026 1.64844 1.64 0.00289 0.00288
5 0.00026 0.00026 1.6487 1.64 0.00028 0.00028
6 0.00002 0.00002 1.64872 1.64 0.00002 0.00002
7 1.6e-06 1.6e-06 1.64872 1.64 1.7e-06 1.7e-06
8 9.7e-08 9.7e-08 1.64872 1.64 1.0e-07 1.0e-07
9 5.4e-09 5.4e-09 1.64872 1.64 5.7e-09 5.7e-09
10 2.7e-10 2.7e-10 1.64872 1.64 2.8e-10 2.8e-10
11 1.2e-11 1.2e-11 1.64872 1.64 1.3e-11 1.3e-11
12 5.1e-13 5.1e-13 1.64872 1.64 5.3e-13 5.3e-13
13 2.0e-14 2.0e-14 1.64872 1.64 2.0e-14 2.0e-14
14 7.0e-16 7.0e-16 1.64872 1.64 7.2e-16 7.2e-16
15 2.3e-17 2.3e-17 1.64872 1.64 2.3e-17 2.3e-17
sumy 1.64872 1.64 1.64872 1.648
błędy bezwzględne/względne 0.00872 0.529% 0.00072 0.044%
Inne z
ródła błędów
Omówiłem już przyczyny występowania błędów wynikające z samej istoty działania komputera. Nic nie wspomniałem
jednak o  czynniku ludzkim (który, nawiasem mówiąc, jest najczęstszą ich przyczyną).
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Błędy obliczeń numerycznych 7
Załóżmy na początek, że mamy program obliczający wartość powierzchni koła o promieniu r (mam nadzieję, że wszyscy
pamiętają, że wynosi ona Ąr2). Błąd jaki program będzie konsekwentnie popełniać będzie wynikać ze skończoności zapisu
rozwinięcia dwójkowego wartości Ą i wartości promienia r oraz z błędów obcięcia popełnianych podczas podnoszenia r
do kwadratu i mnożenia obu wartości. Nic na to nie możemy poradzić (poza oczywiście użyciem maszyny bądz
kompilatora
o dłuższym słowie), ale błędy te będą w rzeczywistości bardzo małe (pamiętajmy, że w przedstawionych wcześniej
przykładach używaliśmy słowa ośmiobitowego  tak prymitywnych komputerów w rzeczywistości nie ma). Co jednak się
stanie jeśli programista za wartość Ą podstawi 3.141 ? Spodziewany błąd będzie wynosić około 0.0005 (!)  będzie więc
wielokrotnie większy niż wszystkie błędy jakie zdoła popełnić komputer... Przykład ten wskazuje dość dobitnie na fakt,
że programista, używając niedokładnych stałych, jest w stanie spowodować błędy o wiele rzędów wielkości większe niż
byłby to w stanie zrobić sam komputer. Jeżeli zależy nam na dokładności i za pomocą wyszukanych sztuczek zdołaliśmy
zmusić kompilator do działania na liczbach o (np.) dwukrotnie większej precyzji niż robi to zazwyczaj  sam z siebie , to
powinniśmy się również upewnić, że zdefiniowane w nim wartości stałych w rodzaju Ą albo e również mają wystarczającą
dokładność. Podobnie wygląda sytuacja z funkcjami zdefiniowanymi w kompilatorze. W starszych wersjach Turbo Pascala
(np.) funkcje trygonometryczne mogły zwracać wartości o podwyższonej dokładności, ale liczyły je zawsze z dokładnością
 normalną (a cyfry znaczące brakujące do podwyższonej dokładności uzupełniały  śmieciami , które akurat znajdowały
się w tym miejscu w pamięci). Można oczywiście dyskutować czy błędy spowodowane w ten sposób są zawinione przez
programistę (który wykazał się nadmiernym zaufaniem do kompilatora) czy przez producenta kompilatora, ale pamiętajmy,
że za błędy w wynikach programu przed zleceniodawcą odpowiadać będzie autor programu, a  opowiadanie o błędach
kompilatora i tak będzie wyglądać na nieudolną próbę  zwalenia winy na kogoś innego.
Podczas omawiania wszystkich dotychczasowych błędów zakładaliśmy, że komputer faktycznie liczy to co powinien.
Tymczasem z
ródła błędów (i to tych największych) pojawiają się już na etapie analizy rozwiązywanego zagadnienia.
Posłużę się tu  nieśmiertelnym przykładem rzutu ukośnego: należy rozwiązać równania ruchu ciała rzuconego ukośnie.
Już na etapie formułowania równań różniczkowych ruchu (czyli nawet jeszcze przed wyborem algorytmu rozwiązywania)
musimy podjąć następujące decyzje:
a) czy uwzględniać opór powietrza czy nie (oczywiście tak  w przeciwnym razie wyniki będą stosowalne co najwyżej
na Księżycu bądz
w komorze próżniowej),
b) w jaki sposób opór powietrza zależy od prędkości: zależność liniowa  zgodna z wzorami Stokesa odnosi się do
cieczy, zależność kwadratowa lepiej pasuje do powietrza (dla niej przynajmniej mamy wyprowadzone wzory, jednak
doświadczenia pokazują, że w grę mogą wchodzić też wyższe potęgi prędkości),
c) w tym miejscu widać, że musimy zdecydować czy dopuszczamy występowanie wiatru czy nie, a jeśli tak to czy będzie
on zawsze wiać dokładnie poziomo czy też jego prędkość może mieć składową pionową,
d) jaki współczynnik proporcjonalności zastosować we wzorze na opór powietrza (zakładając, że wzór już wybraliśmy)
 jego wartość zależy od kształtu ciała (trzebaby więc przyjąć jakąś wartość doświadczalną) na dodatek wiadomo, że sama
wartość tego współczynnika jest funkcją prędkości ciała (względem powietrza  więc z uwzględnieniem wiatru) ,
e) we wzorze na siłę oporu (jakiego byśmy nie przyjęli) występuje gęstość powietrza trzeba więc zdecydować czy
uwzględniamy zmiany siły oporu z wysokością na jaką wzniesie się pocisk czy nie, a na dodatek należałoby się zastanowić
czy będziemy uwzględniać aktualne ciśnienie atmosferyczne i wilgotność powietrza,
f) powinniśmy się zastanowić czy nasz pocisk będzie wykonywać jakieś ruchy obrotowe czy nie  jeśli tak to
należałoby uwzględnić efekt Magnusa, który doda nam nową składową siły działającej na pocisk skierowaną prostopadle
do wektora prędkości pocisku (względem powietrza) i do wektora jego prędkości obrotowej,
g) warto byłoby też zastanowić się jaki wpływ na trajektorię ciała będzie mieć siła Coriolisa (w końcu Ziemia jest
układem nieinercjalnym).
Tego rodzaju rozważania można ciągnąć dalej  pocisk może rozgrzewać się w wyniku tarcia o powietrze (zależnie
od swej prędkości względem powietrza) zatem zmieniać się może jego pole przekroju, zmieni się też prawdopodobnie
lepkość powietrza (bo ono też się rozgrzeje tuż przy ścianie pocisku) itd. itp.. Jak więc widać, już analizując samo zadanie
musimy przyjąć wiele założeń upraszczających, a to powoduje, że właściwie nie można mówić o rozwiązaniach dokładnych
(w jakimś abstrakcyjnym, bezwzględnym sensie). Każde rozwiązanie będzie mieć swój szczególny zakres w jakim można
będzie z grubsza podać oszacowanie jego dokładności. W powyższym przykładzie jeśli np. zaniedbamy zależność gęstości
powietrza od ciśnienia, powinniśmy uprzedzić użytkownika do jakich wysokości dopuszczamy stosowanie naszego
przybliżonego rozwiązania (a sam program powinien  lojalnie informować, że wprowadzono takie parametry wejściowe,
przy których pocisk wznosi się na taką wysokość, na jakiej zastosowane w programie przybliżenia spowodują zbyt wielkie
błędy).
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Błędy obliczeń numerycznych 8
No i jeszcze ostatnia przyczyna  milcząco założyliśmy, że autor programu w ogóle umiał rozwiązać postawiony przed
nim problem i miał świadomość stosowanych przybliżeń. W krańcowym przypadku (może nieco humorystycznym, ale
chyba dobrze ilustrującym problem), ktoś kto znajomość z fizyką zakończył w szkole podstawowej, a miałby akurat
rozwiązać omawiany powyżej problem, przygotowałby program doskonale liczący przypadek bez oporu powietrza i bez
siły Coriolisa i w ogóle nie powiadomiłby użytkownika o istniejących ograniczeniach stosowalności swego programu (bo
po prostu nie byłby ich świadomy).
Podczas zajęć z Metod Numerycznych będziemy zakładać, że naszym jedynym problemem są  psikusy płatane nam
przez metodę (algorytm) rozwiązania. O  psikusach płatanych nam przez komputer, a wynikających z jego cyfrowej
natury, będziemy, pisząc programy, zawsze pamiętać (ale nie będę do nich już wracać w dalszym ciągu wykładu). Natomiast
nie będę przypominać o kwestiach, o których wspomniałem przed chwilą  innymi słowy zakładam, że programista
dysponuje dostateczną wiedzą fizyczną aby problem rozwiązać poprawnie. Cały czas musimy jednak pamiętać (może nawet
zwłaszcza po ukończeniu studiów), że nawet najlepszy algorytm numeryczny nie  uratuje programu, którego założenia
są błędne.
Niestabilność metod obliczeniowych (powodowana przez błędy obcięcia)
Wez
my jako przykład następujące zadanie: Należy obliczyć wartości całek postaci:
1
n
x
yn dx
x 5
0
dla n = 0,1,...,8. Zauważymy, że można zastosować wzór rekurencyjny:
1
yn 5yn 1
n
co wynika z:
1 1 1
n n 1 n 1
x 5x x (x 5) 1
n 1
yn 5yn 1 dx dx x dx
x 5 x 5 n
0 0 0
W przykładzie obliczenia będziemy wykonywać z trzema cyframi po przecinku, a błąd wynikający z obcięcia oznaczymy
przez µ. WykonujÄ…c obliczenia kolejno od y0 bÄ™dziemy otrzymywać:
1
dx 1
yo ln(x 5) ln6 ln5
0
x 5
0
y0 = = 0.182 (bÅ‚Ä…d µ)
y1 = 1  5y0 = 0.090 (bÅ‚Ä…d 5µ)
y2 = ½  5y1 = 0.050 (bÅ‚Ä…d 25µ)
y3 =  5y2 = 0.083 (bÅ‚Ä…d 125µ)
y4 = ź  5y3 =  0.165 (sic!) (bÅ‚Ä…d 625µ)
W tym miejscu zatrzymujemy się, bo przecież z samej postaci obliczanej całki wynika, że dla każdego n powinniśmy
otrzymać wartość dodatnią !
Błąd wynika tu z faktu, że przy każdym kolejnym działaniu nie tylko dodajemy nowy (bardzo mały) błąd związany z
wykonanym mnożeniem i dodawaniem, ale równocześnie, ponieważ poprzedni wynik mnożymy przez 5, pięciokrotnie
powiększamy błąd jakim obarczona była poprzednio obliczona liczba  dla każdej kolejnej błąd początkowy będzie 5 razy
większy (wypisane w kolumnie obok).
Aby tego uniknąć należy wykorzystać wzór rekurencyjny w przeciwnym kierunku:
1 1
yn 1 yn
5n 5
 w tym przypadku błąd popełniany w poprzednich krokach będzie zawsze dzielony przez 5  będzie więc maleć.
Przy takiej metodzie obliczeń problemem będzie przyjęcie wartości początkowej (którą musi być jakieś yn  dla
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Błędy obliczeń numerycznych 9
dostatecznie dużej wartości n). Istnieją dwie możliwości: założenie, że poczynając od jakiejś wartości n wartości szeregu
yn maleją tak wolno, że dwie kolejne będą równe sobie, albo że dla wystarczająco dużego n wartośc yn wynosi zero.
1. Załóżmy, że y10 y9, wówczas:
y9 + 5y9 1/10 y9 1/60 0.017
y8 = 1/45  y9/5 0.019
y7 = 1/40  y8/5 0.021
y6 0.025
y5 0.028
y4 0.034
y3 0.043
y2 0.058
y1 0.088
y0 0.182 (to samo co z obliczenia dokładnego)
2. Załóżmy, że y10 = 0, wówczas dostajemy:
y9 = 0.020
y8 = 0.018
a pozostałe wartości jak wyżej.
Warto zapamiętać, żechoć oba przedstawione wzory rekurencyjne są z matematycznego punktu widzenia poprawne,
to jednak jeden z nich z samej swej natury nie nadaje się do obliczeń numerycznych. Z takimi przypadkami będziemy się
jeszcze spotykać. Istnieją metody pozwalające przewidzieć czy dany algorytm obliczeniowy będzie powodować wzrost
błędów (czyli będzie  niestabilny ) czy też wręcz przeciwnie  będzie ograniczać błędy. Zapoznamy się z nimi póz
niej.
Ćwiczenie: Napisać wzór rekurencyjny dla obliczania całek:
1
n
x
yn dx
4x 1
0
Podać i przeanalizować algorytm działający dobrze i algorytm działający z
le.
Prędkość wykonywania operacji arytmetycznych
Interesujące może być porównanie prędkości wykonywania poszczególnych operacji arytmetycznych. W przypadku
Scilaba warto przy okazji porównać prędkość wykonywania operacji na macierzach z prędkością wykonywania pętli.
Tabela poniżej podaje czasy wymagane do wykonania miliona operacji na przykładowym (moim;) komputerze.
Kolumna  pętla odnosi się do działania na pojedynczym skalarze (powtórzonego oczywiście milion razy), dane z kolumny
 macierz dotyczą jednej operacji za to wykonanej na wektorze zawierającym milion elementów.
operacja pętla macierz krotność
dodawanie 7.859375 0.078125 100.6
odejmowanie 7.921875 0.09375 84.5
mnożenie 7.390625 0.078125 94.6
dzielenie prawostronne 7.28125 0.09375 77.6666667
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Błędy obliczeń numerycznych 10
dzielenie lewostronne 7.328125 0.09375 78.1666667
podnoszenie do dowolnej potęgi 7.859375 0.40625 19.3461538
podnoszenie do kwadratu 7.34375 0.171875 42.7272727
Widać dość wyraz
nie, że warto poznać i stosować operacje macierzowe. W tabeli nie przedstawiłem wyników dla
przypadku gdy w pętli oblicza się kolejne wartości wektora równocześnie go powiększając (ostrzegałem przed taką metodą
na zajęciach)  czasy są ok. dziesięciokrotnie dłuższe niż czasy realizacji pętli (więc też 10000 razy dłuższe od prawidłowo
realizowanych operacji macierzowych  jak widać ostrzeżenie nie było bezpodstawne).
Wyniki przedstawione w powyższej tabeli uzyskano za pomocą prostego programu, którego fragment znajduje się
poniżej (aby samodzielnie przeprowadzić takie testy wystarczy powielić go dla każdego działania, którego prędkość
wykonywania chcemy sprawdzić).
clear
xdel(winsid())
stacksize(10000000)
// Tu nale ałoby pogada o poleceniu "stacksize".
//
rand('seed',.123456);
n = 10000 // liczba powtórze
//
// dodawanie
a = rand(1,1);
b = rand(1,1);
timer();
for i = 1:n
c = a + b;
end
dodawanie_petla = timer()
clear c
aa = ones(1,n)*a;
bb = ones(1,n)*b;
timer();
cc = aa + bb;
dodawanie_macierz = timer()
clear cc
timer();
for i = 1:n
ccc(i) = a + b;
end
dodawanie_petla_macierz = timer()
clear ccc
Schemat Hornera
Inny sposób usprawnienia obliczeń  wzory algebraiczne można zapisywać w różnych postaciach, które mogą lepiej
lub gorzej nadawać się do obliczania. Dobrym przykładem jest schemat Hornera  sposób obliczania wartości wielomianu.
Wielomian postaci:
w(x) = anxn + an 1xn 1 + an 2xn 2 + ... + a1x + ao
należy zapisać następująco:
w(x) = ((...((anx + an 1)x + an 2)x + ... + a1)x + ao)
i obliczać w taki sposób:
w1 = anx + an 1
w2 = xw1 + an 2
.............
wn = xw1 + ao = w(x)
Uzyskujemy w ten sposób znaczną oszczędność czasu obliczeń: potrzeba będzie n mnożeń i n dodawań (zakładając,
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Błędy obliczeń numerycznych 11
że wszystkie współczynniki ai są niezerowe)  podczas gdy obliczając każdą potęgę z osobna należałoby wykonać n2/2
mnożeń i n dodawań.
Pisanie programów
Styl programowania
Komentarze
Istnieje wzór empiryczny:
t tzN
s
gdzie: ts jest czasem straconym z powodu braku komentarzy podczas próby poprawienia lub zmienienia programu, tz jest
czasem zaoszczędzonym dzięki ich nie pisaniu, a N jest liczbą tygodni jakie upłynęły od napisania programu.
Jest to oczywiście żart, ale zawiera sporą dozę prawdy  pod koniec semestru, podczas pisania programów
rozwiązujących równania różniczkowe cząstkowe metodą różnic skończonych konieczne będzie posłużenie się procedurami
rozwiązującymi układy równań algebraicznych (które będą opracowywane na początku semestru). Będzie to najlepszy test
prawdziwości powyższego wzoru.
Powtórne wykorzystywanie procedur i funkcji
Jeśli pisane są jakieś procedury, które mają znaczenie ogólne i można by je wykorzystać w innych programach
(procedury do rozwiązywania układów równań algebraicznych będą pod koniec semestru potrzebne do rozwiązywania
równań różniczkowych, ale dotyczy to również ew. procedur komunikowania się z użytkownikiem itp.) to trzeba je pisać
szczególnie starannie  zaznaczając w komentarzach jakie mają wymagania: tzn. np. jakich stałych globalnych spodziewają
siÄ™ w programie (w zasadzie nie powinny), w jakiej postaci przyjmujÄ… dane i w jakiej postaci zwracajÄ… wyniki, jakie sÄ… ich
ograniczenia (np. maksymalny rozmiar czy ilość danych) itd.
Sposób nazywania zmiennych i stałych
Stosowanie nazw w rodzaju:  c ,  p17 ,  k3 znakomicie utrudnia póz
niejsze zrozumienie programu. Obecnei, z reguły
nazwy mogą mieć niemal dowolną długość, ale do ich rozróżniania wykorzystywanych jest tylko pewna liczba
początkowych znaków. W Turbo Pascalu są to pierwsze 63 znaki, w Scilabie pierwsze 24 znaki  jesto to wystarczająca
liczba by można było posługiwać się nazwami w rodzaju:  wyznacznik ,  pierwiastek_1 (ew.  pierw_1 ),  x_start ,
 x_koniec itd. co z kolei bardzo ułatwia póz
niejszą interpretację zapisanych wyrażeń.
Usuwanie błędów
Właściwie nie istnieje coś takiego jak program napisany od razu bezbłędnie. Popełnianie błędów podczas pisania
programu jest rzeczą normalną i nieuniknioną (możemy oczywiście popełniać ich więcej lub mniej zależnie od tego jak
poważnie podejdziemy do unikania błędów). Istotnym problemem jest jednak przede wszystkim znajdowanie i usuwanie
błędów. W zależności od tego jak napisaliśmy program może to być łatwiejsze lub trudniejsze. Jeśli rozwiązywany problem
wymaga długotrwałych obliczeń warto na początek zacząć pisać program dla mniejszego zakresu, o mniejszej dokładności
itp. W ten sposób można szybciej stwierdzić występowanie błędów. W ostatecznej wersji programu, dla przyspieszenia jego
działania często rezygnujemy z jakichkolwiek komunikatów ekranowych, w wersji początkowej nie powinniśmy z nich
rezygnować, wręcz przeciwnie: należy generować znaczną ich liczbę  w ten sposób będzie wiadomo co program akurat
robi i podczas wykonywania jakiej procedury  padł . Oczywiście wygodnymi narzędziami do poszukiwania błędów są
debuggery wbudowane np. w Borland Pascal czy Borland C++ (ale trzeba się najpierw nauczyć nimi sprawnie posługiwać),
natomiast komunikaty wypisywane przez program podczas pracy, umożliwiając szybkie zidentyfikowanie błędnie
działającej procedury często pozwalają bardzo szybko zidentyfikować i usunąć błąd.
Zawsze warto pamiętać o tym, że błędy w programie są w 99% zawinione przez programistę, a w pozostałym 1%
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej
Metody Numeryczne w Fizyce 1 Błędy obliczeń numerycznych 12
również, tylko, że on o tym jeszcze nie wie. Wiadomo oczywiście, że istnieją błędy kompilatorów, systemów operacyjnych,
a nawet procesorów, jednak prawdopodobieństwo, że wystąpił któryś z tych błędów  niezależnych od programisty jest
znikome. Mówiąc o błędach popełnionych przez programistę miałem na myśli nie tylko trywialną zamianę plusa na minus
albo mnożenia na dzielenie (które to błędy niestety zazwyczaj bardzo trudno znalez
ć, jako że program działa bardzo dobrze,
tylko że daje bezsensowne wyniki) ale także (a może przede wszystkim) błędy polegające na niewłaściwym doborze metody
rozwiązania czy niewłaściwej implementacji algorytmu, a nawet wyborze niewłaściwej metody podejmowania decyzji o
zatrzymaniu programu.
Sporządzanie wykresów
Ten podpunkt pozostał  z dawnych czasów  gdy studenci pisali swoje programy na MNF1 jeszcze w DOS-owych
kompilatorach C lub Turbo Pascal, dających możliwośc wykorztystywania tylko najbardziej podstawowych procedur
graficznych  w rodzaju wyświetlenia pojedynczego punktu bądz
narysowania odcinka pomiędzy punktami o zadanych
współrzędnych. Gdy wykorzystujemy Scilaba wiadomości te nie będą potrzebne, bo można (i należy) posłużyć się
odpowiednimi procedurami Scilaba.
© Andrzej Brozi, Instytut Fizyki Politechniki A
ódzkiej