Całkowanie
Całka jest jednym z najwaŜniejszych pojęć współczesnej analizy matematycznej. Stosuje się ją w
matematyce, fizyce, technice i wielu innych dziedzinach nauki. W matematyce badaniem własności i
obliczaniem wartości całek zajmuje się dział zwany rachunkiem całkowym. Za twórców tego rachunku uwaŜa
się Newtona oraz Leibniza, którzy opracowali teorię i metody związane z pojęciem całki i wprowadzili
terminologię i oznaczenia zbliŜone do stosowanych współcześnie, ukazujące związek rachunku całkowego z
rachunkiem róŜniczkowym.
1.
Całka nieoznaczona
Całką funkcji f(x) będzie kaŜda funkcja F(x), której pochodna jest równa f(x), czyli
F '(x) = f(x)
(1)
Funkcję F(x) nazywamy funkcją pierwotną danej funkcji f(x). Funkcja pierwotna moŜe być znaleziona z
dokładnością do stałej, poniewaŜ stała znika w trakcie wyznaczania pochodnej. Funkcje pierwotne danej
funkcji f(x) tworzą klasę funkcji róŜniących się stałymi:
F(x) + C
(2)
Przykład:
Znajdźmy funkcję pierwotną od f(x) = 2x. Korzystamy ze znanego wzoru na pochodną funkcji:
F '(x) = 2x
poniewaŜ (ax
n
)' = anx
n-1
zatem F '(x) = (x
2
+ C)' = 2x
stąd F(x) = x
2
+ C
Funkcję pierwotną F(x) nazywamy całką nieoznaczoną funkcji f(x) i przedstawiamy:
( )
( )
C
x
F
dx
x
f
+
=
∫
(3)
Zagadnienie obliczania całek nieoznaczonych sprowadza się do znajdowania postaci funkcji pierwotnej co jest
dosyć skomplikowane. W całkowaniu numerycznym głównie zajmujemy się całkami oznaczonymi.
2.
Całka oznaczona
Całka oznaczona funkcji f(x) w przedziale <x
p
,x
k
>
jest liczbą przedstawianą symbolem:
( )
dx
x
f
k
p
x
x
∫
(4)
i obliczaną na róŜne sposoby. Przedział <x
p
,x
k
>
nazywa się przedziałem całkowania. Krańce tego przedziału
nazywa się granicami całkowania. x
p
jest dolną granicą całkowania, a x
k
górną granicą.
Całka oznaczona Newtona-Leibniza
Jeśli znamy przepis na funkcję pierwotną F(x) funkcji podcałkowej f(x), to wartość całki oznaczonej funkcji f(x)
w przedziale <xp,xk> moŜemy obliczyć w prosty sposób wg definicji Newtona-Leibniza:
( )
( )
( )
p
k
x
x
x
F
x
F
dx
x
f
k
p
−
=
∫
(5)
Z punktu widzenia obliczeń komputerowych sposób ten nie jest najlepszy, poniewaŜ wymaga wyznaczenia
funkcji pierwotnej F(x), co jest zadaniem numerycznie skomplikowanym - dopiero w ostatnich latach
rozwiązano problem numerycznego wyznaczania wzoru całki dowolnej, całkowalnej funkcji w postaci
analitycznej. Z metody tej korzystają zaawansowane pakiety matematyczne w stylu MathCAD.
Całka oznaczona Riemanna
Większe zastosowanie w obliczeniach numerycznych ma definicja Riemanna, w której całka oznaczona jest
interpretowana jako suma pól obszarów ograniczonych wykresem funkcji f(x) oraz osią OX. Obszary leŜące
pod osią mają w tej interpretacji pola ujemne.
Przedział całkowania <xp,xk> dzielimy na rozłączne podprzedziały wyznaczając n+1 punktów podziałowych:
x
p
= x
o
< x
1
< ... < x
n-1
< x
n
= x
k
Punkty podziałowe muszą być tak dobrane, aby przy wzroście n do nieskończoności maksymalna odległość
między sąsiednimi punktami malała do zera, czyli:
dla i = 1, 2, ..., n
(
)
0
max
lim
1
=
−
−
∞
→
i
i
n
x
x
Pomiędzy parami sąsiednich punktów x
i-1
i x
i
wybieramy dowolne punkty t
i
spełniające nierówność:
dla i = 1, 2, ..., n
x
i-1
< t
i
< x
i
Całka oznaczona Riemanna jest wtedy granicą sum n prostokątów o podstawie równej (x
i
- x
i-1
) i wysokości
f
(t
i
), czyli
( )
( )(
)
1
1
lim
−
=
∞
→
−
=
∑
∫
i
i
n
i
i
n
x
x
x
x
t
f
dx
x
f
k
p
(6)
Jeśli f(t
i
) jest mniejsze od zera, to pole tego prostokąta zostanie zsumowane ze znakiem minus. Gdy odległości
pomiędzy punktami podziałowymi zbliŜają się do zera, suma pól prostokątów dąŜy do pola obszaru
ograniczonego wykresem funkcji.
RóŜnicę pomiędzy wartością pola pod wykresem funkcji a polem otrzymanym jako suma skończonej ilości
prostokątów nazywamy błędem całkowania. Wraz ze wzrostem liczby prostokątów w całce Riemanna błąd
całkowania dąŜy do zera.
3.
Metody numeryczne całkowania
Metoda prostokątów
W metodzie prostokątów korzystamy z definicji całki oznaczonej Riemanna, w której wartość całki
interpretowana jest jako suma pól obszarów pod wykresem krzywej w zadanym przedziale całkowania <x
p
,x
k
>
.
Sumę tę przybliŜamy przy pomocy sumy pól odpowiednio dobranych prostokątów. Sposób postępowania jest
następujący:
Przedział całkowania <x
p
,x
k
>
dzielimy na n równo odległych punktów x
1
,x
2
,...,x
n
. Punkty te wyznaczamy w
prosty sposób wg wzoru:
dla i = 1,2,...,n
(
)
p
k
p
i
x
x
n
i
x
x
−
+
=
(7)
Obliczamy odległość między dwoma sąsiednimi punktami - będzie to podstawa kaŜdego prostokąta:
n
x
x
dx
p
k
−
=
(8)
Dla kaŜdego wyznaczonego w ten sposób punktu obliczamy wartość funkcji f(x) w tym punkcie:
f
i
= f(x
i
),
dla i = 1,2,...,n
(9)
Obliczamy sumę iloczynów wyznaczonych wartości funkcji przez odległość dx między dwoma sąsiednimi
punktami - da to sumę pól poszczególnych prostokątów ograniczonych wykresem funkcji:
S = f
1
dx + f
2
dx + ... + f
n
dx
(10)
a po wyprowadzeniu wspólnego czynnika przed nawias:
S = dx (f
1
+ f
2
+ ... + f
n
)
(11)
Otrzymana suma jest przybliŜoną wartością całki oznaczonej funkcji f(x) w przedziale <x
p
,x
k
>
.
( )
∑
∫
=
−
+
−
≈
n
i
p
k
p
x
x
p
k
n
x
x
i
x
f
n
x
x
dx
x
f
k
p
1
(12)
Przykład:
Obliczymy ręcznie przybliŜoną wartość całki oznaczonej z funkcji f(x) = sin(x)
w przedziale <0,π>.
Przedział podzielimy na n = 4 punkty:
Odległość między dwoma sąsiednimi punktami wynosi:
Dla kaŜdego z wyznaczonych punktów obliczamy wartość funkcji f(x) = sin(x):
Obliczamy sumę pól prostokątów:
S
= dx (f
1
+ f
2
+ f
3
+ f
4
)
S
= 0,7854 (0,7071 + 1,000 + 0,7071 + 0,0000)
S
= 0,7854 × 2,4142
S
= 1,8961
Dokładna wartość takiej całki oznaczonej wynosi wg tablic:
∫
=
π
0
2
sin xdx
Zatem popełniliśmy błąd równy: 2 - 1,8961 = 0,1039, co stanowi nieco ponad 5% wartości dokładnej. W sumie
nie jest to zły wynik zwaŜywszy na ilość wykonanych przez nas rachunków. Jeśli chcemy zwiększyć
dokładność, to musimy zsumować więcej prostokątów, ale to zostawimy juŜ komputerom.
Metoda trapezów
Metoda prostokątów nie jest zbyt dokładna, poniewaŜ pola uŜytych w niej prostokątów źle
odwzorowują pole pod krzywą. DuŜo lepszym rozwiązaniem jest zastosowanie zamiast nich trapezów o
wysokości dx i podstawach równych odpowiednio wartości funkcji w punktach krańcowych.. Sama zasada nie
zmienia się.
Przedział całkowania <x
p
,x
k
>
dzielimy na n równo odległych punktów x
1
,x
2
,...,x
n
. Punkty te wyznaczamy w
prosty sposób wg wzoru:
dla i = 0,1,2,...,n
(
)
p
k
p
i
x
x
n
i
x
x
−
+
=
(13)
Obliczamy odległość między dwoma sąsiednimi punktami - będzie to wysokość kaŜdego trapezu:
n
x
x
dx
p
k
−
=
(14)
Dla kaŜdego wyznaczonego w ten sposób punktu obliczamy wartość funkcji f(x) w tym punkcie:
f
i
= f(x
i
),
dla i = 0,1,2,...,n
(15)
Pole pod wykresem funkcji przybliŜane jest polami n trapezów. Pole i-tego trapezu obliczamy wg wzoru:
dla i=1,2,...,n
dx
f
f
P
i
i
i
2
1
+
=
−
(16)
PrzybliŜona wartość całki jest sumą pól wszystkich otrzymanych w ten sposób trapezów:
s = P
1
+ P
2
+ ... + P
n
(17)
czyli
(
)
(
)
n
n
n
n
n
n
n
n
n
n
f
f
f
f
f
dx
s
f
f
f
f
f
f
f
f
f
f
dx
s
dx
f
f
dx
f
f
dx
f
f
dx
f
f
dx
f
f
s
+
+
+
+
+
=
+
+
+
+
+
+
+
+
+
+
=
+
+
+
+
+
+
+
+
+
+
=
−
−
−
−
−
−
−
1
2
1
0
1
1
2
3
2
2
1
1
0
1
1
2
3
2
2
1
1
0
2
2
2
2
2
2
2
2
2
2
K
K
K
(18)
+
+
+
+
+
=
−
2
0
1
2
1
n
n
f
f
f
f
f
dx
s
K
(19)
Wyprowadzony na końcu wzór jest podstawą przybliŜonego wyliczania całki w metodzie trapezów.
( )
( )
( )
+
+
−
+
−
≈
∑
∫
−
=
2
1
1
k
p
n
i
p
k
p
p
k
x
x
x
f
x
f
n
x
x
i
x
f
n
x
x
dx
x
f
k
p
(20)
Metoda Simpsona
Metoda ta jest dokładniejsza od opisanych tutaj metod przybliŜonego całkowania. W metodzie
prostokątów całka oznaczona przybliŜana była funkcjami stałymi - liczyliśmy sumę pól prostokątów. W
metodzie trapezów całkę przybliŜaliśmy za pomocą funkcji liniowych - obliczaliśmy sumy pól trapezów. W
metodzie Simpsona stosujemy jako przybliŜenie parabolę - będziemy obliczali sumy wycinków obszarów pod
parabolą.
Zasada jest następująca:
Przedział całkowania <x
p
,x
k
>
dzielimy na n + 1 równo odległych punktów x
o
, x
1
, x
2
,..., x
n
:
dla i = 0,1,2,...,n
(
)
p
k
p
i
x
x
n
i
x
x
−
+
=
(21)
Dla kaŜdych dwóch sąsiednich punktów wyznaczamy punkt środkowy t
i
wg wzoru:
dla i = 1,2,...,n
2
1
i
i
i
x
x
t
+
=
−
(22)
Obliczamy odległość między dwoma sąsiednimi punktami.
n
x
x
dx
p
k
−
=
(22)
Dla kaŜdego wyznaczonego w ten sposób punktu obliczamy wartość funkcji f(x) w tym punkcie:
punkty podziałowe:
dla i = 0,1,2,...,n
f
i
= f(x
i
)
(23)
punkty środkowe:
dla i = 1,2,...,n
f
ti
= f(t
i
)
(24)
W kaŜdym podprzedziale <x
i-1
,x
i
>
przybliŜamy funkcję za pomocą paraboli g(x) o następującej postaci:
dla i = 1,2,...,n
( )
i
i
i
i
i
i
x
x
x
c
x
b
x
a
x
g
,
1
2
−
∈
+
+
=
(25)
Parabola g
i
(x)
musi przechodzić przez punkty: (x
i-1
,f
i-1
), (t
i
,f
ti
), (x
i
,f
i
). Współczynniki a
i
, b
i
i c
i
wyznaczymy
zatem z układu trzech równań:
dla i = 1,2,...,n
i
i
i
i
i
i
ti
i
i
i
i
i
i
i
i
i
i
i
f
c
x
b
x
a
f
c
t
b
t
a
f
c
x
b
x
a
=
+
+
=
+
+
=
+
+
−
−
−
2
2
1
1
2
1
(26)
Po przekształceniach otrzymujemy:
dla i = 1,2,...,n
( )
(
)
(
)
ti
i
i
x
x
p
k
i
f
f
f
n
x
x
dx
x
g
i
i
4
6
1
1
+
+
−
=
−
∫
−
(27)
Wzór ten pozwala wyliczyć pole obszaru pod parabolą aproksymującą funkcję f(x) w przedziale <x
i-1
,x
i
>.
Wartość całej całki otrzymamy sumując te pola, czyli:
( )
(
)
+
+
−
=
∑
∑
∑
∫
=
=
=
−
−
n
i
i
n
i
ti
n
i
i
x
x
p
k
f
f
f
n
x
x
dx
x
f
i
i
1
1
1
1
4
6
1
(28)
4.
Całkowanie w matlab / octave
Problem polega na znalezieniu całki oznaczonej funkcji f(x) na przedziale <a,b>. Jest to łatwe, gdy znana jest
funkcja pierwotna F(x) taka, Ŝe F(x)’=f(x) . Nie zawsze jest to jednak moŜliwe. Metody całkowania
numerycznego (kwadratury) polegają na przybliŜeniu funkcji podcałkowej f na danym przedziale <a,b> lub
jego podprzedziałach przy pomocy innej funkcji, dla której wartość całki jest określona analitycznie. Matlab
stosuje kwadratury Newtona-Cotesa - quad (interpolacja wielomianem drugiego stopnia) i Simpsona -
quadl
(interpolacja wielomianowa – dobierany stopień wielomianu).
Q=quad(f,a,b,tol,trace);
Q=quadl(f,a,b,tol,trace);
f
– łańcuch zawierający nazwę funkcji, funkcja musi być umieszczona w odpowiednim skrypcie, musi
zwracać wektor wartości a jej argumentem jest wektor elementów.
a,b
– przedział całkowania,
tol
– wymagana tolerancja względna, domyślnie 10^(-3)
trace
– parametr opcjonalny, pozwala na rysowanie wykresu z węzłami kwadratury.
Przykład:
function [y] = funkcja_calkowana(x)
%% funkcja
y=sin(x.*x);
end
Wywołanie funkcji całkowania
>> q1 = quad(‘funkcja_calkowana’,0,pi,1e-5,1);
>> ql = quadl(‘funkcja_calkowana’,0,pi,1e-5,1);
Funkcja trapz – całkowanie metodą trapezów
Przykład:
>> x=linspace(0,pi,100);
>> y=sin(x);
>> trapz(x,y)
5.
Ćwiczenia:
5.1.
Zdefiniuj funkcję y=x
2
+2x+1
w pliku f1.m
5.2.
Utwórz skrypt main5.m w którym określ:
•
początek i koniec przedziału całkowania <-1,1> pod zmiennymi x
p
, x
k
•
ilość podprzedziałów całkowania jako 3 pod zmienną n
•
narysuj funkcję podcałkową w przedziale całkowania – kolor czerwony
5.3.
Uzupełnij podprogram f_rectI.m reprezentujący metodę prostokątów według wskazówek zawartych
wewnątrz pliku. Z plikami f_trapI (metoda trapezów) oraz f_SimpI (metoda Simpsona) naleŜy postąpić
analogicznie – wskazówki oraz odniesienia do wzorów określających daną metodę są komentarzami wewnątrz
plików.
5.4.
W programie main5.m wywołaj wymienione funkcje z parametrami w kolejności:
(funkcja całkowana, x
p
, x
k
, n)
5.5.
Porównaj wyniki wszystkich metod dla róŜnych n: np. 3, 7, 10
5.6.
Wykorzystaj funkcje wbudowane (quad, quadl, trapz) do obliczenia całki funkcji f1
5.7.
Porównaj wyniki, sprawdź wpływ parametru tol