Całkowanie numeryczne

Całkowanie numeryczne

Marcin Orchel

1 Wstęp

Całkowanie numeryczne – przybliżone wyliczanie wartości całek oznaczonych.

I (f ) = Z b

a

f (x) dx (1)

Do wszystkich metod. Podział przedziału całkowania na podprzedziały. Przedział cał- kowania: [a, b]. x ₁ , x 2 , . . . , x n – punkty należące do przedziału całkowania, takie, że x _i < x _i+1 , dla i = 1, 2, . . . , n − 1, x ₁ = a, x _n = b. Podprzedziały całkowania: (x ₁ , x ₂ ) , (x ₂ , x ₃ ) , . . . , (x _n−1 , x _n ). Zdefiniujmy h _i := x _i+1 − x _i dla i = 1, 2, . . . , n − 1. Oznaczenie:

y i = f (x _i ).

1.1 Podejście geometryczne 1.1.1 Metoda prostokątów

W każdym podprzedziale wybierany jest punkt środkowy M _i o odciętej m _i = x _i+1 + x _i

2

dla i = 1, 2, . . . , n−1. Dla każdego podprzedziału pole pod wykresem funkcji przybliżane jest prostokątem o bokach długości x _i+1 − x _i oraz f (m _i ). Prostokąt i-ty to prostokąt z podprzedziału (x _i , x _i+1 ). Pole prostokąta i-tego wynosi

P i = h _i f (m i ) Szukana wartość całki oznaczonej wynosi:

I (f ) =

n−1

X

i=1

P _i =

n−1

X

i=1

h _i f (m _i )

Dla jakiego przebiegu funkcji metoda ta daje dobre przybliżenie, a dla jakiego nie?

Oszacowanie górnej granicy błędu dla i-tego przedziału:

R < m _i f (m _i )

Jeśli znana jest maksymalna pochodna w i-tym przedziale:

p _i = max f ⁰ (x)
to oszacowanie można poprawić:

R ≤ 1 2 1 2 h i

1 2 h i p i

R ≤ 1 8 h ² _i p _i

Z innej perspektywy można powiedzieć, że metoda prostokątów polega na interpolacji funkcji f funkcją stałą:

y i (x) = f (m _i ) . 1.1.2 Metoda trapezów

W każdym przedziale konstruowany jest trapez przechodzący przez punkty:

(x _i , 0) , (x _i+1 , 0) , (x _i , y _i ) , (x _i+1 , y _i+1 ) . (2) Pole tego trapezu wynosi:

P _i = (y _i + y _i+1 ) h _i 2 Szukana wartość całki wynosi zatem:

I (f ) =

n−1

X

i=1

P _i =

n−1

X

i=1

(y _i + y _i+1 ) h _i 2

Dla jakiego przebiegu funkcji metoda ta daje gorsze/lepsze przybliżenie od metody prostokątów?

Można również powiedzieć, że metoda trapezów polega na interpolacji funkcji f funk- cją liniową przechodzącą przez punkty (x _i , y _i ) , (x _i+1 , y _i+1 ).

1.1.3 Metoda parabol (Simpsona)

Metoda parabol polega na interpolacji funkcji f funkcją kwadratową w każdym prze- dziale. Dla każdego przedziału (x _i , x i+1 ) konstruujemy funkcję kwadratową przechodzącą przez 3 punkty (x _i , f (x i )), (s _i , f (s i )), (x _i+1 , f (x i+1 )), gdzie s _i jest wartością z przedziału (x _i , x _i+1 ). Funkcja kwadratowa będzie miała postać:

W (x) = a _i x ² + b _i x + c _i

dla i = 1, . . . , n − 1. Współczynniki a _i , b _i oraz c _i można znaleźć rozwiązując następujący układ równań:

y _i = a _i x ² _i + b _i x _i + c _i

f (s i ) = a _i s ² _i + b _i s i + c _i y i+1 = a _i x ² _i+1 + b _i x i+1 + c _i .

Po wyznaczeniu współczynników a _i , b _i oraz c _i należy obliczyć całkę:

I (f ) =

x

Z

x

a i x ² + b _i x + c i dx

=

"

a i

x ³

3 + b _i x ² 2 + c _i x

