Kwadratury numeryczne
Kwadraturami numerycznymi nazywamy wzory służące do przybliżonego wyznaczania wartości całek oznaczonych w obszarze jedno lub wielo wymiarowym. Obliczanie całek oznaczonych w ob- szarze wielowymiarowym sprowadza się do wielokrotnego zastosowania kwadratur dla obszaru jednowymiarowego. Idea postępowania przy całkowaniu numerycznym jest zawsze taka sama, a mianowicie zastępujemy funkcję podcałkową funkcją łatwą do scałkowania analitycznego (wielo- mian) w drodze interpolacji, a następnie funkcję interpolującą całkujemy ściśle (analitycznie). W rezultacie otrzymujemy wzór numerycznego całkowania.
W zależności od sposobu postępowania przy wyborze położeń węzłów interpolacji w przedziale całkowania możemy mieć do czynienia z kwadraturami z:
• zamkniętymi końcami (końce przedziału całkowania wchodzą do wzorów), a węzły są roz- mieszczone równomiernie w całym przedziale całkowania,
• otwartymi końcami (końce przedziału całkowania nie wchodzą do wzorów), a węzły są rozmieszczone w przedziale całkowania nierównomiernie, tak aby zminimalizować błąd kwadratury.
Przykładem kwadratur pierwszego typu są kwadratury Newtona-Cotesa, a kwadratur drugiego typu kwadratury Gaussa.
Na początek rozważmy wyprowadzenie kilku pierwszych kwadratur typu Newtona. Zgodnie z uwagami zapisanymi powyżej, operację wyznaczania całki:
( )
xdx Fb
∫
a , (1)przeprowadzimy dwuetapowo. Najpierw dokonamy interpolacji funkcji F
( )
x wielomianem La- grange’a stosownego stopnia (oznaczonego dalej jako n ), a następnie funkcję interpolującą scałku- jemy. Wzory ogólne w takim przypadku przyjmują następującą postać:( ) ( )
( )( ) ( ) ( )
( )( )
∫ ∑ ∫
∑
=
=
⋅
⋅
≅
⋅
⋅
≅
b a
n j
b a
n j j
n j n
j j
dx x L x F dx x F
x L x F x
F
0
0 . (2)
Po analitycznym wykonaniu całkowania (2)2, dla konkretnej wartości n , pamiętając o równomier- nym rozmieszczeniu węzłów, otrzymamy ogólny wzór całkowania numerycznego całkujący ściśle wielomiany stopnia n . Podamy teraz kilka pierwszych kwadratur typu Newtona:
• interpolacja wielomianowa stopnia 0 (n=0), czyli funkcją stałą:
( ) ( )
( )( )
( ) ( )
( )( ) ( ) ( )
∫ ∑ ∫
∑
− ⋅
=
⋅
⋅
≅
⋅
⋅
≅
→
=
=
= b
a j
b a
j j
j
j j
x a F dx b
x L x F dx x F
x L x F x
F n
0 0
0
0 0 0
0
1 0
, (3)
• interpolacja wielomianowa stopnia 1 (n=1), czyli funkcją liniową:
( ) ( )
( )( )
( ) ( )
( )( ) ( ) ( ) ( ) [ ]
∫ ∑ ∫
∑
+
− ⋅
=
⋅
⋅
≅
⋅
⋅
≅
→
=
=
= b
a j
b a
j j
j j j
x F x a F dx b
x L x F dx x F
x L x F x
F n
1 0
1
0
1 1 1
0
2 1
, (4)
• interpolacja wielomianowa stopnia 2 (n=2), czyli funkcją kwadratową:
( ) ( )
( )( )
( ) ( )
( )( ) ( ) ( ) [ ( ) ( ) ]
∫ ∑ ∫
∑
+
⋅ +
− ⋅
=
⋅
⋅
≅
⋅
⋅
≅
→
=
=
= b
a j
b a
j j
j j j
x F x F x
a F dx b
x L x F dx x F
x L x F x
F n
2 1
0 2
0
2
2 2 0
6 4 2
, (5)
• interpolacja wielomianowa stopnia 3 (n=3), czyli funkcją trzeciego stopnia:
( ) ( )
( )( )
( ) ( )
( )( ) ( ) ( ) [ ( ) ( ) ( ) ]
∫ ∑ ∫
∑
+
⋅ +
⋅ +
− ⋅
=
⋅
⋅
≅
⋅
⋅
≅
→
=
=
= b
a j
b a
j j
j
j j
x F x F x
F x
a F dx b
x L x F dx x F
x L x F x
F n
3 2
1 0
3
0
3
3 3 0
3 8 3
3
, (6)
Kwadratury te noszą kolejno nazwy wzorów prostokątów, trapezów, Simpsona i Newtona. Dokład- ność kwadratury można podnieść zwiększając dokładność interpolacji (patrz rozdział poświęcony interpolacji). W szczególności jest to możliwe przy rezygnacji z warunku równomiernego roz- mieszczenia węzłów interpolacji. I tak, jeżeli węzły interpolacji zostaną przyjęte w miejscach zero- wych wielomianów Legendre’a, będziemy mieli do czynienia z kwadraturami typu Gaussa – Le- gendre’a. Ponieważ wielomiany Legendre’a podobnie jak wielomiany Czebyszewa (patrz rozdział Interpolacja) mają wszystkie miejsca zerowe w przedziale
[
−1 ,1]
, to standardowo wzory te są po- dawane dla takiego właśnie przedziału, i wówczas wyglądają następująco:( ) ∑ ( )
∫
=−
⋅
=
⋅ n
j
j
j F ξ
w ξ
d ξ F
0 1
1
, (7)
gdzie w oznaczają wagi, a j ξ węzły kwadratury. Jeżeli całkowanie ma być dokonane w dowolnym j przedziale
[ ]
a,b konieczna jest zmiana zmiennych:[ ] [ ]
[ ] [ ]
a b
b a x
a b b
x a
b a b
a
−
−
−
= ⋅
⋅ − + +
=
−
−
, 1
, 1
1 , 1 ,
2 2
2 ξ
ξ
, (8)
co prowadzi do:
( ) ( ( ) ) ∑ ( ( ) )
∫ ∫
=−
⋅
− ⋅
=
⋅
− ⋅
=
⋅ n
j
j j
b
a
ξ x F a w
ξ b d ξ x a F dx b
x F
0 1
1 2
2 . (9)
Po przyjęciu położeń węzłów wartości wag wyznacza się w taki sposób, aby kwadratura dawała ścisłe wyniki dla jednomianów możliwie najwyższego stopnia. Ponieważ operacja taka jest dość żmudna i musi być wykonana z dużą dokładnością, praktyczniej jest korzystać z wartości tablico-
wych. Tablice takie można znaleźć w każdym podręczniku dotyczącym metod numerycznych. Tu- taj przytoczymy współczynniki wybranych kwadratur za [1].
Tablica 1. Kwadratury Gaussa rzędu 1 do 9
n i ξi(n) =−ξn(−ni)+1 wi(n) =wn(n−i)+1 1 1 0,0000000000 2,0000000000 2 1 0,5773502692 1,0000000000 3 1
2
0,7745966692 0,0000000000
0,5555555556 0,8888888889 4 1
2
0,8611363116 0,3399810436
0,3478548451 0,6521451549 5
1 2 3
0,9061798459 0,5384693101 0,0000000000
0,2369268851 0,4786286705 0,5688888889
Wyznaczenie wartości wag i położeń węzłów jest również możliwe wprost z warunku ścisłego cał- kowania jednomianów możliwie najwyższego stopnia. Warunek ten pozwala zapisać tyle równań nieliniowych (10), ile jest niewiadomych (wag i położeń węzłów). Niestety ten sposób postępowa- nia, pomimo, że najprostszy koncepcyjnie jest też najbardziej żmudny numerycznie, gdyż prowadzi do układu nieliniowych równań algebraicznych, który musi być rozwiązany z dużą dokładnością, ponieważ dokładność ta decyduje o jakości uzyskanych wyników całkowania.
[ ( ) ] ∑
∫
=+
− +
−
⋅
=
−
− + ⋅
= + ⋅
=
⋅ n
j
k j j k k
k w x
x k dx k
x
0 1 1
1 1 1
1
1 1 1
1 1
1 . (10)
I tak na przykład rozpisując stosowne zależności dla trójpunktowej kwadratury Gaussa, która pro- wadzi do wzorów numerycznego całkowania ścisłych dla wielomianów do stopnia piątego włącznie otrzymujemy następujący układ sześciu nieliniowych równań algebraicznych:
∫
∫
∫
∫
∫
∫
− −
− −
− −
− −
− −
− −
=
⋅ +
⋅ +
⋅
=
⋅
=
⋅
=
⋅ +
⋅ +
⋅
=
⋅
=
⋅
=
⋅ +
⋅ +
⋅
=
⋅
=
⋅
=
⋅ +
⋅ +
⋅
=
⋅
=
⋅
=
⋅ +
⋅ +
⋅
=
⋅
=
⋅
=
⋅ +
⋅ +
⋅
=
=
⋅
1
1
5 2 2 5 1 1 5 0 0 1
1 6 16 5
1
1
25 4 2 2 4 1 1 4 0 0 1
1 5 15 4
1
1
3 2 2 3 1 1 3 0 0 1
1 4 14 3
1
1
23 2 2 2 2 1 1 2 0 0 1
1 3 3 2 1
1
1
2 2 1 1 0 0 1
1 2 2 1 1
1
2 1 0 1
1
0 0 0 2 1 1 1 1
x w x w x w x
dx x
x w x w x w x
dx x
x w x w x w x
dx x
x w x w x w x
dx x
x w x w x w x
dx x
w w w x dx
. (11)
Ze wzorów (10) bądź (11) można natomiast natychmiast wywnioskować, że w przedziale
[
−1 ,1]
węzły i wagi są rozmieszczone symetrycznie względem punktu 0 , co wprost zostało wykorzystane przy konstruowaniu Tablicy 1.
1 J.Szmelter – Metody komputerowe w mechanice, PWN 1980
Możliwości podnoszenia dokładności kwadratury wiążą się oczywiście z:
• zwiększaniem liczby węzłów w przedziale całkowania (podnoszenie stopnia wielomianu in- terpolującego); istotnym ograniczeniem jest tu występowanie efektu Rungego (patrz Inter- polacja), powodujące, że do praktycznego stosowania nadają się wyłącznie kwadratury względnie niskiego (nie wyżej niż czwarty do siódmego) stopnia;
• wykorzystaniem addytywności całki, podzieleniem przedziału całkowania na wiele części i stosowanie w każdej z nich kwadratury niskiego rzędu (tak zwane kwadratury złożone).
Dla praktycznego prześledzenia sposobu postępowania przy numerycznym obliczaniu całek wyzna- czymy wartość wyrażenia:
( ) (
ln( ))
2.545177ln 14
4
1
=
−
⋅
=
∫
x ⋅dx x x x , (11)korzystając z różnych kwadratur. I tak kolejno korzystając ze wzorów (3)-(6) otrzymujemy odpo- wiednio dla kwadratur typu Newtona:
( )
( )
( )
⎪⎪
⎩
⎪⎪
⎨
⎧
=
=
=
=
⋅ +
⋅ +
⋅ +
⋅ +
⋅ +
⋅ +
⋅
=
⋅
−
−
∫
−535590 ,
2
525729 ,
2
079442 ,
2
000000 ,
0
) 0 , 4 ln(
) 0 , 3 ln(
3 ) 0 , 2 ln(
3 ) 0 , 1 ln(
) 0 , 4 ln(
) 5 , 2 ln(
4 ) 0 , 1 ln(
) 0 , 4 ln(
) 0 , 1 ln(
) 0 , 1 ln(
) ln(
81 4 61 4 21 4 13 4
1
dx
x , (12)
a dla kwadratur typu Gaussa po dokonaniu transformacji (8)1 na położeniach węzłów danych w Ta- blicy 1 w celu sprowadzenia ich do przedziału
[ ]
1,4 :{
⎪⎪
⎩
⎪⎪
⎨
⎧
=
⋅ +
=
=
⋅ +
=
=
⋅
−
=
=
⋅
−
=
⎪⎩
⎪⎨
⎧
=
⋅ +
=
=
⋅ +
=
=
⋅
−
=
⎩⎨
⎧
=
⋅ +
=
=
⋅
−
=
=
⋅ +
=
− +
− +
− +
− +
− +
− +
− +
− +
− +
− +
791704 ,
3 861136
, 0
009972 ,
3 339981
, 0
990028 ,
1 339981 ,
0
208296 ,
1 861136
, 0 :
661895 ,
3 774596
, 0
500000 ,
2 000000
, 0
338105 ,
1 774596
, 0 :
3,366025 577350
, 0
1,633974 577350
, : 0
500000 ,
2 000000
, 0 :
21 24 4
) 1 4 ( 4
21 24 4
) 1 4 ( 3
21 24 4
) 1 4 ( 2
21 24 4
) 1 4 ( 1 ) 4 (
21 24 4
) 1 3 ( 3
21 24 4
) 1 3 ( 2
21 24 4
) 1 3 ( 1 ) 3 (
21 24 4
) 1 2 ( 2
21 24 4
) 1 2 ( ) 1 2 (
21 24 4
) 1 1 ( 1 ) 1 (
x x x x x
x x x x
x x x
x x
, (13)
otrzymujemy ostatecznie:
( )
( )
( )
⎪⎪
⎩
⎪⎪
⎨
⎧
=
=
=
=
⋅
⋅ +
⋅ +
⋅ +
⋅
⋅
⋅ +
⋅ +
⋅
⋅
⋅ +
⋅
⋅
⋅
=
⋅
−
−
−
−
∫
545254 .
2
546084 ,
2
557122 ,
2
748872 ,
2
) ln(
) ln(
) ln(
) ln(
) ln(
) ln(
) ln(
) ln(
) ln(
) ln(
) ln(
21 ) 4 4 ( 4 ) 4 ( 4 )
4 ( 3 ) 4 ( 3 )
4 ( 2 ) 4 ( 2 )
4 ( 1 ) 4 ( 1
21 ) 4 3 ( 3 ) 3 ( 3 )
3 ( 2 ) 3 ( 2 )
3 ( 1 ) 3 ( 1
21 ) 4 2 ( 2 ) 2 ( 2 )
2 ( 1 ) 2 ( 1
21 ) 4 1 ( 1 ) 1 ( 1 4
1
w x w
x w
x w
x
w x w
x w
x
w x w
x
w x dx
x .(14)
Zestawmy tak obliczone wartości całki w tablicy w celu porównania dokładności uzyskanych wy- ników. W kolejnych kolumnach Tablicy 2 przedstawiono liczbę węzłów kwadratury n , wyniki cał- kowania kwadraturą typu Newtona i Gaussa oraz błąd względny procentowy tych wyników obli- czony w stosunku do wartości analitycznej całki.
Tablica 2. Wyniki całkowania numerycznego
Kwadratura Błąd
n Newtona Gaussa Newton Gauss
1 0,000000 2,748872 100,000% 8,003%
2 2,079442 2,557122 18,298% 0,469%
3 2,525729 2,546084 0,764% 0,036%
4 2,535590 2,545254 0,377% 0,003%
Jak widać w każdym przypadku przy porównywalnym nakładzie pracy (liczba obliczeń wartości funkcji podcałkowej) zastosowanie kwadratury Gaussa prowadzi do znacząco (od 12 do ponad 120 razy) dokładniejszych wyników.
Rozważmy jeszcze kwestię kwadratur złożonych. Jak już zostało to zaznaczone powyżej, w ich konstrukcji wykorzystujemy addytywność całki, która oznacza, że całkę po pewnym przedziale możemy zastąpić sumą całek po przedziałach których suma daje przedział wyjściowy, czyli:
∫ ∫ ∑
−∫ ∫
=
⋅ +
⋅ +
⋅
=
⋅ +
b a
p a
m i
b p p
p m
i
i
dx x f dx x f dx
x f dx x
f 1 1 1
1
) ( )
( )
( )
( , (15)
gdzie oczywiście a< p1<...< pi <...< pm<b .
W przypadku kwadratur typu Newtona przy stałej długości przedziału całkowania
[
pi,pi+1]
równej h dla (16) i (17), 2⋅h dla (18) i 3⋅h dla (19) (h oznacza tu odległość pomiędzy węzłami kwadra- tury) postępowanie powyższe prowadzi do wzorów:• interpolacja wielomianowa stopnia 0 (n=0), czyli funkcją stałą:
( ) ( )
∫ ∑
=
⋅ +
⋅
≅
b ⋅
a
m i
i h a h F
dx x F
1 0 , (16)
• interpolacja wielomianowa stopnia 1 (n=1), czyli funkcją liniową:
( ) ( ) ( )
∫
⋅ ≅ ⋅ ⋅⎢⎣⎡ + ⋅∑
− + ⋅ + ⎥⎦⎤= b
a
m i
b F i h a F a
h F dx x
F 1
1
) (
2 2
1 , (17)
• interpolacja wielomianowa stopnia 2 (n=2), czyli funkcją kwadratową:
( ) ( ) ( ) ( )
∫ ∑ ∑
⎥⎦
⎢ ⎤
⎣
⎡ + ⋅ + ⋅ + ⋅ + ⋅ +
⋅ ⋅
≅
⋅ −
=
−
= b
a
m i m
i
b F i h a F i
h a F a
h F dx x
F 2
,...
4 , 2 1
,...
3 , 1
) (
2 6 4
2 , (18)
• interpolacja wielomianowa stopnia 3 (n=3), czyli funkcją trzeciego stopnia:
( ) ( ) ( ) ( ) ( ) ( )
∫ ∑ ∑ ∑
⎥⎦
⎢ ⎤
⎣
⎡ + ⋅ + ⋅ + ⋅ + ⋅ + ⋅ + ⋅ +
⋅ ⋅
≅
⋅ −
=
−
=
−
= b
a
m i m
i m
i
b F i h a F i
h a F i
h a F a
h F dx x
F 1
,...
5 , 2 3
,...
6 , 3 2
,...
4 , 1
3 2
8 3
3 , (19)
Postępowanie to w przypadku kwadratur typu Gaussa jest równoznaczne z zastosowaniem prostej kwadratury odpowiedniego rzędu do każdego podprzedziału z osobna i zsumowaniem tak uzyska- nych wyników.
Dla zorientowania się, jaki jest wpływ gęstości podziału przedziału całkowania na dokładność wy- niku finalnego rozwiążemy zadanie (11) tak dobierając liczbę podprzedziałów całkowania o jedna- kowej długości w każdej z metod, aby wynik finalny uzyskać z błędem względnym procentowym nie przekraczającym 0,003%. Wyniki tych obliczeń przedstawiono w Tablicy 3.
Tablica 3. Porównanie efektywności metod całkowania
Newton Gauss
n h ile całka błąd h ile całka błąd
1 0,00011 26850 2,545100 0,003% 0,05000 60 2,545255 0,003%
2 0,03488 87 2,545101 0,003% 0,60000 10 2,545230 0,002%
3 0,30000 11 2,545098 0,003% 1,50000 6 2,545224 0,002%
4 0,20000 13 2,545142 0,003% 3,00000 4 2,545254 0,003%
W kolejnych kolumnach tablicy, dla obydwu metod całkowania, przedstawiono długość poje- dynczego przedziału całkowania (h), całkowitą liczbę punktów, w których należało obliczyć war- tość funkcji podcałkowej (ile), obliczoną wartość całki (całka) oraz błąd względny procentowy uzy- skanego wyniku w stosunku do wartości ścisłej (błąd).
Ponieważ zasadniczym czynnikiem wpływającym na czas całkowania jest liczba obliczeń wartości funkcji podcałkowej, zgodnie z oczekiwaniem okazało się, że najbardziej efektywna z rozważanych powyżej jest czteropunktowa kwadratura Gaussa.
Problemem z którym spotykamy się dość często przy numerycznym całkowaniu jest problem granic niewłaściwych (czyli w − lub ∞ +∞ ). Możemy sobie z nim poradzić na dwa sposoby:
• przez odpowiednią zmianę zmiennych (typu x= ): 1t t dt t f
dx x f
a
b
b a
⎟⋅
⎠
⎜ ⎞
⎝
⋅ ⎛
=
⋅
∫
∫
1
1
1 ) 1
( 2 ; (20)
• zastosowanie kwadratur specjalnych:
∫ ∑
∫ ∑
∫ ∑
=
∞ =
∞
∞
−
−
=
∞ −
⋅
=
⋅
⋅ +
⋅
−
⋅
=
⋅
⋅
⋅
=
⋅
⋅
⋅
n
j j j
n
j j j
x
n
j j j
x
x f w dx x f x x
x f w dx x f e
x f w dx x f e x
0 1
1 0 1
) ( )
( ) 1 ( ) 1 (
) ( )
(
) ( )
(
2
β α
α
, (21)
kwadratury te noszą odpowiednio nazwy Gaussa–Laguerre’a, Gaussa–Hermita i Gaussa–
Jacobiego.
Stosując kwadratury te należy pamiętać o tym, że ich współczynniki (położenia węzłów i wagi) są określone i mogą być używane jedynie wówczas gdy przedział całkowania i funkcja podcałkowa mają postać dokładnie taką jak we wzorach (21).
W przypadku osobliwości całkowalnej (np. typu (22) ) powinniśmy najpierw przez odpowiednią zmianę zmiennych pozbyć się osobliwości i dopiero wtedy zastosować którąś z standardowych kwadratur.
1 0
, ) (
1 ≤ <
− γ γ
a
x . (22)
Stosowne wzory przedstawiono poniżej:
∫ ∫
−
−
−
−
− ⋅ + ⋅
− ⋅
=
⋅
=
−
b a
a b
dt a t f t dx
x f
t a
x γ
γ γ γ
γ
1
11 11 11
) (
0
) 1 (
) 1
( . (23)
Rozpatrzmy takie postępowanie na bardzo prostym przykładzie:
∫
= ⋅∫
⋅⎭⎬
⎫
⎩⎨
⎧
⋅
⋅
=
= =
4 ⋅
0
2
0 2 2
) cos(
2 2 )
cos( t dt
dt t dx
t dx x
x
x . (24)
W przypadku całkowania w obszarze dwuwymiarowym korzystamy z twierdzenia, że całkowanie w takim obszarze po powierzchni czworokąta D może być zastąpione całkowaniem po powierzchni obszaru wzorcowego (kwadrat ∆ ). Warunkiem poprawności tego postępowania jest skonstruowa- nie odwzorowania, w którym obszar D jest przedstawiony jako obraz obszaru ∆ . Stosowne od- wzorowanie można zbudować posługując się wielomianami interpolacyjnymi Lagrange’a w dwóch wymiarach (patrz rozdział dotyczący interpolacji, wzory (68) i następne, w których elementem wzorcowym jest obszar ∆ , a wzory transformacyjne x= u(ξ,η) i y= v(ξ,η) mają postać (69)).
Wówczas, jeżeli tylko spełniona jest zależność:
) 0 , ( ) , (
) , ( ) , (
≠ η
∂ η ξ
∂ ξ
∂ η ξ
∂ ∂η
η ξ
∂ ξ
∂ η ξ
∂
= v v
u u
J , (25)
dla każdego punktu obszaru ∆ , to prawdziwy jest wzór:
( ) ( )
( )
∫∫
∫∫
∆η ξ
⋅
⋅ η ξ η ξ
=
⋅dxdy f u v J d d
y x f
D
, , , )
,
( , (26)
i wtedy całkę w dwóch wymiarach można obliczyć jako złożenie dwóch całek w jednym wymiarze, co ostatecznie prowadzi do zależności:
( ) ( ( ) ( ) )
∫∫ ∑ ∑
= =
η ξ η ξ
⋅
⋅
=
⋅
D
N j
N
j j j j j j j
x
x y
y
y x y x y
x w F x y
w dxdy
y x f
1 1
, , ,
, , (27)
gdzie:
( ) ( ) ( )
η
∂
∂ ξ
∂
∂ ∂η
∂ ξ
∂
∂
⎟⋅
⎠
⎜ ⎞
⎝
⎛ ξ η ⋅ ξ η ⋅
= η
ξ
∑ ∑
=
= y y
x x y N
x N
f F
i i i
i i i
4
1 4
1
, ,
,
, , (28)
(patrz wzory (69) w rozdziale dotyczącym interpolacji).
Warto zauważyć, że ze wzorów (27) wynika możliwość całkowania w kierunkach ξ i η przy uży- ciu kwadratur różnego rzędu (z różną dokładnością). Postępowanie takie może być uzasadnione w przypadku rozwiązywania zadań, w których zmienne podstawowe mają różny charakter (np. prze- strzeń i czas).
Podaną procedurę można łatwo uogólnić na przypadki całkowania w obszarze trój i więcej wymia- rowym.
