Katolicki Uniwersytet Lubelski Jana Pawła II WYDZIAŁ MATEMATYKI, INFORMATYKI I ARCHITEKTURY KRAJOBRAZU KIERUNEK: INFORMATYKA Anna Sagan NR ALBUMU: 125...
17 downloads
25 Views
2MB Size
Katolicki Uniwersytet Lubelski Jana Pawła II WYDZIAŁ MATEMATYKI, INFORMATYKI I ARCHITEKTURY KRAJOBRAZU KIERUNEK: INFORMATYKA
Anna Sagan NR ALBUMU: 125372
ANALIZA PORÓWNAWCZA ALGORYTMÓW NUMERYCZNYCH STOSOWANYCH W OBLICZANIU CAŁEK
Praca magisterska napisana na seminarium Algorytmizacja i modelowanie pod kierunkiem prof. Piotra Matusa
Lublin 2015
Spis treści Wstęp
2
1 Podstawy teoretyczne
3
2 Metody całkowania numerycznego
5
2.1
Metody całkowania wielomianu interpolacyjnego - zagadnienia wstępne . .
5
2.2
Metody Newtona-Cotesa . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
2.3
Metody Gaussa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.4
Ekstrapolacja Richardsona oraz jej zastosowanie w metodzie Romberga . . 22
3 Porównanie wybranych metod numerycznych
24
3.1
Metody całkowania numerycznego Newtona-Cotesa . . . . . . . . . . . . . 24
3.2
Porównanie dokładności metod Newtona-Cotesa i Gaussa . . . . . . . . . . 31
3.3
Metoda Romberga . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.4
Implementacja metody . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
Zakończenie
44
Bibliografia
45
1
Wstęp Całka zajmuje jedno z czołowych miejsc w gronie najważniejszych pojęć analizy matematycznej. Ma bardzo szerokie zastosowanie - wykorzystuje się ją w matematyce, fizyce, technice, a także wielu innych dziedzinach nauki. Obliczanie całki oznaczonej jednak nie zawsze jest łatwe, czasem jest bardzo trudne, a nawet niewykonalne. Szczególnie w sytuacjach, gdy nie możemy wyznaczyć funkcji pierwotnej całkowanej funkcji lub obliczenia są zbyt skomplikowane, szukamy innych sposobów, by uzyskać zadowalające rezultaty. Często wykorzystuje się w tym celu całkowanie numeryczne, czyli przybliżone wyznaczanie wartości całki oznaczonej. Celem niniejszej pracy jest porównanie metod numerycznych, ocenienie ich efektywności. W pierwszym rozdziale zdefiniujemy podstawowe pojęcia związane z całką oznaczoną oraz wymienimy niektóre jej własności. W drugim rozdziale będziemy zajmować się zagadnieniem aproksymacji. Podamy definicje bazy aproksymacyjnej, funkcji aproksymującej oraz błędu aproksymacji, a następnie przedstawimy niektóre rodzaje aproksymacji. Omówimy także zagadnienie interpolacyjne. W tym celu zdefiniujemy wielomian interpolacyjny Lagrange’a oraz wykażemy jego istnienie i jednoznaczność. Ponadto oszacujemy błąd interpolacji Lagrange’a. W trzecim rozdziale przedstawimy ogólne schematy całkowania numerycznego. W szczególności omówimy tzw. kwadratury proste oraz złożone, wśród których wyróżnimy klasyczne typy całkowania numerycznego znane powszechnie jako metody Newtona-Cotesa oraz metody Gaussa. W kolejnym rozdziale przedstawimy wybrane warianty kwadratur różnych typów, podamy też konkretne wzory służące do ich wyznaczania oraz oszacowanie błędu. Ostatni rozdział zawiera implementacje wybranych złożonych metod Newtona - Cotesa w języku programowania C#.
2
Rozdział 1 Podstawy teoretyczne Definicja 1.1. Niech n ∈ Z1 . Dla dowolnego ciągu rosnącego Z0,n 3 k 7→ xn (k) ∈ [a, b] takiego, że xn (0) = a oraz xn (n) = b definiujemy λ(xn ) := max [xn (k + 1) − xn (k)] . k∈Z0,n−1
Liczba λ(xn ) jest średnicą podziału xn , a xn nazywamy podziałem przedziału [a, b]. Definicja 1.2. Niech n ∈ Z1 oraz niech f : [a, b] 7→ R.
Dla dowolnego ciągu
Z0,n 3 k 7→ xn (k) ∈ [a, b] takiego, że xn (0) = a, xn (n) = b oraz dowolnego ciągu Z0,n−1 3 k 7→ tn (k) ∈ [a, b] takiego, że tn (k) ∈ [xn (k), xn (k + 1)] dla k ∈ Z0,n−1 definiujemy σ (xn , tn , f ) :=
n−1 X
f (tn (k)) · [xn (k + 1) − xn (k)] .
k=0
Jest to suma całkowa. Definicja
1.3. Niech
f
:
7→
[a, b]
R
będzie
ograniczona
w
[a, b].
Mówimy, że f jest całkowalna w sensie Riemanna w przedziale [a, b] wtedy i tylko wtedy, gdy dla dowolnego ciągu Z1 3 n 7→ xn takiego, że xn jest rosnący i xn (0) = a, xn (n) = b i lim λ (xn ) = 0 oraz dowolnego ciągu Z1 3 n 7→ tn takiego, że tn (k) ∈ [xn (k), xn (k + 1)] n→∞
dla k ∈ Z0,n−1 istnieje granica I := lim σ (xn , tn , f ) . n→∞
Liczbę I nazywamy całką oznaczoną funkcji f w przedziale [a, b].
3
Dla funkcji całkowalnych f, g : [a, b] 7→ R oraz liczb c ∈ R, d ∈ (a, b) zachodzą poniższe własności: 1. Addytywność Zb
f (x) + g(x)dx =
Zb
f (x)dx +
2. Jednorodność
Zb
cf (x)dx = c
a
g(x)dx
a
a
a
Zb
Zb
f (x)dx
a
3. Addytywność względem przedziału całkowania Zb
f (x)dx =
a
Zd a
f (x)dx +
Zb
f (x)dx
d
Zapiszemy teraz Podstawowy wzór rachunku całkowego, który opisuje związek między całką oznaczoną i nieoznaczoną. Jest on znany jako Twierdzenie Newtona-Leibniza (zob. [1, str.106]). Twierdzenie 1.4 (Newton - Leibniz). Jeśli f jest funkcją ciągłą na przedziale [a, b], to wówczas
Zb
f (x)dx = F (b) − F (a),
a
gdzie F jest dowolną funkcją pierwotną z funkcji f , a więc F 0 = f .
4
Rozdział 2 Metody całkowania numerycznego 2.1 2.1.1
Metody całkowania wielomianu interpolacyjnego - zagadnienia wstępne Aproksymacja
Niech X będzie dowolną unormowaną przestrzenią liniową nad ciałem R, której nośnikiem jest zbiór {f : [a, b] 7→ R} dla pewnych a, b ∈ R, a < b. Definicja 2.1 (baza aproksymacyjna). Niech f ∈ X i niech Xm będzie dowolną m-wymiarową podprzestrzenią liniową przestrzeni X dla m ∈ N.
Dowolną bazę
{ϕ1 , . . . , ϕm } przestrzeni Xm nazywamy bazą aproksymacyjną funkcji f w Xm . Definicja 2.2 (funkcja aproksymująca). Dowolną funkcję F ∈ X nazywamy funkcją aproksymującą daną funkcję f
∈
X względem danej bazy aproksymacyjnej
{ϕk ∈ Xm , k ∈ Z1,m }, gdzie Xm jest pewną podprzestrzenią przestrzeni X, wtedy i tylko wtedy, gdy istnieje funkcja g : Rm 7→ R taka, że F = g(ϕ1 , ϕ2 , . . . , ϕm ). Definicja 2.3 (błąd aproksymacji). Liczbę b := kF − f k nazywamy błędem aproksymacji danej funkcji f ∈ X funkcją aproksymującą F względem danej bazy aproksymacyjnej {ϕk ∈ Xm }. Aproksymacja polega na tym, by dla danej funkcji f ∈ X dobrać bazę {ϕk } oraz funkcję aproksymującą w taki sposób, aby błąd aproksymacji b nie przekraczał danej wartości b0 > 0. Uwaga 2.4. Dla danej funkcji f ∈ X oraz ustalonej bazy aproksymacyjnej zagadnienie aproksymujące może mieć wiele rozwiązań lub może nie mieć żadnego rozwiązania.
5
2.1.1.1
Aproksymacja liniowa
W praktyce wybór funkcji aproksymującej F funkcję f ∈ X względem bazy aproksymacyjnej {ϕk } często zawęża się do pewnej podprzestrzeni Y przestrzeni X. W większości przypadków jest to Y = Xm . Wówczas funkcja aproksymująca jest postaci F =
m X
ak ϕk
k=1
dla pewnych liczb ak ∈ R. Taki przypadek nazywamy aproksymacją liniową (patrz [2, str.74]). 2.1.1.2
Aproksymacja wymierna
Możemy spotkać się z sytuacją, gdy funkcję aproksymującą F poszukujemy w ogólniejszej postaci, czyli
p P
F =
k=1 q P
ak θ k , bk ψk
k=1
gdzie θk , ψk ∈ {ϕk } ⊂ Xm ; ak , bk ∈ R, przy czym m = max {p, q}. W takim przypadku mówimy o tzw. aproksymacji wymiernej (patrz [2, str.74]).
2.1.2
Interpolacja
Rozważmy aproksymację z dodatkowym warunkiem nałożonym na funkcję F aproksymującą f ∈ X w bazie aproksymacyjnej {ϕk } postaci: F (xi ) = f (xi ) dla danych wartości xi ∈ [a, b]. Tak postawiony problem nazywamy zagadnieniem interpolacyjnym. 2.1.2.1
Interpolacja Lagrange’a
Definicja 2.5 (wielomian interpolacyjny Lagrange’a). Niech n ∈ Z0 . Dla dowolnej funkcji [a, b] 3 x 7→ f (x) ∈ R oraz dowolnych punktów x0 , x1 , . . . , xn ∈ [a, b] , takich, że xi 6= xj dla i 6= j,
i, j ∈ Z0,n wielomian Ln (x) :=
n X
f (xi )Φi
(2.1.1)
i=0
nazywamy n-tym
wielomianem
interpolacyjnym
Lagrange’a
funkcji
f
(zob.[3, 276]), gdzie Φi =
n Y
x − xj . j=0 xi − xj j6=i
6
(2.1.2)
Twierdzenie 2.6. Niech n ∈ Z0 . Dla dowolnej funkcji f : [a, b] 7→ R oraz dowolnych punktów x0 , x1 , . . . , xn ∈ [a, b] ,
xi 6= xj dla i 6= j,
i, j ∈ Z0,n istnieje dokładnie jeden
wielomian P : R 7→ R stopnia co najwyżej n taki, że P (xi ) = f (xi ), gdzie i ∈ Z0,n . Twierdzenie 2.7 (oszacowanie błędu interpolacji Lagrange’a). Jeśli f : [a, b] 7→ R jest funkcją należącą do klasy C n+1 ([a, b]), zaś P jest wielomianem stopnia n takim, że f (xi ) = P (xi ), i = 0, 1, . . . , n dla pewnych liczb xi ∈ [a, b], przy czym xi < xj , gdy i < j, to dla dowolnego x ∈ [a, b] istnieje ξx ∈ (a, b) takie, że f (x) − P (x) =
2.1.3
n Y 1 f (n+1) (ξx ) (x − xi ). (n + 1)! i=0
(2.1.3)
Kwadratury
Niech X oznacza dowolną unormowaną przestrzeń liniową nad ciałem R, której nośnikiem jest zbiór {f : [a, b] 7→ R} dla pewnych a, b ∈ R;
a < b. Niech f będzie całkowalna
(w sensie Riemanna) w [a, b]. Oznaczmy I(f ) =
Zb
f (x)dx.
a
Zagadnienie przybliżonego całkowania funkcji f ∈ X możemy sformułować następująco: należy dobrać funkcję F aproksymującą funkcję f w ten sposób, aby |I(f ) − S(f )| < b0 , gdzie S(f ) =
Zb
F (x)dx
a
jest szukanym przybliżeniem całki I(f ) dla ustalonej liczby b0 > 0. Definicja 2.8. Mówimy, że kwadratura numeryczna jest rzędu r, jeśli dla dowolnego wielomianu wn stopnia mniejszego od r (n = 0, 1, · · · , r − 1) daje dokładny wynik (tzn. I(wn ) = S(wn )) oraz istnieje wielomian wr stopnia r taki, że I(wr ) 6= S(wr ). 2.1.3.1
Kwadratury proste
Kwadraturą prostą nazywamy przybliżenie całki oznaczonej na przedziale [a, b] (por. kwadratury złożone). Do przybliżonego obliczania kwadratur prostych z całki I(f ) =
Zb
p(x)f (x)dx,
a
7
gdzie p(x) jest pewną funkcją wagową dodatnią w przedziale [a, b], będziemy używać wzorów postaci I(f ) = S(f ) + R(f ),
(2.1.4)
w których R(f ) jest błędem kwadratury, zaś szukane przybliżenie całki I(f ) jest postaci: S(f ) = I(LN ), gdzie LN =
N X
Φk (x)f (xk )
k=0
jest wielomianem interpolacyjnym Lagrange’a (patrz wzór (2.1.1)), a Φi =
n Y
x − xj j=0 xi − xj j6=i
(patrz wzór (2.1.2)). Oczywiście LN (x0 ) = f (x0 ), LN (x1 ) = f (x1 ), . . . , LN (xN ) = f (xN ). Otrzymujemy więc S(f ) = I
N X
!
Φk (x)f (xk ) =
k=0
Zb a
=
N X
f (xk )
k=0
Zb
p(x)
N X
Φk (x)f (xk )dx =
k=0
p(x)Φk (x)dx.
a
Oznaczając Ak :=
Zb
p(x)Φk (x)dx
(2.1.5)
a
dostajemy S(f ) =
N X
Ak f (xk ),
(2.1.6)
k=0
przy czym Ak nazywamy współczynnikami kwadratury. Twierdzenie 2.9. Kwadratura postaci (2.1.6) jest rzędu co najmniej N + 1. Dowód. Udowodnimy, że kwadratura (2.1.6) jest dokładna dla wielomianu co najwyżej N stopnia. Niech f będzie dowolnym wielomianem stopnia co najwyżej N . Wtedy LN (x) =
N X
f (xk )Φk = f,
k=0
8
a stąd I(f ) =
Zb
p(x)f (x)dx =
a
=
N X k=0
Zb
N X
p(x)
f (xk )Φk dx =
k=0
a
b Z p(x)Φk · f (xk )dx
!
N X
=
Ak f (xk ) = S(f ).
k=0
a
2.1.3.2
Kwadratury złożone
W praktyce stosowanie kwadratur wyższych rzędów jest uciążliwe. Częściej używanym rozwiązaniem jest zastosowanie kwadratur złożonych, czyli podzielenie przedziału całkob−a wania [a, b] na m równych części o długości H = i stosowanie na każdym z nich m kwadratury niskiego rzędu. Wzór ogólny kwadratury złożonej możemy zapisać w postaci: S(f ) =
m−1 X
Si (f|[wi ,wi+1 ] ),
i=0
gdzie wi = a + iH, i = 0, 1, · · · , m − 1, zaś Si (f|[wi ,wi+1 ] ) jest pewną kwadraturą prostą. Podobnie obliczamy błąd dla kwadratur złożonych, sumując błędy z poszczególnych przedziałów, czyli R(f ) =
m−1 X
Ri (f|[wi ,wi+1 ] ),
i=0
gdzie wi = a + iH, i = 0, 1, · · · , m − 1.
2.2 2.2.1
Metody Newtona-Cotesa Kwadratury proste
Jeżeli x0 , x1 , x2 , . . . , xN −1 , xN są równo odległymi węzłami interpolacji funkcji f (x) w przedziale [a, b], czyli a ≤ x0 < x1 < x2 < . . . < xN −1 < xN ≤ b, gdzie xi są elementami dziedziny f , dla których znana jest wartość f (xi ), to przybliżenie całki I(f ) będzie postaci (2.1.6). Węzły w metodach Newtona-Cotesa są ustalone, więc naszym zadaniem pozostaje obliczenie współczynników kwadratury Ak danych wzorem (2.1.5), gdzie waga p(x) = 1.
9
2.2.1.1
Wzory zamknięte
Wzorami zamkniętymi Newtona-Cotesa nazywamy kwadratury, w których końce przedziału są węzłami interpolacji funkcji f (x), czyli a = x0 < x1 < x2 < . . . < xN −1 < xN = b. W takim przypadku h =
b−a oznacza odległość między węzłami interpolacji i x ∈ [a, b] N
można wyrazić w postaci x = a + th, gdzie t ∈ [0, N ]. Wtedy Φk przyjmuje postać Φk (x) =
N Y
t−j . j=0 k − j j6=k
Stąd ostatecznie całkując przez podstawienie (x = a + th) mamy S(f ) =
Zb
LN (x)dx =
N X
f (xk ) · h ·
k=0
a
ZN Y N 0
t−j dt, j=0 k − j
(2.2.1)
j6=k
czyli współczynniki Ak przyjmują postać Ak = h ·
ZN Y N 0
t−j dt. j=0 k − j
(2.2.2)
j6=k
Błąd kwadratur zamkniętych Newtona-Cotesa wyraża się wzorem (por.[2, str.165]):
N hN +3 · f (N +2) (ξ) RN Q · t( (t − i))dt dla N = 2k, k = 0, 1, 2, · · · (N + 2)! i=0 0
N hN +2 · f (N +1) (ξ) RN Q · ( (t − i))dt (N + 1)! 0 i=0
R(f ) =
dla N = 2k + 1, k = 0, 1, 2, · · ·
Przedstawimy teraz kilka wzorów zamkniętych Newtona-Cotesa niskiego rzędu. Wzór trapezów to kwadratura, w której mamy dane dwa węzły (N = 1) będące końcami przedziału [a, b], czyli x0 = a,
x1 = b . Odległość między węzłami wynosi
h = b − a.
10
Rysunek 2.1: Metoda trapezów Wzór jest równy: S(f ) = I(L1 ) =
h b−a (f (a) + f (b)) = (f (a) + f (b)) . 2 2
Błąd wzoru trapezów jest równy
R(f ) = −
(b − a)3 · f (2) (ξ) h3 = − f (2) (ξ). 12 12
Wzór Simpsona (wzór parabol) służy nam do obliczenia kwadratury przy danych trzech węzłach (x0 = a, h=
x1 =
a+b , 2
x2 = b). Odległość między węzłami jest równa
b−a . 2
11
Rysunek 2.2: Metoda Simpsona Mamy więc b−a S(f ) = I(L2 ) = f (a) + 4f 6
!
!
h a+b + f (b) = f (a) + 4f 2 3
!
!
a+b + f (b) 2
Błąd wzoru Simpsona wynosi (b − a)5 · f (4) (ξ) h5 (4) R(f ) = − = − f (ξ) 2880 90
Wzór 3/8 (Simpsona 3/8) oparty jest na znajomości czterech węzłów (N = 3). Odległość między węzłami jest równa h = x2 =
2b + a , 3
b−a , stąd x0 = a, 3
x1 =
x3 = b. Wzór jest równy
b−a S(f ) = I(L3 ) = f (a) + 3f 8 3 = h f (a) + 3f 8
!
2a + b + 3f 3 !
2a + b + 3f 3
12
!
!
2b + a + f (b) = 3 !
!
2b + a + f (b) . 3
2a + b , 3
Błąd wzoru 3/8 jest równy R(f ) =
(b − a)5 · f (4) (ξ) 3(h5 ) (4) =− f (ξ) 6480 80
Wzór Boole’a obliczamy, znając pięć węzłów (N = 4): x0 = a, x2 =
a+b , 2
x3 =
3b + a , 4
3a + b , 4
x1 =
x4 = b, przy odległości między węzłami h =
b−a . 4
Wzór Boole’a jest równy: b−a S(f ) = I(L4 ) = 7f (a) + 32f 90 2 = h 7f (a) + 32f 45
!
!
3a + b + 12f 4 !
3a + b + 12f 4
a+b + 32f 2 !
a+b + 32f 2
!
!
3b + a + 7f (b) = 4 !
!
3b + a + 7f (b) . 4
Błąd wzoru Boole’a wynosi: R(f ) = −
8(h7 ) (6) (b − a)7 · f (6) (ξ) =− f (ξ). 1935360 945
2.2.1.2
Wzory otwarte
Kwadratury Newtona-Cotesa, w których końce przedziału nie są węzłami interpolacji funkcji f (x), nazywamy wzorami otwartymi. Możemy to zapisać następująco: a = x−1 < x0 < x1 < x2 < . . . < xN −1 < xN < xN +1 = b, gdzie x−1 i xN +1 są elementami dziedziny f niebędącymi węzłami interpolacji funkcji f . W tym przypadku h =
b−a oznacza odległość między węzłami interpolacji. Wówczas N +2
x ∈ [a, b] można wyrazić w postaci x = a + (t + 1)h, gdzie t ∈ [−1, N + 1]. Postępując podobnie jak przy wzorach zamkniętych otrzymujemy
S(f ) =
Zb a
LN (x)dx =
N X
f (xk ) · h ·
k=0
czyli współczynniki Ak przyjmują postać
13
N Z+1 Y N −1
t−j dt, j=0 k − j
j6=k
(2.2.3)
Ak = h ·
N Z+1 Y N −1
t−j dt. j=0 k − j
(2.2.4)
j6=k
Błąd kwadratur otwartych Newtona-Cotesa obliczamy ze wzoru (por.[2, str.165]):
R(f ) =
N hN +3 · f (N +2) (ξ) NR+1 Q · t( (t − i))dt dla N = 2k, k = 0, 1, 2, · · · (N + 2)! i=0 −1
N hN +2 · f (N +1) (ξ) NR+1 Q · ( (t − i))dt (N + 1)! −1 i=0
dla N = 2k + 1, k = 0, 1, 2, · · ·
Przyjrzyjmy się kilku wzorom otwartym Newtona - Cotesa niskiego rzędu: Metoda prostokątów (Midpoint rule) to kwadratura, którą obliczamy, znając węzeł x0 =
b−a a+b i mając h = . 2 2
Rysunek 2.3: Metoda prostokątów
14
Wzór prostokątów jest więc równy a+b 2
S(f ) = I(L0 ) = (b − a)f
!
!
= 2hf
a+b , 2
a błąd tego wzoru wynosi
R(f ) =
(b − a)3 · f (2) (ξ) h3 = f (2) (ξ). 24 3
Metoda trapezów pozwala na obliczenia przy znajomości dwóch węzłów x0 = i x1 =
a + 2b b−a . Wtedy N = 1, h = . 3 3
Rysunek 2.4: Metoda trapezów Wzór jest równy b−a S(f ) = I(L1 ) = f 2
3h = f 2
!
2a + b +f 3
!
2a + b +f 3
15
a + 2b 3
a + 2b 3
!!
.
!!
=
2a + b 3
Błąd metody trapezów wynosi R(f ) =
(b − a)3 · f (2) (ξ) 3h3 (2) = f (ξ). 36 4
Uwaga 2.10. Wzór otwarty metody trapezów oraz wzór zamknięty są tej samej postaci S(f ) =
b−a (f0 + f1 ). 2
Metoda Milne’a jest kwadraturą opartą na trzech węzłach: x1 =
a+b , 2
x2 =
3b + a . 4
x0
=
Odległość między nimi jest równa h =
3a + b , 4 b−a . 4
Wzór wynosi !
b−a S(f ) = I(L2 ) = 2f 3 4h = 2f 3
!
3a + b −f 4 !
3a + b −f 4
a+b + 2f 2 !
a+b + 2f 2
3b + a 4
3b + a 4
!!
=
!!
.
Błąd metody Milne’a jest równy R(f ) =
14h5 (4) 7(b − a)5 · f (4) (ξ) = f (ξ). 23040 45
2.2.2
Kwadratury złożone
2.2.2.1
Wzory zamknięte
Złożony wzór trapezów możemy przedstawić w postaci S(f ) =
m−1 X i=0
H 1 1 (f (wi ) + f (wi+1 )) = h · f (w0 ) + f (w1 ) + · · · + f (wm−1 ) + f (wm ) , 2 2 2
gdzie odległość między węzłami h = H =
b−a . m
16
Rysunek 2.5: Złożona metoda trapezów Błąd złożonego wzoru trapezów jest równy R(f ) = −
(b − a)h2 (2) f (ξ). 12
Złożony wzór parabol stosujemy w postaci: S(f ) =
m−1 X i=0
h wi + wi+1 f (wi ) + 4f + f (wi + 1) = 3 2
m−1 m−1 X X h wi + wi+1 f (w0 ) + 4 f = +2 f (wi ) + f (wm ) , 3 2 i=0 i=1
"
gdzie odległość między węzłami jest równa h =
H b−a = . 2 2m
Błąd złożonego wzoru Simpsona wynosi R(f ) = −
(b − a)h4 (4) f (ξ). 180
17
#
2.2.2.2
Wzory otwarte
Złożony wzór prostokątów zapisujemy w następujący sposób S(f ) =
m−1 X i=0
wi + wi+1 2hf 2
=
w0 + w1 w1 + w2 wm−1 + wm =H· f +f + ··· + f 2 2 2 b−a H = . gdzie odległość między węzłami jest równa h = 2 2m
,
Rysunek 2.6: Złożona metoda prostokątów Błąd złożonego wzoru prostokątów jest równy R(f ) =
2.3
(b − a)h2 (2) f (ξ). 6
Metody Gaussa
Rozważane dotąd metody Newtona-Cotesa są kwadraturami z równo odległymi węzłami. Kwadratury Gaussa są natomiast całkami postaci (2.1.6), w których dobierane są nie tylko współczynniki Ak , ale też węzły (nie muszą być równo odległe) - w taki sposób, aby kwadratura była dokładna dla możliwie najwyższego stopnia wielomianu. Ponieważ 18
musimy znaleźć (N +1) współczynników oraz (N +1) węzłów, czyli mamy 2(N +1) stałych do wyznaczenia, należy przypuszczać, iż można je tak dobrać, by kwadratura była rzędu 2(N + 1). Twierdzenie 2.11. Nie istnieje kwadratura postaci (2.1.6) rzędu wyższego niż 2(N + 1). 2 Dowód. Przyjrzyjmy się wielomianowi W (x) = ωN +1 (x) stopnia 2(N + 1), gdzie
ωN +1 (x) := (x − x0 ) · · · (x − xN )
(2.3.1)
dla dowolnie wybranych x0 , x1 , · · · , xN takich, że xi 6= xj , gdy i 6= j. Mamy I(W ) =
Zb
p(x)W (x)dx > 0,
a
gdyż p(x)W (x) > 0 dla x ∈ [a, b]\ {xk : k = 0, 1, · · · , N } oraz S(W ) =
N X
Ak W (xk ) = 0,
k=0
gdyż W (xk ) = 0, k = 0, 1, · · · , N. Istnieją więc wielomiany stopnia 2(N + 1), dla których kwadratura (2.1.6) nie jest dokładna, co należało udowodnić. Twierdzenie 2.12. Kwadratura (2.1.6), o współczynnikach określonych wzorem Ak =
Zb
p(x)Φk (x)dx,
(2.3.2)
a
gdzie p(x) jest wagą, a Φk (x) =
N Y
x − xj , j=0 xk − xj j6=k
jest rzędu 2(N + 1) wtedy i tylko wtedy, gdy xk są pierwiastkami wielomianu PN +1 (x) z ciągu P0 (x), P1 (x), · · · , PN (x), · · · , wielomianów ortogonalnych na przedziale [a, b] względem wagi p, który da się jednoznacznie wyznaczyć przy pewnych dodatkowych założeniach normalizacyjnych. Korzystając z Twierdzenia 2.7 otrzymujemy błąd kwadratury Gaussa postaci: R(f ) =
Rb 1 2 (2N +2) p(x)ωN (ξ)dx = +1 (x)f (2N + 2)! a
(2.3.3) Rb 1 2 = f (2N +2) (ξ) p(x)ωN +1 (x)dx, (2N + 2)! a
przy ξ ∈ (a, b). 19
2.3.1
Kwadratury proste
2.3.1.1
Kwadratury Gaussa-Jacobiego
Kwadratury Gaussa-Jacobiego to kwadratury z wagą p(x) = (1 − x)α (1 + x)β ,
α, β > −1
(2.3.4)
(zob.[7, str.18]), w których węzły są pierwiastkami (N + 1)-ego wielomianu Jacobiego, postaci n (−1)n −α −β d Jn (x; α, β) = n (1 − x) (1 + x) [(1 − x)α+n (1 + x)β+n ]. n 2 · n! dx
Wielomiany Jacobiego tworzą układ ortogonalny (zob.[4, str.58]) w przedziale [-1,1], więc ograniczymy się do tego przedziału. W przypadku przedziału [a, b] dokonujemy podstawienia t= Zb a
gdzie g(x) = f
a+b 2
+
b−a x 2
a+b b−a + x, 2 2
1 b−a Z f (t)dt = g(x)dx, 2
(2.3.5)
−1
.
Współczynniki tej metody oraz błąd kwadratury możemy obliczyć przy pomocy wzorów (patrz [2, str.180]) Ak = −
R(f ) =
(2N + α + β + 4) · Γ(N + α + 2) · Γ(N + β + 2) · 2α+β (N + α + β + 2) · Γ(N + α + β + 2) · JN +1 (xk ; α, β) · JN +2 (xk ; α, β)(N + 2)!
Γ(N + α + 2) · Γ(N + β + 2) · Γ(N + α + β + 2) · (N + 1)!22N +α+β+3 (2N +2) f (ξ), (2N + α + β + 3)[Γ(2N + α + β + 3)]2 (2N + 2)!
gdzie ξ ∈ [−1, 1], a Γ oznacza funkcję gamma Eulera postaci Γ(x) =
Z∞
tx−1 e−t dt,
x > 0.
0
2.3.1.2
Kwadratury Gaussa-Legendre’a
Kwadraturami Gaussa-Legendre’a nazywamy kwadratury z wagą p(x) = 1, w których do wyznaczenia węzłów używamy wielomianów Legendre’a: P0 (x) = 1, 1 P2 (x) = (3x2 − 1), 2
1 P3 (x) = (5x3 − 3x), 2 20
··· ,
Pn (x) =
P1 (x) = x,
1 dn 2 (x − 1)n , 2n · n! dxn
···
Wielomiany Legendre’a, podobnie jak Jacobiego tworzą bazę ortogonalną w przedziale [−1, 1], więc tak jak przy poprzedniej metodzie dla przedziału [a, b] wykonamy podstawienie (2.3.5). Wzór przybliżenia całki kwadraturą Gaussa-Legendre’a ma postać: S(f ) = gdzie tk =
a+b 2
+
b−a xk , 2
N b−a X Ak f (tk ), 2 k=0
a współczynniki określone są wzorem (patrz [2, str.179]): Ak = −
2 , (N + 2)PN +2 (xk )PN0 +1 (xk )
w którym xk są pierwiastkami wielomianu PN +1 (x). Błąd tej metody jest równy (patrz [2, str.179]) R(f ) =
(b − a)2N +3 ((N + 1)!)4 (2N +2) f (ξ), (2N + 3)((2N + 2)!)3
gdzie ξ ∈ (a, b). 2.3.1.3
Kwadratury Gaussa - Czebyszewa
W przypadku, gdy we wzorze (2.3.4) zajdzie α = β = −1/2, czyli waga będzie równa p(x) = (1 − x2 )−1/2 , to węzłami będą pierwiastki (N+1) wielomianu Czebyszewa (patrz [5, str.303]) Tn (x) = cos(n arccos x), a przybliżenie całki tą metodą nazwiemy kwadraturą Gaussa - Czebyszewa. Współczynniki i węzły tej kwadratury wynoszą odpowiednio (por. [2, str.181]): Ak =
π , N +1
xk = cos
(2k + 1)π 2N + 2
Kwadratura ma więc następującą postać: S(f ) =
N π X (2k + 1)π f (cos ). N + 1 k=0 2N + 2
Błąd tej metody jest równy R(f ) =
f (2(N +1)) (ξ) , 22N +1 (2N + 2)! π
21
ξ ∈ (−1, 1).
2.3.2
Kwadratury złożone
2.3.2.1
Złożona kwadratura Gaussa - Legendre’a
rzędu N = 1 ma postać S(f ) = h
m−1 X i=0
"
!
!#
1 1 f a + h(1 − √ + 2i) + f a + h(1 + √ + 2i) 3 3
gdzie odległość między węzłami jest równa h =
,
H b−a = . 2 2m
Błąd tej kwadratury wynosi: (b − a)5 (4) R(f ) = f (ξ), 4320 · m4 gdzie ξ ∈ (a, b).
2.4
Ekstrapolacja Richardsona oraz jej zastosowanie w metodzie Romberga
Ekstrapolacja Richardsona polega na obliczaniu dokładniejszej wartości całki na podstawie dwóch wartości wyznaczonych numerycznie dla różnych wielkości podprzedziałów, na jakie podzielony jest przedział całkowania. Ekstrapolacja odgrywa kluczową rolę w jednej z metod całkowania numerycznego - metodzie Romberga, której wynikiem jest tablica coraz lepszych przybliżeń całki. W tej metodzie wykorzystany jest wzór złożonej metody trapezów. Stosujemy go, dzieląc przedział [a, b] na 2m podprzedziałów równej długości (m ≤ 0). Mamy Tm,0 =
m −1 2X
i=0
H (f (wi ) + f (wi+1 )), 2
gdzie wi = a + iH, i = 0, 1, · · · , m − 1, a odległość między węzłami H =
b−a . 2m
Najprostsze wyrażenie mamy dla m = 0 i jest równe 1 T0,0 = [f (a) + f (b)]. 2 Sumy T1,0 , T2,0 , . . . możemy obliczać rekurencyjnie: m−1
Tm,0
2X 1 = Tm−1,0 + Hf (a + (2i − 1)H). 2 i=1
Korzystając z faktu, że dla m > 0 wartości funkcji f , potrzebne do obliczenia Tm,0 , występują w Tm−1,0 , unikamy wielokrotnego obliczania wartości funkcji f w tych samych punktach. 22
Przybliżenia Tm,0 tworzą pierwszą kolumnę tablicy Romberga: T0,0 T1,0 T1,1 T2,0 T2,1 T2,2 T3,0 T3,1 T3,2 T3,3 ... Kolejne kolumny obliczamy korzystając z następującego wzoru rekurencyjnego: Tm,l = Tm,l−1 +
4l
1 (Tm,l−1 − Tm−1,l−1 ), −1
dla l > 0 (patrz [5, str.478]). Skuteczność ekstrapolacji w tej metodzie zależy od tego, jak dokładna jest lokalna aproksymacja funkcji całkowanej za pomocą wielomianu. Nierzadko się zdarza, że funkcja podcałkowa jest w pewnym punkcie nieskończona lub pewna pochodna niskiego rzędu tej funkcji w jakimś punkcie przedziału całkowania jest nieskończona, co skutkuje tym, że jeden składnik opisanych wzorów, związany z punktem, gdzie np. pochodna funkcji podcałkowej jest nieskończona, wnosi większy błąd niż wszystkie inne składniki. Często warto zbadać możliwość przekształcenia lub modyfikacji danego zadania po to, aby nowa funkcja podcałkowa była bardziej wygodna w całkowaniu numerycznym.
23
Rozdział 3 Porównanie wybranych metod numerycznych W poprzednim rozdziale omówiliśmy najważniejsze rodzaje kwadratur. W tym rozdziale chcemy poszukać odpowiedzi na pytanie: Która z tak wielu kwadratur jest najlepsza do rozwiązania jednego konkretnego zadania? Odpowiedź ta bardzo rzadko może być jednoznaczna i w pełni zadowalająca. Przede wszystkim musimy określić, w jakim sensie kwadratura ma być najlepsza, czy mamy na myśli minimalizację kosztu uzyskania przybliżenia całki z zadaną dokładnością, czy też minimalizację błędu kwadratury o ustalonych węzłach itp. Wydaje się jednak, że poszukiwanie optymalnej kwadratury dla pojedynczego zadania jest niecelowe, gdyż zarówno stopień trudności teoretycznych, jakie musi pokonać osoba poszukująca rozwiązania, jak i koszt ewentualnie potrzebnych obliczeń numerycznych, by wyznaczyć węzły i współczynniki takiej kwadratury są na ogół niewspółmiernie duże do uzyskanych korzyści. Z tych względów skupimy się na kwadraturach optymalnych w klasie funkcji. Wybór stosowanej kwadratury w danej klasie funkcji powinien zależeć zarówno od żądanej dokładności, jak i od dostępnej informacji o całkowanej funkcji - w praktyce informacja ta dotyczy najczęściej jej regularności, np. ciągłości i ograniczoności pochodnych niskiego rzędu. Poniżej zajmiemy się porównywaniem kwadratur, by odnaleźć kwadratury optymalne w sensie minimalizacji reszty (por. [8]).
3.1 3.1.1
Metody całkowania numerycznego Newtona-Cotesa Porównanie metod prostych Newtona-Cotesa
Metoda prostokątów jest obarczona dość dużym błędem, ponieważ prostokąty niezbyt dobrze przybliżają pole pod wykresem funkcji. Zaletą metody jest prosty wzór wyliczania 24
całki. Metoda trapezów jest dokładniejsza od metody prostokątów, gdyż trapezy dużo lepiej przybliżają pole pod wykresem funkcji. W praktyce oznacza to też mniejszą ilość obliczeń w celu uzyskania porównywalnej dokładności wyniku. Parabole przybliżają wykres funkcji z małym błędem, stąd metoda Simpsona jest najdokładniejszą metodą wyznaczania wartości całek oznaczonych z tutaj opisanych. Dokładność okupiona jest nieco skomplikowanym wzorem obliczeniowym. Możliwości podnoszenia dokładności kwadratury wiążą się oczywiście ze zwiększaniem liczby węzłów w przedziale całkowania, czyli podnoszeniem stopnia wielomianu interpolującego. Istotnym ograniczeniem jest tu występowanie efektu Rungego (por.[9]), czyli pogorszenia jakości interpolacji wielomianowej, mimo zwiększenia liczby jej węzłów. Jak opisaliśmy powyżej, początkowo ze wzrostem liczby węzłów przybliżenie poprawia się, jednak po dalszym wzroście węzłów zaczyna się pogarszać, co jest szczególnie widoczne na końcach przedziałów.
Rysunek 3.1: Efekt Rungiego Ponieważ zgodnie z twierdzeniem Weierstrassa istnieje ciąg wielomianów interpolujących coraz większych stopni, które przybliżają jednostajnie funkcje ciągłą, można uważać to za paradoks, iż efekt Rungego ma dokładnie odwrotny wynik. Jest to spowodowane nałożeniem warunku na równoodległość węzłów. Takie zachowanie się wielomianu interpolującego występuje nie tylko przy interpolacji za pomocą wielomianów wysokich stopni przy stałych odległościach węzłów, ale również w sytuacji, jeśli interpolowana funkcja jest 25
nieciągła albo odbiega znacząco od funkcji gładkiej. Aby uniknąć tego efektu, używa się kwadratur względnie niskiego (nie wyżej niż czwarty do siódmego) stopnia lub stosuje się interpolację z węzłami umieszczonymi coraz gęściej na końcach przedziału interpolacji. Najbardziej praktycznym rozwiązaniem jest wykorzystanie addytywności całki, podzielenie przedziału całkowania na wiele części i stosowanie w każdej z nich kwadratury niskiego rzędu (czyli zastosowanie kwadratury złożonej).
3.1.2
Porównanie metod złożonych Newtona-Cotesa
W poprzednim rozdziale opisaliśmy metody proste i złożone Newtona - Cotesa. Mogliśmy zauważyć, że zwiększając liczbę podprzedziałów, możemy dowolnie zmniejszyć błąd całkowania. Ogólnie, jeśli pochodna występująca w reszcie jest ograniczona w przedziale całkowania, to w granicy, kiedy liczba podprzedziałów dąży do nieskończoności, przybliżenie musi być zbieżne do dokładnej wartości całki. Wynika stąd, że stosowanie metod złożonych daje dużo lepsze efekty. Stosując złożone wzory zamknięte Newtona - Cotesa uzyskuje się jeszcze jedną korzyść. Ponieważ końce każdego podprzedziału (z wyjątkiem końców całego przedziału całkowania) są węzłami należącymi do dwóch podprzedziałów, więc całkowita liczba użytych węzłów nie jest równa iloczynowi liczby podprzedziałów m przez liczbę punktów w każdym podprzedziale N + 1, lecz jest od niego mniejsza o m − 1. W przypadku wzorów otwartych, ze względu na to, że punkty końcowe nie są węzłami, nie ma redukcji ostatecznej liczby węzłów. Zaletą wzorów złożonych niskiego rzędu, a w szczególności wzoru trapezów, jest to, że współczynniki mało się wahają, co ma dobry wpływ na właściwości zaokrąglenia. Jeśli całkę aproksymuje się wzorem złożonym, w którym m jest duże, to błąd zaokrąglenia może mieć decydujący wpływ na wynik. Ponieważ oczekiwany błąd zaokrąglenia jest najmniejszy wtedy, gdy współczynniki są równe, więc wzór trapezów ma prawie idealne właściwości zaokrąglenia, do czego wrócimy jeszcze przy metodzie Romberga. Znacznie gorsze właściwości mają wzory Newtona Cotesa wyższych rzędów, w których współczynniki mogą być nawet ujemne. Porównajmy więc metody złożone niskiego rzędu. W tym celu omówimy następujące przykłady: Przykład 1. Mamy za zadanie policzyć całkę Z1 √
1 + x dx.
0
Przyjrzyjmy się funkcji podcałkowej w układzie współrzędnych: 26
Obliczając tę całkę analitycznie otrzymujemy Z1 √ 0
2 √ 1 + x dx = (2 2 − 1) ≈ 1.21895. 3
Zobaczmy, jakie wyniki otrzymamy, stosując poznane przez nas kwadratury złożone Newtona - Cotesa: metodę prostokątów, trapezów oraz parabol. Porównajmy rezultaty dla różnej ilości podprzedziałów m w poniższej tabeli:
Widzimy, że wynik zbliżony do otrzymanego analitycznie otrzymujemy metodą parabol już przy trzech podprzedziałach. Pozostałe metody przy zwiększającej się ilości 27
podprzedziałów dają coraz lepsze przybliżenia naszej całki. Przykład 2. Policzmy całkę Z4 −2
1 − x4 + 2x3 + 2x2 − 5x + 10 dx. 2
Funkcję podcałkową, czyli funkcję wielomianową czwartego stopnia przedstawia następujący wykres:
Analityczne obliczenia całki prowadzą do wyniku Z4 −2
1 2 − x4 + 2x3 + 2x2 − 5x + 10 dx = 92 = 92.4. 2 5
Wyniki numeryczne widzimy w poniższej tabeli:
28
W przypadku tej funkcji obliczenia już nie są tak proste jak w Przykładzie 1. Dokładne wyniki daje, podobnie jak w poprzednim przykładzie, metoda parabol, ale dopiero przy 68 podprzedziałach, z których każdy wynosi ok. 0.08824. Widzimy, że nawet przy długości podprzedziału H = 0.012 żadna z pozostałych metod nie daje jeszcze dokładnego rozwiązania. Przykład 3. Rozpatrzmy całkę zadaną wzorem π
Z2
xsinx dx.
− π2
Powyższą funkcję podcałkową ilustruje rysunek:
29
Po wyliczeniu analitycznie tej całki otrzymujemy, że π
Z2
xsinx dx = 2.
− π2
By obliczyć numerycznie tę całkę skorzystamy z tych samych metod, co poprzednio. Spójrzmy na wyniki:
30
Tak jak w poprzednich przykładach widzimy dużą skuteczność metody parabol, która daje dokładny wynik już przy 12 podprzedziałach. Przyjrzyjmy się dwóm pozostałym metodom. Przy jednym podprzedziale mamy do czynienia z metodami prostymi prostokątów i trapezów. Widzimy, jakim błędem jest obciążona metoda prostokątów, która daje wynik 0. Możemy więc zaobserwować większą skuteczność metody prostej trapezów, co omawialiśmy wcześniej. Jednakże przy większej ilości podprzedziałów widzimy większą efektywność metody prostokątów, która przy łatwiejszym wzorze obliczeniowych daje lepsze rezultaty niż z pozoru skuteczniejsza metoda trapezów. Metoda prostokątów daje dokładny wynik przy 395 podprzedziałach, a metoda trapezów dopiero przy 626 podprzedziałach. Podsumowując widzimy, iż preferowanymi metodami całkowania numerycznego powinny być metoda parabol oraz metoda prostokątów. Druga ma stosunkowo prosty wzór wyliczeniowy. W metodzie parabol wzór jest bardziej skomplikowany, lecz ze względu na jej dokładność wykonamy mniej obliczeń, zatem szybciej uzyskamy wynik i mniejsze będą błędy zaokrągleń.
3.2
Porównanie dokładności metod Newtona-Cotesa i Gaussa
Kwadratur Gaussa rzadko używano w praktycznych obliczeniach przed wprowadzeniem maszyn cyfrowych. Przyczyna tkwiła w aspekcie praktycznym, gdyż w obliczeniach na arytrometrze łatwiej było rachować na liczbach całkowitych lub wymiernych niż na wie31
locyfrowych liczbach występujących we wzorach Gaussa. W takich obliczeniach wartości funkcji często odczytywało się z tablic, co przy prostych argumentach usuwa potrzebę interpolacji. Natomiast w sytuacji, gdy obliczenia wykonuje maszyna cyfrowa, wartości funkcji są obliczane najczęściej za pomocą przybliżeń wielomianowych lub wymiernych. W takim przypadku oraz wtedy, gdy funkcja, którą całkujemy, nie jest stablicowana, ale określona analitycznie, musimy brać pod uwagę nie tylko wzory Newtona-Cotesa, ale i wzory Gaussa, planując wykonanie obliczeń na maszynie cyfrowej. Jeśli mamy zastosować kwadraturę prostą (być może wysokiego rzędu) to dobrze byłoby wybrać wzór Gaussa. Wynika to stąd, ze wzory Gaussa dają większą dokładność przy mniejszej liczbie węzłów niż wzory Newtona-Cotesa - ten sam rząd r = 2N + 2 mają kwadratury Gaussa dla (N + 1) węzłów, a kwadratury Newtona-Cotesa dla (2N + 1) węzłów. Ponadto współczynniki w przypadku kwadratur Gaussa są co do modułu nieco mniejsze. Wobec tego wzory Gaussa nie tylko wymagają obliczenia mniej wartości f (x), lecz także ich reszty są nieco korzystniejsze. W tym przypadku ogólna wyższość wzorów Gaussa jest zatem oczywista. Prócz tego, w przeciwieństwie do wzorów Newtona-Cotesa, stosując dla funkcji ciągłej wzór Gaussa wystarczająco wysokiego rzędu możemy osiągnąć dowolną ustaloną dokładność (patrz [6, str.136]). Stosowanie ciągu wzorów Gaussa przy rosnącym N daje zbieżny ciąg przybliżeń. Mimo to zapamiętywanie w maszynie cyfrowej węzłów i współczynników ciągu wzorów Gaussa jest niewygodne oraz trudno szacować resztę dla wszystkich, nawet niedużych, wartości N , a wzory złożone są skuteczne, więc takie podejście jest w praktyce rzadko stosowane. Jeśli obliczając całkę z funkcji danej analitycznie, używamy wzoru złożonego dla ciągu 2, 4, 8, 16, ... podprzedziałów, to porównując wzory Gaussa i Newtona-Cotesa, nie widzimy już wyższości tych pierwszych. Porównanie reszt pozostaje słuszne. Wtedy jednak we wzorze Newtona-Cotesa węzły przy m podprzedziałach są także węzłami przy 2k m podprzedziałach (k = 1, 2, ...), a nie jest tak we wzorach Gaussa (gdy w każdym podprzedziale używa się co najmniej dwóch węzłów). Dlatego z numerycznego punktu widzenia wzory Newtona-Cotesa są efektywniejsze, gdyż w obliczeniach przy numerycznym całkowaniu większość czasu traci się na obliczanie wartości funkcji f (x). Można udowodnić, że jeśli stosuje się wzór złożony przy 1, 2, 4, ..., 2m podprzedziałach, a w każdym podprzedziale używa się N węzłów, to całkowita liczba różnych węzłów jest równa: • N + (N − 1)(2m − 1) we wzorach zamkniętych Newtona-Cotesa
32
• 2m N we wzorach otwartych Newtona-Cotesa (dla N parzystych) • N + (N + 1)(2m − 1) we wzorach otwartych Newtona-Cotesa (dla N nieparzystych) • N (2m+1 − 1) we wzorach Gaussa.
Powyższa tabela, która zawiera te liczby dla kilku wybranych wartości m i N , potwierdza przewagę wzorów zamkniętych Newtona- Cotesa. W większości przypadków te korzyści obliczeniowe odgrywają większą rolę niż zaleta wyższej dokładności wzorów Gaussa (patrz [6, 138]).
3.3
Metoda Romberga
Wyznaczanie brakujących wartości poza zakresem węzłów jest z reguły obarczone większymi błędami niż uzupełnianie wartości wewnątrz zakresu. Efekty przy wyjściu poza węzły dla wielomianu wysokiego stopnia są podobne do efektu Rungego. Jednakże ekstrapolacja w niewielkiej odległości od węzła może dawać przydatne wartości. Metoda Romberga należy obecnie do najczęściej stosowanych metod całkowania numerycznego, między innymi dlatego, że daje prostą strategię automatycznego wyboru właściwej długości kroku. Zaczyna się od kroku dość dużego, zmniejsza się go dwukrotnie tyle razy, aby ekstrapolacja dała wartości dostatecznie bliskie (por. [3, str.284]). Porównując różne metody całkowania - kwadratury Newtona-Cotesa, metodę ekstrapolacji, metodę Gaussa, stwierdzamy, że przy równym nakładzie obliczeń, zmierzonym 33
jako liczba obliczeń wartości funkcji, najdokładniejsze wyniki otrzymujemy stosując metodę Gaussa. Jeśli wiedzielibyśmy dla danej całki i żądanej dokładności ε, jakie należy przyjąć N , aby aproksymować całkę z dokładnością ε za pomocą wzoru Gaussa do wskaźnika N , to metoda Gaussa z pewnością przewyższałaby inne metody. Niestety w większości przypadków jest to niemożliwe, aby za pomocą wzoru na błąd (2.3.3) znaleźć odpowiednie N , ponieważ najczęściej trudno jest podać rozsądne oszacowanie dla f (2N +2) (ξ) (dla ξ ∈ (a, b)). Aby oszacować błąd kwadratury Gaussa dla wskaźnika N , tę samą całkę aproksymujemy za pomocą dalszej kwadratury Gaussa do wskaźnika N + 1 i badamy, czy tak otrzymane dwie wartości przybliżone różnią się o mniej niż żądaną dokładność ε. Jeśli to nie zachodzi, to powtarzamy postępowanie dla wskaźnika N +2, N +3, . . . tak długo, aż dwie kolejne wartości przybliżone nie zmieniają się istotnie (por. [10, str.121]). Ponieważ w przeciwieństwie do metod ekstrapolacji przy przejściu od N do N + 1 nie można dalej wykorzystywać N obliczonych dotąd wartości funkcji f (xi ), więc zalety metody Gaussa szybko znikają. Przy kilkukrotnym zastosowaniu kwadratur Gaussa dla różnych N nakład obliczeń będzie tak samo duży lub nawet większy niż w metodach ekstrapolacji. Widzimy więc większą efektywność algorytmów ekstrapolacji, dlatego są one w praktyce częściej stosowane. Przykład 4. Policzmy numerycznie całkę Z2 1
1 dx = ln 2 ≈ 0, 69314718. x
Porównamy wyniki obliczeń powyższej całki za pomocą wzorów złożonych: wzoru trapezów, wzoru parabol, wzoru Gaussa-Legendre’a oraz metody Romberga dla podziału na m = 2i , gdzie i = 1, . . . , 7, równych części przedziału [1; 2]. Wyniki numeryczne widzimy w tabeli:
34
Widzimy,że drugi i trzeci wzór złożony jest szybciej zbieżny od pierwszego. Metoda Romberga jest także szybciej zbieżna niż złożony wzór trapezów. Jeżeli wykonujemy obliczenia kolejno dla podziału na m = 21 , 22 , . . . równych części przedziału całkowania, to w przypadku wzorów możemy wykorzystać fakt, że węzły z poprzedniego podziału są także węzłami w następnym podziale, natomiast w przypadku wzorów tak nie jest. Z tego powodu wymienione dwie pierwsze kwadratury wymagają mniejszego nakładu obliczeń. Jeżeli funkcja podcałkowa nie jest dana analitycznie, lecz jest określona za pomocą tablicy (zakładamy, że argumenty są równoodległe, ponieważ najczęściej w praktyce tak jest), to prawie na pewno stosujemy wzory Newtona-Cotesa. Zastosowanie wzorów Gaussa wymagałoby przeprowadzenia interpolacji. Przy wyborze wzoru całkowania kierujemy się takimi samymi przesłankami jak w przypadku funkcji danych analitycznie [2, 186].
35
3.4
Implementacja metody
Zaprezentujemy teraz implementację złożonych metod Newtona-Cotesa wykonaną w języku programowania C#, która została wykorzystana przy porównaniu tych metod w podrozdziale 3.1.2. Wykorzystamy dane z przykładu 3. namespace Calkowanie { class NewtonCotes { static void Main(string[] args) { Console.WriteLine("**ZŁOŻONE METODY NEWTONA-COTESA**"); double a = -0.5 * Math.PI; //lewy koniec przedziału [a,b] double b = 0.5 * Math.PI; //prawy koniec przedziału [a,b] int m = PobierzLiczbePodprzedzialow(); //ilość podprzedziałów m przedziału [a,b] int metoda = PobierzMetode(); //wybrana złożona metoda Newtona - Cotesa if (metoda == 1 || metoda == 3) m = 2 * m; //podstawienie dla większej czytelności zapisu //(por. 2.1.2.2) double H = Math.Abs(b - a) / m; //odległość między kolejnymi węzłami, //(por. wzór (...)) double[] f = WyznaczWartosci(a, b, m, metoda, H); //tablica, w której są przechowywane wartości //funkcji w węzłach ObliczKwadrature(f, H, m, metoda); //oblicza przybliżoną wartość całki w [a,b] }
36
static double Funkcja(double x1) { //zwraca funkcję podcałkową, której wartość numeryczną chcemy otrzymać return x1 * Math.Sin(x1); }
static int PobierzLiczbePodprzedzialow() { /* pobiera od użytkownika ilość podprzedziałów m, na które chcemy podzielić przedział [a,b] */ int m = 0; //ilość podprzedziałów double zmp; //zmienna pomocnicza do { Console.WriteLine("Podaj ilość podprzedziałów: "); zmp = Convert.ToDouble(Console.ReadLine()); m = (int)zmp; if (m < 1 || Math.Abs(zmp) - Math.Abs(m) != 0) //wymagane jest, aby m było liczbą naturalną Console.WriteLine("BŁĘDNA DANA"); } while (m < 1 || (Math.Abs(zmp) - Math.Abs(m) != 0)); return (m); }
static int PobierzMetode() {
//pobiera od użytkownika jedną z trzech dostępnych metod
int m = 0; //numer wybranej metody char typ; //zmienna, która przechowa wpisaną przez //użytkownika daną 37
do { Console.WriteLine("Wybierz metodę: "); Console.WriteLine(" 1 - Złożona metoda prostokątów"); Console.WriteLine(" 2 - Złożona metoda trapezów"); Console.WriteLine("3 - Złożona metoda parabol"); typ = Convert.ToChar(Console.ReadLine()); switch (typ) { case ’1’: m = 1; break; //gdy została wybrana złożona metoda prostokątów case ’2’: m = 2; break; //gdy została wybrana złożona metoda trapezów case ’3’: m = 3; break; //gdy została wybrana złożona metoda parabol default: m = 0; break; //w innym przypadku } } while (m == 0); return m; }
static double[] WyznaczWartosci(double a, double b, int m, int metoda, double H) {
/* wyznacza argumenty, będące węzłami kwadratury, a następnie oblicza
wartości funkcji w węzłach; parametry wejściowe: a,b - końce przedziału [a,b], m - ilość podprzedziałów, metoda - numer wybranej metody, H - odległość między węzłami */ double[] x = new double[m + 1]; //tablica, przechowująca równoodległe węzły kwadratury double[] f = new double[m + 1]; //tablica, przechowująca wartości funkcji podcałkowej 38
//w poszczególnych węzłach if (metoda == 1) { //dla metod otwartych Newtona-Cotesa x[0] = a + H; //pierwszym węzłem jest środek pierwszego podprzedziału for (int i = 2; i <= m; i += 2) {
//odległość między środkami podprzedziałów jest równa //2H, nie bierzemy pod uwagę końców podprzedziałów x[i] = Przyblizenie(x[0] + i * H); //Przyblizenie - funkcja pomocnicza, przybliżająca //liczbę do piątego miejsca po przecinku x[i - 1] = 0;
}
for (int i = 0; i < m; i += 2) { f[i] = Funkcja(x[i]); //Funkcja - funkcja pomocnicza, obliczająca wartość //funkcji podcałkowej dla argumentu x[i] f[i] = Przyblizenie(f[i]); //przybliżenie wartości funkcji if (f[i] < 0) f[i] = Math.Abs(f[i]); //wartość bezwzględna wartości funkcji } } if (metoda == 2 || metoda == 3) {
//dla metod zamkniętych Newtona-Cotesa x[0] = a; //pierwszym węzłem jest lewy koniec przedziału for (int i = 1; i <= m; i++) x[i] = Przyblizenie(x[0] + i * H); //Przyblizenie - funkcja pomocnicza, przybliżająca //liczbę do piątego miejsca po przecinku
39
for (int i = 0; i < m + 1; i++) { f[i] = Funkcja(x[i]); //Funkcja - funkcja pomocnicza, obliczająca wartość //funkcji podcałkowej dla argumentu x[i] f[i] = Przyblizenie(f[i]); //przybliżenie wartości funkcji if (f[i] < 0) f[i] = Math.Abs(f[i]); //wartość bezwzględna wartości funkcji } } return f; //wynikiem są wartości funkcji podcałkowej w węzłach }
static void ObliczKwadrature(double[] f, double H, int m, int metoda) {
/* funkcja oblicza przybliżoną wartość całki dla parametrów: f - wartości funkcji w węzłach, H - odległość między węzłami, m - liczba podprzedziałów, metoda - numer wybranej metody */
Console.Write("Przybliżona wartość całki jest równa: "); switch (metoda) { case 1: Console.WriteLine(Kwadratura1(f, H, m)); break; //funkcja pomocnicza, która oblicza przybliżoną //wartość całki złożoną metodą prostokątów case 2: Console.WriteLine(Kwadratura2(f, H, m)); break; //funkcja pomocnicza, która oblicza przybliżoną //wartość całki złożoną metodą trapezów case 3: Console.WriteLine(Kwadratura3(f, H, m)); break; //funkcja pomocnicza, która oblicza przybliżoną //wartość całki złożoną metodą parabol default: Console.WriteLine("Błąd programu."); break; 40
} Console.WriteLine("Dziękuję za skorzystanie z programu."); System.Threading.Thread.Sleep(10000); }
static double Kwadratura1(double[] f, double H, int m) {
/* oblicza przybliżoną wartość całki złożoną metodą prostokątów; w węzłach, parametry wejściowe: f - tablica przechowująca wartości funkcji, H - odległość między węzłami, m - ilość podprzedziałów; wartość funkcji: S - kwadratura */
double S = 0; //przybliżenie całki (kwadratura) for (int i = 0; i <= m; i += 2) S += 2.0 * H * f[i]; //obliczenie przybliżenia całki przy pomocy //wzoru (...) return Przyblizenie(S); //wynikiem jest przybliżenie S do 5. miejsca //po przecinku }
static double Kwadratura2(double[] f, double H, int m) {
/* oblicza przybliżoną wartość całki złożoną metodą trapezów; parametry wejściowe: f - tablica przechowująca wartości funkcji w węzłach, H - odległość między węzłami, m - ilość podprzedziałów; wartość funkcji: S - kwadratura */
double S = 0; //przybliżenie całki (kwadratura) for (int i = 0; i < m; i++) S += 0.5 * H * (f[i] + f[i + 1]); //obliczenie przybliżenia całki przy pomocy //wzoru (...) 41
return Przyblizenie(S); //wynikiem jest przybliżenie S do 5. miejsca //po przecinku }
static double Kwadratura3(double[] f, double H, int m) {
/* oblicza przybliżoną wartość całki złożoną metodą parabol; parametry wejściowe: f - tablica przechowująca wartości funkcji w węzłach, H - odległość między węzłami, m - ilość podprzedziałów; wartość funkcji: S - kwadratura */
double S = 0; //przybliżenie całki (kwadratura) for (int i = 0; i <= m - 2; i += 2) S += 0.33333333 * H * (f[i] + 4 * f[i + 1] + f[i + 2]); //obliczenie przybliżenia całki przy pomocy //wzoru (...) return Przyblizenie(S); //wynikiem jest przybliżenie S do 5. miejsca //po przecinku }
static double Przyblizenie(double liczba) {
/* funkcja pomocnicza, obliczająca przybliżenie danej liczby do 5.miejsca po przecinku; parametr: liczba - zmienna typu double */
liczba *= 100000.0; if (liczba < 0) liczba = Math.Ceiling(liczba - 0.5); else liczba = Math.Floor(liczba + 0.5); return liczba / 100000.0; } } } 42
Zakończenie ...
43
Bibliografia [1] G. M. Fichtenholz, Rachunek różniczkowy i całkowy. Tom II, Wydawnictwo Naukowe PWN, Warszawa 1980 [2] Z. Fortuna, B. Macukow, J. Wąsowski, Metody numeryczne, Wydawnictwa NaukowoTechniczne, Warszawa 2002 [3] G. Dahlquist, A. Bjorck, Metody numeryczne, Wydawnictwo Naukowe PWN, Warszawa 1987 [4] G. Szegö, Orthogonal Polynomials, American Mathematical Society, 1939 [5] D. Kincaid, W. Cheney, Analiza numeryczna, Wydawnictwa Naukowo-Techniczne, Warszawa 2006 [6] A. Ralston, Wstęp do analizy numerycznej, Wydawnictwo Naukowe PWN, Warszawa 1983 [7] A. H. Stroud, D. Secrest, Gaussian quadrature formulas, Prentice-Hall, London 1966 [8] J. Jankowska, M. Jankowski, Przegląd metod i algorytmów numerycznych. Część 1, Wydawnictwo Naukowo-Techniczne, Warszawa 1981 [9] Philip J. Davis, Interpolation and Approximation, Blaisdell Publishing Company, Waltham, Massachusetts 1963 [10] J. Stoer, Wstęp do metod numerycznych, Wydawnictwo Naukowe PWN, Warszawa 1979
44
