Metody numeryczne
Instytut Sterowania i Systemów Informatycznych
Wydział Elektrotechniki, Informatyki i Telekomunikacji
Uniwersytet Zielonogórski
Elektrotechnika stacjonarne-dzienne pierwszego stopnia z tyt. inżyniera
Informatyka stacjonarne-dzienne drugiego stopnia z tyt. magistra inżyniera
Całkowanie numeryczne
Laboratorium, prowadzący: mgr inż. Błażej Cichy
Rok akademicki 2010/2011
1 Uwagi teoretyczne
1.1 Kwadratury Newtona-Cotesa
Kwadratury Newtona-Cotesa polegają na całkowaniu wielomianu Lagrange a, który
jest reprezentowany następująco:
n
f(x) H" Pn(x) = fkLn,k(x) (1)
k=0
Ogólnie metoda całkowania Newtona-Cotesa posiada następującą postać:
xn xn
f(x)dx H" Pn(x)dx =
x0 x0
xn xn
n n
fkLn,k(x) dx = fkLn,k(x)dx =
(2)
x0 k=0 k=0 x0
xn
n n
= Ln,k(x)dx fk = wkfk
k=0 x0 k=0
Pierwsze cztery kwadratury są następujące:
x1
h
f(x)dx H" (f0 + f1) regua trapezu
2
x0
x2
h
f(x)dx H" (f0 + 4f1 + f2) regua Simpsona
3
x0
(3)
x3
3h
f(x)dx H" (f0 + 3f1 + 3f2 + f3) regua Simpsona3/8
8
x0
x4
2h
f(x)dx H" (7f0 + 32f1 + 12f2 + 32f3 + 7f4) regua Boole a
45
x0
1
Całkowanie numeryczne 2
1.2 Wyprowadzenie metody Simpsona
Rozpoczynamy od zdefiniowania wielomianu Lagrange a (lista zadań o interpolacji)
stopnia drugiego:
(x - x1)(x - x2) (x - x0)(x - x2) (x - x0)(x - x1)
P2(x) = f0 + f1 + f2 (4)
(x0 - x1)(x0 - x2) (x1 - x0)(x1 - x2) (x2 - x0)(x2 - x1)
Całkujemy wielomian wyciągając przed znak całki wartości funkcji:
x2 x2 x2
(x - x1)(x - x2) (x - x0)(x - x2)
P2(x) = f0 dx + f1 dx
(x0 - x1)(x0 - x2) (x1 - x0)(x1 - x2)
x0 x0 x0
(5)
x2
(x - x0)(x - x1)
+f2 dx
(x2 - x0)(x2 - x1)
x0
Aby wykonać całkowanie należy wprowadzić nowe oznaczenia, x = x0 +ht oraz dx = h dt.
Zmieniamy zakres całkowania od 0 to t Ponieważ węzły są równo rozmieszczone, czyli
xk = x0+kh co oznacza, że różnicę elementów xk -xj możemy zastąpić przez (k-j)h oraz
różnicę x - xk przez h(t - k). Po podstawieniu możemy wykonać dalsze przekształcenia:
x2 2 2 2
P2(x) H" f0 h(t-1)h(t-2)h dt + f1 h(t-0)h(t-2)h dt + f2 h(t-0)h(t-1)h dt =
(-h)(-2h) (h)(-h) (2h)(h)
x0 0 0 0
2 2 2
= f0 h (t2 - 3t + 2)dt - f1h (t2 - 2t)dt + f2 h (t2 - t)dt =
2 2
0 0 0 (6)
t=2 t=2 t=2
3t2 t3 t2
= f0 h t3 - + 2t - f1h - t2 + f2 h t3 -
=
2 3 2 3 2 3 2
t=0 t=0 t=0
h
= f0 h 2 - f1h-4 + f2 h 2 = (f0 + 4f1 + f2)
2 3 3 2 3 3
1.3 Złożona metoda trapezów
Znaczące polepszenie numerycznego całkowania uzyskamy, jeśli poszczególne metody
zostaną zastosowane dla podprzedziałów. Załóżmy że mamy pięć punktów: x0, . . . , x4. Gdy
będziemy całkować każdy podprzedział oddzielnie do możemy zastosować np.: metodę
trapezów:
x4 x1 x2 x3 x4
f(x)dx = f(x)dx + f(x)dx + f(x)dx + f(x)dx H"
x0 x0 x1 x2 x3
(7)
h h h h
H" (f0 + f1) + (f1 + f2) + (f2 + f3) + (f3 + f4) =
2 2 2 2
h
= (f0 + 2f1 + 2f2 + 2f3 + f4)
2
T
Błąd En (f) w złożonej metodzie trapezów można wyznaczyć ze wzoru
-h2(b - a)
T
En (f) = f (cn)
12
gdzie cn jest nieznanym punktem w przedziale [a, b], zaÅ› h = (b - a)/n. Natomiast osza-
cowanie błędu dla dużych wartości n określone jest wzorem
-h2
T
En (f) H" (f (b) - f (a)).
12
Całkowanie numeryczne 3
1.4 Złożona metoda Simpsona
Podobnie jak w poprzedniej metodzie również dla metody Simpsona możemy polepszyć
jakość całkowania, jeśli będziemy stosować dla mniejszych przedziałów.
x4 x2 x4
f(x)dx = f(x)dx + f(x)dx H"
x0 x0 x2
(8)
h h
H" (f0 + 4f1 + f2) + (f2 + 4f3 + f4) =
3 3
h
(f0 + 4f1 + 2f2 + 4f3 + f4)
3
S
Błąd En (f) w metodzie Simpsona można wyznaczyć ze wzoru
-h4(b - a)
S
En (f) = f(4)(cn)
180
gdzie cn jest nieznanym punktem w przedziale [a, b], h = (b-a)/n, zaÅ› f(4) oznacza czwartÄ…
pochodną funkcji f. Oszacowanie błędu dla dużych wartości n określone jest wzorem
-h4
S
En (f) H" (f (b) - f (a))
180
2 Metoda Gaussa
2.1 Wstęp
Stosowanie metod opartych na kwadraturach Newtona-Cotesa niesie ze sobÄ… pewien
błąd a ponadto dla osiągnięcia odpowiedniej dokładności wymaga wykonania obliczeń
w wielu punktach. Powstaje więc pytanie: Czy jest możliwe wyznaczenie przybliżonej
wartości całki wykonując obliczenia dla małej liczby punktów? Ponadto, czy metoda taka
może dawać idealne (pomijając błędy numeryczne obliczeń) wyniki dla pewnych klas
funkcji? OdpowiedziÄ… na te pytania jest metoda Gaussa.
2.2 Sformułowanie problemu
Załóżmy, że należy wyznaczyć wartość całki
1
I(f) = f(x)dx
-1
przy pomocy przybliżonej metody
n
I(f) H" In(f) = wjf(xj)
j=1
dodatkowo zakładając, że wagi wj oraz węzły xj będą dobrane tak, aby In(f) = I(f) dla
wszystkich wielomianów stopnia 2n - 1. Błąd musi zatem wynosić zero (a więc I(f) =
In(f)) dla wszystkich wielomianów stopnia niższego lub równego 2n - 1. Konstrukcję
metody przeprowadza się dla wielomianów postaci f(x) = xk dla k = 0, 1, 2, 3, . . . 2n - 1.
Całkowanie numeryczne 4
Pozwala to na wyznaczenie niewiadomych wj oraz xj. Dokładność tej metody całkowania
jest wprost proporcjonalna do wartości n.
Dla metod całkowania definiuje się ponadto stopień precyzji. Metoda całkowania po-
siada stopień precyzji równy k, jeśli daje idealnie dokładny wynik dla wszystkich wie-
lomianów stopnia d" k, zaś dla wielomianu stopnia k + 1, błąd jest niezerowy. Stopień
precyzji metody Gaussa stopnia n jest zatem równy 2n - 1.
2.3 Przypadek n = 1
1
n
I(f) = f(x)dx H" I1(f) = wjf(xj) = w1f(x1) (9)
j=1
-1
Z założenia metoda ma być idealnie dokładna dla wielomianu stopnia 0, a więc f(x) =
1. Podstawiając f(x) = 1 do wzoru (??) i uwzględniając założenie o zerowym błędzie
metody, otrzymuje siÄ™
1
I(f) = (1)dx = I1(f) = w1 " 1
-1
Aby równość była spełniona, należy tak dobrać współczynnik w1, aby
1
(1)dx = w1
-1
Całka po lewej stronie wynosi
1
(1)dx = x|x=1 - x|x=-1 = 1 + 1 = 2
-1
a więc
1
(1)dx = 2 = w1
-1
Współczynnik w1 wynosi więc 2. Pozostało wyznaczenie wartości x1. Aby to uczynić,
możliwe jest wymuszenie, aby błąd dla wielomianu stopnia 1 (a więc f(x) = x) był
zerowy. Oznacza to, że
1
I(f) = f(x)dx = I1(f) = w1 " f(x1)
-1
Podstawiając f(x) = x oraz wcześniej wyznaczony współczynnik w1 = 2
1
I(f) = (x)dx = 2 " x1
-1
Całkowanie numeryczne 5
Całka po lewej stronie wynosi
1
x2 x2 1 1
I(f) = (x)dx = - = - = 0
2 2 2 2
x=1 x=-1
-1
a zatem
I(f) = 0 = 2 " x1
Z tego wynika, że x1 = 0.
PodsumowujÄ…c, dla n = 1 i dowolnej funkcji f zachodzi
1
f(x)dx H" 2f(0)
-1
Metoda jest idealnie dokładna dla wszystkich wielomianów stopnia n d" 1, a więc jej
stopień precyzji wynosi 1. A
atwo pokazać, że dla wielomianów wyższych stopni (np. dla
x2) metoda nie daje idealnie dokładnego wyniku.
2.4 Przypadek n = 2
Ponieważ n = 2, więc obliczenia należy przeprowadzić dla wielomianów stopnia 0, 1, 2, 3.
Podstawiając do wzoru ogólnego na całkowanie numeryczne metodą Gaussa otrzymuje się
1
n
I(f) = f(x)dx H" I2(f) = wjf(xj) = w1f(x1) + w2f(x2) (10)
j=1
-1
Z założenia metoda ma być idealnie dokładna dla wielomianu stopnia 0, a więc f(x) =
1. Podstawiając f(x) = 1 do wzoru (??) i uwzględniając założenie o zerowym błędzie
metody, otrzymuje siÄ™
1
I(f) = (1)dx = I2(f) = w1 " 1 + w2 " 1
-1
upraszczajÄ…c
2 = w1 + w2
Z założenia metoda ma być również idealnie dokładna dla wielomianu stopnia 1, a
więc f(x) = x. Podstawiając f(x) = x do wzoru (??) i uwzględniając założenie o zerowym
błędzie metody, otrzymuje się
1
I(f) = (x)dx = w1 " x1 + w2 " x2
-1
z czego wykonując całkowanie otrzymuje się
0 = w1 " x1 + w2 " x2
Całkowanie numeryczne 6
Z założenia metoda ma być również idealnie dokładna dla wielomianu stopnia 2, a więc
f(x) = x2. Podstawiając f(x) = x2 do wzoru (??) i uwzględniając założenie o zerowym
błędzie metody, otrzymuje się
1
I(f) = (x2)dx = w1 " x2 + w2 " x2
1 2
-1
Całka po lewej stronie wynosi
1
1 1 2
I(f) = (x2)dx = x3 - x3 =
3 3 3
x=1 x=-1
-1
Z tego wynika, że
2
= w1 " x2 + w2 " x2
1 2
3
Z założenia metoda ma być również idealnie dokładna dla wielomianu stopnia 3, a więc
f(x) = x3. Podstawiając f(x) = x2 do wzoru (??) i uwzględniając założenie o zerowym
błędzie metody, otrzymuje się
1
I(f) = (x3)dx = w1 " x3 + w2 " x3
1 2
-1
Całka po lewej stronie wynosi
1
1 1
I(f) = (x3)dx = x4 - x4 = 0
4 4
x=1 x=-1
-1
Z tego wynika, że
0 = w1 " x3 + w2 " x3
1 2
Uwzględnienie wszystkich obliczeń dla f(x) = 1, x, x2, x3 prowadzi do układu równań
2 = w1 + w2
0 = w1x1 + w2x2
2
= w1x2 + w2x2
1 2
3
0 = w1x3 + w2x3
1 2
Można pokazać, że rozwiązaniem tego układu równań są liczby
w1 = 1
w2 = 1
"
3
x1 = -
3
"
3
x2 =
3
oraz
w1 = 1
w2 = 1
"
3
x1 =
3
"
3
x2 = -
3
Całkowanie numeryczne 7
PodsumowujÄ…c, dla n = 2 i dowolnej funkcji f zachodzi
1
" "
3 3
f(x)dx H" f(- ) + f( )
3 3
-1
Metoda jest idealnie dokładna dla wszystkich wielomianów stopnia n d" 3, a więc jej
stopień precyzji wynosi 3. A
atwo pokazać, że dla wielomianów wyższych stopni (np. dla
x4) metoda nie daje idealnie dokładnego wyniku.
2.5 Metoda Gaussa wyższych stopni
Analogicznie do poprzednio opisanych przypadków, dla n > 2 zadanie polega na wy-
znaczeniu wartości x1, x2, . . . , xn oraz w1, w2, . . . , wn tak, aby wyrażenie
1
n
f(x)dx = wjf(xj)
j=1
-1
było prawdziwe dla f(x) = 1, x, x2, x3, . . . , x2n-1 Wykonanie całkowań analogicznych dla
pokazanych poprzednio prowadzi do następującego układu równań
Å„Å‚
ôÅ‚ 2 = w1 + w2 + · · · + wn
ôÅ‚
ôÅ‚
ôÅ‚
ôÅ‚ 0 = w1x1 + w2x2 + · · · + wnxn
ôÅ‚
ôÅ‚
2
ôÅ‚
ôÅ‚ = w1x2 + w2x2 + · · · + wnx2
1 2 n
òÅ‚
3
0 = w1x3 + w2x3 + · · · + wnx3
1 2 n (11)
ôÅ‚
. .
ôÅ‚
. .
ôÅ‚
. .
ôÅ‚
ôÅ‚
ôÅ‚
2
ôÅ‚
= w1x2n-2 + w2x2n-2 + · · · + wnx2n-2
ôÅ‚
1 2 n
2n-1
ôÅ‚
ół
0 = w1x2n-1 + w2x2n-1 + · · · + wnx2n-1
1 2 n
Metoda taka ma stopień precyzji równy 2n - 1. Rozwiązanie takiego równania jest zada-
niem bardzo trudnym. Z tego powodu wartości xk oraz wk często umieszcza się w tabelach.
n xk wk
n xk wk
6 Ä…0.9324695142 0.1713244924
1 0.0 2.0 Ä…0.6612093865 0.3607615730
2 Ä…0.5773502692 1.0 Ä…0.2386191861 0.4679139346
3 Ä…0.7745966692 0.5555555556 7 Ä…0.9491079123 0.1294849662
0.0 0.8888888889 Ä…0.7415311856 0.2797053915
4 Ä…0.8611363116 0.3478548451 Ä…0.4058451514 0.3818300505
Ä…0.3399810436 0.6521451549 0.0 0.4179591837
5 Ä…0.9061798459 0.2369268851 8 Ä…0.9602898565 0.1012285363
Ä…0.5384693101 0.4786286705 Ä…0.7966664774 0.2223810345
0.0 0.5688888889 Ä…0.5255324099 0.3137066459
Ä…0.1834336425 0.3626837834
2.6 Przystosowanie metody Gaussa dla dowolnych przedziałów
Metoda Gaussa nadaje się wyłącznie do wyznaczania całek w przedziale [-1, 1]. Znacz-
nie częściej w zastosowaniach praktycznych konieczne jest jednak wyznaczenie wartości
Całkowanie numeryczne 8
całki
b
f(x)dx
a
StosujÄ…c podstawienie
b + a + t(b - a)
x =
2
dla t " [-1, 1]. Prowadzi to do całki w postaci:
1
b - a b + a + t(b - a)
f dt
2 2
-1
która może zostać wyznaczona przy pomocy metody Gaussa.
3 Zadania
1. Wyznaczyć wzory dla pierwszych czterech rzędów kwadratury Newtona-Cotesa.
2. Wyznaczyć wzory dla metod złożonych dla pierwszych czterech rzędów kwadratury
Newtona-Cotesa.
3. Stosując kwadratury Newtona-Cotesa stopnia 1 - 4 wyznaczyć wartości całek ozna-
czonych z następujących funkcji:
(a) 3x3 - 1
(b) sin(x)
(c) e-x
Przedział całkowania można być dowolny ale wspólny dla wszystkich trzech funk-
cji. Ocenić na podstawie reszty kwadratury oraz wartości wyznaczanej analitycznie
zależność dokładności obliczeń od stopnia kwadratury.
4. Zastosować metody złożone trapezu oraz Simpsona dla wymienionych w poprzednim
punkcie funkcji. Czy wyniki uległy poprawie?
5. Zaimplementować metodę Simpsona
6. Metoda prostokątów polega na zastąpieniu całki z funkcji f całką (polem pod wy-
kresem) z prostokąta o wysokości f(x1) i szerokości (x0 - x1). Zaimplementować tą
metodÄ™.
7. Udowodnić, że metoda Gaussa stopnia 2 daje idealnie dokładny (pomijając błędy
numeryczne) wynik przy całkowaniu wielomianu f(x) = ax2 + bx + c, a = 0.

