Całkowanie numeryczne przy użyciu kwadratur

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

d) Formuła sumacyjna Eulera-Maclaurina (Romberg, Burilsch) 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(x₀), f(x₁),...,f(x_N) definiujemy wielomian interpolacyjny Lagrange'a:

'(x) = L

(x) = X

©

(x)F (x

)

©

(x) = (x ¡ x

)(x ¡ x

) : : : (x ¡ x

)(x ¡ x

) : : : (x ¡ x

)

(x

¡ x

)(x

¡ x

) : : : (x

¡ x

)(x

¡ x

) : : : (x

¡ x

) = Y

j=0 j6=k

x ¡ x

x

¡ x

C =

Z

a

f (x)dx

3 Podstawiamy wielomian interpolacyjny w

miejsce funkcji podcałkowej:

Powyższe wzory definiują tzw. kwadraturę. A_k 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 podcałkowej wielomianem (lub inną funkcją).

Z

a

F (x)dx ¼ Z

a

'(x)dx = X

A

F (x

)

A

= Z

a

©

(x)dx

jf(x) ¡ '(x)j < "; x 2 [a; b]

¯ ¯

¯ ¯

¯ Z

a

F (x)dx ¡ X

A

F (x

)

¯ ¯

¯ ¯

¯ =

= ¯¯ ¯

¯ ¯ Z

a

(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

= Z

a

p(x)©

(x)dx Z

a

F (x)dx =

Z

p(x)f (x)dx ¼

¼

Z

p(x)'(x)dx = X

A

f (x

)

4 Chcemy wyznaczyć wartość całki:

Stosując wzór

Powyższy wzór nosi nazwę kwadratury, a punkty x₁,x₂...,x_N 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

a

p(x)f (x)dx

S(f ) = X

A

f (x

); x 2 [a; b]

E(f ) = I(f ) ¡ S(f)

I(W ) = S(W )

B

= X

jA

j 6 M; N = 1; 2; : : :

Kwadratury Newtona-Cotesa

Rozważamy przypadek z węzłami równoodległymi x_i=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

) = L

(x

); i = 0; 1; 2; : : : ; N L

(x) =

X

f (x

)©

(x)

©

(x) = Y

j=0 j6=k

x ¡ x

x

¡ x

R

(x) = f (x) ¡ L

(x)

= 1

(N + 1)! !

(x)f

(»)

» 2 (a; b)

x = a + ht ©

(t) = Y

j=0 j6=k

t ¡ j

k ¡ j

5 Ostatecznie otrzymujemy:

gdzie:

S(f ) = X

A

f

A

= h ( ¡1)

k!(N ¡ k)!

Z

0

t(t ¡ 1) : : : (t ¡ N) (t ¡ k)

f

= f (a + kh)

h = b ¡ a N

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 C_r 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 A_k

4) Współczynniki A_k zależą od N. W szczególności (wzór na A_k) 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

f

(»); » 2 [a; b]

N

lim

jA

j = 1 Z

a

f (x)dx =

Z

'(x)dx =

=

X

f

Z

a

©

(x)dx

= h X

f

Z

0

'

(t)dt

= h X

f

A

6 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 dokładnie przybliża wielomian N+1=1+1=2 stopnia.

Zatem błąd wyznaczenia przybliżonej wartości całki wynosi

A

= ¡h Z

0

(t ¡ 1)dt = 1 2 h A

= h

Z

tdt = 1 2 h S(f ) = 1

2 h(f

+ f

) h = b ¡ a

E(f ) = 1 2!

Z

a

(x ¡ a)(x ¡ b)f

(»)dx

= ¡ 1

12 h

f

(»); » 2 [a; b]

A

= 4

3 h A

= 1 3 h A

= 1

3 h

h = b ¡ a 2

S(f ) = 1

3 h(f

+ 4f

+ f

)

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

» 2 [a; b]

E(f ) = f

(»

) 4!

Z

(x ¡ a) µ

x ¡ a + b 2

¶

(x ¡ b)dx

= ¡ 1

90 h

f

(»)

7

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 h³ (1/12) f⁽²⁾(» ) trapezów

2 (1/3)h 1 4 1 h⁵ (1/90) f⁽⁴⁾(») parabol

3 (3/8)h 1 3 3 1 h⁵ (3/80) f⁽⁴⁾(») 3/8