# x

x

= a _i 3



x ³ _i+1 − x ³ _i ^ + b _i 2



x ² _i+1 − x ² _i ^ + c _i (x _i+1 − x _i )

= (x _i+1 − x _i ) 6

 2a _i ^ x ² _i+1 + x _i x i+1 + x ² _i ^ + 3b _i (x _i+1 + x _i ) + c _i ^

= (x _i+1 − x _i ) 6

 a i x ² _i+1 + b _i x i+1 + c _i + a _i x ² _i + b _i x i + c _i + a _i x ² _i+1 + 2a _i (x _i x i+1 ) + a _i x ² _i + 2b _i (x _i+1 + x _i ) + 4c _i ^

= (x _i+1 − x _i )

6 a i x ² _i+1 + b _i x i+1 + c _i + a _i x ² _i + b _i x i + c _i + 4 a i

 x i + x _i+1 2

 2

+ b _i x i + x _i+1 2 + c _i

!!

Możemy uprościć powyższą postać podstawiając wartość odciętej m _i punktu środkowego (x _i+1 − x _i )

6 (y _i+1 + y _i + 4f (m _i ))

Znamy jednak wartość dla odciętej s _i , więc zakładamy, że punkt wewnętrzny jest punk- tem środkowym s _i = m _i .

1.2 Podejście analityczne

Interesuje nas zastąpienie funkcji f (x) funkcją g (x) tak aby poniższe równanie było w miarę możliwości spełnione

Z b

a

f (x) dx ≈ Z b

a

g (x) dx . (3)

Interesuje nas szczególny przypadek gdy

b

Z

a

f (x) dx ≈

n

X

i=0

a _i f (x _i ) (4)

gdzie współczynniki a _i nie zależą od f . Wzór ten nazywany jest kwadraturą, punkty x _i

nazywane są węzłami kwadratury, a liczby a _i współczynnikami kwadratury.

1.2.1 Metoda oparta na interpolacji Lagrange’a

Pomysł polega na zastąpieniu funkcji f (x) wzorem interpolacyjnym Lagrange’a, czyli g (x) =

n

X

i=0

y _i L _i (x) (5)

gdzie L _i są wielomianami Lagrange’a. Po podstawieniu otrzymujemy przybliżenie

n

X

i=0

y _i

b

Z

a

L _i (x) dx . (6)

Widzimy, że ta postać jest szczególnym przypadkiem prawej strony (4) gdzie a _i =

Z b a

L _i (x) dx . (7)

Kwadratury Newtona-Cotesa Zobaczmy, jak wygląda powyższa metoda dla węzłów równoodległych. Możemy podstawić do wielomianu interpolacyjnego Lagrange’a

x = a + hs , (8)

gdzie a i h to stałe, a s to nowa zmienna i otrzymujemy L i (x) =

n

Y

k=0 k6=i

x − x _k x _i − x _k =

n

Y

k=0 k6=i

a + hs − a − kh a + ih − a − kh =

n

Y

k=0 k6=i

s − k

i − k (9)

Ponieważ jak widzimy powyżej L _i są niezależne od a i od h, więc możemy zapisać w szczególności

L _i (x) = ϕ _i (s) =

n

Y

k=0 k6=i

s − k

i − k (10)

Stosujemy całkowanie przez podstawienie (8) do (7) zmieniając granice całkowania, po zauważeniu, że dla s = 0 otrzymujemy a, a dla s = n otrzymujemy b

a i = h Z n

0

ϕ i (s) ds . (11)

możemy wprowadzić dodatkowo symbol α i =

Z n 0

ϕ i (s) ds . (12)

Współczynniki α _i nie zależą ani od granic całkowania, ani od funkcji f , tylko od n.

Możemy wyprowadzić wzory na współczynniki dla n = 1 α ₀ =

Z 1 0

ϕ ₀ (s) ds = Z 1

0

s − 1

0 − 1 dx = − 1

2 s ² + s| ¹ ₀ = 1

2 (13)

α 1 = Z 1

0

ϕ 1 (s) ds = s ² 2 | ¹ ₀ = 1

2 (14)

A więc przybliżenie całki jest równe S (f ) = 1

2 h (y ₀ + y ₁ ) . (15)

