Metody Numeryczne II 1
Te materia艂y mo偶na znalez膰 pod adresem:
http://www.phys.uni.torun.pl/~kgrabcze/zajecia/MetNum2.pdf
Bibliografia
" J. Stoer, R. Bulirsch. Wst臋p do analizy numerycznej. PWN 1987.
" G. Dahlquist, A. Bj鰎ck. Metody numeryczne. PWN 1983.
" A. Ralston. Wst臋p do analizy numerycznej. PWN 1971.
" B. P. Demidowicz, I. A. Maron, E. Z. Szuwa艂owa. Metody Numeryczne. PWN 1965.
" Z. Fortuna, B. Macukow, J. W膮sowski. Metody Numeryczne. WNT 1993.
" E. K膮cki, A. Ma艂olepszy, A. Romanowicz. Metody numeryczne dla in偶ynier贸w. WPA
2000.
" J. Krupka, R. Z. Morawski, L. J. Opalski. Metody numeryczne. OWPW 1997.
" A. Smoluk. Metody numeryczne. WAE 1996.
" W. H. Press i in. Numerical Recipes in C (the Art of Scientific Computing).
http://www.phys.uni.torun.pl/nrbook/bookcpdf.html
Metody Numeryczne II 2
Plan
1. Metody numeryczne. B艂臋dy. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2. Interpolacja. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3. Ca艂kowanie numeryczne. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4. Aproksymacja. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
5. Rozwi膮zywanie r贸wna艅 nieliniowych. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
6. Rozwi膮zywanie uk艂ad贸w r贸wna艅 liniowych. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
7. Metody ortogonalizacji. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 204
8. Wyznaczanie warto艣ci i wektor贸w w艂asnych macierzy. . . . . . . . . . . . . . . . . . . . 212
Metody Numeryczne II 3
Metody numeryczne
" Cel oszacowanie dok艂adno艣ci wynik贸w oblicze艅.
" B艂臋dy ograniczaj膮 dok艂adno艣膰.
" Rodzaje b艂臋d贸w:
w danych wej艣ciowych (przyrz膮d贸w pomiarowych itp.),
zaokr膮gle艅 (ograniczenie dok艂adno艣ci oblicze艅),
aproksymacji (uproszczenie problemu, szacowanie a nie wyliczanie),
Przyk艂ad 1:
1 1 1
e = 1 + + + + . . . (1)
1! 2! 3!
mo偶na szacowa膰 licz膮c sko艅czone sumy (powstaje b艂膮d obci臋cia).
Przyk艂ad 2:
ca艂k臋 oznaczon膮 liczymy jako sko艅czon膮 sum臋 p贸l obszar贸w (b艂膮d
dyskretyzacji).
Metody Numeryczne II 4
Reprezentacja liczb
" Systemy liczenia: dw贸jkowy (binarny), 贸semkowy, dziesi臋tny, szesnastkowy.
system cyfry przyk艂ad dziesi臋tnie
dw贸jkowy 0 1 101 1 " 22 + 1 " 20 =
1 " 4 + 1 = 5
贸semkowy 0 1 2 3 4 5 6 7 270 2 " 82 + 7 " 81 =
2 " 64 + 7 " 8 = 184
dziesi臋tny 0 1 2 3 4 309 3 " 102 + 9 " 100 =
5 6 7 8 9 3 " 100 + 9 = 309
szesnastkowy 0 1 2 3 4 5 6 7 10A 1 " 162 + 10 " 160 =
8 9 A B C D E F 1 " 256 + 10 = 266
Metody Numeryczne II 5
Reprezentacja liczb - c.d.
" Reprezentacja sta艂opozycyjna
w praktyce tylko dla liczb ca艂kowitych
" Reprezentacja zmiennopozycyjna (zmiennoprzecinkowa)
mantysa i cecha, liczba = mantysa " 2cecha
18.5 = 0.185 " 102 = 0.100101 " 2101
rozdzielczo艣膰 i sko艅czono艣膰 komputerowych liczb rzeczywistych
typy 4-bajtowe (3 bajty mantysy i 1 bajt cechy), 6-bajtowe, 8-bajtowe
(podw贸jnej precyzji)
nale偶y by膰 ostro偶nym przy dodawaniu ma艂ej liczby zmiennoprzecinkowej
do du偶ej (ma艂a mo偶e znikn膮膰)
warto艣ci, kt贸re nie s膮 poprawnymi liczbami (NAN - not a number)
standard IEEE 754 dobry opis pod
http://research.microsoft.com/~hollasch/cgindex/coding/ieeefloat.html
Metody Numeryczne II 6
Standard IEEE 754
Reprezentacja bitowa
precyzja znak wyk艂adnik cz臋艣膰 u艂amkowa przesuni臋cie
pojedyncza 1 [31] 8 [30-23] 23 [22-00] 127
podw贸jna 1 [63] 11 [62-52] 52 [51-00] 1023
" reprezentacja binarna podstawa r贸wna 2
" znak 0 - liczba dodatnia, 1 - liczba ujemna
" wyk艂adnik i przesuni臋cie rzeczywista cecha = wyk艂adnik - przesuni臋cie
" mantysa
posta膰 znormalizowana liczby w mantysie przecinek po pierwszej
niezerowej cyfrze (1 m < 10)
Przyk艂ad: 1.25 " 102, nie 0.125 " 103 ani 125 " 100
w systemie dw贸jkowym niezerowa znaczy jedynka
pierwszy bit mantysy zawsze r贸wny 1 ukryty, reszta to cz臋艣膰 u艂amkowa
Metody Numeryczne II 7
Standard IEEE 754 przyk艂ady
23.75
" posta膰 binarna: 10111.11 = 24 + 22 + 21 + 20 + 2-1 + 2-2
" po normalizacji: 1.011111 " 24 = 1.011111 " 2100
" znak 0
wyk艂adnik z przesuni臋ciem 4 + 127 = 131 = 10000011
cz臋艣膰 u艂amkowa 011111
" ostatecznie mamy 0 10000011 01111100000000000000000
czyli 01000001 10111110 00000000 00000000
szesnastkowo 0x41BE0000
-0.7
" posta膰 binarna: -0.1011001100110011...
" znak = 0, wyk艂adnik = -1, 127 - 1 = 126 = 01111110
" ostatecznie mamy 1 01111110 01100110011001100110011
czyli 10111111 00110011 00110011 00110011
szesnastkowo 0xBF333333
Metody Numeryczne II 8
Standard IEEE 754 c.d.
Warto艣ci specjalne wyk艂adniki z samych jedynek b膮dz z samych zer
znak wyk艂adnik (w) cz. u艂amkowa (u) warto艣膰
0 00..00 00..00 +0
0 00..00 00..01 11..11 dodatnia zdenorm. 0.u " 2-p+1
0 00..01 11..10 00..00 11..11 dodatnia znorm. 1.u " 2w-p
0 11..11 00..00 +Infinity (niesko艅czono艣膰)
0 11..11 00..01 01..11 SNaN (Signalling Not A Number)
0 11..11 10..00 11..11 QNaN (Quiet Not A Number)
1 00..00 00..00 -0
1 00..00 00..01 11..11 ujemna zdenorm. -0.u " 2-p+1
1 00..01 11..10 00..00 11..11 ujemna znorm. -1.u " 2w-p
1 11..11 00..00 -Infinity (niesko艅czono艣膰)
1 11..11 00..01 01..11 SNaN
1 11..11 10..00 11..11 QNaN
Metody Numeryczne II 9
Standard IEEE 754 c.d.
Operacje specjalne:
Operacja Wynik Operacja Wynik
n / 膮Infinity 0 膮0 / 膮0 NaN
膮Infinity x 膮Infinity 膮Infinity Infinity - Infinity NaN
膮nonzero / 0 膮Infinity 膮Infinity / 膮Infinity NaN
Infinity + Infinity Infinity 膮Infinity x 0 NaN
Zakresy dw贸jkowo:
precyzja zdenormalizowane znormalizowane
pojedyncza 膮2-149 . . . (1 - 2-23) " 2-126 膮2-126 . . . (2 - 2-23) " 2127
podw贸jna 膮2-1074 . . . (1 - 2-52) " 2-1022 膮2-1022 . . . (2 - 2-52) " 21023
Zakresy dziesi臋tnie:
precyzja dziesi臋tnie
pojedyncza 膮 <" 10-44.85 . . . <" 1038.53
podw贸jna 膮 <" 10-323.3 . . . <" 10308.3
Metody Numeryczne II 10
B艂膮d bezwzgl臋dny i wzgl臋dny
Niech x b臋dzie oszacowaniem warto艣ci x (dok艂adnej)
" B艂膮d bezwzgl臋dny
"x = |x - x| (2)
" B艂膮d wzgl臋dny
x - x
x = | | (3)
x
Liczby maszynowe
" zbi贸r liczb maszynowych A - zbi贸r liczb reprezentowalnych w danym formacie
" aproksymacja liczby x liczb膮 maszynow膮 rd(x):
"g"A|x - rd(x)| |x - g| (4)
Metody Numeryczne II 11
B艂臋dy zaokr膮gle艅
" zaokr膮glenie (o ile jest liczb膮 maszynow膮) jest dobr膮 aproksymacj膮 w
systemie o k cyfrach znacz膮cych (gdzie precyzja mantysy to k cyfr) mantys臋
m = 膮1.膮2膮3 . . . zaokr膮glamy do
B
膮1.膮2 . . . 膮k je偶eli 膮k+1 <
2
m = (5)
B
膮1.膮2 . . . 膮k + B-k+1 je偶eli 膮k+1
2
a liczb臋 znormalizowan膮 x = m " Bc do rd(x) = sgn(x)m " Bc.
B艂膮d wzgl臋dny przybli偶enia:
B
" B-k
rd(x) - x 1
2
| | " B-k+1 (6)
x |m| 2
def
1
" dok艂adno艣膰 maszynowa eps = " B-k+1
2
Zatem: rd(x) = x(1 + %E艂), gdzie |%E艂| eps
Dla system贸w binarnych eps = 2-k
Metody Numeryczne II 12
B艂臋dy zaokr膮gle艅 c.d.
" nadmiar i niedomiar cechy nieprawid艂owo艣膰 zatrzymanie oblicze艅, wyj膮tek
na szcz臋艣cie bardzo rzadko si臋 zdarzaj膮
" wyniki dzia艂a艅 arytmetycznych nie musz膮 by膰 liczbami maszynowymi mamy
wi臋c dzia艂ania zmiennopozycyjne aproksymuj膮ce w艂a艣ciwe dzia艂ania, np:
def
x +" y = rd(x + y)
def
x -" y = rd(x - y)
def
x " y = rd(x y)
def
x/"y = rd(x/y)
W贸wczas:
x +" y = (x + y)(1 + %E艂1)
x -" y = (x - y)(1 + %E艂2)
x " y = (x y)(1 + %E艂3)
x/"y = (x/y)(1 + %E艂4),
gdzie |%E艂i| eps
Metody Numeryczne II 13
Problemy arytmetyki zmiennopozycyjnej
eps
" Je偶eli |y| < |x|, to x +" y = x
B
" Dzia艂ania zmiennopozycyjne nie musz膮 by膰 艂膮czne lub rozdzielne.
Przyk艂ad:
a = 2.337125810 - 4, b = 3.3678429102, c = -3.3677811102
a + b + c = 6.4137125810 - 3
a +" (b +" c) = 2.337125810 - 4 +" 6.180000010 - 3 = 6.413712610 - 3
(a +" b) +" c = 3.3678452102 +" -3.3677811102 = 6.410000010 - 3
Oznaczenia:
" fl(wyra偶enie) warto艣膰 wyra偶enia w arytmetyce zmiennopozycyjnej
" dzia艂ania elementarne dzia艂ania arytmetyczne (podstawowe plus pewne
dodatkowe takie jak pierwiastek, sinus itp.)
Metody Numeryczne II 14
Poj臋cie algorytmu
Zadanie: Mamy dan膮 funkcj臋 膯 : D Rm, gdzie D 膮" Rn, n, m " N. Szukamy
metody wyznaczania warto艣ci y = 膯(x).
Odwzorowaniem elementarnym nazywamy odwzorowanie : D 膮" Rk Rl
takie, 偶e jego funkcje sk艂adowe i s膮 dzia艂aniami elementarnymi.
Algorytmem realizuj膮cym funkcj臋 膯 : D Rm (D 膮" Rn, n, m " N) nazywamy
i
taki ci膮g odwzorowa艅 膯(1), . . . , 膯(r), (gdzie 膯(i) : Di Di+1, Di 膮" Rn , D1 = D,
nr+1 = m), 偶e
膯 = 膯(1) 膰% 膰% 膯(r) (7)
Przyk艂ad:
Zadanie 膯(a, b, c) = a + b + c mo偶e by膰 realizowane algorytmem {膯(1), 膯(2)}, gdzie
a + b
膯(1)(a, b, c) = , 膯(2)(d, e) = d + e. (8)
c
Metody Numeryczne II 15
Propagacja b艂臋d贸w
Prze艣ledzmy b艂臋dy zaokr膮gle艅 dla przyk艂adu sumowania (a + b) + c:
fl((a + b) + c) = fl(fl(a + b) + c) =
((a + b)(1 + %E艂1) + c)(1 + %E艂2) =
a + b
(a + b + c) 1 + %E艂1(1 + %E艂2) + %E艂2
a + b + c
a+b
jest wsp贸艂czynnikiem wzmocnienia b艂臋du %E艂1 jego du偶a warto艣膰 powoduje
a+b+c
du偶y b艂膮d wzgl臋dny
Poniewa偶 we wcze艣niejszym przyk艂adzie:
a + b b + c
H" 5 " 104, H" 0.97 (9)
a + b + c a + b + c
to wybra膰 nale偶y metod臋 fl(a + (b + c)).
Metody Numeryczne II 16
Zadanie: Analizujemy naturalny algorytm znajduj膮cy przybli偶enie bn
rozwi膮zania n r贸wnania liniowego
c - a1b1 - - an-1bn-1 - ann = 0, (10)
dla danych c, a1, . . . , an, , b1, . . . , bn-1, an = 0. Oszacowa膰 residuum
r = c - a1b1 - - anbn. (11)
Algorytm: Wyznaczamy
c - a1b1 - - an-1bn-1
bn = fl (12)
an
poprzez:
s0 = c,
sj = fl(sj-1 - ajbj) = [sj-1 - ajbj(1 + 礿)] (1 + 膮j), j = 1, 2, . . . , n - 1,
sn-1 sn-1
bn = fl = (1 + ) .
an an
(13)
Metody Numeryczne II 17
Oszacowanie 1: Z (13) wynikaj膮 (drugie r贸wnanie dla j = 1, 2, . . . , n - 1):
s0 - c = 0,
sj 膮j
sj - (sj-1 - ajbj) = sj - + ajbj礿 = sj - ajbj礿, (14)
1 + 膮j 1 + 膮j
anbn - sn-1 = sn-1.
Dodaj膮c stronami dostajemy:
n n-1
膮j
r = c - aibi = -sj + ajbj礿 - sn-1. (15)
1 + 膮j
i=1 j=1
St膮d oszacowanie:
肱 雠
n-1
0 je偶eli an = 1
eps
砼偞2 |sn-1| +
|r| (|sj| + |ajbj|)艂艂 , 2 =
1 - eps
1 w przec. przyp.
j=1
(16)
Metody Numeryczne II 18
Oszacowanie 2: Z (13):
钆 艂艂
n-1 n-1 n-1
1 +
鹋俢 (1 + 膮k) -
bn = ajbj(1 + 礿) (1 + 膮k) . (17)
an
k=1 j=1 k=j
St膮d
j-1
n-1 n-1
c = ajbj(1 + 礿) (1 + 膮k)-1 + anbn(1 + )-1 (1 + 膮k)-1 (18)
j=1 k=1 k=1
Aatwo pokaza膰 indukcyjnie (ze wzgl臋du na m), 偶e z
m
1 + = (1 + k)膮1, |k| eps, m eps < 1 (19)
k=1
wynika
m eps
|| . (20)
1 - m eps
Metody Numeryczne II 19
Zatem istniej膮 wielko艣ci j takie, 偶e
n-1
c = ajbj(1 + jj) + anbn(1 + (n - 1 + 2 )n),
j=1
(21)
0 je偶eli an = 1|
eps
|j| , 2 =
1 - n eps
1 w przec. przyp.
Ostatecznie dla r = c - a1b1 - - anbn dostajemy:
钆 艂艂
n-1
eps
鹋
|r| j|ajbj| + (n - 1 + 2 )|anbn| . (22)
1 - n eps
j=1
Metody Numeryczne II 20
Wskazniki uwarunkowania
钆 艂艂
膯1
锱 艣艂
.
.
Niech 膯 : D Rm dla otwartego zbioru D 膮" Rn, oraz 膯 = przy czym
鹋
.
膯m
膯i maj膮 ci膮g艂e pochodne na D. Przy za艂o偶eniu pomiar贸w niedok艂adnych x
wielko艣ci x otrzymujemy 艂 = 膯(x) jako przybli偶enie y = 膯(x). Rozwini臋cie w
szereg Taylora z pomini臋ciem wyraz贸w wy偶szych rz臋d贸w daje:
n n
膯i(x) 膯i(x)
.
"yi = 艂 i - yi = 膯i(x) - 膯i(x) = (xj - xj) = "xj (23)
xj xj
j=1 j=1
W zwi膮zku z tym b艂膮d wzgl臋dny:
n
"yi . xj 膯i(x) "xj
= (24)
yi 膯i(x) xj xj
j=1
xj 膯i(x)
Wsp贸艂czynniki nazywamy wskaznikami uwarunkowania zadania.
膯i(x) xj
Metody Numeryczne II 21
Zadania zle i dobrze uwarunkowane
Zadanie nazywamy zle uwarunkowanym je偶eli przynajmniej jeden ze wskaznik贸w
uwarunkowania ma du偶膮 warto艣膰 bezwzgl臋dn膮. W przeciwnym razie zadanie
nazywamy dobrze uwarunkowanym.
" tylko dla niezerowych xj i yi
" mn wskaznik贸w
" inne definicje, np: c jest wskaznikiem uwarunkowania zadania 膯 je艣li
||膯(x) - 膯(x)|| ||x - x||
c (25)
||膯(x)|| ||x||
Metody Numeryczne II 22
B艂臋dy nieuniknione i numeryczna stabilno艣膰
B艂臋dem nieuniknionym zwi膮zanym z y = 膯(x) nazywamy optymalny poziom
b艂臋du wyniku y czyli b艂膮d wynikaj膮cy z zaokr膮gle艅 danych wej艣ciowych (niezale偶ny
od algorytmu). Analiza b艂臋d贸w w kontek艣cie algorytm贸w jako z艂o偶enia
odwzorowa艅 elementarnych prowadzi do oszacowania:
"(0)y = [|D膯(x)| |x| + |y|]eps. (26)
Je艣li b艂臋dy zaokr膮gle艅 algorytmu s膮 w przybli偶eniu r贸wne "(0)y, to m贸wimy, 偶e s膮
one nieszkodliwe, a algorytm nazywamy stabilnym (poprawnym) numerycznie
lub m贸wimy, 偶e zachowuje si臋 dobrze.
Algorytm jest numerycznie u偶yteczniejszy od innego (wyznaczaj膮cego t臋 sam膮
wielko艣膰) dla danego zbioru, gdy ca艂kowity b艂膮d zaokr膮glenia jest mniejszy w
pierwszym ni偶 w drugim.
Uwaga: Stabilno艣膰 numeryczna i wi臋ksza u偶yteczno艣膰 numeryczna s膮
w艂asno艣ciami lokalnymi (zale偶nymi od warto艣ci wej艣cia).
Metody Numeryczne II 23
Arytmetyka przedzia艂owa
" Cel: wyznaczanie g贸rnych ogranicze艅 b艂臋du bezwzgl臋dnego
" Uwzgl臋dnia b艂臋dy zaokr膮gle艅 i b艂臋dy danych
" Obliczenia na przedzia艂ach mo偶liwych warto艣ci liczby zamiast na liczbach: dla
a " R mamy przedzia艂 膰 = [a2 , a2 2 ], gdzie a2 , a2 2 " A.
" Dzia艂ania na przedzia艂ach: dla " " {",
", ", 艢"} wynik
c = 膰 " b " {x y|x " 膰, y " b} (27)
Przyk艂ady:
[a2 , a2 2 ] " [b2 , b2 2 ] = [max{x " A|x a2 + b2 }, min{x " A|x a2 2 + b2 2 }]
[a2 , a2 2 ] " [b2 , b2 2 ] = [max{x " A|x a2 b2 }, min{x " A|x a2 2 b2 2 }]
Uwaga: Oszacowania przedzia艂owe s膮 prawdziwe, ale do艣膰 pesymistyczne
Przyk艂ad: Schemat Hornera a wzory skr贸conego mno偶enia dla
x3 - 3x2 + 3x = (x - 1)3 + 1 (28)
Metody Numeryczne II 24
Arytmetyka przedzia艂owa a statystyka
Pesymizm arytmetyki przedzia艂owej wida膰 ju偶 na przyk艂adzie sumy n liczb.
Niech n = 2 i 膯(x, y) = z = x + y.
Za艂贸偶my, 偶e x = 艂 = [-1, 1].
W贸wczas z = [-2, 2]. Ale rozk艂ad prawdopodobie艅stwa z nie jest jednostajny.
arytm.przedz.
n = 2
n = 10
n = 30
1/n
-n n
Dla n = 30, arytmetyka przedzia艂owa daje zakres oko艂o dwukrotnie szerszy ni偶
przedzia艂 istotnie niezerowego prawdopodobie艅stwa.
Metody Numeryczne II 25
Statystyczna teoria b艂臋d贸w
" B艂臋dem standardowym warto艣ci przybli偶onej nazywamy odchylenie
standardowe tej warto艣ci traktowanej jako zmienna losowa.
" B艂臋dem prawdopodobnym nazywamy tak膮 liczb臋 x, 偶e z
1
prawdopodobie艅stwem mamy |%E艂| < x.
2
Dla rozk艂adu normalnego mamy wi臋c x = 0.6745, gdzie jest
odchyleniem standardowym tego rozk艂adu.
n
W poprzednim przyk艂adzie mamy = , a wi臋c:
3
n = 1 ! H" 0.577
n = 10 ! H" 1.82 = 0.182n
n = 30 ! H" 3.16 H" 0.1n
Metody Numeryczne II 26
Twierdzenie 1 (o b艂臋dzie standardowym) Je偶eli b艂臋dy "x1, . . . , "xn s膮
niezale偶nymi zmiennymi losowymi o warto艣ci oczekiwanej r贸wnej zeru i
odchyleniach standardowych odpowiednio %E艂1, . . . , %E艂n, to b艂膮d standardowy dla
y = y(x1, . . . , xn) dany jest wzorem:
2
n
y
%E艂 H" %E艂2. (29)
i
xi
i=1
Uwaga: Prosta suma wielu zmiennych jest przypadkiem szczeg贸lnym powy偶szego
twierdzenia (pochodne cz膮stkowe s膮 r贸wne 1).
Metody Numeryczne II 27
Zadanie interpolacji
Definicja 1 Niech 艢(x; a0, . . . , an) b臋dzie funkcj膮 jednej zmiennej z n + 1
parametrami. Warto艣ci parametr贸w okre艣laj膮 funkcje rodziny 艢. Zadanie
interpolacyjne dla 艢 przy zadanych m + 1 parach liczb (xi, fi), i = 0, . . . , m,
gdzie xi = xj dla i = j, polega na wyznaczeniu parametr贸w ai tak, aby
艢(xi; a0, . . . , an) = fi, i = 0, . . . , m. (30)
Pary (xi, fi) nazywamy punktami w臋z艂owymi a warto艣ci xi i fi odpowiednio
odci臋t膮 i rz臋dn膮 w臋z艂a.
Metody Numeryczne II 28
Zadania interpolacji
" i. wielomianowa
艢(x; a0, . . . , an) = a0 + a1x + . . . + anxn (31)
" i. trygonometryczna
艢(x; a0, . . . , an) = a0 + a1exi + . . . + anenxi (32)
" i. funkcjami sklejanymi
" i. wymierna
a0 + a1x + . . . + akxk
艢(x; a0, . . . , ak, b0, . . . , bl) = (33)
b0 + b1x + . . . + blxl
" i. eksponencjalna
0 1 n
艢(x; a0, . . . , an; 0, . . . , n) = a0e x + a1e x + . . . + ane x (34)
Metody Numeryczne II 29
Wz贸r interpolacyjny Lagrange a
Oznaczenie: n zbi贸r wszystkich wielomian贸w stopnia nie wiekszego ni偶 n
(postaci P (x) = a0 + a1x + . . . + anxn)
Twierdzenie 2 Dla n+1 dowolnych punkt贸w w臋z艂owych (xi, fi), i = 0, . . . , n,
takich, 偶e xi = xj dla i = j, istnieje tylko jeden wielomian p " n spe艂niaj膮cy
P (xi) = fi, i = 0, . . . , n (35)
Dow贸d:
Jednoznaczno艣ci: Niech P1, P2 " n spe艂niaj膮 warunki (35). W贸wczas
P = P1 - P2 " n jest wielomianem o n + 1 r贸偶nych zerach, a wi臋c jest
to偶samo艣ciowo r贸wny zeru.
Metody Numeryczne II 30
Istnienia: Takim wielomianem jest:
n
P (x) = fiLi(x) (36)
i=0
gdzie Li s膮 tak zwanymi wielomianami Lagrange a zdefiniowanymi jako:
n
x - xj
Li(x) = , (37)
j=0
xi - xj
(j=i)
i spe艂niaj膮cymi nast臋puj膮ce warunki (delty Kroneckera):
1, gdy i = j
Li(xj) = ij = . (38)
0, gdy i = j
Uwagi:
" Niedu偶a przydatno艣膰 w praktyce (n2 - 1 mno偶e艅, 1 dzielenie).
" Metoda przydatna w贸wczas, gdy rozwa偶amy wiele zada艅 interpolacji dla
w臋z艂贸w o tych samych odci臋tych.
Metody Numeryczne II 31
Algorytm Neville a
Idea: Rozwi膮zanie w krokach - od pojedynczych w臋z艂贸w do ca艂ego ich zbioru.
Oznaczenie: Dla danego zbioru punkt贸w w臋z艂owych (xi, fi), i = 0, . . . , n przez
Pi ,...,ik oznaczamy wielomian ze zbioru k taki, 偶e
0
Pi ...ik(xi ) = fi , j = 0, . . . , k. (39)
0 j j
Twierdzenie 3 Wielomiany typu (39) maj膮 nast臋puj膮ce w艂asno艣ci:
1膰% Pi(x) = fi
(x - xi )Pi ...ik(x) - (x - xi )Pi ...ik-1(x)
0 1 k 0
2膰% Pi ...ik(x) =
0
xi - xi
k 0
Dow贸d: W艂asno艣膰 1膰% jest oczywista.
Uwzgl臋dniaj膮c definicj臋 (39) 艂atwo zauwa偶y膰, 偶e prawa strona w艂asno艣ci 2膰% jest dla
x = xi r贸wna fi oraz dla x = xi r贸wna fi .
0 0 k k
Metody Numeryczne II 32
Dla pozosta艂ych xi (dla j = 1, . . . , k - 1) otrzymujemy z prawej strony
j
(xi - x0)fi - (xi - xi )fi
j j j k j
= fi ,
j
xi - xi
k 0
co ostatecznie dowodzi w艂asno艣ci 2膰%.
Tabela ilustruj膮ca algorytm Neville a:
k = 0 1 2 3
x0 f0 = P0(x)
P01(x)
x f1 = P1(x) P012(x)
1
P12(x) P0123(x)
x f2 = P2(x) P123(x)
2
P23(x)
x f3 = P3(x)
3
Uwagi:
" sprawnie wyznacza warto艣ci funkcji dla pojedynczych punkt贸w
" gorzej z wyznaczeniem samego wielomianu interpolacyjnego
Metody Numeryczne II 33
Algorytm Neville a wersja 2
Oznaczenie:
def
Ti+k k = Pi...i+k (40)
Tablica ilustruj膮ca algorytm:
x0 f0 = T00(x)
T11(x)
x1 f1 = T10(x) T22(x)
T21(x) T33(x)
x2 f2 = T20(x) T32(x)
T31(x)
x3 f3 = T30(x)
1膰% Ti0(x) = fi
(x - xi-j)Ti j-1(x) - (x - xi)Ti-1 j-1(x) Ti j-1 - Ti-1 j-1
2膰% Tij(x) = = Ti j-1 +
x-xi-j
xi - xi-j
- 1
x-xi
Metody Numeryczne II 34
Wz贸r interpolacyjny Newtona
Wielomian interpolacyjny P " n, P (xi) = fi dla i = 1, . . . , n mo偶na przedstawi膰
w postaci:
P0...n(x) = a0+a1(x-x0)+a2(x-x0)(x-x1)+ +an(x-x0) (x-xn-1). (41)
Warto艣ci P (x) mo偶na wyznacza膰 na wz贸r schematu Hornera (n - 1 mno偶e艅).
Wsp贸艂czynniki ai mo偶na wyznaczy膰 wprost z:
f0 = P (x0) = a0
f1 = P (x1) = a0 + a1(x1 - x0)
f2 = P (x2) = a0 + a1(x2 - x0) + a2(x2 - x0)(x2 - x1)
. . .
wykonuj膮c n dziele艅 i n(n - 1) mno偶e艅.
Metody Numeryczne II 35
Ale mo偶na te偶 za pomoc膮 n(n + 1)/2 dziele艅. Zauwa偶my, 偶e dla Pi ...ik oraz
0
Pi ...ik-1 istnieje jednoznacznie wsp贸艂czynnik fi ...ik (iloraz r贸偶nicowy k-tego
0 0
rz臋du) taki, 偶e:
Pi ...ik(x) - Pi ...ik-1(x) = fi ...ik(x - xi ) (x - xi ). (42)
0 0 0 0 k-1
Mamy wi臋c nast臋puj膮c膮 posta膰 Newtona cz臋艣ciowego wielomianu Pi ...ik(x):
0
fi +fi i1(x-x0)+fi i1i2(x-x0)(x-x1)+ +fi ...ik(x-x0) (x-xi ). (43)
0 0 0 0 k-1
Twierdzenie 4 Ilorazy r贸偶nicowe fi ...ik s膮 niezale偶ne od permutacji
0
indeks贸w.
Dow贸d: fi ...ik jest wsp贸艂czynnikiem przy najwy偶szej pot臋dze wielomianu Pi ...ik,
0 0
kt贸ry jest jednoznacznie wyznaczony przez w臋z艂y, kt贸re interpoluje.
Metody Numeryczne II 36
Korzystaj膮c z w艂asno艣ci 2膰% z twierdzenia 3 oraz z (42) mamy:
Pi (x) - Pi ...ik-1(x)
(42) ...ik 2膰%
0 0
fi ...ik = =
0
(x - xi ) (x - xi )
0 k-1
(x-xi0 )Pi1...ik (x)-(x-xik )Pi0...ik-1 (x)
- Pi ...ik-1(x)
0
xik -xi0
=
(x - xi ) (x - xi )
0 k-1
(x - xi )Pi ...ik(x) - (x - xi )Pi ...ik-1(x) - (xi - xi )Pi ...ik-1(x)
0 1 k 0 k 0 0
=
(xi - xi )(x - xi ) (x - xi )
k 0 0 k-1
(x - xi )(Pi ...ik(x) - Pi ...ik-1(x))
0 1 0
=
(xi - xi )(x - xi ) (x - xi )
k 0 0 k-1
(Pi ...ik(x) - Pi ...ik-1(x) + Pi ...ik-1(x) - Pi ...ik-1(x)) fi ...ik - fi ...ik-1
(42)
1 1 1 0 1 0
=
(xi - xi )(x - xi ) (x - xi ) xi - xi
k 0 1 k-1 k 0
Metody Numeryczne II 37
Schemat iloraz贸w r贸偶nicowych (tablica ilustruj膮ca metod臋 Newtona ilorazy
liczone na wz贸r metody Neville a):
k = 0 1 2 3
x0 f0
f01
x1 f1 f012
f12 f0123
x2 f2 f123
f23
x3 f3
Wielomian interpolacyjny:
P (x) = f0+f01(x-x0)+f012(x-x0)(x-x1)+ +f0...n(x-x0) (x-xn-1) (44)
Uwagi:
" Dla zminimalizowania b艂臋du przy obliczaniu warto艣ci wielomianu
interpolacyjnego Newtona mo偶na ustali膰 odpowiedni膮 kolejno艣膰 punkt贸w
Metody Numeryczne II 38
w臋z艂owych. Dla danego x dobieramy kolejno艣膰 i0, . . . , in tak, aby
|x - xi | |x - xi | dla k = 1, . . . , n.
k k-1
" Je偶eli wszystkie odci臋te xi w臋z艂贸w s膮 uporz膮dkowane tak, 偶e xi < xi+1, to ze
schematu iloraz贸w r贸偶nicowych mo偶na odczyta膰 wszystkie reprezentacje
Newtona. Dla danego x wyb贸r indeks贸w nast臋puje tak, by ka偶dy nast臋pny by艂
s膮siedni do kt贸rego艣 z poprzednich.
Przyk艂ad: Zb贸r w臋z艂贸w: {(0, 1), (1, 2), (3, 10)}.
P210(x) = 10 + 4(x - 3) + 1(x - 3)(x - 1)
x0 = 0 1
P210(x) = (1(x - 1) + 4)(x - 3) + 10
1
x1 = 1 2 1
P210(2) = 5
4
P120(x) = 2 + 4(x - 1) + 1(x - 1)(x - 3)
x2 = 3 10
P (x) = (1(x - 3) + 4)(x - 1) + 2
120
Metody Numeryczne II 39
Reszta w interpolacji wielomianowej
Twierdzenie 5 Je偶eli funkcja f jest n + 1-krotnie r贸偶niczkowalna, to dla
ka偶dego argumentu x istnieje liczba w najmniejszym przedziale
I[x0, . . . , xn, x] zawieraj膮cym x i odci臋te wszystkich w臋z艂贸w, spe艂niaj膮ca
(x)f(n+1)()
f(x) - P01...n(x) = , (45)
(n + 1)!
gdzie
(x) = (x - x0)(x - x1) (x - xn). (46)
Uwagi:
" Interpolacja wielomianowa nie nadaje si臋 do ekstrapolacji
" Brak zbie偶no艣ci do warto艣ci funkcji interpolowanej przy liczbie punkt贸w
w臋z艂owych rosn膮cej do niesko艅czono艣ci i malej膮cej do zera 艣rednicy
podzia艂贸w.
Metody Numeryczne II 40
Interpolacja wymierna
Zadanie interpolacji w臋z艂贸w (xi, fi) funkcj膮 wymiern膮:
,
P (x) a0 + a1x + + a祒
艢,(x) = = . (47)
Q,(x) b0 + b1x + + bx
" Par臋 (, ) nazywamy stopniem zadania interpolacji wymiernej.
" + + 2 parametry minus czynnik wsp贸lny
" + + 1 w臋z艂贸w (xi, fi), i = 0, . . . , +
Rozwi膮zanie zadania musi by膰 rozwi膮zaniem uk艂adu r贸wna艅 S,:
,
P (xi) - fiQ,(xi) = 0 (48)
czyli
a0 + a1xi + + a祒 - fi(b0 + b1xi + + bx) = 0 (49)
i
i
Metody Numeryczne II 41
Przyk艂ad: Zbi贸r punkt贸w w臋z艂owych {(0, 1), (1, 2), (2, 2)} dla = = 1 prowadzi
do uk艂adu r贸wna艅:
a0 - b0 = 0
a0 + a1 - 2(b0 + b1) = 0
a0 + 2a1 - 2(b0 + 2b1) = 0
Rozwi膮zanie uk艂adu a0 = b0 = 0, a1 = 2, b1 = 1 nie spe艂nia warunk贸w zadania
2x
funkcja wymierna 艢1,1(x) = nie jest okre艣lona w punkcie x = 0. Zatem
x
rozwi膮zanie zadania nie istnieje.
Metody Numeryczne II 42
Definicja 2 Wyra偶enia wymierne:
P1(x) P2(x)
艢1(x) = 艢2(x) = Q1 a" 0, Q2 a" 0
Q1(x) Q2(x)
nazywamy r贸wnowa偶nymi gdy przedstawiaj膮 t臋 sam膮 funkcj臋 wymiern膮 tj.
wtedy, gdy
P1(x)Q2(x) = P2(x)Q1(x)
Wyra偶enie wymierne nazywamy relatywnie pierwszym je艣li jego licznik i
mianownik s膮 relatywnie pierwsze, czyli nie s膮 podzielne przez ten sam
wielomian dodatniego stopnia.
Metody Numeryczne II 43
Twierdzenie 6 Uk艂ad S, ma zawsze rozwi膮zanie nietrywialne. Dla ka偶dego
takiego rozwi膮zania Q, a" 0.
Dow贸d: S, jest jednorodnym uk艂adem + + 1 r贸wna艅 z + + 2
niewiadomymi, wi臋c ma rozwi膮zania nietrywialne (a0, . . . , a, b0, . . . , b).
Gdyby Q, a" 0, czyli " bi = 0, to mieliby艣my
i=0
,
P (xi) = 0, dla i = 0, . . . , + + 1.
, ,
Poniewa偶 P jest stopnia < + + 1, wi臋c mieliby艣my P a" 0, co jest
sprzeczne z za艂o偶eniem.
Metody Numeryczne II 44
Twierdzenie 7 (Jednoznaczno艣膰 rozwi膮zania nietrywialnego) Je偶eli 艢1
oraz 艢2 s膮 dwoma nietrywialnymi rozwi膮zaniami uk艂adu S,, to s膮 one
r贸wnowa偶ne.
P1(x) P2(x)
Dow贸d: Je艣li 艢1(x) = oraz 艢2(x) = , to
Q1(x) Q2(x)
P (x) = P1(x)Q2(x) - P2(x)Q1(x) (50)
ma + + 1 r贸偶nych zer (xi, dla i = 0, . . . , + + 1). Poniewa偶 stopie艅 P jest
niewi臋kszy ni偶 + , to P a" 0, czyli 艢1 i 艢2 s膮 r贸wnowa偶ne.
Wniosek 8 (z twierdze艅 6 i 7) Dla ka偶dego zadania interpolacji wymiernej
A, istnieje jednoznaczna funkcja wymierna 艢,, kt贸ra rozwi膮zuje
odpowiadaj膮cy mu uk艂ad r贸wna艅 S,. Albo funkcja 艢, jest rozwi膮zaniem
zadania A,, albo zadanie to nie jest rozwi膮zalne.
Metody Numeryczne II 45
Oznaczenia:
" Punkty nieosi膮galne zadania A, to w臋z艂y, kt贸re nie s膮 spe艂nione przez
funkcj臋 wymiern膮 艢,.
ozn
" 艢, = wyra偶enie wymierne relatywnie pierwsze, r贸wnowa偶ne 艢,.
Twierdzenie 9 Zadanie A, jest rozwi膮zalne, wtedy i tylko wtedy, gdy dla
ka偶dego rozwi膮zania 艢, uk艂adu S,, 艢, rozwi膮zuje S,.
Dow贸d:
!: Za艂贸偶my, 偶e 艢, nie rozwi膮zuje S,. W贸wczas funkcja wymierna
odpowiadaj膮ca wyra偶eniu 艢, nie rozwi膮zuje zadania A,. St膮d i z
wniosku 8, zadanie A, jest nierozwi膮zalne.
P (x)
!: Niech 艢,(x) = . Rozwa偶my dwa przypadki:
Q(x)
1. "i"{0,...,+}Q(xi) = 0. W贸wczas 艢,(xi) = fi.
2. "i"{0,...,+}Q(xi) = 0. W贸wczas w臋ze艂 (xi, fi) by艂by nieosi膮galny, oraz
mieliby艣my P (xi) = 0. Zatem P i Q zawiera艂yby czynnik x - xi i nie by艂yby
relatywnie pierwsze sprzeczno艣膰!
Metody Numeryczne II 46
Uwagi:
" Zadanie nierozwiazalne mo偶na zredukowa膰 do rozwi膮zalnego poprzez
odpowiednie zmniejszenie liczby w臋z艂贸w.
" Dla zada艅 niedegenerowalnych (tj. takich, 偶e ka偶dy podzbi贸r w臋z艂贸w stanowi
zadanie rozwi膮zalne) Am,n mo偶na rozpatrywa膰 ci膮g wyra偶e艅 wymiernych 艢,
stopni (, ), gdzie m, n.
" W kolejnym kroku zwi臋kszamy stopie艅 licznika b膮dz mianownika.
Metody Numeryczne II 47
Odwrotno艣ci iloraz贸w r贸偶nicowych
xi - xj
def
膯(xi, xj) =
fi - fj
xj - xk
def
膯(xi, xj, xk) =
膯(xi, xj) - 膯(xi, xk)
(51)
.
.
.
xm - xn
def
膯(xi, . . . , xl, xm, xn) =
膯(xi, . . . , xl, xm) - 膯(xi, . . . , xl, xn)
Uwaga: Nie s膮 symetrycznymi funkcjami swoich argument贸w.
0 1 2 3 . . .
0
Metoda wyznacza wyra偶enia wy-
1
mierne odpowiadaj膮ce g艂贸wnej
2
przek膮tnej:
3
.
.
.
Metody Numeryczne II 48
Za艂贸偶my, 偶e wyznaczamy funkcj臋 wymiern膮
n
P (x)
艢n,n = (52)
Qn(x)
tak膮, 偶e 艢n,n(xi) = fi dla i = 0, . . . , 2n. Mamy:
n
P (x0)
n
n n n
P (x) - Qn(x)
P (x) P (x) P (x0)
Qn(x0)
= f0 + - = f0 + =
Qn(x) Qn(x) Qn(x0) Qn(x)
(53)
n-1
P (x) x - x0
f0 + (x - x0) = f0 +
n-1
Qn(x) Qn(x)/P (x)
Dla i = 1, . . . , 2n zachodzi wi臋c:
Qn(xi) xi - x0
= = 膯(x0, xi) (54)
n-1
P (xi) fi - f0
I dalej:
Qn(x) Qn(x) Qn(x1)
= 膯(x0, x1) + - =
n-1 n-1 n-1
P (x) P (x) P (x1)
(55)
x - x1
膯(x0, x1) +
n-1
P (x)/Qn-1(x)
Metody Numeryczne II 49
Kontynuuj膮c, otrzymamy wyra偶enie reprezentowane przez u艂amek 艂a艅cuchowy:
艢n,n(x) = f0 + x - x0/膯(x0, x1) + x - x1/膯(x0, x1, x2) +
(56)
+ x - x2n-1/膯(x0, x1, . . . , x2n)
Uwaga: Kolejne u艂amki cz臋艣ciowe s膮 wyra偶eniami wymiernymi 艢,(x) oraz
艢+1,(x) dla = 0, 1, . . . , n - 1 spe艂niaj膮cymi odpowiednie podzbiory zbioru
w臋z艂贸w zadania:
艢0,0(x) = f0
艢1,0(x) = f0 + x - x0/膯(x0, x1)
(57)
艢1,1(x) = f0 + x - x0/膯(x0, x1) + x - x1/膯(x0, x1, x2)
. . .
Algorytm: Realizacja wzoru 56 z pomoc膮 tablicy odwrotno艣ci iloraz贸w:
0 x0 f0
1 x f 膯(x , x )
1 1 0 1
2 x2 f2 膯(x0, x2) 膯(x0, x1, x2)
3 x3 f3 膯(x0, x3) 膯(x0, x1, x3) 膯(x0, x1, x2, x3)
. . .
Metody Numeryczne II 50
Odwrotne ilorazy r贸偶nicowe
def
(xi) = fi
xi - xi+1
def
(xi, xi+1) =
fi - fi+1
.
.
.
xi - xi+k
def
(xi, . . . , xi+k) = + (xi+1, . . . , xi+k-1)
(xi, . . . , xi+k-1) - (xi+1, . . . , xi+k)
(58)
Uwaga: S膮 symetrycznymi funkcjami swoich argument贸w.
Twierdzenie 10 Dla p = 1, 2, . . ., przyjmuj膮c (x0, . . . , xp-2) = 0 dla p = 1
spe艂nione jest:
膯(x0, . . . , xp) = (x0, . . . , xp) - (x0, . . . , xp-2) (59)
Dow贸d: Indukcja. Dla p = 1 teza jest oczywista.
Metody Numeryczne II 51
Zak艂adamy prawdziwo艣膰 dla p. W贸wczas:
xp - xp+1
膯(x0, . . . , xp+1) =
膯(x0, . . . , xp) - 膯(x0, . . . , xp-1, xp+1)
(60)
xp - xp+1
=
(x0, . . . , xp) - (x0, . . . , xp-1, xp+1)
Z definicji mamy:
xp+1 - xp
(xp+1, x0, . . . , xp)-(x0, . . . , xp-1) = (61)
(xp+1, x0, . . . , xp-1) - (x0, . . . , xp)
Wobec symetryczno艣ci mamy tez臋.
U艂amek 艂a艅cuchowy Thiele a otrzymujemy zast臋puj膮c odwrotno艣ci iloraz贸w
r贸偶nicowych odwrotnymi ilorazami r贸偶nicowymi:
艢n,n(x) = f0 + x - x0/(x0, x1) + x - x1/(x0, x1, x2) - (x0) +
(62)
+ x - x2n-1/(x0, . . . , x2n) - (x0, . . . , x2n-2)
Metody Numeryczne II 52
Uog贸lnienie algorytmu Neville a
Oznaczenia:
,
Ps (x)
def
" 艢,(x) = wyra偶enie wymierne spe艂niaj膮ce warunki: 艢,(xi) = fi
s s
Q,(x)
s
,
dla i = s, s + 1, . . . , s + + . Ps (x) i Q,(x) s膮 wielomianami stopni
s
odpowiednio i .
def
" 膮i = x - xi
Twierdzenie 11 Dla 1 oraz 1 zachodz膮:
艢-1,(x) - 艢-1,(x)
s+1 s
(a) 艢,(x) = 艢-1,(x) +
s s+1
艢-1,(x) - 艢-1,(x)
膮s
s+1 s
1 - - 1
膮s++
艢-1,(x) - 艢-1,-1(x)
s+1 s+1
艢,-1(x) - 艢,-1(x)
s+1 s
(b) 艢,(x) = 艢,-1(x) +
s s+1
艢,-1(x) - 艢,-1(x)
膮s
s+1 s
1 - - 1
膮s++
艢,-1(x) - 艢-1,-1(x)
s+1 s+1
Metody Numeryczne II 53
" Wzor贸w z twierdzenia 11 mo偶na u偶y膰 do
0 1 2 3 . . .
wyznaczenia warto艣ci wyra偶e艅 wymiernych
przy odpowiednio wzrastaj膮cych stopniach.
0
" R贸偶ne 艣cie偶ki to r贸偶ne metody i r贸偶ne final-
1
ne funkcje wymierne.
2
" Dla = 0 i wzrastaj膮cego mamy interpo-
lacj臋 wielomianow膮 Neville a.
3
.
" = 0 i wzrastaj膮ce odpowiadaj膮 interpo-
.
.
lacji wielomianowej dla w臋z艂贸w (xi, 1/fi).
" Do艣wiadczenie pokazuje, 偶e warto realizowa膰 艣cie偶k臋 oznaczon膮 lini膮 ci膮g艂膮.
" W przypadku konkretnej 艣cie偶ki suma + jednoznacznie wyznacza i , co
pozwala upro艣ci膰 oznaczenia.
Metody Numeryczne II 54
Dla 艣cie偶ki oznaczonej ci膮g艂膮 lini膮 mo偶emy wprowadzi膰:
ozn
Tik = 艢,, gdzie i = s + + , k = + (63)
s
Mamy w贸wczas nast臋puj膮cy algorytm:
Ti0 = fi, Ti,-1 = 0
Ti,k-1 - Ti-1,k-1
(64)
Tik = Ti,k-1 +
x - xi-k Ti,k-1 - Ti-1,k-1
1 - - 1
x - xi Ti,k-1 - Ti-1,k-2
Tabela ilustruj膮ca algorytm:
(, ) = (0, 0) (0, 1) (1, 1) (1, 2)
T00 = f0
T0,-1 = 0 T11
T10 = f1 T22
T1,-1 = 0 T21 T33
T20 = f2 T32
T2,-1 = 0 T31
T30 = f3
Metody Numeryczne II 55
Interpolacja wielomianowa a wymierna ctg(2.5膰%) = 22.9037655484
x1 fi = ctg(xi)
1膰% 57.28996163
14.30939911
2膰% 28.63625328 21.47137102
23.85869499 22.36661762
3膰% 19.08113669 23.26186421 22.63519158
21.47137190 23.08281486
4膰% 14.30066626 22.18756808
18.60658719
5膰% 11.43005230
1膰% 57.28996163
22.90760673
2膰% 28.63625328 22.90341624
22.90201805 22.90369573
3膰% 19.08113669 22.90411487 22.90376552
22.91041916 22.90384141
4膰% 14.30066626 22.90201975
22.94418151
5膰% 11.43005230
Metody Numeryczne II 56
Interpolacja funkcjami sklejanymi
def
Definicja 3 Niech " = {x0, x1, . . . , xn}, gdzie a = x0 < x1 < . . . < xn = b
b臋dzie podzia艂em przedzia艂u [a, b]. Funkcj膮 sklejan膮 trzeciego stopnia na
" nazywamy ka偶d膮 funkcj臋 S" : [a, b] R tak膮, 偶e:
" S" " C2[a, b] (dwukrotnie ci膮gle r贸偶niczkowalna na [a, b])
" S" na ka偶dym podprzedziale [xi, xi+1], i = 0, . . . , n - 1 pokrywa si臋 z
wielomianem trzeciego stopnia.
Metody Numeryczne II 57
Oznaczenie: Dla zbioru liczb rzeczywistych Y = {y0, . . . , yn} przez
S"(Y, ) (65)
rozumie膰 b臋dziemy funkcj臋 sklejan膮 S" spe艂niaj膮c膮 warunki S"(Y, xi) = yi dla
i = 0, . . . , n.
Uwaga: S" nie jest wyznaczona jednoznacznie (n funkcji stopnia 3 to 4n
parametr贸w. 2n r贸wna艅 dla punkt贸w w臋z艂owych + 2(n - 1) r贸wno艣ci
pochodnych = 4n - 2 r贸wnania). Pozosta艂e dwa stopnie swobody mo偶na
wyeliminowa膰 np. jednym z nast臋puj膮cych trzech warunk贸w:
2 2 2 2
(a) S"(Y, a) = S"(Y, b) = 0.
(k) (k)
(66)
(b) S" (Y, a) = S" (Y, b) dla k = 0, 1, 2. S"(Y, ) jest okresowa.
2 2 2 2 2 2
(c) S"(Y, a) = y0, S"(Y, b) = yn dla danych y0, yn " R.
Metody Numeryczne II 58
Interpolacja funkcjami sklejanymi
Niech " = {x0, x1, . . . , xn} b臋dzie ustalonym podzia艂em przedzia艂u [a, b] oraz
Y = {y0, y1, . . . , yn} " R
Oznaczenie:
ozn
hj+1 = xj+1 - xj, j = 0, 1, . . . , n - 1 (67)
Definicja 4 Momentami funkcji S"(Y, ) nazywamy warto艣ci jej drugich
pochodnych w w臋z艂ach:
def
2 2
Mj = S"(Y, xj), j = 0, 1, . . . , n (68)
Uwagi:
" Druga pochodna funkcji sklejanej to funkcja kawa艂kami liniowa (liniowa na
przedzia艂ach [xi, xi+1]).
" Funkcje sklejane s膮 jednoznacznie wyznaczone przez swoje momenty.
" Momenty mog膮 by膰 wyznaczone jako rozwi膮zanie uk艂adu r贸wna艅 liniowych.
Metody Numeryczne II 59
Druga pochodna a momenty:
xj+1 - x x - xj
2 2
S"(Y, x) = Mj + Mj+1 dla x " [xj, xj+1] (69)
hj+1 hj+1
Ca艂kuj膮c otrzymujemy
(xj+1 - x)2 (x - xj)2
2
S"(Y, x) = -Mj + Mj+1 + Aj (70)
2hj+1 2hj+1
(xj+1 - x)3 (x - xj)3
S"(Y, x) = Mj + Mj+1 + Aj(x - xj) + Bj (71)
6hj+1 6hj+1
dla x " [xj, xj+1], j = 0, 1, . . . , n - 1, gdzie Aj oraz Bj oznaczaj膮 sta艂e ca艂kowania.
Ze spe艂niania odpowiednich punkt贸w w臋z艂owych otrzymujemy:
h2
Mj j+1 + Bj = yj (72)
6
h2
Mj+1 j+1 + Ajhj+1 + Bj = yj+1 (73)
6
Metody Numeryczne II 60
St膮d:
h2
Bj = yj - Mj j+1 (74)
6
h2
j+1
yj+1 - Mj+1 6 - Bj
yj+1 - yj hj+1
Aj = = - (Mj+1 - Mj) (75)
hj+1 hj+1 6
Mamy wi臋c dla x " [xj, xj + 1]:
S"(Y, x) = 膮j + j(x - xj) + 艂j(x - xj)2 + j(x - xj)3 (76)
gdzie:
2 2
S"(Y, xj) Mj
膮j = S"(Y, xj) = yj, 艂j = = (77)
2 2
Mjhj+1 yj+1 - yj hj+1
2
j = S"(Y, xj) = - + Aj = - (Mj+1 + 2Mj) (78)
2 hj+1 6
2 2 2
S" (Y, x+)
Mj+1 - Mj
j
j = = (79)
6 6hj+1
Metody Numeryczne II 61
Wyznaczanie moment贸w
2
Za艂o偶enie ci膮g艂o艣ci S"(Y, ) w w臋z艂ach wyznacza n - 1 r贸wna艅 dla moment贸w.
Podstawiaj膮c Aj z (75) do (70) otrzymujemy:
(xj+1 - x)2 (x - xj)2 yj+1 - yj hj+1
2
S"(Y, x) = -Mj + Mj+1 + - (Mj+1 - Mj)
2hj+1 2hj+1 hj+1 6
(80)
W szczeg贸lno艣ci dla j = 1, 2, . . . , n - 1 mamy:
yj - yj-1 hj hj
2
S"(Y, x-) = + Mj + Mj-1 (81)
j
hj 3 6
yj+1 - yj hj+1 hj+1
2
S"(Y, x+) = - Mj - Mj+1 (82)
j
hj+1 3 6
A z r贸wno艣ci granic:
hj hj + hj+1 hj+1 yj+1 - yj yj - yj-1
Mj-1 + Mj + Mj+1 = - (83)
6 3 6 hj+1 hj
czyli mamy n - 1 r贸wna艅 dla n + 1 moment贸w.
Metody Numeryczne II 62
Wyznaczanie moment贸w - c.d.
Pozosta艂e 2 r贸wnania z dodatkowych warunk贸w (66).
Przypadek (66.a):
2 2 2 2
S"(Y, a) = M0 = 0 S"(Y, b) = Mn = 0 (84)
Przypadek (66.b): funkcja okresowa, czyli hn+1 = h1, Mn+1 = M1, yn+1 = y1:
2 2 2 2
S"(Y, a) = S"(Y, b) a" M0 = Mn (85)
hn hn + h1 h1 y1 - yn yn - yn-1
2 2
S"(Y, a) = S"(Y, b) a" Mn-1+ Mn+ M1 = -
6 3 6 h1 hn
(86)
Przypadek (66.c):
h1 h1 y1 - y0
2 2 2
S"(Y, a) = y0 a" M0 + M1 = - y0 (87)
3 6 h1
Metody Numeryczne II 63
hn hn yn - yn-1
2 2 2
S"(Y, b) = yn a" Mn-1 + Mn = yn - (88)
6 3 hn
Og贸lna posta膰 dla r贸wna艅 (83), (86), (87), (88):
礿Mj-1 + 2Mj + jMj+1 = dj, j = 1, 2, . . . , n - 1, (89)
gdzie:
hj+1 hj
j = , 礿 = 1 - j = ,
hj + hj+1 hj + hj+1
(90)
6 yj+1 - yj yj - yj-1
dj = -
hj + hj+1 hj+1 hj
Dodatkowo:
Przypadek (66.a):
0 = 0, d0 = 0, 祅 = 0, dn = 0 (91)
Metody Numeryczne II 64
Przypadek (66.c):
6 y1 - y0
2
0 = 1, d0 = - y0
h1 h1
(92)
6 yn - yn-1
2
祅 = 1, dn = yn -
hn hn
Zatem dla przypadk贸w (66.a) i (66.c) otrzymujemy uk艂ad r贸wna艅:
钆 艂艂 钆 艂艂 钆 艂艂
2 0 0 M0 d0
锱 艣艂 锱
1 2 1 M1 艣艂 锱 d1 艣艂
锱 艣艂 锱 艣艂 锱 艣艂
锱 艣艂 锱 艣艂 锱 艣艂
. .
.
. . .
= (93)
锱 艣艂 锱 艣艂 锱 艣艂
.
. .
锱 艣艂 锱 艣艂 锱 艣艂
鹋
祅-1 2 n-1 鹋 Mn-1 鹋 dn-1
祅 2 Mn dn
Metody Numeryczne II 65
Przypadek (66.b):
h1 hn
n = , 祅 = 1 - n =
hn + h1 hn + h1
(94)
6 y1 - yn yn - yn-1
dn = -
hn + h1 h1 hn
Otrzymujemy nast臋puj膮cy uk艂ad r贸wna艅:
钆 艂艂 钆 艂艂 钆 艂艂
2 1 1 M1 d1
锱 艣艂 锱
2 2 2 M2 艣艂 锱 d2 艣艂
锱 艣艂 锱 艣艂 锱 艣艂
锱 艣艂 锱 艣艂 锱 艣艂
. .
.
. . .
= (95)
锱 艣艂 锱 艣艂 锱 艣艂
.
. .
锱 艣艂 锱 艣艂 锱 艣艂
鹋
祅-1 2 n-1 鹋 Mn-1 鹋 dn-1
n 祅 2 Mn dn
oraz M0 = Mn
Metody Numeryczne II 66
Interpolacja funkcjami sklejanymi - podsumowanie
Twierdzenie 12 Uk艂ady r贸wna艅 liniowych 93 i 95 s膮 nieosobliwe dla
dowolnego podzia艂u " przedzia艂u [a, b].
Schemat dowodu: Najpierw dowodzimy w艂asno艣ci:
Az = w ! max |zi| max |wi|.
i i
Nast臋pnie pokazujemy, 偶e osobliwo艣膰 macierzy przeczy tej w艂asno艣ci.
Uwagi:
" Uk艂ady te rozwi膮zujemy metod膮 eliminacji Gaussa (uproszczon膮, bo macierze
tr贸jdiagonalne).
" Bardzo wa偶na w艂asno艣膰! W przeciwie艅stwie do wielomian贸w, interpolowane
funkcje sklejane zbiegaj膮 (przy pewnych s艂abych za艂o偶eniach) do funkcji
interpolowanej.
Metody Numeryczne II 67
Metody ca艂kowania
Definicja
b
2
f(x)dx = F (b) - F (a), gdzie F (x) = f(x) (96)
a
Ca艂kowanie symboliczne i liczenie warto艣ci wprost z definicji cz臋sto jest nierealne
" trudno艣膰 wyznaczenia funkcji pierwotnej
" brak mo偶liwo艣ci wyznaczenia funkcji pierwotnej (np. nie znamy do ko艅ca
samej funkcji f(x))
" wyznaczanie ca艂ek za pomoc膮 dyskretyzacji
Metody Numeryczne II 68
Kwadratury Newtona-Cotesa
" Zast臋pujemy funkcj臋 podca艂kow膮 wielomianem interpolacyjnym.
" Przyjmijmy r贸wnomierny podzia艂 przedzia艂u ca艂kowania [a.b]:
b - a
xi = a + ih, i = 0, . . . , n, gdzie h = , n > 0, n " N (97)
n
Niech Pn b臋dzie wielomianem interpolacyjnym stopnia n takim, 偶e:
Pn(xi) = fi = f(xi) dla i = 0, . . . , n (98)
Ze wzoru interpolacyjnego Lagrange a mamy:
n n
x - xk
Pn(x) = fiLi(x), gdzie Li(x) = (99)
xi - xk
k=0
i=0
k =i
Metody Numeryczne II 69
" Wprowadzamy zmienn膮 t tak膮, 偶e x = a + ht. Mamy wi臋c:
n
t - k
Li(x) = i(t) = (100)
i - k
k=0
k =i
" Ca艂kujemy:
n n n
b b n
Pn(x)dx = fi Li(x)dx = h fi i(t)dt = h fi膮i (101)
a a 0
i=0 i=0 i=0
gdzie wagi
n
膮i = i(t)dt (102)
0
nie zale偶膮 od ca艂kowanej funkcji ani od przedzia艂u ca艂kowania.
Metody Numeryczne II 70
" Dla n = 2 otrzymujemy tzw. kwadratur臋 Simpsona:
2 2
(t - 1)(t - 2) 1
膮0 = dt = (t2 - 3t + 2)dt =
(0 - 1)(0 - 2) 2
0 0
2
1 1 3 1 8 1
= t3 - t2 + 2t = - 6 + 4 =
2 3 2 2 3 3
0
(103)
2 2
(t - 0)(t - 2) 4
膮1 = dt = - (t2 - 2t)dt =
(1 - 0)(1 - 2) 3
0 0
2 2
(t - 0)(t - 1) 1 1
膮2 = dt = (t2 - t)dt =
(2 - 0)(2 - 1) 2 3
0 0
Ostatecznie:
b b
b - a 1
f(x)dx H" Pn(x)dx = (f0 + 4f1 + f2) (104)
2 3
a a
Metody Numeryczne II 71
" w og贸lnej postaci: kwadratury Newtona-Cotesa
n
b
b - a
Pn(x)dx = h fi膮i, fi = f(a + ih), h = (105)
n
a
i=0
n
Zachodzi: 膮i = n
i=0
ozn
Niech s b臋dzie wsp贸lnym mianownikiem wag 膮i. W贸wczas i = s膮i s膮
ca艂kowite oraz:
n n
b
b - a
Pn(x)dx = h fi膮i = ifi (106)
ns
a
i=0 i=0
" reszta kwadratury tj. b艂膮d aproksymacji:
b b
Pn(x)dx - f(x)dx = hp+1Kf(p)(), gdzie " (a, b) (107)
a a
Metody Numeryczne II 72
Zestawienie parametr贸w kwadratur Newtona-Cotesa:
n i ns Reszta Kwadratura
1 1 1 2 h3 1 f(2)() Trapez贸w
12
1
2 1 4 1 6 h f(4)() Simpsona
90
3 1 3 3 1 8 h5 3 f(4)() 3/8
80
4 7 32 12 32 7 90 h7 8 f(6)() Milne a
945
5 19 75 50 50 75 19 288 h7 275 f(6)()
12096
6 41 216 27 272 27 216 41 840 h9 9 f(8)() Weddle a
1400
Dla n > 6 pojawiaj膮 si臋 ujemne wsp贸艂czynniki co stanowi niebezpiecze艅stwo
du偶ych b艂臋d贸w numerycznych.
Uwagi:
def
" rz膮d kwadratury = r gdy ca艂kowanie wielomian贸w stopnia r - 1 jest
bezb艂臋dne
" kwadratura trapez贸w jest rz臋du 2, Simpsona rz臋du 4
Metody Numeryczne II 73
Kwadratury z艂o偶one
b-a
" Ca艂kujemy w [xi, xi+1] dla i = 0, . . . , n - 1, gdzie xi = a + ih, h =
n
" Z艂o偶ona kwadratura trapez贸w
n-1 n-2
1 f(a) + f(b)
def
T (h) = h [f(xi) + f(xi+1)] = h + f(a + ih) (108)
2 2
i=0 i=1
Dla f " C2[a, b] reszta kwadratury w przedziale [xi, xi+1] wynosi
1
h3f(2)(i), gdzie i " [xi, xi+1] (109)
12
Sumuj膮c reszty otrzymujemy:
b
h3 n-1
T (h) - f(x)dx = f(2)(i) =
12
a
i=0
(110)
n-1
b - a 1 b - a
= h2 f(2)(i) = h2f(2)(), gdzie " (a, b)
12 n 12
i=0
Metody Numeryczne II 74
" Z艂o偶ona kwadratura Simpsona (n parzyste)
def
1
S(h) = h[f(a) + 4f(a + h) + 2f(a + 2h) + 4f(a + 3h) + . . .
3
(111)
. . . + 2f(b - 2h) + 4f(b - h) + f(b)]
Reszta kwadratury:
n/2-1
b
h5
S(h) - f(x)dx = f(4)(i) =
90
a
i=0
n/2-1
h4 b - a 2 b - a
= f(4)(i) = h4f(4)(), gdzie " (a, b)
90 2 n 180
i=0
(112)
Uwagi:
" rz膮d kwadratury z艂o偶onej pozostaje taki jak sk艂adanej kwadratury
" reszty zbiegaj膮 do zera z pot臋g膮 h (mimo sumowania si臋 b艂臋d贸w sk艂adanych
kwadratur)
Metody Numeryczne II 75
Z艂o偶ona kwadratura Simpsona w C++:
1 typedef double FunCalk(double);
2
3 double SimpsonZ(FunCalk *f, double a, double b, int n)
4 {
5 if (n%2 || n<4) throw "Niepoprawne n";
6 double h = (b-a)/n;
7 double parts[3];
8 parts[0] = parts[1] = 0;
9 parts[2] = f(a) + f(b);
10 double z = b-h;
11 for (int i=n,k=0; --i; z-=h,k=!k)
12 parts[k] += f(z);
13 return (parts[2]+2*parts[1]+4*parts[0])*h/3;
14 }
Metody Numeryczne II 76
Kwadratury z interpolacji Hermite a
Interpolacja Hermite a
(k)
" zagadnienie interpolacji w臋z艂贸w (xi, yi ), k = 0, . . . , ni - 1, i = 0, . . . , m
m
(k)
(k)
wielomianem P stopnia n = ni - 1 spe艂niaj膮cego P (xi) = yi .
i=0
" rozwi膮zanie jest jednoznaczne
" uog贸lnienie wielomian贸w Lagrange a:
m ni-1
1 gdy i = j '" k =
(k)
P (x) = yi Lik(x), gdzie L()(xj) = (113)
ik
0 w przec. przyp.
i=0 k=0
Metody Numeryczne II 77
Kwadratura dla w臋z艂贸w f(0), f(1), f (0), f (1)
Z uog贸lnionego wzoru Lagrange a dostajemy:
P (t) = f(0)[(t-1)2 +2t(t-1)2]+f(1)[t2 -2t2(t-1)]+f2 (0)t(t-1)2 +f2 (1)t2(t-1)
(114)
Po sca艂kowaniu:
1
1 1
P (t)dt = (f(0) + f(1)) + (f2 (0) - f2 (1)) (115)
2 12
0
Po przekszta艂ceniu zmiennych:
b
1 1
f(x)dx H" M(h) = h(f(a) + f(b)) + h2(f2 (a) - f2 (b)) (116)
2 12
a
gdzie h = b - a.
Reszta kwadratury:
b
1
M(h) - f(x)dx = - h5f(4)(), " (a, b) (117)
720
a
Metody Numeryczne II 78
Wz贸r sumacyjny Eulera-Maclaurina
Twierdzenie 13 Niech g " C2m+2[0, 1], m " N. W贸wczas:
m
1
1 B2l
g(t)dt = (g(0) + g(1)) + (g(2l-1)(0) - g(2l-1)(1))
2 (2l)!
0
l=1 (118)
B2m+2
- g(2m+2)(), " (0, 1)
(2m + 2)!
1 1 1 1
gdzie Bk to liczby Bernoulliego (B2 = , B4 = - , B6 = , B8 = - , . . . ).
6 30 42 30
Schemat dowodu
" Wielomiany i liczby Bernoulliego definicja i w艂asno艣ci.
1
" Stosowne rozpisanie ca艂ki g(t)dt.
0
Metody Numeryczne II 79
Wielomiany i liczby Bernoulliego
Definicja 5 Wielomianami Bernoulliego nazywamy rodzin臋 wielomian贸w
Bi, i = 0, 1, . . . spe艂niaj膮c膮 nast臋puj膮ce warunki:
1
1. B1(x) = x -
2
2
2. Bk+1 = (k + 1)Bk(x), k = 0, 1, . . .
3. B2l+1(0) = B2l+1(1) = 0, l = 1, 2, . . .
1 1
Przyklady:
B0(x) = 1, B1(x) = x - , B2(x) = x2 - x + ,
2 6
(119)
3 1 1
B3(x) = x3 - x2 + x, B4(x) = x4 - 2x3 + x2 -
2 2 30
Definicja 6 Liczby Bernoulliego to warto艣ci Bk = Bk(0), k = 0, 1, . . ..
Metody Numeryczne II 80
Wielomiany Bernoulliego
1
0.8
0.6
0.4
W艂asno艣ci:
0.2
0
" Wielomian Bk jest stopnia k.
-0.2
" (-1)kBk(1 - x) = Bk(x)
-0.4
" Bk(0) = Bk(1) = Bk dla k = 1
-0.6
0.2 0.4 0.6 0.8 1
Metody Numeryczne II 81
Wykorzystanie wielomian贸w Bernoulliego dla dowodu wzoru (118)
Ca艂kowanie przez cz臋艣ci:
1 1
g(t)dt = B1(t)g(t)|1 - B1(t)g2 (t)dt
0
0 0
1
1 1
1 1
B1(t)g2 (t)dt = B2(t)g2 (t) - B2(t)g2 2 (t)dt
2 2
0 0
0
.
.
. (120)
1
1 1
1 1
Bk-1(t)g(k-1)(t)dt = Bk(t)g(k-1)(t) - Bk(t)g(k)(t)dt
k k
0 0
0
Metody Numeryczne II 82
Wz贸r sumacyjny Eulera-Maclaurina
Rozszerzaj膮c wz贸r (118) na n " N przedzia艂贸w otrzymujemy dla g " C2m+2[0, n]:
n-1 m
n
g(0) + g(n) B2l
g(t)dt = + g(i) + (g(2l-1)(0) - g(2l-1)(n))
2 (2l)!
0
i=1 l=1
B2m+2
- ng(2m+2)(), " (0, n) (121)
(2m + 2)!
W og贸lnym przypadku przedzia艂u [a, b] i jego r贸wnomiernego podzia艂u na n cz臋艣ci
b-a
(h = ):
n
m
b
B2l
f(x)dx = T (h) + h2l (f(2l-1)(a) - f(2l-1)(b))
(2l)!
a
l=1
B2m+2
-h2m+2 (b - a)f(2m+2)(), " (a, b), (122)
(2m + 2)!
gdzie T (h) oznacza z艂o偶on膮 kwadratur臋 trapez贸w.
Metody Numeryczne II 83
Kwadratury ekstrapolacyjne
Rozwini臋cie (122) z艂o偶onego wzoru trapez贸w T (h) dla funkcji f w zale偶no艣ci od
d艂ugo艣ci kroku h:
T (h) = 0 + 1h2 + + mh2m + 膮m+1(h)h2m+2, (123)
gdzie
b
B2k
0 = f(t)dt, k = (f(2k-1)(b) - f(2k-1)(a)), k = 1, . . . , m (124)
(2k)!
a
oraz
B2m+2
膮m+1(h) = (b - a)f(2m+2)((h)), (h) " (a, b) (125)
(2m + 2)!
Metody Numeryczne II 84
Kwadratura Romberga
Je艣li f ma wszystkie pochodne w [a, b], to w przej艣ciu do granicy m "
otrzymujemy szereg pot臋gowy:
0 + 1h2 + 2h4 + (126)
Reszta rozwini臋cia maleje z malej膮cym h. Rozwini臋cie mo偶na traktowa膰 jak
wielomian zmiennej h2 przyjmuj膮cy warto艣膰 ca艂ki 0 w h = 0.
b-a
Dla ci膮gu d艂ugo艣ci krok贸w {hi = : i = 0, 1, . . . , m}, gdzie n0 = 1 oraz
ni
{ni : i = 0, 1, . . . , m} jest ci膮giem 艣ci艣le rosn膮cym, wyznaczamy kwadratury
trapez贸w:
Ti0 = T (hi), i = 0, 1, . . . , m, (127)
kt贸re traktujemy jako w臋z艂y zadania interpolacji.
Je艣li Tmm(h) jest wielomianem zmiennej h2 interpoluj膮cym te w臋z艂y, to warto艣膰
ekstrapolowana Tmm(0) jest przybli偶eniem szukanej ca艂ki 0.
Metody Numeryczne II 85
Interpolacja Neville a dla kwadratury Romberga
Cel: Wyznaczy膰 Tmm(0)
Algorytm: Dla indeks贸w i, k takich, 偶e 1 k i m definiujemy Tik(h) jako
wielomian zmiennej h2 stopnia co najwy偶ej k spe艂niaj膮cy:
Tik(hj) = T (hj) dla j = i - k, . . . , i. (128)
Warto艣膰 Tmm(0) wyznaczamy algorytmem Neville a:
Ti,k-1 - Ti-1,k-1
Tik = Ti,k-1 + . (129)
2
hi-k
- 1
hi
W贸wczas:
Tik = Tik(0) (130)
Metody Numeryczne II 86
Tablica algorytmu Neville a dla kwadratury Romberga:
h2 T (h0) = T00
0
T11
h2 T (h1) = T10 T22
1
T T
21 33
h2 T (h2) = T20 T32
2
T31
h2 T (h ) = T
3 30
3
Uwagi:
" Niekt贸re kwadratury dla i = k to kwadratury Newtona-Cotesa:
h0
h1 = ! T11 jest kwadratur膮 Simpsona
2
h1
dodatkowo h2 = ! T22 jest kwadratur膮 Milne a
2
hi-1
" oryginalna metoda Romberga: hi = , i = 1, 2 . . .
2
h2i-2 h2i-2
" ci膮g Bulirsch a: h2i-1 = , h2i = , i = 1, 2 . . .
2 3
Metody Numeryczne II 87
Interpolacja wymierna dla kwadratury Romberga
Dla interpolacji wielomianowej
" trudno okre艣li膰 dok艂adno艣膰,
" cz臋sto warunek stopu to: |Tm-1,m-1 - Tm,m-1| s, gdzie s jest
b
oszacowaniem warto艣ci ca艂ki |f(t)|dt,
a
W interpolacji wymiernej mamy:
Ti0 = T (hi) (131)
Ti,k-1 - Ti-1,k-1
Tik = Ti,k-1 + , 1 k i m (132)
2
hi-k Ti,k-1-Ti-1,k-1
1 - - 1
hi Ti,k-1-Ti-1,k-2
" mo偶na oszacowa膰 reszty
Metody Numeryczne II 88
Kwadratury Gaussa
Zadanie: Wyznaczanie ca艂ek postaci
b
I(f) = (x)f(x)dx, (133)
a
gdzie (x) jest funkcj膮 wagow膮 na (sko艅czonym lub niesko艅czonym) przedziale
[a, b], spe艂niaj膮c膮 nast臋puj膮ce warunki:
" (x) 0
" jest funkcj膮 mierzaln膮 na przedziale [a, b]
b
" dla k = 0, 1, . . . momenty 祂 = xk(x)dx istniej膮 i s膮 sko艅czone
a
b
" dla wielomian贸w s(x) nieujemnych na [a, b] z r贸wno艣ci (x)s(x)dx = 0
a
wynika, 偶e s a" 0.
Metody Numeryczne II 89
Rozpatrywane kwadratury:
n
I(f) = wif(xi). (134)
i=1
" W przeciwie艅stwie do kwadratur Newtona-Cotesa tutaj w臋z艂y nie musz膮 by膰
r贸wnoodleg艂e.
" Jak dobra膰 odci臋te xi oraz wi, 偶eby zmaksymalizowa膰 rz膮d kwadratury?
Metody Numeryczne II 90
Podstawowe definicje:
Definicja 7 Zbiorem unormowanych wielomian贸w rzeczywistych
stopnia j nazywamy
def
j = {p : p(x) = xj + a1xj-1 + + aj} (135)
Definicja 8 Iloczynem skalarnym funkcji f i g nazywamy
b
def
(f, g) = (x)f(x)g(x)dx. (136)
a
Definicja 9 Przez L2[a, b] oznacza膰 b臋dziemy przestrze艅 wszystkich funkcji
okre艣lonych na przedziale [a, b], dla kt贸rych ca艂ka
b
(f, f) = (x)(f(x))2dx (137)
a
jest dobrze okre艣lona i sko艅czona.
Metody Numeryczne II 91
Twierdzenie 14 (Wielomiany ortogonalne) Istniej膮 wielomiany
pj " j, j = 0, 1, . . . spe艂niaj膮ce
(pi, pk) = 0 dla i = k. (138)
Wielomiany te s膮 wyznaczone jednoznacznie wzorami:
p0(x) a" 1 (139)
2
pi+1(x) a" (x - i+1)pi(x) - 艂i+1pi-1(x), (140)
gdzie p-1(x) = 0 oraz
(xpi, pi)
i+1 = dla i 0 (141)
(pi, pi)
0 dla i = 0
2
艂i+1 = (142)
(pi,pi)
dla i 1
(pi-1,pi-1)
Dow贸d: Ortogonalizacja Grama-Schmidta:
Z definicji 0 mamy p0(x) a" 1.
Za艂贸偶my prawdziwo艣膰 tezy dla wszystkich j i. Dowolny wielomian pi+1 " i+1
Metody Numeryczne II 92
ma jednoznaczne przedstawienie
pi+1(x) = (x - i+1)pi(x) + ci-1pi-1(x) + + c0p0(x). (143)
Je艣li pi+1 " i+1, to
" dla j = i mamy:
(pi+1, pj) = (xpi, pi) - i+1(pi, pi) = 0 (144)
" dla j < i mamy:
(pi+1, pj) = (xpi, pj) + cj(pj, pj) = 0 (145)
Z za艂o偶e艅 dla (x) przy podstawieniu s = p2 wynika, 偶e (pk, pk) = 0, wi臋c
k
i+1 i cj s膮 wyznaczone jednoznacznie.
Z za艂o偶enia indukcyjnego dla j < i mamy:
2
pj+1(x) = (x - j+1)pj(x) - 艂j+1pj-1(x) (146)
Metody Numeryczne II 93
St膮d (pj+1, pi) = (xpj, pi) = (xpi, pj), a zatem
2
(pj+1, pi)
-艂j+2 dla j = i - 1
cj = - = , (147)
0 dla j < i - 1
(pj, pj)
czyli
2
pi+1(x) = (x - i+1)pi(x) - 艂i+1pi-1(x). (148)
Wniosek 15 Je艣li p " n-1, to (p, pn) = 0.
Twierdzenie 16 Wielomian pn ma n rzeczywistych i pojedynczych zer -
wszystkie w przedziale [a, b].
Twierdzenie 17 Ci膮g wielomian贸w ortogonalnych jest uk艂adem Czebyszewa
tzn. spe艂nia nast臋puj膮cy warunek Haara:
Dla wzajemnie r贸偶nych argument贸w ti, i = 1, . . . , n, macierz
钆 艂艂
p0(t1) pn-1(t1)
锱 艣艂
. . .
. . .
A = (149)
鹋
. . .
p0(tn) pn-1(tn)
jest nieosobliwa.
Metody Numeryczne II 94
Twierdzenie 18 Niech x1, . . . , xn b臋d膮 zerami wielomianu pn oraz niech
w1, . . . , wn b臋d膮 rozwi膮zaniem nieosobliwego uk艂adu r贸wna艅
n
(p0, p0) dla k = 0,
pk(xi)wi = (150)
0 dla k = 1, . . . , n - 1.
i=1
W贸wczas wi > 0 dla i = 1, . . . , n, oraz dla ka偶dego wielomianu p " 2n-1
n
b
(x)p(x)dx = wip(xi). (151)
a
i=1
Je艣li wi, xi, i = 1, . . . , n spe艂niaj膮 (151) dla wszystkich p " 2n-1, to xi s膮
zerami wielomianu pn, a wi spe艂niaj膮 (150).
Nie istniej膮 takie xi, wi, i = 1, . . . , n, kt贸re spe艂niaj膮 (151) dla wszystkich
p " 2n.
Metody Numeryczne II 95
Twierdzenie 19 Zera wielomianu ortogonalnego pn s膮 warto艣ciami w艂asnymi
macierzy tr贸jdiagonalnej
钆 艂艂
1 艂2 0
锱 艣艂
艂2 2 艂3
锱 艣艂
Jn = . (152)
锱 艣艂
. .
.
. . .
鹋
.
. .
0 艂n n
Twierdzenie 20 (Reszta kwadratury Gaussa) Je偶eli f " C2n[a, b], to dla
pewnego " (a, b)
n
b
f(2n)()
(x)f(x)dx - wif(xi) = (pn, pn). (153)
(2n)!
a
i=1
Metody Numeryczne II 96
Najprostsza kwadratura Gaussa (Gaussa-Legendre a):
Dla funkcji wagowej a" 1 i przedzia艂u [a, b] = [-1, 1], wielomianami ortogonalymi
s膮:
k! dk
pk(x) = (x2 - 1)k, k = 0, 1, . . . (154)
(2k)! dxk
1
0.5
0
-1 -0.5 0.5 1
-0.5
-1
Metody Numeryczne II 97
Inne kwadratury Gaussa:
[a, b] (x) Wielomiany ortogonalne
1
"
[-1, 1] Tn(x) Czebyszewa
1-x2
[0, "] e-x Ln(x) Laguerre a
2
[- ", "] e-x Hn(x) Hermite a
" n " n
x + x2 - 1 + x - x2 - 1
Tn(x) = cos(n arccos(x)) = (155)
2
1
0.5
0
-1 -0.5 0.5 1
-0.5
-1
Metody Numeryczne II 98
Metody aproksymacji
Zadanie aproksymacji funkcji f(x) polega na wyznaczeniu parametr贸w
a0, . . . , am pewnej funkcji F tak, by zminimalizowa膰 ||f(x) - F (a0, . . . , am; x)||.
" przestrze艅 parametr贸w wyznacza przestrze艅 funkcji
" aproksymacja punktowa i integralna
" najcz臋stszy przypadek: przestrze艅 funkcji jest przestrzeni膮 liniow膮 funkcja F
ma posta膰 wielomianu uog贸lnionego
F (x) = a00(x) + a11(x) + + amm(x) (156)
" aproksymacja wymierna
a00(x) + a11(x) + + akk(x)
F (x) = (157)
b00(x) + b11(x) + + bll(x)
Metody Numeryczne II 99
Wyb贸r przestrzeni funkcji i minimalizowanej normy
" Interpolacja - przestrze艅 gwarantuje zerow膮 warto艣膰 normy
" Aproksymacja 艣redniokwadratowa:
1
b
2
||f||2, = (x)f(x)2dx (158)
a
w przypadku dyskretnym
n
||f||2, = (xi)f(xi)2 (159)
i=0
" Aproksymacja jednostajna
||f|| = supx"[a,b]|f(x)| (160)
Metody Numeryczne II 100
Aproksymacja 艣redniokwadratowa wielomianami uog贸lnionymi
Rozpartujemy funkcje postaci
m
F (x) = aii(x). (161)
i=0
0, . . . , m to funkcje bazowe m + 1 wymiarowej przestrzeni liniowej.
R贸偶ne uk艂ady bazowe:
" 1, x, x2, . . . , xm
" rodziny wielomian贸w z interpolacji Lagrange a czy Newtona
" wielomiany ortogonalne
" sin x, cos x, sin 2x, cos 2x, . . . , sin kx, cos kx
Metody Numeryczne II 101
Zadanie: minimalizacja b艂臋du:
钆 艂艂2
n m
鹋俧(xi) -
E(a0, . . . , am) = (xi) ajj(xi) (162)
i=0 j=0
Minimum E w punkcie, gdzie pochodne cz膮stkowe s膮 zerowe
钆 艂艂
n m
E
鹋俧(xi) -
= - 2(xi) ajj(xi) k(xi) k = 0, . . . , m (163)
ak
i=0 j=0
Uk艂ad normalny w postaci macierzowej:
DT DA = DT f (164)
gdzie
钆 艂艂 钆 艂艂 钆 艂艂
0(x0) m(x0) a0 f(x0)
锱 艣艂 锱 艣艂 锱 艣艂
. . . .
. . . .
D = , A = , f = (165)
鹋 鹋 鹋
. . . .
0(xn) m(xn) am f(xn)
Metody Numeryczne II 102
Aproksymacja wielomianowa
Twierdzenie 21 (Weierstrassa) Je偶eli funkcja f(x) jest ci膮g艂a na
sko艅czonym przedziale [a, b], to dla ka偶dego > 0 mo偶na dobra膰 takie n, 偶e
istnieje wielomian P stopnia n, kt贸ry na ca艂ym przedziale [a,b] spe艂nia
nier贸wno艣膰
|f(x) - P (x)| < . (166)
Funkcje bazowe:
i(x) = xi. (167)
Warunki minimum:
钆 艂艂
n m
E
鹋俧(xi) -
= 0 a" ajxj xk = 0 (168)
i
i
ak
i=0 j=0
czyli
n m n
f(xi)xk = aj xj+k, k = 0, . . . , m (169)
i
i
i=0 j=0 i=0
Metody Numeryczne II 103
W prostszej formie:
m
aigik = k, k = 0, . . . , m (170)
i=0
n n
gij = xi+k, k = f(xj)xk (171)
j
j
j=0 j=0
Uwagi:
" Dla m n rozwi膮zanie jest jednoznaczne
" Dla m = n rozwi膮zaniem jest wielomian interpolacyjny
" Dla m = 1 mamy aproksymacj臋 funkcj膮 liniow膮 y = a0 + a1x:
EXY - EXEY cov(X, Y )
a1 = = , a0 = EY - a1EX (172)
EX2 - (EX)2 var(X)
" Zmiana stopnia wielomianu aproksymuj膮cego wymaga ponownego
rozwi膮zania uk艂adu normalnego.
Metody Numeryczne II 104
Problem z艂ego uwarunkowania
Za艂贸偶my, 偶e punkty xi s膮 r贸wno rozmieszczone wewn膮trz przedzia艂u [0, 1]. Dla
du偶ych warto艣ci n:
n
1
n
gik = xi+k H" n xi+k = , i, k = 0, . . . , m (173)
j j
i + k + 1
0
j=1
Macierz wsp贸艂czynnik贸w uk艂adu normalnego:
钆 艂艂
1 1
1
2 m+1
1 1 1
锱 艣艂
2 3 m+2
锱 艣艂
A = (174)
锱 艣艂
. . . .
. . . .
鹋
. . . .
1 1 1
m+1 m+2 2m+1
Macierz odwrotna ma z kolei bardzo du偶e warto艣ci co powoduje bardzo du偶e b艂臋dy
zaokr膮gle艅.
Metody Numeryczne II 105
Aproksymacja za pomoc膮 wielomian贸w ortogonalnych
Definicja 10 Funkcje f(x) i g(x) nazywamy ortogonalnymi na dyskretnym
zbiorze punkt贸w x0, . . . , xn gdy
n
f(xi)g(xi) = 0 (175)
i=0
przy
n n
f2(xi) > 0, g2(xi) > 0. (176)
i=0 i=0
Uwaga: Zastosowanie wielomian贸w ortogonalnych jako bazy przestrzeni
eliminuje problem z艂ego uwarunkowania macierzy uk艂adu normalnego.
Macierz ta staje si臋 diagonaln膮, z elementami przek膮tnej
n
ajj = 2(xi). (177)
j
i=0
Metody Numeryczne II 106
Metoda wielomian贸w ortogonalnych
1. Ka偶dy wielomian jest kombinacj膮 liniow膮 wielomian贸w ortogonalnych.
2. Wielomian aproksymacji to ten o najmniejszym b艂臋dzie kwadratowym.
Niech P0, . . . , Pm b臋d膮 wielomianami ortogonalnymi o stopniach r贸wnych
indeksom oraz Qm " m. W贸wczas istnieje jednoznaczny rozk艂ad
Qm(x) = b0P0(x) + + bmPm(x). (178)
Wsp贸艂czynnik bk 艂atwo wyznaczy膰 mno偶膮c to偶samo艣膰 (178) przez Pk(x) i sumuj膮c
r贸wno艣ci uzyskane dla podstawie艅 x0, . . . , xn. Otrzymujemy:
n n n
2
Qm(xi)Pk(xi) = b0 P0(xi)Pk(xi) + + bk Pk (xi) +
i=0 i=0 i=0
(179)
n
+ bm Pm(xi)Pk(xi)
i=0
Metody Numeryczne II 107
Wobec ortogonalno艣ci rodziny {Pi} mamy:
n n
2
Qm(xi)Pk(xi) = bk Pk (xi). (180)
i=0 i=0
St膮d
n
Qm(xi)Pk(xi)
i=0
bk = , k = 0, . . . , m. (181)
n
2
Pk (xi)
i=0
Szukamy wsp贸艂czynnik贸w b0, . . . , bm minimalizuj膮cych
n
S = [b0P0(xi) + + bmPm(xi) - f(xi)]2
钆傠艂 艂艂
i=0
雠2 肱 雠
(182)
n m m
锱傢艂 bjPj(xi)艂艂 - 2 砼 bjPj(xi)艂艂 f(xi) + f2(xi)艣艂
=
鹋
i=0 j=0 j=0
Metody Numeryczne II 108
W wyniku ortogonalno艣ci {Pi} oraz podstawie艅:
n n
2
sj = Pj (xi), cj = Pj(xi)f(xi), j = 0, . . . , m (183)
i=0 i=0
otrzymujemy
n m n m n
2
S = b2Pj (xi) - 2 bjPj(xi)f(xi) + f2(xi) =
j
i=0 j=0 i=0 j=0 i=0
(184)
m m n m n
= b2sj - 2 bjcj + f2(xi) = b2sj - 2bjcj + f2(xi)
j j
j=0 j=0 i=0 j=0 i=0
Poniewa偶
c2
cj cj c2
j j
b2sj - 2bjcj = sj b2 - 2bj = sj b2 - 2bj + - =
j j j
sj sj s2 sj
j
(185)
2
c2
cj
j
= sj bj - -
sj sj
Metody Numeryczne II 109
b艂膮d kwadratowy
2 m
m n
c2
cj
j
S = sj bj - - + f2(xi) (186)
sj sj i=0
j=0 j=0
osi膮ga minimum dla
cj
bj = . (187)
sj
Zatem wielomian aproksymacyjny ma posta膰
m
cj
Qm(x) = Pj(x). (188)
sj
i=0
Metody Numeryczne II 110
Przypadek r贸wnoodleg艂ych punkt贸w
x-x0
Dla zbioru punkt贸w xi = x0 + ih, i = 0, . . . , n wprowadzamy t = (liczby
h
ca艂kowite) i wyznaczamy odpowiednie unormowane wielomiany ortogonalne:
Pk,n(t) = 1 + b1t + b2t(t - 1) + . . . + bkt(t - 1) (t - k + 1)
(189)
Pk,n(0) = 0 k = 0, 1, . . . , m
Dowolny wielomian Pj,n(t) mo偶na zapisa膰 w postaci:
j
Pj,n(t) = ds(t + s)[s], (190)
s=0
ozn
gdzie ds s膮 sta艂ymi, a wielomiany czynnikowe r[j] = r(r - 1) (r - j + 1).
Warunki ortogonalno艣ci:
j
n n
Pj,n(i)Pk,n(i) = ds(i + s)[s]Pk,n(i) = 0, j = 0, . . . , k - 1, (191)
i=0 i=0 s=0
Metody Numeryczne II 111
wi臋c by je spe艂ni膰 wystarczy by
n
(i + j)[j]Pk,n(i) = 0, j = 0, . . . , k - 1. (192)
i=0
Sumuj膮c
(i + j)[j]Pk,n(i) = (i + j)[j] + b1(i + j)[j+1] + + bk(i + j)[j+k] (193)
po wszystkich i = 0, . . . , n, otrzymujemy
n n n n
(i+ j)[j]Pk,n(i) = (i+j)[j] + b1 (i+ j)[j+1] + + bk (i+j)[j+k] (194)
i=0 i=0 i=0 i=0
Wykorzystuj膮c, 偶e
n
(n + j + 1)[s+1]
(i + j)[s] = , s = j, . . . , j + k, (195)
s + 1
i=0
Metody Numeryczne II 112
i dziel膮c obie strony przez (n + j + 1)[j+1] dostajemy dla j = 0, . . . , k - 1
1 n n[k] Q(j)
+ b1 + + bk = = 0, (196)
1 + j 2 + j k + 1 + j (k + 1 + j)[k]
gdzie Q(j) jest wielomianem stopnia co najwy偶ej k o zerach j = 0, . . . , k - 1.
Przyjmujemy
ak = bkn[k], Q(j) = j(j - 1) (j - k + 1)Ck = Ckj[k], (197)
oraz mno偶ymy (196) przez (1 + j). Dostajemy
1 + j 1 + j Ckj[k]
1 + a1 + + ak = , (198)
2 + j k + 1 + j (2 + j)(3 + j) (k + 1 + j)
sk膮d dla j = -1 wynika
1 2 k
Ck = = (-1)k. (199)
(-1)(-2) (-k)
Metody Numeryczne II 113
Aby wyznaczy膰 as dla s = 1, . . . , k mno偶ymy (196) przez (s+1+j) i przyjmujemy
j = -(s + 1). Otrzymujemy:
(-1)kj[k](s + 1 + j) (-1)k(-s - 1)(-s - 2) (-s - k)
as = =
(k + 1 + j)[k+1] (k - s)(k - s - 1) 1 (-1) (-s)
(-1)2k(k + s)[k] (k + s)[k]s! (k + s)! k!
(200)
= = (-1)s = (-1)s
(k - s)!(-1)ss! s!s!(k - s)! s!k! s!(k - s)!
k + s k
= (-1)s
s s
Zatem ostateczna posta膰 wielomian贸w ortogonalnych (Czebyszewa/Grama) to:
k
k k + s t[s]
Pk,n(t) = (-1)s , k = 0, . . . , m, (201)
s s n[s]
s=0
a wielomian aproksymacyjny ma posta膰
m
cj x - x0
Qm(x) = Pj,n . (202)
sj h
i=0
Metody Numeryczne II 114
Aproksymacja funkcji ci膮g艂ych
Znamy funkcj臋, chcemy j膮 aproksymowa膰 inn膮 funkcj膮.
Dla funkcji aproksymowanej f ci膮g艂ej na [a, b], funkcji wagowej (x) oraz rodziny
funkcji aproksymuj膮cych
F (x) = a00(x) + + ann(x) (203)
minimalizujemy warto艣膰 funkcji
2
n
b b
En = (x) [F (x) - f(x)]2 dx = aii(x) - f(x) dx. (204)
a a
i=0
Rozwi膮zujemy uk艂ad r贸wna艅:
n
b
En
= aii(x) - f(x) j(x)dx
aj
a
i=0
(205)
n
b b
= ai i(x)j(x) - f(x)j(x)dx = 0, j = 0, . . . , n.
a a
i=0
Metody Numeryczne II 115
"
Przyk艂ad: Znalez膰 najlepsz膮 aproksymacj臋 kwadratow膮 funkcji f(x) = x w
przedziale [0, 1] wielomianem stopnia pierwszego Q(x) = a0 + a1x.
Otrzymujemy uk艂ad r贸wna艅:
1 1 1
"
a0 12dx + a1 xdx = xdx
0 0 0
(206)
1 1 1
"
a0 xdx + a1 x2dx = xxdx
0 0 0
Czyli:
1 2
a0 + a1 =
2 3
(207)
1 1 2
a0 + a1 =
2 3 5
Rozwi膮zanie:
4 4 4 4
a0 = , a1 = , Q(x) = + x (208)
15 5 15 5
Metody Numeryczne II 116
Obliczenia mog膮 by膰 mniej uci膮偶liwe przy u偶yciu wielomian贸w ortogonalnych.
Przyk艂ad: Znalez膰 najlepsz膮 aproksymacj臋 kwadratow膮 funkcji f(x) = ex w
przedziale [0, 1] wielomianem stopnia drugiego przy u偶yciu wielomian贸w
ortogonalnych Gaussa-Legendre a (154).
t+1
t+1
2
Podstawiamy x = i przybli偶amy funkcj臋 e wielomianem
2
Q(t) = a0P0(t) + a1P1(t) + a2P2(t):
1 1
t+1
2
a0 12dt = e dt
-1 -1
1 1
t+1
2
a1 t2dt = te dt (209)
-1 -1
2
1 1
t+1
1 1
2
a2 t2 - dt = t2 - e dt
3 3
-1 -1
W rozwi膮zaniu Q(t) podstawiamy t = 2x - 1.
Metody Numeryczne II 117
Metoda wyr贸wnywania
" W przypadku pewnych nieliniowych zale偶no艣ci:
trudno skonstruowa膰 odpowiedni wielomian uniwersalny
E
trudno rozwi膮za膰 wprost uk艂ad r贸wna艅 nieliniowych
ai
" Zastosowanie znajduje metoda wyr贸wnywania (linearyzacji).
" Je艣li zale偶no艣膰 mi臋dzy x a y jest nieliniowa, a istniej膮 odwzorowania
wzajemnie jednoznaczne , takie, 偶e zale偶no艣膰 mi臋dzy X = (x, y) a
Y = (x, y) daje si臋 opisa膰 wielomianem uniwersalnym, to aproksymuje si臋
zale偶no艣膰 mi臋dzy X i Y i minimum przenosi do wyj艣ciowych zmiennych.
" Nieznanej zale偶no艣ci mi臋dzy x a y cz臋sto szuka si臋 eksperymentalnie
(najcz臋艣ciej wizualnie) po艣r贸d nast臋puj膮cych:
b
y = ax + b y = axb y = abx y = a +
x
(210)
1 x
y = y = y = a ln x + b
ax + b ax + b
Metody Numeryczne II 118
W ka偶dym z wymienionych przypadk贸w zale偶no艣ci nieliniowych aproksymujemy
liniow膮 zale偶no艣膰
Y = 膮X + (211)
gdzie stosujemy nast臋puj膮ce podstawienia:
zale偶no艣膰 x i y podstawienia
y = bxa X = log x, Y = log y, = log b
y = bax Y = log y, 膮 = log a, = log b
b
y = a + Y = xy
x
1 1
y = Y =
ax+b y
x x
y = Y =
ax+b y
y = a ln x + b X = log x
Brak specyfikacji podstawienia dla danej zmiennej oznacza odpowiednio:
X = x, Y = y, 膮 = a, = b (212)
Uwaga: Minimalizacja b艂臋du 艣redniokwadratowego dla ww podstawie艅
gwarantuje optymaln膮 warto艣膰 funkcji b艂臋du dla oryginalnych zmiennych tylko
w贸wczas gdy b艂膮d ten jest zerowy.
Metody Numeryczne II 119
Aproksymacja funkcjami sklejanymi
Przypadek r贸wnoodleg艂ych w臋z艂贸w
Dany jest podzia艂 " = {xi : i = 0, . . . , n} przedzia艂u [a, b] taki, 偶e xi = x0 + ih,
b-a
gdzie h = .
n
Twierdzenie 22 (Baza przestrzeni funkcji sklejanych) Niech funkcje
艢3(x), i = -1, 0, . . . , n, n + 1 b臋d膮 okre艣lone wzorem
i
(x - xi-2)3 dla x " [xi-2, xi-1]
h3 + 3h2(x - xi-1) + 3h(x - xi-1)2 - 3(x - xi-1)3 dla x " [xi-1, xi]
1
艢3(x) =
h3 + 3h2(xi+1 - x) + 3h(xi+1 - x)2 - 3(xi+1 - x)3 dla x " [xi, xi+1]
i
h3
(xi+2 - x)3 dla x " [xi+1, xi+2]
0 w pozosta艂ych przyp.
(213)
Funkcje 艢i(x) = 艢3(x) okre艣lone na [a, b] stanowi膮 baz臋 przestrzeni funkcji sklejanych
i
trzeciego stopnia na ".
Metody Numeryczne II 120
Wniosek 23 Ka偶d膮 funkcj臋 sklejan膮 trzeciego stopnia S"(x) mo偶na
przedstawi膰 w postaci:
n+1
S"(x) = ci艢i(x), (214)
i=-1
gdzie ci s膮 liczbami rzeczywistymi.
W艂asno艣ci funkcji bazowych:
Wykres 艢3(x)
i
Warto艣ci 艢3(x) i pochodnych
i
x xi-2 xi-1 xi xi+1 xi+2
艢3(x) 0 1 4 1 0
i
2
3 3
3
0 0 - 0
h h
艢i (x) 2 2
6 12 6
艢3(x) 0 - 0
i
h2 h2 h2
xi-2 xi-1 xi xi+1 xi+2
S膮 klasy C2((- ", "))
Metody Numeryczne II 121
Aproksymacja funkcjami sklejanymi z r贸wnoodleg艂ymi w臋z艂ami
Minimalizujemy:
2
n+1
b
E = f(x) - ci艢3(x) dx (215)
i
a
i=-1
Minimum wyznacza rozwi膮zanie uk艂adu n + 3 r贸wna艅 z n + 3 niewiadomymi:
E
= 0, j = -1, . . . , n + 1 (216)
cj
czyli:
n+1
b b
ci 艢3(x)艢3(x) = f(x)艢3(x), j = -1, . . . , n + 1 (217)
i j j
a a
i=-1
b
Wprowadzaj膮c aij = 艢3(x)艢3(x) otrzymujemy uk艂ad r贸wna艅 o symetrycznej
i j
a
siedmiodiagonalnej macierzy g艂贸wnej.
Metody Numeryczne II 122
Aproksymacja funkcjami sklejanymi dla dyskretnego zbioru punkt贸w
Odci臋te w臋z艂贸w: {zi, i = 0, . . . , m}, gdzie m > n + 3.
2
m n+1
J = f(zk) - ci艢3(zk) (218)
i
k=0 i=-1
Z przyr贸wnania pochodnych cz膮stkowych do zera mamy uk艂ad:
n+1 m
bijci = f(zk)艢3(zk), j = -1, . . . , n + 1 (219)
j
i=-1 k=0
gdzie
m
bij = 艢3(zk)艢3(zk) (220)
i j
k=0
jest symetryczn膮 siedmiodiagonaln膮 macierz膮.
Je艣li wyznacznik uk艂adu jest r贸偶ny od zera, to jego rozwi膮zanie okre艣la minimum
funkcji J.
Metody Numeryczne II 123
Aproksymacja jednostajna
Funkcja b艂臋du:
E = supx"[a,b]|f(x) - F (x)| (221)
" Z twierdzenia Weierstrassa (tw. 21) wynika, 偶e dowoln膮 funkcj臋 f ci膮g艂膮 na
[a, b] mo偶na aproksymowa膰 jednostajnie wielomianem z dowoln膮
dok艂adno艣ci膮.
" Twierdzenie Borela m贸wi, 偶e dla ka偶dej funkcji f(x) na przedziale [a, b] oraz
n " N istnieje wielomian Wn(x) o minimalnym b艂臋dzie aproksymacji
jednostajnej.
" Twierdzenie Czebyszewa m贸wi, 偶e taki wielomian jest tylko jeden.
Brak og贸lnej metody wyznaczania wielomianu najlepszego przybli偶enia
jednostajnego.
Metody Numeryczne II 124
Metoda szereg贸w pot臋gowych
Je偶eli funkcj臋 f(x) mo偶na na [a, b] rozwin膮膰 w szereg Taylora, to przybli偶eniem
mo偶e by膰 odpowiednie obci臋cie szeregu:
f2 (x0) f2 2 (x0) f(n)(x0)
Wn(x) = f(x0)+ (x-x0)+ (x-x0)2 + + (x-x0)n (222)
1! 2! n!
Oszacowanie b艂臋du z reszty w postaci Lagrange a:
f(n+1)(c)
Mn+1
|f(x) - Wn(x)| = (x - x0)n+1 n+1 (223)
(n + 1)! (n + 1)!
dla c le偶膮cego mi臋dzy x i x0, b臋d膮cego promieniem stosownego otoczenia
punktu x0 oraz Mn+1 |f(n+1)(x)| dla x " [a, b].
Odpowiedni膮 dok艂adno艣膰 mo偶na uzyska膰 dobieraj膮c odpowiednio du偶e n.
"
Przyk艂ad: Aproksymowa膰 metod膮 szereg贸w pot臋gowych funkcj臋 x
wielomianem stopnia 2 w przedziale [0,5, 1,5]. Oszacowa膰 b艂膮d.
Metody Numeryczne II 125
Przybli偶enia Pad間o
Dok艂adniejsze od wielomianowych cz臋sto okazuj膮 si臋 by膰 wymierne:
Ln(x)
Rn,k(x) = , (224)
Mk(x)
gdzie
Ln(x) = a0 + a1x + + anxn
(225)
Mk(x) = b0 + b1x + + bkxk
" Parametry okre艣la si臋 tak, by w x0 = 0 warto艣ci funkcji f i Rn,k oraz jak
najwi臋cej ich pochodnych zr贸wna艂o si臋.
" Dla k = 0 uzyskujemy przybli偶enie szeregiem pot臋gowym Maclaurina.
" f(x) przybli偶a si臋 szeregiem Maclaurina i rozwi膮zuje uk艂ad r贸wna艅
wynikaj膮cych z za艂o偶enia r贸wno艣ci warto艣ci funkcji i pochodnych.
" Stosunkowo ma艂e b艂臋dy w pobli偶u zera, wi臋ksze dalej.
Metody Numeryczne II 126
Szeregi Czebyszewa
Przybli偶anie funkcji f(x) sum膮 cz臋艣ciow膮 szeregu
n
1
f(x) H" c0 + cjTj(x), (226)
2
j=1
gdzie
1
2 f(x)Tj(x)
"
cj = dx Tj(x) = cos(j arccos x). (227)
膭
1 - x2
-1
Alternatywnie, funkcj臋 mo偶na aproksymowa膰 funkcj膮 wymiern膮
n
aiTi(x)
i=0
Tn,k(x) = . (228)
k
biTi(x)
i=0
Warto艣ci a bi wyznacza si臋 tak, by w r贸偶nicy f(x) - Tn,k(x) =
"i,
1
c0 + cjTj(x) - Tn,k(x) znika艂y wsp贸艂czynniki przy Tj dla j = 0, . . . , N.
j=1
2
Metody Numeryczne II 127
Rozwi膮zywanie r贸wna艅 nieliniowych
Zadanie: Najcz臋stsze postaci:
1. Znalez膰 pierwiastki rzeczywiste r贸wnania f(x) = 0, funkcji f zmiennej
rzeczywistej.
2. Znalez膰 pierwiastki rzeczywiste i zespolone r贸wnania P (x) = 0 wielomianu P .
Metody iteracyjne
" Ci膮g kolejnych przybli偶e艅 {xi}.
" n-punktowy wz贸r iteracyjny
xi+1 = Fi(xi, xi-1, . . . , xi-n+1). (229)
" Funkcja Fi zwykle zale偶y od warto艣ci funkcji f i od jej pochodnych.
" Funkcja iteracyjna jest stacjonarna gdy jest niezale偶na od indeksu i. W贸wczas
pierwiastki s膮 punktami sta艂ymi F :
膮 = F (膮, . . . , 膮). (230)
Metody Numeryczne II 128
" Niech 膮 b臋dzie szukanym pierwiastkiem. B艂膮d i-tej iteracji i = 膮 - xi. Je艣li
limi " xi = 膮, to m贸wimy, 偶e metoda jest w punkcie 膮 rz臋du p 1 wtw gdy
|xi+1 - 膮| |i+1|
lim = lim = C = 0. (231)
i " i "
|xi - 膮|p |i|p
W贸wczas sta艂膮 C nazywamy sta艂膮 asymptotyczn膮 b艂臋du.
" Efektywno艣膰 metody - odwrotnie proporcjonalna do iloczynu kosztu 艃 jednej
iteracji i liczby iteracji I.
Niech rz臋dy dw贸ch metod b臋d膮 odpowiednio r贸wne p1 i p2. Dla du偶ych i
mamy w przybli偶eniu:
1 2
|(1) | = C1|(1)|p , |(2) | = C2|(2)|p . (232)
i+1 i i+1 i
Niech
Si = - log |(1)|, Ti = - log |(2)|. (233)
i i
Wtedy
Si+1 = - log C1 + p1Si, Ti+1 = - log C2 + p2Ti. (234)
Metody Numeryczne II 129
Rozwi膮zania:
(pi (pi
1 2
Si = S1pi - log C1 -1)/(p1-1), Ti = T1pi - log C2 -1)/(p2-1). (235)
1 2
Zak艂adamy, 偶e S1 = T1, koszt jednej iteracji to odpowiednio 艃1 i 艃2, a liczby
krok贸w potrzebne do uzyskania odpowiedniej dok艂adno艣ci to I1 oraz I2.
Zatem SI i TI s膮 w przybli偶eniu r贸wne. St膮d:
1 2
(pI2
2
C2 -1)/(p2-1)
1 2
S1(pI - pI ) + log = 0. (236)
1 2
(pI1
1
C1 -1)/(p1-1)
Najcz臋艣ciej pierwszy sk艂adnik dominuje, wi臋c mo偶na drugi zaniedba膰. Zatem:
I1 log p1 H" I2 log p2 H" const. (237)
Ostatecznie wskaznik efektywno艣ci mo偶na okre艣li膰 jako:
1 1
艃
E = log p albo E = p . (238)
艃
" Rz膮d metody jest cech膮 lokaln膮.
" Koszt jednej iteracji H" koszt liczenia warto艣ci funkcji i jej pochodnych.
Metody Numeryczne II 130
Metoda interpolacji odwrotnej
Zak艂adamy, 偶e w okolicy pierwiastka 膮 funkcji f istnieje funkcja odwrotna g = f-1
(膮 jest pierwiastkiem pojedynczym).
Niech xi-n, . . . , xi b臋d膮 kolejnymi przybli偶eniami 膮 (b膮dz pewnymi punktami z
otoczenia 膮 dla i = 0). Przyjmujemy
yj = f(xi-j), j = 0, . . . , n. (239)
Mo偶emy interpolowa膰 funkcj臋 g(y) wielomianem interpolacyjnym Lagrange a
n n
h(y) = Lj(y)g(yj) = Lj(y)xi-j. (240)
j=0 j=0
Metody Numeryczne II 131
W贸wczas kolejnym przybli偶eniem xi+1 pierwiastka funkcji f jest:
n n
xi-j(-1)ny0 yj-1yj+1 yn
h(0) = Lj(0)xi-j = .
(yj - y0) (yj - yj-1)(yj - yj+1) (yj - yn)
j=0 j=0
(241)
Ze wzoru Lagrange a:
g(n)()
膮 - xi+1 = (-1)n+1y0 yn, " I[y0, . . . , yn, 0]. (242)
n!
" Mo偶na u偶y膰 dowolnej metody interpolacji.
Metody Numeryczne II 132
Metoda bisekcji
Szukamy zera 膮 w przedziale (a, b), takim, 偶e f(a)f(b) < 0. Przyjmujemy
a0 = a, b0 = b, a tak偶e
ai + bi
xi = , i = 0, 1, . . . , (243)
2
oraz
(ai, xi) o ile f(ai)f(xi) < 0
(ai+1, bi+1) = . (244)
(xi, bi) o ile f(xi)f(bi) < 0
Zak艂adamy, 偶e f(xi) = 0, bo inaczej mamy rozwi膮zanie dok艂adne.
Poniewa偶
1
|xi - xi+1| = (b - a), i = 0, 1, . . . , (245)
2i+1
to
lim xi = 膮. (246)
i "
Metody Numeryczne II 133
Regula falsi
regula linia, falsus fa艂szywy fa艂szywo艣膰 za艂o偶enia liniowo艣ci funkcji
Za艂o偶enia: funkcja f jest ci膮g艂a, f(a)f(b) < 0.
Zak艂adamy (a0, b0) = (a, b). W i-tym kroku prowadzimy ci臋ciw臋 o r贸wnaniu
f(bi) - f(ai)
y - f(ai) = (x - ai) (247)
bi - ai
i przybli偶amy pierwiastek r贸wnania f(x) = 0 pierwiastkiem ci臋ciwy:
f(ai)
xi = ai - (bi - ai). (248)
f(bi) - f(ai)
W kolejnych krokach
(ai, xi) o ile f(ai)f(xi) < 0
(ai+1, bi+1) = . (249)
(xi, bi) o ile f(xi)f(bi) < 0
" na og贸艂 niestacjonarna, stacjonarna np. dla funkcji wypuk艂ych
Metody Numeryczne II 134
Oszacowania b艂臋d贸w metody regula falsi
Z twierdzenia Lagrange a o przyrostach:
f(xn) - f(膮) = f2 (c)(xn - 膮), dla pewnego c " I[xn, 膮] (250)
Poniewa偶 f(膮) = 0, wi臋c
|f(xn)|
|xn - 膮| , m = inf |f2 (x)| (251)
m x"[a,b]
Przy za艂o偶eniu, 偶e f2 i f2 2 nie zmieniaj膮 znaku w [a, b] jeden z ko艅c贸w przedzia艂u
nie zmienia si臋 w kolejnych krokach iteracji. Przyjmijmy f2 (x) > 0, f2 2 (x) > 0. St膮d
bi = b, i = 0, 1, . . . czyli
f(xi)
x-1 = a, xi+1 = xi - (b - xi), i = 0, 1, . . . (252)
f(b) - f(xi)
Mamy wi臋c
f(xi) - f(b)
f(膮) - f(xi) = (xi+1 - xi), (253)
xi - b
Metody Numeryczne II 135
i z tw. Lagrange a
(膮 - xi)f2 (i) = (xi+1 - xi)f2 (xi), i " (xi, 膮), xi " (xi, b). (254)
Po wymno偶eniu i dodaniu obustronnie -xi+1f2 (i) otrzymujemy
|f2 (xi) - f2 (i)|
|膮 - xi+1| = |xi+1 - xi|. (255)
|f2 (i)|
Poniewa偶 i, xi " [a, b] i f2 (x) > 0, wi臋c
|f2 (xi) - f2 (i)| M - m, m = inf |f2 (x)|, M = sup |f2 (x)|. (256)
x"[a,b]
x"[a,b]
Ostatecznie
M - m
|膮 - xi+1| |xi+1 - xi|. (257)
m
Oszacowanie to jest na og贸艂 pesymistyczne.
Dla przybli偶e艅 w niewielkim otoczeniu 膮 mo偶na szacowa膰:
f(xi+1) xi+1 - xi
|膮 - xi+1| H" H" (258)
f (xi+1) f(xi+1) - f(xi) |f(xi+1)|
2
Metody Numeryczne II 136
Metoda siecznych
" Ta sama zasada, co w regula falsi
" Rezygnujemy z za艂o偶enia f(a)f(b) < 0, czasami kosztem zbie偶no艣ci
" Metoda stacjonarna:
xi - xi-1
x0 = a, x1 = b, xi+1 = xi - f(xi) , i = 1, 2, . . . (259)
f(xi) - f(xi-1)
" Istotne znaczenie ma maksymalna graniczna dok艂adno艣膰 (wynikaj膮ca z
przyj臋tej arytmetyki) gdy xi - xi-1 jest tego samego rz臋du co oszacowanie
b艂臋du |膮 - xi|, to kolejne przybli偶enia mog膮 by膰 ca艂kowicie b艂臋dne.
Dodatkowe kryterium stopu: zaburzenie zmniejszania si臋 |f(xi)|.
" Wykrywanie rozbie偶no艣ci r贸偶nica pomi臋dzy kolejnymi przybli偶eniami ro艣nie.
Metody Numeryczne II 137
Metoda Newtona
" Za艂o偶enia: f(a)f(b) < 0, f2 oraz f2 2 nie zmieniaj膮 znaku na [a, b].
" Z tego ko艅ca przedzia艂u x0 " {a, b}, w kt贸rym f(x) ma ten sam znak co
f2 2 (x) prowadzimy styczn膮 do wykresu funkcji y = f(x). Punkt x1 przeci臋cia
stycznej z osi膮 odci臋tych to pierwsze przybli偶enie pierwiastka.
" Prowadzimy styczn膮 w punkcie x1 i wyznaczamy punkt x2 przeci臋cia tej
stycznej z osi膮.
" Uzyskany ci膮g x1, x2, . . . jest zbie偶ny monotonicznie do 膮.
Dow贸d:
1. Ci膮g jest monotoniczny.
2. Zbiega do 膮.
Przypadek f2 (x) > 0, f2 2 (x) > 0:
W贸wczas x0 = b oraz f(x0) > 0. Poka偶emy, 偶e xn+1 < xn dla n = 0, 1, . . ..
R贸wnanie stycznej:
y - f(xn) = f2 (xn)(x - xn). (260)
Metody Numeryczne II 138
Zero stycznej w punkcie
f(xn)
xn+1 = xn - . (261)
f2 (xn)
Ze wzoru Taylora:
1
0 = f(膮) = f(xn) + f2 (xn)(膮 - xn) + f2 2 (c)(膮 - xn)2, c " (膮, xn). (262)
2
St膮d
f(xn) 1 f2 2 (c)
膮 = xn - - (膮 - xn)2. (263)
f2 (xn) 2 f2 (xn)
Z (261):
1 f2 2 (c)
膮 - xn+1 = - (膮 - xn)2 < 0 (264)
2 f2 (xn)
Zatem f(xn+1) > 0 oraz xn+1 < xn, bo
f(xn)
xn+1 - xn = - < 0. (265)
f2 (xn)
Metody Numeryczne II 139
Ci膮g xn : n = 0, 1, . . . jest wi臋c monotoniczny i ograniczony, a zatem zbie偶ny
(do pewnego g). Z (261) przy n " mamy
f(g)
g = g - , (266)
f2 (g)
czyli f(g) = 0, a wi臋c g = 膮.
B艂膮d przybli偶enia xn.
Podobnie jak w metodzie regula falsi (str. 133)
|f(xn)|
|xn - 膮| , m = inf |f2 (x)| (267)
m x"[a,b]
Ze wzoru Taylora, dla pewnego n " I[xn, xn+1]:
1
f(xn+1) = f(xn) + f2 (xn)(xn+1 - xn) + f2 2 (n)(xn+1 - xn)2 (268)
2
Metody Numeryczne II 140
Z (261) f(xn) + f2 (xn)(xn+1 - xn) = 0, wi臋c
1
|f(xn+1)| M(xn+1 - xn)2, M = sup |f2 2 (x)|. (269)
2
x"[a,b]
Zatem
2
1 M M f(xn)
|膮 - xn+1| (xn+1 - xn)2 = . (270)
2 m 2m f2 (xn)
Alternatywne oszacowanie (w bliskim s膮siedztwie 膮):
f(xn)
|膮 - xn| H" . (271)
f (xn)
2
Uwagi:
" Metoda siecznych, to metoda Newtona z przybli偶aniem pochodnej ilorazami
r贸偶nicowymi.
" Znane zastosowanie metody Newtona - szukanie pierwiastk贸w kwadratowych
liczb rzeczywistych jako dodatnich pierwiastk贸w r贸wna艅 x2 - c = 0.
" Metoda Newtona zwykle potrzebuje mniej iteracji ni偶 metoda siecznych, ale
ma wy偶szy koszt jednej iteracji.
Metody Numeryczne II 141
Rz膮d stacjonarnej metody jednopunktowej
Dla metod, w kt贸rych
xn+1 = 艢(xn), (272)
gdzie 艢 jest wystarczaj膮co wiele razy r贸偶niczkowalna w otoczeniu O(膮)
wyznaczanej warto艣ci 膮, je偶eli xi " O(膮) oraz 艢(k)(膮) = 0 dla k = 1, 2, . . . , p - 1 i
艢p(膮) = 0, to zachodz膮:
(xi - 膮)p
xi+1 = 艢(xi) = 艢(膮) + 艢(p)(膮) + R(||xi - 膮||p+1)
p!
(273)
xi+1 - 膮 艢(p)(膮)
lim =
i "
(xi - 膮)p p!
Dla p 2 metoda jest rz臋du p. Metoda jest rz臋du pierwszego gdy p = 1 oraz
|膯2 (膮)| < 1.
Metody Numeryczne II 142
Dla metody Newtona mamy:
f(x)
艢(x) = x - . (274)
f2 (x)
Metoda jest rz臋du 2, bo:
f(x)f2 2 (x) f2 2 (膮)
艢(膮) = 膮, 艢2 (膮) = = 0, 艢2 2 (膮) = (275)
f2 (x)2 x=膮 f2 (膮)
Metody Numeryczne II 143
Modyfikacje metod
" Dotychczasowe rozwa偶ania dotyczy艂y pierwiastk贸w pojedynczych.
" Pierwiastki wielokrotne
Definicja 11 Liczba 膮 jest pierwiastkiem r-krotnym (r 2) r贸wnania
f(x) = 0 wtw gdy jest r - 1-krotnym pierwiastkiem r贸wnania f2 (x) = 0.
" Metoda po艂owienia, regula falsi, metoda siecznych
trac膮 sens dla pierwiastk贸w parzystokrotnych
dla krotno艣ci nieparzystych przy f(a)f(b) < 0 pozostaj膮 s艂uszne, cho膰 dla
siecznych obni偶a si臋 rz膮d.
" Metoda Newtona
pozostaje s艂uszna dla parzystokrotnych o ile w jednostronnym s膮siedztwie
zera spe艂nione s膮 za艂o偶enia sta艂o艣ci znak贸w pochodnych
dla pierwiastk贸w wielokrotnych obni偶ony rz膮d metody tylko liniowa
r-1
zbie偶no艣膰 ze sta艂膮 asymptotyczn膮 C =
r
Metody Numeryczne II 144
Przyk艂ady modyfikacji metod
" Znamy krotno艣膰 pierwiastka:
f(xn)
xn+1 = xn - r (276)
f2 (xn)
" Krotno艣膰 nieznana. Zamiast funkcji f(x) u偶ywamy funkcji
f(x)
u(x) = , (277)
f2 (x)
kt贸ra ma te same pierwiastki co f(x), ale wszystkie pojedyncze. W贸wczas dla
metod regula falsi i siecznych wzory (248) i (259) przyjmuj膮 posta膰
xn - xn-1
xn+1 = xn - u(xn) , (278)
u(xn) - u(xn-1)
a dla metody Newtona mamy:
u(x) f2 2 (xn)
xn+1 = xn - , przy u2 (xn) = 1 - u(xn). (279)
u2 (x) f2 (xn)
Metody Numeryczne II 145
Zera wielomian贸w
Liczba zer wielomianu
Definicja 12 Ci膮giem Sturma dla wielomianu p(x) nazywamy ci膮g
wielomian贸w p0(x), . . . , pm(x) taki, 偶e:
1. p0(x) = p(x).
2. p1(x) = p2 (x).
0
3. dla i = 2, . . . , m + 1 pi(x) jest reszt膮 z dzielenia pi-2 przez pi-1 wzi臋t膮 ze
znakiem przeciwnym.
4. pm+1(x) a" 0
ozn
Oznaczenie: N(z) = liczba zmian znaku w ci膮gu {pi(z) : i = 0, . . . , m, pi(z) = 0}
Twierdzenie 24 (Sturma) Je偶eli ci膮g {pi(x) : i = 0, . . . , m} jest ci膮giem
Sturma dla wielomianu p(x) na przedziale (a, b) oraz p0(a)p0(b) = 0, to liczba
r贸偶nych zer rzeczywistych wielomianu p(x) w przedziale (a, b) jest r贸wna
N(a) - N(b).
Metody Numeryczne II 146
ozn
Oznaczenie: M(z) = liczba zmian znaku w ci膮gu {p(i)(z) : i = 0, 1, . . .}.
Twierdzenie 25 (Fouriera) Je偶eli p(x) " n i p(a)p(b) = 0, to liczba zer
wielomianu p(x) w przedziale [a, b] jest r贸wna M(a) - M(b) lub jest od tej
liczby mniejsza o liczb臋 parzyst膮.
Oznaczenie: Dla wielomianu p(x) = a0xn + . . . + an-1x + an, gdzie a0 = 0,
k
ozn
L(z) = liczba zmian znaku ci膮gu {pk(z) = aizk-i : k = 0, . . . , n}.
i=0
Twierdzenie 26 (Laguerre a) Je偶eli p(x) " n i p(a)p(b) = 0, to liczba zer
wielomianu p(x) w przedziale [a, b] jest r贸wna L(a) - L(b) lub jest od tej liczby
mniejsza o liczb臋 parzyst膮.
Przypadek szczeg贸lny regu艂a Kartezjusza liczba pierwiastk贸w dodatnich
wielomianu jest r贸wna liczbie zmian znaku w ci膮gu wsp贸艂czynnik贸w minus liczba
parzysta.
Metody Numeryczne II 147
Lokalizacja przedzia艂u z zerami
" Zadanie: wyznaczy膰 przedzia艂 zawieraj膮cy wszystkie zera wielomianu.
" Redukcja problemu: wyznaczy膰 kres g贸rny R zbioru zer, bo je艣li
1 -1
p1(x) = xnp , p2(x) = p(-x), p1(x) = xnp (280)
x x
maj膮 kresy g贸rne odpowiednio R1, R2 i R3, to wszystkie dodatnie zera
1 1
wielomianu p(x) le偶膮 w ( , R), a ujemne w (-R2, - )
R1 R3
Twierdzenie 27 (Lagrange a) Niech a0 = 0 i ak (k 1) b臋dzie pierwszym
ujemnym wsp贸艂czynnikiem wielomianu p(x) = a0xn + . . . + an-1x + an.
Wszystkie zera tego wielomianu s膮 mniejsze ni偶
A
k
R = 1 + , (281)
|a0|
gdzie A oznacza maksimum modu艂u ujemnych wsp贸艂czynnik贸w wielomianu.
Je偶eli wszystkie wsp贸艂czynniki wielomianu s膮 nieujemne, to nie ma on zer
dodatnich.
Metody Numeryczne II 148
Wyznaczanie zer wielomian贸w
" Znamy kres g贸rny mo偶na metod膮 Newtona znalez膰 najwi臋kszy pierwiastek.
" Metoda Newtona mo偶e wolno zbiega膰 do najwi臋kszego zera. Metoda
podw贸jnego kroku:
p(xk)
xk+1 = xk - 2 (282)
p2 (xk)
Twierdzenie 28 Niech p(x) " n, n 2 ma wszystkie zera rzeczywiste
膮1 膮2 膮n. Niech 1 b臋dzie najwi臋kszym zerem wielomianu p2 (x)
(膮1 1 膮2). W przypadku n = 2 dodatkowo zak艂adamy 膮1 > 膮2.
Dla ka偶dego z > 膮1 dobrze okre艣lone s膮 liczby
p(z) p(z) p(y)
z2 = z - , y = z - 2 , y2 = y - (283)
p2 (z) p2 (z) p2 (y)
i spe艂niaj膮 nier贸wno艣ci
1 < y, 膮1 y2 z2 . (284)
Metody Numeryczne II 149
Dla ci膮gu przybli偶e艅 {xn : i = 1, 2, . . .} (zgodnych z (282)) mamy dwie
mo偶liwo艣ci:
1. "nxn 膮1 oraz limn " xn = 膮1.
2. "k p(xk )p(x0) < 0 oraz "k
0. W贸wczas ci膮g
0 0 0
p(yk)
y0 = xk , yk+1 = yk - , k = 1, 2, . . . (285)
0
p2 (yk)
zbiega do 膮1.
Deflacja oddzielamy pierwiastek 膮1 od pozosta艂ych wyznaczamy wielomian
stopnia n - 1
p(x)
p1(x) = . (286)
x - 膮1
Najwi臋kszym zerem wielomianu p1(x) jest 膮2, a dobrym punktem startowym jest
ewentualnie przestrzelone xk.
Uwaga na b艂臋dy zaokr膮gle艅 dla zera o najmniejszej warto艣ci bezwzgl臋dnej
wsp贸艂czynniki nale偶y wyznacza膰 od najwy偶szej pot臋gi, dla najwi臋kszej warto艣ci
bezwzgl臋dnej odwrotnie.
Metody Numeryczne II 150
yle uwarunkowane r贸wnania algebraiczne
" pierwiastki wielokrotne s膮 na og贸艂 zle uwarunkowane
" bliskie sobie pierwiastki
" nawet wyraznie rozdzielone pierwiastki:
Przyk艂ad: Wielomian o pierwiastkach rzeczywistych 1,2,. . . ,20:
p(z) = (z - 1)(z - 2) (z - 20) = z20 - 210z19 + + 20! (287)
Je艣li zmniejszy膰 a1 o 2-23 H" 10-10, to r贸wnanie ma pierwiastki:
16.730737466 膮 2.812624894i (288)
Dla schematu Hornera b艂膮d obliczonej warto艣ci p(x) mo偶e by膰 r贸wny
n
E = |(2i + 1)an-ixi|1.06%E艂 (289)
i=0
E
B艂膮d wzgl臋dny pierwiastka mo偶e osi膮ga膰
膮p2 (膮)
Metody Numeryczne II 151
St膮d wskaznik uwarunkowania:
|(2i + 1)an-i膮i| |(2i + 1)an-i膮i|
C = = (290)
|膮p2 (膮)| |ian-i膮i|
Dla badanego wielomianu:
C H" 1.35 1015, dla 膮 = 14. (291)
Gdy pierwiastkami s膮 liczby 膮i = 2i-21, i = 1, 2, . . . , 20, to wskaznik C = 3.83 103.
Miara rozdzielenia wzgl臋dnego: Dla k pierwiastk贸w rzeczywistych
膮i < . . . < 膮k:
k-1
膮i+1 - 膮i
def
s = 膮i = 0 (292)
|膮i|
i=1
W przypadku pierwiastk贸w 1, 2, . . . , 20 mamy s = 8.2 10-18. Dla pierwiastk贸w
膮i = 2i-21, i = 1, 2, . . . , 20 mamy s = 1.
Metody Numeryczne II 152
Algorytm Maehly ego dla wyeliminowania deflacji:
Poniewa偶
p2 (x) p(x)
p2 (x) = - , (293)
1
x - 膮1 (x - 膮1)2
wi臋c
p1(xk) p(xk)
xk+1 = xk - = xk - . (294)
p(xk)
p2 (xk)
1
p2 (xk) -
xk - 膮1
Og贸lnie mamy
p(x)
pj(x) = (295)
(x - 膮1) (x - 膮j)
j
p2 (x) p(x) 1
p2 (x) = - . (296)
j
(x - 膮1) (x - 膮j) (x - 膮1) (x - 膮j) x - 膮i
i=1
St膮d iterujemy
p(xk)
xk+1 = xk - (297)
j p(xk) .
p2 (xk) -
i=1
xk-膮i
Metody Numeryczne II 153
Metoda "2
Dla metod iteracyjnych rz臋du 1 mamy dla du偶ych warto艣ci k
膮 - xk+2 膮 - xk+1
H" H" C. (298)
膮 - xk+1 膮 - xk
St膮d wyznaczamy
xkxk+2 - x2
("xk+1)2
k+1
膮 H" = xk+2 - , (299)
xk+2 - 2xk+1 + xk "2xk
gdzie " jest operatorem r贸偶nicy progresywnej:
"xk = "1xk = xk+1 - xk, "i+1xk = "ixk+1 - "ixk (300)
Uwagi:
" proces "2 Aitkena
" znacznie lepsze przybli偶enie ni偶 przez xk
" tylko dla liniowej zbie偶no艣ci
Metody Numeryczne II 154
Minima funkcji jednej zmiennej
Sposoby znajdowania:
" szukanie rozwi膮za艅 r贸wnania f2 (x) = 0
" aproksymacja wielomianem
" metody podzia艂u
Metody podzia艂u
Lemat 29 Niech f : R R b臋dzie unimodaln膮 na [a, b] tj. ma minimum w
punkcie 膮 " [a, b], oraz f maleje na [a, 膮] oraz ro艣nie na [膮, b]. Aby
zlokalizowa膰 punkt 膮 w przedziale [a2 , b2 ] o mniejszej d艂ugo艣ci ni偶 przedzia艂
[a, b], wystarczy obliczy膰 warto艣膰 funkcji w dw贸ch punktach wewn膮trz [a, b].
Dow贸d: Niech a < t1 < t2 < b. Je艣li f(t1) < f(t2), to 膮 " [a, t2]. W przec. przyp.
膮 " [t1, b].
Metody podzia艂u polegaj膮 na wyznaczaniu ci膮g贸w przedzia艂贸w zst臋puj膮cych
[ai, bi], i = 0, 1, . . . , a0 = a, b0 = b oraz bi - ai n " 0.
Metody Numeryczne II 155
Metoda podzia艂u na trzy cz臋艣ci
2 1 1 2
ti = ai + bi, ti = ai + bi (301)
1 2
3 3 3 3
i
2
bi - ai = (b - a) (302)
3
W ka偶dej iteracji liczymy 2 warto艣ci funkcji.
Metoda po艂owienia
3 1 1 1 1 3
ti = ai + bi, ti = ai + bi, ti = ai + bi (303)
1 2 3
4 4 2 2 4 4
i
1
bi - ai = (b - a) (304)
2
W ka偶dej iteracji liczymy 3 lub 2 warto艣ci funkcji.
Metody Numeryczne II 156
Metoda optymalnych podzia艂贸w Johnsona
Cel: Minimalizacja liczby wylicze艅 warto艣ci funkcji.
Wykorzystujemy liczby Fibonacciego:
F0 = F1 = 1, Fi = Fi-1 + Fi-2, i = 2, 3, . . . (305)
Algorytm: Dla zadanej dok艂adno艣ci
b-a
" c = , a0 = a, b0 = b
" Znajdujemy N takie, 偶e FN-1 < c FN
" Dla i = 1, . . . , N - 2:
FN-i-1 FN-i
ti = ai-1 + (bi-1 - ai-1), ti = ai-1 + (bi-1 - ai-1) (306)
1 2
FN-i+1 FN-i+1
Je偶eli f(ti ) f(ti ), to ai-1 ! a, bi-1 ! ti . W przec. przyp. a ! ti , bi-1 ! b
1 2 2 1
Koszty: Liczymy warto艣ci funkcji w N - 1 punktach otrzymuj膮c ostatecznie
FN-1 FN-2 F2 b - a
bN-2 - aN-2 = (b - a) = 2 2 (307)
FN FN-1 F3 FN
Metody Numeryczne II 157
Metoda z艂otego podzia艂u
Cel:
" Przedzia艂 zmniejsza d艂ugo艣膰 w spos贸b sta艂y.
" Do kolejnych podzia艂贸w wykorzystuje si臋 jeden z punkt贸w z poprzedniego
kroku.
ti oraz ti wybieramy tak aby:
1 2
ti - ai = bi - ti = (bi - ai), " (0, 1)
2 1
(308)
bi - ti = (bi - ti )
2 1
"
5-1
Zatem: 2 + - 1 = 0 czyli = H" 0, 62 stosunek bok贸w z艂otego
2
prostok膮ta.
Metody Numeryczne II 158
Algorytm:
" a0 ! a, b0 ! b, t1 = a + (1 - )(b - a), t1 = b - (1 - )(b - a)
1 2
" Dla kolejnych i:
Je偶eli f(ti ) < f(ti ), to
1 2
ai+1 ! ai, bi+1 ! ti
2
(309)
ti+1 ! ti , ti+1 ! ai + (1 - )(bi - ai)
2 1 1
W przeciwnym przypadku:
ai+1 ! ti , bi+1 ! b
1
(310)
ti+1 ! ti , ti+1 ! bi - (1 - )(bi - ai)
1 2 2
Metody Numeryczne II 159
Por贸wnanie metod Johnsona i z艂otego podzia艂u
Przy za艂o偶eniu K oblicze艅 funkcji:
Metoda Johnsona:
b - a
N = K + 1, |t - 膮| (311)
FK+1
Metoda z艂otego podzia艂u:
b - a 1
|t - 膮| , gdzie Gi = . (312)
2GK-1 i
Por贸wnanie:
2GK-1 < FK+1 < 2GK (313)
Metody Numeryczne II 160
Uk艂ady r贸wna艅 liniowych
Zadanie: rozwi膮za膰 uk艂ad
钆 艂艂 钆 艂艂
a11 . . . a1n b1
锱 艣艂 锱 艣艂
. . .
. . .
Ax = b, A = , b = (314)
鹋 鹋
. . .
an1 . . . ann bn
Zadanie r贸wnowa偶ne znalezieniu macierzy odwrotnej A-1, bo
x = A-1b (315)
Z drugiej strony: i-ta kolumna ai macierzy A-1 jest rozwi膮zaniem uk艂adu
Ax = ei (316)
Metody Numeryczne II 161
Eliminacja Gaussa
Uk艂ad Ax = b przekszta艂camy do
钆 艂艂 钆 艂艂
r11 . . . r1n c1
锱 艣艂 锱 艣艂
. .
.
. . .
Rx = c, R = , c = , (317)
鹋 鹋
.
. .
0 rnn cn
tak by rozwi膮zanie wyznaczy膰 jako:
n
ci - rikxk
k=i+1
xi = , i = n, n - 1, . . . , 1. (318)
rii
Metody Numeryczne II 162
Pierwszy krok: Sprowadzi膰 macierz (A, b) do macierzy (A2 , b2 ):
钆 艂艂
钆 艂艂
a2 a2 . . . a2 b2
11 12 1n 1
a11 . . . a1n b1
锱 艣艂
0 a2 . . . a2 b2
22 2n 2
锱 艣艂 锱 艣艂
. . .
. . .
(A, b) = , (A2 , b2 ) = (319)
锱 艣艂
鹋 . . . .
. . .
. . . .
鹋
. . . .
an1 . . . ann bn
0 a2 . . . a2 b2
n2 nn n
tak, aby uk艂ady Ax = b i A2 x = b2 mia艂y te same rozwi膮zania.
Je艣li a11 = 0, to wystarczy wyliczy膰
ai1 ai1
a2 = aij - a1j , b2 = bi - b1 , i = 2, . . . , n, j = 1, . . . , n. (320)
ij i
a11 a11
W praktyce tzw. wyb贸r elementu g艂贸wnego o maksymalnym module:
" w kolumnie konieczna zamiana odpowiednich wierszy
" w ca艂ej macierzy konieczna zamiana odpowiednich wierszy i kolumn
Rozwi膮zanie pozostaje niezmienione z dok艂adno艣ci膮 do permutacji zmiennych.
Kolejne kroki: To samo co w pierwszym dla macierzy (A2 , b2 ) z poprzedniego
kroku, bez pierwszego wiersza i pierwszej kolumny.
Metody Numeryczne II 163
Eliminacja Gaussa w postaci macierzowej:
(A, b) = (A0, b0) (A1, b1) (An-1, bn-1) = (R, c), (321)
gdzie
(Aj, bj) = GjPj(Aj-1, bj-1)
(322)
(R, c) = Gn-1Pn-1 G1P1(A, b),
1
.
1 0
.
.
.
.
. 0 1
.
1
.
.
Gj = , Pj = (323)
-lj+1,j 1
1 0
.
.
. .
. .
.
.
.
0 -ln,j 0 1
1
Metody Numeryczne II 164
Rozk艂ad na iloczyn macierzy tr贸jk膮tnych
T0 = (A, b).
Dla j = 1, . . . , n:
钆
r11 r12 . . . r1j r1, j + 1 . . . r1n c1 艂艂
. . .
锱 艣艂
. . .
锱 艣艂
21 r22 . . . r2j . . .
锱 艣艂
. . . .
锱 艣艂
. . . .
锱 艣艂
31 32 . . . .
锱 艣艂
. .
锱 艣艂
. .
Tj = (324)
锱
. . rjj rj,j+1 . . . rjn cj 艣艂 .
锱 艣艂
锱 . . 艣艂
. .
锱 艣艂
. . j+1,j aj . . . aj bj
j+1,j+1 j+1,n j+1
锱 艣艂
锱 艣艂
. . . . . .
. . . . . .
鹋
. . . . . .
n,1 n2 . . . nj aj . . . aj bj
nn n
n,j+1
Kolumny z ij s膮 permutacj膮 odpowiednich -lij.
Metody Numeryczne II 165
Otrzymujemy:
LR = Pn-1Pn-2 . . . P1A = PA, (325)
gdzie
钆 艂艂
1 0
钆 艂艂
tn . . . tn
锱 艣艂
11 1n
.
锱倀n . 艣艂
.
锱 艣艂
21 .
锱 艣艂
.
L = , R = . (326)
鹋
.
锱 艣艂
.
.
. .
鹋
.
.
0 tn
nn
tn . . . tn 1
n1 n,n-1
Uwagi:
" Rozwi膮zanie uk艂adu Ax = b: poniewa偶 PAx = LRx = Pb, wi臋c x mo偶na
wyznaczy膰 w dw贸ch krokach:
Lu = Pb, Rx = u. (327)
" Nie ka偶da nieosobliwa macierz ma rozk艂ad tr贸jk膮tny A = LR.
" Z艂o偶ono艣膰 O(n3).
Metody Numeryczne II 166
Schematy zwarte eliminacji Gaussa
" G艂贸wnie dla celu oblicze艅 r臋cznych mniej wynik贸w po艣rednich.
" Eliminacja Gaussa w k-tym kroku wyznacza k-t膮 kolumn臋 L i k-ty wiersz R.
" A = LR jest r贸wnowa偶ne uk艂adowi n2 r贸wna艅 z n(n + 1) niewiadomymi:
min(i,j)
aij = liprpj. (328)
p=1
W k-tym kroku:
k k
akj = lkprpj, j k, aik = liprpk, i > k. (329)
p=1 p=1
Metody Numeryczne II 167
Metoda Doolittle a: lkk = 1, k = 1, . . . , n
k-1
rkj = akj - lkprpj, j = k, . . . , n,
p=1
(330)
k-1
aik - liprpk
p=1
lik = , i = k + 1, . . . , n.
rkk
Metoda Crouta: rkk = 1, k = 1, . . . , n
k-1
lik = aik - liprpk, i = k, . . . , n,
p=1
(331)
k-1
akj - lkprpj
p=1
rkj = , j = k + 1, . . . , n.
lkk
Metody Numeryczne II 168
Algorytm Gaussa Jordana
Cel: Wyznaczenie macierzy odwrotnej A-1. Macierz A realizuje odwzorowanie:
a11x1 + + a1nxn = y1
.
.
(332)
.
an1x1 + + annxn = yn
Wyznaczamy odwzorowanie odwrotne poprzez zamiany zmiennych. Zmienn膮 x1
zamieniamy na takie yr, dla kt贸rego
|ar1| = max |ai1|. (333)
i
Je艣li ar1 = 0, to macierz A jest osobliwa.
Zamieniamy miejscami r贸wnania 1 i r otrzymuj膮c uk艂ad Ax = y, gdzie a11 = ar1
oraz y1 = yr.
Metody Numeryczne II 169
Po zamianie zmiennych x1 i y1 = yr otrzymujemy uk艂ad:
a2 y1 + a2 x2 + + a2 xn = x1
11 12 1n
a2 y1 + a2 x2 + + a2 xn = y2
21 22 2n
, (334)
.
.
.
a2 y1 + a2 x2 + + a2 xn = yn
n1 n2 nn
gdzie dla i, j = 2, . . . , n
1 a1j ai1 ai1a1j
a2 = , a2 = - , a2 = , a2 = aij - . (335)
11 1j i1 ij
a11 a11 a11 a11
W ka偶dym kolejnym kroku k = 2, . . . , n zamieniamy w analogiczny spos贸b
zmienn膮 xk przez jedn膮 ze zmiennych yk, . . . , yn.
Metody Numeryczne II 170
Og贸lnie: konstruujemy ci膮g macierzy:
A = A0 An = A-1, (336)
gdzie dla k = 1, . . . , n macierz Ak = (a2 ) powstaje z Ak-1 = (aij) w nast臋puj膮cy
ij
spos贸b:
1. Wyznaczamy r takie, 偶e
|ark| = max |aik|. (337)
i k
Je艣li ark = 0, to macierz A jest osobliwa Koniec!.
2. Zamieniamy wiersze k i r macierzy Ak-1 otrzymuj膮c macierz A = (aik).
3. Dla i, j = k liczymy:
1 akj aik aikakj
a2 = , a2 = - , a2 = , a2 = aij - . (338)
kk kj ik ij
akk akk akk akk
Metody Numeryczne II 171
Rozk艂ad Choleskiego
Definicja 13 Macierz zespolon膮 A wymiaru n n, nazywamy dodatnio
okre艣lon膮, gdy:
1. A = AH, tzn. A jest macierz膮 hermitowsk膮.
2. Dla wszystkich x " Cn, x = 0 zachodzi xHAx > 0.
Hermitowsk膮 macierz nazywamy dodatnio p贸艂okre艣lon膮 gdy dla
wszystkich x " Cn zachodzi xHAx 0.
Twierdzenie 30 Dla ka偶dej dodatnio okre艣lonej macierzy A wymiaru n n
istnieje jednoznaczna tr贸jk膮tna dolna macierz L wymiaru n n, spe艂niaj膮ca
warunki lkk > 0, k = 1, . . . , n, oraz
A = LLH. (339)
Je偶eli A jest rzeczywista, to r贸wnie偶 L jest rzeczywista.
Metody Numeryczne II 172
R贸wnanie (339) oznacza, 偶e dla k = 1, . . . , n, i = k, . . . , n zachodzi:
aik = li1lk1 + li2lk2 + + linlkn. (340)
Algorytm:
Dla k = 1, . . . , n:
k-1 2
lkk = akk - lkp, (341)
p=1
Dla i = k + 1, . . . , n:
k-1
aik - liplkp
p=1
lik = . (342)
lkk
Uwagi:
" Wykorzystujemy tylko g贸rn膮 tr贸jk膮tn膮 cz臋艣膰 macierzy A.
" Z艂o偶ono艣膰 O(n3).
" n razy liczymy pierwiastek
"
" |lkj| akk, k = 1, . . . , n j = 1, . . . , k.
Metody Numeryczne II 173
B艂臋dy zaokr膮gle艅 w eliminacji Gaussa
Cz臋艣ciowy wyb贸r elementu g艂贸wnego
" zabezpiecza przed zerowym elementem g艂贸wnym w macierzy nieosobliwej
" ma wp艂yw na b艂臋dy numeryczne
Przyk艂ad: Uk艂ad
0.005 1 x 0.5
= (343)
1 1 y 1
4950
ma rozwi膮zanie (5000, ) H" (0.503, 0.497).
9950 9950
Stosujemy 2-cyfrow膮 precyzj臋 oblicze艅.
Je艣li elementem g艂贸wnym jest a11, to mamy:
0.005 1 x 0.5
= , (x, y) = (0, 0.5). (344)
0 -200 y -99
Je艣li wybierzemy a21 jako element g艂贸wny, to otrzymamy:
1 1 x 1
= , (x, y) = (0.5, 0.5). (345)
0 1 y 0.5
Metody Numeryczne II 174
Cz臋艣ciowy wyb贸r elementu g艂贸wnego to nie wszystko:
Z r贸wnowa偶nym do (343) uk艂adem
1 200 x 100
= (346)
1 1 y 1
s膮 te same problemy, cho膰 a11 = a21.
Skalowanie macierzy: Przemno偶enie pierwszego wiersza przez 200, to
zast膮pienie macierzy A przez = DA (i r贸wnocze艣nie b przez b = Db), gdzie
200 0
D = . (347)
0 1
Skalowa膰 mo偶na zar贸wno wiersze jak i kolumny, w贸wczas otrzymujemy uk艂ad:
D1AD2(D-1x) = D1AD2y = D1b, (348)
2
gdzie macierze D1 i D2 s膮 diagonalne.
Metody Numeryczne II 175
Nie ma sposobu skalowania zapewniaj膮cego numeryczn膮 stabilno艣膰 dla
dowolnej macierzy.
Heurystyka: macierze zr贸wnowa偶one, tzn o zbli偶onych sumach warto艣ci
bezwzgl臋dnych element贸w poszczeg贸lnych wierszy.
Najcz臋stsze skalowanie:
1
D2 = I, D1 = diag(s1, . . . , sn), si = . (349)
n
|aik|
k=1
Alternatywa: modyfikacja zasady wyboru elementu g艂贸wnego w kolumnie
(skalowanie przez si).
Metody Numeryczne II 176
B艂臋dy zaokr膮gle艅 dla rozk艂adu tr贸jk膮tnego
Niech A = LR b臋dzie rozk艂adem tr贸jk膮tnym macierzy A wymiaru n n. L i R
liczymy np. wzorami (330).
Wystarczy oceni膰 wyra偶enia typu:
c - a1b1 - . . . - an-1bn-1
bn = fl . (350)
an
Analiza rozchodzenia si臋 b艂臋d贸w (str. 17) pozwala oszacowa膰:
钆 艂艂
n n-1
eps
c -
鹋偞|sn-1| + (|sj| + |ajbj|) ,
ajbj (351)
1 - eps
j=1 j=1
gdzie
0 je艣li an = 1
s0 = c, sj = fl(sj-1 - ajbj), = (352)
1 w przec. przyp.
Metody Numeryczne II 177
Otrzymujemy dla F = (fij) = A - L R oszacowania:
i-1
i
eps
|fij| = |aij - likrkj| (|ak | + |lik||rkj|) dla j i
ij
k=1
1 - eps
k=1
j-1
j
eps
|fij| = |aij - likrkj| |aj-1| + (|ak | + |lik||rkj|) dla j < i,
ij
k=1 ij
1 - eps
k=1
(353)
k
gdzie ak = fl(aij - lisrsj).
ij
s=1
Przyjmuj膮c ak = max |ak |, a = max ai oraz uwzgl臋dniaj膮c, 偶e rij = ai-1, dla
ij
ij
i,j 0i j i zak艂adaj膮c, 偶e |lij| 1 szacujemy:
eps eps
|fij| (a0 + 2a1 + + 2ai-2 + ai-1) 2(i - 1)a dla j i
1 - eps 1 - eps
eps eps
|fij| (a0 + 2a1 + + 2aj-2 + 2aj-1) 2ja dla j < i
1 - eps 1 - eps
(354)
Metody Numeryczne II 178
钆 艂艂
0 0 0 . . . 0 0
锱1 1 1 . . . 艣艂
1 1
锱 艣艂
锱1 2 2 . . . 艣艂
2 2
eps
锱 艣艂
|fij| 2a (355)
锱1 2 3 . . . 艣艂
3 3
锱 艣艂
1 - eps
锱. . . 艣艂
. .
. .
鹋. . .
. . . . .
1 2 3 . . . n - 1 n - 1
Oszacowanie warto艣ci a
Z definicji a0 = maxi,j aij. Dla cz臋艣ciowego wyboru elementu g艂贸wnego mo偶na
wykaza膰, 偶e ak-1 2ka0. Mamy wi臋c (zwykle pesymistyczne, ale czasem
osi膮galne) oszacowanie
a 2n-1a0. (356)
Dla g贸rnej macierzy Hessenberga (aij = 0 dla i > j + 1) mo偶na wykaza膰, 偶e
a (n - 1)a0. (357)
Dla macierzy tr贸jdiagonalnej
a 2a0. (358)
Metody Numeryczne II 179
Pe艂ny wyb贸r elementu g艂贸wnego
" Wilkinson wykaza艂, 偶e
ak f(k + 1)a0, (359)
gdzie
1
1 1
k-1
2 3
f(k) = k213 4 k . (360)
" W praktyce nie znaleziono kontrprzyk艂adu dla oszacowania
ak (k + 1)a0. (361)
" Bardziej kosztowny ni偶 cz臋艣ciowy.
" Niszczy szczeg贸ln膮 struktur臋 macierzy (np. tr贸jdiagonalnych).
Metody Numeryczne II 180
B艂臋dy zaokr膮gle艅 dla uk艂ad贸w o tr贸jk膮tnych macierzach
Rozwi膮zanie uk艂adu Ax = b sprowadzili艣my do dw贸ch uk艂ad贸w tr贸jk膮tnych:
Ly = b oraz Rx = y.
Korzystamy z oszacowania:
n-1
n
eps
c - aibi
i|aibi| + (n + 1 + )|anbn| . (362)
1 - n eps
i=1 i=1
Dla uk艂adu Ly = b dostajemy:
r r
eps
b - lriyi
i|lriyi| + |yr| , (363)
r
1 - n eps
i=1 i=1
Metody Numeryczne II 181
czyli
钆 艂艂
1 0
锱 艣艂
2
eps
锱 艣艂
|b - Ly| (|L|D - I)|y|, D = . (364)
锱 艣艂
.
.
1 - n eps 鹋
.
0 n
A zatem y jest rozwi膮zaniem dok艂adnym nieco zaburzonego uk艂adu:
eps
(L + "L)y = b, |"L| (|L|D - I). (365)
1 - n eps
Podobne wnioskowanie mo偶na przeprowadzi膰 dla uk艂adu Rx = y, co ostatecznie
prowadzi do wniosku, 偶e liczenie rozwi膮zania uk艂adu metod膮 eliminacji Gaussa jest
algorytmem numerycznie stabilnym.
Metody Numeryczne II 182
Normy macierzy
Definicja 14 Norm膮 macierzy A jest liczba ||A||, gdy || || spe艂nia
nast臋puj膮ce warunki:
1膰% ||A|| 0 oraz ||A|| = 0 wtw. gdy A = 0,
2膰% ||cA|| = |c| ||A|| dla dowolnej c " R,
3膰% ||A + B|| ||A|| + ||B|| (nier贸wno艣膰 tr贸jk膮ta),
4膰% ||AB|| ||A|| ||B|| (nier贸wno艣膰 Schwarza).
Z 4膰% wynika, 偶e
||Am|| ||A||m. (366)
Definicja 15 Norma euklidesowa:
m n
||A||E = a2 . (367)
ij
i=1 j=1
Metody Numeryczne II 183
Definicja 16 Norma spektralna:
||A||S = max |i(AAT )|. (368)
i
Definicja 17 Norma euklidesowa wektora w:
"
||w|| = wT w. (369)
Definicja 18 Promie艅 spektralny macierzy:
:(A) = max |i(A)|. (370)
i
Metody Numeryczne II 184
Wa偶ne w艂asno艣ci:
Dla macierzy kwadratowych:
n
||A||2 = T r(AAT ) = i(AAT ). (371)
E
i=1
St膮d
||A||S ||A||E. (372)
Dla macierzy symetrycznych
||A||S = :(A). (373)
Twierdzenie 31 Dla macierzy kwadratowych A stopnia n
||Ax||
max = ||A||S. (374)
x =0
||x||
Metody Numeryczne II 185
Metody iteracyjne rozwi膮zywania
uk艂ad贸w r贸wna艅 liniowych
" Zwykle dok艂adno艣膰 por贸wnywalna z metodami dok艂adnymi uzyskiwana jest w
d艂u偶szym czasie.
" Dla macierzy rzadkich, metody iteracyjne mog膮 by膰 znacznie szybsze.
" Mniejsza z艂o偶ono艣膰 pami臋ciowa.
Jednopunktowy liniowy wz贸r iteracyjny dla uk艂adu r贸wna艅:
xi+1 = Bixi + ci. (375)
Rozwi膮zujemy uk艂ad r贸wna艅 Ax = b. Inaczej
(I + A)x = x + b. (376)
Narzuca si臋 wz贸r iteracyjny
xi+1 = (I + A)xi - b. (377)
Metody Numeryczne II 186
Chcemy, aby rozwi膮zanie dok艂adne xd by艂o punktem sta艂ym przekszta艂cenia (375):
xd = Bixd + ci. (378)
Skoro
xd = A-1b, (379)
wi臋c dostajemy
A-1b = BiA-1b + ci. (380)
Zatem:
ci = (I - Bi)A-1b = Cib. (381)
Warunek zgodno艣ci dla Bi i Ci:
(I - Bi)A-1 = Ci, czyli Bi + CiA = I (382)
Cz臋stsza posta膰 wzoru iteracyjnego (375):
xi+1 = Bixi + Cib. (383)
Metody Numeryczne II 187
Analiza zbie偶no艣ci wzoru (383)
Wprowadzmy:
i = xi - xd. (384)
Zatem:
i+1 = Bixi + Cib - xd = Bixi + Cib - A-1b =
(385)
= Bixi - BiA-1b = Bii.
A wi臋c:
i+1 = Ki1, gdzie Ki = BiBi-1 . . . B1. (386)
Warunek konieczny i dostateczny zbie偶no艣ci {xi} do xd:
lim Kiy = 0. (387)
"
i "
y
Warunek dostateczny:
lim ||i|| = 0. (388)
i "
Metody Numeryczne II 188
Z twierdzenia 31 mamy
||i+1|| ||Ki1||
max = max = ||Ki||S. (389)
1 1
||1|| ||1||
Warunkiem dostatecznym jest wi臋c
lim ||Ki||S = 0. (390)
i "
Dla metod stacjonarnych mamy r贸wnie偶 warunek dostateczny
lim :(Ki) = 0. (391)
i "
Metody Numeryczne II 189
Stacjonarne metody iteracyjne
W (383) Bi = B, Ci = C, oraz
Ki = BiBi-1 . . . B1 = Bi. (392)
Warto艣ci w艂asne Bi s膮 i-tymi pot臋gami warto艣ci w艂asnych macierzy B. St膮d
gwarancj膮 zbie偶no艣ci jest by modu艂y warto艣ci w艂asnych B by艂y mniejsze od 1 (jest
to r贸wnie偶 warunek konieczny).
Z relacji (372) mi臋dzy normami macierzy dostajemy ostrzejszy warunek
dostateczny:
||B||E < 1. (393)
Szybko艣膰 zbie偶no艣ci:
||i+1|| ||Ki||||1|| ||B||i||1||. (394)
Metody Numeryczne II 190
Metoda prostej iteracji
Niech
A = D + L + U, (395)
gdzie D jest diagonalna, L oraz U tr贸jk膮tnymi doln膮 i g贸rn膮 z zerami na
przek膮tnej.
Uk艂ad Ax = b mo偶na przedstawi膰 w postaci:
Dx = -(L + U)x + b. (396)
Metoda prostej iteracji:
xi+1 = -D-1(L + U)xi + D-1b. (397)
Uwagi:
1. Na przek膮tnej macierzy A nie mog膮 wyst臋powa膰 zera. Dla macierzy
nieosobliwych mo偶na poprzestawia膰 kolumny i wiersze.
2. Sprawdzenie czy B = -D-1(L + U) gwarantuje zbie偶no艣膰 mo偶e by膰 trudne.
Metody Numeryczne II 191
Metoda Seidla
Modyfikacja metody prostej iteracji: po obliczeniu dowolnej sk艂adowej wektora
xi+1, korzysta si臋 z niej przy obliczaniu nast臋pnych sk艂adowych.
Oznaczaj膮c j-t膮 sk艂adow膮 wektora xi przez x(j) otrzymujemy:
i
肱 雠
n
1
砼
x(1) = - a1jx(j) - b1艂艂 , (398)
i+1 i
a11 j=2
肱 雠
n
1
砼俛21x(1) +
x(2) = - a2jx(j) - b2艂艂 . (399)
i+1 i+1 i
a22
j=3
Og贸lnie:
肱 雠
r-1 n
1
砼
x(r) = - arjx(j) + arjx(j) - br艂艂 . (400)
i+1 i+1 i
arr j=1
j=r+1
Metody Numeryczne II 192
W postaci macierzowej:
xi+1 = -D-1(Lxi+1 + Uxi) + D-1b. (401)
St膮d:
xi+1 = -(D + L)-1Uxi + (D + L)-1b. (402)
Warunek dostateczny zbie偶no艣ci:
||B|| = ||(D + L)-1U|| < 1. (403)
Twierdzenie 32 Metoda Seidla jest zbie偶na dla dodatnio okre艣lonych
macierzy A dla dowolnego wektora pocz膮tkowego.
Uwaga: Metody iteracyjne rozwi膮zywania uk艂ad贸w r贸wna艅 liniowych s膮 zbie偶ne
liniowo (p. r贸wnanie (385)). Mo偶na wi臋c zastosowa膰 metod臋 "2.
Metody Numeryczne II 193
Relaksacja
Idea: Rozwa偶amy wektor reszt
r = b - Ax. (404)
Modyfikujemy sk艂adowe przybli偶onego rozwi膮zania x tak, aby wyzerowa膰 jedn膮
lub wi臋cej sk艂adowych wektora r.
Uwaga: Metoda Seidla jest metoda relaksacyjn膮.
Metody Numeryczne II 194
Algorytm wygodny w obliczeniach r臋cznych:
1. Wybierz wektor pocz膮tkowy x1.
2. Dla i od 1 w g贸r臋 tak d艂ugo jak trzeba
(a) Oblicz ri = b - Axi.
(L)
(b) Odszukaj sk艂adow膮 wektora ri o najwi臋kszym module (ri ) i element o
najwi臋kszym module z L-tym wierszu macierzy A (aLJ)
(c) Wyznacz wektor xi+1 maj膮cy wszystkie sk艂adowe identyczne ze
sk艂adowymi wektora xi, z wyj膮tkiem
(L)
ri
x(J) = x(J) + (405)
i+1 i
aLJ
Uwagi: ri+1 lepiej wyznacza膰 jako
(L)
ri AL
ri + (x(J) - x(J) )AL = ri - , gdzie AL to L-ta kolumna A. (406)
i i+1
aLJ
Je艣li elementy maksymalne w wierszu s膮 na g艂贸wnej przek膮tnej, to J = L.
Metody Numeryczne II 195
Metoda nadrelaksacji
Szybko艣膰 zbie偶no艣ci zale偶y od promienia spektralnego macierzy B zmniejszenie
promienia spektralnego oznacza przyspieszenie zbie偶no艣ci.
Przekszta艂camy uk艂ad Ax = b tak, by na przek膮tnej macierzy A dosta膰 same
jedynki. Metoda Seidla (400) przybiera posta膰:
r-1 n
(j)
x(r) = - arjx(j) - arjxi + br
i+1 i+1
j=1 j=r+1
(407)
r-1 n
= x(r) - arjx(j) - arjx(j) + br.
i i+1 i
j=1 j=r
Rozwa偶amy
肱 雠
r-1 n
砼
x(r) = x(r) - arjx(j) + arjx(j) - br艂艂 , (408)
i+1 i i+1 i
j=1 j=r
gdzie jest czynnikiem nadrelaksacji.
Metody Numeryczne II 196
W postaci macierzowej:
xi+1 = Bxi + Cb, (409)
gdzie
B = (I + L)-1((1 - )I - U), C = (I + L)-1. (410)
Problem: Jak dobra膰 , by promie艅 :(B) by艂 minimalny?
Og贸lnie: trudne zadanie.
Czasami: np. prosta zale偶no艣膰 mi臋dzy optymalnym a promieniem :(L + U).
Metody Numeryczne II 197
Metody iteracyjne wynikaj膮ce
z minimalizacji formy kwadratowej
Cel: Rozwi膮zanie uk艂adu Ax = b poprzez minimalizacj臋 formy kwadratowej.
Je艣li A nie jest symetryczna i dodatnio okre艣lona, mo偶na minimalizowa膰
R = rT r = (b - Ax)T (b - Ax). (411)
Zwykle ogromna ilo艣膰 oblicze艅.
Je艣li A jest symetryczna i dodatnio okre艣lona, mo偶na minimalizowa膰
1
Q = xT Ax - xT b. (412)
2
Minimum w rozwi膮zaniu Ax = b bo:
1 1
Q(x + "x) - Q(x) = (x + "x)T A(x + "x) - (x + "x)T b - xT Ax - xT b
2 2
1 1
= "xT Ax - "xT b + "xT A"x = "xT A"x.
2 2
Metody Numeryczne II 198
Og贸lna metoda minimalizacji
W kolejnym kroku maj膮c wektor xi, wybieramy kierunek vi i odleg艂o艣膰 膮i:
xi+1 = xi + 膮ivi. (413)
Metoda najwi臋kszego spadku
Kierunek vi wybieramy jako kierunek najszybszego spadku formy Q czyli kierunek
gradientu w punkcie xi:
vi ! "Q = Axi - b = -ri. (414)
Aby wyznaczy膰 膮i liczymy:
1 1 1
2
Q(xi - 膮iri) = - xT ri + 膮irT ri + 膮i rT Ari - xT b. (415)
i i i i
2 2 2
Zatem
Q
= rT ri + 膮irT Ari, (416)
i i
膮i
Metody Numeryczne II 199
czyli minimum w punkcie
rT r
i
膮i = - . (417)
rT Ar
i
Ostatecznie
rT r
i
xi+1 = xi + ri. (418)
rT Ar
i
Uwagi:
" Zbie偶no艣膰 zapewniona, bo Q(xi+1) < Q(xi).
" Zbie偶no艣膰 mo偶e by膰 bardzo wolna dla uk艂ad贸w zle uwarunkowanych.
Metody Numeryczne II 200
Metoda sprz臋偶onych gradient贸w
Idea: Odpowiedniejszy dob贸r kolejnych kierunk贸w.
Niech v1, . . . , vn b臋d膮 baz膮 n wymiarowej przestrzeni euklidesowej i niech x
b臋dzie przybli偶eniem rozwi膮zania dok艂adnego xd. Wtedy:
n
xd - x = 膮jvj. (419)
j=1
Je艣li wektory vj s膮 ortogonalne, to
T
vj (xd - x)
膮j = , j = 1, . . . , n. (420)
T
vj vj
Je艣li wektory vj s膮 A-ortogonalne (lub A-sprz臋偶one), tzn.
T
vj Avi = 0 dla i = j, (421)
to
T T
vj A(xd - x) vj r
膮j = = , j = 1, . . . , n. (422)
T T
vj Avj vj Avj
Metody Numeryczne II 201
Metoda sprz臋偶onych gradient贸w polega na wyznaczaniu wektor贸w vj i
wsp贸艂czynnik贸w 膮j, by ostatecznie wyznaczy膰 xd.
Konstrukcja wektor贸w A-ortogonalnych: za艂贸偶my, 偶e mamy baz臋 wektor贸w
u1, . . . , un. Stosujemy metod臋 ortogonalizacji Grama-Schmidta:
i
v1 = u1, vi+1 = ui+1 + i+1,kvk, i = 1, . . . , n, (423)
k=1
gdzie
T
vk Aui+1
i+1,k = - . (424)
T
vk Avk
Metody Numeryczne II 202
Pozostaje wyznaczy膰 ci膮g wektor贸w {ui}
" Wektory jednostkowe.
" Bardziej wygodny spos贸b: za艂贸偶my, 偶e
i-1
xi = x1 + 膮jvj, (425)
j=1
gdzie 膮j spe艂niaj膮 (422). Zatem
xi+1 = xi + 膮ivi, (426)
ri = b - Axi. (427)
Przyjmujemy
ui = ri, i = 1, . . . , n. (428)
Metody Numeryczne II 203
Rekurencyjne wyznaczanie ri, 膮i oraz vi definiuje metod臋 sprz臋偶onych
gradient贸w:
v1 = r1 = b - Ax1
T
vi ri
膮i =
T
vi Avi
xi+1 = xi + 膮ivi
dla i = 1, . . . , n. (429)
ri+1 = ri - 膮iAvi
T
vi Ari+1
i = -
T
vi Avi
vi+1 = ri+1 + ivi
Uwagi:
" Na u偶ytek vi wyznaczamy tylko jedno i.
" xn+1 = xd metoda sko艅czona, cho膰 b艂臋dy zaokr膮gle艅 psuj膮 t臋 r贸wno艣膰.
" Z艂o偶ono艣膰 O(n3). Troch臋 wi臋cej operacji ni偶 eliminacja Gaussa. Korzystniej
dla macierzy rzadkich, ale wtedy skuteczniejsze prostsze metody (np. Seidla).
Metody Numeryczne II 204
Metoda ortogonalizacji Householdera
Idea: Kolejne transformacje macierzy dobierane tak, by wskaznik uwarunkowania
zadania nadmiernie nie wzr贸s艂.
Analiza prowadzi do wniosku: przekszta艂cenia o macierzach unitarnych zachowuj膮
wskaznik uwarunkowania.
Propozycja Householdera: unitarna macierz transformacji P spe艂niaj膮ca
P = I - 2wwH, gdzie wHw = 1, w " Cn. (430)
Macierz P jest hermitowska:
PH = I - (2wwH)H = I - 2wwH = P, (431)
i unitarna:
PHP = PP = P2 = (I - 2wwH)(I - 2wwH)
(432)
= I - 2wwH - 2wwH + 4wwHwwH = I.
Metody Numeryczne II 205
Wektor w dobierany jest tak, aby pewien wektor x = [x1, . . . , xn]T by艂
przekszta艂cony na kierunek wektora jednostkowego e1:
Px = ke1. (433)
Dok艂adniejsza analiza prowadzi do:
uuH
P = I - , (434)
艂
gdzie
n
= |xi|2, x1 = ei膮|x1|, k = -ei膮,
i=1
钆 艂艂
ei膮(|x1| + )
(435)
锱 艣艂
x2
锱 艣艂
u = x - ke1 = , 艂 = ( + |x1|).
锱 艣艂
.
.
鹋
.
xn
Metody Numeryczne II 206
Algorytm:
Przekszta艂camy macierz A krok po kroku zeruj膮c kolejne kolumny pod przek膮tn膮:
钆 艂艂
r11 . . . r1n
锱 艣艂
.
.
. .
A = A(0) A(1) A(n-1) = R = . (436)
鹋
.
.
0 rnn
W j=tym kroku przekszta艂camy macierz
钆 艂艂
r11 . . . x1,j-1 a(j-1) . . . a(j-1)
1j 1n
锱 艣艂
. . .
.
锱 艣艂
. . . .
.
. . .
锱 艣艂
锱 艣艂
锱 艣艂
0 rj-1,j-1 a(j-1) . . . a(j-1)
D B
j-1,j j-1,n
锱 艣艂
A(j-1) = = . (437)
锱
0 (j-1)
a(j-1) . . . a(j-1) 艣艂
锱 艣艂
jj jn
锱 艣艂
. .
锱 艣艂
. .
0 . .
鹋
a(j-1) . . . a(j-1)
nn
nj
Metody Numeryczne II 207
Szukamy takiej unitarnej macierzy przekszta艂cenia Pj, aby
钆 艂艂
a(j-1)
jj
锱 艣艂
Ij-1 0
艣艂
.
Pj = , Pj 锱 . = ke1 " Cn-j+1. (438)
.
鹋
0 Pj
a(j-1)
nj
Uwagi:
" Mno偶enie przez macierz
ujuH
j
Pj = I - (439)
艂j
realizuje si臋 w nast臋puj膮cy spos贸b:
uH(j-1)
j
H H
Pj(j-1) = (j-1) - ujyj , gdzie yj = . (440)
艂j
2n3
" Algorytm wymaga oko艂o dzia艂a艅.
3
" Wsp贸艂rz臋dne wektora u wpisuje si臋 w macierz A pod przek膮tn膮 brakuje
tylko miejsca na jeden wektor n warto艣ci.
" Kolumny macierzy Q = P-1 = P1 Pn-1 s膮 ortonormalne
Metody Numeryczne II 208
Metoda Householdera przyk艂ad
Rozwi膮za膰 uk艂ad r贸wna艅
钆 艂艂 钆 艂艂
2 4 -1 0
鹋2 -1 2 鹋6
x = (441)
1 -3 1 3
sprowadzaj膮c macierz g艂贸wn膮 do postaci tr贸jk膮tnej metod膮 Householdera.
Metody Numeryczne II 209
Metoda ortogonalizacji Grama-Schmidta
Je艣li macierz A ma rozk艂ad:
A = QR, (442)
gdzie macierz Q ma ortonormalne kolumny, to mo偶na je wyznaczy膰 metod膮
ortogonalizacji Grama-Schmidta.
Dla k-tej kolumny macierzy A = (a1, . . . , an) mamy:
k
ak = rikqi, k = 1, . . . , n. (443)
i=0
Kolumna ak jest kombinacj膮 liniow膮 wektor贸w q1, . . . , qk i odwrotnie: wektor qk
jest kombinacj膮 liniow膮 a1, . . . , ak.
Kolejne qk wyznaczamy rekurencyjnie tak, by zachowywa膰 ortonormalno艣膰.
Metody Numeryczne II 210
Algorytm:
" Wyznaczamy q1 i r11:
a1
r11 ! ||a1||, q1 ! . (444)
r11
" Znamy q1, . . . , qk-1 oraz rij dla j k - 1.
Wyznaczamy wektor
bk ! ak - r1kq1 - . . . - rk-1,kqk-1 (445)
tak, aby by艂 ortogonalny do q1, . . . , qk-1:
H H
qi bk = 0 ! rik ! qi ak, i = 1, . . . , k - 1 (446)
Wyznaczamy qk i rkk:
bk
rkk ! ||bk||, qk ! . (447)
rkk
Metody Numeryczne II 211
Uwagi:
" Dla wszystkich 1 i, j k mamy:
1 dla i = j,
H
qi qj = (448)
0 w przec. przyp.
" Algorytm jest niestabilny numerycznie gdy wektory ai s膮 bliskie liniowej
zale偶no艣ci b艂臋dy numeryczne sprawiaj膮, 偶e bk nie jest ortogonalny do qi.
Reortogonalizacja wektora bk
Obliczamy poprawki w postaci skalar贸w
H
"rik ! qi bk, i = 1, . . . , k - 1, (449)
oraz wektor
bk = bk - "r1kq1 - . . . - "rk-1,kqk-1. (450)
Poprawiamy:
bk
rkk ! || qk ! , rik ! rik + "rik. (451)
bk||,
rkk
Metody Numeryczne II 212
Wyznaczanie warto艣ci i wektor贸w w艂asnych
Definicja 19 Wektor x jest wektorem w艂asnym macierzy A je艣li istnieje
taka liczba , 偶e Ax = x. to warto艣膰 w艂asna macierzy A.
Twierdzenie 33 Liczba jest warto艣ci膮 w艂asn膮 macierzy A wtedy i tylko
wtedy, gdy jest pierwiastkiem wielomianu charakterystycznego
det(A - I) macierzy A.
Definicja 20 Macierze A i B s膮 podobne je艣li istnieje nieosobliwa macierz
podobie艅stwa P tj. taka, 偶e P-1AP = B.
Twierdzenie 34 Macierze podobne maj膮 takie same warto艣ci w艂asne.
Twierdzenie 35 Warto艣ci i wektory w艂asne macierzy symetrycznej s膮
rzeczywiste.
Twierdzenie 36 Warto艣ci w艂asne macierzy A + cI to warto艣ci w艂asne
macierzy A przesuni臋te o c.
Metody Numeryczne II 213
Metoda Kry艂owa
" Szukanie warto艣ci w艂asnych jako pierwiastk贸w wielomianu
charakterystycznego (1931 r.).
" Zagadnienie wyznaczenia wsp贸艂czynnik贸w wielomianu charakterystycznego
jest znacznie gorzej uwarunkowane ni偶 zadanie wyznaczenia warto艣ci
w艂asnych.
" W zasadzie metoda nie do u偶ytku.
Metody Numeryczne II 214
Metoda pot臋gowa
Dla wyznaczania pojedynczych warto艣ci i wektor贸w w艂asnych.
Algorytm:
" i ! 0, wybieramy wektor x0 taki, 偶e ||x0|| " = 1
" Tak d艂ugo jak i jest mniejsze ni偶 za艂o偶ona liczba iteracji:
vi+1 ! Axi, mi+1 ! ||vi+1|| "
Je偶eli mi+1 = 0 to KONIEC!
vi+1
xi+1 ! , i ! i + 1
mi+1
Uwagi:
" Je艣li x2i i " x, to mi i " m.
- -
" Je艣li ||xi+1 - xi|| " i " 0, to Ax = mx.
-
" Je艣li ||xi+1 - xi|| " i " 2, to Ax = -mx.
-
" W og贸lnym przypadku trudno oceni膰 kiedy ci膮g (x2i) jest zbie偶ny.
" Mo偶na wykaza膰 zbie偶no艣膰 np. gdy macierz jest podobna do diagonalnej.
Metody Numeryczne II 215
Metoda QR
Dla macierzy Hessenberga H (hij = 0 dla j < i - 1).
Schemat algorytmu:
" i ! 1, H(1) ! H
" Do uzyskania stosownej zbie偶no艣ci:
Wyznaczamy rozk艂ad H(i) = Q(i)R(i)
H(i+1) ! R(i)Q(i)
i ! i + 1
Uwagi:
" Wszystkie Q(i) oraz H(i) s膮 macierzami Hessenberga.
T
" Wszystkie H(i) s膮 podobne, bo H(i+1) = Q(i) H(i)Q(i).
" Je艣li warto艣ci w艂asne H s膮 rzeczywiste i maj膮 r贸偶ne modu艂y, to H(i) zbiegaj膮
do macierzy tr贸jk膮tnej, a elementy diagonali do warto艣ci w艂asnych H.
" Je艣li w艣r贸d warto艣ci w艂asnych s膮 pary sprz臋偶one (poza tym nie ma warto艣ci o
r贸wnych modu艂ach), to nie wszystkie elementy pod diagonal膮 d膮偶膮 do zera.
Metody Numeryczne II 216
Problemy zbie偶no艣ci
" Analiza szybko艣ci zbie偶no艣ci pokazuje, 偶e zale偶y ona od stosunku modu艂贸w
odpowiednich warto艣ci w艂asnych.
" Dotyczy zar贸wno metody pot臋gowej jak i QR.
" Spos贸b na przyspieszenie: przesuni臋cie B = A - pI.
Metoda pot臋gowa, gdy znamy warto艣ci w艂asne, a szukamy wektor贸w:
1 1
p = (i + n) dla 1, gdzie i = min{k > 1 : k < 1} oraz p = (1 + j)
2 2
dla n, gdzie j = max{k < n : k > n}.
Metoda QR:
H(i) - piI = Q(i)R(i) H(i+1) = R(i)Q(i) + piI,
pi r贸wne otrzymanym w poprzednich iteracjach warto艣ciom w艂asnym.
Sprowadzanie macierzy do postaci Hessenberga
" eliminacja Gaussa
" metoda Householdera
" metoda Givensa
Metody Numeryczne II 217
Szukanie wektor贸w w艂asnych
" Je偶eli P-1AP = B, oraz x jest wektorem w艂asnym macierzy B, to y = Px
jest wektorem w艂asnym macierzy A.
" Je艣li macierz B jest tr贸jk膮tna g贸rna, to warto艣ci w艂asnej i = bii odpowiada
wektor w艂asny x(i) taki, 偶e:
x(i) = 0, j = n, n - 1, . . . , i + 1
j
x(i) = 1
i
i
(452)
bjkx(i)
k
k=j+1
x(i) = - , j = i - 1, . . . , 1
j
bjj - bii
Je艣li B ma wielokrotn膮 warto艣膰 w艂asn膮 , to wektor w艂asny mo偶na
wyliczy膰 tylko dla najmniejszego i takiego, 偶e bii = .
Metody Numeryczne II 218
Metoda QR
" Macierze Q(i) mo偶na wykorzysta膰 do obliczenia wektor贸w w艂asnych.
" Algorytm: Dane: macierz A.
Sprowadzamy macierz A do postaci Hessenberga
P-1AP = H. (453)
Obliczamy ci膮g macierzy H(i) oraz Q(i) algorytmem QR dla macierzy
Hessenberga H. Jednocze艣nie wyznaczamy
Z(1) = P, Z(i+1) = Z(i)Q(i). (454)
Po uznaniu w H(N) zbie偶no艣ci do macierzy tr贸jk膮tnej g贸rnej wyznaczamy
jej wektory w艂asne zgodnie z (452).
Wyznaczamy wektory w艂asne macierzy A:
y(k) = Z(N)x(k). (455)
Metody Numeryczne II 219
Metoda LR
" Algorytm: Dane: macierz A.
Obliczamy ci膮g macierzy A(i) dla i = 1, 2, . . .:
A(1) ! A
(456)
A(i) = L(i)U(i), A(i+1) ! U(i)L(i)
do uzyskania zbie偶no艣ci do macierzy tr贸jk膮tnej.
Jednocze艣nie wyznaczamy:
Z(1) = I, Z(i+1) = Z(i)L(i). (457)
Wyznaczamy wektory w艂asne macierzy A(N) zgodnie z (452).
Wyznaczamy wektory w艂asne macierzy A:
y(k) = Z(N)x(k). (458)
" 3 razy mniej dzia艂a艅 ni偶 w QR ale gorsze uwarunkowanie.
" Stabilna numerycznie dla A symetrycznych i dodatnio okre艣lonych.
" Jak w QR: przesuni臋cia i do postaci Hessenberga (O(n3) O(n2)).
Metody Numeryczne II 220
Metoda iteracji odwrotnej
" Przyspieszenie zbie偶no艣ci metody pot臋gowej
" Algorytm:
i ! 0, wybieramy wektor x0 taki, 偶e ||x0|| " = 1
Tak d艂ugo jak i jest mniejsze ni偶 za艂o偶ona liczba iteracji:
" vi+1 ! (A - kiI)-1xi, mi+1 ! ||vi+1|| "
vi+1
" xi+1 ! , i ! i + 1
mi+1
" Je艣li A ma warto艣ci w艂asne 1, . . . , n takie, 偶e |1| . . . |n|, to
warto艣ciami w艂asnymi (A - kiI)-1 s膮:
1 1 1
, , . . . , . (459)
1 - ki 2 - ki n - ki
" Ci膮g m1, m2, . . . pozwala wyznaczy膰 warto艣膰 w艂asn膮 n jak w metodzie
pot臋gowej.
" Je艣li ki jest bliskie n, to macierz A - kiI jest zle uwarunkowana, co mo偶na
zlekcewa偶y膰 gdy x jest dobrze uwarunkowany i norma ||vi+1|| " jest du偶a.
Wyszukiwarka
Podobne podstrony:
MetNum 2012
MetNum Lab09
MetNum wyk (2)
MetNum Lab03
wi臋cej podobnych podstron