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) Podstawy metod numerycznych 2007/2008 1 / 26
Różniczkowanie i całkowanie numeryczne
1
Różniczkowanie numeryczne
2
Interpolacja w całkowaniu numerycznym
Zastosowanie interpolacji wielomianowej
Wzór trapezów
Wzór Simpsona
Całki z funkcją wagową
3
Kwadratury Gaussa
Wielomiany ortonormalne
Kwadratury Gaussa
Zmiana przedziału całkowania
Zbieżność i błąd
4
Metoda Romberga
Zbieżność metody Romberga
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 2007/2008 2 / 26
Różniczkowanie numeryczne
Różniczkowanie numeryczne
Potrafimy obliczyć wartości funkcji f w zadanych punktach. Czy można

b

stąd uzyskać przybliżenie pochodnej f (c) lub całki f (x) dx?
a
Same wartości f (x0), f (x1), . . . , f (xn) 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

b

niż n  możemy obliczyć f (c) i f (x) dx dokładnie. W praktyce
a
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) Podstawy metod numerycznych 2007/2008 3 / 26
Różniczkowanie numeryczne
Z definicji pochodnej wynika wzór różniczkowania numerycznego
1

f (x) H" [f (x + h) - f (x)] .
h
Dla funkcji liniowej jest on dokładny dla każdego h = 0, dla innych funkcji

jest tak wyjątkowo.

Błąd przybliżenia możemy oszacować ze wzoru Taylora: Jeżeli f jest

ciągła na [x, x + h], a f istnieje na (x, x + h), to
1

f (x + h) = f (x) + hf (x) + h2f (�) ,
2
skąd
1 1

f (x) = [f (x + h) - f (x)] - hf (�) .
h 2

Składnik -1hf (�) nazywamy błędem obcięcia.
2
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 2007/2008 4 / 26
Różniczkowanie numeryczne
Przykład
Ą
Znajdz
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.7106301 ,
0.01
a jej błąd wynosi

1 Ą

hf (�) = 0.005| cos �| d" 0.005| cos | H" 3.535533810-3 .
2 4
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 2007/2008 5 / 26
Różniczkowanie numeryczne
Dokładność przybliżeń pochodnej oceniamy według wykładnika p
w czynniku hp błędu, im większe p, tym lepiej. Poprzedni wzór ma p = 1
i jest dość słaby. Lepszy jest wzór
1

f (x) H" [f (x + h) - f (x - h)] .
2h

Wzór powyższy jest rzędu 2. Istotnie, jeżeli f jest ciągła na
(x - h, x + h), to ze wzoru Taylora
1 1

f (x ą h) = f (x) ą hf (x) + h2f (x) ą f (�ą) ,
2 3!
mamy

1 1

f (x) = [f (x + h) - f (x - h)] - h2 f (�-) + f (�+) ,
2h 12
co, po zastosowaniu tw. o wartości średniej, daje
1 1

f (x) = [f (x + h) - f (x - h)] - h2f (�) .
2h 6
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 2007/2008 6 / 26
Różniczkowanie numeryczne
Przykład
Ą
Znajdz
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
= -0.7070945 ,
2 � 0.01
a jej błąd wynosi

1 0.0001 0.0001 Ą

h2f (�) = | sin(�)| d" | sin( + h)| H" 1.190237310-5 .
6 6 6 4
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 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 znalez
ć dokładnych
wartości całek oznaczonych. Przybliżone wartości całki oznaczonej

b
f (x) dx znajdujemy zastępując funkcję f inną, bliską funkcją g, dla
a
której całkę potrafimy policzyć  stosujemy wzór przybliżony:

b b
f (x) dx H" g(x) dx .
a a
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) Podstawy metod numerycznych 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 x0, x1, . . . , xn i stosujemy
wzór interpolacyjny Lagrange a:
n n

x - xj
p(x) = f (xi)li(x) , gdzie li(x) = .
xi - xj
i=0 j=0
j=i

Stąd

n
b b b

f (x) dx H" p(x) dx = f (xi) li(x) dx .
a a a
i=0
Jest to przykład typowego wzoru całkowania numerycznego, zwanego
kwadraturą:

n
b

f (x) dx H" Aif (xi) .
a
i=0

b
W tym przypadku Ai = li(x) dx (0 d" i d" n). Jeżeli przy tym węzły
a
i
interpolacji są takie, że xi = a + (b - a)n (0 d" i d" n), to kwadraturę tę
nazywamy wzorem Newtona Cotesa.
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 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
x0 = a i xi = b oraz
b - x x - a
l0(x) = , l1(x) = ,
b - a b - a
więc

