• Nie Znaleziono Wyników

Całkowanie numeryczne przy użyciu kwadratur

N/A
N/A
Protected

Academic year: 2021

Share "Całkowanie numeryczne przy użyciu kwadratur"

Copied!
24
0
0

Pełen tekst

(1)

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)

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

j

x

k

¡ x

j

C =

Z

b

a

f (x)dx

(3)

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

b

a

F (x)dx ¼ Z

b

a

'(x)dx = X

N k=0

A

k

F (x

k

)

A

k

= Z

b

a

©

k

(x)dx

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

¯ ¯

¯ ¯

¯ Z

b

a

F (x)dx ¡ X

N k=0

A

k

F (x

k

)

¯ ¯

¯ ¯

¯ =

= ¯¯ ¯

¯ ¯ Z

b

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

0k

= Z

b

a

p(x)©

k

(x)dx Z

b

a

F (x)dx =

Z

b a

p(x)f (x)dx ¼

¼

Z

b a

p(x)'(x)dx = X

N k=0

A

0k

f (x

k

)

(4)

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

b

a

p(x)f (x)dx

S(f ) = X

N k=0

A

k

f (x

k

); x 2 [a; b]

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

I(W ) = S(W )

B

N

= X

N k=0

jA

Nk

j 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=0

f (x

k

k

(x)

©

k

(x) = Y

j=0 j6=k

x ¡ x

j

x

k

¡ x

j

R

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

(5)

Ostatecznie otrzymujemy:

gdzie:

S(f ) = X

N k=0

A

k

f

k

A

k

= h ( ¡1)

N¡k

k!(N ¡ k)!

Z

N

0

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

r

f

(r)

(»); » 2 [a; b]

N

lim

!1

jA

k

j = 1 Z

b

a

f (x)dx =

Z

b a

'(x)dx =

=

X

N k=0

f

k

Z

b

a

©

k

(x)dx

= h X

N k=0

f

k

Z

N

0

'

k

(t)dt

= h X

N k=0

f

k

A

k

(6)

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

1

0

(t ¡ 1)dt = 1 2 h A

1

= h

Z

1 0

tdt = 1 2 h S(f ) = 1

2 h(f

0

+ f

1

) h = b ¡ a

E(f ) = 1 2!

Z

b

a

(x ¡ a)(x ¡ b)f

(2)

(»)dx

= ¡ 1

12 h

3

f

(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

5

f

(4)

(») E(f ) »

Z

b a

(x ¡ a) µ

x ¡ a + b 2

(x ¡ b)dx = 0

(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 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)

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=0

1

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])

»

k

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

3

12

m

X

¡1 k=0

f

(2)

k

)

= ¡ (b ¡ a)

3

12m

2

¢ 1

m

m

X

¡1 k=0

f

(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=0

f

(2)

k

) = f

(2)

(»); » 2 [a; b]

E(f ) = ¡ (b ¡ a)

3

12m

2

f

(2)

(»)

(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

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])

»

k

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

2 m

m=2

X

k=1

f

(4)

k

) = f

(4)

(»)

E(f ) = ¡ (b ¡ a)

5

180m

4

f

(4)

(»)

E(f ) = ¡ h

5

90

m=2

X

k=1

f

(4)

k

)

(10)

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=0

1

k! h

k

f

(k)

(x) f (x ¡ h) =

X

1 k=0

1

k! ( ¡1)

k

h

k

f

(k)

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

(1)

(x) + h

2

2 f

(2)

(x) + h

3

6 f

(3)

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

(1)

(x) + h

2

2 f

(2)

(x) ¡ h

3

6 f

(3)

(x) + : : :

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

(1)

(x) + 2

3! h

3

f

(3)

(x) + 2

5! h

5

f

(5)

(x) + : : :

f

(1)

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

2h ¡

· 1

3! h

2

f

(3)

(x) + 1

5! h

4

f

(5)

(x) + 1

7! h

6

f

(7)

(x) + O(h

8

)

¸

(11)

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

4

h

4

+ b

6

h

6

+ : : : L

h=2;1

= Á

µ h 2

+ a

2

h

2

4 + a

4

h

4

16 + a

6

h

6

64 + : : :

L

h=2;2

= Ã

µ h ¶

+ b

4

h

4

+ b

6

h

6

L

h;1

= Á (h) + a

2

h

2

+ a

4

h

4

+ a

6

h

6

+ : : :

L

h;3

= '(h) + c

6

h

6

+ c

8

h

8

+ : : : L

2

= (4

1

L

h=2;1

¡ L

h;1

)=(4

1

¡ 1)

= 4 3 Á

µ h 2

¡ 1

3 Á(h) ¡ a

4

h

4

4 ¡ 5a

6

h

6

16 ¡ : : :

L

3

= (4

2

L

h=2;2

¡ L

h;2

)=(4

2

¡ 1)

= 16 15 Ã

µ h 2

¡ 1

15 Ã(h) ¡ b

6

h

6

20 Ã(h) ¡ : : :

'(h) = 16 15 Ã

µ h 2

¡ 1

15 Ã(h)

(12)

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

k

4

k

¡ 1 D

n;k¡1

¡ 1

4

k

¡ 1 D

n¡1;k¡1

k = 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=k

A

jk

µ h 2

n

2j

D

0;0

D

1;0

D

1;1

D

2;0

D

2;1

D

2;2

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

D

M;0

D

M;1

D

M;2

: : : D

M;M

(13)

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

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

Ã

n

X

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

2n

X

n

i=1

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

h

n

= b ¡ a 2

