Metody Numeryczne II
1
Te materiały można znaleźć 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¨orck. 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. WPŁ
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. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2. Interpolacja. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3. Całkowanie numeryczne. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4. Aproksymacja. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5. Rozwiązywanie równań nieliniowych. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6. Rozwiązywanie układów równań liniowych. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7. Wyznaczanie wartości i wektorów własnych macierzy. . . . . . . . . . . . . . . . . . . .
8. Rozwiązywanie równań różniczkowych zwyczajnych.
9. Szybka transformata Fouriera.
10. Generatory liczb (pseudo)losowych. Metoda Monte Carlo.
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:
e = 1 +
1
1!
+
1
2!
+
1
3!
+ . . .
(1)
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
∗ 2
2
+ 1
∗ 2
0
=
1
∗ 4 + 1 = 5
ósemkowy
0 1 2 3 4 5 6 7
270
2
∗ 8
2
+ 7
∗ 8
1
=
2
∗ 64 + 7 ∗ 8 = 184
dziesiętny
0 1 2 3 4
309
3
∗ 10
2
+ 9
∗ 10
0
=
5 6 7 8 9
3
∗ 100 + 9 = 309
szesnastkowy
0 1 2 3 4 5 6 7
10A
1
∗ 16
2
+ 10
∗ 16
0
=
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
∗ 10
2
= 0.100101
∗ 2
101
– 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
∗ 10
2
, nie 0.125
∗ 10
3
ani 125
∗ 10
0
– 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 = 2
4
+ 2
2
+ 2
1
+ 2
0
+ 2
−1
+ 2
−2
• po normalizacji: 1.011111 ∗ 2
4
= 1.011111
∗ 2
100
• 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ądź 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
∗ 2
w
−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 ∗ 2
w
−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
n /
±Infinity
0
±Infinity x ±Infinity ±Infinity
±nonzero / 0
±Infinity
Infinity + Infinity
Infinity
Operacja
Wynik
±0 / ±0
NaN
Infinity - Infinity
NaN
±Infinity / ±Infinity
NaN
±Infinity x 0
NaN
Zakresy dwójkowo:
precyzja
zdenormalizowane
znormalizowane
pojedyncza
±2
−149
. . . (1
− 2
−23
)
∗ 2
−126
±2
−126
. . . (2
− 2
−23
)
∗ 2
127
podwójna
±2
−1074
. . . (1
− 2
−52
)
∗ 2
−1022
±2
−1022
. . . (2
− 2
−52
)
∗ 2
1023
Zakresy dziesiętnie:
precyzja
dziesiętnie
pojedyncza
± ∼ 10
−44.85
. . .
∼ 10
38.53
podwójna
± ∼ 10
−323.3
. . .
∼ 10
308.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|
x
(3)
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
˜
m =
α
1
.α
2
. . . α
k
jeżeli α
k+1
<
B
2
α
1
.α
2
. . . α
k
+ B
−k+1
jeżeli α
k+1
B
2
(5)
a liczbę znormalizowaną x = m
∗ B
c
do ˜
rd(x) = sgn(x) ˜
m
∗ B
c
.
Błąd względny przybliżenia:
|
˜
rd(x)
− x
x
|
B
2
∗ B
−k
|m|
1
2
∗ B
−k+1
(6)
• dokładność maszynowa eps
def
=
1
2
∗ B
−k+1
Zatem: ˜
rd(x) = x(1 + ), gdzie
|| 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:
x +
∗
y
def
= rd(x + y)
x
−
∗
y
def
= rd(x
− y)
x
×
∗
y
def
= rd(x
× y)
x/
∗
y
def
= rd(x/y)
Wówczas:
x +
∗
y = (x + y)(1 +
1
)
x
−
∗
y = (x
− y)(1 +
2
)
x
×
∗
y = (x
× y)(1 +
3
)
x/
∗
y = (x/y)(1 +
4
),
gdzie
|
i
| eps
Metody Numeryczne II
13
Problemy arytmetyki zmiennopozycyjnej
• Jeżeli |y| <
eps
B
|x|, to x +
∗
y = x
• Działania zmiennopozycyjne nie muszą być łączne lub rozdzielne.
Przykład:
a = 2.3371258
10
− 4, b = 3.3678429
10
2, c =
−3.3677811
10
2
a + b + c = 6.41371258
10
− 3
a +
∗
(b +
∗
c) = 2.3371258
10
− 4 +
∗
6.1800000
10
− 3 = 6.4137126
10
− 3
(a +
∗
b) +
∗
c = 3.3678452
10
2 +
∗
−3.3677811
10
2 = 6.4100000
10
− 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
→ R
m
, gdzie D
⊆ R
n
, n, m
∈ N. Szukamy
metody wyznaczania wartości y = φ(x).
Odwzorowaniem elementarnym nazywamy odwzorowanie ψ : D
⊆ R
k
→ R
l
takie, że jego funkcje składowe ψ
i
są działaniami elementarnymi.
Algorytmem realizującym funkcję φ : D
→ R
m
(D
⊆ R
n
, n, m
∈ N) nazywamy
taki ciąg odwzorowań φ
(1)
, . . . , φ
(r)
, (gdzie φ
(i)
: D
i
→ D
i+1
, D
i
⊆ R
n
i
, D
1
= D,
n
r+1
= m), że
φ = φ
(1)
◦ · · · ◦ φ
(r)
(7)
Przykład:
Zadanie φ(a, b, c) = a + b + c może być realizowane algorytmem
{φ
(1)
, φ
(2)
}, gdzie
φ
(1)
(a, b, c) =
a + b
c
,
φ
(2)
(d, e) = d + e.
(8)
Metody Numeryczne II
15
Propagacja błędów
Prześledźmy błędy zaokrągleń dla przykładu sumowania (a + b) + c:
f l((a + b) + c) = f l(f l(a + b) + c) =
((a + b)(1 +
1
) + c)(1 +
2
) =
(a + b + c)
1 +
a + b
a + b + c
1
(1 +
2
) +
2
a+b
a+b+c
jest współczynnikiem wzmocnienia błędu
1
– jego duża wartość powoduje
duży błąd względny
Ponieważ we wcześniejszym przykładzie:
a + b
a + b + c
≈ 5 ∗ 10
4
,
b + c
a + b + c
≈ 0.97
(9)
to wybrać należy metodę f l(a + (b + c)).
Metody Numeryczne II
16
Zadanie: Analizujemy „naturalny” algorytm znajdujący przybliżenie b
n
rozwiązania β
n
równania liniowego
c
− a
1
b
1
− · · · − a
n
−1
b
n
−1
− a
n
β
n
= 0,
(10)
dla danych c, a
1
, . . . , a
n
, , b
1
, . . . , b
n
−1
, a
n
= 0. Oszacować residuum
r = c
− a
1
b
1
− · · · − a
n
b
n
.
(11)
Algorytm: Wyznaczamy
b
n
= f l
c
− a
1
b
1
− · · · − a
n
−1
b
n
−1
a
n
(12)
poprzez:
s
0
= c,
s
j
= f l(s
j
−1
− a
j
b
j
) = [s
j
−1
− a
j
b
j
(1 + µ
j
)] (1 + α
j
),
j = 1, 2, . . . , n
− 1,
b
n
= f l
s
n
−1
a
n
= (1 + δ)
s
n
−1
a
n
.
(13)
Metody Numeryczne II
17
Oszacowanie 1: Z (
) wynikają (drugie równanie dla j = 1, 2, . . . , n
− 1):
s
0
− c = 0,
s
j
− (s
j
−1
− a
j
b
j
) = s
j
−
s
j
1 + α
j
+ s
j
b
j
µ
j
= s
j
α
j
1 + α
j
− s
j
b
j
µ
j
,
a
n
b
n
− s
n
−1
= δs
n
−1
.
(14)
Dodając stronami dostajemy:
r = c
−
n
i=1
a
i
b
i
=
n
−1
j=1
−s
j
α
j
1 + α
j
+ a
j
b
j
µ
j
− δs
n
−1
.
(15)
Stąd oszacowanie:
|r|
eps
1
− eps
δ
|s
n
−1
| +
n
−1
j=1
(
|s
j
| + |a
j
b
j
|)
,
δ
=
0
jeżeli a
n
= 1
1
w przec. przyp.
(16)
Metody Numeryczne II
18
Oszacowanie 2: Z (
):
b
n
=
c
n
−1
k=1
(1 + α
k
)
−
n
−1
j=1
a
j
b
j
(1 + µ
j
)
n
−1
k=j
(1 + α
k
)
1 + δ
a
n
.
(17)
Stąd
c =
n
−1
j=1
a
j
b
j
(1 + µ
j
)
j
−1
k=j
(1 + α
k
)
−1
+ a
n
b
n
(1 + δ)
−1
n
−1
k=1
(1 + α
k
)
−1
(18)
Łatwo pokazać indukcyjnie (ze względu na m), że z
1 + σ =
m
k=1
(1 + σ
k
)
±1
,
|σ
k
| eps, m · eps < 1
(19)
wynika
|σ|
m
· eps
1
− m · eps
.
(20)
Metody Numeryczne II
19
Zatem istnieją wielkości ε
j
takie, że
c =
n
−1
j=1
a
j
b
j
(1 + jε
j
) + a
n
b
n
(1 + (n
− 1 + δ
)ε
n
),
|ε
j
|
eps
1
− n · eps
,
δ
=
0
jeżeli a
n
= 1
|
1
w przec. przyp.
(21)
Ostatecznie dla r = c
− a
1
b
1
− · · · − a
n
b
n
dostajemy:
|r|
eps
1
− n · eps
n
−1
j=1
j
|a
j
b
j
| + (n − 1 + δ
)
|a
n
b
n
|
.
(22)
Metody Numeryczne II
20
Wskaźniki uwarunkowania
Niech φ : D
→ R
m
dla otwartego zbioru D
⊆ R
n
, oraz φ =
φ
1
..
.
φ
m
przy czym
φ
i
mają ciągłe pochodne na D. Przy założeniu pomiarów niedokładnych ˜
x
wielkości x otrzymujemy ˜
y = φ(˜
x) jako przybliżenie y = φ(x). Rozwinięcie w
szereg Taylora z pominięciem wyrazów wyższych rzędów daje:
∆y
i
= ˜
y
i
− y
i
= φ
i
(˜
x)
− φ
i
(x)
.
=
n
j=1
(˜
x
j
− x
j
)
δφ
i
(x)
δx
j
=
n
j=1
δφ
i
(x)
δx
j
∆x
j
(23)
W związku z tym błąd względny:
∆y
i
y
i
.
=
n
j=1
x
j
φ
i
(x)
δφ
i
(x)
δx
j
∆x
j
x
j
(24)
Współczynniki
x
j
φ
i
(x)
δφ
i
(x)
δx
j
nazywamy wskaźnikami uwarunkowania zadania.
Metody Numeryczne II
21
Zadania źle i dobrze uwarunkowane
Zadanie nazywamy źle uwarunkowanym jeżeli przynajmniej jeden ze wskaźników
uwarunkowania ma dużą wartość bezwzględną. W przeciwnym razie zadanie
nazywamy dobrze uwarunkowanym.
• tylko dla niezerowych x
j
i y
i
• mn wskaźników
• inne definicje, np: c jest wskaźnikiem uwarunkowania zadania φ jeśli
||φ(˜x) − φ(x)||
||φ(x)||
c
||˜x − x||
||x||
(25)
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ł ˜a = [a
, a
], gdzie a
, a
∈ A.
• Działania na przedziałach: dla ∈ {⊕, , ⊗, } wynik
˜
c = ˜
a
˜b ⊇ {x · y|x ∈ ˜a, y ∈ ˜b}
(27)
Przykłady:
[a
, a
]
⊕ [b
, b
] = [max
{x ∈ A|x a
+ b
}, min{x ∈ A|x a
+ b
}]
[a
, a
]
⊗ [b
, b
] = [max
{x ∈ A|x a
× b
}, min{x ∈ A|x a
× b
}]
Uwaga: Oszacowania przedziałowe są prawdziwe, ale dość pesymistyczne
Przykład: Schemat Hornera a wzory skróconego mnożenia dla
x
3
− 3x
2
+ 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 = ˜
y = [
−1, 1].
Wówczas ˜
z = [
−2, 2]. Ale rozkład prawdopodobieństwa ˜z nie jest jednostajny.
−n
n
1/n
arytm.przedz.
n = 2
n = 10
n = 30
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
prawdopodobieństwem
1
2
mamy
|| < x.
– Dla rozkładu normalnego mamy więc x = 0.6745σ, gdzie σ jest
odchyleniem standardowym tego rozkładu.
– W poprzednim przykładzie mamy σ =
n
3
, a więc:
n = 1
⇒ σ ≈ 0.577
n = 10
⇒ σ ≈ 1.82 = 0.182n
n = 30
⇒ σ ≈ 3.16 ≈ 0.1n
Metody Numeryczne II
26
Twierdzenie 1 (o błędzie standardowym) Jeżeli błędy ∆x
1
, . . . , ∆x
n
są
niezależnymi zmiennymi losowymi o wartości oczekiwanej równej zeru i
odchyleniach standardowych odpowiednio
1
, . . . ,
n
, to błąd standardowy dla
y = y(x
1
, . . . , x
n
) dany jest wzorem:
≈
n
i=1
δy
δx
i
2
2
i
.
(29)
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; a
0
, . . . , a
n
) 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 (x
i
, f
i
), i = 0, . . . , m,
gdzie x
i
= x
j
dla i
= j, polega na wyznaczeniu parametrów a
i
tak, aby
Φ(x
i
; a
0
, . . . , a
n
) = f
i
,
i = 0, . . . , m.
(30)
Pary (x
i
, f
i
) nazywamy punktami węzłowymi a wartości x
i
i f
i
odpowiednio
odciętą i rzędną węzła.
Metody Numeryczne II
28
Zadania interpolacji
• i. wielomianowa
Φ(x
i
; a
0
, . . . , a
n
) = a
0
+ a
1
x + . . . + a
n
x
n
(31)
• i. trygonometryczna
Φ(x
i
; a
0
, . . . , a
n
) = a
0
+ a
1
e
xi
+ . . . + a
n
e
nxi
(32)
• i. funkcjami sklejanymi
• i. wymierna
Φ(x
i
; a
0
, . . . , a
n
, b
0
, . . . , b
n
) =
a
0
+ a
1
x + . . . + a
n
x
n
b
0
+ b
1
x + . . . + b
n
x
n
(33)
• i. eksponencjalna
Φ(x
i
; a
0
, . . . , a
n
; λ
0
, . . . , λ
n
) = a
0
e
λ
0
x
+ a
1
e
λ
1
x
+ . . . + a
n
e
λ
n
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) = a
0
+ a
1
x + . . . + a
n
x
n
)
Twierdzenie 2 Dla n+1 dowolnych punktów węzłowych (x
i
, f
i
), i = 0, . . . , n,
takich, że x
i
= x
j
dla i
= j, istnieje tylko jeden wielomian p ∈ Π
n
spełniający
P (x
i
) = f
i
,
i = 0, . . . , n
(35)
Dowód:
Jednoznaczności: Niech P
1
, P
2
∈ Π
n
spełniają warunki (
). Wówczas
P = P
1
− P
2
∈ Π
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:
P (x) =
n
i=0
f
i
L
i
(x)
(36)
gdzie L
i
są tak zwanymi wielomianami Lagrange’a zdefiniowanymi jako:
L
i
(x) =
n
Π
j=0
(j
=i)
x
− x
j
x
i
− x
j
,
(37)
i spełniającymi następujące warunki (delty Kroneckera):
L
i
(x
j
) = δ
ij
=
1,
gdy i = j
0,
gdy i
= j
.
(38)
Uwagi:
• Nieduża przydatność w praktyce (n
2
− 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 (x
i
, f
i
), i = 0, . . . , n przez
P
i
0
,...,i
k
oznaczamy wielomian ze zbioru Π
k
taki, że
P
i
0
...i
k
(x
i
j
) = f
i
j
,
j = 0, . . . , k.
(39)
Twierdzenie 3 Wielomiany typu (
) mają następujące własności:
1
◦
P
i
(x) = f
i
2
◦
P
i
0
...i
k
(x) =
(x
− x
i
0
)P
i
1
...i
k
(x)
− (x − x
i
k
)P
i
0
...i
k
−1
(x)
x
i
k
− x
i
0
Dowód: Własność 1
◦
jest oczywista.
Uwzględniając definicję (
) łatwo zauważyć, że prawa strona własności 2
◦
jest dla
x = x
i
0
równa f
i
0
oraz dla x = x
i
k
równa f
i
k
.
Metody Numeryczne II
32
Dla pozostałych x
i
j
(dla j = 1, . . . , k
− 1) otrzymujemy z prawej strony
(x
i
j
− x
0
)f
i
j
− (x
i
j
− x
i
k
)f
i
j
x
i
k
− x
i
0
= f
i
j
,
co ostatecznie dowodzi własności 2
◦
.
Tabela ilustrująca algorytm Neville’a:
k = 0
1
2
3
x
0
f
0
= P
0
(x)
P
01
(x)
x
1
f
1
= P
1
(x)
P
012
(x)
P
12
(x)
P
0123
(x)
x
2
f
2
= P
2
(x)
P
123
(x)
P
23
(x)
x
3
f
3
= P
3
(x)
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:
T
i+k k
def
= P
i...i+k
(40)
Tablica ilustrująca algorytm:
x
0
f
0
= T
00
(x)
T
11
(x)
x
1
f
1
= T
10
(x)
T
22
(x)
T
21
(x)
T
33
(x)
x
2
f
2
= T
20
(x)
T
32
(x)
T
31
(x)
x
3
f
3
= T
30
(x)
1
◦
T
i0
(x) = f
i
2
◦
T
ij
(x) =
(x
− x
i
−j
)T
i j
−1
(x)
− (x − x
i
)T
i
−1 j−1
(x)
x
i
− x
i
−j
= T
i j
−1
+
T
i j
−1
− T
i
−1 j−1
x
−x
i
−j
x
−x
i
− 1
Metody Numeryczne II
34
Wzór interpolacyjny Newtona
Wielomian interpolacyjny P
∈ Π
n
, P (x
i
) = f
i
dla i = 1, . . . , n można przedstawić
w postaci:
P
0...n
(x) = a
0
+a
1
(x
−x
0
)+a
2
(x
−x
0
)(x
−x
1
)+
· · ·+a
n
(x
−x
0
)
· · · (x−x
n
−1
). (41)
Wartości P (x) można wyznaczać na wzór schematu Hornera (n
− 1 mnożeń).
Współczynniki a
i
można wyznaczyć wprost z:
f
0
= P (x
0
) = a
0
f
1
= P (x
1
) = a
0
+ a
1
(x
1
− x
0
)
f
2
= P (x
2
) = a
0
+ a
1
(x
2
− x
0
) + a
2
(x
2
− x
0
)(x
2
− x
1
)
. . .
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 P
i
0
...i
k
oraz
P
i
0
...i
k
−1
istnieje jednoznacznie współczynnik f
i
0
...i
k
(iloraz różnicowy k-tego
rzędu) taki, że:
P
i
0
...i
k
(x)
− P
i
0
...i
k
−1
(x) = f
i
0
...i
k
(x
− x
i
0
)
· · · (x − x
i
k
−1
).
(42)
Mamy więc następującą postać Newtona częściowego wielomianu P
i
0
...i
k
(x):
f
i
0
+ f
i
0
i
1
(x
−x
0
) + f
i
0
i
1
i
2
(x
−x
0
)(x
−x
1
) +
· · ·+f
i
0
...i
k
(x
−x
0
)
· · · (x−x
i
k
−1
). (43)
Twierdzenie 4 Ilorazy różnicowe f
i
0
...i
k
są niezależne od permutacji
indeksów.
Dowód: f
i
0
...i
k
jest współczynnikiem przy najwyższej potędze wielomianu P
i
0
...i
k
,
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
oraz z (
) mamy:
f
i
0
...i
k
=
P
i
0
...i
k
(x)
− P
i
0
...i
k
−1
(x)
(x
− x
i
0
)
· · · (x − x
i
k
−1
)
2◦
=
(x
−x
i0
)P
i1...ik
(x)
−(x−x
ik
)P
i0...ik−1
(x)
x
ik
−x
i0
− P
i
0
...i
k
−1
(x)
(x
− x
i
0
)
· · · (x − x
i
k
−1
)
=
(x
− x
i
0
)P
i
1
...i
k
(x)
− (x − x
i
k
)P
i
0
...i
k
−1
(x)
− (x
i
k
− x
i
0
)P
i
0
...i
k
−1
(x)
(x
i
k
− x
i
0
)(x
− x
i
0
)
· · · (x − x
i
k
−1
)
=
(x
− x
i
0
)(P
i
1
...i
k
(x)
− P
i
0
...i
k
−1
(x))
(x
i
k
− x
i
0
)(x
− x
i
0
)
· · · (x − x
i
k
−1
)
=
(P
i
1
...i
k
(x)
− P
i
1
...i
k
−1
(x) + P
i
1
...i
k
−1
(x)
− P
i
0
...i
k
−1
(x))
(x
i
k
− x
i
0
)(x
− x
i
1
)
· · · (x − x
i
k
−1
)
=
f
i
1
...i
k
− f
i
0
...i
k
−1
x
i
k
− x
i
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
x
0
f
0
f
01
x
1
f
1
f
012
f
12
f
0123
x
2
f
2
f
123
f
23
x
3
f
3
Wielomian interpolacyjny:
P (x) = f
0
+f
01
(x
−x
0
)+f
012
(x
−x
0
)(x
−x
1
)+
· · ·+f
0...n
(x
−x
0
)
· · · (x−x
n
−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ść i
0
, . . . , i
n
tak, aby
|x − x
i
k
| |x − x
i
k
−1
| dla k = 1, . . . , n.
• Jeżeli wszystkie odcięte x
i
węzłów są uporządkowane tak, że x
i
< x
i+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)}.
x
0
= 0
1
1
x
1
= 1
2
1
4
x
2
= 3
10
P
210
(x) = 10 + 4(x
− 3) + 1(x − 3)(x − 1)
P
210
(x) = (1(x
− 1) + 4)(x − 3) + 10
P
210
(2) = 5
P
120
(x) = 2 + 4(x
− 1) + 1(x − 1)(x − 3)
P
120
(x) = (1(x
− 3) + 4)(x − 1) + 2
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[x
0
, . . . , x
n
, x] zawierającym x i odcięte wszystkich węzłów, spełniająca
f (x)
− P
01...n
(x) =
ω(x)f
(n+1)
(ξ)
(n + 1)!
,
(45)
gdzie
ω(x) = (x
− x
0
)(x
− x
1
)
· · · (x − x
n
).
(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 (x
i
, f
i
) funkcją wymierną:
Φ
µ,ν
(x) =
P
µ,ν
(x)
Q
µ,ν
(x)
=
a
0
+ a
1
x +
· · · + a
µ
x
µ
b
0
+ b
1
x +
· · · + b
ν
x
ν
.
(47)
• Parę (µ, ν) nazywamy stopniem zadania interpolacji wymiernej.
• µ + ν + 2 parametry minus czynnik wspólny
• µ + ν + 1 węzłów (x
i
, f
i
), i = 0, . . . , µ + ν
Rozwiązanie zadania musi być rozwiązaniem układu równań S
µ,ν
:
P
µ,ν
(x
i
)
− f
i
Q
µ,ν
(x
i
) = 0
(48)
czyli
a
0
+ a
1
x
i
+
· · · + a
µ
x
µ
i
− f
i
(b
0
+ b
1
x
i
+
· · · + b
ν
x
ν
i
) = 0
(49)
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ń:
a
0
− b
0
=
0
a
0
+ a
1
− 2(b
0
+ b
1
)
=
0
a
0
+ 2a
1
− 2(b
0
+ 2b
1
)
=
0
Rozwiązanie układu a
0
= b
0
= 0, a
1
= 2, b
1
= 1 nie spełnia warunków zadania –
funkcja wymierna Φ
1,1
(x) =
2x
x
nie jest określona w punkcie x = 0. Zatem
rozwiązanie zadania nie istnieje.
Metody Numeryczne II
42
Definicja 2 Wyrażenia wymierne:
Φ
1
(x) =
P
1
(x)
Q
1
(x)
Φ
2
(x) =
P
2
(x)
Q
2
(x)
Q
1
≡ 0, Q
2
≡ 0
nazywamy równoważnymi gdy przedstawiają tę samą funkcję wymierną tj.
wtedy, gdy
P
1
(x)Q
2
(x) = P
2
(x)Q
1
(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
µ,ν
≡ 0.
Dowód: S
µ,ν
jest jednorodnym układem µ + ν + 1 równań z µ + ν + 2
niewiadomymi, więc ma rozwiązania nietrywialne (a
0
, . . . , a
µ
, b
0
, . . . , b
ν
).
Gdyby Q
µ,ν
≡ 0, czyli ∀
ν
i=0
b
i
= 0, to mielibyśmy
P
µ,ν
(x
i
) = 0,
dla i = 0, . . . , µ + ν + 1.
Ponieważ P
µ,ν
jest stopnia ν < µ + ν + 1, więc mielibyśmy P
µ,ν
≡ 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.
Dowód: Jeśli Φ
1
(x) =
P
1
(x)
Q
1
(x)
oraz Φ
2
(x) =
P
2
(x)
Q
2
(x)
, to
P (x) = P
1
(x)Q
2
(x)
− P
2
(x)Q
1
(x)
(50)
ma µ + ν + 1 różnych zer (x
i
, dla i = 0, . . . , µ + ν + 1). Ponieważ stopień P jest
niewiększy niż µ + ν, to P
≡ 0, czyli Φ
1
i Φ
2
są równoważne.
Wniosek 8 (z twierdzeń
i
) 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
, zadanie A
µ,ν
jest nierozwiązalne.
⇐: Niech ˜Φ
µ,ν
(x) =
P (x)
Q(x)
. Rozważmy dwa przypadki:
1.
∀
i
∈{0,...,µ+ν}
Q(x
i
)
= 0. Wówczas ˜Φ
µ,ν
(x
i
) = f
i
.
2.
∃
i
∈{0,...,µ+ν}
Q(x
i
) = 0. Wówczas węzeł (x
i
, f
i
) byłby nieosiągalny, oraz
mielibyśmy P (x
i
) = 0. Zatem P i Q zawierałyby czynnik x
− x
i
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) A
m,n
można rozpatrywać ciąg wyrażeń wymiernych Φ
µ,ν
stopni (µ, ν), gdzie µ
m, ν n.
• W kolejnym kroku zwiększamy stopień licznika bądź mianownika.
Metody Numeryczne II
47
Odwrotności ilorazów różnicowych
φ(x
i
, x
j
)
def
=
x
i
− x
j
f
i
− f
j
φ(x
i
, x
j
, x
k
)
def
=
x
j
− x
k
φ(x
i
, x
j
)
− φ(x
i
, x
k
)
..
.
φ(x
i
, . . . , x
l
, x
m
, x
n
)
def
=
x
m
− x
n
φ(x
i
, . . . , x
l
, x
m
)
− φ(x
i
, . . . , x
l
, x
n
)
(51)
Uwaga: Nie są symetrycznymi funkcjami swoich argumentów.
Metoda wyznacza wyrażenia wy-
mierne odpowiadające głównej
przekątnej:
µ
0
1
2
3 . . .
ν
0
1
2
3
..
.
···
Metody Numeryczne II
48
Załóżmy, że wyznaczamy funkcję wymierną
Φ
n,n
=
P
n
(x)
Q
n
(x)
(52)
taką, że Φ
n,n
(x
i
) = f
i
dla i = 0, . . . , 2n. Mamy:
P
n
(x)
Q
n
(x)
= f
0
+
P
n
(x)
Q
n
(x)
−
P
n
(x
0
)
Q
n
(x
0
)
= f
0
+
P
n
(x)
−
P
n
(x
0
)
Q
n
(x
0
)
Q
n
(x)
Q
n
(x)
=
f
0
+ (x
− x
0
)
P
n
−1
(x)
Q
n
(x)
= f
0
+
x
− x
0
Q
n
(x)/P
n
−1
(x)
(53)
Dla i = 1, . . . , 2n zachodzi więc:
Q
n
(x
i
)
P
n
−1
(x
i
)
=
x
i
− x
0
f
i
− f
0
= φ(x
0
, x
i
)
(54)
I dalej:
Q
n
(x)
P
n
−1
(x)
= φ(x
0
, x
1
) +
Q
n
(x)
P
n
−1
(x)
−
Q
n
(x
1
)
P
n
−1
(x
1
)
=
φ(x
0
, x
1
) +
x
− x
1
P
n
−1
(x)/Q
n
−1
(x)
(55)
Metody Numeryczne II
49
Kontynuując, otrzymamy wyrażenie reprezentowane przez ułamek łańcuchowy:
Φ
n,n
(x) = f
0
+ x
− x
0
/φ(x
0
, x
1
) + x
− x
1
/φ(x
0
, x
1
, x
2
) +
· · ·
· · · + x − x
2n
−1
/φ(x
0
, x
1
, . . . , x
2n
)
(56)
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)
=
f
0
Φ
1,0
(x)
=
f
0
+ x
− x
0
/φ(x
0
, x
1
)
Φ
1,1
(x)
=
f
0
+ x
− x
0
/φ(x
0
, x
1
) + x
− x
1
/φ(x
0
, x
1
, x
2
)
. . .
(57)
Algorytm: Realizacja wzoru
z pomocą tablicy odwrotności ilorazów:
0
x
0
f
0
1
x
1
f
1
φ(x
0
, x
1
)
2
x
2
f
2
φ(x
0
, x
2
)
φ(x
0
, x
1
, x
2
)
3
x
3
f
3
φ(x
0
, x
3
)
φ(x
0
, x
1
, x
3
)
φ(x
0
, x
1
, x
2
, x
3
)
. . .
Metody Numeryczne II
50
Odwrotne ilorazy różnicowe
ρ(x
i
)
def
=
f
i
ρ(x
i
, x
i+1
)
def
=
x
i
− x
i+1
f
i
− f
i+1
..
.
ρ(x
i
, . . . , x
i+k
)
def
=
x
i
− x
i+k
ρ(x
i
, . . . , x
i+k
−1
)
− ρ(x
i+1
, . . . , x
i+k
)
+ ρ(x
i+1
, . . . , x
i+k
−1
)
(58)
Uwaga: Są symetrycznymi funkcjami swoich argumentów.
Twierdzenie 10 Dla p = 1, 2, . . ., przyjmując ρ(x
0
, . . . , x
p
−2
) = 0 dla p = 1
spełnione jest:
φ(x
0
, . . . , x
p
) = ρ(x
0
, . . . , x
p
)
− ρ(x
0
, . . . , x
p
−2
)
(59)
Dowód: Indukcja. Dla p = 1 teza jest oczywista.
Metody Numeryczne II
51
Zakładamy prawdziwość dla p. Wówczas:
φ(x
0
, . . . , x
p+1
) =
x
p
− x
p+1
φ(x
0
, . . . , x
p
)
− φ(x
0
, . . . , x
p
−1
, x
p+1
)
=
x
p
− x
p+1
ρ(x
0
, . . . , x
p
)
− ρ(x
0
, . . . , x
p
−1
, x
p+1
)
(60)
Z definicji ρ mamy:
ρ(x
p+1
, x
0
, . . . , x
p
)
−ρ(x
0
, . . . , x
p
−1
) =
x
p+1
− x
p
ρ(x
p+1
, x
0
, . . . , x
p
−1
)
− ρ(x
0
, . . . , x
p
)
(61)
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) = f
0
+ x
− x
0
/ρ(x
0
, x
1
) + x
− x
1
/ρ(x
1
, x
1
, x
2
)
− ρ(x
1
) +
· · ·
· · · + x − x
2n
−1
/ρ(x
1
, . . . , x
2n
)
− ρ(x
1
, . . . , x
2n
−2
)
(62)
Metody Numeryczne II
52
Uogólnienie algorytmu Neville’a
Oznaczenia:
• Φ
µ,ν
s
(x)
def
=
P
µ,ν
s
(x)
Q
µ,ν
s
(x)
– wyrażenie wymierne spełniające warunki: Φ
µ,ν
s
(x
i
) = f
i
dla i = s, s + 1, . . . , s + µ + ν. P
µ,ν
s
(x) i Q
µ,ν
s
(x) są wielomianami stopni
odpowiednio µ i ν.
• α
i
def
= x
− x
i
Twierdzenie 11 Dla µ
1 oraz ν 1 zachodzą:
(a) Φ
µ,ν
s
(x) = Φ
µ
−1,ν
s+1
(x) +
Φ
µ
−1,ν
s+1
(x)
− Φ
µ
−1,ν
s
(x)
α
s
α
s+µ+ν
1
−
Φ
µ
−1,ν
s+1
(x)
− Φ
µ
−1,ν
s
(x)
Φ
µ
−1,ν
s+1
(x)
− Φ
µ
−1,ν−1
s+1
(x)
− 1
(b) Φ
µ,ν
s
(x) = Φ
µ,ν
−1
s+1
(x) +
Φ
µ,ν
−1
s+1
(x)
− Φ
µ,ν
−1
s
(x)
α
s
α
s+µ+ν
1
−
Φ
µ,ν
−1
s+1
(x)
− Φ
µ,ν
−1
s
(x)
Φ
µ,ν
−1
s+1
(x)
− Φ
µ
−1,ν−1
s+1
(x)
− 1
Metody Numeryczne II
53
• Wzorów z twierdzenia
można użyć do
wyznaczenia wartości wyrażeń wymiernych
przy odpowiednio wzrastających stopniach.
• Różne ścieżki to różne metody i różne final-
ne funkcje wymierne.
• Dla ν = 0 i wzrastającego µ mamy interpo-
lację wielomianową Neville’a.
• µ = 0 i wzrastające ν odpowiadają interpo-
lacji wielomianowej dla węzłów (x
1
, 1/f
i
).
µ
0
1
2
3
. . .
ν
0
1
2
3
..
.
···
• 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ć:
T
ik
ozn
= Φ
µ,ν
s
,
gdzie i = s + µ + ν, k = µ + ν
(63)
Mamy wówczas następujący algorytm:
T
i0
= f
i
,
T
i,
−1
= 0
T
ik
= T
i,k
−1
+
T
i,k
−1
− T
i
−1,k−1
x
− x
i
−k
x
− x
i
1
−
T
i,k
−1
− T
i
−1,k−1
T
i,k
−1
− T
i
−1,k−2
− 1
(64)
Tabela ilustrująca algorytm:
(µ, ν) =
(0, 0)
(0, 1)
(1, 1)
(1, 2)
T
00
= f
0
T
0,
−1
= 0
T
11
T
10
= f
1
T
22
T
1,
−1
= 0
T
21
T
33
T
20
= f
2
T
32
T
2,
−1
= 0
T
31
T
30
= f
3
Metody Numeryczne II
55
Interpolacja wielomianowa a wymierna – ctg(2.5
◦
) = 22.9037655484
x
1
f
i
= ctg(x
i
)
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
Definicja 3 Niech ∆
def
=
{x
0
, x
1
, . . . , x
n
}, gdzie a = x
0
< x
1
< . . . < x
n
= 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
∆
∈ C
2
[a, b] (dwukrotnie ciągle różniczkowalna na [a, b])
• S
∆
na każdym podprzedziale [x
i
, x
i+1
], i = 0, . . . , n
− 1 pokrywa się z
wielomianem trzeciego stopnia.
Metody Numeryczne II
57
Oznaczenie: Dla zbioru liczb rzeczywistych Y =
{y
0
, . . . , y
n
} przez
S
∆
(Y,
·)
(65)
rozumieć będziemy funkcję sklejaną S
∆
spełniającą warunki S
∆
(Y, x
i
) = y
i
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:
(a)
S
∆
(Y, a) = S
∆
(Y, b) = 0.
(b)
S
(k)
∆
(Y, a) = S
(k)
∆
(Y, b) dla k = 0, 1, 2. S
∆
(Y,
·) jest okresowa.
(c)
S
∆
(Y, a) = y
0
, S
∆
(Y, b) = y
n
dla danych y
0
, y
n
∈ R.
(66)
Metody Numeryczne II
58
Interpolacja funkcjami sklejanymi
Niech ∆ =
{x
0
, x
1
, . . . , x
n
} będzie ustalonym podziałem przedziału [a, b] oraz
Y =
{y
0
, y
1
, . . . , y
n
} ⊂ R
Oznaczenie:
h
j+1
ozn
= x
j+1
− x
j
,
j = 0, 1, . . . , n
− 1
(67)
Definicja 4 Momentami funkcji S
∆
(Y,
·) nazywamy wartości jej drugich
pochodnych w węzłach:
M
j
def
= S
∆
(Y, x
j
),
j = 0, 1, . . . , n
(68)
Uwagi:
• Druga pochodna funkcji sklejanej to funkcja kawałkami liniowa (liniowa na
przedziałach [x
i
, x
i+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:
S
∆
(Y, x) = M
j
x
j+1
− x
h
j+1
+ M
j+1
x
− x
j
h
j+1
dla x
∈ [x
j
, x
j+1
]
(69)
Całkując otrzymujemy
S
∆
(Y, x)
=
−M
j
(x
j+1
− x)
2
2h
j+1
+ M
j+1
(x
− x
j
)
2
2h
j+1
+ A
j
(70)
S
∆
(Y, x)
=
M
j
(x
j+1
− x)
3
6h
j+1
+ M
j+1
(x
− x
j
)
3
6h
j+1
+ A
j
(x
− x
j
) + B
j
(71)
dla x
∈ [x
j
, x
j+1
], j = 0, 1, . . . , n
− 1, gdzie A
j
oraz B
j
oznaczają stałe całkowania.
Ze spełniania odpowiednich punktów węzłowych otrzymujemy:
M
j
h
2
j+1
6
+ B
j
=
y
j
(72)
M
j+1
h
2
j+1
6
+ A
j
h
j+1
+ B
j
=
y
j+1
(73)
Metody Numeryczne II
60
Stąd:
B
j
=
y
j
− M
j
h
2
j+1
6
(74)
A
j
=
y
j+1
− M
j+1
h
2
j+1
6
− B
j
h
j+1
=
y
j+1
− y
j
h
j+1
−
h
j+1
6
(M
j+1
− M
j
)
(75)
Mamy więc dla x
∈ [x
j
, x
j
+ 1]:
S
∆
(Y, x) = α
j
+ β
j
(x
− x
j
) + γ
j
(x
− x
j
)
2
+ δ
j
(x
− x
j
)
3
(76)
gdzie:
α
j
= S
∆
(Y, x
j
) = y
j
,
γ
j
=
S
∆
(Y, x
j
)
2
=
M
j
2
(77)
β
j
= S
∆
(Y, x
j
) =
−
M
j
h
j+1
2
+ A
j
=
y
j+1
− y
j
h
j+1
−
h
j+1
6
(M
j+1
+ 2M
j
)
(78)
δ
j
=
S
∆
(Y, x
+
j
)
6
=
M
j+1
− M
j
6h
j+1
(79)
Metody Numeryczne II
61
Wyznaczanie momentów
Założenie ciągłości S
∆
(Y,
·) w węzłach wyznacza n − 1 równań dla momentów.
Podstawiając A
j
z (
) otrzymujemy:
S
∆
(Y, x) =
−M
j
(x
j+1
− x)
2
2h
j+1
+ M
j+1
(x
− x
j
)
2
2h
j+1
+
y
j+1
− y
j
h
j+1
−
h
j+1
6
(M
j+1
− M
j
)
(80)
W szczególności dla j = 1, 2, . . . , n
− 1 mamy:
S
∆
(Y, x
−
j
)
=
y
j
− y
j
−1
h
j
+
h
j
3
M
j
−
h
j
6
M
j
−1
(81)
S
∆
(Y, x
+
j
)
=
y
j+1
− y
j
h
j+1
−
h
j+1
3
M
j
−
h
j+1
6
M
j+1
(82)
A z równości granic:
h
j
6
M
j
−1
+
h
j
+ h
j+1
3
M
j
+
h
j+1
6
M
j+1
=
y
j+1
− y
j
h
j+1
−
y
j
− y
j
−1
h
j
(83)
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 (
Przypadek (
S
∆
(Y, a) = M
0
= 0
S
∆
(Y, b) = M
n
= 0
(84)
Przypadek (
.b): funkcja okresowa, czyli h
n+1
= h
1
, M
n+1
= M
1
, y
n+1
= y
1
:
S
∆
(Y, a) = S
∆
(Y, b)
≡ M
0
= M
n
(85)
S
∆
(Y, a) = S
∆
(Y, b)
≡
h
n
6
M
n
−1
+
h
n
+ h
1
3
M
n
+
h
1
6
M
1
=
y
1
− y
n
h
1
−
y
n
− y
n
−1
h
n
(86)
Przypadek (
S
∆
(Y, a) = y
0
≡
h
1
3
M
0
+
h
1
6
M
1
=
y
1
− y
0
h
1
− y
0
(87)
Metody Numeryczne II
63
S
∆
(Y, b) = y
n
≡
h
n
6
M
n
−1
+
h
n
3
M
n
= y
n
−
y
n
− y
n
−1
h
n
(88)
Ogólna postać dla równań (
µ
j
M
j
−1
+ 2M
j
+ λ
j
M
j+1
= d
j
,
j = 1, 2, . . . , n
− 1,
(89)
gdzie:
λ
j
=
h
j+1
h
j
+ h
j+1
,
µ
j
= 1
− λ
j
=
h
j
h
j
+ h
j+1
,
d
j
=
6
h
j
+ h
j+1
y
j+1
− y
j
h
j+1
−
y
j
− y
j
−1
h
j
(90)
Dodatkowo:
Przypadek (
λ
0
= 0,
d
0
= 0,
µ
n
= 0,
d
n
= 0
(91)
Metody Numeryczne II
64
Przypadek (
λ
0
= 1,
d
0
=
6
h
1
y
1
− y
0
h
1
− y
0
µ
n
= 1,
d
n
=
6
h
n
y
n
−
y
n
− y
n
−1
h
n
(92)
Zatem dla przypadków (
.c) otrzymujemy układ równań:
2
λ
0
0
µ
1
2
λ
1
. ..
µ
n
−1
2
λ
n
−1
µ
n
2
M
0
M
1
..
.
M
n
−1
M
n
=
d
0
d
1
..
.
d
n
−1
d
n
(93)
Metody Numeryczne II
65
Przypadek (
λ
n
=
h
1
h
n
+ h
1
,
µ
n
= 1
− λ
n
=
h
n
h
n
+ h
1
d
n
=
6
h
n
+ h
1
y
1
− y
n
h
1
−
y
n
− y
n
−1
h
n
(94)
Otrzymujemy następujący układ równań:
2
λ
1
µ
1
µ
2
2
λ
2
. ..
µ
n
−1
2
λ
n
−1
λ
n
µ
n
2
M
1
M
2
..
.
M
n
−1
M
n
=
d
1
d
2
..
.
d
n
−1
d
n
(95)
oraz M
0
= M
n
Metody Numeryczne II
66
Interpolacja funkcjami sklejanymi - podsumowanie
Twierdzenie 12 Układy równań liniowych
i
są nieosobliwe dla
dowolnego podziału ∆ przedziału [a, b].
Schemat dowodu: Najpierw dowodzimy własności:
Az = w
⇒ max
i
|z
i
| max
i
|w
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
a
f (x)dx = F (b)
− F (a),
gdzie F
(x) = f (x)
(96)
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]:
x
i
= a + ih,
i = 0, . . . , n,
gdzie h =
b
− a
n
, n > 0, n
∈ N
(97)
Niech P
n
będzie wielomianem interpolacyjnym stopnia
n takim, że:
P
n
(x
i
) = f
i
= f (x
i
)
dla i = 0, . . . , n
(98)
Ze wzoru interpolacyjnego Lagrange’a mamy:
P
n
(x) =
n
i=0
f
i
L
i
(x),
gdzie L
i
(x) =
n
k=0
k
=i
x
− x
k
x
i
− x
k
(99)
Metody Numeryczne II
69
• Wprowadzamy zmienną t taką, że x = a + ht. Mamy więc:
L
i
(x) = ϕ
i
(t) =
n
k=0
k
=i
t
− k
i
− k
(100)
• Całkujemy:
b
a
P
n
(x)dx =
n
i=0
f
i
b
a
L
i
(x)dx = h
n
i=0
f
i
n
0
ϕ
i
(t)dt = h
n
i=0
f
i
α
i
(101)
gdzie wagi
α
i
=
n
0
ϕ
i
(t)dt
(102)
nie zależą od całkowanej funkcji ani od przedziału całkowania.
Metody Numeryczne II
70
• Dla n = 2 otrzymujemy tzw. kwadraturę Simpsona:
α
0
=
2
0
(t
− 1)(t − 2)
(0
− 1)(0 − 2)
dt =
1
2
2
0
(t
2
− 3t + 2)dt =
=
1
2
1
3
t
3
−
3
2
t
2
+ 2t
2
0
=
1
2
8
3
− 6 + 4
=
1
3
α
1
=
2
0
(t
− 0)(t − 2)
(1
− 0)(1 − 2)
dt =
−
2
0
(t
2
− 2t)dt =
4
3
α
2
=
2
0
(t
− 0)(t − 1)
(2
− 0)(2 − 1)
dt =
1
2
2
0
(t
2
− t)dt =
1
3
(103)
Ostatecznie:
b
a
f (x)dx
≈
b
a
P
n
(x)dx =
1
3
(f
0
+ 4f
1
+ f
2
)
(104)
Metody Numeryczne II
71
• w ogólnej postaci: kwadratury Newtona-Cotesa
b
a
P
n
(x)dx = h
n
i=0
f
i
α
i
,
f
i
= f (a + ih),
h =
b
− a
n
(105)
– Zachodzi:
n
i=0
α
i
= n
– Niech s będzie wspólnym mianownikiem wag α
i
. Wówczas σ
i
ozn
= sα
i
są
całkowite oraz:
b
a
P
n
(x)dx = h
n
i=0
f
i
α
i
=
b
− a
ns
n
i=0
σ
i
f
i
(106)
• reszta kwadratury tj. błąd aproksymacji:
b
a
P
n
(x)dx
−
b
a
f (x)dx = h
p+1
Kf
(p)
(ξ),
gdzie ξ
∈ (a, b)
(107)
Metody Numeryczne II
72
Zestawienie parametrów kwadratur Newtona-Cotesa:
n
σ
i
ns
Reszta
Kwadratura
1
1
1
2
h
3 1
12
f
(2)
(ξ)
Trapezów
2
1
4
1
6
h
5 1
90
f
(4)
(ξ)
Simpsona
3
1
3
3
1
8
h
5 3
80
f
(4)
(ξ)
3/8
4
7
32
12
32
7
90
h
7 8
945
f
(6)
(ξ)
Milne’a
5
19
75
50
50
75
19
288
h
7 275
12096
f
(6)
(ξ)
—
6
41
216
27
272
27
216
41
840
h
9
9
1400
f
(8)
(ξ)
Weddle’a
Dla n > 6 pojawiają się ujemne współczynniki co stanowi niebezpieczeństwo
dużych błędów numerycznych.
Metody Numeryczne II
73
Kwadratury złożone
• Całkujemy w [x
i
, x
i+1
] dla i = 0, . . . , n
− 1, gdzie x
i
= a + ih, h =
b
−a
n
• Złożona kwadratura trapezów
T (h)
def
=
n
−1
i=0
1
2
h [f (x
i
) + f (x
i+1
)] = h
f (a) + f (b)
2
+
n
−2
i=1
f (a + ih)
(108)
Dla f
∈ C
2
[a, b] reszta kwadratury w przedziale [x
i
, x
i+1
] wynosi
1
12
h
3
f
(2)
(ξ
i
),
gdzie ξ
i
∈ [x
i
, x
i+1
]
(109)
Sumując reszty otrzymujemy:
T (h)
−
b
a
f (x)dx =
h
3
12
n
−1
i=0
f
(2)
(ξ
i
) =
=
b
− a
12
h
2
1
n
n
−1
i=0
f
(2)
(ξ
i
) =
b
− a
12
h
2
f
(2)
(ξ),
gdzie ξ
∈ (a, b)
(110)
Metody Numeryczne II
74
• Złożona kwadratura Simpsona (n parzyste)
S(h)
def
=
1
3
h[f (a) + 4f (a + h) + 2f (a + 2h) + 4f (a + 3h) + . . .
. . . + 2f (b
− 2h) + 4f(b − h) + f(b)]
(111)
Reszta kwadratury:
S(h)
−
b
a
f (x)dx =
h
5
12
n/2
−1
i=0
f
(4)
(ξ
i
) =
=
h
4
90
b
− a
2
2
n
n/2
−1
i=0
f
(4)
(ξ
i
) =
b
− a
180
h
4
f
(4)
(ξ),
gdzie ξ
∈ (a, b)
(112)
Uwagi:
• reszty zbiegają do zera z potęgą h
• rząd kwadratury
def
= 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
75
Kwadratury z interpolacji Hermite’a
Interpolacja Hermite’a
• zagadnienie interpolacji węzłów (x
i
, y
(k)
i
), k = 0, . . . , n
i
− 1, i − 0, . . . , m
wielomianem P stopnia n =
m
i=0
n
i
− 1 spełniającego P
(k)
(x
i
) = y
(k)
i
.
• rozwiązanie jest jednoznaczne
• uogólnienie wielomianów Lagrange’a:
P (x) =
m
i=0
n
i
−1
k=0
y
(k)
i
L
ik
(x),
gdzie L
(σ)
ik
(x) =
1
gdy i = j
∧ k = σ
0
w przec. przyp.
(113)
Metody Numeryczne II
76
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)[t
2
−2t
2
(t
−1)]+f
(0)t(t
−1)
2
+f
(1)t
2
(t
−1)
(114)
Po scałkowaniu:
1
0
P (t)dt =
1
2
(f (0) + f (1)) +
1
12
(f
(0)
− f
(1))
(115)
Po przekształceniu zmiennych:
b
a
f (x)dx
≈ M(h) =
1
2
h(f (a) + f (b)) +
1
12
h
2
(f
(a)
− f
(b))
(116)
gdzie h = b
− a.
Reszta kwadratury:
M (h)
−
b
a
f (x)dx =
−
1
720
h
5
f
(4)
(ξ),
ξ
∈ (a, b)
(117)
Metody Numeryczne II
77
Wzór sumacyjny Eulera-Maclaurina
Twierdzenie 13 Niech g
∈ C
2m+2
[0, 1], m
∈ N. Wówczas:
1
0
g(t)dt =
1
2
(g(0) + g(1)) +
m
l=1
B
2l
(2l)!
(g
(2l
−1)
(0)
− g
(2l
−1)
(1))
−
B
2m+2
(2m + 2)!
g
(2m+2)
(ξ),
ξ
∈ (0, 1)
(118)
gdzie B
k
to liczby Bernoulliego (B
2
=
1
6
, B
4
=
−
1
30
, B
6
=
1
42
, B
8
=
−
1
30
, . . . ).
Schemat dowodu
• Wielomiany i liczby Bernoulliego – definicja i własności.
• Stosowne rozpisanie całki
1
0
g(t)dt.
Metody Numeryczne II
78
Wielomiany i liczby Bernoulliego
Definicja 5 Wielomianami Bernoulliego nazywamy rodzinę wielomianów
B
i
, i = 0, 1, . . . spełniającą następujące warunki:
1. B
1
(x) = x +
1
2
2. B
k+1
= (k + 1)B
k
(x),
k = 0, 1, . . .
3. B
2l+1
(0) = B
2l+1
(1) = 0,
l = 1, 2, . . .
Przyklady:
B
0
(x) = 1,
B
1
(x) = x +
1
2
,
B
2
(x) = x
2
− x +
1
6
,
B
3
(x) = x
3
−
3
2
x
2
+
1
2
x,
B
4
(x) = x
4
− x
3
+ x
2
−
1
30
(119)
Definicja 6 Liczby Bernoulliego to wartości B
k
= B
k
(0), k = 0, 1, . . ..
Wlasnosci:
• Wielomian B
k
jest stopnia k.
• (−1)
k
B
k
(1
− x) = B
k
(x)
• B
k
(0) = B
k
(1) = B
k
dla k
= 1
Metody Numeryczne II
79
Wykorzystanie wielomianów Bernoulliego dla dowodu wzoru (
Całkowanie przez części:
1
0
g(t)dt
=
B
1
(t)g(t)
|
1
0
−
1
0
B
1
(t)g
(t)dt
1
0
B
1
(t)g
(t)dt
=
1
2
B
2
(t)g
(t)
1
0
−
1
2
1
0
B
2
(t)g
(t)dt
..
.
(120)
1
0
B
k
−1
(t)g
(k
−1)
(t)dt
=
1
k
B
k
(t)g
(k
−1)
(t)
1
0
−
1
k
1
0
B
k
(t)g
(k)
(t)dt
Metody Numeryczne II
80
Wzór sumacyjny Eulera-Maclaurina
Rozszerzając wzór (
) na n
∈ N przedziałów otrzymujemy dla g ∈ C
2m+2
[0, n]:
n
0
g(t)dt =
g(0) + g(n)
2
+
n
−1
i=1
g(i) +
m
l=1
B
2l
(2l)!
(g
(2l
−1)
(0)
− g
(2l
−1)
(n))
−
B
2m+2
(2m + 2)!
ng
(2m+2)
(ξ),
ξ
∈ (0, n) (121)
W ogólnym przypadku przedziału [a, b] i jego równomiernego podziału na n części
(h =
b
−a
n
):
b
a
f (x)dx = T (h) +
m
l=1
h
2l
B
2l
(2l)!
(f
(2l
−1)
(a)
− f
(2l
−1)
(b))
−h
2m+2
B
2m+2
(2m + 2)!
(b
− a)f
(2m+2)
(ξ),
ξ
∈ (a, b),
(122)
gdzie T (h) oznacza złożoną kwadraturę trapezów.
Metody Numeryczne II
81
Kwadratury ekstrapolacyjne
Rozwinięcie (
) złożonego wzoru trapezów T (h) dla funkcji f w zależności od
długości kroku h:
T (h) = τ
0
+ τ
1
h
2
+
· · · + τ
m
h
2m
+ α
m+1
(h)h
2m+2
,
(123)
gdzie
τ
0
=
b
a
f (t)dt,
τ
k
=
B
2k
(2k)!
(f
(2k
−1)
(b)
− f
(2k
−1)
(a)), k = 1, . . . , m
(124)
oraz
α
m+1
(h) =
B
2m+2
(2m + 2)!
(b
− a)f
(2m+2)
(ξ(h)),
ξ(h)
∈ (a, b)
(125)
Metody Numeryczne II
82
Kwadratura Romberga
Jeśli f ma wszystkie pochodne w [a, b], to w przejściu do granicy m
→ ∞
otrzymujemy szereg potęgowy:
τ
0
+ τ
1
h
2
+ τ
2
h
4
+
· · ·
(126)
Reszta rozwinięcia maleje z malejącym h. Rozwinięcie można traktować jak
wielomian zmiennej h
2
przyjmujący wartość całki τ
0
w h = 0.
Dla ciągu długości kroków
{h
i
=
b
−a
n
i
: i = 0, 1, . . .
}, gdzie n
0
= 1 oraz
{n
i
: i = 0, 1, . . .
} jest ciągiem ściśle rosnącym, wyznaczamy kwadratury trapezów:
T
i0
= T (h
i
),
i = 0, 1, . . . , m,
(127)
które traktujemy jako węzły zadania interpolacji.
Jeśli ˜
T
mm
(h) jest wielomianem zmiennej h
2
interpolującym te węzły, to wartość
ekstrapolowana ˜
T
mm
(0) jest przybliżeniem szukanej całki τ
0
.
Metody Numeryczne II
83
Interpolacja Neville’a dla kwadratury Romberga
Cel: Wyznaczyć ˜
T
mm
(0)
Algorytm: Dla indeksów i, k takich, że 1
k i m definiujemy ˜
T
ik
(h) jako
wielomian zmiennej h
2
stopnia co najwyżej k spełniający:
˜
T
ik
(h
j
) = T (h
j
)
dla j = i
− k, . . . , i.
(128)
Wartość ˜
T
mm
(0) wyznaczamy algorytmem Neville’a:
T
ik
= T
i,k
−1
+
T
i,k
−1
− T
i
−1,k−1
h
i
−k
h
i
2
− 1
.
(129)
Wówczas:
T
ik
= ˜
T
ik
(0)
(130)
Metody Numeryczne II
84
Tablica algorytmu Neville’a dla kwadratury Romberga:
h
2
0
T (h
0
) = T
00
T
11
h
2
1
T (h
1
) = T
10
T
22
T
21
T
33
h
2
2
T (h
2
) = T
20
T
32
T
31
h
2
3
T (h
3
) = T
30
Uwagi:
• Niektóre kwadratury dla i = k to kwadratury Newtona-Cotesa:
– h
1
=
h
0
2
⇒ T
11
jest kwadraturą Simpsona
– dodatkowo h
2
=
h
1
2
⇒ T
22
jest kwadraturą Milne’a
• oryginalna metoda Romberga: h
i
=
h
i
−1
2
, i = 1, 2 . . .
• ciąg Bulirsch’a: h
2i
−1
=
h
2i
−2
2
, h
2i
=
h
2i
−2
3
, i = 1, 2 . . .
Metody Numeryczne II
85
Interpolacja wymierna dla kwadratury Romberga
Dla interpolacji wielomianowej
• trudno określić dokładność,
• często warunek stopu to: |T
m
−1,m−1
− T
m,m
−1
| εs, gdzie s jest
oszacowaniem wartości całki
b
a
|f(t)|dt,
W interpolacji wymiernej mamy:
T
i0
=
T (h
i
)
(131)
T
ik
=
T
i,k
−1
+
T
i,k
−1
− T
i
−1,k−1
h
i
−k
h
i
2
1
−
T
i,k
−1
−T
i
−1,k−1
T
i,k
−1
−T
i
−1,k−2
− 1
,
1
k i m (132)
• można oszacować reszty
Metody Numeryczne II
86
Kwadratury Gaussa
Zadanie: Wyznaczanie całek postaci
I(f ) =
b
a
ω(x)f (x)dx,
(133)
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]
• dla k = 0, 1, . . . momenty µ
k
=
b
a
x
k
ω(x)dx istnieją i są skończone
• dla wielomianów s(x) nieujemnych na [a, b] z równości
b
a
ω(x)s(x)dx = 0
wynika, że s
≡ 0.
Metody Numeryczne II
87
Rozpatrywane kwadratury:
˜
I(f ) =
n
i=0
w
i
f (x
i
).
(134)
• W przeciwieństwie do kwadratur Newtona-Cotesa tutaj węzły nie muszą być
równoodległe.
• Jak dobrać odcięte x
i
oraz w
i
, żeby zmaksymalizować rząd kwadratury?
Metody Numeryczne II
88
Podstawowe definicje:
Definicja 7 Zbiorem unormowanych wielomianów rzeczywistych
stopnia j nazywamy
Π
j
def
=
{p : p(x) = x
j
+ a
1
x
j
−1
+
· · · + a
j
}
(135)
Definicja 8 Iloczynem skalarnym funkcji f i g nazywamy
(f, g)
def
=
b
a
ω(x)f (x)g(x)dx.
(136)
Definicja 9 Przez L
2
[a, b] oznaczać będziemy przestrzeń wszystkich funkcji
określonych na przedziale [a, b], dla których całka
(f, f ) =
b
a
ω(x)(f (x))
2
dx
(137)
jest dobrze określona i skończona.
Metody Numeryczne II
89
Twierdzenie 14 (Wielomiany ortogonalne) Istnieją wielomiany
p
j
∈ Π, j = 0, 1, . . . spełniające
(p
i
, p
k
) = 0
dla i
= k.
(138)
Wielomiany te są wyznaczone jednoznacznie wzorami:
p
0
(x)
≡ 1
(139)
p
i+1
(x)
≡ (x − δ
i+1
)p
i
(x)
− γ
2
i+1
p
i
−1
(x),
(140)
gdzie p
−1
(x) = 0 oraz
δ
i+1
=
(xp
i
, p
i
)
(p
i
, p
i
)
dla i
0
(141)
γ
2
i+1
=
0
dla i = 0
(p
i
,p
i
)
(p
i
−1
,p
i
−1
)
dla i
1
(142)
Dowód: Ortogonalizacja Grama-Schmidta:
Z definicji Π
0
mamy p
0
(x)
≡ 1.
Załóżmy prawdziwość tezy dla wszystkich j
i. Dowolny wielomian p
i+1
∈ Π
i+1
Metody Numeryczne II
90
ma jednoznaczne przedstawienie
p
i+1
(x) = (x
− δ
i+1
)p
i
(x) + c
i
−1
p
i
−1
(x) +
· · · + c
0
p
0
(x).
(143)
Jeśli p
i+1
∈ Π
i+1
, to
• dla j = i mamy:
(p
i+1
, p
j
) = (xp
i
, p
i
)
− δ
i+1
(p
i
, p
i
) = 0
(144)
• dla j < i mamy:
(p
i+1
, p
j
) = (xp
j
, p
i
) + c
j
(p
j
, p
j
) = 0
(145)
Z założeń dla ω(x) przy podstawieniu s = p
2
j
wynika, że (p
j
, p
j
)
= 0, więc c
j
jest wyznaczone jednoznacznie.
Z założenia indukcyjnego dla j < i mamy:
p
j+1
(x) = (x
− δ
j+1
)p
j
(x)
− γ
2
j+1
p
j
−1
(x)
(146)
Metody Numeryczne II
91
Stąd (p
j+1
, p
i
) = (xp
j
, p
i
), a zatem
c
j
=
−
(p
j+1
, p
i
)
(p
j
, p
j
)
=
−γ
2
j+1
dla j = i
− 1
0
dla j < i
− 1
,
(147)
czyli
p
i+1
(x) = (x
− δ
i+1
)p
i
(x)
− γ
2
i
−1
p
i
−1
(x).
(148)
Wniosek 15 Jeśli p
∈ Π
n
−1
, to (p, p
n
) = 0.
Twierdzenie 16 Wielomian p
n
ma n + 1 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 t
i
, i = 1, . . . , n, macierz
A =
p
0
(t
1
)
· · · p
n
−1
(t
1
)
..
.
..
.
..
.
p
0
(t
n
)
· · · p
n
−1
(t
n
)
(149)
jest nieosobliwa.
Metody Numeryczne II
92
Twierdzenie 18 Niech x
1
, . . . , x
n
będą zerami wielomianu p
n
oraz niech
w
1
, . . . , w
n
będą rozwiązaniem nieosobliwego układu równań
n
i=0
p
k
(x
i
)w
i
=
(p
0
, p
0
)
dla k = 0,
0
dla k = 1, . . . , n
− 1.
(150)
Wówczas w
i
> 0 dla i = 1, . . . , n, oraz dla każdego wielomianu p
∈ Π
2n
−1
b
a
ω(x)p(x)dx =
n
i=1
w
i
p(x
i
).
(151)
Jeśli w
i
, x
i
, i = 1, . . . , n spełniają (
) dla wszystkich p
∈ Π
2n
−1
, to x
i
są
zerami wielomianu p
n
, a w
i
spełniają (
Nie istnieją takie x
i
, w
i
, i = 1, . . . , n, które spełniają (
) dla wszystkich
p
∈ Π
2n
.
Metody Numeryczne II
93
Twierdzenie 19 Zera wielomianu ortogonalnego p
n
są wartościami własnymi
macierzy trójdiagonalnej
J
n
=
δ
1
γ
2
0
γ
2
δ
2
γ
3
..
.
. .. ...
0
γ
n
δ
n
.
(152)
Twierdzenie 20 (Reszta kwadratury Gaussa) Jeżeli f
∈ C
2n
[a, b], to dla
pewnego ξ
∈ (a, b)
b
a
ω(x)f (x)dx
−
n
i=1
w
i
f (x
i
) =
f
(2n)
(ξ)
(2n)!
(p
n
, p
n
).
(153)
Metody Numeryczne II
94
Najprostsza kwadratura Gaussa:
Dla funkcji wagowej ω
≡ 1 i przedziału [a, b] = [−1, 1], wielomianami ortogonalymi
są:
p
k
(x) =
k!
(2k)!
d
k
dx
k
(x
2
− 1)
k
,
k = 0, 1, . . .
(154)
Inne kwadratury Gaussa:
[a, b]
ω(x)
Wielomiany ortogonalne
[
−1, 1]
(1
− x
2
)
−
1
2
T
n
(x) – Czebyszewa
[0,
∞]
e
−x
L
n
(x) – Laguerre’a
[
−∞, ∞]
e
−x
2
H
n
(x) – Hermite’a
Metody Numeryczne II
95
Metody aproksymacji
Zadanie aproksymacji funkcji f (x) polega na wyznaczeniu parametrów
a
0
, . . . , a
m
pewnej funkcji F tak, by zminimalizować
||f(x) − F (a
0
, . . . , a
m
; 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) = a
0
ϕ
0
(x) + a
1
ϕ
1
(x) +
· · · + a
m
ϕ
m
(x)
(155)
• aproksymacja wymierna
F (x) =
a
0
ϕ
0
(x) + a
1
ϕ
1
(x) +
· · · + a
m
ϕ
m
(x)
b
0
ψ
0
(x) + b
1
ψ
1
(x) +
· · · + b
m
ψ
m
(x)
(156)
Metody Numeryczne II
96
Wybór przestrzeni funkcji i minimalizowanej normy
• Interpolacja - przestrzeń gwarantuje zerową wartość normy
• Aproksymacja średniokwadratowa:
||f||
2,ω
=
b
a
ω(x)f (x)
2
dx
1
2
(157)
w przypadku dyskretnym
||f||
2,ω
=
n
i=0
ω(x
i
)f (x
i
)
2
(158)
• Aproksymacja jednostajna
||f|| = sup
x
∈[a,b]
|f(x)|
(159)
Metody Numeryczne II
97
Aproksymacja średniokwadratowa wielomianami uogólnionymi
Rozpartujemy funkcje postaci
F (x) =
m
i=0
a
i
ϕ
i
(x).
(160)
ϕ
0
, . . . , ϕ
m
to funkcje bazowe m + 1 wymiarowej przestrzeni liniowej.
Różne układy bazowe:
• 1, x, x
2
, . . . , x
n
• 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
98
Zadanie: minimalizacja błędu:
E(a
0
, . . . , a
m
) =
n
i=0
ω(x
i
)
f(x
i
)
−
m
j=0
a
j
ϕ
j
(x
i
)
2
(161)
Minimum E w punkcie, gdzie pochodne cząstkowe są zerowe
δE
δa
k
=
−
n
i=0
2ω(x
i
)
f(x
i
)
−
m
j=0
a
j
ϕ
j
(x
i
)
ϕ
k
(x
i
)
k = 0, . . . , m
(162)
Układ normalny w postaci macierzowej:
D
T
DA = D
T
f
(163)
gdzie
D =
ϕ
0
(x
0
)
· · · ϕ
m
(x
0
)
..
.
..
.
ϕ
0
(x
n
)
· · · ϕ
m
(x
n
)
, A =
a
0
..
.
a
m
, f =
f (x
0
)
..
.
f (x
n
)
(164)
Metody Numeryczne II
99
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)| < ε.
(165)
Funkcje bazowe:
ϕ
i
(x) = x
i
.
(166)
Warunki minimum:
δE
δa
k
= 0
≡
n
i=0
f(x
i
)
−
m
j=0
a
j
x
j
i
x
k
i
= 0
(167)
czyli
n
i=0
f (x
i
)x
k
i
=
m
j=0
a
j
n
i=0
x
j+k
i
,
k = 0, . . . , m
(168)
Metody Numeryczne II
100
W prostszej formie:
m
i=0
a
i
g
ik
= ρ
k
,
k = 0, . . . , m
(169)
g
ij
=
n
j=0
x
i+k
j
,
ρ
k
=
n
j=0
f (x
j
)x
k
j
(170)
Uwagi:
• Dla m n rozwiązanie jest jednoznaczne
• Dla m = n rozwiązaniem jest wielomian interpolacyjny
• Dla m = 1 mamy aproksymację funkcją liniową y = a
0
+ a
1
x:
a
1
=
EXY
− EXEY
EX
2
− (EX)
2
,
a
0
= EY
− a
1
EX
(171)
• Zmiana stopnia wielomianu aproksymującego wymaga ponownego
rozwiązania układu normalnego.
Metody Numeryczne II
101
Problem złego uwarunkowania
Załóżmy, że punkty x
i
są równo rozmieszczone wewnątrz przedziału [0, 1]. Dla
dużych wartości m:
g
ik
=
m
j=1
x
i+k
j
≈ m
1
0
x
i+k
j
=
m
i + k + 1
,
i, k = 0, . . . , m
(172)
Macierz współczynników układu normalnego:
A =
1
1
2
· · ·
1
m+1
1
2
1
3
· · ·
1
m+2
..
.
..
.
..
.
..
.
1
m+1
1
m+2
· · ·
1
2m+1
(173)
Macierz odwrotna ma z kolei bardzo duże wartości co powoduje bardzo duże błędy
zaokrągleń.
Metody Numeryczne II
102
Aproksymacja za pomocą wielomianów ortogonalnych
Definicja 10 Funkcje f (x) i g(x) nazywamy ortogonalnymi na dyskretnym
zbiorze punktów x
0
, . . . , x
n
gdy
n
i=0
f (x
i
)g(x
i
) = 0
(174)
przy
n
i=0
f
2
(x
i
) > 0,
n
i=0
g
2
(x
i
) > 0.
(175)
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
a
jj
=
n
i=0
ϕ
2
j
(x
i
).
(176)
Metody Numeryczne II
103
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 P
0
, . . . , P
m
będą wielomianami ortogonalnymi o stopniach równych
indeksom oraz Q
m
∈ Π
m
. Wówczas istnieje jednoznaczny rozkład
Q
m
(x) = b
0
P
0
(x) +
· · · + b
m
P
m
(x).
(177)
Współczynnik b
k
łatwo wyznaczyć mnożąc tożsamość (
) przez P
k
(x) i sumując
równości uzyskane dla podstawień x
0
, . . . , x
n
. Otrzymujemy:
n
i=0
Q
m
(x
i
)P
k
(x
i
) = b
0
n
i=0
P
0
(x
i
)P
k
(x
i
) +
· · · + b
k
n
i=0
P
2
k
(x
i
) +
· · ·
· · · + b
m
n
i=0
P
m
(x
i
)P
k
(x
i
)
(178)
Metody Numeryczne II
104
Wobec ortogonalności rodziny
{P
i
} mamy:
n
i=0
Q
m
(x
i
)P
k
(x
i
) = b
k
n
i=0
P
2
k
(x
i
).
(179)
Stąd
b
k
=
n
i=0
Q
m
(x
i
)P
k
(x
i
)
n
i=0
P
2
k
(x
i
)
,
k = 0, . . . , m.
(180)
Szukamy współczynników b
0
, . . . , b
m
minimalizujących
S =
n
i=0
[b
0
P
0
(x
i
) +
· · · + b
m
P
m
(x
i
)
− f(x
i
)]
2
=
n
i=0
m
j=0
b
j
P
j
(x
i
)
2
− 2
m
j=0
b
j
P
j
(x
i
)
f(x
i
) + f
2
(x
i
)
(181)
Metody Numeryczne II
105
W wyniku ortogonalności
{P
i
} oraz podstawień:
s
j
=
n
i=0
P
2
j
(x
i
),
c
j
=
n
i=0
P
j
(x
i
)f (x
i
),
j = 0, . . . , m
(182)
otrzymujemy
S =
n
i=0
m
j=0
b
2
j
P
2
j
(x
i
)
− 2
n
i=0
m
j=0
b
j
P
j
(x
i
)f (x
i
) +
n
i=0
f
2
(x
i
) =
=
m
j=0
b
2
j
s
j
− 2
m
j=0
b
j
c
j
+
n
i=0
f
2
(x
i
) =
m
j=0
!
b
2
j
s
j
− 2b
j
c
j
"
+
n
i=0
f
2
(x
i
)
(183)
Ponieważ
b
2
j
s
j
− 2b
j
c
j
= s
j
b
2
j
− 2b
j
c
j
s
j
= s
j
b
2
j
− 2b
j
c
j
s
j
+
c
2
j
s
2
j
−
c
2
j
s
j
=
= s
j
b
j
−
c
j
s
j
2
−
c
2
j
s
j
(184)
Metody Numeryczne II
106
błąd kwadratowy
S =
m
j=0
s
j
b
j
−
c
j
s
j
2
−
m
j=0
c
2
j
s
j
+
n
i=0
f
2
(x
i
)
(185)
osiąga minimum dla
b
j
=
c
j
s
j
.
(186)
Zatem wielomian aproksymacyjny ma postać
Q
m
(x) =
m
i=0
c
j
s
j
P
j
(x).
(187)
Metody Numeryczne II
107
Przypadek równoodległych punktów
Dla zbioru punktów x
i
= x
0
+ ih, i = 0, . . . , n wprowadzamy q =
x
−x
0
h
(liczby
całkowite) i wyznaczamy odpowiednie unormowane wielomiany ortogonalne:
P
k,n
(t) = 1 + b
1
t + b
2
t(t
− 1) + . . . + b
k
t(t
− 1) · · · (t − k + 1)
(188)
Dowolny wielomian P
j,n
(t) można zapisać w postaci:
P
j,n
(t) =
j
s=0
d
s
(t + s)
[s]
,
(189)
gdzie d
s
są stałymi, a r
[j]
ozn
= r(r
− 1) · · · (r − j + 1). Warunki ortogonalności:
n
i=0
P
j,n
(i)P
k,n
(i) =
n
i=0
j
s=0
d
s
(t + s)
[s]
P
k,n
(i) = 0,
j = 0, . . . , k
− 1, (190)
więc by je spełnić wystarczy by
n
i=0
(i + j)
[j]
P
k,n
(i) = 0,
j = 0, . . . , k
− 1.
(191)
Metody Numeryczne II
108
Sumując
(i + j)
[j]
P
k,n
(i) = (i + j)
[j]
+ b
1
(i + j)
[j+1]
+
· · · + b
k
(i + j)
[j+k]
(192)
po wszystkich i = 0, . . . , n, otrzymujemy
n
i=0
(i + j)
[j]
P
k,n
(i) =
n
i=0
(i + j)
[j]
+ b
1
n
i=0
(i + j)
[j+1]
+
· · ·+b
k
n
i=0
(i + j)
[j+k]
(193)
Wykorzystując, że
n
i=0
(i + j)
[s]
=
(n + j + 1)
[s+1]
s + 1
,
s = j, . . . , j + k,
(194)
i dzieląc obie strony przez (n + j + 1)
[j+1]
dostajemy dla j = 0, . . . , k
− 1
1
1 + j
+
n
2 + j
b
1
+
· · · +
n
[k]
k + 1 + j
b
k
=
Q(j)
(k + 1 + j)
[k]
= 0,
(195)
gdzie Q(j) jest wielomianem stopnia co najwyżej k o zerach j = 0, . . . , k
− 1.
Metody Numeryczne II
109
Przyjmujemy
a
k
= b
k
n
[k]
,
Q(j) = j(j
− 1) · · · (j − k + 1)C
k
= C
k
j
[k]
,
(196)
oraz mnożymy (
) przez (1 + j). Dostajemy
1 +
1 + j
2 + j
a
1
+
· · · +
1 + j
k + 1 + j
a
k
=
C
k
j
[k]
(2 + j)(3 + j)
· · · (k + 1 + j)
,
(197)
skąd dla j =
−1 wynika
C
k
=
1
· 2 · · · k
(
−1)(−2) · · · (−k)
= (
−1)
k
.
(198)
Metody Numeryczne II
110
Aby wyznaczyć a
s
dla s = 1, . . . , k mnożymy (
) przez (s+1+j) i przyjmujemy
j =
−(s + 1). Otrzymujemy:
a
s
=
(
−1)
k
j
[k]
(s + 1 + j)
(k + 1 + j)
[k+1]
=
(
−1)
k
(
−s − 1)(−s − 2) · · · (−s − k)
(k
− s)(k − s − 1) · · · 1 · (−1) · · · (−s)
=
(
−1)
2k
(k + s)
[k]
(k
− s)!(−1)
s
s!
= (
−1)
s
(k + s)
[k]
s!
s!s!(k
− s)!
= (
−1)
s
(k + s)!
s!k!
k!
s!(k
− s)!
= (
−1)
s
k + s
s
k
s
(199)
Zatem ostateczna postać wielomianów ortogonalnych (Czebyszewa/Grama) to:
P
k,n
(t) =
k
s=0
(
−1)
s
k
s
k + s
s
t
[s]
n
[s]
,
k = 0, . . . , m.
(200)
Metody Numeryczne II
111
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) = a
0
ϕ
0
(x) +
· · · + a
n
ϕ
n
(x)
(201)
minimalizujemy wartość funkcji
E
n
=
b
a
ω(x) [P (x)
− f(x)]
2
dx =
b
a
n
i=0
a
i
ϕ
i
(x)
− f(x)
2
dx.
(202)
Rozwiązujemy układ równań:
δE
n
δa
j
=
b
a
n
i=0
a
i
ϕ
i
(x)
− f(x)
ϕ
j
(x)dx
=
n
i=0
a
i
b
a
ϕ
i
(x)ϕ
j
(x)
−
b
a
f (x)ϕ
j
(x)dx = 0,
j = 0, . . . , n.
(203)
Metody Numeryczne II
112
Przykład: Znaleźć najlepszą aproksymację kwadratową funkcji f (x) =
√
x w
przedziale [0, 1] wielomianem stopnia pierwszego Q(x) = a
0
+ a
1
x.
Otrzymujemy układ równań:
a
0
1
0
1
2
dx + a
1
1
0
xdx =
1
0
√
xdx
a
0
1
0
xdx + a
1
1
0
x
2
dx =
1
0
√
xxdx
(204)
Czyli:
a
0
+
1
2
a
1
=
2
3
1
2
a
0
+
1
3
a
1
=
2
5
(205)
Rozwiązanie:
a
0
=
4
15
,
a
1
=
4
5
,
Q(x) =
4
15
+
4
5
x
(206)
Metody Numeryczne II
113
Metoda wyrównywania
• W przypadku pewnych nieliniowych zależności:
– trudno skonstruować odpowiedni wielomian uniwersalny
– trudno rozwiązać wprost układ równań nieliniowych
δE
δa
i
• 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:
y = ax + b
y = ax
b
y = ab
x
y = a +
b
x
y =
1
ax + b
y =
x
ax + b
y = a ln x + b
(207)
Metody Numeryczne II
114
W każdym z wymienionych przypadków zależności nieliniowych aproksymujemy
liniową zależność
Y = αX + β
(208)
gdzie stosujemy następujące podstawienia:
zależność x i y
podstawienia
y = bx
a
X = log x, Y = log y, α = log a
y = ba
x
Y = log y, α = log a, β = log b
y = a +
b
x
Y = xy
y =
1
ax+b
Y =
1
y
y =
x
ax+b
Y =
x
y
y = a ln x + b
X = log x
Brak specyfikacji podstawienia dla danej zmiennej oznacza odpowiednio:
X = x, Y = y, α = a, β = b
(209)
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
115
Aproksymacja funkcjami sklejanymi
Przypadek równoodległych węzłów
Dany jest podział ∆ =
{x
i
: i = 0, . . . , n
} przedziału [a, b] taki, że x
i
= x
0
+ ih,
gdzie h =
b
−a
n
.
Twierdzenie 22 (Baza przestrzeni funkcji sklejanych) Niech funkcje
Φ
3
i
(x), i =
−1, 0, . . . , n, n + 1 będą określone wzorem
Φ
3
i
(x) =
1
h
3
(x
− x
i
−2
)
3
dla x
∈ [x
i
−2
, x
i
−1
]
h
3
+ 3h
2
(x
− x
i
−1
) + 3h(x
− x
i
−1
)
2
− 3(x − x
i
−1
)
3
dla x
∈ [x
i
−1
, x
i
]
h
3
+ 3h
2
(x
i+1
− x) + 3h(x
i+1
− x)
2
− 3(x
i+1
− x)
3
dla x
∈ [x
i
, x
i+1
]
(x
i+2
− x)
3
dla x
∈ [x
i+1
, x
i+2
]
0
w pozostałych przyp.
(210)
Funkcje Φ
i
(x) = Φ
3
i
(x) określone na [a, b] stanowią bazę przestrzeni funkcji sklejanych
trzeciego stopnia na ∆.
Metody Numeryczne II
116
Wniosek 23 Każdą funkcję sklejaną trzeciego stopnia S
∆
(x) można
przedstawić w postaci:
S
∆
(x) =
n+1
i=
−1
c
i
Φ
i
(x),
(211)
gdzie c
i
są liczbami rzeczywistymi.
Własności funkcji bazowych:
Wykres Φ
3
i
(x)
x
i
−2
x
i
−1
x
i
x
i+1
x
i+2
Wartości Φ
3
i
(x) i pochodnych
x
x
i
−2
x
i
−1
x
i
x
i+1
x
i+2
Φ
3
i
(x)
0
1
4
1
0
!
Φ
3
i
(x)
"
0
3
h
0
−
3
h
0
!
Φ
3
i
(x)
"
0
6
h
2
−
12
h
2
6
h
2
0
Są klasy C
2
((
−∞, ∞))
Metody Numeryczne II
117
Aproksymacja funkcjami sklejanymi z równoodległymi węzłami
Minimalizujemy:
E =
b
a
f (x)
−
n+1
i=
−1
c
i
Φ
3
i
(x)
2
dx
(212)
Minimum wyznacza rozwiązanie układu n + 3 równań z n + 3 niewiadomymi:
δE
δc
j
= 0,
j =
−1, . . . , n + 1
(213)
czyli:
n+1
i=
−1
c
i
b
a
Φ
3
i
(x)Φ
3
j
(x) =
b
a
f (x)Φ
3
j
(x),
j =
−1, . . . , n + 1
(214)
Wprowadzając a
ij
=
b
a
Φ
3
i
(x)Φ
3
j
(x) otrzymujemy układ równań o symetrycznej
pięciodiagonalnej macierzy głównej.
Metody Numeryczne II
118
Aproksymacja funkcjami sklejanymi dla dyskretnego zbioru punktów
Odcięte węzłów:
{z
i
, i = 0, . . . , m
}, gdzie m > n + 3.
J =
m
k=0
f (z
k
)
−
n+1
i=
−1
c
i
Φ
3
i
(z
k
)
2
(215)
Z przyrównania pochodnych cząstkowych do zera mamy układ:
n+1
i=
−1
b
ij
c
i
m
k=0
f (z
k
)Φ
3
j
(z
k
),
j =
−1, . . . , n + 1
(216)
gdzie
b
ij
=
m
k=0
Φ
3
i
(z
k
)Φ
3
j
(z
k
)
(217)
jest symetryczną pięciodiagonalną macierzą.
Jeśli wyznacznik układu jest różny od zera, to jego rozwiązanie określa minimum
funkcji J.
Metody Numeryczne II
119
Aproksymacja jednostajna
Funkcja błędu:
E = sup
x
∈[a,b]
|f(x) − F (x)|
(218)
• Z twierdzenia Weierstrassa (tw.
) 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 W
n
(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
120
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:
W
n
(x) = f (x
0
) +
f
(x
0
)
1!
(x
−x
0
) +
f
(x
0
)
2!
(x
−x
0
)
2
+
· · ·+
f
(n)
(x
0
)
n!
(x
−x
0
)
n
(219)
Oszacowanie błędu z reszty w postaci Lagrange’a:
|f(x) − W
n
(x)
| =
f
(n+1)
(c)
(n + 1)!
(x
− x
0
)
n+1
M
n+1
(n + 1)!
δ
n+1
(220)
dla c leżącego między x i x
0
, δ będącego promieniem stosownego otoczenia
punktu x
0
oraz M
n+1
|f
(n+1)
(x)
| dla x ∈ [a, b].
Metody Numeryczne II
121
Przybliżenia Padégo
Dokładniejsze od wielomianowych często okazują się być wymierne:
R
n,k
(x) =
L
n
(x)
M
k
(x)
,
(221)
gdzie
L
n
(x) = a
0
+ a
1
x +
· · · + a
n
x
n
M
k
(x) = b
0
+ b
1
x +
· · · + b
k
x
k
(222)
• Parametry określa się tak, by w x
0
= 0 wartości funkcji f i R
n,k
oraz jak
najwięcej ich pochodnych zrównało się.
• 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
122
Szeregi Czebyszewa
Przybliżanie funkcji f (x) sumą częściową szeregu
f (x)
≈
1
2
c
0
+
n
j=1
c
j
T
j
(x),
(223)
gdzie
c
j
=
2
π
1
−1
f (x)T
j
(x)
√
1
− x
2
dx
T
j
(x) = cos(j arccos x).
(224)
Alternatywnie, funkcję można aproksymować funkcją wymierną
T
n,k
(x) =
n
i=0
a
i
T
i
(x)
k
i=0
b
i
T
i
(x)
.
(225)
Wartości a
i
, b
i
wyznacza się tak, by w różnicy f (x)
− T
n,k
(x) =
1
2
c
0
+
∞
j=1
c
j
T
j
(x)
− T
n,k
(x) znikały współczynniki przy T
j
dla j = 0, . . . , N .
Metody Numeryczne II
123
Rozwiązywanie równań nieliniowych
Zadanie: Najczęstsze postaci:
1. Znaleźć pierwiastki rzeczywiste równania f (x) = 0, funkcji f zmiennej
rzeczywistej.
2. Znaleźć pierwiastki rzeczywiste i zespolone równania P (x) = 0 wielomianu P .
Metody iteracyjne
• Ciąg kolejnych przybliżeń {x
i
}.
• n-punktowy wzór iteracyjny
x
i+1
= F
i
(x
i
, x
i
−1
, . . . , x
i
−n+1
).
(226)
• Funkcja F
i
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 (α, . . . , α).
(227)
Metody Numeryczne II
124
• Niech α będzie szukanym pierwiastkiem. Błąd i-tej iteracji ε
i
= α
− x
i
. Jeśli
lim
i
→∞
x
i
= α, to mówimy, że metoda jest w punkcie α rzędu p
1 taka, że
lim
i
→∞
|x
i+1
− α|
|x
i
− α|
p
= lim
i
→∞
|ε
i+1
|
|ε
i
|
p
= C
= 0.
(228)
Wówczas stałą C nazywamy stałą asymptotyczną błędu.
• Efektywność metody - iloczyn kosztu ϑ jednej iteracji i liczby iteracji I.
Niech rzędy dwóch metod będą odpowiednio równe p
1
i p
2
. Dla dużych i
mamy w przybliżeniu:
|ε
(1)
i+1
| = C
1
|ε
(1)
i
|
p
1
,
|ε
(2)
i+1
| = C
2
|ε
(2)
i
|
p
2
.
(229)
Niech
S
i
=
− log |ε
(1)
i
|,
T
i
=
− log |ε
(2)
i
|.
(230)
Wtedy
S
i+1
=
− log C
1
+ p
1
S
i
,
T
i+1
=
− log C
2
+ p
2
T
i
.
(231)
Metody Numeryczne II
125
Rozwiązania:
S
i
= S
1
p
i
1
− log C
(p
i
1
−1)/(p
1
−1)
1
,
T
i
= T
1
p
i
2
− log C
(p
i
2
−1)/(p
2
−1)
2
.
(232)
Zakładamy, że S
1
= T
1
, koszt jednej iteracji to odpowiednio ϑ
1
i ϑ
2
, a liczby
kroków potrzebne do uzyskania odpowiedniej dokładności to I
1
oraz I
2
.
Zatem S
I
1
i T
I
2
sa w przybliżeniu równe. Stąd:
S
1
(p
I
1
1
− p
I
2
2
) + log
C
(p
I2
2
−1)/(p
2
−1)
2
C
(p
I1
1
−1)/(p
1
−1)
1
= 0.
(233)
Najczęściej pierwszy składnik dominuje, więc można drugi zaniedbać. Zatem:
I
1
log p
1
≈ I
2
log p
2
≈ const.
(234)
Ostatecznie wskaźnik efektywności można określić jako:
E =
1
ϑ
log p
albo
E = p
1
ϑ
.
(235)
• Rząd metody jest cechą lokalną.
• Koszt jednej iteracji ≈ koszt liczenia wartości funkcji i jej pochodnych.
Metody Numeryczne II
126
Metoda interpolacji odwrotnej
Zakładamy, że w okolicy pierwiastka α funkcji f istnieje funkcja odwrotna g = f
−1
(α jest pierwiastkiem pojedynczym).
Niech x
i
−n
, . . . , x
i
będą kolejnymi przybliżeniami α (bądź pewnymi punktami z
otoczenia α dla i = 0). Przyjmujemy
y
j
= f (x
i
−j
),
j = 0, . . . , n.
(236)
Możemy interpolować funkcję g(y) wielomianem interpolacyjnym Lagrange’a
h(y) =
n
j=0
L
j
(y)g(y
j
) =
n
j=0
L
j
(y)x
i
−j
.
(237)
Metody Numeryczne II
127
Wówczas kolejnym przybliżeniem x
i+1
pierwiastka funkcji f jest:
h(0) =
n
j=0
L
j
(0)x
i
−j
=
n
j=0
x
i
−j
(
−1)
n
y
0
· · · y
n
(y
j
− y
0
)
· · · (y
j
− y
j
−1
)(y
j
− y
j+1
)
· · · (y
j
− y
n
)
.
(238)
Ze wzoru Lagrange’a:
α
− x
i+1
=
g
(n)
(η)
n!
(
−1)
n+1
y
0
· · · y
n
,
η
∈ I[y
0
, . . . , y
n
, 0].
(239)
• Można użyć dowolnej metody interpolacji.
Metody Numeryczne II
128
Metoda bisekcji
Szukamy zera α w przedziale (a, b), takim, że f (a)f (b) < 0. Przyjmujemy
a
0
= a, b
0
= b, a także
x
i
=
a
i
+ b
i
2
,
i = 0, 1, . . . ,
(240)
oraz
(a
i+1
, b
i+1
) =
(a
i
, x
i
)
o ile f (a
i
)f (x
i
) < 0
(x
i
, b
i
)
o ile f (x
i
)f (b
i
) < 0
.
(241)
Zakładamy, że f (x
i
)
= 0, bo inaczej mamy rozwiązanie dokładne.
Ponieważ
|x
i
− x
i+1
| =
1
2
i
(b
− a),
i = 0, 1, . . . ,
(242)
to
lim
i
→∞
x
i
= α.
(243)
Metody Numeryczne II
129
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 (a
0
, b
0
) = (a, b). W i-tym kroku prowadzimy cięciwę o równaniu
y
− f(a
i
) =
f (b
i
)
− f(a
i
)
b
i
− a
i
(x
− a
i
)
(244)
i przybliżamy pierwiastek równania f (x) = 0 pierwiastkiem cięciwy:
x
i
= a
i
−
f (a
i
)
f (b
i
)
− f(a
i
)
(b
i
− a
i
).
(245)
W kolejnych krokach
(a
i+1
, b
i+1
) =
(a
i
, x
i
)
o ile f (a
i
)f (x
i
) < 0
(x
i
, b
i
)
o ile f (x
i
)f (b
i
) < 0
.
(246)
• na ogół niestacjonarna, stacjonarna np. dla funkcji wypukłych
Metody Numeryczne II
130
Oszacowania błędów metody regula falsi
Z twierdzenia Lagrange’a o przyrostach:
f (x
n
)
− f(α) = f
(c)(x
n
− α),
dla pewnego c
∈ I[x
n
, α]
(247)
Ponieważ f (α) = 0, więc
|x
n
− α|
|f(x
n
)
|
m
,
m = inf
x
∈[a,b]
|f
(x)
|
(248)
Przy założeniu, żę f
i f
nie zmieniają znaku w [a, b] jeden z końców przedziału
nie zmienia się w kolejnych krokach iteracji. Przyjmijmy f
(x) > 0, f
(x) > 0. Stąd
b
i
= b, i = 0, 1, . . . czyli
x
−1
= a,
x
i+1
= x
i
−
f (x
i
)
f (b)
− f(x
i
)
(b
− x
i
),
i = 0, 1, . . .
(249)
Mamy więc
f (α)
− f(x
i
) =
f (x
i
)
− f(b)
x
i
− b
(x
i+1
− x
i
),
(250)
Metody Numeryczne II
131
i z tw. Lagrange’a
(α
− x
i
)f
(ξ
i
) = (x
i+1
− x
i
)f
(x
i
),
ξ
i
∈ (x
i
, α), x
i
∈ (x
i
, b).
(251)
Po wymnożeniu i dodaniu obustronnie
−x
i+1
f
(ξ
i
) otrzymujemy
|α − x
i+1
| =
|f
(x
i
)
− f
(ξ
i
)
|
|f
(ξ
i
)
|
|x
i+1
− x
i
|.
(252)
Ponieważ ξ
i
, x
i
∈ [a, b] i f
(x) > 0, więc
|f
(x
i
)
− f
(ξ
i
)
| M − m,
m =
inf
x
∈[a,b]
|f
(x)
|, M = sup
x
∈[a,b]
|f
(x)
|.
(253)
Ostatecznie
|α − x
i+1
|
M
− m
m
|x
i+1
− x
i
|.
(254)
Oszacowanie to jest na ogół pesymistyczne.
Dla przybliżeń w niewielkim otoczeniu α można szacować:
|α − x
i+1
| ≈
f (x
i+1
)
f
(x
i+1
)
≈
x
i+1
− x
i
f (x
i+1
)
− f(x
i
)
|f(x
i+1
)
|
(255)
Metody Numeryczne II
132
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:
x
0
= a, x
1
= b,
x
i+1
= x
i
+
f (x
i
)(x
i
− x
i
−1
)
f (x
i
)
− f(x
i
−1
)
,
i = 1, 2, . . .
(256)
• Istotne znaczenie ma maksymalna graniczna dokładność (wynikająca z
przyjętej arytmetyki) – gdy x
i
− x
i
−1
jest tego samego rzędu co oszacowanie
błędu, to kolejne przybliżenia mogą być całkowicie błędne. Dodatkowe
kryterium stopu: zaburzenie zmniejszania się
|f(x
i
)
|.
• Wykrywanie rozbieżności – różnica pomiędzy kolejnymi przybliżeniami rośnie.
Metody Numeryczne II
133
Metoda Newtona
• Założenia: f(a)f(b) < 0, f
oraz f
nie zmieniają znaku na [a, b].
• Z tego końca przedziału x
0
∈ {a, b}, w którym f(x) ma ten sam znak co
f
(x) prowadzimy styczną do wykresu funkcji y = f (x). Punkt x
1
przecięcia
stycznej z osią odciętych to pierwsze przybliżenie pierwiastka.
• Prowadzimy styczną w punkcie x
1
i wyznaczamy punkt x
2
przecięcia tej
stycznej z osią.
• Uzyskany ciąg x
1
, x
2
, . . . jest zbieżny monotonicznie do α.
Dowód:
1. Ciąg jest monotoniczny.
2. Zbiega do α.
Przypadek f
(x) > 0, f
(x) > 0:
Wówczas x
0
= b oraz f (x
0
) > 0. Pokażemy, że x
n+1
< x
n
dla n = 0, 1, . . ..
Równanie stycznej:
y
− f(x
n
) = f
(x
n
)(x
− x
n
).
(257)
Metody Numeryczne II
134
Zero stycznej w punkcie
x
n+1
= x
n
−
f (x
n
)
f
(x
n
)
.
(258)
Ze wzoru Taylora:
0 = f (α) = f (x
n
) + f
(x
n
)(α
− x
n
) +
1
2
f
(c)(α
− x
n
)
2
,
c
∈ (α, x
n
). (259)
Stąd
α = x
n
−
f (x
n
)
f
(x
n
)
−
1
2
f
(c)
f
(x
n
)
(α
− x
n
)
2
.
(260)
Z (
):
α
− x
n+1
=
−
1
2
f
(c)
f
(x
n
)
(α
− x
n
)
2
< 0
(261)
Zatem f (x
n+1
) > 0 oraz x
n+1
< x
n
, bo
x
n+1
− x
n
=
−
f (x
n
)
f
(x
n
)
< 0.
(262)
Metody Numeryczne II
135
Ciąg x
n
: n = 0, 1, . . . jest więc monotoniczny i ograniczony, a zatem zbieżny
(do pewnego g). Z (
) przy n
→ ∞ mamy
g = g
−
f (g)
f
(g)
,
(263)
czyli f (g) = 0, a więc g = α.
Błąd przybliżenia x
n
.
Podobnie jak w metodzie regula falsi (str.
|x
n
− α|
|f(x
n
)
|
m
,
m = inf
x
∈[a,b]
|f
(x)
|
(264)
Ze wzoru Taylora, dla pewnego ξ
n
∈ I[x
n
, x
n+1
]:
f (x
n+1
) = f (x
n
) + f
(x
n
)(x
n+1
− x
n
) +
1
2
f
(ξ
n
)(x
n+1
− x
n
)
2
(265)
Metody Numeryczne II
136
Z (
) f (x
n
) + f
(x
n
)(x
n+1
− x
n
) = 0, więc
|f(x
n+1
)
|
1
2
M
m
(x
n+1
− x
n
)
2
,
M = sup
x
∈[a,b]
|f
(x)
|.
(266)
Zatem
|α − x
n+1
|
1
2
M (x
n+1
− x
n
)
2
=
M
2m
f (x
n
)
f
(x
n
)
2
.
(267)
Alternatywne oszacowanie (w bliskim sąsiedztwie α):
|α − x
n
| ≈
f (x
n
)
f
(x
n
)
.
(268)
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ń x
2
− c = 0.
• Metoda Newtona zwykle potrzebuje mniej iteracji niż metoda siecznych, ale
ma wyższy koszt jednej iteracji.
Metody Numeryczne II
137
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 f
(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 parzystokrotnych obniżony rząd metody
Metody Numeryczne II
138
Przykłady modyfikacji metod
• Znamy krotność pierwiastka:
x
n+1
= x
n
− r
f (x
n
)
f
(x
n
)
(269)
• Krotność nieznana. Zamiast funkcji f(x) używamy funkcji
u(x) =
f (x)
f
(x)
,
(270)
która ma te same pierwiastki co f (x), ale wszystkie pojedyncze. Wówczas dla
metod regula falsi i siecznych wzory (
) przyjmują postać
x
n+1
= x
n
− u(x
n
)
x
n
− x
n
−1
u(x
n
)
− u(x
n
−1
)
,
(271)
a dla metody Newtona mamy:
x
n+1
= x
n
−
u(x)
u
(x)
,
przy u
(x
n
) = 1
−
f
(x
n
)
f (x
n
)
u(x
n
).
(272)
Metody Numeryczne II
139
Zera wielomianów
Liczba zer wielomianu
Definicja 12 Ciągiem Sturma dla wielomianu p(x) nazywamy ciąg
wielomianów p
0
(x), . . . , p
m
(x) taki, że:
1. p
0
(x) = p(x).
2. p
1
(x) = p
0
(x).
3. dla i = 2, . . . , m + 1 p
i
(x) jest resztą z dzielenia p
i
−2
przez p
i
−1
wziętą ze
znakiem przeciwnym.
4. p
m+1
(x)
≡ 0
Oznaczenie: N (z)
ozn
= liczba zmian znaku w ciągu
{p
i
(z) : i = 0, 1 . . . , p
i
(z)
= 0}.
Twierdzenie 24 (Sturma) Jeżeli ciąg
{p
i
(x) : i = 0, . . . , m
} jest ciągiem
Sturma dla wielomianu p(x) na przedziale (a, b) oraz p
0
(a)p
0
(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
140
Oznaczenie: M (z)
ozn
= 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) = a
0
x
n
+ . . . + a
n
−1
x + a
n
, gdzie a
0
= 0,
L(z)
ozn
= liczba zmian znaku ciągu
{p
k
(z) =
k
i=0
a
i
z
k
−i
: k = 0, . . . , n
}.
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
141
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
p
1
(x) = x
n
p
1
x
,
p
2
(x) = p(
−x),
p
1
(x) = x
n
p
−1
x
(273)
mają kresy górne odpowiednio R
1
, R
2
i R
3
, to wszystkie dodatnie zera
wielomianu p(x) leżą w (
1
R
1
, R), a ujemne w (
−R
2
,
−
1
R
3
)
Twierdzenie 27 (Lagrange’a) Niech a
0
= 0 i a
k
(k
1) będzie pierwszym
ujemnym współczynnikiem wielomianu p(x) = a
0
x
n
+ . . . + a
n
−1
x + a
n
.
Wszystkie zera tego wielomianu są mniejsze niż
R = 1 +
k
#
A
|a
0
|
,
(274)
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
142
Wyznaczanie zer wielomianów
• Znamy kres górny - można metodą Newtona znaleźć największy pierwiastek.
• Metoda Newtona może wolno zbiegać do największego zera. Metoda
podwójnego kroku:
x
k+1
= x
k
− 2
p(x
k
)
p
(x
k
)
(275)
Twierdzenie 28 Niech p(x)
∈ Π
n
, n
2 ma wszystkie zera rzeczywiste
α
1
α
2
· · · α
n
. Niech β
1
będzie największym zerem wielomianu p
(x)
(α
1
β
1
α
2
). W przypadku n = 2 dodatkowo zakładamy α
1
> α
2
.
Dla każdego z > α
1
dobrze określone są liczby
z
= z
−
p(z)
p
(z)
,
y = z
− 2
p(z)
p
(z)
,
y
= y
−
p(y)
p
(y)
(276)
i spełniają nierówności
β
1
< y,
α
1
y
z
.
(277)
Metody Numeryczne II
143
Dla ciągu przybliżeń
{x
n
: i = 1, 2, . . .
} (zgodnych z (
)) mamy:
1.
∀
n
x
n
α
1
oraz lim
n
→∞
x
n
= α
1
.
2.
∃
k
0
p(x
k
0
)p(x
0
) < 0 oraz
∀
k<k
0
p(x
k
)p(x
0
) > 0. Wówczas ciąg
y
0
= x
k
0
,
y
k+1
= y
k
−
p(y
k
)
p
(y
k
)
, k = 1, 2, . . .
(278)
zbiega do α
1
.
Deflacja – „oddzielamy” pierwiastek α
1
od pozostałych – wyznaczamy wielomian
stopnia n
− 1
p
1
(x) =
p(x)
x
− α
1
.
(279)
Największym zerem wielomianu p
1
(x) jest α
2
, a dobrym punktem startowym jest
ewentualnie „przestrzelone” x
k
.
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
144
Źle uwarunkowane równania algebraiczne
• pierwiastki wielokrotne są na ogół źle uwarunkowane
• bliskie sobie pierwiastki
• nawet „wyraźnie rozdzielone” pierwiastki:
Przykład: Wielomian o pierwiastkach rzeczywistych 1,2,. . . ,20:
p(z) = (z
− 1)(z − 2) · · · (z − 20) = z
20
− 210z
19
+
· · · + 20!
(280)
Jeśli zmniejszyć a
1
o 2
−23
≈ 10
−10
, to równanie ma pierwiastki:
16.730737466
± 2.812624894i
(281)
Dla schematu Hornera błąd obliczonej wartości p(x) może być równy
E =
n
i=0
|(2i + 1)a
n
−i
x
i
|1.06
(282)
Błąd względny pierwiastka może osiągać
E
αp
(α)
Metody Numeryczne II
145
Stąd wskaźnik uwarunkowania:
C =
|(2i + 1)a
n
−i
α
i
|
|αp
(α)
|
=
|(2i + 1)a
n
−i
α
i
|
|ia
n
−i
α
i
|
(283)
Dla badanego wielomianu:
C
≈ 1.35 · 10
15
,
dla α = 14.
(284)
Gdy pierwiastkami są liczby α
i
= 2
i
−21
, i = 1, 2, . . . , 20, to wskaźnik C = 3.83
· 10
3
.
Miara rozdzielenia względnego: Dla k pierwiastków rzeczywistych
α
i
< . . . < α
k
:
s =
k
−1
i=1
α
i+1
− α
i
|α
i
|
α
i
= 0
(285)
W przypadku pierwiastków 1, 2, . . . , 20 mamy s = 8.2
· 10
−18
. Dla pierwiastków
α
i
= 2
i
−21
, i = 1, 2, . . . , 20 mamy s = 1.
Metody Numeryczne II
146
Algorytm Maehly’ego dla wyeliminowania deflacji:
Ponieważ
p
1
(x) =
p
(x)
x
− α
1
−
p(x)
(x
− α
1
)
2
,
(286)
więc
x
k+1
= x
k
−
p
1
(x
k
)
p
1
(x
k
)
= x
k
−
p(x
k
)
p
(x
k
)
−
p(x
k
)
x
k
− α
1
.
(287)
Ogólnie mamy
p
j
(x) =
p(x)
(x
− α
1
)
· · · (x − α
j
)
(288)
p
j
(x) =
p
(x)
(x
− α
1
)
· · · (x − α
j
)
−
p(x)
(x
− α
1
)
· · · (x − α
j
)
j
i=1
1
x
− α
i
.
(289)
Stąd iterujemy
x
k+1
= x
k
−
p(x
k
)
p
(x
k
)
−
j
i=1
p(x
k
)
x
k
−α
i
.
(290)
Metody Numeryczne II
147
Metoda δ
2
Dla metod iteracyjnych rzędu 1 mamy dla dużych wartości k
α
− x
k+2
α
− x
k+1
≈
α
− x
k+1
α
− x
k
≈ C.
(291)
Stąd wyznaczamy
α
≈
x
k
x
k+2
− x
2
k+1
x
k+2
− 2x
k+1
+ x
k
= x
k+2
−
(∆x
k+1
)
2
∆
2
x
k
,
(292)
gdzie ∆ jest operatorem różnicy progresywnej:
∆x
k
= ∆
1
x
k
= x
k+1
− x
k
,
∆
i+1
x
k
= ∆
i
x
k+1
− ∆
i
x
k
(293)
Uwagi:
• proces δ
2
Aitkena
• znacznie lepsze przybliżenie niż przez x
k
• tylko dla liniowej zbieżności
Metody Numeryczne II
148
Minima funkcji jednej zmiennej
Sposoby znajdowania:
• szukanie rozwiązań równania f
(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 [a
, b
] 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 < t
1
< t
2
< b. Jeśli f (t
1
) < f (t
2
), to α
∈ [a, t
2
]. W przec. przyp.
α
∈ [t
1
, b].
Metody podziału polegają na wyznaczaniu ciągów przedziałów zstępujących
[a
i
, b
i
], i = 0, 1, . . . , a
0
= a, b
0
= b oraz b
i
− a
i
n
→∞
→ 0.
Metody Numeryczne II
149
Metoda podziału na trzy części
t
i
1
=
2
3
a +
1
3
b,
t
i
2
=
1
3
a +
2
3
b
(294)
b
i
− a
i
=
2
3
i
(b
− a)
(295)
W każdej iteracji liczymy 2 wartości funkcji.
Metoda połowienia
t
i
1
=
3
4
a +
1
4
b,
t
i
2
=
1
2
a +
1
2
b,
t
i
3
=
1
4
a +
3
4
b
(296)
b
i
− a
i
=
1
2
i
(b
− a)
(297)
W każdej iteracji liczymy 3 lub 2 wartości funkcji.
Metody Numeryczne II
150
Metoda optymalnych podziałów – Johnsona
Cel: Minimalizacja liczby wyliczeń wartości funkcji.
Wykorzystujemy liczby Fibonacciego:
F
0
= F
1
= 1,
F
i
= F
i
−1
+ F
i
−2
,
i = 2, 3, . . .
(298)
Algorytm: Dla zadanej dokładności ρ
• c =
b
−a
ρ
, a
0
= a, b
0
= b
• Znajdujemy N takie, że F
N
−1
< c < F
N
• Dla i=1,. . . ,N-2:
t
i
1
= a
i
−1
+
F
N
−i−1
F
N
−i+1
(b
i
−1
− a
i
−1
),
t
i
2
= a
i
−1
+
F
N
−i
F
N
−i+1
(b
i
−1
− a
i
−1
) (299)
Jeżeli f (t
i
1
)
f(t
i
2
), to a
i
−1
← a, b
i
−1
← t
i
2
. W przec. przyp. a
← t
i
1
, b
i
−1
← b
Koszty: Liczymy wartości funkcji w N
− 1 punktach otrzymując ostatecznie
b
N
−2
− a
N
−2
=
F
N
−1
F
N
F
N
−2
F
N
−1
· · ·
F
2
F
3
(b
− a) = 2
b
− a
F
N
2ρ
(300)
Metody Numeryczne II
151
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.
t
i
1
oraz t
i
2
wybieramy tak aby:
t
i
2
− a
i
= b
i
− t
i
1
= τ (b
i
− a
i
),
τ
∈ (0, 1)
b
− t
i
2
= τ (b
− t
i
1
)
(301)
Zatem: τ
2
+ τ
− 1 = 0 czyli τ =
√
5
−1
2
≈ 0, 62 – stosunek boków „złotego”
prostokąta.
Metody Numeryczne II
152
Algorytm:
• a
0
← a, b
0
← b, t
1
1
= a + (1
− τ)(b − a), t
1
2
= b
− (1 − τ)(b − a)
• Dla kolejnych i:
Jeżeli f (t
i
1
) > f (t
i
2
), to
a
i+1
← a
i
,
b
i+1
← t
i
2
t
i+1
2
← t
i
1
,
t
i+1
1
← a
i
+ (1
− τ)(b
i
− a
i
)
(302)
W przeciwnym przypadku:
a
i+1
← t
i
1
,
b
i+1
← b
t
i+1
1
← t
i
2
,
t
i+1
2
← b
i
− (1 − τ)(b
i
− a
i
)
(303)
Metody Numeryczne II
153
Porównanie metod Johnsona i złotego podziału
Przy założeniu K obliczeń funkcji:
Metoda Johnsona:
N = K + 1,
|t − α|
b
− a
F
K+1
(304)
Metoda złotego podziału:
|t − α|
b
− a
2G
K
−1
, gdzie G
i
=
1
τ
i
−1
.
(305)
Porównanie:
2G
K
−1
< F
K+1
< 2G
K
(306)
Metody Numeryczne II
154
Układy równań liniowych
Zadanie: rozwiązać układ
Ax = b,
A =
a
11
. . .
a
1n
..
.
..
.
a
n1
. . .
a
nn
, b =
b
1
..
.
b
n
(307)
Zadanie równoważne znalezieniu macierzy odwrotnej A
−1
, bo
x = A
−1
b
(308)
Z drugiej strony: i-ta kolumna a
i
macierzy A
−1
jest rozwiązaniem układu
Ax = e
i
(309)
Metody Numeryczne II
155
Eliminacja Gaussa
Układ Ax = b przekształcamy do
Rx = c,
R =
r
11
. . .
r
1n
. ..
..
.
0
r
nn
, c =
c
1
..
.
c
n
,
(310)
tak by rozwiązanie wyznaczyć jako:
x
i
=
c
i
−
n
k=i+1
r
ik
x
k
r
ii
,
i = n, n
− 1, . . . , 1.
(311)
Metody Numeryczne II
156
Pierwszy krok: Sprowadzić macierz (A, b) do macierzy (A
, b
):
(A, b) =
a
11
. . .
a
1n
b
1
..
.
..
.
..
.
a
n1
. . .
a
nn
b
n
, (A
, b
) =
a
11
a
12
. . .
a
1n
b
1
0
a
22
. . .
a
2n
b
2
..
.
..
.
..
.
..
.
0
a
n2
. . .
a
nn
b
n
(312)
tak, aby układy Ax = b i A
x = b
miały te same rozwiązania.
Jeśli a
11
= 0, to wystarczy wyliczyć
a
ij
= a
ij
− a
1j
a
i1
a
11
,
i = 2, . . . , n, j = 1, . . . , n.
(313)
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 (A
, b
) bez pierwszego
wiersza i pierwszej kolumny.
Metody Numeryczne II
157
Eliminacja Gaussa w postaci macierzowej:
(A, b) = (A
0
, b
0
)
→ (A
1
, b
1
)
→ · · · → (A
n
−1
, b
n
−1
) = (R, c),
(314)
gdzie
(A
j
, b
j
) = G
j
P
j
(A
j
−1
, b
j
−1
)
(R, c) = G
n
−1
P
n
−1
· · · G
1
P
1
(A, b),
(315)
G
j
=
1
0
. .
.
1
−l
j+1,j
1
..
.
. .
.
0
−l
n,j
0
1
,
P
j
=
1
. .
.
0
1
. .
.
1
0
. .
.
1
(316)
Metody Numeryczne II
158
Rozkład na iloczyn macierzy trójkątnych
T
0
= (A, b).
Dla j = 1, . . . , n:
T
j
=
r
11
r
12
. . .
r
1
j
r
1
, j + 1
. . .
r
1n
c
1
λ
21
r
22
. . .
r
2
j
..
.
..
.
..
.
λ
31
λ
32
..
.
..
.
..
.
..
.
..
.
..
.
r
jj
r
j,j+1
. . .
r
jn
c
j
..
.
..
.
λ
j+1,j
a
j
j+1,j+1
. . .
a
j
j+1,n
b
j
j+1
..
.
..
.
..
.
..
.
..
.
..
.
λ
n,1
λ
n2
. . .
λ
nj
a
j
n,j+1
. . .
a
j
nn
b
j
n
.
(317)
Kolumny z λ
ij
są permutacją odpowiednich l
ij
.
Metody Numeryczne II
159
Otrzymujemy:
LR = P
n
−1
P
n
−2
. . . P
1
A = PA,
(318)
gdzie
L =
1
0
t
n
21
. ..
..
.
. ..
t
n
n1
. . .
t
n
n,n
−1
1
,
R =
t
n
11
. . .
t
n
1n
. ..
0
t
n
nn
.
(319)
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.
(320)
• Nie każda nieosobliwa macierz ma rozkład trójkątny A = LR.
• Złożoność O(n
3
).
Metody Numeryczne II
160
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 n
2
równań z n(n + 1) niewiadomymi:
a
ij
=
min(i,j)
p=1
l
ip
r
pj
.
(321)
W k-tym kroku:
a
kj
=
k
p=1
l
kp
r
pj
,
j
k,
a
ik
=
k
p=1
l
ip
r
pk
,
i > k.
(322)
Metody Numeryczne II
161
Metoda Doolittle’a: l
kk
= 1,
k = 1, . . . , n
r
kj
= a
kj
−
k
−1
p=1
l
kp
r
pj
,
j = k, . . . , n,
l
ik
=
a
ik
−
k
−1
p=1
l
ip
r
pk
r
kk
,
i = k + 1, . . . , n.
(323)
Metoda Crouta: r
kk
= 1,
k = 1, . . . , n
l
ik
= a
ik
−
k
−1
p=1
l
ip
r
pk
,
i = k, . . . , n,
r
kj
=
a
kj
−
k
−1
p=1
l
kp
r
pj
l
kk
,
j = k + 1, . . . , n.
(324)
Metody Numeryczne II
162
Algorytm Gaussa–Jordana
Cel: Wyznaczenie macierzy odwrotnej A
−1
. Macierz A realizuje odwzorowanie:
a
11
x
1
+
· · · + a
1n
x
n
= y
1
..
.
a
n1
x
1
+
· · · + a
nn
x
n
= y
n
(325)
Wyznaczamy odwzorowanie odwrotne poprzez zamiany zmiennych. Zmienną x
1
zamieniamy na takie y
r
, dla którego
|a
r1
| = max
i
|a
i1
|.
(326)
Jeśli a
r1
= 0, to macierz A jest osobliwa.
Zamieniamy miejscami równania 1 i r otrzymując układ Ax = y, gdzie a
11
= a
r1
oraz y
1
= y
r
.
Metody Numeryczne II
163
Po zamianie otrzymujemy układ:
a
11
y
1
+ a
12
x
2
+
· · · + a
1n
x
n
= x
1
a
21
y
1
+ a
22
x
2
+
· · · + a
2n
x
n
= y
2
..
.
a
n1
y
1
+ a
n2
x
2
+
· · · + a
nn
x
n
= y
n
,
(327)
gdzie dla i, j = 2, . . . , n
a
11
=
1
a
11
,
a
1j
=
−
a
1j
a
11
,
a
i1
=
a
i1
a
11
,
a
ij
= a
ij
−
a
i1
a
1j
a
11
.
(328)
W każdym kolejnym kroku k = 2, . . . , n zamieniamy w analogiczny sposób
zmienną x
k
przez jedną ze zmiennych y
k
, . . . , y
n
.
Metody Numeryczne II
164
Ogólnie: konstruujemy ciąg macierzy:
A = A
0
→ · · · → A
n
= A
−1
,
(329)
gdzie dla k = 1, . . . , n macierz A
k
= (a
ij
) powstaje z A
k
−1
= (a
ij
) w następujący
sposób:
1. Wyznaczamy r takie, że
|a
rk
| = max
i
k
|a
ik
|.
(330)
Jeśli a
rk
= 0, to macierz A jest osobliwa – Koniec!.
2. Zamieniamy wiersze k i r macierzy A
k
−1
otrzymując macierz A = (a
ik
).
3. Dla i, j
= k liczymy:
a
kk
=
1
a
kk
,
a
kj
=
−
a
kj
a
kk
,
a
ik
=
a
ik
a
kk
,
a
ij
= a
ij
−
a
ik
a
kj
a
kk
.
(331)
Metody Numeryczne II
165
Rozkład Choleskiego
Definicja 13 Macierz zespolona A wymiaru n
× n, nazywamy dodatnio
określoną, gdy:
1. A = A
H
, tzn. A jest macierzą hermitowską.
2. Dla wszystkich x
∈ C
n
, x
= 0 zachodzi x
H
Ax > 0.
Hermitowską macierz nazywamy dodatnio półokreśloną gdy dla
wszystkich x
∈ C
n
zachodzi x
H
Ax
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 l
kk
> 0, k = 1, . . . , n, oraz
A = LL
H
.
(332)
Jeżeli A jest rzeczywista, to również L jest rzeczywista.
Metody Numeryczne II
166
Równanie (
) oznacza, że dla k = 1, . . . , n, i = k, . . . , n zachodzi:
a
ik
= l
i1
l
k1
+ l
i2
l
k2
+
· · · + l
in
l
kn
.
(333)
Algorytm:
Dla k = 1, . . . , n:
l
kk
=
$
a
kk
−
k
−1
p=1
l
2
kp
,
(334)
Dla i = k + 1, . . . , n:
l
ik
=
a
ik
−
k
−1
p=1
l
ip
l
kp
l
kk
.
(335)
Uwagi:
• Wykorzystujemy tylko górną trójkątną część macierzy A.
• Złożoność O(n
3
).
• n razy liczymy pierwiastek
• |l
kj
| √a
kk
,
k = 1, . . . , n
j = 1, . . . , k.
Metody Numeryczne II
167
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
1
1
x
y
=
0.5
1
(336)
ma rozwiązanie (
5000
9950
,
4950
9950
)
≈ (0.503, 0.497).
Stosujemy 2-cyfrową precyzję obliczeń.
Jeśli elementem głównym jest a
11
, to mamy:
0.005
1
0
−200
x
y
=
0.5
−99
,
(x, y) = (0, 0.5).
(337)
Jeśli wybierzemy a
21
jako element główny, to otrzymamy:
1
1
0
1
x
y
=
0.1
0.5
,
(x, y) = (0.5, 0.5).
(338)
Metody Numeryczne II
168
Częściowy wybór elementu głównego to nie wszystko:
Z równoważnym do (
) układem
1
200
1
1
x
y
=
100
1
(339)
są te same problemy, choć a
11
= a
21
.
Skalowanie macierzy: Przemnożenie pierwszego wiersza przez 200, to
zastąpienie macierzy A przez ˜
A = DA (i równocześnie b przez ˜
b = Db), gdzie
D =
200
0
0
1
.
(340)
Skalować można zarówno wiersze jak i kolumny, wówczas otrzymujemy układ:
D
1
AD
2
(D
−1
2
x) = D
1
AD
2
y = D
1
b,
(341)
gdzie macierze D
1
i D
2
są diagonalne.
Metody Numeryczne II
169
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:
D
2
= I,
D
1
= diag(s
1
, . . . , s
n
),
s
i
=
1
n
k=1
|a
ik
|
.
(342)
Alternatywa: modyfikacja zasady wyboru elementu głównego w kolumnie
(skalowanie przez s
i
).
Metody Numeryczne II
170
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 (
).
Wystarczy ocenić wyrażenia typu:
b
n
= f l
c
− a
1
b
1
− . . . − a
n
−1
b
n
−1
a
n
.
(343)
Analiza rozchodzenia się błędów pozwala oszacować:
c
−
n
j=1
a
j
b
j
eps
1
− eps
δ|s
n
−1
| +
n
−1
j=1
(
|s
j
| + |a
j
b
j
|)
,
(344)
gdzie
s
0
= c,
s
j
= f l(s
j
−1
− a
j
b
j
),
δ =
0
jeśli a
n
= 1
1
w przec. przyp.
(345)
Metody Numeryczne II
171
Otrzymujemy dla F = (f
ij
) = A
− L R oszacowania:
|f
ij
| = |a
ij
−
i
k=1
l
ik
r
kj
|
eps
1
− eps
i
−1
k=1
(
|a
k
ij
| + |l
ik
||r
kj
|)
dla j
i
|f
ij
| = |a
ij
−
j
k=1
l
ik
r
kj
|
eps
1
− eps
|a
j
−1
ij
| +
j
−1
k=1
(
|a
k
ij
| + |l
ik
||r
kj
|)
dla j < i,
(346)
gdzie a
k
ij
= f l(a
ij
−
k
s=1
l
is
r
sj
).
Przyjmując a
k
= max
i,j
|a
k
ij
|, a = max
0<i<n
−1
a
i
oraz uwzględniając, że r
ij
= a
i
−1
ij
, dla
i
j i zakładając, że |l
ij
| 1 szacujemy:
|f
ij
|
eps
1
− eps
(a
0
+ 2a
1
+
· · · + 2a
i
−2
+ a
i
−1
)
2(i − 1)a
eps
1
− eps
dla j
i
|f
ij
|
eps
1
− eps
(a
0
+ 2a
1
+
· · · + 2a
j
−2
+ 2a
j
−1
)
2ja
eps
1
− eps
dla j < i
(347)
Metody Numeryczne II
172
|f
ij
| 2a
eps
1
− eps
0
0
0
. . .
0
0
1
1
1
. . .
1
1
1
2
2
. . .
2
2
1
2
3
. . .
3
3
..
.
..
.
..
.
..
.
..
.
1
2
3
. . .
n
− 1 n − 1
(348)
Oszacowanie wartości a
Z definicji a
0
= max
i,j
a
i
j. Dla częściowego wyboru elementu głównego można
wykazać, że a
k
−1
2
k
a
0
. Mamy więc (zwykle pesymistyczne, ale czasem
osiągalne) oszacowanie
a
2
n
−1
a
0
.
(349)
Dla górnej macierzy Hessenberga (a
ij
= 0 dla i > j + 1) można wykazać, że
a
(n − 1)a
0
.
(350)
Dla macierzy trójdiagonalnej
a
2a
0
.
(351)
Metody Numeryczne II
173
Pełny wybór elementu głównego
• Wilkinson wykazał, że
a
k
f(k + 1)a
0
,
(352)
gdzie
f (k) =
$
k2
1
3
1
2
4
1
3
· · · k
1
k
−1
.
(353)
• W praktyce nie znaleziono kontrprzykładu dla oszacowania
a
k
(k + 1)a
0
.
(354)
• Bardziej kosztowny niż częściowy.
• Niszczy szczególną strukturę macierzy (np. trójdiagonalnych).
Metody Numeryczne II
174
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:
c
−
n
i=1
a
i
b
i
eps
1
− n · eps
n
−1
i=1
i
|a
i
b
i
| + (n + 1 + δ)|a
n
b
n
|
.
(355)
Dla układu Ly = b dostajemy:
b
r
−
r
i=1
l
ri
y
i
eps
1
− n · eps
r
i=1
i
|l
ri
y
i
| + |y
r
|
,
(356)
Metody Numeryczne II
175
czyli
|b − Ly|
eps
1
− n · eps
(
|L|D − I)|y|,
D =
1
0
2
. ..
0
n
.
(357)
A zatem y jest rozwiązaniem dokładnym nieco zaburzonego układu:
(L + ∆L)y = b,
|∆L|
eps
1
− n · eps
(
|L|D − I).
(358)
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
176
Metoda ortogonalizacji Householdera
Idea: Kolejne transformacje macierzy dobierane tak, by wskaźnik uwarunkowania
zadania nadmiernie nie wzrósł.
Analiza prowadzi do wniosku: przekształcenia o macierzach unitarnych zachowują
wskaźnik uwarunkowania.
Propozycja Householdera: unitarna macierz transformacji P spełniająca
P
− I − 2ww
H
,
gdzie w
H
w = 1, w
∈ C
n
.
(359)
Macierz P jest hermitowska:
P
H
= I
− (2ww
H
)
H
= I
− 2ww
H
= P,
(360)
i unitarna:
P
H
P = PP = P
2
= (I
− 2ww
H
)(I
− 2ww
H
)
= I
− 2ww
H
− 2ww
H
+ 4ww
H
ww
H
= I.
(361)
Metody Numeryczne II
177
Wektor w dobierany jest tak, aby pewien wektor x = [x
i
, . . . , x
n
]
T
był
przekształcony na kierunek wektora jednostkowego e
1
:
Px = ke
1
.
(362)
Dokładniejsza analiza prowadzi do:
P = I
−
uu
H
γ
,
(363)
gdzie
σ =
n
i=1
|x
i
|
2
,
x
1
= e
iα
|x
1
|, k = −σe
iα
,
u = x
− ke
1
=
e
iα
(
|x
1
| + σ)
x
2
..
.
x
n
,
γ = σ(σ +
|x
1
|).
(364)
Metody Numeryczne II
178
Algorytm:
Przekształcamy macierz A krok po kroku zerując kolejne kolumny pod przekątną:
A = A
(0)
→ A
(1)
→ · · · → A
(n
−1)
= R =
r
11
. . .
r
1n
. ..
..
.
0
r
nn
.
(365)
W j=tym kroku przekształcamy macierz
A
(j
−1)
=
r
11
. . .
x
1,j
−1
a
(j
−1)
1j
. . .
a
(j
−1)
1n
. ..
..
.
..
.
..
.
0
r
j
−1,j−1
a
(j
−1)
j
−1,j
. . .
a
(j
−1)
j
−1,n
a
(j
−1)
jj
. . .
a
(j
−1)
jn
0
..
.
..
.
a
(j
−1)
nj
. . .
a
(j
−1)
nn
=
D
B
0
˜
A
(j
−1)
. (366)
Metody Numeryczne II
179
Szukamy takiej unitarnej macierzy przekształcenia P
j
, aby
P
j
=
I
j
−1
0
0
˜
P
j
,
˜
P
j
a
(j
−1)
jj
..
.
a
(j
−1)
nj
= ke
1
∈ C
n
−j+1
.
(367)
Uwagi:
• Mnożenie przez macierz
˜
P
j
= I
−
u
j
u
H
j
γ
j
(368)
realizuje się w następujący sposób:
˜
P
j
˜
A
(j
−1)
= ˜
A
(j
−1)
− u
j
y
H
j
,
gdzie y
H
j
=
u
H
j
˜
A
(j
−1)
γ
j
.
(369)
• Algorytm wymaga około
2n
3
3
działań.
• 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
= P
1
· · · P
n
−1
są ortonormalne
Metody Numeryczne II
180
Metoda ortogonalizacji Grama-Schmidta
Jeśli macierz A ma rozkład:
A = QR,
(370)
gdzie macierz Q ma ortonormalne kolumny, to można je wyznaczyć metodą
ortogonalizacji Grama-Schmidta.
Dla k-tej kolumny macierzy A = (a
1
, . . . , a
n
) mamy:
a
k
=
k
i=0
r
ik
q
i
,
k = 1, . . . , n.
(371)
Kolumna a
k
jest kombinacją liniową wektorów q
1
, . . . , q
k
i odwrotnie: wektor q
k
jest kombinacją liniową a
1
, . . . , a
k
.
Kolejne q
k
wyznaczamy rekurencyjnie tak, by zachowywać ortonormalność.
Metody Numeryczne II
181
Algorytm:
• Wyznaczamy q
1
i r
11
:
r
11
← ||a
1
||,
q
1
←
a
1
r
11
.
(372)
• Znamy q
1
, . . . , q
k
−1
oraz r
ij
dla j
k − 1.
– Wyznaczamy wektor
b
k
← a
k
− r
1k
q
1
− . . . − r
k
−1,k
q
k
−1
(373)
tak, aby był ortogonalny do q
1
, . . . , q
k
−1
:
q
H
i
b
k
= 0
⇒ r
ik
← q
H
i
a
k
,
i = 1, . . . , k
− 1
(374)
– Wyznaczamy q
k
i r
kk
:
r
kk
← ||b
k
||,
q
k
←
b
k
r
kk
.
(375)
Metody Numeryczne II
182
Uwagi:
• Dla wszystkich 1 i, j k mamy:
q
H
i
q
j
=
1
dla i = j,
0
w przec. przyp.
(376)
• Algorytm jest niestabilny numerycznie gdy wektory a
i
są bliskie liniowej
zależności – błędy numeryczne sprawiają, że b
k
nie jest ortogonalny do q
i
.
Reortogonalizacja wektora b
k
Obliczamy poprawki w postaci skalarów
∆r
ik
← q
H
i
b
k
,
i = 1, . . . , k
− 1,
(377)
oraz wektor
˜
b
k
= b
k
− ∆r
1k
q
1
− . . . − ∆r
k
−1,k
q
k
−1
.
(378)
Poprawiamy:
r
kk
← ||˜b
k
||,
q
k
←
˜
b
k
r
kk
,
r
ik
← r
ik
+ ∆r
ik
.
(379)
Metody Numeryczne II
183
Wyznaczanie wartości i wektorów własnych
Definicja 14 Wektor x jest wektorem własnym macierzy A jeśli istnieje
taka liczba λ, że Ax = λx.
Twierdzenie 31 Liczba λ jest wartością własną macierzy A wtedy i tylko
wtedy, gdy jest pierwiastkiem wielomianu charakterystycznego
det(A
− λI) macierzy A.
Definicja 15 Macierze A i B są podobne jeśli istnieje nieosobliwa macierz
podobieństwa P tj. taka, że P
−1
AP = B.
Twierdzenie 32 Macierze podobne mają takie same wartości własne.
Twierdzenie 33 Wartości i wektory własne macierzy symetrycznej są
rzeczywiste.
Twierdzenie 34 Wartości własne macierzy A + cI to wartości własne
macierzy A przesunięte o c.
Metody Numeryczne II
184
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
185
Metoda potęgowa
Dla wyznaczania pojedynczych wartości i wektorów własnych.
Algorytm:
• i ← 0, wybieramy wektor x
0
taki, że
||x
0
||
∞
= 1
• Tak długo jak i jest mniejsze niż założona liczba iteracji:
– v
i+1
← Ax
i
,
m
i+1
← ||v
i+1
||
∞
– Jeżeli m
i+1
= 0 to KONIEC!
– x
i+1
←
v
i+1
m
i+1
,
i
← i + 1
Uwagi:
• Jeśli x
2i
i
→∞
−→ x, to m
i
i
→∞
−→ m.
• Jeśli ||x
i+1
− x
i
||
∞
i
→∞
−→ 0, to Ax = mx.
• Jeśli ||x
i+1
− x
i
||
∞
i
→∞
−→ 2, to Ax = −mx.
• W ogólnym przypadku trudno ocenić kiedy ciąg (x
2i
) jest zbieżny.
• Można wykazać zbieżność np. gdy macierz jest podobna do diagonalnej.
Metody Numeryczne II
186
Metoda QR
Dla macierzy Hessenberga H (h
ij
= 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.
• Wszystkie H
(i)
są podobne, bo H
(i+1)
= Q
(i)T
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
187
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: p =
1
2
(λ
i
+ λ
n
) dla λ
1
oraz p =
1
2
(λ
1
+ λ
j
) dla λ
n
.
– Metoda QR: p równe otrzymanym w poprzednich iteracjach wartościom
własnym.
Sprowadzanie macierzy do postaci Hessenberga
• eliminacja Gaussa
• metoda Householdera
• metoda Givensa
Metody Numeryczne II
188
Szukanie wektorów własnych
• Jeżeli P
−1
AP = B, oraz x jest wektorem własnym macierzy B, to y = Bx
jest wektorem własnym macierzy A.
• Jeśli macierz B jest trójkątna górna, to wartości własnej λ
i
= b
ii
odpowiada
wektor własny x taki, że:
x
(i)
j
= 0,
j = n, n
− 1, . . . , i + 1
x
(i)
i
= 1
x
(i)
j
=
−
i
k=j+1
b
jk
x
(i)
k
b
jj
− b
ii
,
j = i
− 1, . . . , 1
(380)
– Jeśli B ma wielokrotną wartość własną λ, to wektor własny można
wyliczyć tylko dla najmniejszego i takiego, że b
ii
= λ.
Metody Numeryczne II
189
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
−1
AP = H.
(381)
– 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)
.
(382)
– Po uznaniu w H
(N )
zbieżności do macierzy trójkątnej górnej wyznaczamy
jej wektory własne zgodnie z (
).
– Wyznaczamy wektory własne macierzy A:
y
(k)
= Z
(N )
x
(k)
.
(383)
Metody Numeryczne II
190
Metoda LR
• Algorytm: Dane: macierz A.
– Obliczamy ciąg macierzy A
(i)
dla i = 1, 2, . . .:
A
(1)
← A
A
(i)
= L
(i)
U
(i)
,
A
(i+1)
← U
(i)
L
(i)
(384)
do uzyskania zbieżności do macierzy trójkątnej.
Jednocześnie wyznaczamy:
Z
(1)
= I,
Z
(i+1)
= Z
(i)
L
(i)
.
(385)
– Wyznaczamy wektory własne macierzy A
(N )
zgodnie z (
).
– Wyznaczamy wektory własne macierzy A:
y
(k)
= Z
(N )
x
(k)
.
(386)
• 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(n
3
)
→ O(n
2
)).
Metody Numeryczne II
191
Metoda iteracji odwrotnej
• Przyspieszenie znieżności metody potęgowej
• Algorytm:
– i
← 0, wybieramy wektor x
0
taki, że
||x
0
||
∞
= 1
– Tak długo jak i jest mniejsze niż założona liczba iteracji:
∗ v
i+1
← (A − k
i
I)
−1
x
i
,
m
i+1
← ||v
i+1
||
∞
∗ x
i+1
←
v
i+1
m
i+1
,
i
← i + 1
• Jeśli A ma wartości własne λ
1
. . . λ
n
, to wartościami własnymi A
− k
i
I
są:
1
λ
1
− k
i
,
1
λ
2
− k
i
, . . . ,
1
λ
n
− k
i
.
(387)
• Ciąg m
1
, m
2
, . . . pozwala wyznaczyć wartość własną λ
n
jak w metodzie
potęgowej.
• Jeśli k
i
jest bliskie λ
n
, to macierz A
− k
i
I
jest źle uwarunkowana, co można
zlekceważyć gdy x jest dobrze uwarunkowany i norma
||v
i+1
||
∞
jest duża.
Wyszukiwarka
Podobne podstrony:
EDI wyk
Wyk ad 5 6(1)
zaaw wyk ad5a 11 12
Wyk 02 Pneumatyczne elementy
Automatyka (wyk 3i4) Przel zawory reg
Wyk ECiUL#1 2013
wyk II
Wyk 07 Osprz t Koparki
budownictwo stany skupenia wyk 3
6 wykad WiĄzania chemiczne[F]
Wyk ECiUL#9S 2013
Wyk ad II
zaaw wyk ad6
więcej podobnych podstron