b b
b - a b - a
A0 = l0(x) dx = , A1 = l1(x) dx = .
2 2
a a
Ostatecznie otrzymujemy tzw. wzór trapezów

b
b - a
f (x) dx H" [f (a) + f (b)] .
2
a
Wzór ten jest dokładny jeżeli f jest funkcją liniową, w pozostałych
przypadkach jego błąd wynosi
(b - a)3

- f (�) , gdzie � " (a, b) .
12
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 2007/2008 10 / 26
Interpolacja w całkowaniu numerycznym Wzór trapezów
Jeżeli podzielimy przedział [a, b] na podprzedziały punktami
a = x0 < x1 < � � � < xn-1 < xn = b ,
a następnie zastosujemy do każdego z nich wzór trapezów, to
otrzymujemy złożony wzór trapezów:

n n
b xi

xi - xi-1
f (x) dx = f (x) dx H" [f (xi-1) + f (xi)] .
2
a xi-1
i=1 i=1
Wzór jest dokładny, jeżeli wykres funkcji jest linią łamaną, której
wierzchołki mają odcięte xi.
Złożony wzór trapezów upraszcza się, jeżeli wszystkie podprzedziały są
b-a
jednakowo długie, tj. dla xi = a + ih, gdzie h = :
n


n-1
b

h
f (x) dx H" f (x) + 2 f (a + ih) + f (b) .
2
a
i=1

Jeżeli f jest ciągła na (a, b), to błąd tego wzoru wynosi
(b - a)3

- f (�) , gdzie � " (a, b) .
12n2
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 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:


b
b - a
f (x) dx H" f (a) + 4f (a+b ) + f (b) .
2
6
a
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
(4)
- f � , gdzie � " (a, b) .
2880
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 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 xi = a + ih (0 d" i d" 2n)
b-a
i h = , mamy
2n


n n-1
b

h
f (x) dx H" f (x0) + 4 f (x2i-1) + 2 f (x2i) + f (x2n) .
3
a
i=1 i=1
Błąd tego wzoru wynosi
(b - a)5
(4)
- f (�) , gdzie � " (a, b) .
2880n4
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 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

n
b

f (x)w(x) dx H" Aif (xi) ,
a
i=0
gdzie w jest nieujemną funkcją wagową (wagą). W tym przypadku

b
Ai = li(x)w(x) dx .
a
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 2007/2008 14 / 26
Kwadratury Gaussa
Kwadratury Gaussa
Wiemy jak tworzyć kwadratury

n
b

f (x)w(x) dx H" Aif (xi)
a
i=0
(w jest ustaloną nieujemną funkcją wagową), dokładne dla każdego
wielomianu f " n. We wzorach tych układ węzłów xi jest dowolny
i określa jednoznacznie współczynniki Ai.
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 Ai
byłyby równe:

n
b

f (x)w(x) dx H" A f (xi) .
a
i=0
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) Podstawy metod numerycznych 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

b
f , g = f (x)g(x)w(x) dx .
a
Układ wielomianów pn (n = 0, 1, . . .) jest ortogonalny na przedziale [a, b]
z wagą w, jeżeli stopień każdego z nich jest równy n oraz
pm, pn = 0 dla m = n .