4 (4/90)h 7 32 12 32 7 h⁷ (8/945) f⁽⁶⁾(») Milne'a 5 (5/288)h 19 75 50 50 75 19 h⁷ (275/12096) f⁽⁶⁾(») --- 6 (6/840)h 41 216 27 272 27 216 41 h⁹ (9/1400) f⁽⁸⁾(») 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

2 h(f

+ f

)

= h µ 1

2 f

+ f

+ : : : + f

+ 1 2 f

¶

f

= f (a + k ¢ h) f 2 C

([a; b])

»

2 (a + kh; a + (k + 1)h) E(f ) = ¡ h

12

m

X

f

(»

)

= ¡ (b ¡ a)

12m

¢ 1

m

m

X

f

(»

)

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

f

(»

) = f

(»); » 2 [a; b]

E(f ) = ¡ (b ¡ a)

12m

f

(»)

9 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

+ 4f

+ f

)

= h

3 [f

+ f

+ 2(f

+ f

+ : : : + f

) + 4(f

+ f

+ : : : + f

)]

f 2 C

([a; b])

»

2 (a + 2(k ¡ 1)h; a + 2kh)

2 m

m=2

X

k=1

f

(»

) = f

(»)

E(f ) = ¡ (b ¡ a)

180m

f

(»)

E(f ) = ¡ h

90

m=2

X

k=1

f

(»

)

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! h

f

(x) f (x ¡ h) =

X

1

k! ( ¡1)

h

f

(x) f (x + h) = f (x) + hf

(x) + h

2 f

(x) + h

6 f

(x) + : : : f (x ¡ h) = f(x) ¡ hf

(x) + h

2 f

(x) ¡ h

6 f

(x) + : : :

f (x + h) ¡ f(x ¡ h) = 2hf

(x) + 2

3! h

f

(x) + 2

5! h

f

(x) + : : :

f

(x) = f (x + h) ¡ f(x ¡ h)

2h ¡

· 1

3! h

f

(x) + 1

5! h

f

(x) + 1

7! h

f

(x) + O(h

)

¸

11 Powyższa formuła w postaci ogólnej

co można interpretować jako przybliżenie f⁽¹⁾(x).

Za h podstawimy h/2

i obliczmy różnicę

Zatem L₁ przybliża f⁽¹⁾(x) z dokładnością O(h⁴) (wyrazów rzędu h⁴).

Dokonujemy podstawienia

w L₁

Ã(h) = 4 3 Á

µ h 2

¶

¡ 1

3 Á(h)

'(h) = 15 16 Ã

µ h 2

¶

¡ 1

15 Ã(h)

Podstawiając do L₂

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(h^2M).

(o ile h<<1).

L

= Ã µ h

2

¶

+ b

h

16 + b

h

64 L

= Á(h) + a

h

+ a

h

+ a

h

+ : : :

L

= Á µ h

2

¶

+ a

h

4 + a

h

16 + a

h

64 + : : : L

= (L

¡ 4L

)=3

= 4 3 Á

µ h 2

¶

¡ 1

3 Á(h) ¡ a

h

4 ¡ 5a

h

16 ¡ : : :

L

= Ã(h) + b

h

+ b

h

+ : : :

L

= (L

¡ 4L

)=3

= 15 16 Ã

µ h 2

¶

¡ 1

15 Ã(h) ¡ b

h

20 Ã(h) ¡ : : :

L

= '(h) + c

h

+ c

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

= L + O(h

) D

= L + O(h

) D

= L + O(h

) D

= L + O(h

)

: : : : : : : : : : : :

D

= L + O(h

); h ! 0 D

= 4

4

¡ 1 D

¡ 1

4

¡ 1 D

k = 1; 2; : : : ; M

n = k; k + 1; : : : ; M D

= Á

µ h 2

¶

; 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 D_m,m.

Kolejne kroki algorytmu można zapisać w postaci tablicy

D

= L + X

A

µ h 2

¶

D

D

D

D

D

D

.. . .. . .. . . ..

D

D

D

: : : D

13 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 S

= h

X

f (a + ih) ¡ f (a) + f (b) 2

= b ¡ a n

X

f µ

a + i b ¡ a n

¶

¡ f (a) + f (b) 2

x 2 [0; 1]

S

= 1

2 f (0) + 1

2 f (1) S

= 1

4 f (0) + 1 2

· f

µ 1 2

¶¸

