Całkowanie numeryczne przy użyciu kwadratur
Plan wykładu:
1. Kwadratury Newtona-Cotesa a) wzory: trapezów, parabol etc.
b) kwadratury złożone 2. Ekstrapolacja
a) Ekstrapolacja Richardsona b) Metoda Romberga
c) Metody adaptacyjne
3. Kwadratury Gaussa a)Gaussa-Legendre'a b)Gaussa-Hermitte'a c) Gaussa-Laguerre'a
4. Całkowanie funkcji wielu zmiennych
2 Całkowanie numeryczne oznacza zastosowanie metod numerycznych w celu
wyznaczenia przybliżonej wartości całki oznaczonej.
Skoro funkcję podcałkową możemy interpolować
to wielomian interpolacyjny można wykorzystać do całkowania.
Dla danego ciągu wartości funkcji podcałkowej f(x0), f(x1),...,f(xN) definiujemy wielomian interpolacyjny Lagrange'a:
'(x) = L
N(x) = X
N k=0©
k(x)F (x
k)
©
k(x) = (x ¡ x
0)(x ¡ x
1) : : : (x ¡ x
k¡1)(x ¡ x
k+1) : : : (x ¡ x
N)
(x
k¡ x
0)(x
k¡ x
1) : : : (x
k¡ x
k¡1)(x
k¡ x
k+1) : : : (x
k¡ x
N) = Y
j=0 j6=k
x ¡ x
jx
k¡ x
jC =
Z
ba
f (x)dx
Podstawiamy wielomian interpolacyjny w miejsce funkcji podcałkowej:
Powyższe wzory definiują tzw. kwadraturę. Ak są współczynnikami kwadratur.
Jeśli spełniony jest warunek
To wtedy
Dokładność wyznaczonej wartości całki jest ograniczona dokładnością przybliżenia funkcji
Z
ba
F (x)dx ¼ Z
ba
'(x)dx = X
N k=0A
kF (x
k)
A
k= Z
ba
©
k(x)dx
jf(x) ¡ '(x)j < "; x 2 [a; b]
¯ ¯
¯ ¯
¯ Z
ba
F (x)dx ¡ X
N k=0A
kF (x
k)
¯ ¯
¯ ¯
¯ =
= ¯¯ ¯
¯ ¯ Z
ba
(F (x) ¡ '(x))dx ¯¯ ¯
¯ ¯ 6 "(b ¡ a)
Jeśli funkcja podcałkowa posiada osobliwości (np. jest nieograniczona, lub przedział całkowania jest
nieskończony) wówczas powyższy schemat
całkowania ulega modyfikacji – funkcję podcałkową zastępujemy iloczynem funkcji wagowej i nowej gładkiej funkcji:
Funkcja wagowa p(x) zawiera wszystkie osobliwości funkcji F(x) lub jej dobór wynika z zastosowanych wielomianów ortogonalnych:
F (x) = p(x)f (x)
A
0k= Z
ba
p(x)©
k(x)dx Z
ba
F (x)dx =
Z
b ap(x)f (x)dx ¼
¼
Z
b ap(x)'(x)dx = X
N k=0A
0kf (x
k)
4 Chcemy wyznaczyć wartość całki:
Stosując wzór
Powyższy wzór nosi nazwę kwadratury, a punkty x1,x2...,xN węzłami kwadratury.
Błąd przybliżenia całki kwadraturą (błąd metody):
Kryterium dokładności kwadratury można przyjąć zgodność I(W) z S(W), gdy W jest wielomianem.
Wówczas mówimy że dana kwadaratura jest rzędu r (r ≥ 1 ) jeśli
dla wszystkich wielomianów stopnia mniejszego niż r.
Kwadratura jest zbieżna dla każdej funkcji f∈ C([a,b]) wtedy gdy:
1) Jest ona zbieżna dla każdego wielomianu 2) Istnieje liczba M niezależna od N taka że
I(f ) = Z
ba
p(x)f (x)dx
S(f ) = X
N k=0A
kf (x
k); x 2 [a; b]
E(f ) = I(f ) ¡ S(f)
I(W ) = S(W )
B
N= X
N k=0jA
Nkj 6 M; N = 1; 2; : : :
Kwadratury Newtona-Cotesa
Rozważamy przypadek z węzłami równoodległymi xi=a+ih, i=0,1,2,...,N. Jeśli końce przedziału są również węzłami wówczas kwardatury noszą nazwę kwadratur zamkniętych.
Przybliżamy funkcję podcałkową wielomianem Lagrange'a stopnia conajwyżej N
Błąd przybliżenia (interpolacji)
Wprowadzamy nową zmienną t
f (x
i) = L
N(x
i); i = 0; 1; 2; : : : ; N L
N(x) =
X
N k=0f (x
k)©
k(x)
©
k(x) = Y
j=0 j6=k
x ¡ x
jx
k¡ x
jR
N +1(x) = f (x) ¡ L
N(x)
= 1
(N + 1)! !
N +1(x)f
(N +1)(»)
» 2 (a; b)
x = a + ht ©
k(t) = Y
j=0 j6=k
t ¡ j
k ¡ j
Ostatecznie otrzymujemy:
gdzie:
S(f ) = X
N k=0A
kf
kA
k= h ( ¡1)
N¡kk!(N ¡ k)!
Z
N0
t(t ¡ 1) : : : (t ¡ N) (t ¡ k)
f
k= f (a + kh) b ¡ a
Własności:
1) Gdy N jest nieparzyste wówczas kwadratura jest rzędu (N+1) (dokładna dla wielomianów stopnia N), dla parzystego N rząd kwadratury wynosi (N+2)
2) Jeżeli funkcja podcałkowa jest r-krotnie różniczkowalna, wówczas błąd metody można przedstawić w postaci:
współczynnik Cr nie zależy od f.
3) Dla dużych N oszacowanie błędu jest trudne ze względu na pochodne wysokich rzędów lub ze względu na numeryczne kasowanie się
współczynników Ak
4) Współczynniki Ak zależą od N. W szczególności (wzór na Ak) zachodzi
dlatego metoda kwadratur Newtona-Cotesa nie jest zbieżna w klasie funkcji ciągłych.
W praktyce przedział całkowania dzieli się na m
podprzedziałów. W każdym podprzedziale określa się N (N=1,2,3) i przeprowadza całkowanie. Taka
procedura prowadzi do uzyskania kwadratur złożonych.
E(f ) = C
rf
(r)(»); » 2 [a; b]
N
lim
!1jA
kj = 1 Z
ba
f (x)dx =
Z
b a'(x)dx =
=
X
N k=0f
kZ
ba
©
k(x)dx
= h X
N k=0f
kZ
N0
'
k(t)dt
= h X
N k=0f
kA
k6 Kwadratury dla N=1,2,...,6 (całkowanie w
podprzedziale)
N=1 (wzór trapezów)
Ze wzoru na błąd interpolacji wynika, że kwadratura jest N+1=2 rzędu, a dokładnie przybliża wielomian N=1 stopnia.
Zatem błąd wyznaczenia przybliżonej wartości całki wynosi
A
0= ¡h Z
10
(t ¡ 1)dt = 1 2 h A
1= h
Z
1 0tdt = 1 2 h S(f ) = 1
2 h(f
0+ f
1) h = b ¡ a
E(f ) = 1 2!
Z
ba
(x ¡ a)(x ¡ b)f
(2)(»)dx
= ¡ 1
12 h
3f
(2)(»); » 2 [a; b]
A
1= 4
3 h A
2= 1 3 h A
0= 1
3 h
h = b ¡ a 2
S(f ) = 1
3 h(f
0+ 4f
1+ f
2)
N=2 (wzór parabol – Simpsona)
Ponieważ N jest parzyste więc kwadratura jest
dokładna dla wielomianów stopnia N+1 i jest rzędu N+2. Dlaczego? Zgodnie z wzorem na błąd wzoru interpolacyjnego dostajemy
z powodu nieparzystości funkcji podcałkowej.
Dodajmy więc dodatkowy węzeł w x=(a+b)/2, który nie zmienia warunku interpolacji. Wówczas stopień wielomianu czynnikowego rośnie o 1:
(funkcja podcałkowa teraz jest parzysta)
» 2 [a; b]
E(f ) = f
(4)(»
1) 4!
Z
b a(x ¡ a) µ
x ¡ a + b 2
¶
2(x ¡ b)dx
= ¡ 1
90 h
5f
(4)(») E(f ) »
Z
b a(x ¡ a) µ
x ¡ a + b 2
¶
(x ¡ b)dx = 0
N w A0/w A1/w A2/w A3/w A4/w A5/w A6/w błąd wzór
1 (1/2)h 1 1 h3 (1/12) f(2)(» ) trapezów
2 (1/3)h 1 4 1 h5 (1/90) f(4)(») parabol
3 (3/8)h 1 3 3 1 h5 (3/80) f(4)(») 3/8
4 (4/90)h 7 32 12 32 7 h7 (8/945) f(6)(») Milne'a 5 (5/288)h 19 75 50 50 75 19 h7 (275/12096) f(6)(») --- 6 (6/840)h 41 216 27 272 27 216 41 h9 (9/1400) f(8)(») Weddle'a
Kwadratury złożone Newtona-Cotesa
Kwadratury wyższych rzędów są rzadko stosowane. Natomiast błąd kwadratur niższych rzędów jest proporcjonalny do długości przedziału całkowania w odpowiedniej potędze. Zatem niski rząd kwadratury może nie zapewiniać wymaganej dokładności.
Problemu tego można uniknąć, dzieląc przedział całkowania na m podprzedziałów, w których przeprowadza się całkowanie kwadaraturami niższych rzędów a wyniki całkowania sumuje się.
8 Wzór złożony trapezów
Przedział całkowania dzieli się na m poprzedziałów:
Gdzie
Zakładamy że
Błąd złożonego wzoru trapezów
h = b ¡ a m S(f ) =
m
X
¡1 k=01
2 h(f
k+ f
k+1)
= h µ 1
2 f
0+ f
1+ : : : + f
m¡1+ 1 2 f
m¶
f
k= f (a + k ¢ h) f 2 C
2([a; b])
»
k2 (a + kh; a + (k + 1)h) E(f ) = ¡ h
312
m
X
¡1 k=0f
(2)(»
k)
= ¡ (b ¡ a)
312m
2¢ 1
m
m
X
¡1 k=0f
(2)(»
k)
Wyraz
jest średnią arytmetyczną wartości drugiej pochodnej w przedziale całkowania.
Można więc zapisać
Błąd zależy od 3 potęgi długości przedziału. Ale zwiększając m można istotnie ograniczyć jego wartość.
1 m
m
X
¡1 k=0f
(2)(»
k) = f
(2)(»); » 2 [a; b]
E(f ) = ¡ (b ¡ a)
312m
2f
(2)(»)
Wzór złożony parabol.
Przedział całkowania [a,b] dzielimy na m podprzedziałów (m jest parzyste).
W podprzedziałach [a,a+2h],..., [a+(m-1)h,b]
stosuje się wzór parabol a wyniki cząstkowe sumuje:
Zakładamy
Ze względu na ciągłość pochodnej istnieje taki punkt że:
Wówczas błąd złożonego wzoru parabol wyraża się wzorem
S(f ) = h 3
m=2
X
k=1
(f
2k¡2+ 4f
2k¡1+ f
2k)
= h
3 [f
0+ f
m+ 2(f
2+ f
4+ : : : + f
m¡2) + 4(f
1+ f
3+ : : : + f
m¡1)]
f 2 C
4([a; b])
»
k2 (a + 2(k ¡ 1)h; a + 2kh)
2 m
m=2
X
k=1
f
(4)(»
k) = f
(4)(»)
E(f ) = ¡ (b ¡ a)
5180m
4f
(4)(»)
E(f ) = ¡ h
590
m=2
X
k=1
f
(4)(»
k)
10 Ekstrapolacja Richardsona
(przypadek dla różniczkowania ale zastosowanie ogólne) Rozwijamy funkcję f(x) w szereg Taylora w otoczeniu punktów
i odejmujemy od siebie oba wyrażenia
a następnie przegrupowujemy wyrazy by obliczyć pierwszą pochodną
x § h f (x + h) =
X
1 k=01
k! h
kf
(k)(x) f (x ¡ h) =
X
1 k=01
k! ( ¡1)
kh
kf
(k)(x) f (x + h) = f (x) + hf
(1)(x) + h
22 f
(2)(x) + h
36 f
(3)(x) + : : : f (x ¡ h) = f(x) ¡ hf
(1)(x) + h
22 f
(2)(x) ¡ h
36 f
(3)(x) + : : :
f (x + h) ¡ f(x ¡ h) = 2hf
(1)(x) + 2
3! h
3f
(3)(x) + 2
5! h
5f
(5)(x) + : : :
f
(1)(x) = f (x + h) ¡ f(x ¡ h)
2h ¡
· 1
3! h
2f
(3)(x) + 1
5! h
4f
(5)(x) + 1
7! h
6f
(7)(x) + O(h
8)
¸
Powyższa formuła w postaci ogólnej
co można interpretować jako przybliżenie f(1)(x).
Za h podstawimy h/2
i obliczmy różnicę
Zatem L1 przybliża f(1)(x) z dokładnością O(h4) (wyrazów rzędu h4).
Dokonujemy podstawienia
w L2
Ã(h) = 4 3 Á
µ h 2
¶
¡ 1
3 Á(h)
Podstawiając do L2
Otrzymujemy
Powtarzając M-krotnie powyższy proces dostaniemy coraz lepsze przybliżenie pierwszej pochodnej tzn.
dokładność jej przybliżenia jest na poziomie O(h2M).
(o ile h<<1).
L
h;2= Ã(h) + b
4h
4+ b
6h
6+ : : : L
h=2;1= Á
µ h 2
¶
+ a
2h
24 + a
4h
416 + a
6h
664 + : : :
L
h=2;2= Ã
µ h ¶
+ b
4h
4+ b
6h
6L
h;1= Á (h) + a
2h
2+ a
4h
4+ a
6h
6+ : : :
L
h;3= '(h) + c
6h
6+ c
8h
8+ : : : L
2= (4
1L
h=2;1¡ L
h;1)=(4
1¡ 1)
= 4 3 Á
µ h 2
¶
¡ 1
3 Á(h) ¡ a
4h
44 ¡ 5a
6h
616 ¡ : : :
L
3= (4
2L
h=2;2¡ L
h;2)=(4
2¡ 1)
= 16 15 Ã
µ h 2
¶
¡ 1
15 Ã(h) ¡ b
6h
620 Ã(h) ¡ : : :
'(h) = 16 15 Ã
µ h 2
¶
¡ 1
15 Ã(h)
12 Algorytm dla powyższej procedury jest następujący
1. Wybieramy h i liczymy
2. Następnie obliczamy
Obliczając rekurencyjnie wyrazy wg 2 dostajemy przybliżenia
D
n;0= L + O(h
2) D
n;1= L + O(h
4) D
n;2= L + O(h
6) D
n;3= L + O(h
8)
: : : : : : : : : : : :
D
n;k¡1= L + O(h
2k); h ! 0 D
n;k= 4
k4
k¡ 1 D
n;k¡1¡ 1
4
k¡ 1 D
n¡1;k¡1k = 1; 2; : : : ; M
n = k; k + 1; : : : ; M D
n;0= Á
µ h 2
n¶
; n = 0; 1; 2; : : : ; M
Algorytm ten definiuje tzw. ekstrapolację Richardsona. Generalnie jest to proces
rekurencyjnego wyznaczania pewnej wielkości (pochodnej, całki), co można zdefiniować przy pomocy wzoru
co w połączeniu z pkt. 2 daje szukane przybliżenie Dm,m.
Kolejne kroki algorytmu można zapisać w postaci tablicy
D
n;k¡1= L + X
1 j=kA
jkµ h 2
n¶
2jD
0;0D
1;0D
1;1D
2;0D
2;1D
2;2.. . .. . .. . . ..
D
M;0D
M;1D
M;2: : : D
M;M13 Metoda Romberga
Korzytamy z wzoru trapezów
Jeśli
to dla kolejnych wartości n dostajemy poniższy ciąg przybliżeń wartości całki
h = b ¡ a n
x 2 [0; 1]
Łatwo zauważyć, że do obliczenia T2n można wykorzystać już obliczone Tn
co ogólnie dla przedziału całkowania [a,b] można zapisać jako
S
n= h
Ã
nX
i=0
f (a + ih) ¡ f (a) + f (b) 2
!
S
0= 1
2 f (0) + 1
2 f (1) S
2= 1
4 f (0) + 1 2
½ f
µ 1 2
¶¾ + 1
4 f (1) S
4= 1
8 f (0) + 1 4
½ f
µ 1 4
¶ + f
µ 1 2
¶ + f
µ 3 4
¶¾ + 1
8 f (1) S
6= 1
16 f (0) + 1 8
½ f
µ 1 8
¶ + f
µ 1 4
¶ + f
µ 3 8
¶ + f
µ 1 2
¶ µ 5 ¶ µ
3 ¶ µ
7 ¶¾
1
S
2n= 1
2 S
2(n¡1)+ h
2nX
ni=1
f (a + (2i ¡ 1)h)
h
n= b ¡ a 2
nS
2= 1
2 S
0+ 1 2
½ f
µ 1 2
¶¾
S
4= 1
2 S
2+ 1 4
½ f
µ 1 4
¶ + f
µ 3 4
¶¾
S
6= 1
2 S
4+ 1 8
½ f
µ 1 8
¶ + f
µ 3 8
¶ + f
µ 5 8
¶ + f
µ 7 8
¶¾
14 W metodzie Romberga zakładamy, że odległość
między (n+1) węzłami wynosi
Do obliczenia całki wykorzystujemy rekurencyjną formułę z wzorem trapezów
Wartości kolejnych przybliżeń można uporządkować w postaci tablicy podobnie jak w przypadku
ekstrapolacji Richardsona.
Obliczenia przerywa się gdy spełniony jest warunek
lub po osiągnięciu zadanej liczby iteracji k.
Metoda Romberga jest przykładem kwadratury adaptacyjnej.
h
n= b ¡ a 2
njR
k;k¡ R
k¡1;k¡1j · "; " 2 R
Metody adaptacyjne
Liczymy numerycznie całkę np. wzorem parabol
Dzielimy przedział [a,b] na n podprzedziałów i stosujemy wzór parabol w każdym z nich
gdzie: ei jest lokalnym błędem przybliżenia wartości całki w i-tym podprzedziale [xi-1,xi].
Załóżmy że jego wartość możemy oszacować zgodnie z poniższym wzorem
Z
ba
f (x)dx = S(a; b) ¡ (b ¡ a)
590 f
(4)(»)
» 2 [a; b]
Z
ba
f (x)dx = X
ni=1
(S
i+ e
i)
je
ij · " x
i¡ x
i¡1b ¡ a R
0;0= 1
2 (b ¡ a) [f(a) + f (b)]
R
n;0= 1
2 R
n¡1;0+ b ¡ a 2
n2
X
n¡1 i=1f µ
a + (2i ¡ 1) b ¡ a 2
n¶
R
n;m= R
n;m¡1+ 4
mR
n;m¡1¡ R
n¡1;m¡14
m¡ 1
S(a; b) = b ¡ a 3
½
f (a) + 4f
µ a + b 2
¶
+ f (b)
¾
wówczas oszacowanie błędu całkowitego od góry jest następujące
co pozwala oszacować wartość bezwzględną błędu całki wyznaczonej numerycznie
Wniosek: przy założonej wartości , odpowiednio niski poziom błędu wartości całki osiągniemy zwiększając liczbę węzłów całkowania.
¯ ¯
¯ ¯
¯ X
ni=1
e
i¯ ¯
¯ ¯
¯ · X
ni=1
je
ij · "
b ¡ a X
ni=1
(x
i¡ x
i¡1) = "
¯ ¯
¯ ¯
¯ Z
ba
f (x)dx ¡ X
ni=1
S
i¯ ¯
¯ ¯
¯ · "
16 Kwadratury Gaussa
Nadal rozpatrujemy kwadratury typu:
ale nieco zmienimy metodologię postępowania.
Ustalamy funkcję wagową p(x)
oraz liczbę węzłów (N+1). Szukamy:
a)położenia węzłów b)współczynników Ak
tak aby rząd kwadratury był jak najwyższy.
Kwadratura tego typu nosi nazwę kwadratury Gaussa. Do wyznaczenia kwadratur Gaussa używa się wielomianów ortogonalnych.
Ciąg wielomianów
Nazywamy ortogonalnymi w przedziale [a,b] jeśli zachodzi pomiędzy nimi związek:
S(f ) = X
N k=0A
kf (x
k)
A
k= Z
ba
p(x)©
k(x)dx
Tw.1. Wielomiany ortogonalne mają tylko pierwiastki rzeczywiste, leżące w przedziale [a,b].
Tw.2. Nie istnieje kwadratura Gaussa rzędu wyższego niż 2(N+1). Kwadratura Gaussa jest rzędu 2(N+1) wtedy i tylko wtedy, gdy węzły xk są pierwiastkami wielomianu PN+1(x).
Tw. 3. Wszystkie współczynniki Ak w kwadraturach Gaussa są dodatnie.
Dlaczego rząd kwadratury Gaussa jest tak wysoki?
Musimy ustalić położenia N+1 węzłów oraz
współczynniki kombinacji liniowej N+1 wielomianów ortogonalnych. Daje to 2n+2 rząd.
Metoda kwadratur Gaussa jest zbieżna do każdej funkcji ciągłej w [a,b]. Kwadratury te są dokładne dla wielomianów stopnia 2N+1.
('
r; '
s) = Z
ba
p(x)'
r(x)'
s(x)dx = 0
r 6= s
f'
n(x) g = f'
0(x); '
1(x); : : : ; '
N(x) g
17 Korzystamy z tożsamości Christoffela-Darboux
βk – współczynnik stojący w wielomianie ϕk przy zmiennej w najwyższej potędze
Podstawmy za y zero wielomianu n-tego stopnia
Po wykonaniu mnożenia a następnie całkowania otrzymamy
X
n k=0'
k(x)'
k(y)
°
k= '
n+1(x)'
n(y) ¡ '
n(x)'
n+1(y)
®
n°
n(x ¡ y)
Korzystamy teraz z definicji wielomianu interpolacyjnego Lagrange'a
Wybieramy oczywiście przypadek taki że:
oraz korzystamy z faktu
f (x) = X
n j=0f (x
j)l
j(x)
l
j(x) = !
n(x)
(x ¡ a
j)!
n0(a
j)
y = d
j nX
¡1k=0
'
k(x)'
k(d
j)
°
k= ¡ '
n(x)'
n+1(d
j)
®
n°
n(x ¡ d
j) = ¢ p(x)'
0(x)
'
0(d
j)
°
0°
0= ¡ '
n+1(d
j)
®
n°
nZ
b ap(x) '
0(x)'
n(x) x ¡ d
jdx
!
n(x) = '
n(x)
'
0(x) = 1
®
k= ¯
k+1¯
k1 = ¡ '
n+1(d
j)
®
n°
nZ
b ap(x) '
n(x) x ¡ d
jdx
= ¡ '
n+1(d
j)'
0n(d
j)
®
n°
nZ
b ap(x)l
j(x)dx
= ¡ '
n+1(d
j)'
0n(d
j) A
°
k= Z
b ap(x)'
2k(x)dx
19 Współczynniki Ak:
Błąd kwadratury:
Węzły xk stanowią pierwiastki wielomianu PN+1(x).
(jak je znaleźć? => metody poszukiwania zer wielomianów)
Dla kwadratur niskiego rzędu węzły i współczynniki Ak są stablicowane. Aby zastosować wzory z przedziału [-1,1] w przedziale [a,b] należy dokonać transformacji liniowej zmiennej niezależnej:
P
0(x) = 1 P
1(x) = x
P
2(x) = 3x
2¡ 1 2
P
3(x) = 5x
3¡ 3x
22
A
k= ¡ 2
(N + 2)P
N +2(x
k)P
N +10(x
k)
E(f ) = 2
2N +3((N + 1)!)
4(2N + 3)((2N + 2)!)
3f
(2N +3)(»)
¡1 < » < 1
t = a + b
2 + b ¡ a 2 x
g(x) = f
µ a + b
2 + b ¡ a 2 x
¶ x 2 [¡1; 1]; t 2 [a; b]
Z
ba
f (t)dt = b ¡ a 2
Z
1¡1
g(x)dx
W praktyce nie używa się kwadratur wysokiego rzędu. O wiele lepszym rozwiązaniem jest
zastosowanie kwadratur złożonych tj. kwadratur niskiego rzędu w każdym podprzedziale a
wyniki sumuje się.
N k xk Ak
1 0, 1 (-/+)0.577350 1
2 0, 2
1 (-/+)0.774597
0 5/9
8/9 3 0, 3
1, 2 (-/+)0.861136
(-/+)0.339981 0.347855 0.652145
4 0,4 1, 3 2
(-/+)0.906180 (-/+)0.538469
0
0.236927 0.478629 0.568889
Z
ba
f (t)dt ¼ S(f) = b ¡ a 2
X
N k=0A
kf (t
k)
t
k= a + b
2 + b ¡ a 2 x
kPoszukiwanie zer wielomianów Legandre'a Mając ustalony stopień wielomianu (n) i numer jego zera (k), najpierw wyznaczamy jego przybliżone położenie (Tricomi)
gdzie:
a następnie iteracyjnie poprawiamy je metodą Newtona z zadaną dokładnością.
x
k=
½
1 ¡ n ¡ 1
8n
3¡ 1 384n
4µ
39 ¡ 28 sin
2Á
k¶¾
£cos(Á
k) + O(n
¡5)
Á
k= ¼(k ¡ 1=4)
n + 1=2
21 Kwadratury dla przedziału jedno- i obustronnie
nieskończonego
Kwadratura Gaussa-Laguerre'a
Ciąg wielomianów ortogonalnych stanowią wielomiany Laguerre'a:
[a; b] = [0; 1) p(x) = e
¡xL
n(x) = ( ¡1)
ne
xd
ndx
n(x
ne
¡x) L
0(x) = 1
L
1(x) = 1 ¡ x
L
2(x) = x
2¡ 4x + 2 2
L
3(x) = ¡x
3+ 9x
2¡ 18x + 6 6
Relacja rekurencyjna dla stowarzyszonych wielomianów Laguerre'a
wówczas funkcja wagowa ma postać:
(n + 1)L
®n+1= (2n + 1 ¡ x)L
®n¡ nL
®n¡1p(x) = x
®e
¡xWęzły xk są pierwiastkami wielomianu LN+1(x).
Wzór całkowania:
Kwadratura Gaussa-Hermite'a
A
k= ((N + 1)!)
2L
0N +1(x
k)L
N +2(x
k)
E(f ) = ((N + 1)!)
2(2N + 2)! f
(2N +2)(´)
´ 2 (0; 1)
Z
10
e
¡xf (x)dx ¼ S(f) = X
N k=0A
kf (x
k)
p(x) = e
¡x2(a; b) = ( ¡1; 1)
Ciąg wielomianów ortogonalnych stanowią wielomiany Hermite'a
Relacja rekurencyjna
H
n(x) = ( ¡1)
ne
x2d
ndx
ne
¡x2H
0(x) = 1 H
1(x) = 2x
H
2(x) = 4x
2¡ 2
H
3(x) = 8x
3¡ 12x
H
n+1= 2xH
n¡ 2nH
n¡124 Całkowanie funkcji wielu zmiennych
Przy całkowaniu funkcji wielu zmiennych pojawiają się problemy:
1) Konstrukcja wielomianów interpolacyjnych jest możliwa tylko dla odpowiednio położonych węzłów i regularnych obszarów całkowania
2) Czas obliczeń rośnie bardzo szybko wraz z liczbą zmiennych.
W praktyce liczba zmiennych nie przekracza 4.
Zakładamy, że obszar całkowania można opisać układem nierówności:
Szukamy wartości całki wielokrotnej:
½ R
MI(f ) = Z
: : : Z
| {z }
f (x
1; x
2; : : : ; x
M)dx
1: : : dx
Ma
1· x
1· b
1a
2(x
1) · x
2· b
2(x
1) : : : : : : : : : : : : : : : : : : : : : : : :
a
M(x
1; x
1; : : : ; x
M¡1) · x
M· b
M(x
1; x
2; : : : ; x
M¡1)
25 Wartość całki wielokrotnej oblicza się poprzez M-krotne zastosowanie
kwadratur jednowymiarowych.
Przykład dla dwóch wymiarów.
Po złożeniu obu kwadratur otrzymujemy:
I(f ) =
Z
b1a1
dx
1Z
b2(x1)a2(x1)
dx
2: : :
Z
bM(x1;x2;:::;xM¡1)aM(x1;x2;:::;xM¡1)
f (x
1; x
2; : : : ; x
M)dx
MI(f ) =
Z
b1a1
g(x
1)dx
1g(x
1) =
Z
b2a2
f (x
1; x
2)dx
2I
N1(g) =
N1
X
n=0
A
ng(x
1;n) g(x
1;n) ¼ I
N2;n(f
n) =
N
X
2;nº=0
B
º;nf (x
1;n; x
2;º)
I(f ) = ZZ
f (x
1; x
2)dx
1dx
2=
N1
X
n=0
N
X
2;nº=0
A
nB
º;nf (x
1;n; x
2;º) + R
N1(g) +
+
N
X
2;nA R (f )
26 gdzie: -reszta kwadratury
-reszta kwadratury
Uwagi:
1) Przedział całkowania po zmiennej x2 może się zmieniać wraz z wartością x1 2) Liczba węzłów kwadratur może być różna dla każdego węzła x1,n 3) Liczba użytych węzłów
Jeśli liczba w każdej kwadraturze byłaby jednakowa i równa (N+1) wówczas obliczenie wartości całki w M wymiarowej przestrzeni wiązałoby się z wykonaniem (N+1)M obliczeń.
Przykład.
Jeśli N=10 i M=10 wówczas (N+1)M >25.9 109
Przy dużej liczbie wymiarów (M>4) lepiej jest posługiwać się znacznie wydajniejszą metodą Monte Carlo.
N1
X
n=0