A więc otrzymujemy wzór trapezów. Dla n = 2 a ₀ =

Z 2 0

ϕ ₀ (s) ds = . . . = 1

3 (16)

a ₁ = Z 2

0

ϕ ₁ (s) ds = . . . = 4

3 (17)

a ₂ = Z 2

0

ϕ ₂ (s) ds = . . . = 1

3 (18)

Zauważmy, że suma wszystkich współczynników jest równa n. Przybliżenie całki jest równe

S (f ) = 1

3 h (y ₀ + 4y ₁ + y ₂ ) . (19)

Jest to wzór parabol (Simpsona).

Dla n = 3 otrzymujemy wzór 3/8 Newtona (wzór trzech ósmych).

S (f ) = 3

8 h (y 0 + 3y ₁ + 3y ₂ + y ₃ ) . (20) Przykład: Obliczyć wartość całki dokładną i przybliżoną metodami trapezów, Simp- sona i 3/8 Newtona

Z 1 0

dx

1 + x (21)

Metoda trapezów

Z 1 0

dx 1 + x = 1

2 1

 1 + 1

2



= 0.75 . (22)

Metoda Simpsona

Z 1 0

dx

1 + x = 0.6944 . (23)

Metoda 3/8 (trzech ósmych) Newtona Z 1

0

dx

1 + x = 0.6937 . (24)

Wartość dokładna

Z 1 0

dx

1 + x = ln 2 ≈ 0, 6931 . (25)

Dla funkcji f , która jest na [a, b] dowolnie wiele razy różniczkowalna, błąd kwadratury wynosi

Z b a

g n (x) dx − Z b

a

f (x) dx = h ^p+1 Kf ^(p) (ξ) (26) gdzie ξ ∈ (a, b), p i K zależą od n.

Dla n = 1, dla wzoru trapezów otrzymujemy h ³ 1

12 f ⁽²⁾ (ξ) (27)

Dla n = 2 dla wzoru Simpsona otrzymujemy h ⁵ 1

90 f ⁽⁴⁾ (ξ) (28)

Dla n = 3 dla wzoru trzech ósmych otrzymujemy h ⁵ 3

80 f ⁽⁴⁾ (ξ) (29)

Możemy stosować wzory Newtona-Cotesa dla podprzedziałów, a następnie zsumować wartości całek z podprzedziałów. Dla n = 1 dla całego przedziału

I (h) :=

n−1

X

i=0

h

2 (f (x _i ) + f (x _i+1 )) (30)

= h

 f (a)

2 + f (a + h) + f (a + 2h) + . . . + f (b − h) + f (b) 2



(31) Jest to złożony wzór trapezów. Błąd kwadratury dla każdego podprzedziału wynosi

I i − Z x

x

f (x) dx = h ³

12 f ⁽²⁾ (ξ _i ) (32)

gdzie ξ _i ∈ (x _i , x i+1 ). Dla ograniczenia drugiej pochodnej ograniczona 

 f ⁽²⁾ (x)  

 ≤ m (33)

dla x ∈ [a, b], błąd wynosi 

I (h) − Z b

a

f (x) dx     

=     

n−1

X

i=0

h ³

12 f ⁽²⁾ (ξ _i )     

≤ h ³

12 nm = b − a

12 h ² m (34)

1.2.2 Metoda Gaussa

Po angielsku Gaussian quadrature, http://en.wikipedia.org/wiki/Gaussian_quadrature.

Dla węzłów nierównoodległych, możemy otrzymać lepsze przybliżenie. 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 dla wielomianów n-tego stopnia. W metodzie Gaussa, możemy otrzymać dokładne przybliżenie całki dla wielomianów stopnia 2n. A więc dla n = 2 otrzymamy dokładne przybliżenie całki dla wielomianów 4 stopnia.

Mamy dany przedział skończony [a, b], a < b oraz dla x ∈ [a, b] dodatnia ciągła funkcja wagowa ω(x), n ≥ 1 to liczba naturalna. Zadanie polega na aproksymacji całki

I (f ) = Z b

a

w (x) f (x) (35)

wartością

I (f ) = ˜

n

X

i=1

w _i f (x _i ) (36)