+ 1

4 f (1) S

= 1

8 f (0) + 1 4

· f

µ 1 4

¶ + f

µ 1 2

¶ + f

µ 3 4

¶¸

+ 1

8 f (1) S

= 1

16 f (0) + 1 8

· f

µ 1 8

¶ + f

µ 1 4

¶ + f

µ 3 8

¶ + f

µ 1 2

¶

+f µ 5

8

¶ + f

µ 3 4

¶ + f

µ 7 8

¶¸

+ 1

16 f (1)

Łatwo zauważyć, że do obliczenia T_2n można wykorzystać już obliczone T_n

co ogólnie dla przedziału całkowania [a,b] można zapisać jako

S

= 1

2 S

+ 1 2

· f

µ 1 2

¶¸

S

= 1

2 S

+ 1 4

· f

µ 1 4

¶ + f

µ 3 4

¶¸

S

= 1

2 S

+ 1 8

· f

µ 1 8

¶ + f

µ 3 8

¶ + f

µ 5 8

¶ + f

µ 7 8

¶¸

S

= 1

2 S

+ h X

i=1

f (a + (2i ¡ 1)h)

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

= b ¡ a 2

jR

¡ R

j · "; " 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: e_i jest lokalnym błędem przybliżenia wartości całki w i-tym podprzedziale [x_i-1,x_i].

Załóżmy że jego wartość możemy oszacować zgodnie z poniższym wzorem

Z

a

f (x)dx = S(a; b) ¡ (b ¡ a)

90 f

(»)

S(a; b) = b ¡ a 6

·

f (a) + 4f

µ a + b 2

¶

+ f (b)

¸

» 2 [a; b]

Z

a

f (x)dx = X

i=1

(S

+ e

)

je

j · " x

¡ x

b ¡ a R

= 1

2 (b ¡ a) [f(a) + f (b)]

R

= 1

2 R

+ b ¡ a 2

2

X

f µ

a + (2i ¡ 1) b ¡ a 2

¶

R

= R

+ R

¡ R

4

¡ 1

15 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

i=1

e

¯ ¯

¯ ¯

¯ · X

i=1

je

j · "

b ¡ a X

i=1

(x

¡ x

) = "

¯ ¯

¯ ¯

¯ Z

a

f (x)dx ¡ X

i=1

S

¯ ¯

¯ ¯

¯ · "

16 Formuła sumacyjna Eulera-Maclaurina

Jest metodą ekstrapolacyjną (można z niej uzyskać np. metodę Romberga)

Dla funkcji

gdzie: współczynniki B_2l oraz B_2m+2 są liczbami Bernoulliego.

Dowód – Stoer, Burilsch „Introduction to ...”

Liczby te wyznacza się korzystając z wielomianów Bernouliego – to wartości wielomianu dla x=0

B

(x) = (k + 1)B

(x) B

(x) = 1; B

(x) = x ¡ 1

2 ; B

(x) = x

¡ x + 1 6 B

(x) = x

¡ 3

2 x

+ 1

2 x; B

(x) = x

¡ 2x

+ x

¡ 1 30 B

= B

(0) ) B

= 1

6 ; B

= ¡ 1

30 ; B

= 1

42 ; B

= ¡ 1

30 ; : : : f (x) 2 C

; x 2 [0; 1]

Z

0

f (x)dx = f (0)

2 + f (1) 2 +

X

B

(2l)! (f

(0) ¡ f

(1))

¡ B

(2m + 2)! g

(»); 0 < » < 1

17 Zastępując funkcję g(x) jej wartościami określonymi na siatce

równoodległych węzłów

W powyższym wzorze rozpoznajemy wzór złożony trapezów. Po przegrupowaniu wyrazów dostajemy związek formuły Eulera- Maclaurina z kwadraturą:

h = b ¡ a N x

= a + ih

i = 0; 1; : : : ; N Z

0

f (x)dx =

X

Z

xi¡1

f (x)dx = f (x

)

2 + f (x

) + : : : + f (x

) + f (x

) 2 +

X

h

B

(2l)! (f

(x

) ¡ f

(x

))

¡ h

B

(2m + 2)! (b ¡ a)f

(»); x

< » < x

S

= f (x

)

2 + f (x

+ h) + : : : + f (x

+ h(N ¡ 1)) + f (x

+ hN ) 2

=

Z

x0

f (x)dx + X

