dr hab. Piotr Fronczak
Metody numeryczne
Wykład nr 7
Całkowanie numeryczne
Całkowanie numeryczne to przybliżone obliczanie całek oznaczonych.
Metody całkowania numerycznego polegają na przybliżeniu całki za pomocą odpowiedniej sumy ważonej wartości całkowanej funkcji w kilku punktach. Aby uzyskać dokładniejsze przybliżenie dzieli się przedział całkowania na
niewielkie fragmenty. Ostateczny wynik jest sumą oszacowań całek w poszczególnych podprzedziałach.
I. Metoda prostokątów
b
a
b
a
a b a f dx a f dx
x f f
I( ) ( ) ( ) ( )( )
b
a
b
a
a b b f dx b f dx
x f f
I( ) ( ) ( ) ( )( )
Podzielmy przedział
[a,b]
naN
podprzedziałów.) )(
( ...
) )(
( )
)(
( )
( )
(
1 2 1 2 3 2 N N 1 Nb
a
x x
x f x
x x f x
x x f dx x f f
I
Ni
i i
i
x x
x f f
I
1
1
)
)(
( )
(
Gdy podprzedziały są równe
Ni
x
if h f
I
1
) ( )
(
Pewnym udoskonaleniem metody kwadratów jest wybór wartości funkcji w środku przedziału.
ba
a b b
f a b dx
f a f
I ( )
2 ) 2
(
N
i
i
i
x
f x h
f I
1
1
) 2 (
W przypadku wielu podprzedziałów:
II. Metoda trapezów
x0 x1 x
f(x) L(x)
Przybliżamy funkcję podcałkową prostą.
0 1
0 1
0 1 1 0
0 1
( ) ( ) ( )
a x , b x , , d dx ; h
0
( ) (1 ) ( ) ( ) ( )
1
x x x x
L x f x f x
x x x x
let x a h b a
b a x a
L f a f b
x b
0 1
0 1
0 1 1 0
0 1
( ) ( ) ( )
a x , b x , , d dx ; h
0
( ) (1 ) ( ) ( ) ( )
1
x x x x
L x f x f x
x x x x
let x a h b a
b a x a
L f a f b
x b
0 1
0 1
0 1 1 0
0 1
( ) ( ) ( )
a x , b x , , d dx ; h
0
( ) (1 ) ( ) ( ) ( )
1
x x x x
L x f x f x
x x x x
let x a h b a
b a x a
L f a f b
x b
1 0
1 1
0 0
1 1
2 2
0 0
( ) ( ) ( )
( ) (1 ) ( )
( ) ( ) ( ) ( ) ( )
2 2 2
b b
a
f x dx
aL x dx h L d
f a h d f b h d
f a h f b h h f a f b
W przypadku wielu podprzedziałów:
Ni
i
i
f x
x h f
f I
1
1
) (
) 2 (
)
(
II. Metoda Simpsona 1/3
x0 x1 x
f(x)
x2
h h
L(x)
Przybliżamy funkcję podcałkową parabolą.
1 x
x
0 x
x
1 x
x
h d dx
h ,
x , x
2 a h b
2 b x a
, b x
, a x
let
) x ( ) f x x
)(
x x
(
) x x )(
x x (
) x ( ) f x x
)(
x x
(
) x x )(
x x ) (
x ( ) f x x
)(
x x
(
) x x )(
x x ) (
x ( L
2 1 0
1 1 2
0
2 1
2 0
2
1 0
1 2
1 0
1
2 0
0 2
0 1
0
2 1
1 x
x
0 x
x
1 x
x
h d dx
h ,
x , x
2 a h b
2 b x a
, b x
, a x
let
) x ( ) f x x
)(
x x
(
) x x )(
x x (
) x ( ) f x x
)(
x x
(
) x x )(
x x ) (
x ( ) f x x
)(
x x
(
) x x )(
x x ) (
x ( L
2 1 0
1 1 2
0
2 1
2 0
2
1 0
1 2
1 0
1
2 0
0 2
0 1
0
2 1
1 x
x
0 x
x
1 x
x
h d dx
h ,
x , x
2 a h b
2 b x a
, b x
, a x
let
) x ( ) f x x
)(
x x
(
) x x )(
x x (
) x ( ) f x x
)(
x x
(
) x x )(
x x ) (
x ( ) f x x
)(
x x
(
) x x )(
x x ) (
x ( L
2 1 0
1 1 2
0
2 1
2 0
2
1 0
1 2
1 0
1
2 0
0 2
0 1
0
2 1
1 x
x
0 x
x
1 x
x
h d dx
h ,
x , x
2 a h b
2 b x a
, b x
, a x
let
) x ( ) f x x
)(
x x
(
) x x )(
x x (
) x ( ) f x x
)(
x x
(
) x x )(
x x ) (
x ( ) f x x
)(
x x
(
) x x )(
x x ) (
x ( L
2 1 0
1 1 2
0
2 1
2 0
2
1 0
1 2
1 0
1
2 0
0 2
0 1
0
2 1
) x ( 2 f
) 1 ) (
x ( f ) 1
( ) x ( 2 f
) 1 ) (
(
L
0
2 1
21
1 2
3 2
1
1 3
1 1
1 2
3 0
1 2 1
1 0
2 1
1 0 1
1 1 b
a
2 ) ξ 3
( ξ 2 ) h f(x
3 ) (ξ ξ
)h f(x
2 ) ξ 3
( ξ 2 ) h f(x
)dξ ξ(ξ 1
2 ) h f(x ξ )dξ
1 ( )h f(x
)dξ ξ(ξ 1
2 ) h f(x dξ
) ( L h
f(x)dx
) x ( 2 f
) 1 ) (
x ( f ) 1
( ) x ( 2 f
) 1 ) (
(
L
0
2 1
2 f(x ) 4f(x ) f(x )
3
f(x)dx h
0 1 2b
a
W przypadku
N
podprzedziałów o szerokościh = (b – a) / N
:UWAGA: Ponieważ wielomian kwadratowy jest określony poprzez trzy punkty (a zatem dwa przedziały), to liczba przedziałów musi być parzysta.
( ) 4 ( ) ( )
) 3 ( ),
( )
(
1 16 , 4 , 2
N i i i ii
i
h f x f x f x
f I gdzie
f I f
I
) ( )
( 2
) ( 4
) 3 (
) (
1
7 , 5 , 3 6
, 4 , 2
b f x
f x
f a
h f f
I
N
j
j N
i
i
Jest to ważona suma wartości funkcji w punktach definiujących podprzedziały.
II. Metoda Simpsona 3/8
x0 x1 x
f(x)
x2
h h
L(x)
x3 h
Przybliżamy funkcję podcałkową wielomianem trzeciego stopnia.
) x ( ) f x x
)(
x x
)(
x x
(
) x x )(
x x )(
x x (
) x ( ) f x x
)(
x x
)(
x x
(
) x x )(
x x )(
x x (
) x ( ) f x x
)(
x x
)(
x x
(
) x x )(
x x )(
x x (
) x ( ) f x x
)(
x x
)(
x x
(
) x x )(
x x )(
x x ) (
x ( L
3 2
3 1
3 0
3
2 1
0
2 3
2 1
2 0
2
3 1
0
1 3
1 2
1 0
1
3 2
0
0 3
0 2
0 1
0
3 2
1
f ( x ) 3 f ( x ) 3 f ( x ) f ( x )
8 h 3
3 a - h b
; L(x)dx f(x)dx
3 2
1 0
b a b
a
W przypadku
N
podprzedziałów o szerokościh = (b – a) / N
:UWAGA: Ponieważ wielomian kubiczny jest określony poprzez cztery punkty (a zatem trzy przedziały), to liczba podprzedziałów musi być podzielna przez 3.
( ) (
) 2 ( ) ( )
3 ) 8 (
) 3 (
2
10 , 7 , 4 1
8 , 5 , 2
1
f x f b
x f x
f a
h f f
I
N
j
j N
i
i i
Przykład: I
4xe dx
0
x
2 (= 5216.93)
metoda prostokątów
metoda trapezów
metoda Simpsona 1/3
metoda Simpsona 3/8
N = 1 436 23847 15898 6819
-91.64% 357.19% 204.79% 30.73%
N = 6 4749 6179 5331 5430
-8.95% 18.46% 2.20% 4.10%
N = 12 5094 5468 5225 5235
-2.34% 4.83% 0.17% 0.36%
N = 18 5162 5326 5218 5220
-1.04% 2.11% 0.04% 0.08%
Oszacowanie błędów
Metoda prostokątów
b
a
b
a
a b a f dx a f dx
x f f
I( ) ( ) ( ) ( )( )
) )(
( ' ) ( )
(x f a f x a
f
'( )( )22 ) 1 )(
( )
)(
( ' ) ( )
(x dx f a f x a dx f a b a f b a
f
b
a b
a
b
a
a b a f dx x f
E ( ) ( )( )
Błąd wynosi
Rozwijając
f(x)
w szereg Taylora blisko punktua
otrzymujemygdzie jest punktem między
a
ib
.Całkując obie strony tego równania otrzymujemy
)2
)(
( 2 '
1 f b a
E
Zatem błąd wyniesie
Zależy on od szerokości przedziału całkowania i wartości pierwszych pochodnych
f(x)
w tym przedziale.Błąd można znacząco zredukować dzieląc przedział na mniejsze podprzedziały.
)2
)(
( 2 '
1 f b a
E
Ponieważ
) 2
( 2 '
1 f h
Ei i to
N
i
i i
i h f
E E
1
2 '( )
2
1
Całkowity błąd
N f f
N
i
i 1
) ( ' '
Przyjmując, że średnio
2 ' ) (b a h f E
N a h (b ) oraz
) (h
O
)
2)(
2 ( ) 1 )(
( )
( )
( a f x f x a x f a x
f
ba
a b f
a b b f dx x f a
b a
f ( )( )
36 ) 1 )(
( )
( 2
) )(
(
b
a
a b f
a b b f a
f dx
x
f ( )( )
312 ) 1
))(
( )
( 2 ( ) 1
(
Metoda trapezów
Rozwijając w szereg Taylora
b
a
b
a
a b f
dx x f x a
af b
af dx
x f a
b a
f ( )( )
36 ) 1
( )
( )
( )
( )
)(
(
Całkując
ba b
a b
a
dx x f x
xf dx
x f
x ( ) ( ) ( )
Pamiętając, że
N f f
N
i
i 1
) ( '' ''
Przyjmując, że średnio
N a h (b ) oraz
b a
a b f
a b b f a f dx
x
f ( )( )3
12 ) 1 ))(
( )
( 2( ) 1
(
) 24 (
)
(
2 2h O h
a f
E b
Metoda prostokątów (udoskonalona)
) 180 (
)
(
4 4h O h
a f
E b
IV
Metoda Simpsona 1/3
) 80 (
)
(
4 4h O h
a f
E b
IV
Metoda Simpsona 3/8
) 12 (
)
(
2 2h O h
a f
E b
Otrzymujemy
Ekstrapolacja Richardsona
...
'' '
)
( 1 2
I h Khk K hk K hk I
) ( )
( 1
I h Khk O hk I
) ( )
2 / ( ) 2 /
( 1
I h K h k O hk I
Idea: zredukować błąd metody numerycznej z
O(h
k)
doO(h
k+1)
Rozwińmy całkę jako funkcję
h
w szereg wokół punktuh=0
.K,K’,K’’
niewiadome – reprezentują wartości błędówI
– rzeczywista wartość całkiI(h)
– numeryczna wartość całkiRedukując
h
otrzymujemy inne równanie naI
:Pamiętajmy, że
O(h
k+1)
w obu równaniach jest inne.(2) oraz
(1)
) 2 (
) 2 / ( ) ( )
2 / ( ) 2 / (
) ( )
(
1 1
1
k k
k k
k k k
h Kh O
h I h
O h
K h
I I
h O Kh h
I I
) ( ) ( ) 2 / ( 2
) ( )
( ) ( )
2 / ( 2 1
2
2
2
1
1 1
k k
k k
k k
k k
k k
h O h I h
I
h O Kh h
I h
O Kh h
I )I
- (
I
I (2) (1)
Zatem mamy dwa równania na
I
Wyeliminujmy wiodący błąd
K
:) 1 (
2
) ( ) 2 / (
2 1
kI h k I h O hk I
Czyli wyeliminowaliśmy błąd
Kh
k i zredukowaliśmy rząd błędu zk
dok+1
.Przykład: metoda trapezów
2
12 ) ) (
( b a f h
h I
I
2 1 1
)
( h Kh I
I
2 2 2
)
( h Kh I
I
K
jest takie same (średnia wartość drugiej pochodnej pof
nie zależy odh
).2
2 1
2 2
2 1 1
1
) ( )
(
h h
h h I
h h I I
Eliminując
K
otrzymujemy:Dla
h
1= 2h
23
) ( )
(
4 I h
2I h
1I
=> Błąd
O(h
4)
...
)
(
1 2
2 4
3 6
I h K h K h K h
Można pokazać, że
I
4 1 1
)
( h Kh I
I
4 2 2
)
( h Kh I
I
4
2 1
2 4
2 1 1
1
) ( )
(
h h
h h I
h h I I
Eliminując
K
otrzymujemy:Dla
h
1= 2h
215
) ( )
(
16 I h
2I h
1I
BłądO(h
6)
Gdy mamy całkę obliczoną numerycznie z dokładnością
O(h
4)
postępujemy analogicznie:Uogólnienie:
Jeśli
I
n jest przybliżeniem wartości całki przy użyciun
podprzedziałów (krokh
), aI
2n jest przybliżeniem wartości całki przy użyciu2n
podprzedziałów (krokh/2
) Jeśli oba przybliżenia mają ten sam błąd rzęduh
p, to nowe przybliżenie całki1 2
2
2
pI
pn I
nI
będzie miało błąd rzędu
h
p+2.3, 2,
1, k
;
1 4
4
1, ,1
, k
k j k
j k k
j
I I I
255 256
63 64
15 16
3 4
16 /
8 /
4 /
2 /
) ( )
( )
( )
( )
(
4 3
2 1
0
3 , 3
, 1 2
, 2
, 1 1
, 1
, 1 0
, 0
, 1 0
, 4
1 , 3 0
, 3
2 , 2 1
, 2 0
, 2
3 , 1 2
, 1 1
, 1 0
, 1
4 , 0 3
, 0 2
, 0 1
, 0 0
, 0
10 8
6 4
2
j j
j j
j j
j
j
I I I I I I I
I I
h
I I
h
I I
I h
I I
I I
h
I I
I I
I h
h O h
O h
O h
O h
O
k k
k k
k
Metoda Romberga
Metoda ta zwiększa dokładność oszacowania całki poprzez sukcesywne zastosowanie ekstrapolacji Richardsona.
926477 .
5216 dx
xe I
40
x
2
% 00050 .
0
% 00168 .
0
% 0053 .
0
% 0527 .
0
% 66 . 2
95 . 5355 25
. 0
68 . 5219 76
. 5764 5
. 0
20 . 5217 75
. 5256 79
. 7288 1
01 . 5217 14
. 5229 98
. 5670 2
. 12142 2
95 . 5216 84
. 5224 68
. 5499 41
. 8240 7
. 23847 4
) ( )
( )
( )
( )
(
4 3
2 1
0 .
10 8
6 4
2
hh h h h
h O h
O h
O h
O h
O
k k
k k
k
trapezów Met
Przykład:
Błąd przy całkowaniu podstawową metodą trapezów.
By uzyskać błąd rzędu 10-8 należy
przyjąć krok 10-8 (czyli 108 podprzedziałów)
Kolejne przybliżenia metodą Romberga
O(h
2), O(h
4) i O(h
6)
.Ten sam błąd można uzyskać przy
znacznie mniejszej liczbie podprzedziałów.
Kwadratura Gaussa na przedziale [-1, 1]
Dotychczas wszystkie metody wykorzystywały do obliczenia całki ważone wartości równo rozmieszczonych punktów. Położenia punktów były stałe dla danej metody.
W metodzie zwanej kwadraturą Gaussa nie są rozłożone równomiernie (ponadto punkty krańcowe również nie są brane pod uwagę). Położenie punktów oraz
odpowiadające im wagi są dobierane tak, by zminimalizować błąd.
Wybieramy (
c
1, c
2, x
1, x
2) tak, by metoda zwróciła dokładną wartość całki, gdy funkcja ma postaćf(x) = x
0, x
1, x
2, x
3) ( )
( )
( )
( )
(
1 1 2 21
1 1
n n
i n
i
i
f x c f x c f x c f x
c dx
x
f
) f(x c
) f(x c
f(x)dx :
n
2 2
1 1
1
2
1
x
2x
1-1 1
Ogólna postać kwadratury Gaussa
Dokładna całka dla
f = x
0, x
1, x
2, x
3– Cztery równania z czterema niewadomymi
) f(x c
) f(x c
f(x)dx :
n
1 1 1 2 22
1
3 1
3 1 1 1
0 3 2 0 2 1
1
2 1 2 1
3 2 2 3
1 1
1 1
3 3
2 2 2 2
1 1
1 1
2 2
2 2 1
1
1 1
2 1
1 1
x x c c
x c x
c dx
x x
f
x c x
c dx
x x
f
x c x
c xdx
x f
c c
dx f
3 ) ( 1 3 )
( 1 )
1
(
1
f x dx f f
I
-1 0 1 0.5
1.0
f(x) = x
23
1
3 1
Kwadratura Gaussa daje dokładną wartość,
gdy funkcja
f(x) = 1, x, x
2, x
3, lub liniowa kombinacja powyższych.-1 0 1
0.5 1.0
3
1
3 1
4 6829 .
1 )
cos(
1
1
dx x
676 . 3 1
cos 1 3
cos 1
Błąd 4.2%) ( )
( )
( )
( :
3
1 1 1 2 2 3 31
f x dx c f x c f x c f x
n
x
3x
1-1 x
21
Aby zwiększyć dokładność obliczeń należy wziąć więcej wyrazów w kwadraturze:
Wybieramy (