a więc musimy wybrać odpowiednio 2n liczb: w _i , x _i , i = 1, . . . , n. Likwidujemy założenie aby węzły musiały być równoodległe. Zadanie polega na takim dobraniu w _i oraz x _i aby błąd aproksymacji

I (f ) − I (f ) ˜ (37)

znikał dla wielomianów możliwie dużego stopnia. Czyli gdy podstawimy za f wielomian.

Okazuje się, że dla każdego n = 1, 2, . . . istnieją jednoznacznie określone liczby w _i , x _i gdzie i = 1, . . . , n takie, że ˜ I(f ) = I(f ) dla wszystkich wielomianów p stopnia ≤ 2n − 1, oraz w _i > 0 i x i ∈ (a, b) dla i = 1, . . . , n.

Oznaczmy

Π ¯ _j := ⁿ p|p (x) = x ^j + a ₁ x ^j−1 + . . . + a _j ^o (38) zbiór wszystkich wielomianów unormowanych, oraz

Π _j := {p|p wielomian stopnia ≤ j} (39) Twierdzenie 1. Istnieją jednoznacznie określone wielomiany p j ∈ ˜ Π _j , j = 0, 1, . . . takie, że

(p _i , p k ) = 0 (40)

dla i 6= k. Wielomiany te spełniają związek rekurencyjny

p 0 (x) = 1 (41)

p _i+1 (x) = (x − δ _i+1 ) p _i (x) − γ _i+1 ² p _i−1 (x) (42) dla i ≥ 0, p −1 (x) ≡ 0 oraz

δ i+1 := (xp _i , p i ) / (p _i , p i ) (43) dla i ≥ 0, oraz

γ ² _i+1 :=