l=1

h

B

(2l)! (f

(x

) ¡ f

(x

))

¡ h

B

(b ¡ a)f

(»); x

< » < x

18 Związek formuły E-L z metodą Romberga

Jeśli pochodna f^(2m+2) jest ciągła w [a,b] oraz istnieje taka liczba L, że

(czyli pochodna nie posiada osobliwości) to wówczas istnieje także taka liczba M_m+1, że zachodzi:

Sposoby podziału przedziału całkowania:

a) Romberga

b) Burilscha - liczba wykonywanych operacji nie rośnie tak szybko jak w metodzie Romberga

S(h) = ¿

+ ¿

h

+ ¿

h

+ : : : + ¿

h

+ ®

(h)h

¿

= B

(2k)! (f

(b) ¡ f

(a)); k = 1; 2; : : : ; m

®

(h) = B

(2m + 2)! (b ¡ a)f

(»(h)); a < » < b

j®

(h) j · M

¿

= Z

a

f (x)dx

h

= b ¡ a

2

; i = 0; 1; 2; : : :

h

= b ¡ a; h

= h

2 ; h

= h

3 ; : : : ; h

h

2 ; i = 3; 4; : : :

S

= S

+ S

¡ S

h

i

jf

(x) j · L; x 2 [a; b] ¡ 1

S

¼ ¿

19 Wzór na S(h) jest rozwinięciem asymptotycznym w h jeśli współczynniki

τ_k k ≤ m nie zależą od h. Jeśli m → ∞ wówczas prawa strona w formule E-M staje się szeregiem nieskończonym

Uwaga: dla h>0 nieskończony szereg jest rozbieżny, ale dla m skończonego (małe) wzór pracuje dobrze, tj. zwiększajac m minimalizujemy wszystkie wyrazy

poza τ₀ (czyli poszukiwaną wartością całki).

Przykład. Obliczyć wartość całki

¿

+ ¿

h

+ ¿

h

+ : : :

I = Z

0

t

dt h

= 1; h

= 1

2 ; 1 4

S

= 0:500000 S

= 0:265625 S

= 0:192383

S

= 0:187500 S

= 0:167969

S

= 0:166667(= 1=6)

20 Formuła Eulera-Maclaurina (met. Romberga) bazuje na

wielomianowej interpolacji funkcji podcałkowej. Okazuje się że znacznie lepszym przybliżeniem są ilorazy wielomianowe (funkcje wymierne).

Definicja funkcji S_ik następująca

z warunkami

Także wzór rekurencyjny ulega modyfikacji

S

(h) = p

+ p

h

+ : : : + p

h

q

+ q

h

+ : : : + q

h

¹ + º = k; ¹ = º _ ¹ = º ¡ 1

S

= S

+ S

¡ S

³

hi

´

³

1 ¡

´

¡ 1 1 · k · i · m

T

= 0; i = 0; 1; : : : ; m ¡ 1

21 Kwadratury Gaussa

Nadal rozpatrujemy kwadratury typu:

ale nieco zmienimy metodologię postępowania.

Ustalamy funkcję wagową p(x) oraz liczbę węzłów N. Szukamy:

a) położenia węzłów b) współczynników A_k

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

A

f (x

)

A

= Z

a

p(x)©

(x)dx

(P

(x)) = P

(x); P

(x); : : : ; P

(x)

(P

; P

) = Z

a

p(x)P

(x)P

(x)dx = 0 r 6= s

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 tytlko wtedy, gdy węzły x_k są pierwiastkami wielomianu P_N+1(x).

Tw. 3. Wszystkie współczynniki A_k w kwadraturach Gaussa są dodatnie.

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.

