Całkowanie numeryczne
Piotr Modliński
20 października 2010
Spis treści
1
Wstęp
1
2
Definicje i oznaczenia
1
3
Kwadratury proste Newtona-Cotesa
2
3.1
Formuła prostokątów . . . . . . . . .
2
3.2
Formuła trapezów . . . . . . . . . . .
2
4
Kwadratury złożone Newtona-Cotesa
2
4.1
Złożona formuła prostokątów
. . . .
3
4.1.1
Oszacowanie błędu . . . . . .
3
4.2
Złożona formuła trapezów . . . . . .
3
4.2.1
Oszacowanie błędu . . . . . .
3
5
Uwagi na koniec
4
6
Zadania kontrolne
4
6.1
Zadanie 1 . . . . . . . . . . . . . . .
4
6.2
Zadanie 2 . . . . . . . . . . . . . . .
4
1
Wstęp
Wielokrotnie w naukach technicznych występuje ko-
nieczność obliczenia całki oznaczonej. O ile jest
to rozwiązanie bardzo wygodne z matematycznego
punktu widzenia, wpisuje się idealnie w analizę funk-
cji, można snuć teoretyczne rozważania, o tyle jest
często bardzo trudne do wykonania w praktyce. Z
definicji obliczenie wartości całki oznaczonej wy-
maga wyznaczenia funkcji pierwotnej, a następnie
obliczenie jej wartości na granicy przedziałów. Wła-
śnie to określenie postaci funkcji pierwotnej jest czę-
sto największym problemem. Najprostszą sytuacją
jest taka, w której po prostu nie potrafimy tej funk-
cji wyznaczyć. Zdarza się również, że nie jest to
wina naszej niekompetencji, ale funkcja pierwotna
jest w ogóle niemożliwa do wyznaczenia – np. nie
ma postaci analitycznej. Może zdarzyć się, że ana-
lityczne obliczenia już po wyznaczeniu funkcji pier-
wotnej dają bardzo duży błąd numeryczny, zatem są
z punktu widzenia dalszego przetwarzania nieprzy-
datne. Wreszcie ostatnią grupę stanowią funkcje ta-
blicowane, w których postać funkcji podcałkowej jest
aproksymowana za pomocą listy wartości pobieranej
np. z pomiarów. W tej sytuacji także nie ma możli-
wości określenia postaci analitycznej.
Zamiast robić to analitycznie, w przypadku całki
oznaczonej można wykorzystać do obliczeń kompu-
ter. W sposób numeryczny możliwe jest wyznacze-
nie przybliżonej wartości całki oznaczonej na zada-
nym przedziale jeśli tylko potrafimy określić war-
tość funkcji podcałkowej dla konkretnych wartości.
W dalszej części zajmiemy się przykładowymi me-
todami określania całki, a także oszacowaniem do-
kładności takiego przybliżenia.
2
Definicje i oznaczenia
Kwadraturą (czy kwadraturą numeryczną) okre-
ślamy przedstawioną niżej metodę obliczania całek
oznaczonych.
Przyjmijmy następujące oznaczenia:
• I(f ) =
R
b
a
f (x)dx – całka oznaczona
• S(f ) =
P
n
i=0
f (x
i
)A
i
– formuła liniowa (kombi-
nacja liniowa wartości funkcji w punktach x
i
),
gdzie A
i
jest tzw. współczynnikiem formuły li-
niowej, oraz x
i
∈< a, b >
• R(f ) = I(f ) − S(f ) – błąd formuły
Można zauważyć wówczas, że I(f ) ≈ S(f ), a przy-
bliżenie jest tym dokładniejsze, im większa jest
liczba n.
Liczbę n nazywamy rzędem kwadratury – kwadra-
tura jest dokładna dla wszystkich wielomianów stop-
nia mniejszego od n, oraz istnieje wielomian stopnia
n, dla którego nie jest dokładna.
1
3
Kwadratury
proste
Newtona-Cotesa
Tzw. kwadratury Newtona-Cotesa polegają na przy-
bliżeniu funkcji podcałkowej za pomocą wielomianu
określonego stopnia, czyli:
Z
b
a
f (x)dx ≈
Z
b
a
W (x)dx
gdzie W (x) jest tzw. wielomianem interpolacyjnym
Lagrange’a funkcji f opartym na równoodległych
węzłach x
i
= a + ih, gdzie i = 0, 1, . . . , n, zaś
h =
b−a
n
.
Dla zastosowań praktycznych (nie wchodząc zbyt
głęboko w teorię) wykorzystywanych na zajęciach
będziemy używali tylko dwóch kwadratur:
• Kwadraturę stopnia 0 – tzw. formułę prostoką-
tów
• Kwadraturę stopnia 1 – tzw. formułę trapezów
Wykorzystywane są także kwadratury wyższych rzę-
dów (np. rzędu 2 – tzw. formuła parabol wykorzy-
stująca 3 kolejne punkty), oraz metody Monte-Carlo
(opierające się o losowe próbkowanie i szacowanie
pola pod i nad wykresem funkcji), jednak nie bę-
dziemy się nimi dalej zajmować skupiając się jedynie
na dobrym zrozumieniu dwóch najprostszych.
3.1
Formuła prostokątów
Funkcję przybliżamy wielomianem stopnia zerowego,
a więc funkcją stałą. Zgodnie z definicją wykorzystu-
jemy do tego tylko jeden węzeł x
0
= a, toteż nasze
przybliżenie ma postać:
Z
b
a
f (x)dx ≈
Z
b
a
f (a)dx = f (a) · (b − a)
Schematycznie przedstawiono to na rysunku 1, gdzie
zakreskowano dokładną wartość całki, natomiast
przybliżenie oznaczono kolorem szarym.
Oczywiście nic nie stoi na przeszkodzie, aby wy-
boru punktu węzłowego dokonać w inny sposób –
np. biorąc środek przedziału, bądź dowolny inny.
3.2
Formuła trapezów
Nieznacznie bardziej złożoną kwadraturą jest tzw.
formuła trapezów. W tym przypadku wielomian in-
terpolacyjny jest wielomianem stopnia 1, czyli do-
wolną prostą. Ponieważ mamy do czynienia z wie-
lomianem stopnia 1, potrzebujemy dwóch punktów,
x
b
f(x)
W(x)=f(a)
y
a
Rysunek 1: Kwadratura prosta stopnia 0 – formuła
prostokątów
by go wyznaczyć
1
. Przyjmujemy więc następujące
punkty węzłowe:
x
0
= a ⇒ f (x
0
) = f (a)
x
1
= b ⇒ f (x
1
) = f (b)
Wyznaczamy zatem funkcję liniową i całkujemy ją
na przedziale < a, b >. Powstaje trapez (stąd nazwa
formuły trapezów), którego pole można wyznaczyć
prosto ze wzoru:
Z
b
a
f (x)dx ≈
Z
b
a
W (x)dx =
f (a) + f (b)
2
(b − a)
Schematycznie rozwiązanie przedstawione zostało na
rysunku 2, gdzie zakreskowano dokładną wartość
całki, natomiast przybliżenie oznaczono kolorem sza-
rym.
x
b
f(x)
W(x)
y
a
Rysunek 2: Kwadratura prosta stopnia 1 – formuła
trapezów
4
Kwadratury
złożone
Newtona-Cotesa
Jak nietrudno zauważyć, metody te są tym skutecz-
niejsze, im granice całkowania są położone bliżej sie-
bie. Wykorzystując powyższe spostrzeżenie można
1
nie będę tego dowodził, jeśli ktoś jest zainteresowany, od-
syłam do własności układu równań Cramera
2
wykazać, że jeśli podzielimy przedział < a, b > na
m 0 równych części, otrzymamy następującą za-
leżność:
R(f )
m
X
i=1
R
i
(f )
gdzie:
R
i
(f ) =
R
b
i
a
i
f (x)dx −
R
b
i
a
i
W (x)dx
a
i
= a + (i − 1)h
b
i
= a + ih
h =
b−a
m
lub bardziej intuicyjnie – suma błędów zrobionych
na wielu małych przedziałach jest mniejsza niż błąd
zrobiony na jednym dużym przedziale (będącym
sumą małych).
Powyższe spostrzeżenia wykorzystać można do
budowy znacznie dokładniejszych kwadratur – tzw.
kwadratur złożonych. Podobnie, jak przy omawianiu
kwadratur prostych, tak i tutaj mamy do czynienia
z różnymi możliwymi podejściami, jednak skupimy
się na dwóch najprostszych – wielomianach stałych
i liniowych.
4.1
Złożona formuła prostokątów
Dla każdego podprzedziału (długości h) przybliżamy
funkcję pewną stałą – wartością funkcji podcałkowej
liczonej w początku przedziału (lub w dowolnym in-
nym ustalonym punkcie). Mamy wówczas:
S(f ) = h ·
m
X
i=1
f (a
i
)
Poglądowo przedstawiona została na rysunku 3.
x
b
f(x)
y
a
Rysunek 3: Kwadratura złożona stopnia 0 – formuła
prostokątów
4.1.1
Oszacowanie błędu
Zakładając, że funkcja podcałkowa jest klasy C
1
(czyli ma ciągłą pierwszą pochodną) na przedziale
< a, b > można oszacować, że
|R(f )| ¬
M
1
(b − a)
2
· h
gdzie:
M
1
= max
x∈<a,b>
|f
0
(x)|
Widać wyraźnie, że zagęszczając podział, popeł-
niany błąd zbiega do zera (lim
h→0
R(f ) = 0), co po-
twierdza intuicyjne spostrzeżenia, które poczynili-
śmy wyżej. Jeżeli nasza funkcja ma jeszcze lepsze
własności, tj. f ∈ C
2
< a, b > (ma ciągłą drugą
pochodną na rozważanym przedziale), można osza-
cować wielkość błędu znacznie dokładniej, ale poniż-
sze oszacowanie jest prawdziwe jedynie dla wyboru
punktów węzłowych pośrodku przedziału:
|R(f )| ¬
M
2
(b − a)
24
· h
2
gdzie:
M
2
= max
x∈<a,b>
|f
00
(x)|
4.2
Złożona formuła trapezów
Dokładna analogia do złożonej formuły prostokątów
– dla każdego podprzedziału przybliżamy funkcję
podcałkową funkcją liniową wyznaczaną na podsta-
wie wartości funkcji podcałkowej na granicach pod-
przedziału. Można wyznaczyć następujące przybli-
żenie:
S(f ) = h ·
f (a) + f (b)
2
+
m
X
i=2
f (a
i
)
!
Interpretacja graficzna przedstawiona została na ry-
sunku 4.
x
b
f(x)
y
a
Rysunek 4: Kwadratura złożona stopnia 1 – formuła
trapezów
4.2.1
Oszacowanie błędu
Można wykazać, że jeśli f ∈ C
2
< a, b > to błąd
przybliżenia spełnia następującą nierówność:
|R(f )| ¬
M
2
(b − a)
12
· h
2
3
5
Uwagi na koniec
Zbieżność
Z powyższych rozważań (w szczególności z osza-
cowania błędów) wynika, że ciąg kwadratur
Newtona-Cotesa S
m
(f ), gdzie m jest liczbą pod-
przedziałów, jest zbieżny dla wszystkich funkcji
f ∈< a, b > dla m → ∞ (czyli h → 0), a więc:
lim
m→∞
S
m
(f ) = I(f ) =
Z
b
a
f (x)dx
Warunek stopu
W praktyce często trudne jest szacowanie ilo-
ści iteracji niezbędnych do osiągnięcia żądanej
dokładności na podstawie przedstawionych teo-
retycznych oszacowań, lub jest ona silnie za-
wyżona (zwłaszcza jeśli nie umiemy dokładnie
ograniczyć pochodnych). Aby poradzić sobie z
tym problemem stosuje się często iteracyjne wy-
konywanie coraz dokładniejszych przybliżeń i
badanie różnic kolejnych, tj. dąży się do tego,
by:
S
m
(f ) − S
m+1
(f )
S
m+1
(f )
¬
gdzie jest pożądaną dokładnością. Podejście
to wystarcza dla całkowania znakomitej więk-
szości funkcji wykorzystywanych w praktyce.
W każdym kolejnym przybliżeniu wykorzystuje
się często wielokrotność liczby używanych pod-
przedziałów, co umożliwia wykorzystanie w ob-
liczeniach już wykorzystanych wartości.
6
Zadania kontrolne
6.1
Zadanie 1
Proszę obliczyć złożoną formułę prostokątów (dwu-
krotnie – raz z wyborem punktu węzłowego na po-
czątku, drugi raz na środku podprzedziału), oraz
trapezów, a następnie porównać wyniki obliczeń z
wartościami dokładnymi dla następujących zesta-
wów funkcji i przedziałów:
1. f (x) = x
2
, x ∈< 0, 1 >
2. f (x) = sin(x), x ∈< 0, π >
3. f (x) = e
x
, x ∈< 1, 2 >
4. f (x) = x, x ∈< 0, 1 >
Dla jakich przypadków która z kwadratur jest do-
kładna (dokładniejsza) i dlaczego? Proszę porównać
wyniki obliczeń z teoretycznymi oszacowaniami błę-
dów.
6.2
Zadanie 2
Dana jest funkcja f (x) = 0, 001x
2
e
x
. Zdefiniujmy
całkę oznaczoną postaci I(f ) =
R
4
1
f (x)dx. Proszę
obliczyć wartość całki za pomocą złożonej kwadra-
tury trapezów stosując minimalną liczbę potrzeb-
nych przedziałów tak, aby:
1. wykorzystując teoretyczne oszacowanie błędu
mieć gwarancję, że błąd kwadratury |R(f )| ¬
= 0, 01
2. Błąd kwadratury spełniał rzeczywiście zależ-
ność |R(f )| ¬ = 0, 01 (na podstawie porów-
nania z wartością rzeczywistą wyznaczoną ana-
litycznie, nie na podstawie oszacowań!)
3. Spełniony został warunek stopu w postaci
|S
m
(f ) − S
m+1
(f )| ¬ = 0, 01, gdzie S
m
(f )
oznacza wynik kwadratury dla podziału prze-
działu całkowania na m równych części
Proszę porównać uzyskaną liczbę przedziałów dla
badanych trzech metod definiowania kryterium
stopu.
4