Mariusz PYRZ
SIMR (PW), Instytut Pojazdów
Metody numeryczne w mechanice
Całkowanie numeryczne
8.
Mariusz PYRZ
SIMR (PW), Instytut Pojazdów
Metody numeryczne w mechanice
Całkowanie numeryczne
8.
Przybli
Ŝ
one obliczanie całek oznaczonych
Aby obliczyć wartość całki oznaczonej
gdzie f: R->R spełnia niezbędne warunki róŜniczkowania i całkowania
przybliŜamy funkcję f(x) w przedziale [a,b] za pomocą wielomianu P
N
(x)
łatwego do scałkowania i na jego podstawie obliczamy przybliŜoną wartość I
a
( )
b
b
I
f x dx
=
∫
MNM Całkowanie numeryczne M.Pyrz 04-2014
2
Aby uzyskać dokładniejsze przybliŜenie rozwiązania
(ale uniknąć zwiększania stopnia wielomianu interpolacyjnego P
N
(x) w [a,b]),
dzielimy przedział [a,b] na podprzedziały, w kaŜdym z nich obliczamy
oddzielnie wartość całki i następnie sumujemy wszystkie oszacowania.
( )
( )
b
b
a
N
b
b
I
P x dx
f x dx
=
≈
∫
∫
Obliczanie całki w przedziale [a,b]
Zastosowanie:
obliczanie powierzchni, objętości obszarów o znanym brzegu, charakterystyk przekrojów
(środek cięŜkości, moment bezwładności), itp.
3
Rys. 8.1
(
) /
h
b
a
N
= −
MNM Całkowanie numeryczne M.Pyrz 04-2014
Metoda prostok
ą
tów
Oszacowaniem wartości całki w przedziale [a,b] jest pole prostokąta
( )
( )(
)
b
a
a
f x dx
I
f a b
a
≈
=
−
∫
( )(
)
a
I
f b b
a
=
−
f(a) – lewa granica przedziału
f(b) - prawa granica przedziału
4
Rys. 8.2
(
)(
)
2
a
a
b
I
f
b
a
+
=
−
f((a+b)/2) – środek przedziału
MNM Całkowanie numeryczne M.Pyrz 04-2014
Numeryczne obliczanie całki - kwadratura
Dla skończonej liczby punktów (węzłów)
moŜna wyznaczyć wielomian interpolacyjny (przechodzący przez te węzły)
Wartość całki (dla wielomianu) równa jest sumie waŜonej wartości całkowanej
0
1
2
...
N
a
x
x
x
x
b
=
< <
< <
=
[ , ]
0,1,...,
k
x
a b
k
N
∈
=
0
( )
( ) ( )
N
N
i
i
P x
L x f x
=
∑
5
funkcji liczonej w kilku odpowiednich punktach (nazywanej takŜe kwadraturą)
A
k
są stałymi niezaleŜnymi od f(x):
( )
b
k
k
b
A
L x dx
=
∫
0
ˆ
( )
(
),
[ , ]
b
N
N
N
k
k
k
k
a
I
P x dx
A f x
x
a b
=
=
=
∈
∑
∫
MNM Całkowanie numeryczne M.Pyrz 04-2014
Dokładno
ść
w obliczaniu całek
Współczynniki A
k
i węzły kwadratury x
k
powinny być tak dobrane,
aby zminimalizować błąd aproksymacji
ˆ
N
N
R
I
I
=
−
Aby kwadratura była zbieŜna, musi zachodzić
0
N
N
R
→∞
→
6
MNM Całkowanie numeryczne M.Pyrz 04-2014
Obliczanie przybli
Ŝ
onej warto
ś
ci całek
Metody Newtona-Cotesa
wykorzystują interpolację funkcji f(x) za pomocą
wielomianów Lagranga L
N
, budowanych na węzłach równoodległych
Metoda trapezów (najprostsza z metod Newtona-Cotesa) dokonuje interpolacji funkcji f(x) za
0
0,
( )
(
)
N
N
j
N
n
n
j
j n
n
j
x
x
L
x
f x
x
x
=
=
≠
−
=
−
∑
∏
,
(
) /
0,1,...,
n
x
a
nh
h
b a
N
n
N
= +
= −
=
7
Metoda trapezów (najprostsza z metod Newtona-Cotesa) dokonuje interpolacji funkcji f(x) za
pomocą wielomianu pierwszego stopnia (tj. odcinka prostej).
Metody Gaussa
stosują interpolacje budowane na węzłach
ξ
i
rozłoŜonych
nierównomiernie i dobieranych tak, aby anulować błąd aproksymacji R
N
.
MNM Całkowanie numeryczne M.Pyrz 04-2014
1
1
1
( )
( )
n
i
i
i
I
f x dx
f
w
ξ
+
=
−
=
=
∑
∫
Metoda trapezów
Przedział całkowania [a,b] dzielony jest na N
podprzedziałów o długości h=(b-a)/N.
W kaŜdym z podprzedziałów krzywa f(x) jest
zastępowana odcinkiem prostej.
Otrzymana aproksymacja jest rzędu O(h
2
) co
oznacza, Ŝe róŜnica między dokładną wartością całki
a wartością obliczoną jest proporcjonalna do h
2
).
8
a wartością obliczoną jest proporcjonalna do h
2
).
f x dx
x
x
i
i
( )
+
z
1
[
]
1
ˆ
(
)
( )
2
i
i
i
h
I
f x
f x
−
=
+
MNM Całkowanie numeryczne M.Pyrz 04-2014
Algorytm
Ustalić h=(b-a)/N oraz x
i
=a+ih (mamy w sumie N+1 punktów x
i
, x
0
=a, x
N
=b ).
W podprzedziale [x
i
, x
i+1
] zastępujemy wartość całki
polem trapezu
Metoda trapezów
1
1
1
ˆ
( )
( )
( )
2
(
)
2
b
N
N
i
N
N
i
i
a
h
f x dx
I
R
f a
f b
f a
ih
R
−
=
=
=
+
=
+
+
+
+
∑
∑
∫
1
2
1
ˆ
ˆ
ˆ
ˆ
( )
...
b
N
N
i
i
a
f x dx
I
I
I
I
=
≈ + + +
=
∑
∫
9
Błąd
3
[ , ]
max
( )
12
N
x
a b
h
R
N
f
x
∈
′′
≤
MNM Całkowanie numeryczne M.Pyrz 04-2014
Metoda Simpsona
Przedział całkowania [a,b] dzielony jest na N
podprzedziałów o długości h=(b-a)/N.
W kaŜdym podprzedziale krzywa f(x) jest
zastępowana wielomianem interpolującym drugiego
stopnia (czyli łukiem paraboli).
Uzyskana aproksymacja jest rzędu O(h
4
) (czyli
róŜnica między wartością dokładną całki a wartością
obliczoną jest proporcjonalna do h
4
).
10
obliczoną jest proporcjonalna do h
4
).
2
2
2
( )
i
i
x
x
f x dx
+
∫
[
]
2
2
2
1
2
2
ˆ
(
)
4 (
)
(
)
3
i
i
i
i
h
I
f x
f x
f x
+
+
=
+
+
MNM Całkowanie numeryczne M.Pyrz 04-2014
Algorytm
Ustalić h=(b-a)/2N oraz x
i
=a+ih, (mamy teraz 2N+1 punktów x
i
, x
0
=a, x
2N
=b ).
W podprzedziale [x
2i
, x
2i+2
] zastępujemy wartość całki
polem „pod wielomianem”
Metoda Simpsona
1
/ 2 1
/ 2 1
2
0
1
0
ˆ
( )
( )
( )
2
(
2 )
4
(
(2
1) )
3
b
N
N
N
i
N
N
i
i
i
a
h
f x dx
I
R
f a
f b
f a
ih
f a
i
h
R
−
−
−
=
=
=
=
+
==
+
+
+
+
+
+
+
∑
∑
∑
∫
1
0
2
4
2(
1)
2
0
ˆ
ˆ
ˆ
ˆ
ˆ
( )
...
b
N
N
i
i
a
f x dx
I
I
I
I
I
−
−
=
≈ + + + +
=
∑
∫
11
5
(4)
[ , ]
max
( )
180
N
x
a b
h
R
N
f
x
∈
≤
MNM Całkowanie numeryczne M.Pyrz 04-2014
Błąd
Metoda Boole’a
Przedział [a,b] dzielony jest na podprzedziały. W kaŜdym z nich
funkcja f(x) jest zastępowana wielomianem 4 stopnia.
Uzyskana aproksymacja jest rzędu O(h
6
)
Algorytm
Ustalić h=(b-a)/4N, x
i
=a+ih.
W kaŜdym podprzedziale [x
4i
, x
4i+4
] wartość całki
f x dx
x
i
( )
4
4
+
z
12
W kaŜdym podprzedziale [x
4i
, x
4i+4
] wartość całki
jest przybliŜana wzorem
f x dx
x
i
( )
4
z
[
]
4
4
4
1
4
2
4
3
4
4
2
ˆ
7 (
) 32 (
) 12 (
) 32 (
) 7 (
)
45
i
i
i
i
i
i
h
I
f x
f x
f x
f x
f x
+
+
+
+
=
+
+
+
+
1
4
0
ˆ
( )
b
N
i
i
a
f x dx
I
−
=
≈
∑
∫
MNM Całkowanie numeryczne M.Pyrz 04-2014
Oszacowanie bł
ę
du
Oszacowanie wartości błędu całki w przedziale [a,b]:
Dla [a,b] podzielonego na N podprzedziałów:
(błąd liczymy rozwijając f(x) w szereg Taylora, całkujemy Ŝeby obliczyć I, odejmujemy I
a
,
wynik wyraŜamy w funkcji h – szerokości podprzedziału)
a
E
I
I
= −
1
N
i
i
E
E
=
=
∑
13
Metoda prostokątów
Metoda prostokątów (f w środku przedz.)
Metoda trapezów
Metoda Simpsona
2
( )
2
b
a
E
f h
O h
−
′
=
=
2
2
(
)
24
b
a
E
f h
O h
−
′′
=
=
2
2
(
)
12
b
a
E
f h
O h
−
′′
=
=
4
4
(
)
180
IV
b
a
E
f
h
O h
−
=
=
MNM Całkowanie numeryczne M.Pyrz 04-2014
Metoda Romberga
Metoda ta pozwala na ustalenie dokładności wyniku i prowadzenie obliczeń
iteracyjnych aŜ do osiągnięcia Ŝądanej dokładności. Algorytm metody
bazuje na oryginalnym wykorzystaniu metody trapezów.
•T
0
(h)
– oszacowanie całki uzyskane na podstawie metody trapezów o kroku h
•T
0
(h/2)
– dokładniejsze oszacowanie wartości całki I otrzymane dla kroku h/2
•T (h/4)
– oszacowanie całki I uzyskane przy kroku h/4
14
•T
0
(h/4)
– oszacowanie całki I uzyskane przy kroku h/4
•T
0
(h/8) …
•T
0
(h/2
n
)
- ciąg jest zbieŜny do dokładnej wartości I.
MoŜna wykazać, Ŝe moŜna utworzyć drugi ciąg T
1
(h/2
n-1
) w którym element
T
1
(h/2
k
) stanowi lepsze przybliŜenie całki niŜ T
0
(h/2
k
).
MNM Całkowanie numeryczne M.Pyrz 04-2014
Metoda Romberga
Otrzymuje się
W analogiczny sposób moŜna otrzymać trzeci ciąg
h
h
FG IJ FG IJ
T
h
T
h
T
h
k
n
k
k
k
1
0
1
0
2
4
2
2
4 1
0
1
F
HG
I
KJ
=
F
HG
I
KJ
−
F
HG
I
KJ
−
∈
−
+
[ ,
]
15
W sformułowaniu ogólnym
T
h
T
h
T
h
k
n
k
k
k
2
2
1
1
1
2
2
4
2
2
4
1
0
2
F
HG
I
KJ
=
F
HG
I
KJ
−
F
HG
I
KJ
−
∈
−
+
[ ,
]
T
h
T
h
T
h
k
n
p
p
k
p
p
k
p
k
p
2
4
2
2
4
1
0
1
1
1
F
HG
I
KJ
=
F
HG
I
KJ
−
F
HG
I
KJ
−
∈
−
−
+
−
[ ,
]
MNM Całkowanie numeryczne M.Pyrz 04-2014
Metoda Romberga
Grupujemy elementy w następującej tablicy
T h
T
h
T
h
T
h
T
h
T h
T
h
T
h
T
h
T h
T
h
T
h
n
n
n
n
n
n
0
0
0
2
0
1
0
1
1
1
2
1
1
2
2
2
2
2
2
2
2
2
2
2
2
2
b g
b g
b g
F
HG
I
KJ
F
HG
I
KJ
F
HG
I
KJ
F
HG
I
KJ
F
HG
I
KJ
F
HG
I
KJ
F
HG
I
KJ
F
HG
I
KJ
F
HG
I
KJ
−
−
−
−
−
⋯
⋯
⋯
16
MoŜna wykazać, Ŝe ciąg T
k
(h) zbiega się szybciej do I niŜ ciąg T
0
(h/2
k
).
T
h
T
h
T h
n
n
n
n
2
2
2
2
1
1
2
2
2
b g
b g
b g
HG KJ
HG KJ
F
HG
I
KJ
−
−
−
⋯
⋯
MNM Całkowanie numeryczne M.Pyrz 04-2014
Metoda Romberga
Algorytm
Obliczyć metodą trapezów T
0
(h), T
0
(h/2),
Na ich podstawie wyznaczyć T
1
(h) i następnie obliczyć T
0
(h/4),
Na podstawie znanych wartości wyznaczyć T
1
(h/2), T
2
(h), itd.
17
W ten sposób na podstawie tablicy rzędu n otrzymuje się tablicę rzędu (n+1)
obliczając T
0
(h/2
n+1
) i na tej podstawie wyznaczając elementy T
k
(h/2
n+1-k
).
Obliczenia kończy się gdy róŜnica między T
n+1
(h) i T
n
(h) staje się mniejsza
niŜ załoŜona wcześniej wartość.
MNM Całkowanie numeryczne M.Pyrz 04-2014
Metoda Romberga
( )
( )
0
0
0
0
0
2
3
4
1
1
1
1
2
3
2
2
2
2
2
2
2
h
h
h
h
T h
T
T
T
T
h
h
h
T h
T
T
T
⋯
⋯
18
Kryterium zatrzymania
( )
( )
1
n
n
T
h
T h
ε
+
−
≤
( )
( )
2
2
2
2
3
3
2
2
2
...
h
h
T
T
T h
h
T h
T
⋯
⋯
MNM Całkowanie numeryczne M.Pyrz 04-2014
Metody Gaussa
Całkowanie numeryczne
Idea: wybrać połoŜenie N węzłów x
k
(i odpowiadających im współczynników
wagowych A
k
) w taki sposób, aby anulować błąd R
N
aproksymacji Î
N
.
We wzorze na całkowanie występuje 2N niewiadomych (x
k
i A
k
). Zatem
moŜemy spróbować je określić dla wielomianu interpolacyjnego stopnia
1
ˆ
( )
(
)
b
N
N
N
k
k
N
k
a
I
f x dx
I
R
A f x
R
=
=
≈
+
=
+
∑
∫
19
moŜemy spróbować je określić dla wielomianu interpolacyjnego stopnia
k=2N-1 (w którym występuje 2N współczynników).
Kwadratury metody Gaussa wykorzystywane są do obliczania całek w
granicach [-1,1], (-
∞
, +
∞
), [0,
∞
) …
W dalszej części rozwaŜany będzie przypadek [a,b]=[-1,1] tj.
MNM Całkowanie numeryczne M.Pyrz 04-2014
1
1
( )
I
f x dx
+
−
=
∫
Metoda punktów Gaussa
Stosowana do całkowania wielomianów w granicach [-1,+1]
(duŜe zastosowanie w Metodzie Elementów Skończonych).
Całka obliczana jest jako suma wartości wielomianu wyznaczanych w
tzw. punktach całkowania, mnoŜonych przez współczynniki wagowe
1
1
( )
( )
n
i
i
i
I
f x dx
f
w
ξ
+
=
=
=
∑
∫
20
ξ
i
– połoŜenie i-ego punktu Gaussa,
w
i
– współczynnik (waga),
n – liczba punktów całkowania.
Stosując
n
punktów całkowania moŜna otrzymać dokładny wynik dla
wielomianu stopnia
2n-1
(lub niŜszego)
(na przykład do całkowania wielomianu stopnia 3,
potrzeba co najmniej 2 punktów Gaussa).
1
1
i
=
−
∑
∫
MNM Całkowanie numeryczne M.Pyrz 04-2014
Metoda punktów Gaussa
PołoŜenie punktów całkowania (punktów Gaussa)
- symetryczne rozmieszczenie względem środka przedziału całkowania
- współczynniki (wagi) punktów połoŜonych symetrycznie są jednakowe
PołoŜenie
ξ
i
i wagi w
i
punktów Gaussa
21
PołoŜenie
ξ
i
i wagi w
i
punktów Gaussa
Liczba punktów
PołoŜenie punktów
Współczynniki wagowe
1
ξ
1
=0.0
w
1
=2.0
2
ξ
1
=
ξ
2
=
±
0.57735
w
1
=w
2
=1.0
3
ξ
1
=
ξ
3
=
±
0.774596
ξ
2
=0.0
w
1
=w
3
=0.555556
w
2
=0.888889
4
ξ
1
=
ξ
4
=
±
0.86113
ξ
2
=
ξ
3
=
±
0.33998
w
1
=w
4
=0.34785
w
2
=w
3
=0.65214
MNM Całkowanie numeryczne M.Pyrz 04-2014
Całkowanie wielomianu - przykład
Przypadek n=2. Wielomian stopnia 3 do scałkowania
Całkowanie dokładne prowadzi do wyniku
JeŜeli zastosujemy 2 punkty całkowania
ξ
1
=+k,
ξ
2
=-k , będziemy
stosowali wzór
f
a
a
a
a
( )
ξ
ξ
ξ
ξ
= +
+
+
1
2
3
2
4
3
I
f
a
a
=
=
+
−
+
z
( )
ξ
2
2
3
1
3
1
1
22
wynika stąd błąd obliczeń
Błąd zniknie (dla dowolnych wartości a
1
i a
3
) jeŜeli spełnimy zaleŜności:
czyli
2
1
3
(
)
( )
2 (
)
a
I
wf
k
wf k
w a
a k
=
− +
=
+
2
1
3
1
2 (1
)
2
(
)
3
a
I
I
a
w
a
wk
−
=
−
+
−
1
0
1
3
0
2
− =
−
=
w
wk
w
k
=
= ±
= ±
1 0
1
3
0 57735
.
.
MNM Całkowanie numeryczne M.Pyrz 04-2014
Całki wielokrotne
Analogiczne sformułowanie
Dobór punktów w zaleŜności od stopnia
wielomianu względem poszczegόlnych
zmiennych
Carl Friedrich Gauss
(1777-1855)
1 1 1
p
n
m
+ + +
∑∑∑
∫ ∫ ∫
23
MNM Całkowanie numeryczne M.Pyrz 04-2014
(1777-1855)
1
1
1
1 1 1
( , , )
( ,
,
)
i
j
k
i
j
k
i
j
k
I
f x y z dxdydz
f
w w w
ξ η ζ
=
=
=
− − −
=
=
∑∑∑
∫ ∫ ∫