10.3 Kwadratury Gaussa
Dla dowolnej n, niech x∗0, x∗1, . . . , x∗nbędą (różnymi) zerami (n+1)-szego wielomianu pn+1w ciągu wielomianów ortogonalnych w przestrzeni uni-tarnej L2,ρ(a, b), tzn.
pn+1(x) = (x − x∗0)(x − x∗1) · · · (x − x∗n). (10.1) Ponieważ x∗j leżą w przedziale (a, b), możemy mówić o kwadraturze interpolacyjnej opartej na tych węzłach.
Definicja 10.3 Kwadraturę interpolacyjną QGSn opartą na zerach wie-lomianu ortogonalnego pn+1 nazywamy kwadraturą Gaussa.
Okazuje się, że właśnie kwadratury Gaussa mają najwyższy rząd.
Dokładniej, mamy następujące twierdzenie.
Twierdzenie 10.2 Kwadratura Gaussa QGSn ma najwyższy rząd spo-śród wszystkich kwadratur opartych na węzłach o łącznej krotności n+1, oraz
rz(QGSn ) = 2n + 2.
Dowód Wobec Lematu 10.1(ii) wystarczy pokazać, że kwadratura QGSn jest dokładna dla każdego wielomianu stopnia nie większego niż 2n + 1.
Niech f ∈ Π2n+1. Niech wf ∈ Πn będzie wielomianem interpolującym f w węzłach-zerach x∗0, . . . , x∗n wielomianu pn+1. Jeśli deg f ≤ n to oczywiście wf = f i QGSn (f) = Sρ(wf) = Sρ(f). Jeśli zaś
n+ 1 ≤ deg f ≤ 2n + 1,
to f − wf jest wielomianem stopnia tego samego co f i zeruje się w x∗j, 0 ≤ j ≤ n. Stąd
f(x) − wf(x) = (x − x∗0)(x − x∗1) · · · (x − x∗n)g(x), gdzie g jest wielomianem,
deg g = deg f − (n + 1) ≤ (2n + 1) − (n + 1) = n.
124 ROZDZIAŁ 10. CAŁKOWANIE A APROKSYMACJA Korzystając z (10.1) i faktu, że pn+1jest prostopadły do Πn, ostatecznie otrzymujemy
Sρ(f) − QGSn (f) = Z b
a (f(x) − wf(x))ρ(x) dx
= Z b
a (x − x∗0) · · · (x − x∗n)g(x)ρ(x) dx
= Z b
a pn+1(x)g(x)ρ(x) = hpn+1, gi = 0, co kończy dowód. 2
Zajmiemy się teraz błędem kwadratur Gaussa. Pokażemy, że również ze względu na błąd mają one dobre własności.
Twierdzenie 10.3 Jeśli f ∈ C2n+2([a, b]) to błąd kwadratury Gaussa QGSn wyraża się wzorem
Sρ(f) − QGSn (f) = kpn+1k2f(2n+2)(ξ) (2n + 2)!, gdzie ξ∈ [a, b]. Stąd, w szczególności,
max
f ∈FM2n+1([a,b])|Sρ(f) − QGSn (f)| = Mkpn+1k2 (2n + 2)! .
Dowód Niech wf ∈ Πnbędzie wielomianem interpolującym f w zerach x∗j wielomianu pn+1. Niech ˜wf ∈ Π2n+1 będzie z kolei wielomianem (Hermite’a) interpolującym f w dwukrotnych węzłach x∗j, tzn. takim, że ˜wf(x∗j) = f(x∗j) i ˜wf′(x∗j) = f′(x∗j), 0 ≤ j ≤ n. Ponieważ rz(QGSn ) = 2n + 2 to QGSn ( ˜wf) = Sρ( ˜wf), a stąd i ze wzoru na błąd interpolacji Hermite’a mamy
Sρ(f) − QGSn (f) = Sρ(f) − QGSn (wf) = Sρ(f) − QGSn ( ˜wf)
= Sρ(f) − Sρ( ˜wf) = Z b
a (f(x) − ˜wf(x))ρ(x) dx
= Z b
a (x − x∗0)2· · · (x − x∗n)2f(x∗0, . . . , x∗n, x)ρ(x) dx
= Z b
a p2n+1(x)ρ(x)f(x∗0, . . . , x∗n, x) dx.
10.3. KWADRATURY GAUSSA 125 Ponieważ funkcja p2n+1(x)ρ(x) jest prawie wszędzie dodatnia, możemy teraz zastosować twierdzenie o wartości średniej, aby ostatecznie otrzy-mać
Sρ(f) − QGSn (f) = f(2n+2)(ξ) (2n + 2)!
Z b
a p2n+1(x)ρ(x) dx
= f(2n+2)(ξ)
(2n + 2)! kpn+1k2, co kończy dowód. 2
Na końcu, zwrócimy jeszcze uwagę na inną, bardzo ważną własność kwadratur Gaussa, a mianowicie, że ich współczynniki są dodatnie.
Rzeczywiście, zapisując
QGSn (f) = Xn
j=0
ajf(xj) i podstawiając
fj(x) = (x − x∗0)2· · · (x − x∗j−1)2(x − x∗j+1)2· · · (x − x∗n)2 mamy, że fj ∈ Π2n i fj jest prawie wszędzie dodatnia. Stąd
0 < Sρ(fj) = QGSn (fj) = ajfj(x∗j)
i aj > 0, bo fj(x∗j) > 0. Przypomnijmy, że dodatniość współczynni-ków kwadratury ma duże znaczenie przy numerycznym ich obliczaniu, zwłaszcza gdy funkcja podcałkowa f ma stały znak, zob. Rozdział 2.5.2.
Mimo niewątpliwych zalet kwadratur Gaussa, ich stosowalność ogra-niczają trudności w wyliczeniu pierwiastków wielomianów ortogonal-nych, gdy stopień wielomianu jest duży. Wyjątkiem są tutaj kwadra-tury interpolacyjne oparte na zerach wielomianów Czebyszewa, zob. U.
10.3.
Uwagi i uzupełnienia
U. 10.1 Pokażemy teraz, że wielomiany ortogonalne {pk}k≥0 w danej prze-strzeni L2,ρ(a, b) spełniają następującą formulę trójczłonową. Załóżmy dla
126 ROZDZIAŁ 10. CAŁKOWANIE A APROKSYMACJA uproszczenia, że współczynnik przy xk w wielomianie pk jest dla każdego k jednością. Wtedy istnieją liczby βk (dla k ≥ 1) i γk > 0 (dla k≥ 2) takie, że
p0(x) = 1,
p1(x) = (x − β1), (10.2)
pk(x) = (x − βk)pk−1(x) − γkpk−2(x), k≥ 2.
Aby to pokazać zauważmy, że pk można dla k ≥ 1 przedstawić w postaci rozwinięcia
pk(x) = (x − ck−1)pk−1(x) +
k−2X
j=0
cjpj(x).
Mnożąc skalarnie obie strony tego równania przez ps dla 0 ≤ s ≤ k − 3, otrzymujemy
0 = hpk(x), ps(x)i = h(x − ck−1)pk−1(x), ps(x)i + cshps(x), ps(x)i.
Wobec tego, że (x−ck−1)ps(x) jest wielomianem stopnia mniejszego niż k−1, mamy
h(x − ck−1)pk−1(x), ps(x)i = hpk−1(x), (x − ck−1)ps(x)i = 0,
a stąd cshps(x), ps(x)i = 0 i cs = 0. Możemy więc napisać, że dla k = 1 mamy p1(x) = (x − β1), a dla k ≥ 2,
pk(x) = (x − βk)pk−1(x) − γkpk−2(x), (10.3) gdzie βk = ck−1 i γk = ck−2. Aby jawnie wyznaczyć βk i γk, pomnożymy skalarnie obie strony równania (10.3) kolejno przez pk−1i pk−2. Otrzymujemy
0 = hpk(x), pk−1i
= h(x − βk)pk−1(x), pk−1(x)i − γkhpk−2(x), pk−1(x)i
= hxpk−1(x), pk−1(x)i − βkhpk−1(x), pk−1(x)i, czyli
βk = hxpk−1(x), pk−1(x)i hpk−1(x), pk−1(x)i , oraz
0 = hpk(x), pk−2(x)i
= h(x − βk)pk−1(x), pk−2(x)i − γkhpk−2(x), pk−2(x)i
= hpk−1(x), xpk−2(x)i − γkhpk−2(x), pk−2(x)i,
10.3. KWADRATURY GAUSSA 127 a stąd i z równości hpk−1(x), xpk−2(x)i = hpk−1(x), pk−1(x)i,
γk = hpk−1(x), pk−1(x)i hpk−2(x), pk−2(x)i.
Zauważmy, że z formuły trójczłonowej wynika w szczególności algorytm wyznaczenia ciągu wielomianów ortogonalnych. Wystarczy bowiem wyzna-czać kolejne współczynniki βki γk(obliczając odpowiednie iloczyny skalarne hxpk(x), pk(x)i i hpk(x), pk(x)i) i stosować wzór rekurencyjny (10.2). Do-dajmy jeszcze, że w obliczeniach numerycznych najlepiej jest przechowywać informację o ciągu {pk}k≥0 po prostu w postaci liczb βk i γk.
U. 10.2 Wygodnie jest posłużyć się wielomianami ortogonalnymi w przy-padku, gdy chcemy znaleźć najlepszą aproksymację danej funkcji f wielomia-nem ustalonego stopnia n, i błąd mierzymy w normie przestrzeni L2,ρ(a, b).
Rzeczywiście, jak wiadomo, najlepszą aproksymacją dla f w przestrzeni Πn
jest jej rzut prostopadły na Πn. Ponieważ n + 1 początkowych wielomianów ortogonalnych pk, 0 ≤ k ≤ n, tworzy bazę w Πn, rzut ten wyraża się wzorem
wf∗ = Xn k=0
hf, pki hpk, pkipk.
U. 10.3 Zachodzi następujące twierdzenie Łuzina. Niech przedział [a, b] bę-dzie skończony. Niech Qn(f) =Pnj=0a(n)j f (x(n)j ) będzie takim ciągiem kwa-dratur, że:
(i) wszystkie współczynniki a(n)j są dodatnie,
(ii) rząd kwadratur Qn rośnie do nieskończoności gdy n → ∞.
Wtedy ciąg Qn(f) zbiega do Rabf (x)dx dla każdej funkcji ciągłej f .
Zauważmy, że twierdzenie to stosuje się do ciągu kwadratur Gaussa QGSn , ale nie do ciągu kwadratur Newtona-Cotesa QN Cm , ponieważ w tych ostatnich pojawiają się dla dużych n współczynniki ujemne.
W praktyce, najczęściej stosuje się ciąg kwadratur interpolacyjnych opar-tych na zerach kolejnych wielomianów Czebyszewa, ponieważ zera te dane są jawnie i “zagęszczają się”, tzn. zera wielomianu Czebyszewa Tksą też zerami wielomianu T2k. Powstające w ten sposób kwadratury noszą nazwę formuł Clanshow-Curtis’a. Są one w pewnym sensie uniwersalne, bowiem posiadają optymalną szybkość zbieżności n−(r+1) w klasach FMr ([a, b]) dla dowolnych r i M .
128 ROZDZIAŁ 10. CAŁKOWANIE A APROKSYMACJA U. 10.4 Jeśli przedział całkowania jest skończony i waga jest jednostkowa to kwadratury Gaussa (a dokładniej kwadratury Legendre’a) można użyć do tworzenia kwadratur złożonych ¯QGSr,k, gdzie k oznacza liczbę podprzedzia-łów. Łatwo widać, że dla f ∈ FM2r+1([a, b]), błąd takiej kwadratury można
czyli jest on porównywalny do błędu “zwykłej” złożonej kwadratury interpo-lacyjnej. Jednak złożona kwadratura Legendre’a korzysta z dwa razy mniej-szej liczby węzłów.
U. 10.5 Błąd złożonej kwadratury Legendre’a w klasie FM2r+1([a, b]) można podać dokładnie. Wystarczy wykorzystać wzory z Ćw. 10.4 i 10.5, aby otrzy-mać
(a) jeśli kwadratura jest dokładna dla dowolnych n + 1 wielomianów tworzą-cych bazę w Πn, to jest ona rzędu co najmniej n + 1, oraz
(b) jeśli kwadratura jest rzędu n + 1 to jest ona niedokładna dla każdego wielomianu stopnia dokładnie n + 1.
Ćw. 10.2 Uzasadnić, że kwadratura prostokątów jest kwadraturą Legen-dre’a, natomiast żadna z kwadratur Newtona-Cotesa QN Cn nie jest kwadra-turą Gaussa.
Ćw. 10.3 Załóżmy, że dane są liczby βk i γk definiujące ciąg wielomia-nów ortogonalnych przez formułę trójczłonową. Zaproponować ekonomiczny (tzn. o koszcie proporcjonalnym do n) algorytm obliczania wartości n-tego wielomianu ortogonalnego w danym punkcie x, wykorzystujący formułę trój-członową.
10.3. KWADRATURY GAUSSA 129
Wskazówka. Wykorzystać fakt, że dla n-tego wielomianu Legendre’a mamy R1
będzie kwadraturą interpolacyjną opartą na zerach (n+1)-ezego wielomianu Legendre’a. Niech −∞ < a < b < ∞. Pokazać, że wtedy kwadratura
Q˜GSn (f) = b− a
jest kwadraturą Gaussa opartą na n + 1 węzłach, dla całki na przedziale (a, b) z wagą jednostkową. Ponadto, jeśli f ∈ C(2n+2)([a, b]), to ortogonalnego Legendre’a (tzn. na przedziale [−1, 1] z wagą 1). Niech dla 0 ≤ j ≤ n,
wj = Z 1
−1
(x − x0) · · · (x − xj−1)(x − xj+1) · · · (x − xn) (xj− x0) · · · (xj − xj−1)(xj− xj+1) · · · (xj− xn)dx.
Pokazać, że jeśli f i g są wielomianami stopnia nie większego niż n, to ich iloczyn skalarny w L2,1(−1, 1),
hf, gi =
130 ROZDZIAŁ 10. CAŁKOWANIE A APROKSYMACJA
Rozdział 11
Iteracje dla równań liniowych
Algorytmy rozwiązywania układów równań liniowych postaci A~x = ~b,
gdzie A jest nieosobliwą macierzą rzeczywistą n × n, a ~b jest wektorem rzeczywistym w Rn, które rozpatrywaliśmy w Rozdziałach 3, 4 i 5, na-leżą do grupy algorytmów dokładnych albo bezpośrednich. To znaczy, że po wykonaniu skończonej liczby dopuszczalnych operacji elementarnych dostajemy w arytmetyce idealnej dokładne rozwiązanie
~x∗ = A−1~b.
W tym rozdziale zajmiemy się algorytmami iteracyjnymi rozwiązy-wania układów równań liniowych. Polegają one na tym, że, startując z pewnego przybliżenia początkowego ~x0, konstruuje się ciąg kolejnych przybliżeń
~xk = Φk(A,~b; ~x0), k = 1, 2, . . . , które w granicy osiągają rozwiązanie dokładne,
k→∞lim ~xk = ~x∗.
131
132 ROZDZIAŁ 11. ITERACJE DLA RÓWNAŃ LINIOWYCH
11.1 Kiedy stosujemy iteracje?
Jasne jest, że algorytmy iteracyjne stosujemy wtedy, gdy są one konku-rencyjne w stosunku do algorytmów bezpośrednich. Dlatego przekształ-cenia Φk należy wybierać tak, aby kolejne przybliżenia można było ła-two obliczać i jednocześnie kolejne błędy k~xk− ~x∗k szybko zbiegały do zera.
Zwykle zakłada się również, że dokładne rozwiązanie ~x∗ jest punk-tem stałym przekształcenia Φk(A,~b; ·). Wtedy kolejne błędy spełniają zależność
~xk− ~x∗ = Φk(A,~b; ~x0) − Φk(A,~b; ~x∗).
Jeśli teraz Φk(A,~b; ·) są lipschitzowskie ze stałymi mk <+∞, tzn. dla pewnej normy wektorowej k · k mamy
kΦk(~x) − Φk(~y)k ≤ mkk~x − ~yk, ∀~x, ~y, to
k~xk− ~x∗k ≤ mkk~x0− ~x∗k.
Warunek limk→∞mk = 0 jest więc dostateczny na to, aby metoda była zbieżna dla dowolnego przybliżenia początkowego ~x0, przy czym szyb-kość zbieżności zależy od tego, jak szybko mk maleją do zera. Dla więk-szości stosowanych metod Φk jest funkcją liniową błędu początkowego, tzn.
Φk(A,~b; ~x0− ~x∗) = Mk(~x0− ~x∗),
gdzie Mk jest pewną macierzą. Wtedy jako mk można przyjąć normę tej macierzy,
mk = kMkk = sup
k~xk=1kMk~xk.
Dla ilustracji, rozpatrzmy ogólną metodę iteracji prostej, w której
~xk = B~xk−1 + ~c, (11.1) dla pewnej macierzy B wymiaru n × n i wektora ~c ∈ Rn. W tym przypadku
~xk− ~x∗ = Bk(~x0− ~x∗),
11.1. KIEDY STOSUJEMY ITERACJE? 133 a stąd i z nierówności kBkk ≤ kBkk, mamy
k~xk− ~x∗k ≤ kBkkk~x0− ~x∗k.
Warunkiem dostatecznym zbieżności iteracji prostych jest więc kBk <
1. Mówimy, że metoda jest zbieżna liniowo z ilorazem kBk.
Przykład 11.1 Rozkładając macierz A = (ai,j)ni,j=1 na sumę A = D + C,
gdzie D jest macierzą diagonalną składającą się z wyrazów stojących na głównej przekątnej macierzy A, układ A~x = ~b jest równoważny układowi
D~x = −C~x + ~b,
a stąd (o ile na przekątnej macierzy A nie mamy zera) otrzymujemy metodę iteracyjną
~xk = B~xk−1 + ~c,
gdzie B = −D−1C i ~c = D−1~b, zwaną metodą Jacobiego.
W metodzie Jacobiego warunek dostateczny zbieżności, kBk < 1, jest spełniony wtedy, gdy macierz A ma dominującą przekątną, tzn.
gdy
2|ai,i| >
Xn j=1
|ai,j|, 1 ≤ i ≤ n. (11.2) Rzeczywiście, ponieważ wyraz (i, j) macierzy D−1C wynosi 0 dla i = j i ai,j/ai,i dla i 6= j, to
kD−1Ck∞ = max
1≤i≤n
Xn j=1,j6=i
|ai,j|/|ai,i|
= max
1≤i≤n
Xn j=1
|ai,j|/|ai,i| − 1 < 1, przy czym ostatnia nierówność wynika z (11.2).
Inne przykłady iteracji prostych podane są w U. 11.3 i Ćw. 11.2.
Zastanówmy się teraz nad złożonością metod iteracyjnych. Ponieważ możemy jedynie znaleźć pewne przybliżenie rozwiązania dokładnego ~x∗,
134 ROZDZIAŁ 11. ITERACJE DLA RÓWNAŃ LINIOWYCH przez złożoność metody będziemy rozumieli koszt kombinatoryczny ob-liczenia ~xkz zadaną dokładnością ε > 0. Dla uproszczenia założymy, że medoda jest zbieżna liniowo z ilorazem m. Zauważmy, że aby zreduko-wać błąd początkowy do ε > 0, wystarczy wykonać k = k(ε) iteracji, gdzie k spełnia
mkk~x0− ~x∗k ≤ ε, czyli
k ≥ log(1/ε) − log(1/k~x0 − ~x∗k)
log(1/m) .
Liczba ta zależy więc w istotny sposób od błędu początkowego i (przede wszystkim) od stałej Lipschitza m, natomiast zależność od dokładności εi wymiaru n układu jest dużo mniej istotna. Zakładając, że koszt jed-nej iteracji wynosi c = c(n) (zwykle c(n) jest tym mniejszy, im mniejsza jest liczba niezerowych elementów macierzy A), złożoność metody jest proporcjonalna do
c(n) log(1/ε) log(1/m).
Stąd oczywisty wniosek, że metody iteracyjne warto stosować zamiast metod bezpośrednich w przypadku gdy
• wymiar n układu A~x = ~b jest “duży”, oraz
• macierz A układu jest “rozrzedzona”, tzn. ma stosunkowo nie-wielką liczbę elementów niezerowych, np. proporcjonalną do n.
Układy o tych własnościach powstają często przy numerycznym roz-wiązywaniu równań różniczkowych cząstkowych.
Zaletą metod iteracyjnych jest również ich prostota, przez co są one łatwe do zaprogramowania.