5. Całkowanie numeryczne
W tej części wykładu omówimy najważniejsze metody numerycznego obliczania całek funkcji rzeczywistych jednej
zmiennej. Najpierw przypomnimy de
�nicję całki oznaczonej:
łR b �
DEFINICJA 26. Niech f będzie funkcjąokreśloną i ograniczonąna[a; b]. Całką oznaczonąRiemanna f(x)dx
a
funkcji f na [a; b] nazywamy granicę sum
n
X
Sn = (xi+1 Ą xi) f (~i) , (53)
x
i=0
gdzie a = x0 przy n !1, a xi 2 [xi; xi+1] są dowolnymi punktami pośrednimi.
~
Można pokazać, że funkcja f ograniczona na przedziale [a; b] jest na nim całkowalna w sensie Riemanna, tzn.
R
b
f(x)dx = limn!1 Sn istnieje i jest skończona, wtedy i tylko wtedy, gdy punkty nieciągłości funkcji f w [a; b]
a
tworzą zbiór miary zero.
Funkcjonał liniowy, którym jest całka, będziemy oznaczali przez I, tzn.
Z
b
I(f) = f(x)dx .
a
Wcałkowaniu numerycznym powyższą całkę przybliża się funkcjonałami Q postaci
n0 n1 nk
X X X
Q(f) = Ai;0f (xi;0) + Ai;1f (xi;1) +::: + Ai;kf (xi;k) ,
i=0 i=0 i=0
które nazywamy kwadraturami liniowymi, natomiast wspólczynniki Ai;j i punkty xi;j odpowiednio współczynnikami i
węzłami kwadratury Q. Resztą kwadratury będziemy nazywać wartość
R(f) =I(f) Ą Q(f) . (54)
Najczęściej stosowanymi są kwadratury wykorzystujące tylko wartości funkcji f, tzn.
n
X
Q(f) = Aif (xi) , (55)
i=0
które można uważać za uogólnienie sum częściowych Riemanna (53). Takie kwadratury mają 2n +3 parametrów:
liczba n, węzły x0; x1; :::; xn i współczynniki A0; A1; :::; An.
Wybór konkretnych wartości parametrów odbywaćsiębędzie na podstawie różnych żądańzwiązanych z całkowaniem
numerycznym, np. ograniczony wybór węzłów, łatwość wyznaczenia współczynników (równość wszystkich Ai) lub
węzłów (np. równoodległe), minimalizacja wartości bezwzględnej reszty lub dokładność dla wielomianów możliwie
wysokiego stopnia przy ustalonej liczbie węzłów.
Dalej będziemy rozważać również nieco ogólniejszą postać całek, tzn.
Z
b
Ip(f) = p(x)f(x)dx ,
a
gdzie funkcja wagowa p jest nieujemna na [a; b], zeruje się w nim w skończonej liczbie punktów i jest całkowalna, tzn.
R
b
p(x)dx < 1.
a
5.1. Zbieżność ciągu kwadratur
Rozpatrzmy ciąg kwadratur Qn, n =0; 1; :::
n
ł �
X
Qn(f) = A(n)f x(n) (56)
i i
i=0
R
b
przybliżających całkę Ip(f) = p(x)f(x)dx.
a
Zbadamy warunki zbieżności ustalonego ciągu kwadratur do całek funkcji ciągłych na [a; b]. Załóżmy, że dane są
1
nieskończone macierze trójkątne dolne węzłów x(j) i współczynników A(j)
i i
x(0) A(0)
0 0
x(1) x(1) A(1) A(1)
0 1 0 1
oraz
x(2) x(2) x(2) A(2) A(2) A(2)
0 1 2 0 1 2
::: ::: ::: ::: ::: ::: ::: :::
de
�niujące ciąg (56).
TWIERDZENIE 27. Ciąg kwadratur (56) jest zbieżny dla dowolnych funkcji f ciągłychnaprzedziale [a; b], tzn.
Z
n
ł � b
X
lim A(n)f x(n) = p(x)f(x)dx
i i
n!1
a
i=0
wtedy i tylko wtedy, gdy
(a) ciąg (56) jest zbieżny dla dowolnego wielomianu,
(b) istnieje stała K taka, że dla n =0; 1; ::: zachodzi nierówność
n
Ż Ż
X
ŻA Ż
(n)
Ż Ż � K .
i
i=0
DOWÓD pomijamy.
Jako wniosek z powyższego można otrzymać następujący:
WNIOSEK 28. Jeśli wszystkie współczynniki A(n) ciagu kwadratur (56) są nieujemne, to warunkiem koniecznym i
i
dostatecznym zbieżności tego ciągu dla dowolnych funkcji ciągłych na [a; b] jest jego zbieżność dla dowolnych wielomi-
anów.
5.2. Kwadratury interpolacyjne
DEFINICJA 29. Kwadraturami intrepolacyjnym nazywamy kwadratury otrzymane przez całkowanie wielomi-
anów Hermite a (w szczególności Lagrange a) interpolujących funkcję podcałkową
Pnf.
Założmy, że xi, i = 0; 1; :::; n są ustalonymi �i-krotnymi węzłami, n + �i = m, wówczas wielomian
i=0
interpolacyjny Hermite a Hm jest określony jednoznacznie warunkami
(s)
Hm (xi) =f(s) (xi) , i =0; 1; :::; n, s =0; 1; :::; �i
i spełnia równanie
f(x) =Hm(x) +Rm(x) ,
gdzie Rm(x) jest resztą wzoru interpolacyjnego.
Całkując powyższe równanie otrzymamy
Z Z Z
b b b
p(x)f(x)dx = p(x)Hm(x)dx + p(x)Rm(x)dx ,
a a a
zatem kwadratury interpolacyjne są postaci
Q(f) =Ip (Hm) ,
a dla wyrażenia ich reszty można skorzystać z różnych przedstawień reszty wzoru interpolacyjnego Hermite a (patrz
J.M.JANKOWSCY - Przegląd metod i algorytmów numerycznych. Część 1) - jeśli f 2 Cm+1 [a; b], to
Z
b
1
Rm(f) = p(x)!� (x)f(m+1) (x) dx , (57)
m+1
(m +1)!
a
Qn
gdzie !� (x) = (x Ą xi)�i.
m+1 i=o
Najczęściej używa się kwadratur z węzłami jednokrotnymi, czyli gdy wielomian interpolacyjny Hermite a upraszcza
sie do wielomianu interpolacyjnego Lagrange a. Wówczas
n n
X Y
x Ą xi
f(x) = f (xk) + Rn(x) .
xk Ą xi
k=0 i=0;i6=k
2
Całkując otrzymamy kwadraturę postaci (55)
n
X
Q(f) = Akf (xk) ,
k=0
gdzie xi są ustalonymi węzłami, a współczynniki Ak sa określone wzorem
Z
n
b
Y
x Ą xi
Ak = p(x) dx .(A)
xk Ą xi
a
i=0;i6=k
PRZYKA
AD 1: Funkcję f przybliżymy wielomianem interpolacyjnym Lagrange a opartym na węzłach a i b.
x Ą b x Ą a
f(x) ź W1(x) =f(a) + f(b) .
a Ą b b Ą a
Całkując wielomian W1 otrzymujemy kwadraturę
Z � ś
b
x Ą b x Ą a b Ą a b Ą a
Q(f) =I (W1) = f(a) + f(b) dx = f(a) + f(b) . (58)
a Ą b b Ą a 2 2
a
Jest to tzw. wzór trapezów, gdyż przybliżenie to można geometrycznie zinterpretować jako przybliżenie pola trapezu
krzywoliniowego polem trapezu o podstawach f(a) i f(b) oraz wysokości b Ą a.
PRZYKA
AD 2: Niech S będzie funkcją sklejaną stopnia pierwszego z węzłami a = x0 < x1 < ::: < xn = b
interpolującą funkcję f. Wkażdym z przedziałów [xi; xi+1], i =0; 1; :::; n Ą 1 funkcja S jest określona wzorem
f (xi+1) Ą f (xi)
Si(x) =f (xi) + (x Ą xi) .
xi+1 Ą xi
Całkując S otrzymujemy kwadraturę postaci
Z � ś
nĄ1 nĄ1
xi+1
X X
xi+1 Ą xi xi+1 Ą xi
Q(f) =I (S) = S(x)dx = f (xi) + f (xi+1) . (59)
2 2
xi
i=0 i=0
Dla funkcji f dodatniej na [a; b] takie przybliżenie odpowiada geometrycznie przybliżaniu pola trapezu krzywolin-
iowego sumą pól trapezów o wysokościach xi+1 Ą xi oraz podstawach f (xi) i f (xi+1).
UWAGA: W myśl przyjętego określenia kwadratury interpolacyjnej, kwadratura określona wzorem (59) dla n >1
nie jest kwadraturą interpolacyjną.
5.2.1. Kwadratury Newtona-Cotesa
DEFINICJA 30. Kwadraturami Newtona-Cotesa są nazywane kwadratury Q(f) =I (Wn), gdzie Wn jest wielo-
mianem interpolacyjnym Lagrange a funkcji f opartym na równoodległych węzłach
xk = a + kh; i =0; 1; :::; n ,
gdzie h =(b Ą a)=n.
UWAGA: Kwadratury oparte na powyższych węzłach są nazywane zamkniętymi. Rozważa się również otwarte
kwadratury Newtona-Cotesa oparte na węzłach
2k +1
xk = a + h0; k =0; 1; :::; n Ą 1
2
h a+b
gdzie h0 =(b Ą a)=(n +1). Np. kwadratura prostokątów (n =0), która posiada tylko1 węzeł x0 = a + = :
2 2
� ś
Z Z Z
b b b
a + b
f(x)dx ź W0(x)dx = f (x0) dx =(b Ą a)f ,
2
a a a
(bĄa)3
a jej błąd wynosi R(f) = f00(), 2 (a; b).
24
Zapiszmy teraz wzór interpolacyjny Lagrange a postaci (14) dla węzłów równoodległych xk = a + kh, stosując
podstawienie x = a + qh, a mianowicie
n n
X Y
q Ą i
Wn(x) =Wn(a + qh) = f (xk) .
k Ą i
k=0 i=0;i6=k
3
Całkując prawąstronępowyższego wzoru, otrzymujemy kwadraturę(55) zewspółczynnikami Ak określonymi wzorem
Z
n
n
Y
q Ą i
Ak = h dq . (60)
k Ą i
0
i=0;i6=k
UWAGA: A
atwo sprawdzić, że Ak = AnĄk.
Przedstawimy teraz najprostsze kwadratury interpolacyjne Newtona-Cotesa i obliczymy ich reszty. Będziemy
musieli wykorzystać twierdzenie o wartości średniej dla całek. Przypomnijmy je bez dowodu w postaci poniższego
LEMAT 31. Jeżeli w przedziale [a; b] funkcja f jest ciągła, a funkcja g całkowalna i nieujemna (lub niedodatnia),
to istnieje wewnątrz tego przedziałutaki punkt , że
Z Z
b b
f(x)g(x)dx = f () g(x)dx .
a a
PRZYKA
AD 1: n =1 - węzłami kwadratury są wówczas krańce przedziału całkowania, tj. x0 = a, x1 = b. Na
podstawie wzoru (60) możemy obliczyć
Z
1
q Ą 1 b Ą a
A0 = (b Ą a) dq =
0 Ą 1 2
Z0
1
q Ą 0 b Ą a
A1 = (b Ą a) dq =
1 Ą 0 2
0
Zatem kwadratura Newtona-Cotesa przy n =1 jest postaci
b Ą a
Q(f) =I (W1) = (f(a) +f(b)) (61)
2
i pokrywa się z wzorem trapezów (58):
Aby obliczyć resztę kwadratury trapezów, skorzystamy ze wzoru (57). Zakładając, że f 2 C2 [a; b], otrzymamy
stąd i z lematu 31, że
Z Z
b b
1 f00 ()
R(f) = (x Ą a)(x Ą b)f00(x)dx = (x Ą a)(x Ą b)dx =
2! 2!
a a
(b Ą a)3
= Ą f00 () , gdzie 2 (a; b) :
12
PRZYKA
AD 2: n =2 - węzłami kwadratury są x0 = a, x1 =(a+b)=2, x2 = b. Na podstawie wzoru (63) możemy
obliczyć
b Ą a
A0 =
6
4(b Ą a)
A1 =
6
A2 = A0
Zatem kwadratura Newtona-Cotesa przy n =2 jest postaci
� � ś ś
b Ą a a + b
Q(f) =I (W2) = f(a) +4f + f(b) . (62)
6 2
Powyższa kwadratura jest nazywana wzorem parabol lub wzorem Simpsona.
Zbadajmy teraz błąd kwadratury parabol. Załóżmy, że f 2 C4 [a; b] i zde
�niujmy wielomian
� ś
a + b
H3(x) =W2(x) +K (x Ą a) x Ą (x Ą b) , (63)
2
Ąa+bó Ąa+bó
0
wktórymstałą K dobieramy tak, by H3 2 = f0 2 . Zatem H3 jest wielomianem interpolacyjnym Hermite a
funkcji f opartym na węzłach a,a+b i b okrotnościach odpowiednio 1,2 i 1, spełniającym równość
2
� ś2
f(4) (x) a + b
f(x) =H3(x) + (x Ą a) x Ą (x Ą b) . (64)
4! 2
R Ą ó
b
a+b
Ale ponieważ (x Ą a) x Ą (x Ą b) dx =0, więc ze wzoru (63) wynikają związki Q(f) =I (H3) =I (L2).
a 2
4
Korzystając z równości (64) i tw. o wartości średniej (lemat 31) uzyskamy
R(f) =I(f) Ą Q(f) =I(f) Ą I (H3) =
Z � ś2
b
1 a + b
= f(4) (x)(x Ą a) x Ą (x Ą b) dx =
4! 2
a
� ś5
1 b Ą a
= Ą f(4) () , 2 (a; b) .
90 2
W poniższej tabeli podane zostaną współczynniki kwadratur Newtona-Cotesa oraz wartości stałych określających
resztę dla funkcji f 2 Cm [a; b], zgodnie z nastepującymi wzorami:
x(n) = xk = a + kh
k
n n
ł �
X X
(n)
Q(f) = Akf (xk) =Qn(f) =(b Ą a) Bk f x(n)
k
k=0 k=0
Rn(f) = Ąchjf(m) () , 2 (a; b) .
(n) (n)
i Bk = BnĄk Rn(f)
n 0 1 2 3 c j m
1 1
1 3 2
2 12
1 4 1
2 5 4
6 6 90
1 3 3
3 5 4
8 8 80
7 32 12 8
4 7 6
90 90 90 945
19 75 50 275
5 7 6
288 288 288 12096
41 216 27 272 9
6 9 8
840 840 840 840 1400
(n)
UWAGA: (i) Między współczynnikami Ak określonymi dla danego n wzorem (60) a współczynnikami Bk za-
chodzi następujący związek
Ak
(n)
Bk = .
b Ą a
(ii) Tylko dla n � 7 i n =9 wszystkie współczynniki kwadratur Newtona-Cotesa są dodatnie. Dla dowolnego n
n
X
(n)
(b Ą a) Bk = Qn(1) = I(1) = b Ą a ,
k=0
(n)
zatem suma współczynników Bk wynosi 1. J�sli dla dużych n część z nich jest ujemna, to suma modułów tych
współczynników jest większa od jedności. Przykładowo
10
Ż Ż
X
ŻB Ż
(10)
dla n = 10 Ż Ż ź 3:1
k
k=0
15
Ż Ż
X
ŻB Ż
(15)
dla n = 15 Ż Ż ź 8:3
k
k=0
20
Ż Ż
X
ŻB Ż
(20)
dla n = 20 Ż ź 560
k
k=0
i bardzo szybko rośnie, czyli nie jest spełnione założenie (b) twierdz
nia 27, co oznacza, że istnieją funkcje ciągłe, dla
których ciąg kwadratur Newtona-Cotesa fQn(f)g nie jest zbieżny do całki I(f).
(iii) Kwadratura dla n =3 jest nazywana kwadratrurą 3=8, zaś dla n =4 - kwadraturą Bodego.
5.2.2. Złożone kwadratury Newtona-Cotesa
Wzwiązku z poprzednią uwagą (ii) w praktyce nie stosuje się kwadratur Newtona-Cotesa wyższego rzędu. Kon-
struuje się natomiast kwadratury, nazywane złożonymi kwadraturami Newtona-Cotesa, dzieląc przedział całkowania
[a; b] na N równych podprzedziałów [xi; xi+1] długości h = (b Ą a)=N i stosując na każdym z nich kwadratury
Newtona-Cotesa niskiego rzędu.
5
PRZYKA
AD 3: Złożona kwadratura trapezów. Dla i =0; 1; :::; N Ą 1 mamy
Z
xi+1
h h3
f(x)dx = (f (xi) +f (xi+1)) Ą f00 (i) ,
2 12
xi
gdzie h =(b Ą a)=N oraz i 2 (xi; xi+1).
Stąd
� !
Z Z
NĄ1 NĄ1
b xk+1
X X
h
f(x)dx = f(x)dx = f(a) +2 f(a + kh) +f(b) Ą
2
a xk
k=0 k=1
X
h3 NĄ1
Ą f00 (i) .
12
k=0
Jeśli jf00(x)j �M dla x 2 [a; b], to
� ś
Ą ó
h3 h2 1
jR(f)j � NM = (b Ą a) M = O h2 = O .
12 12 N2
Zaś dla f 2 C2 [a; b]
NĄ1
X
N min f00 (k) � f00 (k) � N max f00 (k) ,
k k
k=0
a zatem
NĄ1
X
f00 (k) =Nf00 ()
k=0
dla pewnego 2 [a; b] i wówczas zachodzi równość
� !
Z
NĄ1
b
X
h h2
f(x)dx = f(a) +2 f(a + kh) +f(b) Ą (b Ą a)f00 () . (65)
2 12
a
k=1
UWAGA: Zauważmy, że złożona kwadratura trapezów (65) jest szczególnym przypadkiem kwadratury (59).
PRZYKA
AD 4: Złożona kwadratura parabol. Jeśli f 2 C4 [a; b], to
Z � � ś ś
xk+1
h xk + xk+1 1
f(x)dx = f (xk) +4f + f (xk) Ą h5f(4) (k) , k 2 (xk; xk+1) .
6 2 90
xk
Stąd
� !
Z � ś
NĄ1 N
b
X X
h h
f(x)dx = f(a) +2 f(a + kh) +4 f a +(2k Ą 1) + f(b) + O(h4) (66)
6 2
a
k=1 k=1
3
Innym przykładem złożonej kwadratury Newtona-Cotesa jest złożona kwadratura (dla n =3)
8
� !
Z PNĄ1
b
Ą ó
h
Ł Ą óń
f(x)dx = PNf(a) +2 i=1 f(aó+ ih)Ą+f(b)+ + O h4 (67)
8 +3 f a +(3i Ą 2)h + f a +(3i Ą 1)h
a
i=1 3 3
Można dowieść następującego twierdzenia:
TWIERDZENIE 32. Dla dowolnego n ciąg złożonych kwadratur Newtona-Cotesa QN(f) jest zbieżny dla wszyst-
kich funkcji f 2 C [a; b] przy N !1(czyli h ! 0)
Z
b
QN(f) ! f(x)dx:
a
DOWÓD pomijamy.
DEFINICJA 33. Mówimy, że kwadratura Q jest rzędu n, jeśli jest dokładna dla wszystkich wielomianów w stopnia
mniejszego od n, tzn. Q (w) =I (w) oraz istnieje wielomian wn stopnia n, dla którego Q (wn) = I (wn).
6
Ponieważ dla funkcji f będącymi wielomianami stopnia niewiększego od m zachodzi f(m+1) =0, zatem bezpośred-
nioz wzoru(57) mamynastępujący
LEMAT 34. Kwadratury interpolacyjne oparte na węzłach o łącznej krotności m +1sąrzędu co najmniej m +1.
6
Z przykładów 1 i 2 wynika, że kwadratury Newtona-Cotesa dla n = 1 są co najmniej rzędu 2, a dla n = 2 co
najmniej rzędu 4. Uogólniając, zachodzi następujące twierdzenie
TWIERDZENIE 35. Kwadratury Newtona-Cotesa oparte na n +1węzłach sąrzędu
�
n +2 dla n parzystych
n +1 dla n nieparzystych
DOWÓD pomijamy.
UWAGA: (i) Z wniosku 28 wynika zbieżność kwadratur interpolacyjnych o nieujemnych współczynnikach.
(ii) Załóżmy, że całkę I(f) przybliżamy dowolną kwadraturą Q(f) =Q(f; a; b) rzędu n. Dla funkcji f 2 Cn [a; b]
Ż
Żf(n)(x)Ż � Mn dla x 2 [a; b] mamy wtedy oszacowanie reszty postaci
Ż
takich, że
jI(f) Ą Q(f)j �MnKn (b Ą a)n+1 ,
gdzie stała Kn nie zależy od b Ą a.
bĄa
Dzieląc przedział całkowania na N równych podprzedziałów [xk; xk+1] długości xk+1 Ą xk = i przybliżając na
N
R
xk+1
każdym z nich całkę f(x)dx kwadraturą Q(f; xk; xk+1) konstruujemy kwadraturę złożoną
xk
NĄ1
X
QN(f) = Q (f; xk; xk+1) ,
k=0
której oszacowanie reszty ma postać
� śn
b Ą a
jI(f) Ą QN(f)j �MnKn(b Ą a) ; N =1; 2; 3; :::
N
Jeśli liczbę podprzedziałów będziemy kolejno np. podwajali (czyli N =2j), to
MnKn(b Ą a)n+1
jI(f) Ą Q2j (f)j � .
(2n)j
Przy założeniu, że stałe Mn i Kn są tego samego rzędu wielkości dla różnych n, wynika stąd, że wyższy rząd
kwadratury zapewnia lepszą dokładność.
(iii) W każdym kroku oblicza się wartości kolejnych kwadratur złożonych QN0, QN1, QN2, :::, za każdym razem
dzieląc poprzednie podprzedziałycałkowania na m części, czyli N0 =1; Nj = mj. Wielkość m należy tak dobrać, aby
przy wyznaczaniu nowej kwadratury można było skorzyst�c z poprzednio obliczonych wartości funkcji podcałkowej
(czyli aby punkty, które były węzłami dla Nj pozostaly nimi dla Nj+1). A
atwo stwierdzić, że m =2 dla kwadratury
trapezów i parabol, m =3 dla kwadratury 3=8, itd. Obliczenia wykonuje się tak długo, aż dwa kolejne przybliżenia
całki QNj i QNj+1 będą się różniłymniej niż o zadaną dokładność ".
(iv) Funkcje podcałkowe mają często zupełnie inny rząd wielkości w różnych podprzedziałach przedziału całkowania.
Wtedy w każdym podprzedziale stosuje się niezależnie kwadraturę i każdą z całek
Z Z Z Z
c1 c2 b b
+ +::: + =
a c1 ck a
oblicza się oddzielnie.
5.2.3. Przyspieszanie zbieżności ciągu kwadratur interpolacyjnych
Koncepcję analogiczną do ekstrapolacji Richardsona w przypadku różniczkowania numerycznego można zastosować
równieżprzycałkowaniu numerycznym, np. za pomocą kwadratur Newtona-Cotesa. Metodę przyspieszania zbieżności
kwadratur nazywamy metodą Romberga.
Oznaczmy przez R(n) złożoną kwadraturę trapezów przybliżającącałkę I(f) przy użyciu n podprzedziałów. Zatem
� !
� ś
nĄ1
X
b Ą a f(a) b Ą a f(b)
R(n) = + f a + k + :
n 2 n 2
k=1
Stosując wielokrotnie ekstrapolację (podobnie jak przy różniczkowaniu), otrzymamy metodę Romberga, w której kon-
struuje się trójkątną tablicę za pomocą następującego wzoru:
R[j; 0] = R(2j) (68)
4lR [j; l Ą 1] Ą R [j Ą 1; l Ą 1]
R [j; l] = (1)
4l Ą 1
7
Można sformułować:
TWIERDZENIE 36. Jeśli f 2 C [a; b], to każda kolumna w tablicy Romberga jest zbieżna do całki z funkcji f na
przedziale [a; b], tzn. dlakażdego l
Z
b
lim R [j; l] = f(x)dx:
j!1
a
DOWÓD wynika z wniosku 28, gdyż współczynniki dowolnej z kwadratur Romberga są dodatnie.
UWAGA: (i) A
atwo sprawdzić, że kwadratury R[j; 1] tworzące drugą kolumnę tablicy Romberga dla kwadratury
trapezów, są złozonymi kwadraturami parabol. Jednak elementy nastepnych kolumn nie są już innymi złożonymi
kwadraturami Newtona-Cotesa.
(ii) Dla innych kwadratur metoda Romberga wyrażają s/
e innymi wzorami. Jeśli ogólną post�c metody Romberga
zapiszemy jako:
alR [j; l Ą 1] Ą R [j Ą 1; l Ą 1]
R [j; l] = ; (69)
al Ą 1
to współczynnik al dla poszczególnych kwadratur wynosi:
4l+1 dla kwadratury parabol (a)
9l+1 dla kwadratury 3=8 (b)
16l+2 dla kwadratury Bodego (c)
25l+2 dla kwadratury przy n =5 (d)
36l+3 dla kwadratury przy n =6 (e)
PRZYKA
AD: Jaka wygląda złożona formułatrapezówdlan =1; 2; 4; 8, gdy przedziałem całkowania jest [0; 1] ?
1 1
R(1) = f(0) + f(1)
2 2
� ś
1 1 1 1
R(2) = f(0) + f + f(1)
4 2 2 4
� � ś � ś � śś
1 1 1 1 3 1
R(4) = f(0) + f + f + f + f(1)
8 4 4 2 4 8
� � ś � ś � ś � ś � ś � ś � śś
1 1 1 1 3 1 5 3 7 1
R(8) = f(0) + f + f + f + f + f + f + f + f(1)
16 8 8 4 8 2 8 4 8 16
Zpowyższego przykładu widać, że aby obliczyć R(2n), można wykorzystać obliczone R(n). Wówczas trzeba będzie
tylko obliczyć składniki nieobecne w R(2n), tzn. wartości funkcji w węzłach, które nie były używane w R(n). Mi-
anowicie:
� ś
1 1 1
R(2) = R(1) + f
2 2 2
� � ś � śś
1 1 1 3
R(4) = R(2) + f + f
2 4 4 4
� � ś � ś � ś � śś
1 1 1 3 5 7
R(8) = R(4) + f + f + f + f
2 8 8 8 8 8
Ogólna formuła na dowolnym przedziale [a; b] jest następująca:
n
X
1
R(2n) = R(n) +h f (a +(2k Ą 1) h) ;
2
k=1
bĄa
gdzie h = .
2n
Przyjmijmy oznaczenia:
1
R [0; 0] = (b Ą a)(f(a) +f(b)) ;
2
� ś
2nĄ1
X
1 b Ą a b Ą a
R [n; 0] = R [n Ą 1; 0] + f a +(2i Ą 1) ;
2 2n 2n
i=1
8
czyli kolejne R [0; 0], R [1; 0], R [2; 0], ... oznaczająwartości rekurencyjnej kwadratury trapezów, które zostaną znalezione
unikając powtórnego obliczania tych samych wartości funkcji. Zatem dodatkowo, oprócz stosowania ekstrapolacji,
obliczanie niektórych kwadratur Newtona-Cotesa można łatwo przyspieszyć stosując powyższą koncepcję.
UWAGA: Oczywiście metodę Romberga dla dowolnej kwadratury można stosować bez rekurencji. Nie powoduje
to jedynie skrócenia czasu obliczeń.
5.3. Obliczanie całek niewłaściwych i całek osobliwych
W tej części zajmiemy się pewnymi niezależnymi od stosowanej kwadratury metodami całkowania, w sytuacji gdy
przedział całkowania nie jest skończony (całki niewłaściwe) lub funkcje podcałkowe nie są dostatecznie regularne (całki
osobliwe).
Zacznijmy od przypomnienia, że dla skończonego a całkę niewłaściwą w przedziale [a; +1) de
�niuje się jako
granicę całek oznaczonych w [a; b] przy b ! +1, tzn.
Z Z
+1 b
f(x)dx = lim f(x)dx: (70)
b!+1
a a
Załóżmy, że funkcja f jest określona w przedziale skończonym [a; b] oraz dla dowolnego c 2 (a; b) jest całkowalna
w [c; b] i nieograniczona w [a; c]. Punkt a nazywamy wtedy punktem osobliwym funkcji f, a przez całkę w [a; b]
rozumiemy granicę całek oznaczonych w [c; b] przy c ! a+, tzn.
Z Z
b b
f(x)dx = lim f(x)dx: (71)
c!a+ c
a
Analogicznie de
�niuje się całkę niewłaściwą w przedziale (Ą1; a] oraz całkę osobliwą w [a; b] z punktem osobliwym
b. W dalszych rozważaniach zakładamy, że powyższe granice istnieją i są skończone.
UWAGA: Założenie, że punktem osobliwym jest kraniec przedziału całkowania nie zmniejsza ogólności, gdyż
zawsze możemy przejść do sumy całek na podprzedziałach, w których punktami osobliwymi będą krańce.
Podstawowe metody numerycznego obliczania całek niewłaściwych i osobliwych można podzielić na następujące
rodzaje:
(a) Dokonanie zamiany zmiennej, przekształecenie funkcji podcałkowej, całkowanie przez części, itp. sposoby znane z
analizy, które pozwolą na przejście do całki na skończonym przedziale lub całki bez osobliwości.
(b) Korzystanie wprost z de
�nicji (70) lub (71) za pomocą rozbicia przedziałucałkowania
Z Z Z
+1 b +1
f(x)dx = f(x)dx + f(x)dx = I1(f) +I2(f)
a a b
Z Z Z
b b c
f(x)dx = f(x)dx + f(x)dx = I1(f) +I2(f)
a c a
tak, by składnik I2(f) dałosię obliczyć w specjalny sposób albo był dostatecznie mały, by można go zaniedbać. Innymi
słowy: należy znalez
ć dostatecznie dużą liczbę b w pierwszym przypadku lub liczbę c dostatecznie bliską a wdrugim
"
przypadku, by zachodziła nierównósć jI2(f)j < . Wówczas, jeśli całkę I1(f) przybliżymy za pomocą dowolnej
2
"
kwadratury z dokładnością , to otrzymamy, że zadaną dokładność ".
2
(c) Zastosowanie następującego algorytmu: dla całki osobliwej podstawiamy
Z
aj
b Ą a
aj = a + ; j =1; 2; ::: oraz Ij := f(x)dx;
j
aj+1
zaś dla całki niewłaściwej
Z
aj+1
Ij := f(x)dx;
aj
gdzie fajg jest pewnym ciągiem rozbieżnym i obliczamy I = I1 + I2 + :::.
(d) Przedstawienie funkcji f w postaci iloczynu p ó g, aby funkcja p mogła bżc traktowana jako funkcja wagowa,
a w przypadku całek osobliwych zawierała również osobliwósci funkcji podcałkowej i zastosowanie odpowiedniej
kwadratury Gaussa dla całki funkcji g z wagą p (patrz dalsze części wykładu).
PRZYKA
AD 1: Połączenie metod (a) i (b) dla całki osobliwej
Z ź
2
ln (sin x) dx
0
9
Ąsin xó
Korzystamy z tożsamości ln (sin x) =ln x +ln i otrzymujemy
x
Z ź Z ź Z ź � ś
2 2 2
sin x
ln (sin x) dx = ln xdx + ln dx = I1 + I2:
x
0 0 0
Całka I1 jest nadal osobliwa w x =0, ale łatwo ją obliczyć z de
�nicji (71)
Żź
Z ź Z ź
2
2 2 Ż
ź ź ź
ln xdx = lim ln xdx = lim (x ln x Ą x)Ż = ln Ą ;
Ż
c!0+ c c!0+ 2 2 2
0
c
zaś numeryczne obliczenie całki I2 za pomocą dowolnej kwadratury nie jest już problemem.
PRZYKA
AD 2: Różne użyteczne transformacje całek osobliwych i niewłaściwych
Z Z
1 ź
f(x)
p dx = f (cos t) dt
1 Ą x2
Ą1 0
Z Z
1 1
1
xĄ n f(x)dx = n tnĄ2f (tn) dt dla f 2 C [0; 1] i n � 2
0 0
Z Z ź
1
2
Ą ó
f(x)
p dx = 2 f sin2 t dt
x(1 Ą x)
0 0
PRZYKA
AD 3: Całkowanie przez części dla całki, która nie wydaje się być kłopotliwa
Z
1
p
x + x3dx
0
0
gdyż funkcja podcałkowa jest ciągła w [0; 1], ale f posiada osobliwość w x =0, czyli mamy tylko f 2 C [0; 1], a to
zbyt mało np. przy stosowaniu kwadratur Newtona-Cotesa. A zatem
Ż1 1 5
Z Z Z
1 1
p p p
Ż
1 2 3 2 x2
x + x3dx = x2 1+x2dx = x2 1+x2Ż Ą p dx
Ż
3 3
1+x2
0 0 0
0
i otrzymujemy całkę, dla której f 2 C2 [0; 1], więc można ją swobodnie przybliżać kwadraturą trapezów.
Czasami trzeba zastosowaćcałkowanie przez części kilkakrotnie, zanim zastosuje sięmetodę numeryczną, np. dla całek
Z
1
ex
p :
x
0
Po pierwszym calkowaniu przez części otrzymujemy funkcję podcałkową tylko klasy C [0; 1]
PRZYKA
AD 4: Całki podwójnie niewłaściwe, np.
Z
+1
2
eĄx dx,
Ą1
całkuje się numerycznie na przedziale [ĄR1; R2] tak wybranym, aby funkcja podcałkowa i jej pochodne niskich rzędów
były dostatecznie małe dla x �ĄR1 i x � R2. Okazuje s/
e, że nawet wzór trapezów daje zaskakująco dobrą dokład-
ność przy obc/
eciu przedziału całkowania do [Ą4; 4]. Spowodowane jest to faktem, że wartości funkcji podcałkowej
1
dla x = ż4 są mniejsze od Ł 10Ą6 (1: 125 4 Ł 10Ą7). Oszacujmy bład obcięcia przedziałucałkowania
2
Z Z Z
+1 +1 +1
2 1 1 1 1 1
jRj =2 eĄx dx =2 eĄt tĄ 2 dt < 2 ó ó 16Ą 2 eĄtdt = eĄ16 =2: 813 4 Ł 10Ą8:
2 2 4
4 16 16
10
5.4. Kwadratury o ustalonej liczbie węzłów
W tej części zbadamy, jaki jest maksymalny rząd kwadratury opartej na ustalonej liczbie węzłów. Najpierw jednak
omówimy niezbędne własności wielomianów ortogonalnych, które bedą nam potrzebne dalej.
5.4.1. Wielomiany ortogonalne
DEFINICJA 37. Ciąg P0; P1; :::; Pn, gdzie Pi; i = 0; 1; :::; n (n � 1), jest wielomianem stopnia dokładnie i,
nazywamy ciągiem wielomianów ortogonalnych na przedziale [a; b] z wagą p, jeśli tworzą one układ ortogonalny w
przestrzeni L2 [a; b], tzn.
p
Z
b
hPk; Pli = Pk(x)Pl(x)p(x)dx =0 dla k = l:
6
a
Jeżeli, prócz tego hPi; Pii =1 dla każdego i, tociąg taki nazywamy ortonormalnym.
UWAGA: Dla konstrukcji wielomianów ortogonalnych wystarczy zortogonalizować metodą Grama-Schmidta np.
ciag 1; x; :::; xn.
TWIERDZENIE 38. (Regułatrójczłonowa) Wielomiany ortogonalne Pi; i =0; 1; :::; n spełniają zależność rekuren-
cyjną
Pi(x) = (�ix + Żi) PiĄ1(x) +�iPiĄ2(x); i =1; 2; :::; n (72)
PĄ1(x) = 0; Po(x) =a0
ai
gdzie �i = =0, Żi oraz �i =0 dane sąwzorami
6 6
aiĄ1
hxPiĄ1; PiĄ1i hPiĄ1; PiĄ1i
Żi = Ą oraz �i = Ą ;
hPiĄ1; PiĄ1i hPiĄ2; PiĄ2i
ai zaś oznacza współczynnik wielomianu Pi przy najwyższej potędze.
DOWÓD pomijamy.
Dalsze własności (podane bez dowodów) są ważne dla ich zastosowańwcałkowaniu nuemrycznym.
LEMAT 39. Dowolny wielomian Qj stopnia j niższego od i jest ortogonalny do i-tego wielomianu ortogonalnego
Pi, tzn. hQj; Pii =0 dla j TWIERDZENIE 40. Pierwiastki wielomianów ortogonalnych w przestrzeni L2 [a; b] są rzeczywiste, jednokrotne i
p
leżąwewnątrz przedziału [a; b].
Przedstawimy teraz kilka przykładów ciągów wielomianów ortogonalnych na różnych przedziałach i z różnymi
funkcjami wagowymi.
PRZYKA
AD 1: Wielomiany Legendre a postaci:
1 di Ą ói
P0(x) =1; Pi(x) = x2 Ą 1 dla i =1; 2; :::
2ii! dxi
są ortogonalne na przedziale [Ą1; 1] z funkcją wagową p(x) =1.
Współczynniki we wzorze (72) są równe:
2i Ą 1 i Ą 1 (2i)!
�i = ; =0; �i = Ą ; ai = :
%2ł
i i 2i(i!)2
Ponadto
Z
1
2
kPik2 = (Pi(x))2 dx = .
2i +1
Ą1
PRZYKA
AD 2: Wielomiany Hermite a postaci:
2 2
di
H0(x) =1; Hi(x) =(Ą1)iex eĄx
dxi
2
są ortogonalne na przedziale (Ą1; +1) z funkcja wagową p(x) =eĄx .
Współczynniki we wzorze (72) są równe:
�i =2; =0; �i = Ą(2i Ą 2); ai =2i:
%2ł
Ponadto
Z
+1
2 p
kHik2 = eĄx (Hi(x))2 dx = ź2ii!.
Ą1
11
PRZYKA
AD 3: Wielomiany Laguerre a postaci:
di Ą ó
L0(x) =1; Li(x) =(Ą1)iex xieĄx
dxi
są ortogonalne na przedziale [0; +1) z funkcją wagową p(x) =eĄx.
Współczynniki we wzorze (72) są równe:
i Ą 1
�i = Ą1; =0; �i = Ą ; ai =(Ą1)i:
%2ł
i
Ponadto
kLik2 =(i!)2 .
PRZYKA
AD 4: Wielomiany Czebyszewa postaci:
Ą p ói Ą p ói
x + x2 Ą 1 + x Ą x2 Ą 1
Ti(x) = dla i =0; 1; :::
2
1
p
są ortogonalne na przedziale [Ą1; 1] z funkcją wagową p(x) = .
1Ąx2
Współczynniki we wzorze (72) są równe:
�i =2; =0; �i = Ą1; ai =2iĄ1:
%2ł
Ponadto
ź
kTik2 = dla i =0:
6
2
5.4.2. Kwadratury Gaussa
Niech Q będzie kwadraturą opartą na dowolnyP węzłach x0; x1; :::; xn o dowolnych krotnościach �0 +1; �1 +
ch
n
1; :::; �n +1, których suma jest ustalona, tzn. n + �i = m, tzn.
i=o
�0 �1 �n
X X X
Q(f) = Ak;0f(k) (x0) + Ak;1f(k) (x1) +::: + Ak;nf(k) (xn) : (73)
k=0 k=0 k=0
W klasie takich kwadratur szukamy kwadratury o maksymalnym rzędzie przybliżającej całkę Ip(f).
Z lematu 34 wiemy, że rząd kwadratury interpolacyjnej opartej na węzłach o łącznej krotności m+1 może byćrówny
co najmniej m +1. Można się spodziewać, że przy odpowiednio wybranych węzłach (niekoniecznie równoodległych)
maksymalny rząd kwadratury będzie istotnie większy od m +1. Okazuje s/
e, że zachodzi następujące:
Pk
TWIERDZENIE 41. Jeśli kwadratura Q postaci (73) jest rzędu r, gdzie r > �i = m +1, to jest ona
i=0
kwadraturą interpolacyjną.
DOWÓD: Niech Hm będzie wielomianem interpolacyjnym Hermite a funkcji f opartym na węzłach xi okrotnoś-
(j
ciach �i; i =0; 1; :::; n. Zatem Hm)(xi) =f(j)(xi) dla j =0; 1; :::; �i; i =0; 1; :::; n, a stąd otrzymujemy równość
Q(f) =Q(Hm). Ale z założenia kwadratura Q jest rzędu r >m+1, czyli jest dokładna dla wielomianu Hm stopnia
co najwyżej m, tzn. Q (Hm) =Ip (Hm). Tymsamym Q(f) =Ip (Hm), co kończy dowód. Ą
Z powŹzszego twierdzenia wynika, że szukając kwadratury o maksymalnym rzędzie można ograniczyć s/
e do
kwadratur interpolacyjnych. Niech Q(f) = Ip(Hm) będzie dowolną kwadraturą interpolacyjną opartą na węzłach
o łacznej krotności m +1. Korzystając z postaci reszty interpolacji Hermite a można zapisać
f(x) =Hm(x) +!� (x)g(x)
m+1
(patrz wzór (57)), a więc węzłysą pierwiastkami wielomianu !� . Całkując obustronnie, otrzymamy
m+1
Ip(f) =Q(Hm) +R(f);
gdzie reszta kwadratury
Z
b
R(f) = p(x)!� (x)g(x)dx
m+1
a
� �
może być interpretowana jako iloczyn skalarny !� ; g w przestrzeni L2 [a; b].
m+1 p
Jeśli f jest wielomianem stopnia l, to g jest również wielomianem i to stopnia l Ą m Ą 1. Znalezienie kwadratury
o maksymalnym rzędzie, czyli dokładnej dla wielomianów możliwie wysokiego stopnia, sprowadza się zatem do wyz-
naczenia wielomianu !� ortogonalnego do wielomianów jak najwyższego stopnia. Z lem.39 wiemy, że rozwiązaniem
m+1
tego zadania jest m +1-szy wielomian ortogonalny na przedziale [a; b] z wagą p. Jest on prostopadły do wszystkich
12
� �
wielomianów g stopnia mniejszego od m +1 (tzn. !� ; g =0), a stąd l Ą m Ą 1 m+1
Ą ó2
Biorąc zaś wielomian !� stopnia 2m +2 dostajemy
m+1
łĄ ó2� łĄ ó2�
0 =Q !� = Ip !� > 0.
6
m+1 m+1
Tym samym udowodniliśmy:
TWIERDZENIE 42. Kwadraturą postaci (73) o maksymalnym rzędzie (równym 2m +2) jest kwadartura interpo-
lacyjna, której węzłami sąpierwiastki m +1-szego wielomianu ortogonalnego na przedziale [a; b] z wagą p.
Kwadratury takie są nazywane kwadraturami Gaussa. Ich węzły są rzeczywiste, jednokrotne i leżą wewnątrz
przedziałucałkowania [a; b] (tw.40). Kwadratury Gaussa są zatem postaci (55), a ich współczynniki można wyznaczać
z ogólnego wzoru (A) dla kwadratur interpolacyjnych. Można udowodnić (korzystając z formuły Christofela-Darboux),
że
am+1 kPmk2
Ak = : (74)
0
am Pm+1(xk)Pm(xk)
Można również pokazać, że
1
Ak =
Pm ł ~ �2 ;
Pi(xk)
i=0
~
gdzie Pi sa wielomianami ortonormalnymi, skąd wynika, że współczynniki kwadratur Gaussa są dodatnie. Na mocy
wn.28 prawdziwe jest zatem następujące:
TWIERDZENIE 43. Dla dowolnego przedziału [a; b] i funkcji wagowej p ciąg kwadratur Gaussa jest zbieżny do
całki Ip(f) dla wszystkich funkcji f ciągłych na [a; b].
Dla przedstawienia reszty kwadratury Gaussa można skorzystać z wzoru (57), otrzymamy wówczas:
LEMAT 44. Jeśli f 2 C2n+2 [a; b], toresztakwadratury Gaussajest równa
Z
b
Ą ó2
f(2n+2) ()
R(f) = p(x) !� dx; 2 (a; b) : (75)
n+1
(2n +2)!
a
UWAGA: (i) Oczywiście, aby móc stosować kwadratury Gaussa musimy znać pierwiastki odpowiednich wielomi-
anów ortogonalnych. Dla wielu funkcji wagowych p wartości węzłów i współczynników kwadratur Gaussa sa stabli-
cowane. Z drugiej strony koszt wyznaczenia samych węzłów może byćtakduży, że stosowanie kwadratur Gaussa może
stać się nieopłacalne.
(ii) W przypadku funkcji wagowej symetrycznej względem środka przedziałucałkowania, węzły kwadratury Gaussa są
równiez symetryczne oraz zachodzi równość współczynników Ak = AmĄk.
(iii) Podobnie, jak w przypadku kwadratur Newtona-Cotesa można konstruwać złożone kwadratury Gaussa, dzieląc
przedział całkowania na podprzedziały i stosując w każdym z nich kwadraturęR ustalonego stopnia.
Gaussa
1
PRZYKA
AD 1: Skonstruujmy kwadraturę Gaussa z 2 węzłami dla całki f(x)dx. Ponieważ p(x) = 1, w/
ec
Ą1
1 1 1
p
odpowiednim wielomianem ortogonalnym jest wielomian Legendre a P2(x) = (3x2 Ą1), skąd x0 = Ąp3, x1 = .
2
3
Współczynniki można łatwo wyznaczyć stosując kwadraturę
Z
1
f(x)dx = A0f(x0) +A1f(x1)
Ą1
do funkcji f(x) =1 i f(x) =x, wówczas otrzymamy układ równań
A0 + A1 = 2
1 1
p p
Ą A0 + A1 = 0
3 3
którego rozwiązaniami są A0 = A1 =1. Ostatecznie otrzymujemy kwadraturę Gaussa-Legendre a postaci
Z � ś � ś
1
1 1
p p
f(x)dx ź f Ą + f ;
3 3
Ą1
1 1
która jest rzędu 4. Np. eĄ p3 + ep3 = 2: 342 7, prosta kwadratura Simpsona (przy użyciu 3 wartości funkcji pod-
Ą ó
całkowej) daje eĄ1 +4e0 + e1 =3 = 2: 362 1, zaś dokładny wynik I =2 sinh (1) =2: 350 4.
13
PRZYKA
AD 2. Obliczmy całkę
Z
1
g(x)
p dx
1 Ą x2
Ą1
Ąp óĄ1
z osobliwościami w obu krańcach przedziału całkowania. Można je usunąć traktując funkcję p(x) = 1 Ą x2
jako wagę. Naturalne będzie przybliżenie całki za pomocą kwadratury Gaussa-Czebyszewa, której węzłami są xi =
2i+1 ź
cos ź, zas współczynniki Ai = . Zatem
2n+2 n+1
� ś
n
X
ź 2k +1
Q(f) = g cos ź
n +1k=0 2n +2
oraz
ź g(2n+2)()
R(f) = ; 2 (Ą1; 1) :
22n+1 (2n +2)!
UWAGA: Węzły i współczynniki kwadratur Gaussa-Laguerre a i Gaussa-Hermite a są stablicowane i mozna je
znalez
ć np. w [FMW].
LEMAT 45. Przy liniowej zamianie zmiennej x = �t + Ż, gdzie
b Ą a ad Ą bc
� = ; Ż =
d Ą c d Ą c
Pn
kwadratura Q(f) = Aif(xi) spełniająca równanie
i=0
Z
b
p(x)f(x)dx = Q(f) +R(f) (76)
a
Pn
~
przekształca się wkwadraturę Q(g) = Big(ti), spełniająca zależność
i=0
Z
d
~ ~
v(t)g(t)dt = Q(g) + R(g); (77)
c
gdzie
1 1
ti = (xi Ą Ż) ; Bi = Ai dla i =0; 1; :::; n
� �
1
~
oraz R(g) = R(f).
�
~
UWAGA: Jeśli kwadratura Q jest rzędu n, to kwadratura Q, otrzymana przez liniową zamianę zmiennej, jest również
rzędu n.
R
3
PRZYKA
AD 3: Skonstruujmy kwadraturę Gaussa opartą na 4 węzłach dla całki g(t)dt, a więc znów z funkcją
2
wagową p(t) =1. Użyjemy wielomianów ortogonalnych Legendre a, a nastepnie dokonamy liniowej zamiany zmien-
nych. Można obliczyć
5 3 35 30 3
P3(x) = x3 Ą x oraz P4(x) = x4 Ą x2 + :
2 2 8 8 8
Węzłami kwadratury będą
s
p
15 + 2 30
Ąx0 = x3 =
35
s
p
15 Ą 2 30
Ąx1 = x2 =
35
Z wzoru (74) można wyznaczyć współczynniki
35 2
8 7
Ai = ;
5 0
P4(xi)P3(xi)
2
14
skąd
49
A0 = A3 = Ą p ó
6 18 + 30
49
A1 = A2 = Ą p ó
6 18 Ą 30
W ten sposób budujemy kwadraturę pomocniczą
Z
3
1
X
f(x)dx = Aif(xi) +R(f):
Ą1
iĄ0
Ą1 ó
1
Następnie podstawiamy t = (x +5) i przyjmując f(x) =g (x +5) otrzymujemy szukaną kwadraturę Gaussa-
2 2
Legendre a, dla której
Z
3
3
X
~
g(t)dt = Big(ti) + R(g);
2
i=0
gdzie
1 1 1
~
ti = (xi +5); Bi = Ai oraz R(g) = R(f):
2 2 2
5.5. Błędy obliczeń zmiennopozycyjnych przy obliczaniu kwadratur
Załóżmy, że dla wszystkich węzłów xi wartości funkcji są obliczane (lub dane) z błędem nie większym niż ąi, tzn.
jeśli yi są znanymi przybliżeniami nieznanych dokładnych wartości yi = f (xi), to
~
yi � yi + ąi oraz jąij �ą dla i =0; 1; :::n:
~
Wynika stąd następujące:
n n n
X X X
~
~
Q(f) Ą Q(f) = Aiyi (1 + ei) Ą Aiyi � [Ai (yi + ąi)(1+ei) Ą Aiyi] =
i=0 i=0 i=0
n n
X X
ź Aiąi + Aiyiei = E1(f) +E2(f);
i=0 i=0
gdzie jeij � 2Ąt (n Ą i +3), a t oznacza liczbę cyfr mantysy (wielkość ta wynika ze sposobu obliczania wartości
kwadratury - kolejno dosumowuje się iloczyny) - patrz [JD].
Błąd E1(f) związany z przybliżonym obliczaniem wartości funkcji ma oszacowanie
n
X
jE1(f)j �ą jAij ;
i=0
Wynika stąd, że ograniczoność sumy wartósci bezwzględnych współczynników kwadratury jest istotna nie tylko ze
względu na zbieżność kwadratur (tw.27), ale również ze względu na wielkość błędu wytworzonego przy obliczaniu
kwadratury w arytmetyce zmiennopozycyjnej.
Błąd E2(f) spełnia nierówność
n
X
jE2(f)j �2Ąt (n Ą i +3) jAiyij :
i=0
Oznacza to, że wzrastająca liczba węzłów n powoduje wzrost wytworzonego błęduconajwyżej proporcjonalnie do n:
p
W rzeczywistości wzrost ten jest raczej rzędu n, a praktycznie jego wpływ mozna zauważyć jedynie przy kwadrat-
urach o bardzo dużej liczbie węzłów - można stosować s/
e wtedy algorytm Gilla-M�llera sumowania z poprawkami.
Wowczas
n
X
jE2(f)j �3 ó 2Ąt jAiyij :
i=0
15
Podsumowując, nawet jeśli reszta kwadratury R(f) dąży do zera przy n !1, to trzeba się liczyć z tym, że błędy
operacji zmiennopozycyjnych nie pozwolą na uzyskanie dokładności wiekszej niż
n n
X X
jAiąf(xi)j +2Ąt jAif(xi)j :
i=0 i=0
Jest to istotne ze względu na zadawaną żadaną dokładność w procedurach numerycznego całkowania.
UWAGA: Nie można ogólnie stwierdzić, że w konkretnych przypadkach mniej dokładny kwadratura przy tym
samym doborze węzłów nie może dać lepszych wyników. Patrz poniższe przykłady.
PRZYKA
AD 1. Porównajmy dokładnośc różnych kwadratur opartych na 3 węzłach na przykładzie obliczenia całki
Z
1
p p
2
I = 2+xdx =2 3 Ą ź 2: 797 434 948:
3
Ą1
Stosując złożony wzór trapezów otrzymamy
Ąp p p ó
1 1
I ź 2 Ą 1+2 2+0+ 2+1 = ó 5: 560 477 932 = 2: 780 238 966;
2 2
zastosowanie wzoru parabol daje wynik
Ąp p p ó
1 1
I ź 2 Ą 1+4 2+0+ 2+1 = ó 8: 388 905 057 = 2: 796 301 686;
3 3
a według kwadratury Gaussa-Legendre a obliczymy
Ąp p ó p
I ź 0:555566 2 Ą 0:774597 + 2+0:774597 +0:888889 2+0=2: 797 491 94:
Zatem w tym przypadku najbardziej dokadny wynik daje kwadratura Gaussa-Legendre a.
PRZYKA
AD 2. Całka funkcji f(x) =Ą8+45x2 Ą 25x4 na [Ą1; 1] wynosi I =4. Dla h =1 wzór trapezów daje
wartość dokładną
1 1
I1 = f(Ą1) + f(0) + f(1) = 4;
2 2
zaś wzór parabol daje nawet złyznakcałki
1 8
I2 = [f(Ą1) + 4f(0) + f(1)] = Ą :
3 3
Dokładność kwadratur z ustalona liczbą węzłów zależy w sposób istotny od ich rozmieszczenia. Ogólnie, gdy funkcja
podcalkowa f(x) ma w przedziale całkowania dużą liczbę miejsc zerowych lub ekstremów, to dokładność kwadratur
znacznie się zmniejsza (ze względu na to, że pochodne wyższych rzędów funkcji podcalkowej przyjmują wówczas duże
wartości). Powinno sie wtedy wybierać h o wiele mniejsze niż najmniejsza odległość pomiędzy sąsiednimi mijescami
zerowymi funkcji f i jej pochodnej f0.
5.6. Całki wielokrotne
Numeryczne obliczanie całek wielokrotnych postaci
Z Z Z
In = ::: f (x1; x2; :::; xn) dx1dx2:::dxn;
�
gdzie � jest ustalonym obszarem w przestrzeni Rn, jest znacznie trudniejsze niż całek jednokrotnych, gdyż obszarów
różnego typu jest wówczas nieskończenie wiele (w przeciwieństwie do R1, gdzie mamy tylko trzy podstawowe typy
przedziałów całkowania, tj. [a; b], [a; 1) oraz (Ą1; +1)).
W najprostszym przypadku, w którym całka wielokrotna daje się wyrazić jakocałka iterowana
Z Z Z
b1 b2(x1) bn(x1;x2;:::;xnĄ1)
In = dx1 dx2::: f (x1; x2; :::; xn) dxn;
a1 a2(x1) an(x1;x2;:::;xnĄ1)
n-krotnie stosuje się kwadratury jednowymiarowe.
R R
PRZYKA
AD: Rozpatrzmy całkę podwójną f(x; y)dxdy, gdzie
�
� = f(x; y) : a1 � x � b1; a2(x) � y � b2(x)g :
16
Jeśli funkcje f(x; y), a2(x) i b2(x) sa ciągłe, to
Z Z Z Z Z
b1 b2(x) b1
f(x; y)dxdy = dx f(x; y)dy = g(x)dx;
a1 a2(x) a1
�
gdzie
Z
b2(x)
g(x) = f(x; y)dy:
a2(x)
R
Pn1
b1
Całkę g(x)dx możemy przybliżać dowolną z omówionych kwadratur Q1(g) = Aig(xi). Potrzebne wartości
a1 i=0
R
Pn2
b2(xi)
g(xi) = f(xi; y)dy mogąbyć aproksymowane inną(lubtąsamą) kwadraturą postaci Q2(') = Bj'(yj).
a2(xi) j=0
Dla każedgo xi przedział całkowania [a2(xi); b2(xi)] jest na ogół różny, a więc węzły yj zależą od i, tzn. yj = yij.
Oznaczając fi(y) =f(xi; y) otrzymujemy
n2
X
Q2(fi) = Bjf(xi; yij):
j=0
I dalej
Z Z
n1 n2
X X
f(x; y)dxdy = Ai Bjf(xi; yij) +R(f);
i=0 j=0
�
gdzie
n2
X
R(f) =R1(g) + AiR2(fi);
i=0
a R1 i R2 są odpowiednio resztami kwadratur Q1 i Q2.
UWAGA: W przypadku wielowymiarowym nie ma możliwości dowolnego wyboru węzłów, co stanowi istotne
utrudnienie. Poza tym wraz ze wzrostem wymiaru przestrzeni wystepuje szybki wzrost liczby węzłów potrzebnych
do uzyskania nawet niewielkiej dokładności przybliżenia całki. Niestety już wielomiany ortogonalne dwóch zmien-
nych mają pierwiastki zespolone lub rzeczywiste leżace poza obszarem całkowania, więc nie można skonstruować
wielowymiarowych kwadratur Gaussa dla dowolnego obszaru. W takim przypadku dzieli się obszar całkowania na
mniejsze podstawowe (n-wymiarowy prostopadłościan, kula lub sympleks) albo stosuje zamianę zmiennych.
17