Każdy wielomian jest określony z dokładnością do stałego czynnika.
Czynniki te wybieramy zazwyczaj tak, aby pn, pn = 1, albo tak aby pn
był wielomianem standardowym, tj. współczynnik w pn(x) przy xn jest
równy 1.
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 2007/2008 16 / 26
Kwadratury Gaussa Wielomiany ortonormalne
Twierdzenie
Ciąg {pn} wielomianów ortogonalnych standardowych w przedziale [a, b]
z wagą w można wyznaczyć ze wzorów:
p0(x) = 1 ,
p1(x) = x - a1 ,
pn(x) = (x - an)pn-1(x) - bnpn-2(x) (n e" 2) ,
gdzie
xpn-1, pn-1 xpn-1, pn-2
an = , bn = .
pn-1, pn-1 pn-2, pn-2
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 2007/2008 17 / 26
Kwadratury Gaussa Wielomiany ortonormalne
Wielomiany Legendre a
Wielomiany Legendre a są ortogonalne w przedziale [-1, 1] z wagą 1.
P0(x) = 1 , P1(x) = x ,
1 3
P2(x) = x2 - , P3(x) = x3 - x ,
3 5
6 3 10 5
P4(x) = x4 - x2 + , P5(x) = x5 - x3 + x .
7 35 9 21
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 2007/2008 18 / 26
Kwadratury Gaussa Wielomiany ortonormalne
Wielomiany Czebyszewa
1
"
Wielomiany Czebyszewa są ortogonalne w przedziale [-1, 1] z wagą .
1-x2
T0(x) = 1 , T1(x) = x ,
T2(x) = 2x2 - 1 , T3(x) = 4x3 - 3x ,
T4(x) = 8x4 - 8x2 + 1 , T5(x) = 16x5 - 20x3 + 5x ,
T6(x) = 32x6 - 48x4 + 18x2 - 1 .
Dla n e" 2 można je określić wzorem rekurencyjnym
Tn(x) = 2xTn-1(x) - Tn-2(x) .
Twierdzenie
W przedziale [-1, 1] wielomiany Czebyszewa wyrażają się wzorem
Tn(x) = cos(n arccos x) .
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 2007/2008 19 / 26
Kwadratury Gaussa Kwadratury Gaussa
Kwadratury Gaussa
Twierdzenie
Jeżeli węzły x0, x1, . . . , xn są zerami (n + 1)-szego wielomianu
ortogonalnego pn+1 na przedziale [a, b] z wagą w, to kwadratura Gaussa

n
b

f (x)w(x) dx H" Aif (xi)
a
i=0
o współczynnikach

n
b

x - xj
Ai = w(x) dx
xi - xj
a
j=0
j=i

jest dokładna dla każdej funkcji f " 2n+1.
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 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) Podstawy metod numerycznych 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ć

n
d

f (t) dt H" Aif (ti) .
c
i=0
Aby otrzymać wzór dla przedziału [a, b] dokonujemy podstawienia
(b - a)t + ad - bc
x = (t) = ,
d - c
co daje

n
b d

b - a b - a
f (x) dx = f ((t)) dt H" Aif ((b-a)ti+ad-bc ) .
d-c
d - c d - c
a c
i=0
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 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

n
b

lim Anif (xni) = f (x)w(x) dx .
n"
a
i=0
Twierdzenie
2n
Jeżeli f " C [a, b], to dla kwadratury Gaussa z n węzłami zachodzi

n-1
(2n)
b b

f (�)
f (x)w(x) dx = Aif (xi) + q2(x)w(x) dx ,
(2n)!
a a
i=0
n-1
gdzie � " (a, b) i q(x) = (x - xi).
i=0
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 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 2n podprzedziałów równej
długości (n e" 0). Przyjmujemy:

2n-1

b - a
R(n, 0) = hn f (a) + f (a + ihn) + f (b) gdzie hn = .
2n
i=1
W najprostszym przypadku, dla n = 0 mamy prosty wzór trapezów:
1
R(0, 0) = [f (a) + f (b)] .
2
R(1, 0), R(2, 0), . . . obliczamy rekurencyjnie, aby uniknąć obliczania
wartości funkcji f w tych samych punktach:
2n-1

1
R(n, 0) = R(n - 1, 0) + hn f (a + (2i - 1)hn) .
2
i=1
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 2007/2008 24 / 26
Metoda Romberga
Następne serie przybliżeń całki (dla m > 0) obliczamy ze wzorów:
4mR(n, m - 1) - R(n - 1, m - 1)
R(n, m) = .
4m - 1
metoda Romberga
Require: a, b, f (x), M

b
Ensure: R(n, m) H" f (x) dx
a
h �! b - a
R(0, 0) �! (b - a)(f (a) - f (b))/2
for n �! 1, 2, . . . , M do
h �! h/2
2
n-1
R(n, 0) �! R(n - 1, 0)/2 + h f (a + (2i - 1)h)
i=1
for m �! 1, 2, . . . , n do
R(n, m) �! (4mR(n, m - 1) - R(n - 1, m - 1))/(4m - 1)
end for
end for
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 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
2m
dostatecznie wiele pochodnych. Jeżeli f " C [a, b], to R(n, m) przybliża
całkę z błędem rzędu O(h2m), gdzie h = (b - a)/2n. Metoda jest także
sensowna dla funkcji, które są jedynie ciągłe:
Twierdzenie
Jeżeli f " C [a, b], to

b
lim R(n, m) = f (x) dx .
n"
a
dr K. A. Smoliński (WSInf) Podstawy metod numerycznych 2007/2008 26 / 26