8. Wyznaczyć błąd metody Gaussa stopnia 2 przy całkowaniu funkcji f(x) = x4.
9. Zaimplementować metodę Gaussa stopnia 1, 2, 3 i 8 dla dowolnego przedziału cał-
kowania.
Całkowanie numeryczne 9
10. Wyznaczyć korzystając z metody Gaussa stopnia 2, 3 oraz 8 wartość całki I(f) =
1
1
dx. Wskazówka: Przed wykonywaniem obliczeń zmodyfikować granice całko-
1+x
0
wania.
11. Korzystając z wyników zadań poprzednich porównać dokładność poznanych metod.
Literatura
[1] Bjärck Ake i Dahlquist Germund. Metody numeryczne. PWN, Warszawa, 1987.
[2] Jerzy Brzózka i Lech Dorobczyński. Programowanie w MATLAB. Warszawa,
Wydanie I, 1998.
[3] Zenon Fortuna, Bohdan Macukow i Janusz WÄ…sowski. Metody numeryczne. WNT,
Warszawa, 1995.
[4] Jerzy Klamka i in. Metody numeryczne. Politechnika ÅšlÄ…ska, Gliwice, 1998.
[5] David Kincaid i Ward Cheney. Analiza numeryczna. WNT, Warszawa, 2006.
[6] Anna Kamińska i Beata Pańczyk. Matlab. Ćwiczenia z . . . , Przykłady i zadania.
Warszawa, Wydanie I, 2002.
[7] Wanat Kazimierz. Algorytmy numeryczne. Helion, Gliwice, 1994.
[8] Bogumiła Mrozek i Zbigniew Mrozek. MATLAB i Simulink. Poradnik użytkownika.
Wydanie II, 2004.
[9] Jurij Povstenko. Wprowadzenie do metod numerycznych. Akademicka Oficyna
Wydawnicza EXIT, Warszawa, Wydanie drugie poprawione i uzupełnione, 2005.
[10] Rudra Pratap. MATLAB 7 dla naukowców i inżynierów. PWN, 2007.
[11] Wiesława Regel. Wykresy i obiekty graficzne w MATLAB. Warszawa, Wydanie I,
2003.
[12] Marcin Stachurski. Metody numeryczne w programie Matlab. Warszawa, Wydanie I,
2003.