( 0 dla i ≥ 0

(p _i , p _i ) / (p _i−1 , p _i−1 ) dla i ≥ 1 (44)

Wielomiany p _i ze względu na pierwszy warunek są wielomianami ortogonalnymi z funkcją wagową ω(x). Dla przypomnienia

(p _i , p _k ) = Z b

a

w (x) p _i (x) p _k (x) dx (45)

Przykład. Wyliczenie wielomianu p ₁ (x) dla w(x) ≡ 1, wtedy i = 0 p ₁ (x) =



x − (xp ₀ (x) , p ₀ (x)) (p ₀ (x) , p ₀ (x))



p ₀ (x) − 0 =



x − (x, 1) (1, 1)



1 = x −

1 2 − ¹ ₂

2 = x (46) Dla wszystkich wielomianów p ∈ Π _n−1 zachodzi równość (p, p _n ) = 0.

Twierdzenie 2. Zera x i dla i = 1, . . . , n wielomianu p _n są rzeczywiste, pojedyncze i leżą w przedziale (a, b).

Twierdzenie 3. Dla dowolnych liczb t _i dla i = 1, . . . , n takich, że t ₁ < t ₂ < . . . < t _n macierz kwadratowa n × n

A :=







p ₀ (t ₁ ) . . . p ₀ (t _n ) . . . . . . . . . p n (t ₁ ) p n−1 (t _n )





 (47)

jest nieosobliwa.

Twierdzenie 4. Niech liczby w ₁ , . . . , w _n i zera x ₁ , . . . , x _n wielomianu p _n (x) będą roz- wiązaniem układu równań liniowych

n

X

i=1

p _k (x _i ) w _i =

( (p ₀ , p ₀ ) jesli k = 0

0 jesli k = 1, 2, . . . , n − 1 (48) Wtedy w _i > 0 dla i = 1, 2, . . . , n jak również

Z b a

ω (x) p (x) dx =

n

X

i=1

w _i p (x _i ) (49)

dla wszystkich wielomianów p ∈ Π _2n−1 .

Jeśli na odwrót, dla pewnych liczb w _i , x _i dla i = 1, . . . , n prawdziwa jest równość (49) dla wszystkich wielomianów p ∈ Π _2n−1 , to wówczas x _i są zerami wielomianu p _n , a w i spełniają układ równań (48).

Nie istnieją liczby x _i , w _i dla i = 1, . . . , n takie, żeby zachodziło (49) dla wszystkich wielomianów p ∈ Π _2n .

Dla funkcji wagowej w(x) :≡ 1 i przedziału [−1, 1] wielomianami ortogonalnymi są wtedy wielomiany Legendre’a (kwadratura Gaussa-Legendre’a)

p _k (x) := k!

(2k)!

d ^k dx ^k

 x ² − 1 ^{ k} (50)

dla k = 0, 1, . . ..

Przykład. Obliczyć współczynniki w _i x i dla całki Z 1

−1 x ³ + 2x ² + x + 1 (51)

Wartość dokładna całki to 10/3. Możemy wyliczyć 2n − 1 = 3, czyli n = 2 Układ równań to

p ₀ (x ₁ ) w ₁ + p ₀ (x ₂ ) w ₂ = (p ₀ , p ₀ ) (52) p ₁ (x ₁ ) w ₁ + p ₁ (x ₂ ) w ₂ = 0 (53) Wielomiany Legendre’a to

p ₀ (x) = 1 (54)

p 1 (x) = x (55)

p ₂ (x) = 1 2

 3x ² − 1 ^ (56)

Wartości x ₁ i x ₂ to zera wielomianu p ₂ (x), a zatem x ₁ = −1/ √

3, x ₂ = 1/ √ 3. A zatem układ równań to

1w ₁ + 1w ₂ = 2 (57)

− 1

√ 3 w 1 + 1

√ 3 w 2 = 0 (58)

Rozwiązanie: w ₁ = 1, w ₂ = 1. A więc ostatecznie Z 1

−1

x ³ + 2x ² + x + 1 = f (x ₁ ) + f (x ₂ ) =



− 1

√ 3

 3

+ 2



− 1

√ 3

 2

− 1

√ 3 + 1 (59)

+

 1

√ 3

 3

+ 2

 1

√ 3

 2

+ 1

√ 3 + 1 = 10/3 (60)

Alternatywnie możemy wyliczyć wartości współczynników ze wzoru

a _k = − 2

(n + 2) p _n+2 (x _k ) p ⁰ _n+1 (x _k ) (61) gdzie x _k są zerami wielomianu p _n+1 (x).

Innymi możliwymi wielomianami ortogonalnymi są

1. wielomiany Czebyszewa T _n (x), przedział [−1, 1] z funkcją wagową w (x) = ^ 1 − x ^{2  −}

1

2

(62)

2. wielomiany Laguerre’a L _n (x), przedział [0, ∞] z funkcją wagową w(x) = e ^−x

3. wielomiany Hermite’a H _n (x), przedział [−∞, ∞] z funkcją wagową w(x) = e ^−x

Konwersja przedziału [a, b] do [−1, 1] za pomocą całkowania przed podstawienie.

Mamy po zróżniczkowaniu wzoru na x względem t dx = b − a

2 dt (63)

Po podstawieniu powyższego wzoru i wzoru na x względem t, oraz zastąpieniu granicy całkowania (też po podstawieniu) otrzymujemy

Z b a

f (x) dx = b − a 2

Z 1

−1 f

 b − a

2 t + a + b 2



dt (64)

Dla n = 1. Mamy P ₂ (x) = 1/2(3x ² −1), P ₃ (x) = 1/2(5x ³ −3x). Węzły są pierwiastka- mi wielomianu P ₂ (x), czyli x = ±1/ √

3. Zaś oba współczynniki są równe 1, przykładowe wyprowadzenie

a 0 = −2

3/2(5x ³ ₀ − 3x ₀ ) ∗ 3x ₀ = 2 3/2(−5/ √

3 ³ + 3/ √

3) ∗ 3/ √ 3

(65)

= 2

3/2(−5/(3 √

3) + 3/ √

3) ∗ 3/ √

3 = 2

3/2(4/(3 √

3)) ∗ 3/ √

3 = 2

6(1/( √

3)) ∗ 1/ √ 3 (66)

= 2/2 = 1 (67)

Przykład: Obliczyć całkę

Z 1 0

dx

1 + x (68)

Na początku zamieniamy problem na przedział [−1, 1] i mamy całkę 1

2 Z 1

−1

dt

1/2t + 3/2 ≈ 0.69315 (69)

Obliczamy t ₀ = −1/ √

3, podobnie t ₁ = 1/ √

3 Przybliżenie całki jest równe a ₀ f (t ₀ )+a ₁ f (t ₁ ) = 1

−1/ ^ 2 √

3 ^ + 3/2

+ 1

1/ ^ 2 √

3 ^ + 3/2

= 0.825542+0.559073 = 1.384615 (70) Po podzieleniu przez 2 otrzymujemy 0.6923075.

Przykład. Chcemy znaleźć wartość całki Z 1

−1

x ³ + 1dx = 2 (71)

A zatem

f (x) ≡ x ³ + 1 (72)

ω (x) ≡ 1 (73)

Jest to wielomian stopnia 3, a zatem 2n − 1 = 3, czyli n = 2. A zatem szukamy Z 1

−1

1 ^ x ³ + 1 ^ dx =

2

X

i=1

w i

 x ³ _i + 1 ^ = w ₁ ^ x ³ ₁ + 1 ^ + w ₂ ^ x ³ ₂ + 1 ^ (74)

Musimy podstawić odpowiednie wartości miejsc zerowych wielomianu Legendre’a 2 stop- nia, oraz współczynników w _i rozwiązując układ równań. Miejsca zerowe to x ₁ = −1/ √

3, x 2 = 1/ √

3. A więc układ równań będzie miał postać:

p ₀ (x ₁ ) w ₁ + p ₀ (x ₂ ) w ₂ = (p ₀ , p ₀ ) (75) p ₁ (x ₁ ) w ₁ + p ₁ (x ₂ ) w ₂ = 0 (76) Wiemy, że p ₀ (x) = 1, p ₁ (x) = x. Po podstawieniu

1 · w ₁ + 1 · w ₂ = Z 1

−1

1dx = 2 (77)

x 1 w 1 + x ₂ w 2 = 0 (78)

Czyli

w 1 = 2 − w ₂ (79)

x ₁ (2 − w ₂ ) + x ₂ w ₂ = 0 (80)

w 2 = −2x ₁

x 2 − x ₁ = 1 (81)

w 1 = 1 (82)

A więc po podstawieniu otrzymujemy Z 1

−1

1 ^ x ³ + 1 ^ dx = ^ x ³ ₁ + 1 ^ + ^ x ³ ₂ + 1 ^ = 2 (83) Obliczenie błędu.

Twierdzenie 5. Jeśli funkcja f ma na przedziale [a, b] ciągłą pochodną rzędu 2n, to

Z b a

w (x) f (x) dx −

n

X

i=1

w _i f (x _i ) = f ⁽²ⁿ⁾ (ξ)

(2n)! (p _n , p _n ) (84)

gdzie ξ ∈ (a, b).

Kwadratury Gaussa-Czebyszewa Po angielsku Chebyshev-Gauss quadrature, http:

//en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature. Szczególny przy- padek, gdy wartości współczynników a _i są stałymi. Kwadratury Gaussa-Czebyszewa dla przedziału [−1, 1]

a _k = π

n + 1 (85)

x _k = cos (2k + 1) π

2n + 2 (86)

Są to miejsca zerowe wielomianów Czebyszewa T _n (x), z funkcją wagową

√ 1

1 − x ² (87)

Przykładowe pierwiastki dla n = 1: x ₀ = 1/ √

2, x ₁ = −1/ √

2. Dla n = 2, x ₀ = √ 3/2, x 1 = 0, x ₂ = − √

3/2.

Przykład: Obliczyć całkę

Z 1 0

dx

1 + x (88)

Na początku zamieniamy problem na przedział [−1, 1] i mamy całkę 1

2 Z 1

−1

dt

1/2t + 3/2 (89)

Następnie musimy mieć funkcję wagową pod całką:

1 2

Z 1

−1

√ 1 1 − t ²

√

1 − t ² dt

1/2t + 3/2 (90)

Obliczamy t ₀ = 1/ √

2, podobnie t ₁ = −1/ √

2. I otrzymujemy:

π 4

q 1 − t ² ₀ 1/2t ₀ + 3/2 + π

4 q

1 − t ² ₁

1/2t ₁ + 3/2 = 0.7840 (91)

Dla n = 2 mamy t ₀ = √

3/2, t ₁ = 0, t ₂ = − √

3/2 i otrzymujemy π

6 q

1 − t ² ₀ 1/2t ₀ + 3/2 + π

6 q

1 − t ² ₁ 1/2t ₁ + 3/2 + π

6 q

1 − t ² ₂

1/2t ₂ + 3/2 = 0.7299

Kwadratura Czebyszewa Po angielsku Chebyshev quadrature, http://mathworld.

wolfram.com/ChebyshevQuadrature.html. Wszystkie współczynniki a _i są stałymi i są sobie równe oraz funkcja wagowa równa 1. Rozpatrywany przedział [−1, 1].

Przyjmujemy, że

a = 2

n (92)

Węzły wyznaczamy za pomocą wzoru (4) tak aby był spełniony dla funkcji x, x ² , . . . , x ⁿ , czyli

Z 1

−1

f (x) dx = 2 n

n

X

i=1

f (x _i )

A więc otrzymujemy układ równań. Z układu równań możemy wyznaczyć węzły x _i . Dla n = 3, układ równań na węzły jest następujący, całka funkcji x jest równa 0, zatem

t 1 + t ₂ + t ₃ = 0 (93)

Całka funkcji x ² jest równa 2/3 zatem

t ² ₁ + t ² ₂ + t ² ₃ = 1 (94)

Całka funkcji x ³ jest równa 0 zatem

t ³ ₁ + t ³ ₂ + t ³ ₃ = 0 (95)

Rozwiązanie t ₁ = −1/ √

2, t ₂ = 0, t ₃ = 1/ √ 2.

Przykład obliczyć całkę

Z 1 0

dx

1 + x (96)

dla n = 3. Najpierw przekształcamy całkę do przedziału (−1, 1) 1

2 Z 1

−1

dt

1/2t + 3/2 (97)

I otrzymujemy 1 2

Z 1

−1

dt

1/2t + 3/2 = 1

3 (f (x ₁ ) + f (x ₂ ) + f (x ₃ )) = 0.6928 (98) 1.3 Całkowanie adaptacyjne

Dla każdego przedziału sprawdzany jest błąd całkowania. Błąd całkowania jest spraw- dzany w ten sposób, że każdy przedział dzielony jest na dwa przedziały. Sprawdzana jest różnica między przybliżoną wartością całki w danym przedziale i sumą przybliżo- nych wartości całek w dwóch przedziałach. Jeśli ta różnica jest duża oznacza, to że w danym przedziale przybliżona wartość całki nie jest dokładna, wtedy przedział ten dzielony jest na dwa przedziały. Sprawdzanie jest powtarzane dla każdego przedziału.

2 Zadania

2.1 Zadania na 3.0

• Wyprowadzić wartość dokładną i przybliżone kwadraturami Newtona-Cotesa dla n = 1, 2, 3, gdzie n to stopień wielomianu interpolacyjnego całki

Z 2 1

 x ² + 3x ^ dx . (99)

Obliczyć błędy względne procentowe. Naszkicować wykres funkcji podcałkowej i jej przybliżeń wraz z punktami.

2.2 Zadania na 4.0

• Obliczyć całkę

Z 1 0

2 (1 + x)

p x (1 − x) dx (100)

za pomocą kwadratury Gaussa-Czebyszewa dla n = 2. Obliczyć błędy względne procentowe. Naszkicować wykres funkcji podcałkowej wraz z punktami.

• Za pomocą kwadratury Czebyszewa obliczyć całkę Z 2

0

dx

1 + x ² (101)

Obliczyć błędy względne procentowe. Naszkicować wykres funkcji podcałkowej i jej przybliżeń wraz z punktami.

2.3 Zadania na 5.0

• Obliczyć całkę

Z 1 0

2 (1 + x)

p x (1 − x) dx (102)

kwadraturami Gaussa-Legendre’a dla n = 1.

A Tabela wielomianów Legendre’a

1 x

1

2 3x ² − 1
^{−1 √}

3 , ^{√ 1}

3 1

2 5x ³ − 3x
0, − ^{q 3} ₅ , ^{q 3} ₅

1

8 35x ⁴ − 30x ² + 3
−0.33998, 0.33998, −0.86114, 0.86114

(103)