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