Metody numeryczne dla fizyków
Krzysztof Golec–Biernat
(1 maja 2015)
Wersja robocza nie do dystrybucji
Kraków
2007-15
Spis treści
1 Liczby i ich reprezentacje 7
1.1 Systemy liczbowe . . . 7
1.2 System dwójkowy . . . 7
1.3 Konwersja do systemu dwójkowego . . . 9
1.3.1 Konwersja liczby całkowitej . . . 9
1.3.2 Konwersja liczby ułamkowej . . . 10
1.4 Reprezentacja liczb całkowitych . . . 10
1.5 Reprezentacja zmiennoprzecinkowa . . . 12
1.5.1 Podział bitów . . . 12
2 Błe¸dy i ich źródła 14 2.1 Dokładność zapisu . . . 14
2.2 Błąd bezwględny i względny . . . 15
2.2.1 Przyklad . . . 16
2.3 Źródła bł¸edów . . . 16
2.3.1 Błędy wejściowe. . . 16
2.3.2 Bł¸edy obci¸ecia. . . 17
2.3.3 Bł¸edy zaokr¸agleń. . . 18
2.3.4 Przyklad . . . 18
2.4 Odejmowanie małych liczb . . . 19
3 Algorytmy 21 3.1 Algorytmy stabilne . . . 21
3.2 Algorytmy niestabilne . . . 22
3.3 Złożoność obliczeniowa . . . 24
4 Macierze 25
4.1 Układ równań liniowych . . . 25
4.2 Metoda eliminacji Gaussa . . . 26
4.3 Rozkład LU macierzy . . . 28
4.4 Metoda Doolittle’a . . . 30
4.5 Wyznacznik macierzy . . . 31
5 Interpolacja 33 5.1 Wielomian interpolacyjny Lagrange’a . . . 33
5.2 Interpolacja Lagrange’a znanej funkcji . . . 35
5.3 Interpolacja przy pomocy funkcji sklejanych . . . 36
6 Interpolacja optymalna 38 6.1 Wielomiany Czebyszewa . . . 38
6.2 Optymalny wybór w¸ezłów interpolacji . . . 39
6.3 Wzory dla dowolnego przedziału . . . 40
7 Aproksymacja 42 7.1 Regresja liniowa . . . 42
7.2 Aproksymacja średniokwadratowa . . . 43
8 Przykłady aproksymacji 46 8.1 Aproksymacja Czebyszewa . . . 46
8.2 Aproksymacja Czebyszewa w dowolnym przedziale . . . 48
8.3 Aproksymacja trygonometryczna . . . 49
8.4 Wzory dla dowolnego okresu . . . 51
9 Różniczkowanie 52 9.1 Metoda z aproksymacj¸a . . . 52
9.2 Metody z rozwini¸eciem Taylora . . . 52
9.2.1 Wi¸eksza dokładność . . . 53
9.3 Wyższe pochodne . . . 54
10 Zera funkcji 55
10.1 Metoda połowienia . . . 55
10.2 Metoda Newtona . . . 57
10.3 Metoda siecznych . . . 59
10.4 Metoda fałszywej prostej (falsi) . . . 60
11 Całkowanie 62 11.1 Kwadratury . . . 62
11.2 Metoda prostokatów . . . 63
11.3 Metoda trapezów . . . 63
11.4 Metoda parabol Simpsona . . . 64
11.5 Bł¸ad przybliżeń . . . 65
12 Kwadratury Gaussa 66 12.1 Rz¸ad kwadratury . . . 67
12.2 Wielomiany ortogonalne . . . 68
12.3 Kwadratura Gaussa . . . 69
12.4 Współczynniki kwadratury Gaussa . . . 70
12.5 Przykład kwadratury Gaussa . . . 71
13 Równania różniczkowe 72 13.1 Równania zwyczajne pierwszego rz¸edu . . . 72
13.2 Układ równań różniczkowych . . . 73
13.3 Równania wyższych rz¸edów . . . 73
13.4 Metody Eulera . . . 74
13.5 Metoda przewidź i popraw . . . 75
13.6 Stabilność numeryczna rozwi¸azań . . . 76
13.7 Równania typu stiff . . . 77
14 Metoda Rungego-Kutty 78 14.1 Metoda drugiego rz¸edu . . . 78
14.2 Metoda czwartego rz¸edu . . . 80
15 Metody Monte Carlo 81 15.1 Zmienna losowa i rozkład prawdopodobieństwa . . . 81
15.2 Przykłady rozkładów prawdopodobieństwa . . . 83
15.3 Generowanie liczb losowych . . . 84
15.4 Całkowanie metod¸a Monte Carlo . . . 85
Wykład 1
Liczby i ich reprezentacje
1.1 Systemy liczbowe
Przypomnijmy zapis liczby całkowitej w systemie dziesie
‘tnym, na przykład (1234)10 = 4 · 100+ 3 · 101+ 2 · 103+ 1 · 104.
Do zapisu używamy dziesie
‘ć cyfr {0, 1, 2, . . . , 9}, a 10 to podstawa systemu liczbowego. Wartość cyfry zależy od pozycji w liczbie.
W ogólności, dowolna
‘ liczbe
‘ całkowita
‘ można zapisać w systemie pozycyj- nym o podstawie p przy pomocy cyfr ci∈ {0, 1, . . . , p − 1} w następujacy sposób (ck. . . c1c0)p = c0· p0+ c1· p + . . . + ck· pk (1.1) Podobnie, dowolna
‘ liczbe
‘ ułamkowa
‘ z przedziału (0, 1) możemy zapisać przy pomocy (na ogół nieskończonego) rozwinie
‘cia w ujemnych potęgach podstawy p,
(0. c−1c−2. . . c−n)p = c−1· p−1+ c−2· p−2+ . . . + c−n· p−n (1.2)
1.2 System dwójkowy
Szczególnie ważna
‘role
‘ odgrywa system dwójkowy o podstawie p = 2 i cyfrach ci∈ {0, 1}. Znaczenie tego systemu dla komputerów wynika z faktu, że pod- stawowa
‘ jednostka
‘ informacji jest bit przyjmuja
‘cy jedna
‘ z dwóch możliwych wartości: 0 lub 1. Przyjęło się nazywać
1 bajt = 8 bitów
Przykładowy zapis liczby całkowitej w systemie dwójkowym to (11010)2= 0 · 20+ 1 · 21+ 0 · 22+ 1 · 23+ 1 · 24 = 2 + 8 + 16 = (26)10.
Podobnie dla liczby ułamkowej
(0.1101)2 = 1 ·12+ 1 ·14+ 0 ·18+ 1 ·161 =1316
10.
Załóżmy, że dysponujemy n bitami do zapisu liczby całkowitej. Pozwala- ją one zapisać 2n liczb. Najmniejsza
‘ liczba
‘ całkowita
‘ dodatnia
‘ zapisana
‘ przy
pomocy n bitów jest zero
(00 . . . 00)2= 0 , (1.3)
natomiast najwie
‘ksza
‘liczba
‘ jest
(11 . . . 11)2 = 1 · 20+ 1 · 21+ . . . + 1 · 2n−1 = 2n− 1 . (1.4) Wykorzystaliśmy tutaj wzór na sume
‘ skończonej liczby wyrazów szeregu geo- metrycznego o ilorazie a 6= 1:
1 + a + a2+ . . . + aN −1=1 − aN
1 − a (1.5)
Podsumowuja
‘c, dysponuja
‘c n bitami możemy zapisać 2n liczb całkowitych z przedziału
[0, 2n− 1] . (1.6)
Podobnie, przy pomocy n bitów możemy zapisać 2n− 1 liczb ułamkowych z przedziału (0, 1). Najmniejsza
‘ z nich jest
(0.00 . . . 01)2 = 1 · 2−n = 1
2n (1.7)
natomiast najwie
‘ksza to
(0.11 . . . 11)2 = 1 · 2−1+ 1 · 2−2+ . . . + 1 · 2−n= 1 − 1
2n. (1.8) Podsumowując, przy pomocy n bitów zapiszemy (2n− 1) (brak zera) liczb z przedziału
1
2n, 1 − 1 2n
, (1.9)
Wygodnie jest zapamie
‘tać naste
‘puja
‘ca tabelke
‘ ułatwiaja
‘ca
‘ zmiane
‘ zapisu liczb w systemie dwójkowym na system dziesie
‘tny.
n 2n 2−n
0 1 1
1 2 1/2
2 4 1/4
3 8 1/8
4 16 1/16
5 32 1/32
6 64 1/64
7 128 1/128
8 256 1/256
9 512 1/512
10 1024 1/1024 11 2048 1/2048 12 5096 1/5096
1.3 Konwersja do systemu dwójkowego
1.3.1 Konwersja liczby całkowitej Konwersje
‘ odwrotna
‘z systemu dziesiętnego na dwójkowy dla liczb całkowitych realizuje sie
‘ zauważaja
‘c, że zachodzi
x = c0+ c1· 21+ . . . + ck· 2k = c0+ 2 · (c1+ . . . + ck· 2k−1) (1.10) Tak wie
‘c c0 to reszta z dzielenia x przez 2. Podobnie c1 jest reszta
‘ z dzielenia (c1+ . . . + ck· 2k−1) przez 2, itd.
Przykładowo, konwertuja
‘c liczbę 30 otrzymamy 30 = 2 · 15 + 0 15 = 2 · 7 + 1
7 = 2 · 3 + 1 3 = 2 · 1 + 1 1 = 2 · 0 + 1 Odczytuja
‘c reszty od dołu do góry otrzymujemy:
30 = (11110)2. (1.11)
1.3.2 Konwersja liczby ułamkowej Dla liczby ułamkowej,
x = c−1· 2−1+ c−2· 2−2+ . . . + c−n· 2−n (1.12)
po pomnożeniu przez dwa, otrzymujemy
2x = c−1+ (c−2· 2−1+ . . . + c−n· 2−n+1)
| {z }
r
(1.13)
Jak łatwo zauważyć c−1 jest cze
‘ścia
‘ całkowita
‘ liczby 2x (równa
‘ 0 lub 1) nato- miast r < 1 jest cze
‘ścia
‘ ułamkowa
‘. Aby wydobyć c−2 mnożymy wtedy r przez dwa, itd.
Przykładowo, dla 0.375 otrzymujemy
0.375 × 2 = 0. 750 0.750 × 2 = 1. 500 0.500 × 2 = 1. 000 Odczytuja
‘c cze
‘ści całkowite liczb po prawej stronie równości od góry do dołu, dostajemy
0.375 = (0.011)2. (1.14)
Otrzymana reprezentacja dwójkowa jest skończona. Rozwinie
‘cie dwójkowe ułam- ków może być jednak nieskończone, podobnie jak w systemie dziesie
‘tnym.
Ła‘cza
‘c dwie przedstawione metody konwersji, możemy znaleźć reprezentacje
‘
dwójkowa
‘ dowolnej liczby rzeczywistej, na przykład:
30.375 = 30 + 0.375 = (11110.011)2 (1.15)
1.4 Reprezentacja liczb całkowitych
Załóżmy, że przeznaczamy n bitów do zapisu liczby całkowitej w komputerze.
Pierwszy z nich przeznaczamy na znak liczby:
(0) = (+) (1) = (−) (1.16)
Pozostałe (n − 1) bitów posłuży do zapisu liczb całkowitych C z przedziału C ∈ [ −2n−1+ 1 , 2n−1− 1 ] . (1.17) Zauważmy, że poste
‘puja
‘c w ten sposób zero jest reprezentowane dwukrotnie ze wzgle
‘du na znak:
0 = (0) 0 . . . 0
| {z }
n−1
oraz 0 = (1) 0 . . . 0
| {z }
n−1
(1.18)
Przyjmuja
‘c tylko pierwszy sposób zapisu unikamy tej niejednozanczności zwal- niaja
‘c jednocześnie miejsce dla dodatkowej liczby, na przykład 2n−1. Tak wie
‘c, przy pomocy n bitów można zapisać 2nliczb całkowitych z przedziału
C ∈ [ −2n−1+ 1 , 2n−1] (1.19)
Podkreślmy, że liczby całkowite są reprezentowane dokładnie.
W praktyce, jednoznaczną reprezentację zera można zrealizować w następu- jący sposób. Przeznaczmy n bitów na zapis 2n nieujemnych liczb całkowitych C0 z przedziału
C0∈ [ 0 , 2n− 1 ] . (1.20) Przesuńmy następnie wszystkie liczby całkowite odejmując od każdej z nich liczbę E = 2n−1− 1. Otrzymujemy liczby
C = C0− E (1.21)
z przedziału (1.19). Zauważmy, że zero nie jest już reprezentowane przez pierw- szy ciąg bitów w reprezentacji (1.18). Jego reprezentacja odpowiada teraz dwój- kowej reprezentacji liczby E.
Współczesne procesory przeznaczają w pojedynczej precyzji 4 bajty czyli 32 bity do zapisu liczby całkowitej lub 8 bajtów czyli 64 bity w podwójnej precyzji. Poniższa tabelka pokazuje zakres liczb całkowitych w tych przypad- kach
Precyzja Liczba bitów Zakres Uwagi
Pojedyncza 32 [ −231+ 1 , 231] 231∼ 109
Podwójna 64 [ −263+ 1 , 263] 263∼ 1019
1.5 Reprezentacja zmiennoprzecinkowa
Przedyskutujemy teraz zapis liczb rzeczywistych w komputerze. Najbardziej efektywny zapis wykorzystuje reprezentacje
‘ zmiennoprzecinkowa
‘. W repre- zentacji tej nie ustala sie
‘ maksymalnej liczby cyfr po przecinku przy zapisie liczby.
Przykładowo, rozważmy liczbe
‘ x = 0.00045. Możemy ja zapisać w naste
‘pu-
ja‘cy sposób
x = 45 · 10−5 = 4.5 · 10−4 = 0.45 · 10−3. Liczba stoja
‘ca przez pote
‘ga
‘ dziesia
‘tki to mantysa M , natomiast wykładnik pote‘gi to cecha C. Dowolna
‘ liczbe
‘ rzeczywista
‘można wie
‘c zapisać w postaci
x = ± M · 10C. (1.22)
Przyjmujemy konwencje
‘, że mantysa należy do przedziału
M ∈ [1, 10) , (1.23)
natomiast cecha C jest liczbą całkowitą (dodatnią ujemna lub równą zeru). Dla przykładowej liczby x: M = 4.5 i C = −4.
Przyjmuja
‘c za podstawe
‘ p = 2 dostajemy analogiczny zapis w systemie dwój- kowym
x = ± M · 2C, (1.24)
gdzie tym razem dla mantysy zachodzi:
M ∈ [1, 2) (1.25)
1.5.1 Podział bitów
Reprezentacja (1.24) jest podstawa
‘ zapisu liczb rzeczywistych w komputerze.
W tym celu dzielimy n bitów przeznaczonych do zapisu liczby rzeczywistej,
x = S × M × 2C, (1.26)
w naste
‘puja
‘cy sposób
n = 1 + nM+ nC. (1.27)
W powyższej sumie
• pierwszy bit przeznaczamy na znak liczby S = ±1,
• naste
‘pne nM bitów służa
‘do zapisu cze‘ści ułamkowej F mantysy
M = 1 + F . (1.28)
Ze względu na zawsze występująca jedynkę nie musimy jej kodować. Część ułamkowa mantysy należy więc do przedziału
F ∈
1
2nM, 1 − 1 2nM
(1.29)
• pozostałe nC bitów przeznaczamy na zapis cechy C wraz ze znakiem.
Zgodnie z warunkiem (1.19)
C ∈ [ −2nC−1+ 1 , 2nC−1] . (1.30) Mantysę C można zapisać wykorzystując przesuniecie C = C0− E, gdzie E = 2nC−1− 1. Kodujemy wtedy liczbę całkowitą C0∈ [ 0 , 2nC− 1 ].
Oczywistym jest, że prawie wszystkie liczby rzeczywiste są reprezentowane w przybliżeniu.
Współczesne procesory reprezentują liczby zmiennoprzecinkowe w konwencji zgodnej ze standardem IEEE Standard 754-1985. Podział bitów dla cechy i mantysy dla tego standardu podany jest w tabelce.
Precyzja Liczba bitów nM nC Zakres C
Pojedyncza 32 23 8 [ -127 , 128 ]
Podwójna 64 52 11 [ -1023 , 1024 ]
Poniżej podajemy kilka zapisów liczb w pojedynczej precyzji.
0( 0 . . . 0
| {z }
23 zera
)(00000000) = + 1.0 × 20= 1.0
0( 0 . . . 0
| {z }
23 zera
)(00000001) = + 1.0 × 21= 2.0
1(1 0 . . . 0
| {z }
22 zera
)(10000001) = − (1.0 +12) × 2−1= −0.75
Wykład 2
Błe ¸dy i ich źródła
2.1 Dokładność zapisu
Miara
‘ przybliżonego zapisu liczby rzeczywistej jest epsilon maszynowy, zdefi- niowany w naste
‘puja
‘cy sposób.
Rozważmy najmniejszą liczbe
‘ możliwa
‘ do zapisu w komputerze wie
‘ksza od 1 i oznaczny ja
‘ przez x. Epsilon maszynowy to różnica:
m= x − 1 . (2.1)
W przykładzie z podziałem 8 = 1 + 3 + 4 bitów najmniejsza
‘ zapisana
‘liczba
‘
wie‘ksza
‘od 1 jest
(0)001(0)000 = (1 + 0 · 2−1+ 0 · 2−2+ 1 · 2−3) · 20 = 1 +18 (2.2) i sta
‘d m= 2−3= 0.125.
W ogólności, przeznaczaja
‘c nM bitów do zapisu mantysy (z wyła
‘czeniem znaku), zachodzi
m = 1
2nM (2.3)
Liczba
p = nM (2.4)
nazywana jest precyzja‘i jest liczba
‘ bitów przeznaczonych do zapisu mantysy z wyłaczeniem znaku.
Zwie‘kszaja
‘c liczbe
‘ bitów mantysy nM zwie
‘kszamy precyzje
‘. Odbywa sie to kosztem zmniejszenia maksymalnej liczby możliwej do zapisania, gdyż przy‘
ustalonej całkowitej liczbie bitów n, mniej bitów można przeznaczyć do zapisu cechy:
nC= n − 1 − nM. (2.5)
Jest to ilustracja wymienności pomie
‘dzy dokładnościa‘ a zakresem reprezento- wanych liczb zmiennoprzecinkowych w komputerze.
Tabelka poniżej precyzję i epsilon maszynowy dla używanych powszechnie precyzji zapisu liczby rzeczywistej.
Precyzja Liczba bitów nM p m
Pojedyncza 32 23 23 2−23∼ 10−7
Podwójna 64 52 52 2−52∼ 10−15
2.2 Błąd bezwględny i względny
Jak widzimy, na ogół dowolona liczba zmiennoprzecinkowa x nie może być re- prezentowana w komputerze dokładnie.
Jeżeli najbliższa
‘liczba
‘reprezentującą x w komputerze jest x0, to bła‘d bez- wgle‘dny zapisu liczby x jest zdefiniowany jako różnica
δx = x0− x (2.6)
Bła‘d wzgle‘dny to stosunek błe
‘du bezwzgle
‘dnego do wartości dokładnej liczby
= x0− x
x (2.7)
Z powyższego wzoru wynika. że liczbę przybliżoną x0 można traktować ja- ko zaburzenie liczby dokładnej x przez czynnik multiplikatywny zawierający wartość błędu względnego ,
x0= x (1 + ). (2.8)
Można pokazć, że moduł błe
‘du względnego zapisu liczby zmiennoprzecinkowej jest mniejszy od epsilona maszynowego
|| ¬ m (2.9)
2.2.1 Przyklad Zapiszmy liczbe
‘ 0.48 = 12/25 przy pomocy 8 = 1 + 3 + 4 bitów. Epsilon maszy- nowy w tym przypadku to
m=18= 0.125 . (2.10)
W reprezentacji cecha-mantysa otrzymujemy
0.48 = 1.92/4 = (1.92) · 2−2. Mantysa zapisana w systemie dwójkowym to
1.92 = (1.111010 . . .)2. (2.11) W przyje
‘tej reprezentacji pozostawiamy pierwsze trzy bity w rozwinie
‘ciu dwój- kowym mantysy, sta
‘d
0.48 ≈ (0)111(1)010 = (1 + 1 · 2−1+ 2−2+ 2−3) · 2−2=1532= 0.46875 . Jest to ilustracja zaokra‘glania, polegaja
‘cego na odrzuceniu cyfr w rozwinie
‘- ciu dwójkowym mantysy wykraczaja
‘cych poza przyje
‘ta
‘ liczbe
‘ bitów (w tym przypadku trzech). Powstały w ten sposób bła
‘d wzgle
‘dny to
|| =
15 32−1225
12 25
= 964 ≈ 0.04 . Jak widać || < m, zgodnie z warunkiem (2.9).
2.3 Źródła błe ¸dów
Przy obliczeniach komputerowych mamy do czynienia z trzema rodzajami błe
‘dów.
2.3.1 Błędy wejściowe.
Sa‘ one spowodowane tym, że dane wprowadzane do komputera sa
‘ obarczone błe‘dem. Moga
‘to być dane eksperymentalne, które ze swej natury sa
‘obarczone błe‘dem. Najcze
‘ściej jednak błe
‘dy te wynikaja
‘ z zaokra
‘glania liczb rzeczywi- stych, reprezentowanych w komputerze przez słowa binarne o skończonej dłu- gości.
2.3.2 Błe¸dy obcie¸cia.
Błe‘dy obcie
‘cia powstaja
‘ podczas obliczeń wymagaja
‘cych nieskończonej liczby działań, np. podczas nieskończonych sumowań lub ścisłych przejść granicznych.
Dla przykładu, obliczaja
‘c wartość funkcji przy pomocy rozwinie
‘cia Taylora, f (x) = f (x0) +f0(x0)
1! h +f00(x0)
2! h2+ . . . +f(N )(x0)
N ! hN+ RN +1(ξ) (2.12) gdzie h = x − x0, popełniamy bła
‘d zadany przez odrzucana
‘reszte
‘
RN +1(ξ) = f(N +1)(ξ)
(N + 1)! hN +1, (2.13)
gdzie punkt ξ ∈ [min{x0, x}, max{x0, x}].
W szczególności, dla kilku podstawowych funkcji używamy naste
‘puja
‘cych
przybliżonych rozwinie
‘ć wokół x0= 0:
exp x ≈ 1 + x +x2
2! + . . . +xN
N ! (−∞ < x < ∞) (2.14)
sin x ≈ x −x3 3! +x5
5!+ . . . + (−1)N x2N +1
(2N + 1)! (−∞ < x < ∞) (2.15) cos x ≈ 1 −x2
2! +x4
4!+ . . . + (−1)N x2N
(2N )! (−∞ < x < ∞) (2.16) 1
1 − x ≈ 1 + x + x2+ . . . + xN, (−1 < x < 1) (2.17) Ponadto, do obliczenia logarytmu dowolnej dodatniej liczby rzeczywistej dodat- niej stosujemy wzór
ln1 + x
1 − x ≈ 2 x + x3 3 + x5
5 + . . . + x2N +1 2N + 1
!
, (−1 < x < 1) (2.18) Zauważmy bowiem, że podstawiaja
‘c x ∈ (−1, 1), liczba w =1 + x
1 − x (2.19)
przebiega zakres wartości dziedziny logarytmu: w ∈ (0, +∞). Sta
‘d metoda by licza
‘c ln w najpierw odwrócić powyższa
‘ relacje
‘: x =(w − 1)
(w + 1), (2.20)
a naste
‘pnie podstawić ja
‘ po prawej stronie wzoru (2.18).
2.3.3 Błe¸dy zaokra¸gleń.
To błe
‘dy wynikaja
‘ce z zaokra
‘gleń wykonywanych w trakcie wielokrotnych dzia- łań arytmetycznych.
Oznaczmy przez fl(x) wynik obliczeń numerycznych różnia
‘cy sie
‘ od dokład- nej wartości x na skutek błe
‘dów zaokra
‘gleń. Zgodnie z (2.7) dla pojedynczej liczby mamy:
fl(x) = x (1 + ) . (2.21)
W przypadku działań arytmetycznych obowia
‘zuje:
Lemat Wilkinsona Błe‘dy zaokra
‘gleń powstaja
‘ce przy wykonywaniu działań zmiennoprzecinkowych
sa‘równoważne zaste
‘pczemu zaburzaniu liczb, na których wykonujemy działania.
W przypadku pojedynczych działań zachodzi
fl(x1± x2) = x1(1 + 1) ± x2(1 + 2) (2.22) fl(x1· x2) = x1(1 + 3) · x2 = x1· x2(1 + 3) (2.23)
fl(x1/x2) = x1(1 + 4)/x2 = x1/(x2(1 + 5)) , (2.24) Wszystkie zaburzenia można wybrać tak by ich moduł był mniejszy od epsilona maszynowego
|i| < m. (2.25)
2.3.4 Przyklad
Zilustrujmy lemat Wilkinsona dodaja
‘c do siebie dwie liczby w reprezentacji jednobajtowej z podziałem 8 = 1 + 3 + 4,
0.48 + 0.24 =1225+1250= 0.72 . Dostajemy
0.48 = 1.92 · 2−2 ≈ (0)111(1)010 = (1 + 1 · 2−1+ 2−2+ 2−3) · 2−2=1532 0.24 = 1.92 · 2−3 ≈ (0)111(1)110 = (1 + 1 · 2−1+ 2−2+ 2−3) · 2−3=1564. gdyż
1.92 = (1.111010 . . .)2≈ (1.111)2. Po dodaniu numerycznym otrzymujemy
fl(0.48 + 0.24) =1532+1564=4564= 1.40625 · 2−1≈ (0)011(1)100 =1116= 0.6875 , gdyż
1.40625 = (1.01101)2≈ (1.011)2.
Zgodnie z lematem Wilkinsona, ten sam wynik można otrzymać zaburzaja
‘c dodawane do siebie liczby
12
25(1 + 1) +1250(1 + 2) =1116, ska‘d wynika relacja
21+ 2= −1396. Wybór liczb i nie wie
‘c jest jednoznaczny. Można je jednak wybrać tak by zachodziło: |i| ¬ m= 1/8. Na przykład: 1= −1/16 oraz 2= −1/96.
2.4 Odejmowanie małych liczb
Obliczmy bła
‘d wzgle
‘dny różnicy dwóch liczb
|| =
fl(x1− x2) − (x1− x2) x1− x2 .
(2.26) Zgodnie z lematem Wilkinsona
fl(x1− x2) = x1(1 + 1) − x2(1 + 2) . (2.27) Podstawiając do pierwszego wzoru, otrzymujemy
|| =
x1· 1− x2· 2 x1− x2 .
. (2.28)
Załóżmy, ze dwie liczby różnią się mało od siebie
x1= x2(1 + δ) , δ → 0 (2.29)
Podstawiając powyżej, znajdujemy
|| ' |1− 2|
|δ| → ∞ dla δ → 0 . (2.30)
Błąd względny odejmowania może być bardzo duży dla liczb niewiele różniących się od siebie. Należy wie‘c unikać odejmowania mało różnia‘cych sie‘ od siebie liczb.
Nieograniczony wzrost (2.30) nie jest sprzeczny z warunkiem (2.25) w le- macie Wilkinsona, gdyż epsilony w tym lemacie to zaburzenia indywidualnych liczb, a nie bła
‘d wzgle
‘dny (2.26) wyniku działań arytmetycznych.
Jako przykład, obliczmy jeden z pierwiastków równania kwadratowego dla b2 4ac,
x =−b +√
b2− 4ac
2a . (2.31)
W tym celu należy wykorzystać równoważna
‘ postać x =−b +√
b2− 4ac
2a ·(−b −√
b2− 4ac) (−b −√
b2− 4ac) = −2c b +√
b2− 4ac, (2.32) gdyz dodajemy wtedy dwie wielkości tego samego rze
‘du w mianowniku, zamiast odejmować dwie bliskie liczby w liczniku.
Wykład 3
Algorytmy
Omówimy po krótce problem stabilności i efektywności algorytmów. Zagad- nienie to dotyczy wpływu błe
‘dów zaokra
‘gleń lub błe
‘dów obcie
‘ć na stabilność algorytmów numerycznych. Wynikaja
‘ce sta
‘d problemy przedstawimy na przy- kładzie algorytmów realizowanych przy pomocy relacji rekurencyjnych.
3.1 Algorytmy stabilne
Przykładem algorytmu stabilnego jest relacja rekurencyjna służa
‘ca do oblicze- nia kolejnych przybliżeń liczby√
2:
xn+1=1 2
xn+ 2
xn
, n 0 (3.1)
Zaczynaja
‘z x0= 2 otrzymujemy kolejne przybliżenia
x1= 1.5 , x2= 1.416667 , x3= 1.414216 Zakładaja
‘c, że istnieje granica
‘ tego cia
‘gu otrzymujemy jako granice
‘ x∗=√
2, gdyż z relacji (3.1) wynika
x∗=x∗
2 + 1 x∗
=> x∗=
√
2 . (3.2)
Istnienie granicy wynika z faktu, że bład wzgle
‘dny kolejnych przybliżeń szybko zmierza do zera. Zaczynaja
‘c od n = 1, otrzymujmy x1=
√
2 (1 + ) . (3.3)
gdzie < 1. W kolejnym przybliżeniu otrzymujemy x2 = x1
2 + 1 x1 = 1
√ 2
1 + + 1 1 +
= 1
√2
n1 + + (1 − + 2+ . . .)o≈√ 2
( 1 +2
2 )
(3.4)
Teraz bła
‘d wzgle
‘dny jest kwadratem poprzedniego. W naste
‘pnych krokach bład ten maleje jak kolejne pote
‘gi 2ktak, że wystarczy niewiele iteracji do osia
‘gniecia z góry założonej dokładności Algorytm znajdowania pierwiastka jest wie
‘c szyb- ko zbieżny (efektywny).
Przykładem mało efektywnego algorytmu jest sposób obliczania liczby ln 2 przy pomocy rozwinie
‘cia w szereg Taylora ln(1 + x) ≈ x −x2
2 +x3
3 + . . . + (−1)N −1xN
N , (3.5)
dla x = 1. Bła
‘d bezwgle
‘dny to |RN +1| < 1/(N + 1) i chca
‘c na przykład osia
‘gna
‘c dokładność równa
‘10−8, musimy wysumować N ≈ 108 wyrazów. Znacznie lepiej jest tutaj skorzystać ze wzoru (2.18)
ln1 + x
1 − x≈ 2 x + x3 3 +x5
5 + . . . + x2N +1 2N + 1
!
(−1 < x < 1) dla x = 1/3, który jest dużo szybciej zbieżny ze wzgle
‘du na szybko maleja
‘ce kolejne nieparzyste pote
‘gi tej liczby.
3.2 Algorytmy niestabilne
Rozważmy trójczłonowa
‘ relacje
‘ rekurencyjna
‘
xn+2 = xn− xn+1, n 0 (3.6)
z warunkiem pocza
‘tkowym
x0= 1 , x1=
√5 − 1
2 . (3.7)
Łatwo sie przekonać, że rozwia
‘zaniem tej rekurencji jest xn=
√ 5 − 1
2
!n
. (3.8)
W oczywisty sposób xn→ 0 dla n → ∞, gdyż (√
5 − 1)/2 ≈ 0.618 < 1 Obliczaja
‘c jednak kolejne wartości xn z warunku (3.6) dostajemy od pewnego n szyb- ko rosna
‘ce wartości. Jest to spowodowane nieuniknionym błe
‘dem zaokra
‘glenia liczby niewymiernej (√
5 − 1)/2.
Aby to zrozumieć, rozważmy relacje
‘ rekurencyjna
‘ (3.6) podstawiaja
‘c roz- wia‘zanie próbne w postaci
xn= λn, (3.9)
Otrzymujemy wtedy równanie
λ2+ λ − 1 = 0 . (3.10)
Ma ono dwa rozwia
‘zania λ1=
√5 − 1
2 , λ2= −
√5 + 1
2 , (3.11)
spełniaja
‘ce warunki: λ1< 1 oraz |λ2| ≈ 1.618 > 1. Sta
‘d ogólne rozwia
‘zanie re- kurencji w postaci
xn = A (λ1)n+ B (λ2)n. (3.12) Współczynniki A i B wyznaczymy korzytaja
‘c z warunków pocza
‘tkowych (3.7) A + B = 1
A λ1+ B λ2 = λ1+ , (3.13) gdzie dodaliśmy mały czynnik 6= 0 po prawej stronie. Zdaje on sprawe
‘ z tego, że liczba λ1 jest zawsze reprezentowana w komputerze w sposób przybliżony ze wzgle
‘du na bła
‘d zaokra
‘glenia nieskończonego rozwinie
‘cia dwójkowego tej liczby. Rozwia
‘zaniem układu równań (3.13) jest A = 1 +
√5, B = −
√5. (3.14)
Sta‘d ostateczne rozwia
‘zanie rekurencji xn =
1 +
√ 5
(λ1)n−
√
5(λ2)n. (3.15) Tak wiec, niezależnie od tego jak mały jest bła
‘d , drugi wyraz be
‘dzie do- minował w sumie od pewnej wartości n, gdyż |(λ2)n| → ∞ dla n → ∞. Relacja rekurencyjna (3.6) z warunkiem pocz¸atkowym (3.7) jest wie
‘c niestabilna ze wzgle
‘du na małe zaburzenia wartości pocza
‘tkowych.
Łatwo sprawdzić, że relacja rekurencyjna (3.6) z warunkiem
x0= 1 , x1= λ2, (3.16)
który prowadzi do rozwia
‘zania xn= (λ2)n, jest stabilna ze wzgle
‘du na małe zaburzenia warunku pocza
‘tkowego.
3.3 Złożoność obliczeniowa
Złożoność obliczeniowa jest określona przez liczbe
‘ kroków (działań), które na- leży wykonać w danym algorytmie by osia
‘gna
‘ć zamierzony wynik.
Jako przykład, rozważmy problem obliczenia wartości wielomianu stopnia n,
Wn(x) = a0+ a1x + a2x2+ . . . + anxn, (3.17)
dla pewnej wartości x. Licza
‘c "naiwnie" mamy do wykonania n dodawań oraz naste
‘puja
‘ca
‘liczbe
‘ mnożeń
1 + 2 + . . . + (n − 1) + n = n(n + 1)
2 (3.18)
W sumie liczba działań zależy kwadratowo od stopnia wielomianu n.
Dla wielomianu stopnia 10 wykonujemy wie
‘c 55 ∼ O(102) działań, nato- miast dla wielomianu stopnia 100 to liczba 5050 ∼ O(104). Dobrze by było zna- leźć algorytm, w którym liczba wykonanych działań byłaby znacza
‘co mniejsza, chociażby ze wgle
‘du ma ryzyko kumulacji błe
‘dów zaokra
‘gleń.
Takim algorytmem jest algorytm Hornera. Obliczmy wartość wielomianu zapisuja
‘c go w naste
‘puja
‘cy sposób:
Wn(x) = (. . . ((anx + an−1) x + an−2) x + an−3) x + . . . + a1) x + a0) . (3.19) Tym razem wyste
‘puje tu n mnożeń oraz n dodawań. Liczba działań zależy wie
‘c liniowo od n. Na przykład, dla wielomianu stopnia 100 wykonujemy tylko 200 działań!
Algorytmy, w których liczba kroków do wykonania by osia
‘gna
‘c cel rośnie szybciej niż pote
‘gowo z wymiarem problemu n (w naszym przykładzie określo- nym przez stopień wielomianu) nie jest praktycznie obliczalny. Takim wzrostem jest wzrost wykładniczy en. Już dla n = 10 mamy e10≈ 104 kroków, podczas gdy dla n = 100 liczba kroków to e100≈ 1043!
Przykładem takiego problemu jest rozkład danej liczby naturalnej na ilo- czyn liczb pierwszych. Wszystkie algorytmy implementowane na komputerze klasycznym zależa
‘ wykładniczo od długości liczby n, która określa wymiar problemu. Dopiero kwantowy algorytm Schora, implementowany na komputrze kwantowym, pozwala rozwia
‘zać ten problem w przy pomocy około (ln n)3 licz- by kroków. Problem tylko w tym, że jak dota
‘d nie skonstruowano stabilnie działaja
‘cego komputera kwantowego.
Wykład 4
Macierze
4.1 Układ równań liniowych
Macierze pojawiaja
‘sie
‘ naturalnie w problemie rozwia
‘zywania układów równań liniowych
a11x1+ a12x2 + ... + a1nxn = y1 a21x1+ a22x2 + ... + a2nxn = y2
...
an1x1 + an2x2+ ... + annxn = yn. (4.1) Układ ten można zapisać w postaci macierzowej
a11 a12 . . . a1n
a21 a22 . . . a2n ... ... ... an1 an2 . . . ann
x1
x2 ... xn
=
y1
y2 ... yn
. (4.2)
lub w skróconej formie.
A x = y . (4.3)
Tak wie
‘c A jest macierza
‘ kwadratowa
‘ o wymiarze n × n. Rozwia
‘zanie układu (4.1) istnieje i jest ono jednoznaczne jeśli wyznacznik
det A 6= 0 . (4.4)
Istnieje wtedy macierz odwrotna A−1 i rozwia
‘zanie jest zadane przez
x = A−1y . (4.5)
Metoda ta wymaga wie
‘c znajomości metod odwracania macierzy A.
Alternatywnie, można skorzystać ze wzorów Cramera x1 = det A1
det A , x2 = det A2
det A , . . . xn = det An
det A . (4.6)
Macierze Ai w liczniku powstały poprzez zamiane
‘ i-tej kolumny w macierzy A przez wektor y. Przykładowo
A1 =
y1 a12 . . . a1n y2 a22 . . . a2n
... ... ... yn an2 . . . ann
. (4.7)
Z punktu widzenia metod numerycznych wzory Cramera można zastosować jedynie dla macierzy o małych wymiarach. Obliczenie wyznacznika macierzy o wymiarze n × n wymaga bowiem wykonanie n! mnożeń oraz n! dodawań.
We wzorach (4.6) musimy policzyń (n + 1) wyznaczników macierzy, sta
‘d cał- kowita liczba działań potrzebnych do znalezienia rozwia
‘zania wynosi 2 (n + 1)!.
Korzystaja
‘c ze wzoru Stirlinga, słusznego dla dużych n, n! ≈ nne−n√
2π n ∼ en ln n−n, (4.8)
widzimy, że liczba działań w metodzie Cramera rośnie wykładniczo z n i metoda szybko staje sie
‘ bezużyteczna. Przykładowo, rozwia
‘zuja
‘c układ złożony z 15 równań musimy wykonać
2 · 16! ≈ 1013 (4.9)
działań.
Przykład ten ilustruje zagadnienie złożoności obliczeniowej. Wszystkie algorytmy, w których liczba działań do wykonania rośnie wykładniczo z wymia- rem problemu szybko staja
‘ sie
‘ bezużyteczne. Aby rozwia
‘zać problem musimy znaleźć bardziej wydajne algorytmy, w których liczba działań zależy na przy- kład pote
‘gowo od n.
4.2 Metoda eliminacji Gaussa
Metoda eliminacji Gaussa pozwala znaleźć rozwia
‘zanie układu równań linio- wych przy pote
‘gowo zależnej od n liczbie wykonanych działań. Pozwala ona sprowadzić układ (4.1) do równoważnej postaci z macierza
‘ trójka
‘tna
‘, w której wszystkie składowe poniżej (lub powyżej) diagonali sa
‘ równe zeru
a110 a120 . . . a1(n−1)0 a1n0 0 a220 . . . a2(n−1)0 a2n0 ... ... ... ... 0 0 . . . a(n−1)(n−1)0 a(n−1)n0
0 0 . . . 0 ann0
x1
x2 ... xn−1
xn
=
y10 y20 ... yn−10
yn0
. (4.10)
Wyznacznik macierzy trójka
‘tnej jest iloczynem liczb na diagonali det A0 =
n
Y
i=1
aii0 . (4.11)
Warunkiem istnienia rozwia
‘zania układu (4.10) jest by det A06= 0. Sta
‘d wszyst- kie wyrazy na diagonali sa
‘ rózne od zera. Łatwo już rowia
‘zać taki układ.
Zaczynaja
‘c od ostatniego równania dostajemy xn = 1
ann0 yn0
xn−1 = 1
a(n−1)(n−1)0
yn−10 − a(n−1)n0 xn
. . . xi = 1
aii0
yi0−
n
X
k=i+1
aik0 xk
. (4.12)
Rozważmy równanie (4.1):
A x = y .
W pierwszym kroku metody eliminacji Gaussa mnożymy pierwsze równanie w tym układzie przez ai1/a11 i odejmujemy kolejno od i-tych równań, gdzie i = 2, 3, . . . , n. W wyniku tego otrzymamy układ równań:
A(1)x = y(1), (4.13)
lub w jawnej postaci
a(1)11 x1+ a(1)12 x2+ ... + a(1)1nxn = y1(1) a(1)22 x2+ ... + a(1)2nxn = y2(1) ...
a(1)n2x2+ ... + a(1)nnxn = yn(1). (4.14) Wyeliminowaliśmy w ten sposób pierwsza
‘ niewiadoma
‘ w równaniach o nume- rach i 2.
Mnożymy naste
‘pnie drugie równanie w układzie (4.14) przez a(1)i2 /a(1)22 i odej- mujemy od i-tego równania dla i = 3, 4, . . . , n. Eliminujemy wie
‘c druga
‘ niewia- doma‘ x2 z równań o numerach i 3, otrzymuja
‘c nowy układ:
A(2)x = y(2), (4.15)
lub
a(2)11 x1 + a(2)12 x2+ a(2)13 x3+ ... + a(2)1nxn = y(2)1 a(2)22 x2+ a(2)23 x3+ ... + a(2)2nxn = y(2)2
...
a(2)n2x3+ ... + a(2)nnxn = y(2)n . (4.16) Procedure
‘ te
‘ kontynuujemy aż do chwili otrzymania układu równań z macierza trójka
‘tna
‘:
A(n−1)x = y(n−1), , (4.17)
lub
a(n−1)11 x1 + a(n−1)12 x2+ a(n−1)13 x3+ ... + a(n−1)1n xn = y1(n−1) a(n−1)22 x2+ a(n−1)23 x3+ ... + a(n−1)2n xn = y2(n−1)
...
a(n−1)nn xn = yn(n−1)(4.18). Naste
‘pnie rozwia
‘zujemy ten układ korzystaja
‘c ze wzorów (4.12).
Do wyznaczenia ta
‘metoda
‘ rozwia
‘zania wykonujemy, odpowiednio
1
3n3+ n2−13n , 13n3+12n2−56n (4.19) mnożeń i dzieleń oraz dodawań. Całkowita liczba działań zależy wie
‘c pote
‘gowo
od n.
4.3 Rozkład LU macierzy
Bardzo pożytecznym dla zastosowań numerycznych jest rozkład danej macierzy A na iloczyn dwóch macierzy trójka
‘tnych L i U:
A = L U (4.20)
takich, że
a11 a12 . . . a1n a21 a22 . . . a2n
... ... ... an1 an2 . . . ann
=
1 0 . . . 0 l21 1 . . . 0 ... ... ... ln1 ln2 . . . 1
u11 u12 . . . u1n 0 u22 . . . u2n
... ... ... 0 0 . . . unn
. (4.21)
W takim przypadku wyznaczenie rozwia
‘zania układu równań
A x = L(U x) = y (4.22)
polega na rozwia
‘zaniu dwóch układów równań liniowych z macierzami trójka
‘t- nymi
L w = y , U x = w , (4.23)
stosuja
‘c metode
‘ z poprzedniego rozdziału.
Metoda eliminacji Gaussa prowadzi do rozkładu LU, gdyż kolejne jej kroki sa‘ równoważne operacji mnożenia przez odpowiednie macierze L(i):
A(n−1) = L(n−1)L(n−2). . . L(1)A , (4.24) gdzie A(n−1) jest macierza
‘trójka
‘tna
‘ typu U. Podobnie dla wektora y:
y(n−1) = L(n−1)L(n−2). . . L(1)y . (4.25)