n

S

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)

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

n

jR

k;k

¡ R

k¡1;k¡1

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: 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

b

a

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

5

90 f

(4)

(»)

» 2 [a; b]

Z

b

a

f (x)dx = X

n

i=1

(S

i

+ e

i

)

je

i

j · " x

i

¡ x

i¡1

b ¡ a R

0;0

= 1

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

R

n;0

= 1

2 R

n¡1;0

+ b ¡ a 2

n

2

X

n¡1 i=1

f µ

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

n

R

n;m

= R

n;m¡1

+ 4

m

R

n;m¡1

¡ R

n¡1;m¡1

4

m

¡ 1

S(a; b) = b ¡ a 3

½

f (a) + 4f

µ a + b 2

+ f (b)

¾

(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

n

i=1

e

i

¯ ¯

¯ ¯

¯ · X

n

i=1

je

i

j · "

b ¡ a X

n

i=1

(x

i

¡ x

i¡1

) = "

¯ ¯

¯ ¯

¯ Z

b

a

f (x)dx ¡ X

n

i=1

S

i

¯ ¯

¯ ¯

¯ · "

(16)

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=0

A

k

f (x

k

)

A

k

= Z

b

a

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

b

a

p(x)'

r

(x)'

s

(x)dx = 0

r 6= s

f'

n

(x) g = f'

0

(x); '

1

(x); : : : ; '

N

(x) g

(17)

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=0

f (x

j

)l

j

(x)

l

j

(x) = !

n

(x)

(x ¡ a

j

)!

n0

(a

j

)

y = d

j n

X

¡1

k=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

°

n

Z

b a

p(x) '

0

(x)'

n

(x) x ¡ d

j

dx

!

n

(x) = '

n

(x)

'

0

(x) = 1

®

k

= ¯

k+1

¯

k

1 = ¡ '

n+1

(d

j

)

®

n

°

n

Z

b a

p(x) '

n

(x) x ¡ d

j

dx

= ¡ '

n+1

(d

j

)'

0n

(d

j

)

®

n

°

n

Z

b a

p(x)l

j

(x)dx

= ¡ '

n+1

(d

j

)'

0n

(d

j

) A

°

k

= Z

b a

p(x)'

2k

(x)dx

(18)

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

2

2

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

3

f

(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

b

a

f (t)dt = b ¡ a 2

Z

1

¡1

g(x)dx

(19)

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

b

a

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

X

N k=0

A

k

f (t

k

)

t

k

= a + b

2 + b ¡ a 2 x

k

Poszukiwanie 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

(20)

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

¡x

L

n

(x) = ( ¡1)

n

e

x

d

n

dx

n

(x

n

e

¡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¡1

p(x) = x

®

e

¡x

(21)

Węzły xk są pierwiastkami wielomianu LN+1(x).

Wzór całkowania:

Kwadratura Gaussa-Hermite'a

A

k

= ((N + 1)!)

2

L

0N +1

(x

k

)L

N +2

(x

k

)

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

2

(2N + 2)! f

(2N +2)

(´)

´ 2 (0; 1)

Z

1

0

e

¡x

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

N k=0

A

k

f (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)

n

e

x2

d

n

dx

n

e

¡x2

H

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¡1

(22)

24 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

M

I(f ) = Z

: : : Z

| {z }

­

f (x

1

; x

2

; : : : ; x

M

)dx

1

: : : dx

M

a

1

· x

1

· b

1

a

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

)

(23)

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

b1

a1

dx

1

Z

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

M

I(f ) =

Z

b1

a1

g(x

1

)dx

1

g(x

1

) =

Z

b2

a2

f (x

1

; x

2

)dx

2

I

N1

(g) =

N1

X

n=0

A

n

g(x

1;n

) g(x

1;n

) ¼ I

N2;n

(f

n

) =

N

X

2;n

º=0

B

º;n

f (x

1;n

; x

2;º

)

I(f ) = ZZ

­

f (x

1

; x

2

)dx

1

dx

2

=

N1

X

n=0

N

X

2;n

º=0

A

n

B

º;n

f (x

1;n

; x

2;º

) + R

N1

(g) +

+

N

X

2;n

A R (f )

(24)

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

(N

2;n

+ 1) R

N1

(g) I

N1

(g)

R

N2;n

(f

n

) I

N2;n

(f

n

)

I

N2;n

(f

n

)

Cytaty

Powiązane dokumenty

Wartość pierwszej komórki w kolumnie C dajemy równą 0, a następnie liczymy w kolejnej komórce pole paska, według wzoru na pole trapezu.. Formułę tę kopiujemy

W poprzedniej metodzie, dla n = 2 otrzymywaliśmy dokładne przybliżenie całki dla wielomianów dru- giego stopnia, w ogólności można otrzymać dokładne przybliżenie całki

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

3.. W sprawozdaniu należy dodatkowo: a) przedyskutować dokładność oszacowania wartości całki ze względu na stopień wielomianu podcałkowego i liczbę użytych węzłów,

Funkcje gaussowskie stanowią w naszym przypadku funkcje wagowe, a właściwą funkcją podcałkową jest wyraz 1/r 12... Wyniki zapisać

[r]

Niestety, dla innych całek takiej kontroli najczęściej nie mamy – gdybyśmy nie znali wyniku analitycznego opisywanej przykładowej funkcji, we wniosku końcowym

Duży otwór – o rozmiarze większym od długości fali – traktuje się jako wiele stykających się z sobą małych otworków i mówi: każdy z tych fikcyjnych małych otworków