Podstawy metod numerycznych
Wykład 5
dr K. A. Smoliński
Wyższa Szkoła Informatyki
rok akademicki 2007/2008
dr K. A. Smoliński (WSInf)
2007/2008
1 / 26
Podstawy metod numerycznych
Wykład 5
dr K. A. Smoliński
Wyższa Szkoła Informatyki
rok akademicki 2007/2008
dr K. A. Smoliński (WSInf)
2007/2008
1 / 26
Różniczkowanie i całkowanie numeryczne
1
2
Interpolacja w całkowaniu numerycznym
Zastosowanie interpolacji wielomianowej
Wzór trapezów
Wzór Simpsona
Całki z funkcją wagową
3
Wielomiany ortonormalne
Kwadratury Gaussa
Zmiana przedziału całkowania
Zbieżność i błąd
4
dr K. A. Smoliński (WSInf)
2007/2008
2 / 26
Różniczkowanie numeryczne
Różniczkowanie numeryczne
Potrafimy obliczyć wartości funkcji f w zadanych punktach. Czy można
stąd uzyskać przybliżenie pochodnej f
0
(c) lub całki
R
b
a
f (x) dx?
Same wartości f (x
0
), f (x
1
), . . . , f (x
n
) nie mówią wiele o funkcji, nawet
jeżeli wiemy, że jest ona rzeczywista i ciągła. Jednoznacznie wyznaczają ją
tylko wtedy, gdy wiemy o niej, że jest wielomianem stopnia nie wyższego
niż n — możemy obliczyć f
0
(c) i
R
b
a
f (x) dx dokładnie. W praktyce
sytuacja wygląda tak, że informacje o funkcji f nie określają jej całkowicie,
ale poozwalają oszacować błąd obliczonych przybliżeń wartości
pochodnych lub całki.
dr K. A. Smoliński (WSInf)
2007/2008
3 / 26
Różniczkowanie numeryczne
Z definicji pochodnej wynika wzór różniczkowania numerycznego
f
0
(x) ≈
1
h
[f (x + h) − f (x)]
.
Dla funkcji liniowej jest on dokładny dla każdego h 6= 0, dla innych funkcji
jest tak wyjątkowo.
Błąd przybliżenia możemy oszacować ze wzoru Taylora: Jeżeli f
0
jest
ciągła na [x, x + h], a f
00
istnieje na (x, x + h), to
f (x + h) = f (x) + hf
0
(x) +
1
2
h
2
f
00
(ξ) ,
skąd
f
0
(x) =
1
h
[f (x + h) − f (x)] −
1
2
hf
00
(ξ)
.
Składnik −
1
2
hf
00
(ξ) nazywamy
błędem obcięcia
.
dr K. A. Smoliński (WSInf)
2007/2008
4 / 26
Różniczkowanie numeryczne
Przykład
Znajdź przybliżoną wartość pochodnej funkcji f (x) = cos x w punkcie
π
4
,
przyjmując h = 0.01.
Przybliżeniem szukanej wartości jest
0.7000005 − 0.7071068
0.01
=
−0.7106301
,
a jej błąd wynosi
1
2
hf
00
(ξ)
= 0.005| cos ξ| ≤ 0.005| cos
π
4
| ≈
3.5355338
10
−3
.
dr K. A. Smoliński (WSInf)
2007/2008
5 / 26
Różniczkowanie numeryczne
Dokładność przybliżeń pochodnej oceniamy według wykładnika p
w czynniku h
p
błędu, im większe p, tym lepiej. Poprzedni wzór ma p = 1
i jest dość słaby. Lepszy jest wzór
f
0
(x) ≈
1
2h
[f (x + h) − f (x − h)]
.
Wzór powyższy jest rzędu 2. Istotnie, jeżeli f
000
jest ciągła na
(x − h, x + h), to ze wzoru Taylora
f (x ± h) = f (x) ± hf
0
(x) +
1
2
h
2
f
00
(x) ±
1
3!
f
000
(ξ
±
) ,
mamy
f
0
(x) =
1
2h
[f (x + h) − f (x − h)] −
1
12
h
2
f
000
(ξ
−
) + f
000
(ξ
+
)
,
co, po zastosowaniu tw. o wartości średniej, daje
f
0
(x) =
1
2h
[f (x + h) − f (x − h)] −
1
6
h
2
f
000
(ξ)
.
dr K. A. Smoliński (WSInf)
2007/2008
6 / 26
Różniczkowanie numeryczne
Przykład
Znajdź przybliżoną wartość pochodnej funkcji f (x) = cos x w punkcie
π
4
,
przyjmując h = 0.01.
Przybliżeniem szukanej wartości jest
0.7000005 − 0.7141424
2 · 0.01
=
−0.7070945
,
a jej błąd wynosi
1
6
h
2
f
000
(ξ)
=
0.0001
6
| sin(ξ)| ≤
0.0001
6
| sin(
π
4
+ h)| ≈
1.1902373
10
−5
.
dr K. A. Smoliński (WSInf)
2007/2008
7 / 26
Interpolacja w całkowaniu numerycznym
Interpolacja w całkowaniu numerycznym
Całki nieoznaczone z wielu funkcji (w tym bardzo prostych) nie wyrażają
się przez funkcje elementarne. Nie możemy więc znaleźć dokładnych
wartości całek oznaczonych. Przybliżone wartości całki oznaczonej
R
b
a
f (x) dx znajdujemy zastępując funkcję f inną, bliską funkcją g, dla
której całkę potrafimy policzyć — stosujemy wzór przybliżony:
Z
b
a
f (x) dx ≈
Z
b
a
g(x) dx .
Dobrą funkcją g może być wielomian (łatwo obliczyć całkę) interpolujący
funkcję f w danych węzłach, co zazwyczaj zapewnia bliskość obu funkcji.
dr K. A. Smoliński (WSInf)
2007/2008
8 / 26
Interpolacja w całkowaniu numerycznym
Zastosowanie interpolacji wielomianowej
Zastosowanie interpolacji wielomianowej
W przedziale całkowania [a, b] wybieramy węzły x
0
, x
1
, . . . , x
n
i stosujemy
wzór interpolacyjny Lagrange’a:
p(x) =
n
X
i=0
f (x
i
)l
i
(x) ,
gdzie
l
i
(x) =
n
Y
j=0
j6=i
x − x
j
x
i
− x
j
.
Stąd
Z
b
a
f (x) dx ≈
Z
b
a
p(x) dx =
n
X
i=0
f (x
i
)
Z
b
a
l
i
(x) dx .
Jest to przykład typowego wzoru całkowania numerycznego, zwanego
kwadraturą
:
Z
b
a
f (x) dx ≈
n
X
i=0
A
i
f (x
i
)
.
W tym przypadku A
i
=
R
b
a
l
i
(x) dx (0 ≤ i ≤ n). Jeżeli przy tym węzły
interpolacji są takie, że
x
i
= a + (b − a)
i
n
(0 ≤ i ≤ n), to kwadraturę tę
nazywamy
wzorem Newtona–Cotesa
.
dr K. A. Smoliński (WSInf)
2007/2008
9 / 26
Interpolacja w całkowaniu numerycznym
Wzór trapezów
Wzór trapezów
Najprostszy z wzorów Newtona–Cotesa otrzymujemy dla n = 1; wówczas
x
0
= a i x
i
= b oraz
l
0
(x) =
b − x
b − a
,
l
1
(x) =
x − a
b − a
,
więc
A
0
=
Z
b
a
l
0
(x) dx =
b − a
2
,
A
1
=
Z
b
a
l
1
(x) dx =
b − a
2
.
Ostatecznie otrzymujemy tzw.
wzór trapezów
Z
b
a
f (x) dx ≈
b − a
2
[f (a) + f (b)]
.
Wzór ten jest dokładny jeżeli f jest funkcją liniową, w pozostałych
przypadkach jego błąd wynosi
−
(b − a)
3
12
f
00
(ξ)
,
gdzie
ξ ∈ (a, b) .
dr K. A. Smoliński (WSInf)
2007/2008
10 / 26
Interpolacja w całkowaniu numerycznym
Wzór trapezów
Jeżeli podzielimy przedział [a, b] na podprzedziały punktami
a = x
0
< x
1
< · · · < x
n−1
< x
n
= b ,
a następnie zastosujemy do każdego z nich wzór trapezów, to
otrzymujemy
złożony wzór trapezów
:
Z
b
a
f (x) dx
=
n
X
i=1
Z
x
i
x
i−1
f (x) dx
≈
n
X
i=1
x
i
− x
i−1
2
[f (x
i−1
) + f (x
i
)]
.
Wzór jest dokładny, jeżeli wykres funkcji jest linią łamaną, której
wierzchołki mają odcięte x
i
.
Złożony wzór trapezów upraszcza się, jeżeli wszystkie podprzedziały są
jednakowo długie, tj. dla
x
i
= a + ih
, gdzie
h =
b−a
n
:
Z
b
a
f (x) dx ≈
h
2
"
f (x) + 2
n−1
X
i=1
f (a + ih) + f (b)
#
.
Jeżeli f
00
jest ciągła na (a, b), to błąd tego wzoru wynosi
−
(b − a)
3
12n
2
f
00
(ξ)
,
gdzie
ξ ∈ (a, b) .
dr K. A. Smoliński (WSInf)
2007/2008
11 / 26
Interpolacja w całkowaniu numerycznym
Wzór Simpsona
Wzór Simpsona
Wzór Newtona–Cotesa dla n = 2 na przedziale [a, b] nosi nazwę
wzoru
Simpsona
:
Z
b
a
f (x) dx ≈
b − a
6
h
f (a) + 4f (
a+b
2
) + f (b)
i
.
Wzór ten jest dokładny nie tylko dla funkcji kwadratowej, ale także dla
funkcji sześciennej! Błąd jest dany przez
−
(b − a)
5
2880
f
(4)
ξ
,
gdzie
ξ ∈ (a, b) .
dr K. A. Smoliński (WSInf)
2007/2008
12 / 26
Interpolacja w całkowaniu numerycznym
Wzór Simpsona
Złożony wzór Simpsona
otrzymujemy stosując ostatni wzór do n
podprzedziałów jednakowej długości. Oznaczając
x
i
= a + ih
(0 ≤ i ≤ 2n)
i
h =
b−a
2n
, mamy
Z
b
a
f (x) dx ≈
h
3
"
f (x
0
) + 4
n
X
i=1
f (x
2i−1
) + 2
n−1
X
i=1
f (x
2i
) + f (x
2n
)
#
.
Błąd tego wzoru wynosi
−
(b − a)
5
2880n
4
f
(4)
(ξ)
,
gdzie
ξ ∈ (a, b) .
dr K. A. Smoliński (WSInf)
2007/2008
13 / 26
Interpolacja w całkowaniu numerycznym
Całki z funkcją wagową
Całki z funkcją wagową
Procedurą prowadzącą do wzorów Newtona–Cotesa można również
otrzymać ogólniejsze wzory przybliżone postaci
Z
b
a
f (x)w(x) dx ≈
n
X
i=0
A
i
f (x
i
) ,
gdzie w jest nieujemną
funkcją wagową
(
wagą
). W tym przypadku
A
i
=
Z
b
a
l
i
(x)w(x) dx .
dr K. A. Smoliński (WSInf)
2007/2008
14 / 26
Kwadratury Gaussa
Kwadratury Gaussa
Wiemy jak tworzyć kwadratury
Z
b
a
f (x)w(x) dx ≈
n
X
i=0
A
i
f (x
i
)
(w jest ustaloną nieujemną funkcją wagową), dokładne dla każdego
wielomianu f ∈ Π
n
. We wzorach tych układ węzłów x
i
jest dowolny
i określa jednoznacznie współczynniki A
i
.
Powstaje pytanie, czy pewne układy węzłów są lepsze od innych? Np.
wygodny byłby wzór całkowania taki, że wszystkie współczynniki A
i
byłyby równe:
Z
b
a
f (x)w(x) dx ≈ A
n
X
i=0
f (x
i
) .
Gdy funkcja wagowa jest stała (w(x) = 1), wzór taki, zwany
kwadraturą
Czebyszewa
, istnieje tylko dla n = 0, 1, . . . , 6 i n = 8.
dr K. A. Smoliński (WSInf)
2007/2008
15 / 26
Kwadratury Gaussa
Wielomiany ortonormalne
Wielomiany ortonormalne
Dla dwóch funkcji f i g określonych na [a, b] określamy ich
iloczyn
skalarny
jako
hf , gi =
Z
b
a
f (x)g(x)w(x) dx
.
Układ wielomianów p
n
(n = 0, 1, . . .) jest
ortogonalny na przedziale [a, b]
z wagą w
, jeżeli stopień każdego z nich jest równy n oraz
hp
m
, p
n
i = 0
dla
m 6= n
.
Każdy wielomian jest określony z dokładnością do stałego czynnika.
Czynniki te wybieramy zazwyczaj tak, aby hp
n
, p
n
i = 1, albo tak aby p
n
był wielomianem standardowym, tj. współczynnik w p
n
(x) przy x
n
jest
równy 1.
dr K. A. Smoliński (WSInf)
2007/2008
16 / 26
Kwadratury Gaussa
Wielomiany ortonormalne
Twierdzenie
Ciąg {p
n
} wielomianów ortogonalnych standardowych w przedziale [a, b]
z wagą w można wyznaczyć ze wzorów:
p
0
(x) = 1 ,
p
1
(x) = x − a
1
,
p
n
(x) = (x − a
n
)p
n−1
(x) − b
n
p
n−2
(x)
(n ≥ 2) ,
gdzie
a
n
=
hxp
n−1
, p
n−1
i
hp
n−1
, p
n−1
i
,
b
n
=
hxp
n−1
, p
n−2
i
hp
n−2
, p
n−2
i
.
dr K. A. Smoliński (WSInf)
2007/2008
17 / 26
Kwadratury Gaussa
Wielomiany ortonormalne
Wielomiany Legendre’a
Wielomiany Legendre’a
są ortogonalne w przedziale [−1, 1] z wagą 1.
P
0
(x) = 1 ,
P
1
(x) = x ,
P
2
(x) = x
2
−
1
3
,
P
3
(x) = x
3
−
3
5
x ,
P
4
(x) = x
4
−
6
7
x
2
+
3
35
,
P
5
(x) = x
5
−
10
9
x
3
+
5
21
x .
dr K. A. Smoliński (WSInf)
2007/2008
18 / 26
Kwadratury Gaussa
Wielomiany ortonormalne
Wielomiany Czebyszewa
Wielomiany Czebyszewa
są ortogonalne w przedziale [−1, 1] z wagą
1
√
1−x
2
.
T
0
(x) = 1 ,
T
1
(x) = x ,
T
2
(x) = 2x
2
− 1 ,
T
3
(x) = 4x
3
− 3x ,
T
4
(x) = 8x
4
− 8x
2
+ 1 ,
T
5
(x) = 16x
5
− 20x
3
+ 5x ,
T
6
(x) = 32x
6
− 48x
4
+ 18x
2
− 1 .
Dla n ≥ 2 można je określić wzorem rekurencyjnym
T
n
(x) = 2xT
n−1
(x) − T
n−2
(x) .
Twierdzenie
W przedziale [−1, 1] wielomiany Czebyszewa wyrażają się wzorem
T
n
(x) = cos(n arccos x) .
dr K. A. Smoliński (WSInf)
2007/2008
19 / 26
Kwadratury Gaussa
Kwadratury Gaussa
Kwadratury Gaussa
Twierdzenie
Jeżeli węzły x
0
, x
1
, . . . , x
n
są zerami (n + 1)-szego wielomianu
ortogonalnego p
n+1
na przedziale [a, b] z wagą w, to
kwadratura Gaussa
Z
b
a
f (x)w(x) dx ≈
n
X
i=0
A
i
f (x
i
)
o współczynnikach
A
i
=
Z
b
a
w(x)
n
Y
j=0
j6=i
x − x
j
x
i
− x
j
dx
jest dokładna dla każdej funkcji f ∈ Π
2n+1
.
dr K. A. Smoliński (WSInf)
2007/2008
20 / 26
Kwadratury Gaussa
Kwadratury Gaussa
Twierdzenie
Jeżeli funkcja f ∈ C [a, b] jest ortogonalna w tym przedziale z wagą w
względem wszystkich wielomianów klasy Π
n
, to w (a, b) zmienia znak co
najmniej n + 1 razy.
dr K. A. Smoliński (WSInf)
2007/2008
21 / 26
Kwadratury Gaussa
Zmiana przedziału całkowania
Zmiana przedziału całkowania
Wzór całkowania w pewnym przedziale można dostosować do innego
przedziału przez liniowe przekształcenie zmiennej przy zachowaniu jego
dokładności dla wielomianów pewnego stopnia. Niech znany wzór
całkowania ma postać
Z
d
c
f (t) dt ≈
n
X
i=0
A
i
f (t
i
) .
Aby otrzymać wzór dla przedziału [a, b] dokonujemy podstawienia
x = λ(t) =
(b − a)t + ad − bc
d − c
,
co daje
Z
b
a
f (x) dx
=
b − a
d − c
Z
d
c
f (λ(t)) dt
≈
b − a
d − c
n
X
i=0
A
i
f (
(b−a)t
i
+ad−bc
d−c
)
.
dr K. A. Smoliński (WSInf)
2007/2008
22 / 26
Kwadratury Gaussa
Zbieżność i błąd
Zbieżność i błąd
Twierdzenie (Stieltjes)
Jeżeli f ∈ C [a, b], to
lim
n→∞
n
X
i=0
A
ni
f (x
ni
) =
Z
b
a
f (x)w(x) dx .
Twierdzenie
Jeżeli f ∈ C
2n
[a, b], to dla kwadratury Gaussa z n węzłami zachodzi
Z
b
a
f (x)w(x) dx =
n−1
X
i=0
A
i
f (x
i
) +
f
(2n)
(ξ)
(2n)!
Z
b
a
q
2
(x)w(x) dx ,
gdzie ξ ∈ (a, b) i q(x) =
Q
n−1
i=0
(x − x
i
).
dr K. A. Smoliński (WSInf)
2007/2008
23 / 26
Metoda Romberga
Metoda Romberga
Dla
metody Romberga
punktem wyjścia jest złożony wzór trapezów,
stosowany do przedziału [a, b] podzielonego na 2
n
podprzedziałów równej
długości (n ≥ 0). Przyjmujemy:
R(n, 0) = h
n
"
f (a) +
2
n
−1
X
i=1
f (a + ih
n
) + f (b)
#
gdzie h
n
=
b − a
2
n
.
W najprostszym przypadku, dla n = 0 mamy prosty wzór trapezów:
R(0, 0) =
1
2
[f (a) + f (b)] .
R(1, 0), R(2, 0), . . . obliczamy rekurencyjnie, aby uniknąć obliczania
wartości funkcji f w tych samych punktach:
R(n, 0) =
1
2
R(n − 1, 0) + h
n
2
n−1
X
i=1
f (a + (2i − 1)h
n
) .
dr K. A. Smoliński (WSInf)
2007/2008
24 / 26
Metoda Romberga
Następne serie przybliżeń całki (dla m > 0) obliczamy ze wzorów:
R(n, m) =
4
m
R(n, m − 1) − R(n − 1, m − 1)
4
m
− 1
.
metoda Romberga
Require: a, b, f (x), M
Ensure: R(n, m) ≈
R
b
a
f (x) dx
h ← b − a
R(0, 0) ← (b − a)(f (a) − f (b))/2
for n ← 1, 2, . . . , M do
h ← h/2
R(n, 0) ← R(n − 1, 0)/2 + h
P
2
n−1
i=1
f (a + (2i − 1)h)
for m ← 1, 2, . . . , n do
R(n, m) ← (4
m
R(n, m − 1) − R(n − 1, m − 1))/(4
m
− 1)
end for
end for
dr K. A. Smoliński (WSInf)
2007/2008
25 / 26
Metoda Romberga
Zbieżność metody Romberga
Zbieżność metody Romberga
Zbieżność metody Romberga jest szczególnie dobra, gdy f ma
dostatecznie wiele pochodnych. Jeżeli f ∈ C
2m
[a, b], to R(n, m) przybliża
całkę z błędem rzędu O(h
2m
), gdzie h = (b − a)/2
n
. Metoda jest także
sensowna dla funkcji, które są jedynie ciągłe:
Twierdzenie
Jeżeli f ∈ C [a, b], to
lim
n→∞
R(n, m) =
Z
b
a
f (x) dx .
dr K. A. Smoliński (WSInf)
2007/2008
26 / 26