Kwadratura dla przedziału skończonego (Gaussa-Legendre'a).

Dla tego typu kwadratury przyjmujemy:

W tym przedziale ciąg wielomianów ortogonalnych tworzą wielomiany Legendre'a

p(x) = 1

[a; b] = [ ¡1; 1]

P

(x) = 1 2

¢ n!

d

dx

(x

¡ 1)

22 Współczynniki A_k:

Błąd kwadratury:

Węzły x_k stanowią pierwiastki wielomianu P_N+1(x).

(jak je znaleźć? => metody poszukiwania zer wielomianów)

Dla kwadratur niskiego rzędu węzły i współczynniki A_k są stablicowane. Aby zastosować wzory z przedziału [-1,1] w przedziale [a,b] należy dokonać transformacji liniowej zmiennej niezależnej:

P

(x) = 1 P

(x) = x

P

(x) = 3x

¡ 1 2

P

(x) = 5x

¡ 3x

2

A

= ¡ 2

(N + 2)P

(x

)P

(x

)

E(f ) = 2

((N + 1)!)

(2N + 3)((2N + 2)!)

f

(»)

¡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

a

f (t)dt = b ¡ a 2

Z

¡1

g(x)dx

23 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

a

f (t)dt ¼ S(f) = b ¡ a 2

X

A

f (t

)

t

= a + b

2 + b ¡ a 2 x

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

L

(x) = ( ¡1)

e

d

dx

(x

e

) L

(x) = 1

L

(x) = 1 ¡ x

L

(x) = x

¡ 4x + 2 2

L

(x) = ¡x

+ 9x

¡ 18x + 6

6

24 Węzły x_k są pierwiastkami wielomianu L_N+1(x).

Wzór całkowania:

Kwadratura Gaussa-Hermite'a

A

= ((N + 1)!)

L

(x

)L

(x

)

E(f ) = ((N + 1)!)

(2N + 2)! f

(´)

´ 2 (0; 1)

Z

0

e

f (x)dx ¼ S(f) = X

A

f (x

)

p(x) = e

(a; b) = ( ¡1; 1)

Ciąg wielomianów ortogonalnych stanowią wielomiany Hermite'a

H

(x) = ( ¡1)

e

d

dx

e

H

(x) = 1 H

(x) = 2x

H

(x) = 4x

¡ 2

H

(x) = 8x

¡ 12x

25 x_k są zerami wielomianu H_N+1

Wzór przybliżonego całkowania:

A

= 2

(N + 1)!

H

(x

)H

(x

)

E(f ) = (N + 1)! p

¼

2

(2N + 2)! f

(´)

´ 2 (¡1; 1)

Z

¡1

e

f (x)dx ¼ S(f) = X

A

f (x

)

Uwagi końcowe:

1) Kwadratury Gaussa są dokładniejsze od kwadratur Newtona-Cotesa przy uwzględnieniu tej samej liczby węzłów

2) Kwadratury Gaussa mają rząd r=2N+2 dla (N+1) węzłów, podczas gdy kwadratury NC osiągają ten rząd dla (2N+1) węzłów

3) Po ustaleniu rzędu kwadratury stosuje się wzory złożone dla coraz mniejszego kroku całkowania do momentu braku zmian w kolejnym przybliżeniu 4) Całkowanie stablicowanej funkcji podcałkowej

lepiej wykonać przy użyciu kwadratur Newtona- Cotes'a (użycie kwadratur Gaussa może wymagać dodatkowej interpolacji)

26 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

I(f ) = Z

: : : Z

| {z }

­

f (x

; x

; : : : ; x

)dx

: : : dx

a

· x

· b

a

(x

) · x

· b

(x

) : : : : : : : : : : : : : : : : : : : : : : : :

a

(x

; x

; : : : ; x

) · x

· b

(x

; x

; : : : ; x

)

27 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

a1

dx

Z

a2(x1)

dx

: : :

Z

aM(x1;x2;:::;xM¡1)

f (x

; x

; : : : ; x

)dx

I(f ) =

Z

a1

g(x

)dx

g(x

) =

Z

a2

f (x

; x

)dx

I

(g) =

N1

X

n=0

A

g(x

) g(x

) ¼ I

(f

) =

N

X

º=0

B

f (x

; x

)

I(f ) = Z Z

­

f (x

; x

)dx

dx

=

N1

X

n=0

N

X

º=0

B

f (x

; x

) + R

(g) +

+

N

X

A R (f )

28 gdzie: -reszta kwadratury

-reszta kwadratury

Uwagi:

1) Przedział całkowania po zmiennej x₂ może się zmieniać wraz z wartością x₁ 2) Liczba węzłów kwadratur może być różna dla każdego węzła x_1,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=1 i M=10 wówczas (N+1)^M=(1+1)¹⁰=1048576

Przy dużej liczbie wymiarów (M>4) lepiej jest posługiwać się znacznie wydajniejszą metodą Monte Carlo.

N1

X

n=0

(N

+ 1) R

(g) I

(g)

R

(f

) I

(f

)

I

(f

)