Metody numeryczne dla fizyk ´
ow
Krzysztof Golec–Biernat
(June 5, 2007)
Wersja robocza nie do dystrybucji
Krak ´ow
2007
2
Spis tre´sci
1
Liczby i ich reprezentacje
6
1.1
Systemy liczbowe . . . . . . . . . . . . . . . . . . . . . . . . . .
6
1.2
Konwersja z systemu dziesie¸tnego na dw ´ojkowy . . . . . . . . . .
8
1.3
Zapis liczb całkowitych . . . . . . . . . . . . . . . . . . . . . . .
9
1.4
Zapis liczby rzeczywistych . . . . . . . . . . . . . . . . . . . . .
9
1.5
Dokładno´s´c zapisu . . . . . . . . . . . . . . . . . . . . . . . . .
11
2
Błe¸dy
12
2.1
Bła¸d zapisu . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.2
´
Zr´odła błe¸d´ow . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
3
Algorytmy
17
3.1
Algorytmy stabilne . . . . . . . . . . . . . . . . . . . . . . . . .
17
3.2
Algorytmy niestabilne
. . . . . . . . . . . . . . . . . . . . . . .
18
3.3
Zło˙zono´s´c obliczeniowa
. . . . . . . . . . . . . . . . . . . . . .
20
4
Macierze
22
4.1
Układ r´owna´n liniowych . . . . . . . . . . . . . . . . . . . . . .
22
4.2
Metoda eliminacji Gaussa . . . . . . . . . . . . . . . . . . . . . .
24
4.3
Rozkład LU macierzy . . . . . . . . . . . . . . . . . . . . . . . .
26
4.4
Metoda Doolittle’a . . . . . . . . . . . . . . . . . . . . . . . . .
27
4.5
Wyznacznik macierzy . . . . . . . . . . . . . . . . . . . . . . . .
29
3
5
Interpolacja
30
5.1
Wielomian interpolacyjny Lagrange’a . . . . . . . . . . . . . . .
30
5.2
Interpolacja Lagrange’a znanej funkcji . . . . . . . . . . . . . . .
32
5.3
Wielomiany Czebyszewa . . . . . . . . . . . . . . . . . . . . . .
33
5.4
Optymalny wyb´or we¸zł´ow interpolacji . . . . . . . . . . . . . . .
35
5.5
Wzory dla dowolnego przedziału . . . . . . . . . . . . . . . . . .
36
5.6
Interpolacja przy pomocy funkcji sklejanych . . . . . . . . . . . .
36
6
Aproksymacja
38
6.1
Regresja liniowa . . . . . . . . . . . . . . . . . . . . . . . . . . .
38
6.2
Aproksymacja ´sredniokwadratowa . . . . . . . . . . . . . . . . .
40
6.3
Aproksymacja Czebyszewa . . . . . . . . . . . . . . . . . . . . .
41
6.4
Aproksymacja Czebyszewa w dowolnym przedziale . . . . . . . .
43
6.5
Aproksymacja trygonometryczna . . . . . . . . . . . . . . . . . .
44
6.6
Wzory dla dowolnego okresu . . . . . . . . . . . . . . . . . . . .
46
7
R ´o˙zniczkowanie
48
7.1
Metoda z aproksymacja¸ . . . . . . . . . . . . . . . . . . . . . . .
48
7.2
Metody z rozwinie¸ciem Taylora
. . . . . . . . . . . . . . . . . .
48
7.3
Wie¸ksza dokładno´s´c
. . . . . . . . . . . . . . . . . . . . . . . .
49
7.4
Wy˙zsze pochodne . . . . . . . . . . . . . . . . . . . . . . . . . .
50
8
Całkowanie
51
8.1
Kwadratury . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
8.2
Metoda prostokat´ow
. . . . . . . . . . . . . . . . . . . . . . . .
52
8.3
Metoda trapez´ow . . . . . . . . . . . . . . . . . . . . . . . . . .
52
8.4
Metoda parabol Simpsona
. . . . . . . . . . . . . . . . . . . . .
53
8.5
Bła¸d przybli˙ze´n . . . . . . . . . . . . . . . . . . . . . . . . . . .
54
4
9
Kwadratury Gaussa
56
9.1
Rza¸d kwadratury . . . . . . . . . . . . . . . . . . . . . . . . . .
57
9.2
Wielomiany ortogonalne . . . . . . . . . . . . . . . . . . . . . .
58
9.3
Kwadratura Gaussa . . . . . . . . . . . . . . . . . . . . . . . . .
59
9.4
Wsp´ołczynniki kwadratury Gaussa . . . . . . . . . . . . . . . . .
60
9.5
Przykład kwadratury Gaussa . . . . . . . . . . . . . . . . . . . .
61
10 Zera funkcji
62
10.1 Metoda połowienia . . . . . . . . . . . . . . . . . . . . . . . . .
62
10.2 Metoda Newtona . . . . . . . . . . . . . . . . . . . . . . . . . .
64
10.3 Metoda siecznych . . . . . . . . . . . . . . . . . . . . . . . . . .
67
10.4 Metoda fałszywej prostej (falsi)
. . . . . . . . . . . . . . . . . .
67
11 R ´ownania r´o˙zniczkowe
70
11.1 R ´ownania zwyczajne pierwszego rze¸du
. . . . . . . . . . . . . .
70
11.2 Układ r´owna´n r´o˙zniczkowych
. . . . . . . . . . . . . . . . . . .
71
11.3 R ´ownania wy˙zszych rze¸d´ow . . . . . . . . . . . . . . . . . . . .
71
11.4 Metody Eulera
. . . . . . . . . . . . . . . . . . . . . . . . . . .
73
11.5 Metoda przewid´z i popraw . . . . . . . . . . . . . . . . . . . . .
74
11.6 Stabilno´s´c numeryczna rozwia¸za´n . . . . . . . . . . . . . . . . .
74
11.7 R ´ownania typu stiff . . . . . . . . . . . . . . . . . . . . . . . . .
75
12 Metoda Rungego-Kutty
77
12.1 Metoda drugiego rze¸du . . . . . . . . . . . . . . . . . . . . . . .
77
12.2 Metoda czwartego rze¸du . . . . . . . . . . . . . . . . . . . . . .
79
13 Metody Monte Carlo
80
13.1 Zmienna losowa i rozkład prawdopodobie´nstwa . . . . . . . . . .
80
13.2 Przykłady rozkład´ow prawdopodobie´nstwa
. . . . . . . . . . . .
82
13.3 Generowanie liczb losowych . . . . . . . . . . . . . . . . . . . .
83
13.4 Całkowanie metoda¸ Monte Carlo . . . . . . . . . . . . . . . . . .
84
5
Rozdział 1
Liczby i ich reprezentacje
1.1
Systemy liczbowe
Przypomnijmy zapis liczby całkowitej w systemie dziesie
‘
tnym, na przykład
(1234)
10
= 4 · 10
0
+ 3 · 10
1
+ 2 · 10
3
+ 1 · 10
4
.
Do zapisu u˙zywamy dziesie
‘
´c cyfr
{0,1,2,... ,9}, a 10 to podstawa systemu licz-
bowego. Warto´s´c cyfry zale˙zy od pozycji w liczbie.
W og´olno´sci, dowolna
‘
liczbe
‘
całkowita
‘
mo˙zna zapisa´c w systemie pozycyjnym
o podstawie p przy pomocy cyfr c
i
∈ {0,1,... , p − 1}:
(c
k
. . . c
1
c
0
)
p
= c
0
· p
0
+ c
1
· p + ... + c
k
· p
k
(1.1)
Szczeg´olnie wa˙zna
‘
role
‘
odgrywa system dw ´ojkowy o podstawie p
= 2 i cyfrach
c
i
∈ {0,1}. Przykładowo
(11010)
2
= 0 · 2
0
+ 1 · 2
1
+ 0 · 2
2
+ 1 · 2
3
+ 1 · 2
4
= 2 + 8 + 16 = (26)
10
Znaczenie tego systemu dla maszyn cyfrowych wynika z faktu, ˙ze podstawowa
‘
jednostka
‘
informacji jest bit przyjmuja
‘
cy jedna
‘
z dw ´och mo˙zliwych warto´sci: 0
lub 1.
6
Najmniejsza
‘
liczba
‘
całkowita
‘
dodatnia
‘
zapisana
‘
przy pomocy n bit´ow jest zero
0
= (00 . . . 00), natomiast najwie
‘
ksza
‘
liczba
‘
jest
(11 . . . 11
|
{z
}
n
−razy
)
2
= 1 · 2
0
+ 1 · 2
1
+ . . . + 1 · 2
n
−1
= 2
n
− 1.
(1.2)
Wykorzystali´smy tutaj wz´or na sume
‘
sko´nczonej liczby wyraz´ow w szeregu geom-
etrycznym o ilorazie a
6= 1:
1
+ a + a
2
+ . . . + a
N
−1
=
1
− a
N
1
− a
(1.3)
Podsumowuja
‘
c, dysponuja
‘
c n bitami mo˙zemy zapisa´c wszystkie liczby całkowite
dodatnie od 0 do 2
n
− 1.
Podobnie dowolna
‘
liczbe
‘
ułamkowa
‘
z przedziału
(0, 1) mo˙zemy zapisa´c w sys-
temie dw ´ojkowych przy pomocy (na og´oł niesko´nczonego) rozwinie
‘
cia
(0. c
−1
c
−2
. . . c
−n
)
2
= c
−1
· 2
−1
+ c
−2
· 2
−2
+ . . . + c
−n
· 2
−n
(1.4)
Przykładowo
(0.1101)
2
= 1 ·
1
2
+ 1 ·
1
4
+ 0 ·
1
8
+ 1 ·
1
16
=
13
16
.
Najwie
‘
ksza
‘
liczba
‘
ułamkowa
‘
zapisana
‘
przy pomocy n bit´ow jest
(0. 11 . . . 11
|
{z
}
n
−razy
)
2
= 1 · 2
−1
+ 1 · 2
−2
+ . . . + 1 · 2
−n
= 1 −
1
2
n
.
(1.5)
Wygodnie jest zapamie
‘
ta´c naste
‘
puja
‘
ca tabelke
‘
ułatwiaja
‘
ca
‘
zamiane
‘
liczby z sys-
temu dw ´ojkowego na dziesie
‘
tny.
n
2
n
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
7
1.2
Konwersja z systemu dziesie¸tnego na dw´ojkowy
Konwersje
‘
odwrotna
‘
(10
→ 2) dla liczb całkowitych realizuje sie
‘
zauwa˙zaja
‘
c,
x
= c
0
+ c
1
· 2
1
+ . . . + c
k
· 2
k
= c
0
+ 2 · (c
1
+ . . . + c
k
· 2
k
−1
)
(1.6)
Tak wie
‘
c c
0
to reszta z dzielenia x przez 2. Podobnie c
1
jest reszta
‘
z dzielenia
(c
1
+ . . . + c
k
· 2
k
−1
) przez 2, itd.
Przykładowo, konwertuja
‘
c 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´ory otrzymujemy:
(31)
10
= (11110)
2
.
Podobnie, dla liczb ułamkowych x
∈ (0,1), po pomno˙zeniu przez dwa otrzy-
mamy
2x
= c
−1
+ (c
−2
· 2
−1
+ . . . + c
−n
· 2
−n+1
)
|
{z
}
r
(1.7)
Jak łatwo zauwa˙zy´c c
−1
jest cze
‘
´scia
‘
całkowita
‘
liczby 2x (r´owna
‘
0 lub 1) natomiast
r
< 1 jest cze
‘
´scia
‘
ułamkowa
‘
. Mno˙zymy wtedy r przez dwa by wydoby´c c
−2
, 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
‘
´sci całkowite liczb po prawej stronie r´owno´sci od g´ory do dołu,
dostajemy
(0.375)
10
= (0.011)
2
. Otrzymana reprezentacja dw ´ojkowa jest sko´n-
czona. Na og´oł jednak rozwinie
‘
cie dw ´ojkowe ułamk´ow jest niesko´nczone, podob-
nie jak w systemie dziesie
‘
tnym.
Ła
‘
cza
‘
c dwie przedstawione metody konwersji 10
→ 2, mo˙zemy znale´z´c repre-
zentacje
‘
dw ´ojkowa
‘
dowolnej liczby rzeczywistej, na przykład:
(30.375)
10
= 30 + 0.375 = (11110.011)
2
(1.8)
8
1.3
Zapis liczb całkowitych
Zał´o˙zmy, ˙ze mamy do dyspozycji n bit´ow. Pierwszy z nich przeznaczamy na znak
liczby:
(0) dla znaku dodatniego i (1) dla znaku ujemnego. Pozostałe (n−1) bit´ow
słu˙zy do zapisu liczb całkowitych. z przedziału
(−2
n
−1
+ 1 , 2
n
−1
−1). Zauwa˙zmy,
˙ze poste
‘
puja
‘
c w ten spos´ob zero jest reprezentowane dwukrotnie ze wzgle
‘
du na
znak: 0
= (0)0 . . . 0 oraz 0 = (1)0 . . . 0. Przyjmuja
‘
c tylko pierwszy spso´ob zapisu
unikamy tej niejednozanczno´sci zwalniaja
‘
c jednocze´snie miejsce dla dodatkowej
liczby, na przykład r´ownej
−2
n
−1
. Tak wie
‘
c przy pomocy n bit´ow mo˙zna zapisa´c
2
n
liczb całkowitych z przedziału
(−2
n
−1
, 2
n
−1
− 1).
(1.9)
Przykładowo, dysponuja
‘
c 4 bajtami czyli 32 bitami (zapis w pojedynczej precyzji),
mo˙zemy zapisa´c wszystkie liczby całkowite pomie
‘
dzy:
(−2
31
, 2
31
− 1) = (−2 147 483 648, 2 147 483 647)
Operacje arytmetyczne na tych liczbach moga
‘
prowadzi´c do przekroczenia za-
kresu. Chwilowym wyj´sciem z problemu jest zwie
‘
kszenie liczby bajt´ow przezna-
czonych do zapisu jednej liczby całkowitej, np. do 8 bajt´ow czyli 64 bit´ow (zapis
w podw ´ojnej precyzji). Wtedy zakres wynosi:
(−2
63
, 2
63
− 1).
1.4
Zapis liczby rzeczywistych
Najbardziej efektywny zapis liczb rzeczywistych wykorzystuje reprezentacje
‘
zmie-
nnoprzecinkowa
‘
. W reprezentacji tej nie ustala sie
‘
maksymalnej liczby cyfr po
przecinku przy zapisie liczby.
Przykładowo, rozwa˙zmy liczbe
‘
x
= 0.00045. Mo˙zemy ja zapisa´c w naste
‘
puja
‘
cy
spos´ob
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˙zna wie
‘
c zapisa´c w postaci
x
= M · 10
C
.
(1.10)
9
Przyjmujemy konwencje
‘
, ˙ze mantysa zawiera sie
‘
w przedziale:
1
≤ |M| < 10.
(1.11)
Dla przykładowej liczby: M
= 4.5 i C = −4.
Przyjmuja
‘
c za podstawe
‘
p
= 2 dostajemy analogiczny zapis w systemie dw´oj-
kowym
x
= M · 2
C
,
(1.12)
gdzie tym razem dla mantysy zachodzi:
1
≤ |M| < 2.
(1.13)
Reprezentacja (1.12) jest podstawa
‘
zapisu liczb rzeczywistych w komputerze.
Dzielimy n bit´ow przeznaczonych do zapisu liczby rzeczywistej na n
M
bit´ow słu-
˙za
‘
cych do zapisu mantysy i n
C
bit´ow dla zapisu cechy, łacznie ze znakiem,
n
= n
M
+ n
C
.
(1.14)
Przykładowo, dla podziału 8
= 4 + 4 najmniejsza liczba dodatnia to
(0)000
| {z }
M
(1)111
| {z }
C
= (1 + 0 · 2
−0
+ 0 · 2
−1
+ 0 · 2
−2
) · 2
−(1·2
0
+1·2
1
+1·2
2
)
= 1 · 2
−7
=
1
128
,
gdzie liczby odczytujemy od lewa do prawa. Pierwsze bity mantysy i cechy (w
nawiasie) sa
‘
zarezerwowane na znak:
(0) = (+) oraz (1) = (−). Zauwa˙zmy, ˙ze
ze wzgle
‘
du na warunek M
≥ 1 nie musimy rezerwowa´c jednego bitu dla zapisu
jedynki w mantysie, pojawia sie
‘
ona z definicji. Liczba 1 jest reprezentowana
przez
(0)000(0)000 = 1 · 2
0
= 1 ,
natomiast liczba 0 nie jest reprezentowana i na jej zapis musi by´c przeznaczone
dodatkowe osiem bit´ow.
Najwie
‘
ksza liczba zapisana w ten spos´ob to
(0)111(0)111 = (1 + 1 · 2
−1
+ 1 · 2
−2
+ 1 · 2
−3
) · 2
(1·2
0
+1·2
1
+1·2
2
)
=
15
8
· 2
7
= 240 .
10
Okre´slona ona jest gł´ownie przez liczbe
‘
bit´ow cechy n
c
. Im wie
‘
ksza n
c
, tym
wie
‘
ksza
‘
liczbe
‘
maksymalna
‘
mo˙zna zapisa´c.
Podsumowuja
‘
c, przy podziale 8
= 4 + 4 bit´ow mo˙zemy zapisa´c 2
8
= 256 liczb
rzeczywistych z przedziału
−240, −
1
128
∪
1
128
, 240
.
Oczywistym jest, ˙ze prawie wszystkie liczby sa
‘
reprezentowane z przybli˙zeniem.
1.5
Dokładno´s´c zapisu
Miara
‘
przybli˙zonego zapisu liczby rzeczywistej jest epsilon maszynowy, zdefinio-
wany w naste
‘
puja
‘
cy spos´ob.
Rozwa˙zmy najmniejsza liczbe
‘
mo˙zliwa
‘
do zapisu w komputerze wie
‘
ksza od 1
i oznaczny ja
‘
przez x. Epsilon maszynowy to r´o˙znica:
ε
m
= x − 1.
(1.15)
W naszym przykładzie z podziałem 8
= 4 + 4 bit´ow najmniejsza
‘
zapisana
‘
liczba
‘
wie
‘
ksza
‘
od 1 to
(0)001(0)000 = (1 + 0 · 2
−1
+ 0 · 2
−2
+ 1 · 2
−3
) · 2
0
= 1 +
1
8
(1.16)
i sta
‘
d
ε
m
= 2
−3
= 0.125.
W og´olno´sci, przeznaczaja
‘
c n
M
bit´ow do zapisu mantysy (z czego jeden bit jest
po´swie
‘
cony na znak), zachodzi
ε
m
=
1
2
(n
M
−1)
(1.17)
Liczba p
= n
M
− 1 nazywana jest precyzja
‘
i jest liczba
‘
bit´ow przeznaczonych do
zapisu mantysy z wyłaczeniem znaku.
Zwie
‘
kszaja
‘
c liczbe
‘
bit´ow mantysy n
M
zwie
‘
kszamy precyzje
‘
. Odbywa sie
‘
to
kosztem zmniejszenia maksymalnej liczby mo˙zliwej do zapisania, gdy˙z przy usta-
lonej całkowitej liczbie bit´ow n, mniej bit´ow mo˙zna przeznaczy´c do zapisu cechy:
n
C
= n − n
M
.
(1.18)
Jest to ilustracja wymienno´sci pomie
‘
dzy dokładno´scia
‘
a zakresem reprezentowa-
nych liczb zmiennoprzecinkowych w komputerze.
11
Rozdział 2
Błe¸dy
2.1
Bła¸d zapisu
Na og´oł dowolona liczba zmiennoprzecinkowa x nie mo˙ze by´c reprezentowana w
komputerze dokładnie. Je˙zeli najbli˙zsza
‘
liczba
‘
zapisana
‘
dokładnie jest x
′
, to bła
‘
d
bezwgle
‘
dny jest zdefiniowany jako r´o˙znica
δ
x
= x
′
− x
(2.1)
Bła
‘
d wzgle
‘
dny to stosunek błe
‘
du bezwzgle
‘
dnego do warto´sci dokładnej liczby
ε
=
x
′
− x
x
(2.2)
Z ostatniego wzoru wynika
x
′
= x (1 +
ε
).
(2.3)
Mo˙zna pokaz´c, ˙ze moduł błe
‘
du wzgle
‘
dnego jest mniejszy od epsilona maszyno-
wego
|
ε
| ≤
ε
m
.
(2.4)
Przyklad
Chcemy zapisa´c liczbe
‘
0
.48 = 12/25 przy pomocy 8 = 4 + 4 bit´ow. Zapiszmy
ja
‘
w reprezentacji cecha-mantysa
0
.48 = 1.92/4 = (1.92) · 2
−2
.
12
Mantysa zapisana w systemie dw ´ojkowym to
1
.92 = (1.111010 . . .)
2
.
(2.5)
W przyje
‘
tej reprezentacji pozostawiamy pierwsze trzy bity w rozwinie
‘
ciu dw ´ojko-
wym mantysy, sta
‘
d
0
.48 ≈ (0)111(1)010 = (1 + 1 · 2
−1
+ 2
−2
+ 2
−3
) · 2
−2
=
15
32
= 0.46875 .
Jest to ilustracja zaokra
‘
glania, polegaja
‘
cego na odrzuceniu cyfr w rozwinie
‘
ciu
dw ´ojkowym mantysy wykraczaja
‘
cych poza przyje
‘
ta
‘
liczbe
‘
bit´ow (w tym przy-
padku trzech). Powstały w ten spos´ob bła
‘
d wzgle
‘
dny to
|
ε
| =
15
32
−
12
25
12
25
=
4
96
≈ 0.04.
Jak wida´c
|
ε
| < 1/8 =
ε
m
, zgodnie z warunkiem (2.4).
2.2
´
Zr´odła błe¸d´ow
Przy obliczeniach komputerowych mamy do czynienia z trzema rodzajami błe
‘
d´ow.
Błe
‘
dy wej´sciowe.
Sa
‘
one spowodowane tym, ˙ze dane wprowadzane do komputera sa
‘
obarczone
błe
‘
dem. Moga
‘
to by´c dane eksperymentalne, kt´ore ze swej natury sa
‘
obarczone
błe
‘
dem. Najcze
‘
´sciej jednak błe
‘
dy te wynikaja
‘
z zaokra
‘
glania liczb rzeczywistych,
reprezentowanych w komputerze przez słowa binarne o sko´nczonej długo´sci.
Błe
‘
dy obcie
‘
cia.
Błe
‘
dy obcie
‘
cia powstaja
‘
podczas oblicze´n wymagaja
‘
cych niesko´nczonej liczby
działa´n, np. podczas niesko´nczonych sumowa´n lub ´scisłych przej´s´c granicznych.
Dla przykładu, obliczaja
‘
c warto´s´c funkcji przy pomocy rozwinie
‘
cia Taylora,
f
(x) = f (x
0
) +
f
′
(x
0
)
1!
h
+
f
′′
(x
0
)
2!
h
2
+ . . . +
f
(N)
(x
0
)
N!
h
N
+ R
N
+1
(
ξ
)
(2.6)
gdzie h
= x − x
0
, popełniamy bła
‘
d zadany przez odrzucana
‘
reszte
‘
R
N
+1
(
ξ
) =
f
(N+1)
(
ξ
)
(N + 1)!
h
N
+1
,
(2.7)
13
gdzie punkt
ξ
∈ [min{x
0
, x},max{x
0
, x}].
W szczeg´olno´sci, dla kilku podstawowych funkcji u˙zywamy naste
‘
puja
‘
cych
przybli˙zonych rozwinie
‘
´c wok´oł x
0
= 0:
exp x
≈ 1 + x +
x
2
2!
+ . . . +
x
N
N!
(−
∞
< x <
∞
)
(2.8)
sin x
≈ x −
x
3
3!
+
x
5
5!
+ . . . + (−1)
N
x
2N
+1
(2N + 1)!
(−
∞
< x <
∞
)
(2.9)
cos x
≈ 1 −
x
2
2!
+
x
4
4!
+ . . . + (−1)
N
x
2N
(2N)!
(−
∞
< x <
∞
)
(2.10)
1
1
− x
≈ 1 + x + x
2
+ . . . + x
N
,
(0 < x < 1)
(2.11)
ln
1
+ x
1
− x
= 2
x
+
x
3
3
+
x
5
5
+ . . . +
x
2N
+1
2N
+ 1
(−1 < x < 1)
(2.12)
Z ostatniego wzoru korzystamy obliczaja
‘
c logarytm dowolnej dodatniej liczby
rzeczywistej. Zauwa˙zmy bowiem, ˙ze podstawiaja
‘
c x
∈ (−1,1) liczba
w
=
1
+ x
1
− x
(2.13)
przebiega zakres warto´sci dziedziny logarytmu: w
∈ (0,+
∞
). Sta
‘
d metoda by
licza
‘
c ln w najpierw odwr´oci´c powy˙zsza
‘
relacje
‘
:
x
=
(w − 1)
(w + 1)
,
(2.14)
a naste
‘
pnie podstawi´c ja
‘
po prawej stronie wzoru (2.12).
Błe
‘
dy zaokra
‘
gle´n.
To błe
‘
dy wynikaja
‘
ce z zaokra
‘
gle´n wykonywanych w trakcie wielokrotnych
działa´n arytmetycznych. Je´sli przez fl
(x) oznaczymy wynik oblicze´n numerycz-
nych r´o˙znia
‘
cy sie
‘
od dokładnej warto´sci x na skutek błe
‘
d´ow zaokra
‘
gle´n, to zgodnie
z (2.2) dla pojedynczej liczby mamy:
fl
(x) = x (1 +
ε
) .
(2.15)
W przypadku wykonywania działa´n arytmetycznych obowia
‘
zuje:
14
Lemat Wilkinsona
Błe
‘
dy zaokra
‘
gle´n powstaja
‘
ce przy wykonywaniu działa´n zmiennoprzecinko-
wych sa
‘
r´ownowa˙zne zaste
‘
pczemu zaburzaniu liczb, na kt´orych wykonujemy dzia-
łania. W przypadku pojedynczych działa´n zachodzi
fl
(x
1
± x
2
) = x
1
(1 +
ε
1
) ± x
2
(1 +
ε
2
)
(2.16)
fl
(x
1
· x
2
) = x
1
(1 +
ε
3
) · x
2
= x
1
· x
2
(1 +
ε
3
)
(2.17)
fl
(x
1
/x
2
) = x
1
(1 +
ε
4
)/x
2
= x
1
/(x
2
(1 +
ε
5
)) ,
(2.18)
Wszystkie zaburzenia mo˙zna tak dobra´c by zachodziło
|
ε
i
| <
ε
m
.
(2.19)
Przyklad
Zilustrujmy lemat Wilkinsona dodaja
‘
c do siebie dwie liczby w reprezentacji
8
= 4 + 4 bit´ow:
0
.48 + 0.24 =
12
25
+
12
50
= 0.72 .
Na podstawie przykładu z poprzedniego rozdziału:
0
.48 = 1.92 · 2
−2
≈ (0)111(1)010 = (1 + 1 · 2
−1
+ 2
−2
+ 2
−3
) · 2
−2
=
15
32
0
.24 = 1.92 · 2
−3
≈ (0)111(1)110 = (1 + 1 · 2
−1
+ 2
−2
+ 2
−3
) · 2
−3
=
15
64
.
Po dodaniu numerycznym otrzymujemy
fl
(0.48+0.24) =
15
32
+
15
64
=
45
64
= 1.40625·2
−1
≈ (0)011(1)100 =
11
16
= 0.6875 ,
gdy˙z
(1.40625)
10
= (1.01101)
2
. Zgodnie z lematem Wilkinsona, ten sam
wynik mo˙zna otrzyma´c zaburzaja
‘
c dodawane do siebie liczby
12
25
(1 +
ε
1
) +
12
50
(1 +
ε
2
) =
11
16
,
ska
‘
d wynika relacja
2
ε
1
+
ε
2
= −
13
96
.
Wyb´or liczb
ε
i
nie wie
‘
c jest jednoznaczny. Mo˙zna je jednak wybra´c tak by
zachodziło:
|
ε
i
| ≤
ε
m
= 1/8. Na przykład:
ε
1
= −1/16 oraz
ε
2
= −1/96.
15
Obliczmy bła
‘
d wzgle
‘
dny r´o˙znicy dw ´och liczb. Zgodnie z lematem Wilkinsona
|
ε
| =
fl
(x
1
− x
2
) − (x
1
− x
2
)
x
1
− x
2
=
x
1
·
ε
1
− x
2
·
ε
2
x
1
− x
2
.
(2.20)
Bła
‘
d wzgle
‘
dny mo˙ze by´c dowolnie du˙zy w sytuacji gdy odejmowane liczby niewiele
sie
‘
od siebie r´o˙znia
‘
. Na przkład dla x
1
= x
2
(1 +
δ
) otrzymujemy
|
ε
| ≃
|
ε
1
−
ε
2
|
|
δ
|
→
∞
dla
δ
→ 0.
(2.21)
Nale˙zy wie
‘
c unika´c odejmowania mało r´o˙znia
‘
cych sie
‘
od siebie liczb.
Nieograniczony wzrost (2.21) nie jest sprzeczny z warunkiem (2.19) w lemacie
Wilkinsona, gdy˙z epsilony w tym lemacie to zaburzenia indywidualnych liczb, a
nie bła
‘
d wzgle
‘
dny (2.20) wyniku działa´n arytmetycznych.
Przyklad
Obliczaja
‘
c jeden z pierwiastk´ow r´ownania kwadratowego dla b
2
≫ 4ac,
x
=
−b +
√
b
2
− 4ac
2a
,
(2.22)
nale˙zy wykorzysta´c r´ownowa˙zna
‘
posta´c
x
=
−b +
√
b
2
− 4ac
2a
·
(−b −
√
b
2
− 4ac)
(−b −
√
b
2
− 4ac)
=
−2c
b
+
√
b
2
− 4ac
.
(2.23)
Dodajemy wtedy dwie wielko´sci tego samego rze
‘
du w mianowniku, zamiast
odejmowa´c dwie bliskie liczby w liczniku.
16
Rozdział 3
Algorytmy
Om ´owimy po kr´otce problem stabilno´sci i efektywno´sci algorytm ´ow. Zagadnie-
nie to dotyczy wpływu błe
‘
d´ow zaokra
‘
gle´n lub błe
‘
d´ow obcie
‘
´c na stabilno´s´c algo-
rytm ´ow numerycznych. Wynikaja
‘
ce sta
‘
d problemy przedstawimy na przykładzie
algorytm ´ow realizowanych przy pomocy relacji rekurencyjnych.
3.1
Algorytmy stabilne
Przykładem algorytmu stabilnego jest relacja rekurencyjna słu˙za
‘
ca do obliczenia
kolejnych przybli˙ze´n liczby
√
2:
x
n
+1
=
1
2
x
n
+
2
x
n
,
n
≥ 0
(3.1)
Zaczynaja
‘
z x
0
= 2 otrzymujemy kolejne przybli˙zenia
x
1
= 1.5 ,
x
2
= 1.416667 ,
x
3
= 1.414216
Zakładaja
‘
c, ˙ze istnieje granica
‘
tego cia
‘
gu otrzymujemy jako granice
‘
x
∗
=
√
2, gdy˙z
z relacji (3.1) wynika
x
∗
=
x
∗
2
+
1
x
∗
=>
x
∗
=
√
2
.
(3.2)
Istnienie granicy wynika z faktu, ˙ze bład wzgle
‘
dny kolejnych przybli˙ze´n szy-
bko zmierza do zera. Zaczynaja
‘
c od n
= 1, otrzymujmy
x
1
=
√
2
(1 +
ε
) .
(3.3)
17
gdzie
ε
< 1. W kolejnym przybli˙zeniu otrzymujemy
x
2
=
x
1
2
+
1
x
1
=
1
√
2
1
+
ε
+
1
1
+
ε
=
1
√
2
1
+
ε
+ (1 −
ε
+
ε
2
+ . . .)
≈
√
2
(
1
+
ε
√
2
2
)
(3.4)
Teraz bła
‘
d wzgle
‘
dny jest kwadratem poprzedniego. Widzimy, ˙ze algorytm ten
jest bardzo efektywny, gdy˙z wystarczy niewiele krok´ow by uzyska´c bardzo dobra
‘
dokładno´s´c.
Przykładem mało efektywnego algorytmu jest spos´ob obliczania liczby ln 2
przy pomocy rozwinie
‘
cia w szereg Taylora
ln
(1 + x) ≈ 1 −
x
2
+
x
3
3
+ . . . + (−1)
N
−1
x
N
N
,
(3.5)
dla x
= 1. Bła
‘
d bezwgle
‘
dny to
|R
N
+1
| < 1/(N + 1) i chca
‘
c na przykład osia
‘
gna
‘
c
dokładno´s´c r´owna
‘
10
−8
, musimy wysumowa´c N
≈ 10
8
wyraz´ow. Znacznie lepiej
jest tutaj skorzysta´c ze wzoru (2.12) dla x
= 1/3, kt´ory jest du˙zo szybciej zbie˙zny
ze wzgle
‘
du na szybko maleja
‘
ce kolejne nieparzyste pote
‘
gi tej liczby.
3.2
Algorytmy niestabilne
Rozwa˙zmy liczbe
‘
φ
=
√
5
− 1
2
≈ 0.618033989.
(3.6)
Chcemy obliczy´c kolejne pote
‘
gi
φ
n
. Poniewa˙z
φ
< 1 to
lim
n
→
∞
φ
n
= 0
(3.7)
i kolejne warto´sci
φ
n
szybko maleja
‘
. Spr´obujmy uzyska´c ten wynik w inny spos´ob,
wykonuja
‘
c jedynie dodawania. Poniewa˙z
φ
2
= 1 −
φ
(3.8)
to po pomno˙zeniu obu stron przez
φ
n
, otrzymujemy naste
‘
puja
‘
ca
‘
relacje
‘
rekuren-
cyjna
‘
φ
n
+2
=
φ
n
−
φ
n
+1
.
(3.9)
18
Obliczaja
‘
c przy jej pomocy kolejne warto´sci
φ
n
dostajemy jednak od pewnego n
szybko rosna
‘
ce warto´sci, co jest sprzeczne z warunkiem (3.7). Jest to spowodowa-
ne nieuniknionym błe
‘
dem zaokra
‘
glenia liczby niewymiernej
φ
.
Aby to zrozumie´c, rozwa˙zmy og´olna
‘
relacje
‘
rekurencyjna
‘
F
n
+2
= F
n
− F
n
+1
(3.10)
z warunkami pocza
‘
tkowymi F
0
= 1 oraz F
1
=
φ
. Podstawiaja
‘
c rozwia
‘
zanie pr´obne
w postaci
F
n
=
λ
n
,
(3.11)
otrzymujemy r´ownanie
λ
2
+
λ
− 1 = 0.
(3.12)
Ma ono dwa rozwia
‘
zania
λ
−
=
√
5
− 1
2
,
λ
+
= −
√
5
+ 1
2
.
(3.13)
Zauwa˙zmy, ˙ze
|
λ
−
| < 1 natomiast |
λ
+
| > 1. Og´olne rozwia
‘
zanie rekurencji to
F
n
= A
λ
n
−
+ B
λ
n
+
,
(3.14)
gdzie wsp´ołczynniki A i B wyznaczymy korzytaja
‘
c z warunk´ow pocza
‘
tkowych
F
0
= A + B = 1
F
1
= A
λ
−
+ B
λ
+
=
φ
.
(3.15)
Rozwia
‘
zaniem jest
A
=
φ
−
λ
+
√
5
,
B
=
λ
−
−
φ
√
5
(3.16)
i sta
‘
d ostateczna posta´c rozwia
‘
zania rekurencji (3.10):
F
n
=
φ
−
λ
+
√
5
λ
n
−
+
λ
−
−
φ
√
5
λ
n
+
(3.17)
Zauwa˙zmy, ˙ze
φ
jako liczba niewymierna jest reprezentowana zawsze z pewnym
przybli˙zeniem
φ
=
λ
−
+
ε
.
(3.18)
Wtedy kolejny wyraz cia
‘
gu F
n
przyjmuje posta´c
F
n
=
1
+
ε
√
5
λ
n
−
−
ε
√
5
λ
n
+
.
(3.19)
19
Ze wzgle
‘
du na warunek
|
λ
+
| > 1, drugie wyra˙zenie dominuje dla du˙zych warto´sci
n, prowadza
‘
c do rozbie˙znego cia
‘
gu warto´sci F
n
. Gdyby
φ
=
λ
−
to wtedy zgodnie
z oczekiwaniem
F
n
=
λ
n
−
=
√
5
− 1
2
!
n
.
(3.20)
Jest to jednak niemo˙zliwe w praktyce.
3.3
Zło˙zono´s´c obliczeniowa
Zło˙zono´s´c obliczeniowa jest okre´slona przez liczbe
‘
krok´ow (działa´n), kt´ore nale˙zy
wykona´c w danym algorytmie by osia
‘
gna
‘
´c zamierzony wynik.
Jako przykład, rozwa˙zmy problem obliczenia warto´sci wielomianu stopnia n,
W
n
(x) = a
0
+ a
1
x
+ a
2
x
2
+ . . . + a
n
x
n
,
(3.21)
dla pewnej warto´sci x. Licza
‘
c ”naiwnie” mamy do wykonania n dodawa´n oraz
naste
‘
puja
‘
ca
‘
liczbe
‘
mno˙ze´n
(n + 1) + n + (n − 1) + ... + 1 =
(n + 1)(n + 2)
2
(3.22)
W sumie liczba działa´n zale˙zy kwadratowo od stopnia wielomianu n.
Dla wielomianu stopnia 10 wykonujemy wie
‘
c 76
∼ O(10
2
) działa´n, natomiast
dla wielomianu stopnia 100 to liczba 5251
∼ O(10
4
). Dobrze by było znale´z´c algo-
rytm, w kt´orym liczba wykonanych działa´n byłaby znacza
‘
co mniejsza, chocia˙zby
ze wgle
‘
du ma ryzyko kumulacji błe
‘
d´ow zaokra
‘
gle´n.
Takim algorytmem jest algorytm Hornera. Obliczmy warto´s´c wielomianu za-
pisuja
‘
c go w naste
‘
puja
‘
cy spos´ob:
W
n
(x) = (. . . ((a
n
x
+ a
n
−1
) x + a
n
−2
) x + a
n
−3
) x + . . . + a
1
) x + a
0
) .
(3.23)
Tym razem wyste
‘
puje tu n mno˙ze´n oraz n dodawa´n. Liczba działa´n zale˙zy wie
‘
c
liniowo od n. Na przykład, dla wielomianu stopnia 100 wykonujemy tylko 200
działa´n!
Algorytmy, w kt´orych liczba krok´ow do wykonania by osia
‘
gna
‘
c cel ro´snie
szybciej ni˙z pote
‘
gowo z wymiarem problemu n (w naszym przykładzie okre´slonym
przez stopie´n wielomianu) nie jest praktycznie obliczalny. Takim wzrostem jest
20
wzrost wykładniczy e
n
. Ju˙z dla n
= 10 mamy e
10
≈ 10
4
krok´ow, podczas gdy dla
n
= 100 liczba krok´ow to e
100
≈ 10
43
!
Przykładem takiego problemu jest rozkład danej liczby naturalnej na iloczyn
liczb pierwszych. Wszystkie algorytmy implementowane na komputerze klasycz-
nym zale˙za
‘
wykładniczo od długo´sci liczby n, kt´ora okre´sla wymiar problemu.
Dopiero kwantowy algorytm Schora, implementowany na komputrze kwantowym,
pozwala rozwia
‘
za´c ten problem w przy pomocy około
(ln n)
3
liczby krok´ow. Prob-
lem tylko w tym, ˙ze jak dota
‘
d nie skonstruowano stabilnie działaja
‘
cego komputera
kwantowego.
21
Rozdział 4
Macierze
4.1
Układ r´owna ´n liniowych
Macierze pojawiaja
‘
sie
‘
naturalnie w problemie rozwia
‘
zywania układ´ow r´owna´n
liniowych
a
11
x
1
+ a
12
x
2
+ ... + a
1n
x
n
= y
1
a
21
x
1
+ a
22
x
2
+ ... + a
2n
x
n
= y
2
...
a
n1
x
1
+ a
n2
x
2
+ ... + a
nn
x
n
= y
n
.
(4.1)
Układ ten mo˙zna zapisa´c w postaci macierzowej
a
11
a
12
. . . a
1n
a
21
a
22
. . . a
2n
..
.
..
.
..
.
a
n1
a
n2
. . . a
nn
x
1
x
2
..
.
x
n
=
y
1
y
2
..
.
y
n
.
(4.2)
lub w skr´oconej 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´sli wyznacznik
det A
6= 0.
(4.4)
22
Istnieje wtedy macierz odwrotna A
−1
i rozwia
‘
zanie jest zadane przez
x
= A
−1
y
.
(4.5)
Metoda ta wymaga wie
‘
c znajomo´sci metod odwracania macierzy A.
Alternatywnie, mo˙zna skorzysta´c ze wzor´ow Cramera
x
1
=
det A
1
det A
,
x
2
=
det A
2
det A
,
. . . x
n
=
det A
n
det A
.
(4.6)
Macierze A
i
w liczniku powstały poprzez zamiane
‘
i-tej kolumny w macierzy A
przez wektor y. Przykładowo
A
1
=
y
1
a
12
. . . a
1n
y
2
a
22
. . . a
2n
..
.
..
.
..
.
y
n
a
n2
. . . a
nn
.
(4.7)
Z punktu widzenia metod numerycznych wzory Cramera mo˙zna zastosowa´c je-
dynie dla macierzy o małych wymiarach. Obliczenie wyznacznika macierzy o
wymiarze n
×n wymaga bowiem wykonanie n! mno˙ze´n oraz n! dodawa´n. We wzo-
rach (4.6) musimy policzy´n
(n + 1) wyznacznik´ow macierzy, sta
‘
d całkowita liczba
działa´n potrzebnych do znalezienia rozwia
‘
zania wynosi 2
(n + 1)!. Korzystaja
‘
c ze
wzoru Stirlinga, słusznego dla du˙zych n,
n!
≈ n
n
e
−n
√
2
π
n
∼ e
n ln n
,
(4.8)
widzimy, ˙ze liczba działa´n w metodzie Cramera ro´snie wykładniczo z n. Metoda
szybko staje sie
‘
wie
‘
c bezu˙zyteczna. Przykładowo, rozwia
‘
zuja
‘
c układ zło˙zony z 15
r´owna´n musimy wykona´c
2
· 16! ≈ 10
13
(4.9)
działa´n.
Przykład ten ilustruje zagadnienie zło˙zono´sci obliczeniowej. Wszystkie algo-
rytmy, w kt´orych liczba działa´n do wykonania ro´snie wykładniczo z wymiarem
problemu szybko staja
‘
sie
‘
bezu˙zyteczne. Aby rozwia
‘
za´c problem musimy znale´z´c
bardziej wydajne algorytmy, w kt´orych liczba działa´n zale˙zy na przykład pote
‘
gowo
od n.
23
4.2
Metoda eliminacji Gaussa
Metoda eliminacji Gaussa pozwala znale´z´c rozwia
‘
zanie układu r´owna´n liniowych
przy pote
‘
gowo zale˙znej od n liczbie wykonanych działa´n. Pozwala ona sprowadzi´c
układ (4.1) do r´ownowa˙znej postaci z macierza
‘
tr´ojka
‘
tna
‘
, w kt´orej wszystkie skła-
dowe poni˙zej (lub powy˙zej) diagonali sa
‘
r´owne zeru
a
′
11
a
′
12
. . .
a
′
1
(n−1)
a
′
1n
0
a
′
22
. . .
a
′
2
(n−1)
a
′
2n
..
.
..
.
..
.
..
.
0
0
. . . a
′
(n−1)(n−1)
a
′
(n−1)n
0
0
. . .
0
a
′
nn
x
1
x
2
..
.
x
n
−1
x
n
=
y
′
1
y
′
2
..
.
y
′
n
−1
y
′
n
.
(4.10)
Wyznacznik macierzy tr´ojka
‘
tnej jest iloczynem liczb na diagonali
det A
′
=
n
∏
i
=1
a
′
ii
.
(4.11)
Warunkiem istnienia rozwia
‘
zania układu (4.10) jest by det A
′
6= 0. Sta
‘
d wszystkie
wyrazy na diagonali sa
‘
r´ozne od zera. Łatwo ju˙z rowia
‘
za´c taki układ. Zaczynaja
‘
c
od ostatniego r´ownania dostajemy
x
n
=
1
a
′
nn
y
′
n
x
n
−1
=
1
a
′
(n−1)(n−1)
y
′
n
−1
− a
′
(n−1)n
x
n
. . .
x
i
=
1
a
′
ii
y
′
i
−
n
∑
k
=i+1
a
′
ik
x
k
!
.
(4.12)
Rozwa˙zmy r´ownanie (4.1):
A x
= y .
W pierwszym kroku metody eliminacji Gaussa mno˙zymy pierwsze r´ownanie w
tym układzie przez a
i1
/a
11
i odejmujemy kolejno od i-tych r´owna´n, gdzie i
=
2
, 3, . . . , n. W wyniku tego otrzymamy układ r´owna´n:
A
(1)
x
= y
(1)
,
(4.13)
24
lub w jawnej postaci
a
(1)
11
x
1
+ a
(1)
12
x
2
+ ... + a
(1)
1n
x
n
= y
(1)
1
a
(1)
22
x
2
+ ... + a
(1)
2n
x
n
= y
(1)
2
...
a
(1)
n2
x
2
+ ... + a
(1)
nn
x
n
= y
(1)
n
.
(4.14)
Wyeliminowali´smy w ten spos´ob pierwsza
‘
niewiadoma
‘
w r´ownaniach o numerach
i
≥ 2.
Mno˙zymy naste
‘
pnie drugie r´ownanie w układzie (4.14) przez a
(1)
i2
/a
(1)
22
i odej-
mujemy od i-tego r´ownania dla i
= 3, 4, . . . , n. Eliminujemy wie
‘
c druga
‘
niewiado-
ma
‘
x
2
z r´owna´n o numerach i
≥ 3, otrzymuja
‘
c nowy układ:
A
(2)
x
= y
(2)
,
(4.15)
lub
a
(2)
11
x
1
+ a
(2)
12
x
2
+ a
(2)
13
x
3
+ ... + a
(2)
1n
x
n
= y
(2)
1
a
(2)
22
x
2
+ a
(2)
23
x
3
+ ... + a
(2)
2n
x
n
= y
(2)
2
...
a
(2)
n2
x
3
+ ... + a
(2)
nn
x
n
= y
(2)
n
.
(4.16)
Procedure
‘
te
‘
kontynuujemy a˙z do chwili otrzymania układu r´owna´n z macierza
tr´ojka
‘
tna
‘
:
A
(n−1)
x
= y
(n−1)
, ,
(4.17)
lub
a
(n−1)
11
x
1
+ a
(n−1)
12
x
2
+ a
(n−1)
13
x
3
+ ... + a
(n−1)
1n
x
n
= y
(n−1)
1
a
(n−1)
22
x
2
+ a
(n−1)
23
x
3
+ ... + a
(n−1)
2n
x
n
= y
(n−1)
2
...
a
(n−1)
nn
x
n
= y
(n−1)
n
. (4.18)
Naste
‘
pnie rozwia
‘
zujemy ten układ korzystaja
‘
c ze wzor´ow (4.12).
Do wyznaczenia ta
‘
metoda
‘
rozwia
‘
zania wykonujemy, odpowiednio
1
3
n
3
+ n
2
−
1
3
n
,
1
3
n
3
+
1
2
n
2
−
5
6
n
(4.19)
mno˙ze´n i dziele´n oraz dodawa´n. Całkowita liczba działa´n zale˙zy wie
‘
c pote
‘
gowo
od n.
25
4.3
Rozkład LU macierzy
Bardzo po˙zytecznym dla zastosowa´n numerycznych jest rozkład danej macierzy A
na iloczyn dw ´och macierzy tr´ojka
‘
tnych L i U:
A
= L U
(4.20)
takich, ˙ze
a
11
a
12
. . . a
1n
a
21
a
22
. . . a
2n
..
.
..
.
..
.
a
n1
a
n2
. . . a
nn
=
1
0
. . . 0
l
21
1
. . . 0
..
.
..
.
..
.
l
n1
l
n2
. . . 1
u
11
u
12
. . . u
1n
0
u
22
. . . u
2n
..
.
..
.
..
.
0
0
. . . u
nn
. (4.21)
W takim przypadku wyznaczenie rozwia
‘
zania układu r´owna´n
A x
= L(U x) = y
(4.22)
polega na rozwia
‘
zaniu dw ´och układ´ow r´owna´n liniowych z macierzami tr´ojka
‘
tny-
mi
L w
= y ,
U x
= w ,
(4.23)
stosuja
‘
c metode
‘
z poprzedniego rozdziału.
Metoda eliminacji Gaussa prowadzi do rozkładu LU, gdy˙z kolejne jej kroki sa
‘
r´ownowa˙zne operacji mno˙zenia 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´ojka
‘
tna
‘
typu U. Podobnie dla wektora y:
y
(n−1)
= L
(n−1)
L
(n−2)
. . . L
(1)
y
.
(4.25)
Kolejne macierze L to
L
(1)
=
1
0
0
. . . 0
−l
21
1
0
. . . 0
−l
31
0
1
. . . 0
..
.
..
.
..
.
..
.
−l
n1
0
0
. . . 1
l
i1
= a
i1
/a
11
,
i
> 1
(4.26)
26
naste
‘
pnie
L
(2)
=
1
0
0
. . . 0
0
1
0
. . . 0
0
−l
32
1
. . . 0
..
.
..
.
..
.
..
.
0
−l
n2
0
. . . 1
l
i2
= a
(1)
i2
/a
(1)
12
,
i
> 2
(4.27)
gdzie a
(1)
ik
to wsp´ołczynniki macierzy
A
(1)
= L
(1)
A
,
(4.28)
w kt´orej wyzerowane sa
‘
elementy pierwszej kolumny z wyja
‘
tkiem a
(1)
11
. Poste
‘
puja
‘
c
w ten spos´ob a˙z do macierzy L
(n−1)
otrzymujemy wz´or (4.24).
Macierze L
(i)
sa
‘
nieosobliwe, istnieje wie
‘
c dla ka˙zdej z nich macierz odwrotna
(L
(i)
)
−1
, kt´ora powstaje poprzez zamiane
‘
wyste
‘
puja
‘
cych w L
(i)
minus´ow na plusy.
Odwracaja
‘
c wie
‘
c wz´or (4.24), znajdujemy
A
= (L
(1)
)
−1
(L
(2)
)
−1
. . . (L
(n−1)
)
−1
|
{z
}
L
(A
(n−1)
)
|
{z
}
U
.
(4.29)
Otrzymujemy w ten spos´ob rozkład (4.21), gdy˙z A
(n−1)
jest macierza
‘
tr´ojka
‘
tna
‘
typu U , natomiast iloczyn kolejnych macierzy odwrotnych to macierz typu L:
L
=
1
0
. . . 0
l
21
1
. . . 0
..
.
..
.
..
.
l
n1
l
n2
. . . 1
.
(4.30)
4.4
Metoda Doolittle’a
Wsp´ołczynniki macierzy L i U mo˙zna tak˙ze znale´z´c bezpo´srednio traktuja
‘
c (4.21)
jako układ n
2
r´owna´n na n
2
poszukiwanych wsp´ołczynnik´ow l
i j
(i > j) oraz u
i j
(i ≤ j).
Rozwa˙zmy dla uproszczenia n
= 3
a
11
a
12
a
13
a
21
a
22
a
23
a
31
a
32
a
33
=
1
0
0
l
21
1
0
l
31
l
32
1
u
11
u
12
u
13
0
u
22
u
23
0
0
u
33
.
(4.31)
27
Mno˙za
‘
c pierwszy wiersz macierzy L przez kolumny macierzy U, otrzymujemy
u
11
= a
11
u
12
= a
12
u
13
= a
13
.
(4.32)
Mno˙za
‘
c drugi wiersz macierzy L przez kolumny macierzy U, znajdujemy
l
21
u
11
= a
21
,
l
21
u
12
+ u
22
= a
22
,
l
21
u
13
+ u
23
= a
23
,
ska
‘
d mo˙zna wyznaczy´c
l
21
= a
21
/u
11
u
22
= a
22
− l
21
u
12
u
23
= a
23
− l
21
u
13
.
(4.33)
Mno˙za
‘
c trzeci wiersz macierzy L przez kolumny macierzy U, dostajemy
l
31
u
11
= a
31
,
l
31
u
12
+ l
32
u
22
= a
32
,
l
31
u
13
+ l
32
u
23
+ u
33
= a
33
,
co prowadzi do
l
31
= a
31
/u
11
l
32
= (a
32
− l
31
u
12
)/u
22
u
33
= a
33
− l
31
u
13
− l
32
u
23
.
(4.34)
Łatwo sta
‘
d wydedukowa´c wz´or dla macierzy o dowolnym wymiarze n.
l
i j
=
1
u
j j
a
i j
−
j
−1
∑
k
=1
l
ik
u
k j
!
,
i
> j
(4.35)
u
i j
= a
i j
−
i
−1
∑
k
=1
l
ik
u
k j
,
i
≤ j
(4.36)
Zauwa˙zmy, ˙ze obliczaja
‘
c l
i j
wykorzystujemy elementy macierzy L poprzedzaja
‘
ce
ten element w wierszu oraz elementy macierzy U w wierszach powy˙zej. Podobnie,
dla obliczenia u
i j
wykorzystujemy obliczone ju˙z elementy macierzy L w wierszu
oraz elementy macierzy U w wierszach powy˙zej. W ten spos´ob powy˙zsze wzory
definiuja
‘
algorytm rozkładu LU , w kt´orym obliczamy kolejne elementy macierzy
L i U poruszaja
‘
c sie
‘
wzdłu˙z wierszy od lewa do prawa, z g´ory na d´oł.
Przedstawiona metoda nazywa sie
‘
metoda
‘
Doolittle’a. Rozwia
‘
zuja
‘
c przy jej
pomocy układ r´owna´n (4.1) wykonuje sie
‘
tyle samo działa´n algebraicznych co w
metodzie eliminacji Gaussa. Zale˙zno´s´c od n jest wie
‘
c pote
‘
gowa.
28
4.5
Wyznacznik macierzy
Rozkładu LU umo˙zliwia tak˙ze policzenie wyznacznika macierzy A, korzystaja
‘
c z
własno´sci wyznacznika iloczynu macierzy:
det A
= det L · det U
(4.37)
oraz z faktu, ˙ze wyznacznik macierzy tr´ojka
‘
tnej jest iloczynem liczb na diagonali.
Tak wie
‘
c ostatecznie otrzymujemy
det A
=
n
∏
i
=1
u
ii
(4.38)
29
Rozdział 5
Interpolacja
5.1
Wielomian interpolacyjny Lagrange’a
Poszukujemy og´olnej zale˙zno´sci funkcyjnej w sytuacji gdy znamy te
‘
zale˙zno´s´c w
(n + 1) punktach:
x
0
x
1
... x
n
y
0
y
1
... y
n
.
(5.1)
Punkty x
i
nazywamy we
‘
złami interpolacji, przy czym x
i
< x
j
dla i
< j. Problem ten
ma jednoznaczne rozwia
‘
zanie je´sli poszukiwana
‘
funkcja
‘
jest wielomian stopnia n
W
n
(x) = a
0
+ a
1
x
+ ... + a
n
x
n
.
(5.2)
Wsp´ołczynniki wielomianu a
i
musza
‘
by´c tak dobrane by wielomian W
n
przechodził
przez ka˙zdy punkt
(x
i
, y
i
):
y
i
= W
n
(x
i
)
i
= 0, 1, . . . , n .
(5.3)
Warunek ten prowadzi do układu liniowego
(n + 1) r´owna´n na (n + 1) niewia-
domych wsp´ołczynnik´ow wielomianu a
i
:
a
0
+ a
1
x
0
+ ... + a
n
x
n
0
= y
0
a
0
+ a
1
x
1
+ ... + a
n
x
n
1
= y
1
...
a
0
+ a
1
x
n
+ ... + a
n
x
n
n
= y
n
,
(5.4)
30
lub zapisuja
‘
c w postaci macierzowej
1
x
0
x
2
0
. . . x
n
0
1
x
1
x
2
1
. . . x
n
1
..
.
..
.
..
.
..
.
1
x
n
x
2
n
. . . x
n
n
a
0
a
1
..
.
a
n
=
y
0
y
1
..
.
y
n
.
(5.5)
Jednoznaczne rozwia
‘
znie powy˙zszego układu r´owna´n istnieje gdy˙z wyznacznik
macierzy gł´ownej (zwany wyznacznikiem Vandermonde’a) jest r´o˙zny od zera
det
=
∏
(i, j) i> j
(x
i
− x
j
) 6= 0.
(5.6)
Istnieje wie
‘
c dokładnie jeden wielomian interpoluja
‘
cy (5.2) przechodza
‘
cy przez
punkty (5.1), zwany wielomianem Lagrange’a.
Aby znale´z´c jawna
‘
posta´c wielomianu Lagrange’a nale˙zy rozwia
‘
za´c układ r´ow-
na´n (5.5). Dzie
‘
ki warunkowi jednoznaczno´sci alternatywna
‘
metoda
‘
jest poszuki-
wanie wielomianu w postaci
W
n
(x) = y
0
L
0
(x) + y
1
L
1
(x) + ... + y
n
L
n
(x) ,
(5.7)
gdzie L
i
(x) sa
‘
wielomianami stopnia n spełniaja
‘
cymi naste
‘
puja
‘
cy warunek
L
i
(x
j
) =
δ
i j
=
1
dla i
= j
0
dla i
6= j .
(5.8)
Spełniony jest wtedy warunek (5.3) nało˙zony na wielomian Lagrange’a:
W
n
(x
j
) =
n
∑
i
=0
y
i
L
i
(x
j
) =
n
∑
i
=0
y
i
δ
i j
= y
j
.
(5.9)
Zgodnie z warunkiem (5.8) ka˙zdy wielomian L
i
(x) zeruje sie
‘
we wszystkich
we
‘
złach z wyja
‘
tkiem x
= x
i
. Jest wie
‘
c proporcjonalny do iloczynu wszystkich
jednomian´ow
(x − x
j
) z wyła
‘
czeniem czynnika
(x − x
i
):
L
i
(x) = a (x − x
0
) ... (x − x
i
−1
) (x − x
i
+1
) ... (x − x
n
) .
(5.10)
Wsp´ołczynnik a wyznaczamy z warunku L
i
(x
i
) = 1, co prowadzi do wzoru
L
i
(x) =
(x − x
0
) ... (x − x
i
−1
) (x − x
i
+1
) ... (x − x
n
)
(x
i
− x
0
) ... (x
i
− x
i
−1
) (x
i
− x
i
+1
) ... (x
i
− x
n
)
(5.11)
31
Wyra˙zenia (5.7) i (5.11) tworza
‘
ostateczna
‘
posta´c wielomianu interpolacyjnego
Lagrange’a.
Przyklad
Rozwa˙zmy interpolacje
‘
zale˙zno´sci funkcyjnej okre´slonej w trzech we
‘
złach
x
0
, x
1
, x
2
. Wielomian Lagrange’a przyjmuje posta´c
W
2
(x) = y
0
(x − x
1
)(x − x
2
)
(x
0
− x
1
)(x
0
− x
2
)
+ y
1
(x − x
0
)(x − x
2
)
(x
1
− x
0
)(x
1
− x
2
)
+ y
2
(x − x
0
)(x − x
1
)
(x
2
− x
0
)(x
2
− x
1
)
Podsumowuja
‘
c, interpolacja przy pomocy wielomianu Lagrange’a jest jednoz-
nacznym rozwia
‘
zaniem problem znalezienia funkcji, przechodza
‘
cej przez zadane
punkty. Poniewa˙z prawdziwa zale˙zno´s´c funkcyjna nie jest znana, pytanie o bła
‘
d
interpolacji nie ma sensu.
5.2
Interpolacja Lagrange’a znanej funkcji
Cze
‘
sto w zastosowaniach numerycznych mamy do czynienia z problemem przy-
bli˙zenia znanej funkcji y
= f (x) w przedziale [a, b] przy pomocy wyra˙zenia wyko-
rzystuja
‘
cego znajomo´s´c tej funkcji w sko´nczonej liczbie punkt´ow z tego przedziału.
Wybierzmy w tym celu
(n + 1) we
‘
zł´ow
a
= x
0
< x
1
< . . . < x
n
= b
(5.12)
i znajd´zmy odpowiadaja
‘
ce im warto´sci funkcji
f
(x
0
)
f
(x
1
)
...
f
(x
n
) .
(5.13)
Przybli˙zmy naste
‘
pnie funkcje
‘
przy pomocy wielomianu Lagrange’a (5.7):
W
n
(x) =
n
∑
i
=0
f
(x
i
) L
i
(x) .
(5.14)
Zakładaja
‘
c, ˙ze istnieje
(n + 1) pochodna f mo˙zna pokaza´c (dow´od w [?]), ˙ze
f
(x) −W
n
(x) =
f
(n+1)
(
ξ
)
(n + 1)!
ω
n
+1
(x)
(5.15)
32
gdzie
ξ
∈ [a,b], natomiast
ω
n
+1
(x) to wielomian stopnia (n + 1):
ω
n
+1
(x) = (x − x
0
)(x − x
1
) . . . (x − x
n
) .
(5.16)
Ze wzoru tego wynika, ˙ze maksymalny bła
‘
d interpolacji Lagrange’a jest ograni-
czony przez
| f (x) −W
n
(x)| ≤
M
n
+1
(n + 1)!
|
ω
n
+1
(x)|
(5.17)
gdzie M
n
+1
to kres g´orny
M
n
+1
= sup
x
∈[a,b]
| f
(n+1)
(x)|
(5.18)
Zauwa˙zmy, ˙ze dla wielomian´ow stopnia n pochodne rze
‘
du
≥ (n + 1) znikaja
‘
. Sta
‘
d
M
n
+1
= 0 i bła
‘
d interpolacji Lagrange’a wynosi zero. Prawdziwe sa
‘
zatem
Twierdzenia
Dla ka˙zdego wielomianu stopnia n, interpolacja Lagrange’a z liczba
‘
we
‘
zł´ow
≥ (n + 1) jest dokładna.
Interpolacja Lagrange’a przy pomocy
(n + 1) we
‘
zł´ow jest dokładna dla ka˙z-
dego wielomianu stopnia
≤ n.
Powstaje pytanie, jak dobra´c we
‘
zły interpolacji by prawa strona wyra˙zenia
(5.17) była jak najmniejsza. Nale˙zy w tym celu znale´z´c najmniejsze oszacowanie
kresu g´ornego warto´sci
|
ω
(n+1)
(x)| w przedziale [a,b]. Odpowied´z na to pytanie
otrzymujemy rozwa˙zaja
‘
c wielomiany Czebyszewa.
5.3
Wielomiany Czebyszewa
Wielomian Czebyszewa stopnia n jest zdefiniowany przy pomocy wzoru
T
n
(x) = cos(n arccos x)
(5.19)
dla n
= 0, 1, 2, . . . oraz x ∈ [−1,1]. Z definicji tej wynika wa˙zna relacja
|T
n
(x)| ≤ 1.
(5.20)
33
Dwa pierwsze wielomiany Czebyszewa to
T
0
(x) = 1
T
1
(x) = x .
(5.21)
Kolejne wielomiany mo˙zna wyliczy´c korzystaja
‘
c z relacji rekurencyjnej słusznej
dla warto´sci stopnia wielomianu n
≥ 1:
T
n
+1
(x) = 2x T
n
(x) − T
n
−1
(x) ,
(5.22)
Dow ´od
Podstawiaja
‘
c
φ
= arccos x, rozwa˙zmy
T
n
+1
(x) = cos(n + 1)
φ
= cos(n
φ
) cos
φ
− sin(n
φ
) sin
φ
= x T
n
(x) − sin(n
φ
) sin
φ
.
(5.23)
Wykorzystuja
‘
c naste
‘
pnie to˙zsamo´s´c trygonometryczna
‘
sin
α
sin
β
=
1
2
{cos(
α
−
β
) − cos(
α
+
β
)}
(5.24)
otrzymamy
T
n
+1
= x T
n
−
1
2
{cos((n − 1)
φ
) − cos((n + 1)
φ
)}
= x T
n
−
1
2
(T
n
−1
− T
n
+1
) .
(5.25)
Sta
‘
d wynika ju˙z relacja (5.22).
Tak wie
‘
c naste
‘
pne wielomiany Czebyszewa to
T
2
(x) = 2x
2
− 1
T
3
(x) = 4x
3
− 3x
T
4
(x) = 8x
4
− 8x
2
+ 1 .
(5.26)
Wida´c sta
‘
d, ˙ze wielomian T
n
jest funkcja
‘
o okre´slonej parzysto´sci, dodatniej dla
parzystych n i ujemna
‘
dla nieparzystych warto´sci stopnia wielomianu. Ponadto,
rozwa˙zaja
‘
c wsp´ołczynnik przy najwy˙zszej pote
‘
dze zauwa˙zamy, ˙ze
T
n
(x) = 2
n
−1
x
n
+ ··· .
(5.27)
34
Ka˙zdy wielomian Czebyszewa stopnia n ma dokładnie n pierwiastk´ow rzeczy-
wistch. Aby je znale´z´c rozwia
‘
zujemy r´ownanie
T
n
(x) = cos(n arccos x) = 0 ,
(5.28)
otrzymuja
‘
c naste
‘
puja
‘
ce pierwiastki dla k
= 0, 1, . . . , n − 1
x
k
= cos
π
(k + 1/2)
n
(5.29)
Dowolny wielomian Czebyszewa mo˙zna wie
‘
c zapisa´c w postaci
T
n
(x) = 2
n
−1
(x − x
0
)(x − x
1
) . . . (x − x
n
−1
) = 2
n
−1
ω
n
(x) ,
(5.30)
gdzie wsp´ołczynnik 2
n
−1
wynika z uwagi (5.27).
5.4
Optymalny wyb´or we¸zł´ow interpolacji
Rozwa˙zmy najpierw funkcje
‘
f
(x) okre´slona
‘
na przedziale
[−1,1]. Oszacujemy
maksymalna
‘
warto´s´c kresu g´ornego warto´sci
|
ω
(n+1)
(x)| w tym przedziale. W tym
celu wybierzmy
(n + 1) we
‘
zł´ow interpolacji be
‘
da
‘
cych zerami wielomianu Czeby-
szewa T
n
+1
(x):
x
k
= cos
π
(k + 1/2)
(n + 1)
k
= 0, 1, . . . n .
(5.31)
Na podstawie relacji (5.30) oraz własno´sci (5.20) otrzymujemy dla tak wybranych
we
‘
zł´ow relacje
‘
|
ω
n
+1
(x)| =
1
2
n
|T
n
+1
(x)| ≤
1
2
n
.
(5.32)
Wz´or (5.17) przyjmuje wie
‘
c teraz posta´c
| f (x) −W
n
(x)| ≤
M
n
+1
2
n
(n + 1)!
(5.33)
gdzie M
n
+1
to kres g´orny warto´sci modułu
(n + 1) pochodnej funkcji f :
M
n
+1
=
sup
x
∈[−1,1]
| f
(n+1)
(x)|.
(5.34)
Otrzymali´smy w ten spos´ob najmniejsze oszacowanie maksymalnego błe
‘
du inter-
polacji Lagrange’a.
35
5.5
Wzory dla dowolnego przedziału
Rozwa˙zmy naste
‘
pnie funkcje
‘
f
(y) okre´slona
‘
w dowolnym przedziale
[a, b]. Op-
tymalnie wybranymi we
‘
złami sa
‘
teraz obrazy zer (5.31) poprzez transformacje
‘
liniowa
‘
y
k
=
1
2
{(b − a)x
k
+ (a + b)} ,
(5.35)
Łatwo sprawdzi´c, wszystkie obrazy we
‘
zł´ow Czebyszewa le˙za
‘
w przedziale
[a, b].
Zastosujemy teraz wzory z poprzedniego rozdziału do y
∈ [a,b]
ω
n
+1
(y) =
n
∏
k
=0
(y − y
k
) =
(b − a)
n
+1
2
n
+1
n
∏
k
=0
(x − x
k
)
=
b
− a
2
n
+1
ω
n
+1
(x) .
(5.36)
Wykorzystuja
‘
c naste
‘
pnie relacje
‘
(5.32) otrzymujemy
|
ω
n
+1
(y)| ≤
b
− a
2
n
+1
1
2
n
.
(5.37)
Sta
‘
d najmniejsze oszacowanie błe
‘
du maksymalnego interpolacji Lagrange’a
| f (y) −W
n
(y)| ≤
M
n
+1
2
n
(n + 1)!
b
− a
2
n
+1
(5.38)
gdzie M
n
+1
to kres g´orny warto´sci modułu
(n + 1) pochodnej funkcji f :
M
n
+1
= sup
y
∈[a,b]
| f
(n+1)
(y)|.
(5.39)
5.6
Interpolacja przy pomocy funkcji sklejanych
Dla ka˙zdego układu
(n + 1) we
‘
zł´ow w przedziale
[a, b] istnieje funkcja cia
‘
gła
w tym przedziale, dla kt´orej metoda interpolacyjna nie jest zbie˙zna. Dla takiej
funkcji zwie
‘
kszenie liczby we
‘
zł´ow pogarsza dokładno´s´c interpolacji (zwłaszcza
bli˙zej ko´nc´ow przedziału). Przykładem jest funkcja f
(x) = |x| interpolowana przy
pomocy r´ownoodległych we
‘
zł´ow.
36
W interpolacji przy pomocy funkcji sklejanych mo˙zna osia
‘
gna
‘
´c bardzo dobra
‘
dokładno´s´c przy pomocy małej liczby we
‘
zł´ow. Rozwa˙zmy
(n + 1) we
‘
zł´ow oraz n
odpowiadaja
‘
cych im przedział´ow:
[x
0
, x
1
] ∪ [x
1
, x
2
] ∪ ... ∪ [x
n
−1
, x
n
] ,
(5.40)
numerowanych indeksem prawego ko´nca przedziału.
W ka˙zdym przedziale funkcja f
(x) jest interpolowana przy pomocy wielomia-
nu trzeciego stopnia:
P
i
(x) = a
0 i
+ a
1 i
x
+ a
2 i
x
2
+ a
3 i
x
3
,
(5.41)
gdzie i
= 1, 2, . . . n. Otrzymujemy wie
‘
c n wielomian´ow o 4n wsp´ołczynnikach do
wyznaczenia. Nakładamy na nie naste
‘
puja
‘
ce 4n warunk´ow, kt´ore prowadza
‘
do
jednoznacznego rozwia
‘
zania.
1. Na brzegach ka˙zdego przedziału spełnione jest 2n warunk´ow interpolacji:
P
i
(x
i
−1
) = f (x
i
−1
)
P
i
(x
i
) = f (x
i
) ,
(5.42)
gdzie i
= 1, 2, . . . n.
2. W punktach wewne
‘
trznych
{x
1
, . . . , x
n
−1
} spełnione jest 2(n − 1) warunk´ow
cia
‘
gło´sci dla pierwszych i drugich pochodnych:
P
′
i
(x
i
) = P
′
i
+1
(x
i
)
P
′′
i
(x
i
) = P
′′
i
+1
(x
i
) ,
(5.43)
gdzie i
= 1, 2, . . . (n − 1).
3. W punktach zewne
‘
trznych
{x
0
, x
n
} ˙za
‘
damy dla drugich pochodnych:
P
′′
1
(x
0
) = f
′′
(x
0
)
P
′′
n
(x
n
) = f
′′
(x
n
) .
(5.44)
Sta
‘
d dwa warunki.
W sumie otrzymujemy
{2n + 2(n − 1) + 2} = 4n poszukiwanych warunk´ow.
37
Rozdział 6
Aproksymacja
W wielu przypadkach przybli˙zanie zale˙zno´sci funkcyjnej poprzez wielomian prze-
chodza
‘
cy przez wszystkie znane punkty nie jest dobra
‘
metoda
‘
, w szczeg´olno´sci,
gdy punkty sa
‘
obarczone błe
‘
dem lub gdy jest ich bardzo du˙zo. W tym ostat-
nim przypadku stopie´n wielomianu Lagrange’a byłby bardzo wysoki co prowadzi´c
mo˙ze do niestabilno´sci numerycznych.
Znacznie lepsza
‘
metoda
‘
jest wtedy aproksymacja przy pomocy wzgle
‘
dnie pro-
stej funkcji, tak by była ona jak najmniej “odległa” od aproksymowanej funkcji.
Metode
‘
te
‘
zilustrujemy na pocza
‘
tek w najprostszym przypadku, gdy dobrym przy-
bli˙zeniem jest zało˙zenie o liniowej zale˙zno´sci funkcyjnej.
6.1
Regresja liniowa
Zał´o˙zmy, ˙ze mamy do czynienia z n we
‘
złami
{x
0
, x
1
. . . x
n
−1
} i odpowiadaja
‘
cymi
im warto´sciami
{y
0
, y
1
. . . y
n
−1
}. Zał´o˙zmy, ˙ze dobra
‘
poszukiwana
‘
zale˙zno´scia
‘
jest
funkcja liniowa
y
= ax + b ,
(6.1)
gdzie parametry a i b pozostaja
‘
do wyznaczenia. Utw ´orzmy w tym celu funkcje
‘
S
(a, b) =
n
−1
∑
k
=0
{y
k
− (ax
k
+ b)}
2
(6.2)
38
be
‘
da
‘
ca
‘
suma
‘
kwadrat´ow odległo´sci punkt´ow
(x
k
, y
k
) od prostej (6.1), mierzonych
wzdłu˙z osi y. Dobierzmy naste
‘
pnie parametry a i b tak, by warto´s´c tej funkcji była
jak najmniejsza. R ´o˙zniczkuja
‘
c, otrzymamy
∂
S
∂
a
=
n
−1
∑
k
=0
(−2x
k
)(y
k
− ax
k
− b) = 0
∂
S
∂
b
=
n
−1
∑
k
=0
(−2)(y
k
− ax
k
− b) = 0,
(6.3)
Sta
‘
d dostajemy układ r´owna´n na wsp´ołczynniki a i b:
a
(
∑
k
x
2
k
) + b (
∑
k
x
k
) =
∑
k
x
k
y
k
a
(
∑
k
x
k
) + b (
∑
k
1
) =
∑
k
y
k
.
(6.4)
lub w postaci macierzowej
∑
k
x
2
k
∑
k
x
k
∑
k
x
k
n
a
b
=
∑
k
x
k
y
k
∑
k
y
k
.
(6.5)
Wykorzystali´smy przy tym obserwacje
‘
∑
n
−1
k
=0
1
= n. Rozwia
‘
zanie zadane jest poprzez
wzory Cramera
a
=
n
(
∑
k
x
k
y
k
) − (
∑
k
x
k
)(
∑
k
y
k
)
n
(
∑
k
x
2
k
) − (
∑
k
x
k
)
2
b
=
(
∑
k
x
2
k
)(
∑
k
y
k
) − (
∑
k
x
k
y
k
)(
∑
k
x
k
)
n
(
∑
k
x
2
k
) − (
∑
k
x
k
)
2
.
(6.6)
Uog´olnienie tej metody polega na rozwa˙zeniu zale˙zno´sci wielomianowej
y
= a
0
+ a
1
x
+ . . . + a
M
x
M
(6.7)
i minimalizacje
‘
ze wzgle
‘
du na wsp´ołczynniki tego wielomianu naste
‘
puja
‘
cej r´o˙znicy
S
(a
i
) =
n
−1
∑
k
=0
y
k
− (a
0
+ a
1
x
+ . . . + a
M
x
M
)
2
.
(6.8)
Zwr´o´cmy uwage
‘
, ˙ze na og´oł liczba punkt´ow n, w kt´orych znamy aproksymowana
‘
funkcje
‘
, jest du˙zo wieksza od stopnia wielomianu aproksymuja
‘
cego M.
39
6.2
Aproksymacja ´sredniokwadratowa
Przedstawiona w poprzednim rozdziale metoda jest przykładem aproksymacji ´sre-
dniokwadratowej. W og´olno´sci, znaja
‘
c n we
‘
zł´ow
{x
0
, x
1
. . . x
n
−1
} oraz odpowiada-
ja
‘
cych im warto´sci
{ f (x
0
), f (x
1
), . . . f (x
n
−1
)}, nieznanej na og´oł funkcji, chcemy
ja
‘
aproksymowa´c przy pomocy wyra˙zenia
F
(x) = c
0
φ
0
(x) + c
1
φ
1
(x) + . . . + c
M
φ
M
(x) ,
(6.9)
gdzie
{
φ
0
(x),
φ
1
(x), . . .
φ
M
(x)}
(6.10)
jest układem dogodnie wybranych funkcji bazowych. W przykładzie (6.8) funk-
cjami bazowymi sa
‘
jednomiany
φ
I
(x) = x
I
. Zauwa˙zmy, ˙ze liczba we
‘
zł´ow n jest
niezale˙zna od liczby funkcji bazowych
(M + 1).
Wsp´ołczynniki c
i
sa
‘
tak dobrane by zminimalizowa´c odległo´s´c pomie
‘
dzy funk-
cja
‘
dokładna
‘
f
(x) i przybli˙zona
‘
F
(x), zdefiniowana
‘
poprzez:
S
(c
i
) =
n
−1
∑
k
=0
{ f (x
k
) − (c
0
φ
0
(x
k
) + c
1
φ
1
(x
k
) + . . . + c
M
φ
M
(x
k
))}
2
.
(6.11)
R ´o˙zniczkuja
‘
c po parametrach otrzymujemy dla ka˙zdego parametru c
i
:
∂
S
∂
c
i
=
n
−1
∑
k
=0
(−2)
φ
i
(x
k
)
f
(x
k
) − c
0
φ
0
(x
k
) − c
1
φ
1
(x
k
) − ... − c
M
φ
M
(x
k
)
= 0 .
Sta
‘
d naste
‘
puja
‘
cy układ r´owna´n na poszukiwane wsp´ołczynniki
(
φ
0
,
φ
0
) c
0
+ (
φ
0
,
φ
1
) c
1
+ ··· + (
φ
0
,
φ
M
) c
M
= (
φ
0
, f )
(
φ
1
,
φ
0
) c
0
+ (
φ
1
,
φ
1
) c
1
+ ··· + (
φ
1
,
φ
M
) c
M
= (
φ
1
, f )
...
(
φ
n
,
φ
0
) c
0
+ (
φ
n
,
φ
1
) c
1
+ ··· + (
φ
n
,
φ
M
) c
M
= (
φ
n
, f ) ,
(6.12)
gdzie
(·, ·) to peudo-iloczyn skalarny
1
, zdefiniowany przy pomocy dyskretnego
zbioru we
‘
zł´ow:
(
φ
I
,
φ
J
) =
n
−1
∑
k
=0
φ
I
(x
k
)
φ
J
(x
k
)
(6.13)
1
Pseudo-, gdy˙z nie jest spełniona własno´s´c iloczynu skalarnego
(
φ
I
,
φ
J
) = 0 =>
φ
I
(x) = 0. W
takim wypadku funkcje
φ
I
(x) musza
‘
sie
‘
zerowa´c tylko w we
‘
złach x
k
.
40
Rozwia
‘
zanie układu (6.12) jest szczeg´olnie proste, gdy układ funkcji bazowych
jest ortogonalny:
(
φ
I
,
φ
J
) ∼
δ
IJ
.
(6.14)
W r´ownaniach (6.12) pozostaja
‘
tylko składowe diagonalne i sta
‘
d
c
I
=
(
φ
I
, f )
(
φ
I
,
φ
I
)
,
I
= 0, 1, . . . , M .
(6.15)
Ostatecznie funkcja aproksymuja
‘
ca to
F
(x) =
M
∑
I
=0
(
φ
I
, f )
(
φ
I
,
φ
I
)
φ
I
(x)
(6.16)
W naste
‘
pnych rozdziałach przedstawimy przykłady dw ´och szczeg´olnie wa˙znych
realizacji powy˙zszego rozwia
‘
zania. Rozwa˙zymy aproksymacje
‘
przy pomocy wielo-
mian´ow Czebyszewa oraz funkcji trygonometrycznych.
6.3
Aproksymacja Czebyszewa
Wybierzmy n we
‘
zł´ow
{x
0
, x
1
. . . x
n
−1
} be
‘
da
‘
cych zerami wielomianu Czebyszewa:
T
n
(x
k
) = 0
(6.17)
Zgodnie ze wzorem (5.29) we
‘
zły to
x
k
= cos
π
(2k + 1)
2n
,
k
= 0, 1, . . . (n − 1).
(6.18)
Wielomiany Chebyszewa
{T
0
, T
1
, . . . , T
n
−1
} sa
‘
ortogonalne wzgle
‘
dem peudo-ilo-
czynu skalarnego (6.13) z powy˙zszymi we
‘
złami.
(T
I
, T
J
) =
n
−1
∑
k
=0
T
I
(x
k
) T
J
(x
k
) =
n
dla
I
= J = 0
n
/2
dla
I
= J 6= 0
0
dla
I
6= J
(6.19)
41
Dow ´od
Rozwa˙zmy peudo-iloczyn skalarny dla 0
≤ I,J ≤ (n − 1):
n
−1
∑
k
=0
T
I
(x
k
) T
J
(x
k
) =
n
−1
∑
k
=0
cos
I
π
(2k + 1)
2n
cos
J
π
(2k + 1)
2n
.
(6.20)
Dla I
= J = 0 otrzymujemy sume
‘
r´owna
‘
n.
Rozwa˙zmy naste
‘
pnie przypadek I
6= J. Korzystaja
‘
c ze wzoru
cos
α
cos
β
=
1
2
{cos(
α
+
β
) + cos(
α
−
β
)}
(6.21)
otrzymujemy
(T
I
, T
J
) =
1
2
n
−1
∑
k
=0
cos
(2k + 1)
π
(I − J)
2n
+ cos
(2k + 1)
π
(I + J)
2n
.
(6.22)
Policzmy naste
‘
pnie dla dowolnego ka
‘
ta
φ
6= 0:
n
−1
∑
k
=0
cos
(2k + 1)
φ
= Re
n
−1
∑
k
=0
e
i
(2k+1)
φ
= Re
n
e
i
φ
(1 + e
2i
φ
+ . . . + e
2i
(n−1)
φ
)
o
= Re
e
i
φ
1
− e
2i n
φ
1
− e
2i
φ
=
sin
(2n
φ
)
2 sin
φ
.
(6.23)
Podstawiaja
‘
c
φ
±
=
π
(I ± J)/2n 6= 0, otrzymujemy dla wzoru (6.22)
(T
I
, T
J
) =
sin
π
(I − J)
4 sin
φ
+
+
sin
π
(I + J)
4 sin
φ
−
.
(6.24)
Wyra˙zenie to znika dla I
6= J ze wzgle
‘
du na zerowanie sie
‘
sinus´ow w liczniku.
Sta
‘
d warunek ortogonalno´sci dla wielomian´ow Czebyszewa.
Dla I
= J 6= 0 r´ownanie (6.22) przyjmuje posta´c
(T
I
, T
I
) =
1
2
n
−1
∑
k
=0
1
+ cos(2k + 1)
π
I
n
=
n
2
(6.25)
gdy˙z z r´ownania (6.23) dla
φ
=
π
I
/n wynika, ˙ze suma cosinus´ow daje zero.
42
Mo˙zemy wie
‘
c aproksymowa´c przy pomocy wielomian´ow Czebyszewa funkcje
‘
f
(x) okre´slona
‘
w przedziale
[−1,1]:
f
(x) ≃
1
2
c
0
+
n
−1
∑
J
=0
c
J
T
J
(x) .
(6.26)
Wsp´ołczynniki c
J
sa
‘
wyliczone ze wzoru (6.15), w kt´orym
φ
J
= T
J
:
c
J
=
2
n
n
−1
∑
k
=0
T
J
(x
k
) f (x
k
)
(6.27)
gdzie dla J
= 0 podstawili´smy dwukrotnie wie
‘
ksza
‘
warto´s´c c
0
ni˙z ta wynikaja
‘
ca z
normalizacji:
(T
0
, T
0
) = n. Sta
‘
d wsp´olczynnik 1
/2 we wzorze (6.26).
W sumie (6.26) mo˙zna zachowa´c jedynie M
≤ (n − 1) składnik´ow przy nie-
zmienionych wsp´ołczynnikach c
J
. W wie
‘
kszo´sci przypadk´ow staja
‘
sie
‘
one coraz
mniejsze dla rosna
‘
cych warto´sci wska´znika, natomiast ka˙zdy wielomian Czeby-
szewa jest ograniczony warunkiem (5.20). Nie popełniamy w ten spos´ob du˙zego
błe
‘
du odrzucaja
‘
c wyrazy z J
> M. Sta
‘
d ostateczny wz´or
f
(x) ≃
1
2
c
0
+
M
<n
∑
J
=0
c
J
T
J
(x)
(6.28)
Podsumowuja
‘
c, kluczowym punktem w aproksymacji Czebyszewa jest zna-
jomo´s´c funkcji f w we
‘
złach Chebyszewa. Na tej podstawie konstruuje sie
‘
wsp´oł-
czynniki (6.27), a naste
‘
pnie przybli˙zenie (6.28).
6.4
Aproksymacja Czebyszewa w dowolnym przedziale
Rozwa˙zmy aproksymacje
‘
funkcji f
(y) okre´slona
‘
dla y
∈ [a,b]. Niech
y
= y(x) ,
x
∈ [−1,1]
(6.29)
be
‘
dzie dowolna
‘
transformacja
‘
bijektywna
‘
odwzorowuja
‘
ca
‘
[−1,1] → [a,b]. Trans-
formacja odwrotna to
x
= y
−1
(y) ,
y
∈ [a,b].
(6.30)
43
Przyklad
Dla transformacji liniowej
y
= y(x) =
1
2
{(b − a)x + (a + b)}
x
∈ [−1,1],
(6.31)
transformacja odwrotna to
x
= y
−1
(x) =
2y
− a − b
b
− a
y
∈ [a,b].
(6.32)
Zdefiniujmy nowa
‘
funkcje
‘
okre´slona
‘
dla x
∈ [−1,1] wzorem:
˜f(x) = f (y(x)),
(6.33)
a naste
‘
pnie zastosujmy do niej wz´or aproksymacyjny (6.28):
˜
f
(x) ≈
1
2
c
0
+
M
<n
∑
J
=0
c
J
T
J
(x) .
(6.34)
Podstawiaja
‘
c po prawej stronie relacje
‘
odwrotna
‘
x
= y
−1
(y), otrzymujemy apro-
ksymacje
‘
funkcji f
(y):
f
(y) ≈
1
2
c
0
+
M
<n
∑
J
=0
c
J
T
J
(y
−1
(y))
(6.35)
gdzie wsp´ołczynniki c
J
dla J
= 0, 1 . . . (n − 1) sa
‘
zadane przez
c
J
=
2
n
n
−1
∑
k
=0
T
J
(x
k
) f (y
k
)
(6.36)
gdzie y
k
= y(x
k
).
6.5
Aproksymacja trygonometryczna
Aproksymacje
‘
te
‘
stosujemy, gdy mamy do czynienia z funkcja
‘
okresowa
‘
f
(x) o
okresie 2
π
. Zał´o˙zmy, ˙ze znamy te
‘
funkcje
‘
w parzystej liczbie 2n r´ownoodległych
punkt´ow z przedziału
[0, 2
π
]:
x
k
= k
π
n
,
k
= 0, 1, . . . , (2n − 1)
(6.37)
44
Wła´sciwym układem funkcji bazowych dla aproksymacji (6.9) sa
‘
wtedy funkcje
trygonometryczne
{1, sin x, cos x, sin 2x, cos 2x ... sin(n − 1)x, cos(n − 1)x}
(6.38)
Zbi´or tych funkcji jest ortogonalny wzgle
‘
dem peudo-iloczynu skalarnego (6.13) z
punktami (6.37).
Dow ´od
Zachodzi bowiem dla 1
≤ I,J ≤ (n − 1):
2n
−1
∑
k
=0
e
iIx
k
=
2n
−1
∑
k
=0
{cos(Ix
k
) + i sin(Ix
k
)}
=
2n
−1
∑
k
=0
e
(iI
π
/n)k
=
1
− e
i2
π
I
1
− e
iI
π
/n
= 0.
(6.39)
Sta
‘
d relacje ortogonalno´sci
2n
−1
∑
k
=0
sin
(Ix
k
) · 1 = 0
2n
−1
∑
k
=0
cos
(Ix
k
) · 1 = 0
2n
−1
∑
k
=0
1
· 1 = 2n.
Ponadto
2n
−1
∑
k
=0
sin
(Ix
k
) · sin(Jx
k
) =
1
2
2n
−1
∑
k
=0
{cos(I − J)x
k
− cos(I + J)x
k
} = n
δ
IJ
2n
−1
∑
k
=0
cos
(Ix
k
) · cos(Jx
k
) =
1
2
2n
−1
∑
k
=0
{cos(I − J)x
k
+ cos(I + J)x
k
} = n
δ
IJ
2n
−1
∑
k
=0
cos
(Ix
k
) · sin(Jx
k
) =
1
2
2n
−1
∑
k
=0
{sin(I − J)x
k
+ sin(I + J)x
k
} = 0.
45
Sta
‘
d, aproksymuja
‘
c funkcje
‘
f
(x) przy pomocy funkcji trygonometrycznych,
otrzymujemy
f
(x) ≈
1
2
a
0
+
M
<n
∑
J
=1
{a
J
cos
(Jx) + b
J
sin
(Jx)} .
(6.40)
Wsp´ołczynniki a
J
i b
J
dla J
= 0, 1 . . . (n − 1) mo˙zna wyliczy´c ze wzoru (6.15):
a
J
=
1
n
2n
−1
∑
k
=0
f
(x
k
) cos(Jx
k
)
b
J
=
1
n
2n
−1
∑
k
=0
f
(x
k
) sin(Jx
k
) .
(6.41)
6.6
Wzory dla dowolnego okresu
Zastosujmy powy˙zsza
‘
aproksymacje
‘
do funkcji czasu o okresie T
f
(t) = f (t + T ) .
(6.42)
Zdefiniujmy nowa
‘
funkcje
‘
˜
f
(x) o okresie 2
π
:
˜f(x) = f (t)
(6.43)
przy pomocy transformacji:
t
= x
T
2
π
.
(6.44)
Aproksymujemy funkcje
‘
˜
f
(x) za pomoca
‘
wzoru (6.40), a naste
‘
pnie podstawiamy
po prawej stronie transformacje
‘
odwrotna
‘
:
x
=
2
π
T
t
.
(6.45)
Otrzymujemy w ten spos´ob dla funkcji f
(t):
f
(t) ≈
1
2
a
0
+
M
<n
∑
J
=1
a
J
cos
2
π
J
T
t
+ b
J
sin
2
π
J
T
t
(6.46)
46
Wsp´ołczynniki a
J
i b
J
sa
‘
teraz zadane wzorami:
a
J
=
1
n
2n
−1
∑
k
=0
f
(t
k
) cos(J x
k
)
b
J
=
1
n
2n
−1
∑
k
=0
f
(t
k
) sin(J x
k
) ,
(6.47)
w kt´orych wielko´sci
t
k
= x
k
T
2
π
= k
T
2n
(6.48)
dla k
= 0, 1, . . . (2n − 1), sa
‘
obrazami we
‘
zł´ow (6.37) poprzez transformacje
‘
(6.44).
47
Rozdział 7
R´o˙zniczkowanie
7.1
Metoda z aproksymacja¸
W metodzie tej najpierw aproksymujemy funkcje
‘
, kt´orej pochodna
‘
chcemy znale´z´c
przy pomocy jednej z metod opisanych w poprzednim rozdziale:
f
(x) ≃
M
∑
i
=0
c
i
φ
i
(x) ,
(7.1)
gdzie
φ
i
to znane i r´o˙zniczkowalne funkcje bazowe, np. jednomiany x
i
lub wielo-
miany Czebyszewa T
i
(x). Przyjmujemy, ˙ze pochodna tej funkcji to
f
′
(x) ≃
M
∑
i
=0
c
i
φ
′
i
(x)
(7.2)
Podobnie poste
‘
pujemy przy obliczniu wy˙zszych pochodnych.
7.2
Metody z rozwinie¸ciem Taylora
Ta metoda r´o˙zniczkowania wykorzystuja
‘
rozwinie
‘
cie Taylora funkcji. Zakladaja
‘
c,
˙ze rozwinie
‘
cie takie istnieje w otoczeniu punktu x, mamy
f
(x + h) = f (x) + f
′
(x) h + f
′′
(x)
h
2
2!
+ f
′′′
(x)
h
3
3!
+
O
(h
4
)
(7.3)
48
gdzie symbol
O
(h
4
) oznacza reszte
‘
rze
‘
du h
4
, tzn
lim
h
→0
O
(h
4
)
h
4
= const .
(7.4)
Aby obliczy´c pierwsza
‘
pochodna
‘
wykorzystujemy wz´or zapisany z dokładno-
´scia do członu liniowego w h:
f
(x + h) = f (x) + f
′
(x) h +
O
(h
2
) .,
Sta
‘
d wynika
f
′
(x) =
f
(x + h) − f (x)
h
+
O
(h) .
(7.5)
Aby poprawi´c dokładno´s´c oblicznej w ten spos´ob pochodnej wykorzystujemy dwa
rozwinie
‘
cia Taylora zapisane z dokładno´scia
‘
h
2
:
f
(x + h) = f (x) + f
′
(x) h + f
′′
(x)
h
2
2!
+
O
(h
3
)
f
(x − h) = f (x) − f
′
(x) h + f
′′
(x)
h
2
2!
−
O
(h
3
) .
(7.6)
Po odje
‘
ciu stronami wyra˙zenia z parzystymi potegami h upraszczaja
‘
sie
‘
i sta
‘
d
dostajemy
f
′
(x) =
f
(x + h) − f (x − h)
2h
+
O
(h
2
)
(7.7)
Znaja
‘
c wie
‘
c warto´sc funkcji w punktach sa
‘
siednich
(x − h) oraz (x + h), poprawia-
my dokładno´s´c pierwszej pochodnej.
Dodaja
‘
c stronami wyra˙zenia (7.6) otrzymujemy wz´or na druga
‘
pochodna
‘
f
′′
(x) =
f
(x + h) − 2 f (x) + f (x − h)
h
2
+
O
(h
2
)
(7.8)
Rza
‘
d reszty wynika z faktu kasowanie sie
‘
wyra˙ze´n z nieparzystymi pote
‘
gami h
przy dodawaniu stronami.
7.3
Wie¸ksza dokładno´s´c
Lepsza
‘
dokładno´s´c obliczanych pochodnych mo˙zna uzyska´c rozwa˙zaja
‘
c rozwi-
nie
‘
cia Taylora zapisane z dokładno´scia do pia
‘
tej pote
‘
gi h dla warto´sci funkcji
49
w czterech punktach: f
(x − 2h), f (x − h), f (x + h), f (x + 2h). Jako po˙zyteczne
´cwiczenie nale˙zy udowodni´c, ˙ze w takim przypadku dla pierwszej pochodnej otrzy-
mujemy
f
′
(x) =
f
(x − 2h) − 8 f (x − h) + 8 f (x + h) − f (x + 2h)
12h
+
O
(h
4
)
(7.9)
Natomiast druga pochodna jest dana przez wyra˙zenie
f
′′
(x) =
− f (x − 2h) + 16 f (x − h) − 30 f (x) + 16 f (x + h) − f (x + 2h)
12h
2
+
O
(h
4
)
(7.10)
7.4
Wy˙zsze pochodne
W podobny spos´ob mo˙zna wyprowadzi´c wzory na wy˙zsze pochodne, korzystaja
‘
c
z wyprowadzonych wcze´sniej formuł dla ni˙zszych pochodnych. Na przykład, aby
obliczy´c trzecia
‘
pochodna
‘
odejmujemy od siebie dwa rozwinie
‘
cia
f
(x + h) = f (x) + f
′
(x) h + f
′′
(x)
h
2
2!
+ f
′′′
(x)
h
3
3!
+
O
(h
4
)
f
(x − h) = f (x) − f
′
(x) h + f
′′
(x)
h
2
2!
− f
′′′
(x)
h
3
3!
+
O
(h
4
) ,
(7.11)
otrzymuja
‘
c
f
′′′
(x) = 3
f
(x + h) − f (x − h)
h
3
− 6
f
′
(x)
h
2
+
O
(h
2
) .
(7.12)
Zwr´o´cmy uwage
‘
, ˙ze podstawienie wzoru (7.7) w miejsce pierwszej pochodnej w
powy˙zszym wzorze daje reszte
‘
O
(1), niezale˙zna
‘
od odchylenia h. Jeste´smy wie
‘
c
zmuszeni do u˙zycia dokładniejszego wzoru (7.9), prowadza
‘
cego do
f
′′′
(x) =
− f (x − 2h) + 2 f (x − h) − 2 f (x + h) + f (x + 2h)
2h
3
+
O
(h
2
)
(7.13)
50
Rozdział 8
Całkowanie
8.1
Kwadratury
Rozwa˙zmy jednowymiarowa
‘
całke
‘
na przedziale
[a, b ]
I
( f ) =
b
Z
a
f
(x) dx
(8.1)
Podzielmy przedział całkowania na N r´ownych odcink´ow o długo´sci
h
=
b
− a
N
(8.2)
i wyznaczonych przez kolejne N
+ 1 punkt´ow
x
k
= a + kh
k
= 0, 1, . . . , N .
(8.3)
Zauwa˙zmy, ˙ze x
0
= a oraz x
N
= b sa
‘
ko´ncami przedział´ow. Wtedy
I
( f ) =
N
∑
i
=1
x
i
Z
x
i
−1
f
(x) dx .
(8.4)
Nale˙zy wie
‘
c przybli˙zy´c całke
‘
w ka˙zdym z podprzedział´ow
[x
i
−1
, x
i
]. Otrzymany
wynik nazywa sie
‘
kwadratura
‘
i ma og´olna
‘
posta´c:
I
( f ) ≃
N
∑
i
=0
A
i
f
(x
i
)
(8.5)
A
i
to wsp´ołczynniki (wagi) kwadratury, natomiast punkty x
i
to jej we
‘
zły.
51
8.2
Metoda prostokat´ow
W metodzie prostoka
‘
t´ow
x
i
Z
x
i
−1
f
(x) dx ≃ h f (x
i
)
(8.6)
i wtedy, wprowadza
‘
c oznacznie f
i
≡ f (x
i
),
b
Z
a
f
(x) dx ≃ h( f
1
+ f
1
+ . . . + f
N
)
(8.7)
Otrzymany wynik odpowiada przybli˙zeniu funkcji f w ka˙zdym z podprzedzia-
łow poprzez funkcje
‘
stała
‘
be
‘
a
‘
ca
‘
warto´scia
‘
funkcji w prawym ko´ncu ka˙zdego pod-
przedziału.
Mo˙zemy r´ownie˙z wykorzysta´c warto´sci funkcji w lewych ko´ncach i wtedy
b
Z
a
f
(x) dx ≃ h( f
0
+ f
1
+ . . . + f
N
−1
)
(8.8)
Warto´sci funkcji f
i
moga
‘
by´c r´ownie˙z wzie
‘
te w ´srodkach podprzedział´ow
f
i
= f
x
i
+ x
i
+1
2
.
(8.9)
8.3
Metoda trapez´ow
W metodzie trapez´ow otrzymujemy
x
i
Z
x
i
−1
f
(x) dx ≃
h
2
{ f
i
−1
+ f
i
}
(8.10)
i wtedy
b
Z
a
f
(x) dx ≃
h
2
f
0
+ 2( f
1
+ . . . + f
N
−1
) + f
N
(8.11)
52
Prawa strona wzoru (8.10) to pole trapezu gdy obie warto´sci funkcji f
i
−1
i f
i
sa
‘
dodatnie, lub minus pole, gdy obie sa
‘
ujemne. W przypadku, gdy warto´sci r´o˙znia
‘
sie
‘
znakiem przedstawiona interpretacja geometryczna nie jest ju˙z prawdziwa. Mo-
˙zna jednak sformułowa´c problem og´olnie korzystaja
‘
c z liniowej interpolacji La-
grange’a w ka˙zdym z podprzedziałow
(x
i
−1
, x
i
):
f
(x) ≃ f
i
−1
x
− x
i
x
i
−1
− x
i
+ f
i
x
− x
i
−1
x
i
− x
i
−1
.
(8.12)
Całkuja
‘
c bowiem (8.12) i pamie
‘
taja
‘
c, ˙ze x
i
− x
i
−1
= h, otrzymujemy wz´or (8.10)
niezale˙znie od znaku warto´sci funkcji w kra´ncach podprzedział´ow:
1
h
x
i
Z
x
i
−1
{ f
i
(x − x
i
−1
) − f
i
−1
(x − x
i
)} dx =
h
2
{ f
i
−1
+ f
i
}.
Łatwo uog´olni´c te
‘
metode
‘
, przybli˙zaja
‘
c funkcje
‘
podcałkowa
‘
przy pomocy wie
‘
kszej
liczby punkt´ow. Przykładem jest metoda Simpsona.
8.4
Metoda parabol Simpsona
W metodzie tej dzielimy przedział
[a, b ] na parzysta
‘
liczbe
‘
2N przedziałow o
r´ownej długo´sci:
h
=
b
− a
2N
.
(8.13)
W ka˙zdej sa
‘
siedniej parze przedział´ow wyznaczonych przez
(x
2i
−2
, x
2i
−1
, x
2i
) sto-
sujemy interpolacje
‘
Lagrange’a:
f
(x) ≃ f
2i
−2
(x − x
2i
−1
)(x − x
2i
)
(x
2i
−2
− x
2i
−1
)(x
2i
−2
− x
2i
)
+
f
2i
−1
(x − x
2i
−2
)(x − x
2i
)
(x
2i
−1
− x
2i
−2
)(x
2i
−1
− x
2i
)
+
f
2i
(x − x
2i
−2
)(x − x
2i
−1
)
(x
2i
− x
2i
−2
)(x
2i
− x
2i
−1
)
.
(8.14)
Pamie
‘
taja
‘
c, ˙ze
(x
2i
− x
2i
−1
) = (x
2i
−1
− x
2i
−2
) = h
(8.15)
53
otrzymujemy po wykonaniu całkowania
x
2i
Z
x
2i
−2
f
(x) dx ≃
h
3
{ f
2i
−2
+ 4 f
2i
−1
+ f
2i
}.
(8.16)
Sta
‘
d przybli˙zony wz´or dla warto´sci całki:
b
Z
a
f
(x) dx ≃
h
3
f
0
+ 4( f
1
+ . . . + f
2N
−1
) + 2( f
2
+ . . . + f
2N
−2
) + f
2N
(8.17)
8.5
Bła¸d przybli˙ze ´n
Zdefiniujmy bła
‘
d obliczenia całki, rozumiany jako r´o˙znice
‘
pomie
‘
dzy warto´scia
‘
dokładna
‘
a przybli˙zona
‘
ε
≡ I( f ) −
N
∑
i
=0
A
i
f
(x
i
) .
(8.18)
Mo˙zna pokaza´c (podre
‘
cznik [?]), ˙ze bła
‘
d ten jest ograniczony od g´ory.
W metodzie prostoka
‘
t´ow:
|
ε
| < (b − a)
h
2
sup
x
∈[a,b]
| f
′
(x)|
(8.19)
Całkowanie ta
‘
metoda
‘
jest wie
‘
c dokładne dla funkcji stałej – wielomianu stopnia
zerowego.
W metodzie trapez´ow:
|
ε
| < (b − a)
h
2
12
sup
x
∈[a,b]
| f
(2)
(x)|
(8.20)
Tym razem całkowanie jest dokładne dla wszystkich wielomian´ow stopnia
≤ 1, dla
kt´orych znika druga pochodna.
W metodzie parabol:
|
ε
| < (b − a)
h
4
180
sup
x
∈[a,b]
| f
(4)
(x)|
(8.21)
54
Metoda ta jest dokładna dla ka˙zdego wielomianu stopnia
≤ 3, dla kt´orego znika
czwarta pochodna.
Najwy˙zszy stopie´n wielomianu, dla kt´orego metoda całkowania nie jest do-
kładna nazywamy rze
‘
dem metody. Tak wie
‘
c rza
‘
d przedstawionych kwadratur
wynosi odpowiednio: 1
, 2, 4. Poje
‘
cie to odgrywa podstawowa
‘
role
‘
przy konstrukcji
kwadratur Gaussa.
Zauwa˙zmy na koniec, ˙ze przy tej samej długo´sci podprzedział´ow h
< 1, naj-
lepsza parametrycznie jest metoda parabol Simpsona.
55
Rozdział 9
Kwadratury Gaussa
Naszym celem jest znalezienie og´olnej metody przybli˙zonego obliczania całek
I
( f ) =
b
Z
a
f
(x) w(x) dx
(9.1)
Granice całkowania moga
‘
by´c r´owne
±
∞
. Dodatnia funkcja w
(x) ≥ 0 jest nazy-
wana waga
‘
. Mo˙ze ona zawiera´c całkowalne osobliwo´sci, kt´ore wyła
‘
czyli´smy z
funkcji podcałkowej. Na przykład,
w
(x) =
1
√
1
− x
2
,
(9.2)
dla całek w przedziale
[−1,1] lub
w
(x) = e
−x
2
.
(9.3)
dla całek niewła´sciwych z granicami całkowania
±
∞
. Ostatnia waga zapewnia
zbie˙zno´s´c całki dla szerokiej klasy funkcji f .
Wybierzmy N we
‘
zł´ow w przedziale całkowania:
{x
0
, x
1
, . . . x
N
−1
}.
(9.4)
Posłu˙za
‘
one do skonstruowania interpolacji Lagrange’a funkcji f :
f
(x) ≃
N
−1
∑
k
=0
f
(x
k
) L
k
(x) ,
L
k
(x) =
N
−1
∏
i
=0
x
− x
i
x
k
− x
i
′
,
(9.5)
56
gdzie symbol
( )
′
oznacza, ˙ze w iloczynie pomijamy składnik z i
= k. Podstawiaja
‘
c
do całki otrzymamy wz´or na kwadrature
‘
:
I
( f ) ≃
N
−1
∑
k
=0
A
k
f
(x
k
) ,
A
k
=
b
Z
a
L
k
(x) w(x) dx .
(9.6)
Wsp´ołczynniki kwadratury A
k
zale˙za
‘
od wyboru we
‘
zł´ow poprzez wielomiany L
k
(x).
Zauwa˙zmy, ˙ze powy˙zsza kwadratura jest dokładna dla wielomian´ow stopnia
< N, gdy˙z interpolacja Lagrange’a jest dokładna dla ka˙zdego z tych wielomian´ow.
9.1
Rza¸d kwadratury
Rza¸d kwadratury jest miara
‘
jej dokładno´sci.
Definicja
M´owimy, ˙ze kwadratura jest rze
‘
du r
= N je´sli jest dokładna dla wszystkich
wielomian´ow stopnia
< N oraz istnieje wielomian stopnia N, dla kt´orego
kwadratura nie jest dokładna.
Innymi słowy, rza
‘
d kwadratury jest okre´slony przez najni˙zszy stopie´n wielomianu,
dla kt´orego kwadratura nie jest dokładna.
Tak wie
‘
c dla kwadratury (9.6) zachodzi: r
≥ N. Łatwo udowodni´c, ˙ze dla
ka˙zdej kwadratury z N we
‘
złami r
≤ 2N. Istnieje bowiem wielomian stopnia 2N:
W
(x) = [(x − x
0
)(x − x
1
) . . . (x − x
N
−1
)]
2
≥ 0,
(9.7)
dla kt´orego ka˙zda kwadratura nie jest dokładna
Z
b
a
W
(x) w(x) dx > 0 ,
N
−1
∑
k
=0
A
k
W
(x
k
) = 0 .
(9.8)
Ostatecznie, rza
‘
d kwadratury (9.6) zawarty jest w przedziale
N
≤ r ≤ 2N
(9.9)
We
‘
zły
{x
k
} w kwadraturze Gaussa sa
‘
tak wybrane by jej rza
‘
d był r´owny maksy-
malnej warto´sci 2N. Mo˙zna to osia
‘
gna
‘
´c przy pomocy wielomian´ow ortogonal-
nych.
57
9.2
Wielomiany ortogonalne
Wielomiany okre´slone na odcinku
[a, b]
{P
0
(x) , P
1
(x) , . . . P
N
(x) , . . .}
(9.10)
tworza
‘
układ ortogonalny z waga
‘
w
(x) ≥ 0, je´sli dla i 6= j zachodzi
P
i
|P
j
≡
b
Z
a
P
i
(x) P
j
(x) w(x) dx = 0
(9.11)
Wielomiany ortogonalne maja
‘
bardzo wa˙zna
‘
własno´s´c wyra˙zona
‘
w naste
‘
puja
‘
-
cym twierdzeniu (dow ´od w podre
‘
czniku [?]).
Twierdzenie
Ka˙zdy wielomian ortogonalny P
n
(x) ma dokładnie n jednokrotnych pier-
wiastk´ow w przedziale
[a, b].
Poni˙zej w tabelce podajemy przykłady wielomian´ow ortogonalnych, poda-
ja
‘
c wage
‘
w
(x) oraz przedział [a, b] w iloczynie skalarnym (9.11), a tak˙ze relacje
‘
rekurencyjna
‘
, przy pomocy kt´orej mo˙zna obliczy´c kolejne wielomiany.
Wielomiany
P
n
w
(x)
[a, b]
Rekurencja
Legendre’a
P
n
1
[−1,1]
( j + 1)P
j
+1
= (2 j + 1)xP
j
− jP
j
−1
Czebyszewa
T
n
1
/
√
1
− x
2
[−1,1]
T
j
+1
= 2xT
j
− T
j
−1
Hermite’a
H
n
e
−x
2
[−
∞
,
∞
]
H
j
+1
= 2xH
j
− 2 jH
j
−1
Laguerre’a
L
n
e
−x
[0,
∞
]
jL
j
+1
= (2 j + 1 − x))L
j
− jL
(
α
)
j
−1
Przyklad
Dla wielomian´ow Legendre’a zachodzi
1
Z
−1
P
n
(x) P
m
(x) dx =
2
2n
+ 1
δ
nm
,
(9.12)
a kolejne wielomiany konstruowane przy pomocy relacji rekurencyjnej to
P
0
(x) = 1 ,
P
1
(x) = x ,
P
2
(x) =
1
2
(3x
2
− 1).
(9.13)
58
9.3
Kwadratura Gaussa
N poszukiwanych we
‘
zł´ow w kwadraturze Gaussa (9.6),
I
( f ) ≃
N
−1
∑
k
=0
A
k
f
(x
k
) ,
A
k
=
b
Z
a
L
k
(x) w(x) dx,
jest zadane przez zera wielomianu ortogonalnego P
N
z waga
‘
w
(x):
P
N
(x
k
) = 0
(9.14)
Dla takich we
‘
zł´ow rza
‘
d kwadratury Gaussa jest maksymalny i wynosi 2N.
Dow ´od
Poka˙zemy, ˙ze kwadratura Gaussa jest dokładna dla dowolnego wielomianu
W
(x) stopnia < 2N. Dziela
‘
c go przez P
N
(x) otrzymamy
W
(x) = W
1
(x) P
N
(x) + W
2
(x) ,
(9.15)
gdzie W
1
, W
2
sa
‘
wielomianami stopnia
< N. Mo˙zemy wtedy zapisa´c W
1
jako
kombinacje
‘
liniowa
‘
wielomian´ow ortogonalnych P
0
, . . . , P
N
−1
i wtedy
W
(x) =
N
−1
∑
i
=0
c
i
P
i
(x) P
N
(x) + W
2
(x)
(9.16)
Całkuja
‘
c z waga
‘
w
(x), otrzymamy
Z
b
a
W
(x) w(x) dx =
N
−1
∑
i
=0
c
i
Z
b
a
P
i
(x) P
N
(x) w(x) dx +
Z
b
a
W
2
(x) w(x) dx .
Suma całek znika na mocy ortogonalno´sci układu wielomian´ow P
i
, tak wie
‘
c
Z
b
a
W
(x) w(x) dx =
Z
b
a
W
2
(x) w(x) dx .
(9.17)
Kwadratura Gaussa jest dokładna dla wielomian´ow stopnia
< N, sta
‘
d
Z
b
a
W
2
(x) w(x) dx =
N
−1
∑
k
=0
A
k
W
2
(x
k
)
(9.18)
59
Wyliczaja
‘
c W
2
(x) z relacji (9.15) i podstawiaja
‘
c po prawej stronie, otrzy-
mamy
N
−1
∑
k
=0
A
k
W
2
(x
k
) =
N
−1
∑
k
=0
A
k
W
(x
k
) −
N
−1
∑
k
=0
A
k
W
1
(x
k
) P
N
(x
k
)
|
{z
}
0
.
(9.19)
Ostatnia suma jest r´owna zeru, gdy˙z P
N
(x
k
) = 0. Tak wie
‘
c udowodnili´smy,
˙ze kwadratura Gaussa jest dokładna, tzn. dla dowolnego wielomianu W
(x)
stopnia
< 2N zachodzi
Z
b
a
W
(x) w(x) dx =
N
−1
∑
k
=0
A
k
W
(x
k
) .
(9.20)
9.4
Wsp´ołczynniki kwadratury Gaussa
Wsp´ołczynniki A
k
kwadratury Gaussa sa
‘
zadane przez [?]:
A
k
=
hP
N
−1
|P
N
−1
i
P
N
−1
(x
k
) P
′
N
(x
k
)
(9.21)
gdzie k
= 0, 1, . . . (N − 1), natomiast P
′
N
to pochodna wielomianu P
N
.
Twierdzenie
Wszystkie wsp´ołczynniki kwadratury Gaussa sa
‘
dodatnie.
Dow ´od
Rozwa˙zmy kwadrat wielomianu L
i
(x) z interpolacji Lagrange’a (9.5), stop-
nia 2N. Kwadratura Gaussa jest dokładna dla tego wielomianu i sta
‘
d
Z
b
a
w
(x) [L
i
(x)]
2
dx
=
N
−1
∑
k
=0
A
k
[L
i
(x
k
)]
2
=
N
−1
∑
k
=0
A
k
δ
ik
= A
i
.
(9.22)
Całka po lewej stronie jest dodatnia i sta
‘
d A
i
> 0.
60
9.5
Przykład kwadratury Gaussa
Skonstruujemy kwadrature
‘
Gaussa-Legendre’a rze
‘
du r
= 4. Mamy wtedy do czy-
nienia z dwoma we
‘
złami
{x
0
, x
1
} i wtedy
1
Z
−1
f
(x) dx ≈ A
0
f
(x
0
) + A
1
f
(x
1
) .
(9.23)
We
‘
zły sa
‘
zerami wielomianu Legendre’a P
2
:
P
2
(x) =
1
2
(3x
2
− 1) = 0
(9.24)
i rozwia
‘
zaniami sa
‘
:
x
0
= −
1
√
3
,
x
1
=
1
√
3
.
(9.25)
Wsp´ołczynniki A
i
mo˙zna obliczy´c ze wzoru (9.21). Pro´sciej jest wykorzysta´c fakt,
˙ze kwadratura (9.23) jest dokładna dla wielomian´ow stopnia
< 4, w szczeg´olno´sci
dla W
0
(x) = 1 oraz W
1
(x) = x. Wtedy
A
0
+ A
1
=
Z
1
−1
dx
= 2
−
A
0
√
3
+
A
1
√
3
=
Z
1
−1
x dx
= 0 .
Sta
‘
d wsp´ołczynniki kwadratury A
1
= A
2
= 1 i ostatecznie
1
Z
−1
f
(x) dx ≈ f (−
1
√
3
) + f (
1
√
3
) .
(9.26)
61
Rozdział 10
Zera funkcji
Nie ka˙zde r´ownanie (lub układ r´owna´n) mo˙zna rozwia
‘
za´c dokładnie. Na przykład,
nie mo˙zna poda´c og´olnych wzor´ow na rozwia
‘
zanie dowolnego r´ownanie alge-
braicznego stopnia wy˙zszego ni˙z cztery.
Rozpatrzmy zagadnienie znajdowania pierwiastka r´ownania
f
( ¯x) = 0
(10.1)
w pewnym przedziale. Zakładamy, ˙ze jest to pierwiastek jednokrotny, tzn. w jego
otoczeniu
f
(x) ∼ (x − ¯x).
(10.2)
Prezentowane poni˙zej metody mo˙zna stosowa´c gdy znamy przedział, w kt´orym
znajduje sie
‘
pierwiastek. Nale˙zy wie
‘
c wcze´sniej okre´sli´c taki przedział, na przykład
rysuja
‘
c wykres funkcji y
= f (x).
10.1
Metoda połowienia
Zakładamy, ˙ze funkcja f
(x) jest cia
‘
gła w przedziale
[a, b]. Zgodnie z twierdze-
niem Bolzano-Cauchego, je´sli na ko´ncach tego przedziału warto´sci funkcji maja
‘
przeciwne znaki,
f
(a) f (b) < 0 ,
(10.3)
to wewna
‘
trz przedziału znajduje sie
‘
co najmniej jeden pierwiastek r´ownania (10.1).
62
Zgodnie z uwaga
‘
z poprzedniego rozdziału zakładamy, ˙ze jest to pierwiastek
jednokrotny. W pierwszym kroku dzielimy przedział
[a, b] na połowe
‘
x
0
=
a
+ b
2
.
(10.4)
Je˙zeli f
(x
0
) = 0 to x
0
jest poszukiwanym pierwiastkiem. W przeciwnym wypadku
pierwiastek le˙zy w tym z przedział´ow
[a, x
0
]
lub
[x
0
, b] ,
na ko´ncach kt´orego funkcja f ma przeciwny znak. Dzielimy ten przedział na
połowe
‘
otrzymuja
‘
c x
1
. Zauwa˙zmy, ˙ze długo´s´c nowego przedziału to
∆
1
=
1
2
(b − a).
(10.5)
W wyniku wielokrotnego zastosowania tej procedury otrzymujemy pierwiastek ¯
x
lub cia
‘
g punkt´ow x
n
, be
‘
da
‘
cych ´srodkiem przedziału o długo´sci
∆
n
=
1
2
n
(b − a)
(10.6)
Po dostatecznie du˙zej liczbie krok´ow długo´s´c takiego przedziału jest dowolnie
mała. Zadaja
‘
c wie
‘
c dokładno´s´c
ε
, przerywamy procedure
‘
przy n takim, ˙ze
∆
n
≤ 2
ε
.
(10.7)
Poszukiwanym pierwiastkiem jest wtedy
¯
x
= x
n
±
ε
,
(10.8)
Przyklad
Rozwia
‘
˙zemy r´ownanie exp
(−x) = x. W tym celu poszukajmy zer funkcji
f
(x) = x − exp(−x).
(10.9)
Zero znajduje sie
‘
przedziale
[0, 1] gdy˙z f (0) = −1 i f (1) ≈ 0.63. Kolejne
warto´sci ´srodk´ow przedziałow x
n
wraz z długo´sciami
∆
n
/2 to
63
n
x
n
∆
n
/2
0
0.50000
0.50000
1
0.75000
0.25000
2
0.62500
0.12500
3
0.56250
0.06250
4
0.59375
0.03125
5
0.57812
0.01562
6
0.57031
0.00781
7
0.56641
0.00391
8
0.56836
0.00195
W 18 kroku otrzymujemy liczbe
‘
0.56714, kt´ora ju˙z nie ulega zmianie przy
zadanej liczbie pie
‘
ciu cyfr po przecinku (dokładno´s´c
ε
< 10
−5
).
10.2
Metoda Newtona
Zał´o˙zmy, ˙ze f jest klasy C
2
w przedziale
[a, b] na ko´ncach, kt´orego zmienia znak.
Ponadto, niech pochodne f
′
oraz f
′′
maja
‘
stały znak w całym przedziale.
Metoda Newtona obejmuje cztery przypadki be
‘
da
‘
ce kombinacja
‘
naste
‘
puja
‘
cych
warunk´ow. Funkcja f jest
• rosna
‘
ca oraz wypukła ku dołowi
( f
′
> 0, f
′′
> 0)
• rosna
‘
ca oraz wypukła ku g´orze
( f
′
> 0), f
′′
< 0)
• maleja
‘
ca oraz wypukła ku dołowi
( f
′
< 0, f
′′
> 0)
• maleja
‘
ca oraz wypukła ku g´orze
( f
′
< 0, f
′′
< 0).
Przyjmuja
‘
c dla ustalenia uwagi, ˙ze obie pochodne sa
‘
dodatnie, wystawmy
styczna
‘
do wykresu funkcji w punkcie x
0
= b, w kt´orym f (x
0
) > 0, patrzy rysunek
10.1,
y
− f (x
0
) = f
′
(x
0
)(x − x
0
) .
(10.10)
Kłada
‘
c y
= 0 otrzymujemy punkt przecie
‘
cia stycznej z osia
‘
x:
x
1
= x
0
−
f
(x
0
)
f
′
(x
0
)
.
(10.11)
Udowodnimy, ˙ze ¯
x
< x
1
< x
0
, gdzie ¯
x jest poszukiwanym pierwiastkiem.
64
Rys. 10.1: Ilustracja metody Newtona
Dow ´od
G ´orne ograniczenie wynika ze wzoru (10.11), gdy˙z f
(x
0
) oraz f
′
(x
0
) sa
‘
wie
‘
ksze od zera. W dowodzie dolnego ograniczenia rozwa˙zmy rozwinie
‘
cia
Taylora funkcji f wok´oł punktu x
0
:
f
(x) = f (x
0
) + f
′
(x
0
)(x − x
0
) +
1
2
f
′′
(
ξ
)(x − x
0
)
2
,
gdzie
ξ
∈ ( ¯x,x
0
). Kłada
‘
c x
= ¯x i wykorzystuja
‘
c r´owno´s´c f
( ¯x) = 0, wyliczymy
¯
x
= x
0
−
f
(x
0
)
f
′
(x
0
)
|
{z
}
x
1
−
1
2
f
′′
(
ξ
)
f
′
(x
0
)
( ¯x − x
0
)
2
.
Sta
‘
d na podstawie zało˙zenia f
′
, f
′′
> 0 otrzymujemy nasza
‘
teze
‘
¯
x
− x
1
= −
1
2
f
′′
(
ξ
)
f
′
(x
0
)
( ¯x − x
0
)
2
< 0 .
Powt´orzmy procedure
‘
, wystawiaja
‘
c styczna
‘
w punkcie
(x
1
, f (x
1
)):
y
− f (x
1
) = f
′
(x
1
)(x − x
1
) .
(10.12)
Przecina ona o´s x w punkcie
x
2
= x
1
−
f
(x
1
)
f
′
(x
1
)
.
(10.13)
65
Punkt ten spełnia warunek ¯
x
< x
2
< x
1
. Dolny warunek udowodnimy podobnie
jak dla punktu x
1
. Dow ´od g´ornego warunku polega na pokazaniu, ˙ze f
(x
1
) > 0.
W tym celu skorzystamy z twierdzenia Lagrange’a, kt´ore m ´owi, ˙ze istnieje punkt
ζ
∈ ( ¯x,x
1
), dla kt´orego zachodzi
f
(x
1
) − f ( ¯x) = f
′
(
ζ
)(x
1
− ¯x).
Ze wzgle
‘
du na warunki f
( ¯x) = 0 oraz f
′
(
ζ
) > 0 dostajemy wie
‘
c f
(x
1
) > 0.
Kolejne kroki procedury prowadza
‘
do relacji rekurencyjnej definiuja
‘
cej kolejne
wyrazy cia
‘
gu przybli˙ze´n x
n
poszukiwanego zera funkcji:
x
n
+1
= x
n
−
f
(x
n
)
f
′
(x
n
)
.
(10.14)
Otrzymany cia
‘
g punkt´ow jest maleja
‘
cy i ograniczony od dołu. Z twierdzenia
Cauchego wynika, ˙ze istnieje granica tego cia
‘
gu g. Tak wie
‘
c z relacji (10.14)
wynika r´owno´s´c
g
= g −
f
(g)
f
′
(g)
.
(10.15)
Sta
‘
d wnioski
f
(g) = 0
=>
g
= ¯x.
(10.16)
Procedura Newtona jest wie
‘
c zbie˙zna do poszukiwanego zera funkcji. W praktyce
procedure
‘
przerywamy gdy
|x
n
+1
− x
n
| <
ε
.
(10.17)
Przyklad
Metoda Newtona jest znacznie szybciej zbie˙zna ni˙z metoda połowienia. Dla
przykładu z poprzedniego rozdziału otrzymujemy wynik po czterech krokach
(pamie
‘
tamy, ˙ze x
0
= 1):
n
x
n
|x
n
− x
n
−1
|
1
0.53788
0.46212
2
0.56699
0.02910
3
0.56714
0.00016
4
0.56714
0.00000
66
10.3
Metoda siecznych
W metodzie Newtona konieczna jest znajomo´s´c pochodnej f
′
. W metodzie sie-
cznych unikamy tego warunku przybli˙zaja
‘
c pochodna
‘
w wyra˙zeniu (10.14) przez
iloraz r´o˙znicowy
f
′
(x
n
) ≃
f
(x
n
) − f (x
n
−1
)
x
n
− x
n
−1
.
(10.18)
W ten spos´ob otrzymujemy
x
n
+1
= x
n
− f (x
n
)
x
n
− x
n
−1
f
(x
n
) − f (x
n
−1
)
.
(10.19)
Zatem do wyznaczenia
(n + 1) przybli˙zenia pierwiastka ¯x wykorzystuje sie
‘
punkty
(x
n
, f (x
n
)) oraz (x
n
−1
, f (x
n
−1
)), przez kt´ore przeprowadza sie
‘
sieczna
‘
y
− f (x
n
) =
f
(x
n
) − f (x
n
−1
)
x
n
− x
n
−1
(x − x
n
) .
(10.20)
Przecinaja ona o´s x w punkcie x
n
+1
zadanym wzorem (10.19). Aby rozpocza
‘
c
procedure
‘
konieczne jest wie
‘
c wybranie dw´och punkt´ow startowych x
0
i x
1
. W
omawianej przez nas przykładzie wybieramy x
0
= b oraz x
1
= x
0
−
∆
z mała war-
to´scia
‘
∆
≪ x
0
.
Metoda siecznych mo˙ze nie by´c zbie˙zna, na przykład, gdy pocza
‘
tkowe przy-
bli˙zenia nie le˙za
‘
dostatecznie blisko szukanego pierwiastka. Jako dodatkowe kry-
terium przerwania iteracji, opr´ocz warto´sci r´o˙znic
|x
n
+1
− x
n
|, nale˙zy przyja
‘
´c war-
to´sci
| f (x
n
)|, tak by tworzyły one cia
‘
g maleja
‘
cy w ko´ncowej fazie oblicze´n.
Wracaja
‘
c do przykładu, w metodzie siecznych otrzymujemy wynik 0
.56714 z
błe
‘
dem mniejszym ni˙z 10
−5
po pie
‘
ciu krokach
10.4
Metoda fałszywej prostej (falsi)
W metodzie falsi (fałszywej prostej) znajomo´s´c f
′
tak˙ze nie jest potrzebna. Przy
przyje
‘
tych zało˙zeniach ( f
′
, f
′′
> 0) przeprowadzamy w pierwsza
‘
sieczna
‘
przez
punkty
(b, f (b)) oraz (a, f (a)):
y
− f (a) =
f
(b) − f (a)
b
− a
(x − a).
(10.21)
67
Rys. 10.2: Ilustracja metody falsi
Przecina ona o´s x w punkcie x
1
, patrz rysunek 10.2,
x
1
= a − f (a)
b
− a
f
(b) − f (a)
.
(10.22)
Punkt
(b, f (b)) jest punktem stałym wszystkich cie
‘
ciw i w naste
‘
pnym przybli˙zeniu
przeprowadzamy cie
‘
ciwe
‘
przez punkt
(x
1
, f (x
1
), otrzymuja
‘
c x
2
jako punkt jej prze-
cie
‘
cia z osia
‘
x. W og´olno´sci przeprowadzamy cie
‘
ciwy przez punkty
(x
n
, f (x
n
)):
y
− f (x
n
) =
f
(b) − f (x
n
)
b
− x
n
(x − x
n
) .
(10.23)
Przecinaja
‘
one o´s x w punkcie:
x
n
+1
= x
n
− f (x
n
)
b
− x
n
f
(b) − f (x
n
)
(10.24)
Mo˙zna pokaza´c, ˙ze metoda ta prowadzi do cia
‘
gu warto´sci x
n
zbie˙znych do granicy
be
‘
da
‘
cej jednokrotnym pierwiastkiem r´ownania f
( ¯x) = 0.
Przyklad
W naszym przykładzie otrzymujemy w metodzie falsi:
68
n
x
n
|x
n
− x
i
−1
|
1
0.61270
0.61270
2
0.56384
0.04886
3
0.56739
0.00355
4
0.56713
0.00026
5
0.56714
0.00002
6
0.56714
0.00000
Jak wida´c metoda falsi jest stosunkowo wolno zbie˙zna.
69
Rozdział 11
R´ownania r´o˙zniczkowe
11.1
R´ownania zwyczajne pierwszego rze¸du
R ´ownanie r´o˙zniczkowe zwyczajne rze
‘
du pierwszego ma posta´c
dy
dt
= F(t, y) .
(11.1)
Funkcja y
= y(t) jest rozwia
‘
zaniem je´sli po podstawieniu do (11.1) otrzymujemy
to˙zsamo´s´c.
dy
(t)
dt
= F(t, y(t)) .
(11.2)
Rozwia
‘
zanie jest wyznaczone jednoznacznie poprzez zadanie warunku pocza
‘
tko-
wego:
y
(t
0
) = y
0
.
(11.3)
Przyklad
Rozwia
‘
˙zmy r´ownanie
dy
dt
= y
(11.4)
z warunkiem pocza
‘
tkowym y
(t
0
) = y
0
. Rozwia
‘
zaniem og´olnym jest funkcja
y
(t) = C exp(t) ,
(11.5)
co łatwo sprawdzi´c poprzez podstawienie do r´ownania. Dla chwili pocza
‘
tkowej
t
0
otrzymujemy
y
(t
0
) = C exp(t
0
) = y
0
.
(11.6)
70
Sta
‘
d wynika, ˙ze
C
= y
0
exp
(−t
0
)
(11.7)
i rozwia
‘
zaniem spełniaja
‘
cym warunek pocza
‘
tkowy jest
y
(t) = y
0
exp
(t − t
0
) .
(11.8)
11.2
Układ r´owna ´n r´o˙zniczkowych
Układ r´owna´n r´o˙zniczkowych 1. rze
‘
du ma naste
‘
puja
‘
ca
‘
posta´c
dy
1
dt
= F
1
(t, y
1
, y
2
, . . . y
n
)
dy
2
dt
= F
2
(t, y
1
, y
2
, . . . y
n
)
...
dy
n
dt
= F
n
(t, y
1
, y
2
, . . . y
n
) .
(11.9)
Wprowadzaja
‘
c notacje
‘
wektorowa
‘
y
= (y
1
, y
2
, . . . y
n
) ,
F
= (F
1
, F
2
, . . . F
n
)
(11.10)
układ ten mo˙zemy zapisa´c w formie analogocznej do r´ownania (11.1)
dy
dt
= F(t, y) .
(11.11)
Rozwia
‘
zanie układu r´ownan´n y
= y(t) jest wyznaczone jednoznacznie poprzez
warunek pocza
‘
tkowy
y
(t
0
) = (y
1
(t
0
), y
2
(t
0
), . . . , y
n
(t
0
)) ≡ y
0
.
(11.12)
11.3
R´ownania wy˙zszych rze¸d´ow
R ´ownanie r´o˙zniczkowe zwyczajne rze
‘
du n ma posta´c
d
n
y
dt
n
= F(t, y, y
(1)
, y
(2)
, . . . y
(n−1)
) ,
(11.13)
71
gdzie y
(k)
jest k-ta
‘
pochodna
‘
. Rozwia
‘
zanie y
= y(t) jest wyznaczone jednoznacznie
przez warunki pocza
‘
tkowe
y
(t
0
) = y
10
,
y
(1)
(t
0
) = y
20
,
. . .
y
(n−1)
(t
0
) = y
n0
.
(11.14)
Ka˙zde r´ownanie r´o˙zniczkowe zwyczajne rze
‘
du n mo˙zna zapisa´c jako układ n
r´owna´n 1. rze
‘
du. Wprowad´zmy bowiem oznaczenia
y
1
= y ,
y
2
= y
(1)
,
. . .
y
n
= y
(n−1)
.
(11.15)
Wtedy r´ownanie (11.13) mo˙zna zapisa´c w r´ownowa˙znej formie układu r´owna´n 1.
rze
‘
du:
dy
1
dt
= y
2
dy
2
dt
= y
3
...
dy
n
dt
= F(t, y
1
, y
2
, . . . y
n
) .
(11.16)
Z obserwacji tej wynika wniosek, ˙ze r´ownanie (11.1) jest podstawowym obiektem
zainteresowa´n. Rozwinie
‘
cie metod przybli˙zonego rozwia
‘
zywania tego r´ownania
pozwala na znalezienie rozwia
‘
za´n dla układu r´owna´n (11.11), a tym samym dla
r´owna´n wy˙zszych rze
‘
d´ow (11.13).
Przyklad
Rozwa˙zmy układ trzech r´owna´n Newtona (masa m
= 1) drugiego rze
‘
du
d
2
r
dt
2
= F
t
, r,
dr
dt
.
(11.17)
Definiuja
‘
c v
= dr/dt dostajemu układ sze´sciu r´owna´n 1. rze
‘
du na poszuki-
wane funkcje
(r(t), v(t)):
dr
dt
= v
dv
dt
= F(t, r, v) .
(11.18)
72
11.4
Metody Eulera
Prezentowane tutaj metody znajdowania przybli˙zonych rozwia
‘
za´n r´ownania (11.1)
bazuja
‘
na dyskretyzacji przedziału czasu
[t
0
, T ], w kt´orym chcemy znale´z´c rozwia
‘
-
zanie
t
0
< t
1
. . . < t
n
−1
< t
n
. . . < T .
(11.19)
Zał´o˙zmy, ˙ze punkty czasowe sa
‘
r´ownoodległe:
t
n
− t
n
−1
=
τ
.,
(11.20)
Omawiane metody dostarczaja
‘
rekurencji wia
‘
˙za
‘
cej warto´s´c rozwia
‘
zania y
n
+1
w
chwili t
n
+1
z warto´sciami rozwia
‘
zania y
n
w chwili wcze´sniejszej t
n
:
y
n
+1
= f (t
n
, y
n
) .
(11.21)
Warto´s´c y
0
= y(t
0
) jest zadana przez warunek pocza
‘
tkowy (11.3). Bardzo wa˙znym
zagadnieniem tak okre´slonych metod jest numeryczna stabilno´s´c rozwia
‘
zania dla
du˙zych czas´ow T .
Punktem wyj´scia metod Eulera jest r´ownanie otrzymane po scałkowaniu po
czasie obu stron r´ownania (11.1) w przedziale
(t
n
, t
n
+1
):
y
n
+1
= y
n
+
t
n
+1
Z
t
n
F
(t, y(t)) dt .
(11.22)
Poniewa˙z nie znamy y
(t) w funkcji podcałkowej musimy zastosowa´c metode
‘
przy-
bli˙zona
‘
obliczenia tej całki. Wykorzystuja
‘
c wz´or (8.6) na całkowanie metoda
‘
pros-
toka
‘
t´ow,
t
n
+1
Z
t
n
dt F
(t, y(t)) ≃
τ
F
(t
n
, y
n
) .,
(11.23)
otrzymujemy relacje
‘
rekurencyjna
‘
w metodzie Eulera:
y
n
+1
= y
n
+
τ
F
(t
n
, y
n
)
(11.24)
Zauwa˙zmy, ˙ze mogliby´smy przybli˙zy´c całke
‘
(11.23) przez warto´s´c funkcji dla
g´ornej granicy:
y
n
+1
= y
n
+
τ
F
(t
n
+1
, y
n
+1
) .
(11.25)
Pojawia sie
‘
jednak w tym momencie problem, gdy˙z nieznane y
n
+1
wyste
‘
puje po
obu stronach r´ownania. Nale˙załoby wie
‘
c u˙zy´c dodatkowo metody poszukiwania
zer z rozdziału 10. Jest to do´s´c niepraktyczne, chocia˙z r´ownanie (11.25) mo˙ze by´c
lepsze ze wzgle
‘
du na stabilno´s´c numeryczna
‘
rozwia
‘
zania. Zagadnienie to zostanie
om ´owione w naste
‘
pnym rozdziale.
73
11.5
Metoda przewid´z i popraw
Przybli˙zmy całke
‘
w r´ownaniu (11.22) stosuja
‘
c metoda
‘
trapez´ow (8.10):
t
n
+1
Z
t
n
F
(t, y(t)) dt ≃
τ
2
(F(t
n
, y
n
) + F(t
n
+1
, y
n
+1
)) .
(11.26)
Sta
‘
d
y
n
+1
= y
n
+
τ
2
{F(t
n
, y
n
) + F(t
n
+1
, y
n
+1
)} .
(11.27)
Podobnie jak we wzorze (11.25), nieznana warto´s´c y
n
+1
znajduje sie
‘
po obu stronach
r´ownania rekurencyjnego. Nale˙zy wie
‘
c wykorzysta´c jedna
‘
z metod znajdowania
zer funkcji z rozdziału 10.
Alternatywna
‘
metoda
‘
jest dwustopniowa procedura. Najpierw okre´slamy prze-
widywa
‘
warto´s´c y
n
+1
korzystaja
‘
c z metody Eulera (11.24):
¯
y
n
+1
= y
n
+
τ
F
(t
n
, y
n
)
(11.28)
a naste
‘
pnie obliczamy poprawiona
‘
warto´s´c:
y
n
+1
= y
n
+
τ
2
{F(t
n
, y
n
) + F(t
n
+1
, ¯y
n
+1
)}
(11.29)
Sta
‘
d nazwa metody przewid´z i popraw.
11.6
Stabilno´s´c numeryczna rozwia¸za ´n
Rozwa˙zmy r´ownanie r´o˙zniczkowe
dy
dt
= −
λ
y
,
λ
> 0 .
(11.30)
Rozwia
‘
zanie dokładne z warunkiem pocza
‘
tkowym y
(0) = y
0
to
y
(t) = y
0
e
−
λ
t
.
(11.31)
Da
‘
˙zy ono do zera dla t
→
∞
.
Szukaja
‘
c rozwia
‘
zania metoda
‘
Eulera (11.24), otrzymujemy
y
n
+1
= y
n
−
λτ
y
n
= (1 −
λτ
) y
n
,
(11.32)
74
co prowadzi do naste
‘
puja
‘
cej relacji dla kolejnych chwil t
n
= n
τ
:
y
n
= (1 −
λτ
)
n
y
0
.
(11.33)
Powy˙zsze rozwia
‘
znie dobrze przybli˙za rozwia
‘
zanie dokładne (11.31) gdy zachodzi
|1 −
λτ
| < 1.
(11.34)
Warunek ten oznacza, ˙ze krok czasowu musi by´c dostatecznie mały
τ
<
2
λ
.
(11.35)
W przeciwnym przypadku
|y
n
| →
∞
dla rosna
‘
cych n.
Najprostszym rozwia
‘
zaniem dla tego problemu (opr´ocz dobrania
τ
≪ 2/
λ
) jest
skorzystanie z drugiego wzoru Eulera (11.25). Wtedy otrzymujemy dla naszego
r´ownania
y
n
+1
= y
n
−
λτ
y
n
+1
,
(11.36)
co prowadzi do
y
n
+1
=
y
n
1
+
λτ
=
1
1
+
λτ
n
+1
y
0
.
(11.37)
Warto´s´c 1
/(1 +
λτ
) < 1 i metoda ta jest stabilna niezale˙znie od kroku
τ
oraz
warto´sci
λ
.
11.7
R´ownania typu stiff
Rozwa˙zmy układ dw ´och r´owna´n r´o˙zniczkowych
dy
dt
= −Ay
(11.38)
z macierza
‘
A o dodatnich warto´sciach własnych
λ
1
oraz
λ
2
. Rozwia
‘
zanie jest suma
‘
dw ´och rozwia
‘
za´n:
y
(t) = e
−
λ
1
t
y
1
+ e
−
λ
2
t
y
2
.
(11.39)
Je˙zeli
λ
2
≫
λ
1
to drugie rozwia
‘
zanie szybko da
‘
zy do zera i przestaje by´c istotne.
W numerycznej realizacji drugie rozwia
‘
zanie mo˙ze jednak odgrywa´c kluczowa
‘
role
‘
, prowadza
‘
c do niestabilnego zachowania całego rozwia
‘
zania. Stosuja
‘
c metode
‘
Eulera (11.24), otrzymujemy
y
n
= (1 −
λ
1
τ
)
n
y
1
+ (1 −
λ
2
τ
)
n
y
2
.
(11.40)
75
Je´sli krok czasowy
τ
> 2/
λ
2
to wyra˙zenie
|1−
λ
2
τ
| > 1, co prowadzi do rosna
‘
cych
i oscyluja
‘
cych warto´sci y
n
gdy n
→
∞
. Nawet dla
τ
≤ 2/
λ
2
rozwia
‘
zanie zachowuje
charakter oscylacyjny dla niezbyt długich czas´ow.
Podsumowuja
‘
c, drugie rozwia
‘
zanie okre´sla krok czasowy przy kt´orym rozwia
‘
-
zanie (11.40) jest stabilne numerycznie
τ
≪ 2/
λ
2
.
(11.41)
Cena
‘
jest zwykle
‘
mała warto´s´c
τ
, czyli du˙za liczba krok´ow n
= t/
τ
by osia
‘
gnac
ko´ncowe t.
Przyklad
Je˙zeli
y
(t) = e
−t
y
1
+ e
−1000 t
y
2
,
(11.42)
to krok
τ
≪ 2 · 10
−3
, a całkowita liczba krok´ow n
≫ 10
3
dla t
∼ 1.
76
Rozdział 12
Metoda Rungego-Kutty
Zastosujmy rozwinie
‘
cie Taylora do rozwia
‘
zania r´ownania (11.1):
y
(t
n
+1
) = y(t
n
+
τ
) = y(t
n
) +
τ
dy
(t
n
)
dt
+
τ
2
2
d
2
y
(t
n
)
dt
2
+ . . . .
(12.1)
Pierwsza pochodna to
dy
(t
n
)
dt
= F(t
n
, y
n
) .
(12.2)
Natomiast druga pochodna to
d
2
y
dt
2
=
dF
(t, y)
dt
=
∂
F
∂
t
+
∂
F
∂
y
dy
dt
=
∂
F
∂
t
+
∂
F
∂
y
F
.
(12.3)
Sta
‘
d wz´or (12.1) przyjmuje posta´c
y
n
+1
= y
n
+
τ
F
(t
n
, y
n
) +
τ
2
2
{
∂
t
F
+ F
∂
y
F
}(t
n
, y
n
) +
O
(
τ
3
) ,
(12.4)
Jest to wz´or ´scisły i mo˙zmy go wykorzysta´c do znajdowania rozwia
‘
zania. Wada
‘
tej metody jest konieczno´s´c znajomo´sci pochodnych funkcji F .
12.1
Metoda drugiego rze¸du
Metoda Rungego-Kutty pozwala unikna
‘
´c problemu pochodnych po prawej stronie
wzoru (12.4) poprzez wzie
‘
cie warto´sci funkcji F po prawej stronie w odpowiednio
dobranych punktach pomie
‘
dzy t
n
i t
n
+1
.
77
W metodzie drugiego rze
‘
du wybiera sie
‘
dwa takie punkty:
k
1
=
τ
F
(t
n
, y
n
)
k
2
=
τ
F
(t
n
+ a
τ
, y
n
+ a k
1
)
y
n
+1
= y
n
+
α
1
k
1
+
α
2
k
2
,
(12.5)
gdzie wsp´ołczynniki
α
1
,
α
2
i a sa
‘
tak dobrane by spełni´c r´ownanie (12.4) po roz-
winie
‘
ciu (12.5) w szereg w zmiennej
τ
. Rozwijaja
‘
c k
2
z dokładno´scia
‘
O
(
τ
3
), otrzy-
mujemy:
k
2
=
τ
{F + a
τ ∂
t
F
+ a k
1
∂
y
F
}
=
τ
{F + a
τ ∂
t
F
+ a
τ
F
∂
y
F
}
=
τ
F
+ a
τ
2
(
∂
t
F
+ F
∂
y
F
) ,
(12.6)
gdzie wszystkie funkcje po prawej stronie sa
‘
wzie
‘
te w punkcie
(t
n
, y
n
). Podstawia-
ja
‘
c k
1
i k
2
do wzoru (12.5), otrzymujemy
y
n
+1
= y
n
+
α
1
τ
F
+
α
2
τ
F
+ a
τ
2
(
∂
t
F
+ F
∂
y
F
)
= y
n
+ (
α
1
+
α
2
)
τ
F
+
α
2
a
τ
2
(
∂
t
F
+ F
∂
y
F
) .
(12.7)
Zgodno´s´c z r´ownaniem (12.4) wymaga spełnienia warunk´ow
α
1
+
α
2
= 1
α
2
a
=
1
2
.
(12.8)
Jednym z mo˙zliwych wybor´ow jest
α
1
=
α
2
=
1
2
a
= 1 .
(12.9)
Otrzymujemy w ten spos´ob wzory dla metody Rungego-Kutty drugiego rze
‘
du:
k
1
=
τ
F
(t
n
, y
n
)
(12.10)
k
2
=
τ
F
(t
n
+
τ
, y
n
+ k
1
)
(12.11)
y
n
+1
= y
n
+
1
2
(k
1
+ k
2
)
(12.12)
78
Jak łatwo zauwa˙zy´c wzory te sa
‘
identyczne ze wzorami w metodzie przewid´z
i popraw, gdy˙z dla argument´ow funkcji f we wzorze (12.11), zachodzi
τ
n
+
τ
=
τ
n
+1
,
y
n
+ k
1
= y
n
+
τ
F
(t
n
, y
n
) = ¯y
n
+1
i sta
‘
d otrzymujemy wz´or (11.29)
y
n
+1
= y
n
+
τ
2
{F(t
n
, y
n
) + F(t
n
+1
, ¯y
n
+1
)}
Dla układu r´owna´n r´o˙zniczkowych (11.11) otrzymujemy wzory analogiczne
do powy˙zszych, zapisane w formie wektorowej
k
1
=
τ
F
(t
n
, y
n
)
k
2
=
τ
F
(t
n
+
τ
, y
n
+ k
1
)
y
n
+1
= y
n
+
1
2
(k
1
+ k
2
)
(12.13)
12.2
Metoda czwartego rze¸du
W metodzie czwartego rze
‘
du otrzymuje sie
‘
naste
‘
puja
‘
ce wzory dla układu r´owna´n
r´o˙zniczkowych:
k
1
=
τ
F
(t
n
, y
n
)
k
2
=
τ
F
(t
n
+
1
2
τ
, y
n
+
1
2
k
1
)
k
3
=
τ
F
(t
n
+
1
2
τ
, y
n
+
1
2
k
2
)
k
4
=
τ
F
(t
n
+
τ
, y
n
+ k
3
)
y
n
+1
= y
n
+
1
6
(k
1
+ 2k
2
+ 2k
3
+ k
4
) .
(12.14)
Relacja (12.14) jest zgodna z rozwie
‘
ciem Taylora (12.1) rozwia
‘
zania r´ownania
r´o˙zniczkowego a˙z do wyraz´ow rze
‘
du
τ
4
.
79
Rozdział 13
Metody Monte Carlo
Podstawowym elementem metod Monte Carlo sa
‘
liczby losowe, kt´ore sa
‘
genero-
wane z pewnym prawdopodobie´nstwem. Metody te stosuje sie
‘
w takich zagad-
nieniach jak
• symulacja proces´ow o charakterze stochastycznym, w kt´orych wynik po-
jawia sie
‘
z okre´slonym prawdopodobie´nstwem,
• obliczanie całek. Zaleta
‘
metod Monte Carlo jest mo˙zliwo´s´c obliczania wie-
lowymiarowych całek po skomplikowanych obszarach w wielowymiarowej
przestrzeni. Na og´oł inne metody całkowania zawodza
‘
w takich przypad-
kach.
W og´olno´sci, bardziej og´olnym poje
‘
ciem jest zmienna losowa. Warto´sciami
tej zmiennej sa
‘
wła´snie liczby losowe.
13.1
Zmienna losowa i rozkład prawdopodobie ´nstwa
Zmienna losowa X przyjmuje warto´sci liczbowe z pewnym prawdopodobie ´n-
stwem.
Innymi słowy, zmienna losowa jest w pełni okre´slona, gdy podamy zakres
warto´sci liczbowych, kt´ore mo˙ze przyjmowa´c oraz dodatnia
‘
funkcje
‘
P
(X = x),
przy pomocy kt´orej okre´slamy prawdopodbie´nstwo, ˙ze zmienna losowa X przyj-
muje warto´s´c liczbowa
‘
x.
80
Przyklad
Zmienna
‘
losowa
‘
jest wynik rzutu kostka
‘
. Przyjmuje ona sze´s´c warto´sci:
x
∈ {1,2,3,4,5,6}.
(13.1)
Zakladaja
‘
c, ˙ze kostka jest idealna okre´slamy prawdopodobie´nstwo, i˙z zmien-
na losowa przyjmie dozwolone warto´sci jako
P
(X = x) =
1
6
.
(13.2)
Przykład ten jest ilustracja
‘
zmiennej losowej o warto´sciach dyskretnych.
W przypadku gdy zmienna losowa przyjmuje warto´sci cia
‘
głe x
∈
ℜ
, praw-
dopodbie´nstwo, ˙ze warto´sci zmiennej losowej X sa
‘
z przedziału
[a, b] jest okre´slone
przez
P
(a ≤ X ≤ b) =
b
Z
a
p
(x) dx
(13.3)
Dodatnia funkcja p
(x) jest nazywana ge
‘
sto´scia
‘
prawdopodobie ´nstwa zmiennej
losowej X . Jest ona unormowana do jedynki
∞
Z
−
∞
p
(x) dx = 1 .
(13.4)
Rozkłady zmiennej losowej sa
‘
charakteryzowane przez takie wielko´sci jak,
• warto´s´c ´srednia zmiennej losowej
hXi =
∞
Z
−
∞
x p
(x) dx .
(13.5)
• wariancja zmiennnej losowej
σ
2
X
=
∞
Z
−
∞
(x − hXi)
2
p
(x) dx =
X
2
− hXi
2
.
(13.6)
Wariancja jest miara
‘
rozrzutu warto´sci zmiennej losowej wok´oł jej warto´sci
´sredniej
• w og´olno´sci definiujemy momenty zmiennej losowej
hX
n
i =
∞
Z
−
∞
x
n
p
(x) dx ,
n
= 1, 2, . . .
(13.7)
Gdy powy˙zsze całki nie sa
‘
zbie˙zne momenty nie istnieja
‘
.
81
13.2
Przykłady rozkład´ow prawdopodobie ´nstwa
Przykładami rozkład´ow zmiennych losowych o cia
‘
głych warto´sciach sa
‘
.
• Rozkład jednostajny:
p
(x) =
1
x
∈ [0,1]
0
poza
.
(13.8)
Warto´s´c ´srednia
hXi = 1/2, natomiast wariancja
σ
2
X
= 1/12.
• Rozkład normalny Gaussa z parametrami µ oraz
σ
2
:
p
(x) =
1
√
2
πσ
2
e
−
(x−µ)
2
2
σ
2
(13.9)
Warto´s´c ´srednia
hXi = µ, natomiast wariancja
σ
2
X
=
σ
2
.
Przykładami rozkład´ow zmiennych losowych o warto´sciach dyskretnych sa
‘
.
• Rozkład dwumianowy.
Rozwa˙zmy N niezale˙znych pr´ob (np. rzut moneta
‘
). Pytamy jakie jest praw-
dopodobie´nstwo odniesienia n sukces´ow, je´sli prawdopodbie´nstwo pojedyn-
czego sukcesu wynosi q. Zmienna losowa X opisuje liczbe
‘
sukces´ow przyj-
muja
‘
c warto´sci n
= 0, 1, . . . N. Jej rozkład prawdopodobie´nstwa to
P
(X = n) =
N
n
q
n
(1 − q)
N
−n
.
(13.10)
Warto´s´c ´srednia
hXi = Nq, natomiast dyspersja
σ
2
X
= Nq(1 − q).
• Rozkład Poissona.
Zmienna losowa X przyjmuje warto´sci n
= 0, 1, 2, . . . . Jej rozkład praw-
dopodobie´nstwa jest dany przez
P
(X = n) =
µ
n
n!
e
−µ
,
µ
> 0 .
(13.11)
Warto´s´c ´srednia
hXi = µ, a dyspersja
σ
2
X
= µ. Rozkład Poissona mo˙zna
otrzyma´c z rozkładu dwumianowego w granicy: N
→
∞
, q → 0 takiej, ˙ze
Nq
= µ = const. Opisuje on wie
‘
c liczbe
‘
sukces´ow w bardzo du˙zej pr´obie,
gdy prawdopodobie´nstwo pojedynczego sukcesu jest małe.
82
13.3
Generowanie liczb losowych
Generowanie liczb losowych o zadanym rozkładzie prawdopodbie´nstwa to cen-
tralny element metod Monte Carlo. Wykorzystuje sie
‘
do tego celu generatory liczb
losowych o rozkładzie jednostajnym na odcinku
[0, 1], dla kt´orch ge
‘
sto´s´c praw-
dopodobie´nstwa jest zadana przez wz´or
p
(x) =
1
x
∈ [0,1]
0
poza
[0, 1] .
(13.12)
Generatory takie sa
‘
oferowane w ramach bibliotek program ´ow komputero-
wych, na przykład CERNLIB. W praktyce generowane liczby nie sa
‘
w pełni lo-
sowe. Na przykład, powtarzaja
‘
sie
‘
w ramach do´s´c długiego cyklu. Moga
‘
tak˙ze
istnie´c korelacje mie
‘
dzy nimi.
Tym niemniej dobry generator powinien spełnia´c trzy podstawowe kryteria.
• Mie´c bardzo długi okres powtarzalno´sci.
Na przykład dla komputera 32-bitowego powinien mie´c okres bliski
2
31
− 1 = 2 147 483 647,
gdy˙z zakres liczb całkowitych dla tego komputera to
[−2
31
, 2
31
− 1]
• Charakteryzowa´c sie
‘
dobra
‘
losowo´scia
‘
. Oznacza to, ˙ze korelacje pomie
‘
dzy
wszystkimi, kolejno generowanymi liczbami powinny by´c mo˙zliwie małe.
• By´c szybki.
Przykładem dobrego generatora jest generator, w kt´orym kolejne liczby losowe
otrzymuje sie
‘
z relacji
x
n
+1
= (ax
n
+ b) mod c .
(13.13)
Liczby a
, b, c nazywa sie
‘
magicznymi, od ich wyboru zale˙zy jako´s´c generatora.
Jednym z dobrych zestaw ´ow tych liczb jest a
= 7
5
= 16 807, b = 0, c = 2
31
− 1.
Przyklad
Przyjmijmy a
= 3, b = 0, c = 2
3
− 1 = 7. Zaczynaja
‘
c od 1 otrzymujemy
kolejne liczby, generowane zgodnie z reguła
‘
x
n
+1
= 3x
n
mod 7:
1 3 2 6 4 5 1
. . .
Dziela
‘
c je przez 7 znajdujemy sze´s´c liczb z przedziału
[0, 1].
Dysponuja
‘
c liczbami o rozkładzie jednostajnym mo˙zemy otrzyma´c liczby lo-
sowe o dowolnych innych rozkładach, np. gausowskim.
83
13.4
Całkowanie metoda¸ Monte Carlo
Podstawowy wz´or na oblicznie całki w metodzie Monte Carlo to
b
Z
a
f
(x) dx ≈ V h f i ± V
s
h f
2
i − h f i
2
N
,
(13.14)
gdzie V
= b − a oraz
h f i =
1
N
N
∑
i
=1
f
(x
i
) ,
f
2
=
1
N
N
∑
i
=1
f
2
(x
i
) .
(13.15)
W powy˙zszych wzorach x
i
sa
‘
liczbami losowymi o rozkładzie jednostajnym na
odcinku
[a, b], natomiast N jest liczba
‘
wygenerowanych liczb losowych. Stała
V wynika z konieczno´sci normalizacji rozkładu jednostajnego do jedynki, gdy˙z
przymuja
‘
c f
≡ 1 powinni´smy otrzyma´c
b
Z
a
dx
= b − a.
(13.16)
Drugie wyra˙zenie po prawej stronie (13.14) jest estymacja
‘
błe
‘
du wyniku
Wz´or (13.14) jest słuszny tak˙ze dla całkowania w n-wymiarowej przestrzeni z
punktami x
= (x
1
, x
2
, . . . , x
n
),
Z
W
f
(x) dx ≈ V h f i ± V
s
h f
2
i − h f i
2
N
.
(13.17)
We wzorze tym V jest n-wymiarowa
‘
obje
‘
to´scia
‘
przestrzeni W , po kt´orej całkujemy
V
=
Z
W
dx
1
dx
2
. . . dx
n
,
(13.18)
natomiast
h f i =
1
N
N
∑
i
=1
f
(x
i
) ,
f
2
=
1
N
N
∑
i
=1
f
2
(x
i
) ,
(13.19)
gdzie tym razem x
i
= (x
1i
, x
2i
, . . . , x
ni
) to wielowymiarowe liczby losowe o rozkła-
dzie jednostajnym w obszarze całkowania.
84
W praktyce obszar W mo˙ze by´c na tyle skomplikowany, ˙ze nie jest łatwo
pr´obkowa´c go przy pomocy rozkładu jednostajnego. Wybieramy wtedy obszar
U zawieraja
‘
cy W ,
W
⊂ U
(13.20)
w kt´orym łatwo wygenerowa´c liczby losowe o rozkładzie jednostajnym. Dodat-
kowo definiujemy nowa
‘
funkcje
‘
˜
f r´owna
‘
funkcji f w obszarze W i zero poza nim.
Wtedy
Z
U
˜f(x)dx =
Z
W
f
(x) dx .
(13.21)
Wyb´or obszaru U wpływa na wielko´s´c błe
‘
du w metodzie Monte Carlo. Bła
‘
d jest
tym wie
‘
kszy im wie
‘
ksza jest r´o˙znica mie
‘
dzy obszarami U i W . Dobrze jest wie
‘
c
wybra´c U tak, by liczba wylosowanych punkt´ow poza obszarem W była jak naj-
mniejsza.
Przyklad
Opisana
‘
metoda
‘
mo˙zemy policzy´c całke
‘
R
b
a
f
(x) dx z dodatniej funkcji f ,
definiuja
‘
c funkcje
‘
dw ´och zmiennych
F
(x, y) =
1
x
≤ y = f (x)
0
x
> y
(13.22)
okre´slona
‘
w kwadracie U
= [a, b] × [0,max( f )]. Wtedy
b
Z
a
f
(x) dx =
Z
U
F
(x, y) dx dy ≈ V
1
N
∑
x
i
≤y
i
1
,
(13.23)
gdzie V jest polem kwadratu U , po kt´orym całkujemy. Warto´s´c całki jest
wie
‘
c proporcjonalna do stosunku liczby punkt´ow le˙zacych pod wykresem
funkcji y
= f (x) do całkowitej liczby N wygenerowanych punkt´ow o rozkła-
dzie jednostajnym w obszarze całkowania.
85
Wyszukiwarka
Podobne podstrony:
K Golec Biernat Metody Numeryczne dla Fizyków
Notka, Informatyka, Informatyka, Informatyka. Metody numeryczne, Kosma Z - Metody numeryczne dla zas
Metody numeryczne dla inżynierów (notatki do wykładów)
Errata, Informatyka, Informatyka, Informatyka. Metody numeryczne, Kosma Z - Metody numeryczne dla za
5 b, Informatyka, Informatyka, Informatyka. Metody numeryczne, Kosma Z - Metody numeryczne dla zasto
4 l pom, Informatyka, Informatyka, Informatyka. Metody numeryczne, Kosma Z - Metody numeryczne dla z
7 e, Informatyka, Informatyka, Informatyka. Metody numeryczne, Kosma Z - Metody numeryczne dla zasto
Nowoczesne metody antykoncepcji dla kobiet i mezczyzn
Metody numeryczne w6
metoda siecznych, Elektrotechnika, SEM3, Metody numeryczne, egzamin metody numeryczn
MN energetyka zadania od wykładowcy 09-05-14, STARE, Metody Numeryczne, Część wykładowa Sem IV
METODA BAIRSTOWA, Politechnika, Lab. Metody numeryczne
testMNłatwy0708, WI ZUT studia, Metody numeryczne, Metody Numeryczne - Ćwiczenia
Metodyka galwanizacji, Fizjoterapia, Fizyko
Metody numeryczne Metoda węzłowa
Metody numeryczne, wstep
metody numeryczne w4
więcej podobnych podstron