2013-11-23
METODY NUMERYCZNE
Wykład 5.
Całkowanie numeryczne
dr hab. inż. Katarzyna Zakrzewska, prof. AGH
1
Met.Numer. wykład 5
Plan
" Wzór trapezów
" Złożony wzór trapezów
" Metoda ekstrapolacji Richardsona
" Metoda Romberga
" Metoda Simpsona  wzór parabol
" Metoda Gaussa
2
Met.Numer. wykład 5
Całkowanie numeryczne - idea
b
Y
f(x)
I = f ( x )dx
+"
a
Całkę można
b
przybliżyć sumą f (x)dx
+"
a
n
S = f (ci )"xi
"
i=1
ab X
xi d" ci d" xi+1
3
Met.Numer. wykład 5
1
2013-11-23
Kwadratury Newtona-Cotesa
Kwadratura Newtona  Cotesa należy do metod z
ustalonymi węzłami, polega na tym, że funkcja f(x) jest
interpolowana wielomianem (np. Lagrange a)
f ( x ) H" fn( x )
gdzie:
fn( x ) = a0 + a1x + ...+ an-1xn-1 + anxn
Wówczas całka z funkcji f(x) może być przybliżana
całką z funkcji interpolującej fn(x) n-tego stopnia
b b
I = f (x) H" fn (x)
+" +"
a a
4
Met.Numer. wykład 5
Wzór trapezów
Zakładamy n = 1 czyli
f1(x) = a0 + a1x
b b b
I = f (x) H" f1(x) = (a0 + a1x)dx
+" +" +"
a a a
b2 - a2
= a0(b - a) + a1
2
Szukamy współczynników a0 i a1
Zakładamy, że prosta, która przybliża funkcję f(x) przechodzi
przez punkty (a,f(a)) i (b,f(b)). Czyli:
f (a)b - f (b)a
a0 =
f (a) = f1(a) = a0 + a1a
b - a
f (b) - f (a)
f (b) = f1(b) = a0 + a1b
a1 =
b - a
5
Met.Numer. wykład 5
Wzór trapezów
b
b
f (a) + f (b)
ń#
f (x)dx = (b - a)Ą# f (x)dx H" pole trapezu
+"
ó# Ą# +"
2
a Ł# Ś#
a
Y
f(x)
fI(x)
f(b)
f(a)
ab X
6
Met.Numer. wykład 5
2
2013-11-23
Wzór trapezów
Przykład 1:
Prędkość rakiety w przedziale czasu od t1=8 s do t2=30 s jest
opisana wzorem:
140000
ń#
v(t) = 2000lnĄ# - 9.8t
ó#140000 - 2100t Ą#
Ł# Ś#
a) Przy pomocy metody trapezów znajdz
przemieszczenie rakiety
w tym przedziale czasu czyli "x=x(t2)-x(t1)
b) Wyznacz błąd względny odniesiony do wartości dokładnej (ang.
true relative error)
7
Met.Numer. wykład 5
Wzór trapezów
f (a) + f (b)
ń#
a)
a = 8 s
I H" (b - a)Ą# b = 30 s
ó# Ą#
2
Ł# Ś#
140000
ń#
f (t) = 2000lnĄ# - 9.8t
ó#140000 - 2100t Ą#
Ł# Ś#
Ą# 140000 ń#
= 177.27 m / s
f (8) = 2000lnó# - 9.8(8)
Ł#140000 - 2100(8)Ą#
Ś#
Ą# 140000 ń#
f (30) = 2000lnó# - 9.8(30) = 901.67 m / s
Ł#140000 - 2100(30)Ą#
Ś#
177.27 + 901.67
ń#
I = (30 -8)Ą# = 11868 m
ó# Ą#
2
Ł# Ś#
8
Met.Numer. wykład 5
Wzór trapezów
Wartość dokładna
b)
30
�# 140000 ś#
ń#
"x = ś#2000lnĄ# - 9.8t ź#dt
=11061m
+"
ś# ó#140000 Ą# ź#
Ś#
8 Ł# - 2100t
�# #
Błąd względny:
11061 - 11868
"t = � 100 = 7.2959 %
11061
9
Met.Numer. wykład 5
3
2013-11-23
Złożony wzór trapezów
Błąd, jak pokazuje
Y
poprzedni przykład, f(x)
jest zbyt duży. Można
zaproponować podział
przedziału całkowania
na n segmentów,
każdy o długości h.
b - a
h =
n
b b
a - a - a b - a
b X
a + 2 a + 3
a +
4 4
dla n=4 4
b a+h a+2h a+3h a+4h
I = f (x)dx = f (x)dx + f (x)dx + f (x)dx + f (x)dx
+" +" +"+" +"
a a a+h a+2h a+3h
10
Met.Numer. wykład 5
Złożony wzór trapezów
b a+( n-1 )h b
a+h a+2h
f ( x )dx = f ( x )dx + f ( x )dx + ... + f ( x )dx + f ( x )dx
+" +" +" +" +"
a a a+h a+( n-2 )h a+( n-1 )h
f (a) + f (a + h) f (a + h) + f (a + 2h)
ń# ń#
= hĄ# + hĄ# + ...
ó# Ą# ó# Ą#
2 2
Ł# Ś# Ł# Ś#
f (a + (n -1)h) + f (b)
ń#
... + [b - (a + (n -1)h]Ą#
ó# Ą#
2
Ł# Ś#
b
n-1
b - a
Ą#
f ( x )dx = f ( a ) + 2ż# f ( a + ih )�# + f ( b )ń#
+" �#" Ź#
ó# Ą#
2n
a Ł# �#i=1 �# Ś#
11
Met.Numer. wykład 5
Złożony wzór trapezów
Przykład 2:
Prędkość rakiety w przedziale czasu od t1=8 s do t2=30 s jest
opisana wzorem:
140000
ń#
v(t) = 2000lnĄ# - 9.8t
ó#140000 - 2100t Ą#
Ł# Ś#
a) Przy pomocy metody trapezów znajdz
przemieszczenie rakiety
w tym przedziale czasu czyli "x=x(t2)-x(t1) przyjmując n=2
b) Wyznacz błąd względny odniesiony do wartości dokładnej (true
relative error)
12
Met.Numer. wykład 5
4
2013-11-23
Złożony wzór trapezów
n-1
b - a Ą#
a)
I = f (a) + 2ż# f (a + ih)�# + f (b)ń#
"
�# Ź#
ó# Ą#
2n �#i=1 �#
Ł# Ś#
b - a 30 - 8
n = 2 a = 8 s b = 30 s h = = = 11s
n 2
2-1
30 - 8 Ą#
I = f (8) + 2ż# f (a + ih)�# + f (30)ń#
"
�# Ź#
ó# Ą#
2(2)
�#i=1 �#
Ł# Ś#
22 22
= [f (8) + 2 f (19) + f (30)] = [177.27 + 2( 484.75 ) + 901.67]
4 4
= 11266 m
13
Met.Numer. wykład 5
Złożony wzór trapezów
b)wartość dokładna wynosi:
30
140000
�# Ą# ń# ś#dt
x =
+"ś#2000lnó#140000- 2100tĄ# - 9.8t ź# =11061m
8 �# Ł# Ś# #
Błąd względny to:
11061-11266
"t = �100 =1.8534%
11061
14
Met.Numer. wykład 5
Złożony wzór trapezów
n "x Et "t % "a %
1 11868 -807 7.296 ---
2 11266 -205 1.853 5.343
3 11153 -91.4 0.8265 1.019
4 11113 -51.5 0.4655 0.3594
5 11094 -33.0 0.2981 0.1669
6 11084 -22.9 0.2070 0.09082
7 11078 -16.8 0.1521 0.05482
8 11074 -12.9 0.1165 0.03560
15
Met.Numer. wykład 5
5
2013-11-23
Analiza błędu dla wzoru trapezów
Błąd bezwzględny metody z pojedynczym segmentem
( b - a )3
Et = f " ( ś ), a < ś < b
12
ś
gdzie jest punktem w [a,b]
Błąd w metodzie złożonej (wielosegmentowej) jest sumą błędów dla
każdego segmentu. Błąd pojedynczego segmentu wynosi:
[( a + h ) - a]3 f " ( ś1 ), a < ś1 < a + h
E1 =
12
h3
= f " ( ś1 )
12
16
Met.Numer. wykład 5
Analiza błędu dla wzoru trapezów
Podobnie:
3
[(a + ih) - (a + (i -1)h)]
Ei = f "(śi ), a + (i -1)h < śi < a + ih
12
h3
= f " ( śi )
12
dla n:
[b - {a + ( n - 1)h}]3 f " ( ś ), a + ( n - 1)h < ś < b
En =
n n
12
h3
= f " ( ś )
n
12
17
Met.Numer. wykład 5
Analiza błędu dla wzoru trapezów
Całkowity błąd w metodzie złożonej jest sumą błędów pojedynczych segmentów:
n
f " ( śi )
"
n
n
( b - a )3 i=1
Et = Ei h3 =
= f " ( śi )
" "
12
12n2 n
i=1 i=1
n
f " ( ś )
"
i
Wyrażenie
i=1
n
jest przybliżoną średnią wartością drugiej pochodnej w przedziale a < x < b
1
Et " ą
n2
18
Met.Numer. wykład 5
6
2013-11-23
Analiza błędu dla wzoru trapezów
30
�# 140000 ś#
ń#
Poniższa tablica dla całki
- 9.8t dt
ź#
+"ś# 2000lnĄ#
ó#140000 - 2100t Ą#
�# Ł# Ś# #
8
w funkcji liczby segmentów n. Widać, że gdy liczba segmentów jest
podwajana, błąd Et zmniejsza się w przybliżeniu czterokrotnie.
n Value
Et "t % "a %
2 11266 -205 1.854 5.343
4 11113 -51.5 0.4655 0.3594
8 11074 -12.9 0.1165 0.03560
16 11065 -3.22 0.02913 0.00401
19
Met.Numer. wykład 5
Całkowanie metodą Romberga
Metoda Romberga jest rozszerzeniem metody
trapezów i daje lepsze przybliżenie całki poprzez
zasadniczą redukcję błędu (true error)
20
Met.Numer. wykład 5
Ekstrapolacja Richardsona
Błąd (Et true error) w n-segmentowym wzorze trapezów wynosi
C
Et E"
n2
gdzie C jest współczynnikiem proporcjonalności
Stąd:
Et = TV - In
wartość dokładna (true
wartość przybliżona np.
value)
wyliczona ze wzoru trapezów
(approximate value)
21
Met.Numer. wykład 5
7
2013-11-23
Ekstrapolacja Richardsona
Można pokazać, że C
E" TV - I2n
(2n)2
gdy segment zostaje zmniejszony o połowę
C
Ze wzorów:
E" TV - In
2
(n)
C
E" TV - I2n
(2n)2
I2n - In
TV E" I2n +
otrzymujemy:
3
22
Met.Numer. wykład 5
Ekstrapolacja Richardsona
Przykład 3:
Prędkość rakiety w przedziale czasu od t1=8 s do t2=30 s jest
opisana wzorem:
140000
ń#
v(t) = 2000lnĄ# - 9.8t
ó#140000 - 2100t Ą#
Ł# Ś#
a) Przy pomocy ekstrapolacji Richardsona znajdz
przemieszczenie
rakiety w tym przedziale czasu czyli "x=x(t2)-x(t1)
b) Wyznacz błąd względny odniesiony do wartości dokładnej (true
relative error)
Przyjąć n=2
23
Met.Numer. wykład 5
Tabela wyników ze złożonego wzoru trapezów
do n=8 segmentów
n "x Et "t % "a %
1 11868 -807 7.296 ---
2 11266 -205 1.853 5.343
3 11153 -91.4 0.8265 1.019
4 11113 -51.5 0.4655 0.3594
5 11094 -33.0 0.2981 0.1669
6 11084 -22.9 0.2070 0.09082
7 11078 -16.8 0.1521 0.05482
8 11074 -12.9 0.1165 0.03560
24
Met.Numer. wykład 5
8
2013-11-23
Ekstrapolacja Richardsona
a)
I2 = 11266m I4 = 11113m
I2n - In dla n=2
TV E" I2n +
3
I4 - I2 11113 -11266
TV E" I4 + = 11113 +
3 3
= 11062m
25
Met.Numer. wykład 5
Ekstrapolacja Richardsona
b)
Wartość dokładna to:
30
�# 140000 ś#dt = 11061 m
ń#
x = 2000 lnĄ# - 9.8t
ś# ź#
+"
ó#140000 - 2100t Ą#
8 �# Ł# Ś# #
Stąd
=
Et = 11061-11062 -1 m
26
Met.Numer. wykład 5
Ekstrapolacja Richardsona
Błąd względny
c)
11061-11062
"t = �100 = 0.00904%
11061
Porównanie wyników z metodą złożoną trapezów
.
n
"x (m) "t % "x (m) "t %
wzór ekstrapolacja
wzór trapezów ekstrapolacja
trapezów Richardsona
Richardsona
1 11868 7.296 -- --
2 11266 1.854 11065 0.03616
4 11113 0.4655 11062 0.009041
8 11074 0.1165 11061 0.0000
27
Met.Numer. wykład 5
9
2013-11-23
Metoda Romberga
Całkowanie metodą Romberga stosuje ten sam wzór
co ekstrapolacja Richardsona. Jednakże, metoda
Romberga jest to algorytm rekurencyjny.
Przypomnijmy:
I2n - In
TV E" I2n +
3
Można zapisać
I2n - In I2n - In
(I2n )R = I2n + = I2n +
3 42-1 -1
Wartość dokładna TV jest zastąpiona przez wynik całkowania
metodą Richardsona
(I2n )R
Znak jest zastąpiony przez znak równości.
E"
28
Met.Numer. wykład 5
Metoda Romberga
Estymowana wartość dokładna wynosi:
TV E" (I2n )R + Ch4
gdzie Ch4 jest wartością błędu przybliżenia
Następna wartość całki (podwajając liczbę segmentów)
wynosi:
I4n - I2n
(I4n )R = I4n +
3
Estymowana wartość dokładna wynosi teraz:
(I4n) -(I2n )
R R
TV E" (I4n) +
R
15
(I4n) - (I2n)
R R
= (I4n )R +
43-1 -1
29
Met.Numer. wykład 5
Metoda Romberga
Ogólne wyrażenie w metodzie Romberga
Ik-1, j+1 - Ik-1, j
Ik , j = Ik-1, j+1 + ,k e" 2
4k -1 -1
Wskaz
nik k reprezentuje rząd ekstrapolacji
k=1 odpowiada wartościom uzyskanym ze wzoru trapezów
k=2 odpowiada wartościom uzyskanym z błędem O(h2)
Wskaz
nik j reprezentuje dokładność; j+1 daje całkę
wyznaczoną dokładniej niż j
30
Met.Numer. wykład 5
10
2013-11-23
Metoda Romberga
Przykład 4:
Prędkość rakiety w przedziale czasu od t1=8 s do t2=30 s jest
opisana wzorem:
140000
ń#
v(t) = 2000lnĄ# - 9.8t
ó#140000 - 2100t Ą#
Ł# Ś#
a) Przy pomocy wzoru Romberga znajdz
przemieszczenie rakiety
w tym przedziale czasu czyli "x=x(t2)-x(t1)
b) Wyznacz błąd względny odniesiony do wartości dokładnej (true
relative error)
Przyjąć n=1, 2, 4, 8
31
Met.Numer. wykład 5
Tabela wyników ze złożonego wzoru trapezów
do n=8 segmentów
n "x Et "t % "a %
1 11868 -807 7.296 ---
2 11266 -205 1.853 5.343
3 11153 -91.4 0.8265 1.019
4 11113 -51.5 0.4655 0.3594
5 11094 -33.0 0.2981 0.1669
6 11084 -22.9 0.2070 0.09082
7 11078 -16.8 0.1521 0.05482
8 11074 -12.9 0.1165 0.03560
32
Met.Numer. wykład 5
Metoda Romberga
Na podstawie tabeli odczytujemy:
I1,1 = 11868 I1,2 = 11266
I1,4 = 11074
I1,3 = 11113
W pierwszym przybliżeniu:
I1,2 - I1,1
I2,1 = I1,2 +
3
11266 -11868
= 11266 +
3
33
Met.Numer. wykład 5
11
2013-11-23
Metoda Romberga
Podobnie,
I1,4 - I1,3
I1,3 - I1,2
I2,3 = I1,4 +
I2,2 = I1,3 +
3
3
11113 -11266
11074 -11113
= 11113 +
= 11074 +
3
3
= 11062 = 11061
34
Met.Numer. wykład 5
Metoda Romberga
W drugim przybliżeniu,
I2,2 - I2,1
I3,1 = I2,2 +
15
11062 -11065
= 11062 +
15
= 11062
I2,3 - I2,2
I3,2 = I2,3 +
Podobnie,
15
11061 -11062
= 11061 +
15
= 11061
35
Met.Numer. wykład 5
Metoda Romberga
Dla trzeciego rzędu,
I3,2 - I3,1
I4,1 = I3,2 +
63
11061 -11062
= 11061 +
63
= 11061m
36
Met.Numer. wykład 5
12
2013-11-23
Metoda Romberga
Rząd 1 Rząd 2 Rząd 3
1-segment 11868
11065
2-segment 11262 11062
11062 11061
4-segment 11113 11061
11061
8-segment 11074
Poprawione wartości całki metodą Romberga
37
Met.Numer. wykład 5
Metoda Simpsona
Metoda trapezów była oparta na przybliżeniu funkcji
podcałkowej f(x) wielomianem stopnia pierwszego.
W metodzie Simpsona wykorzystuje się wielomiany
stopnia drugiego, jest to tzw. metoda parabol.
b b
I = f ( x )dx H" f2( x )dx
+" +"
a a
gdzie:
f2( x ) = a0 + a1x + a2x2
38
Met.Numer. wykład 5
Metoda Simpsona
Równanie paraboli dla 3
Y
f(x)
punktów:
( a, f ( a )),
�# a + b a + b ś#,
�# ś#
, f
ś# ś# ź#
ź#
2 2
�# �# # # fI(x)
( b, f ( b ))
f ( a ) = f2( a ) = a0 + a1a + a2a2
ab X
2
a + b a + b a + b a + b
�# ś# �#
f = f2 �# ś# = a0 + a1ś# ś# + a2 �# ś#
ś# ź# ś# ź# ź# ś# ź#
2 2 2 2
�# # �# # �# # �# #
f ( b ) = f2( b ) = a0 + a1b + a2b2
39
Met.Numer. wykład 5
13
2013-11-23
Metoda Simpsona
Wyznaczone współczynniki funkcji f2(x) to:
a + b
�# ś#
a2 f ( b ) + abf ( b ) - 4abf + abf ( a ) + b2 f ( a )
ś# ź#
2
�# #
a0 =
a2 - 2ab + b2
a + b a + b
�# ś# �# ś#
af ( a ) - 4af + 3af ( b ) + 3bf ( a ) - 4bf + bf ( b )
ś# ź# ś# ź#
2 2
�# # �# #
a1 = -
a2 - 2ab + b2
a + b
�# ś#
2�# f ( a ) - 2 f + f ( b )ś#
ś# ś# ź# ź#
2
�# �# # #
a2 =
a2 - 2ab + b2
40
Met.Numer. wykład 5
Metoda Simpsona
Wówczas:
b
I H" f2( x )dx
+"
a
b
= (a0 + a1x + a2 x2)dx
+"
a
b
Ą# x2 x3 ń#
= x + a1 + a2 Ą#
0
ó#a 2 3
Ł# Ś#a
b2 - a2 b3 - a3
= a0( b - a ) + a1 + a2
2 3
41
Met.Numer. wykład 5
Metoda Simpsona
b
b - a a + b
Ą# �# ś#
f2( x )dx = f ( a ) + 4 f + f ( b )ń#
ś# ź#
+"
ó# Ą#
6 2
a Ł# �# # Ś#
b - a
h =
2
Co daje:
b
h Ą# a + b ń#
�# ś#
f2 (x)dx = f (a) + 4 f + f (b)Ą#
ś# ź#
+" ó#
3 2
�# #
Ł# Ś#
a
wzór parabol
42
Met.Numer. wykład 5
14
2013-11-23
Metoda Simpsona
w wersji złożonej
b x2 x4
f(x)
f ( x )dx = f ( x )dx + f ( x )dx + .....
+" +" +"
a x0 x2
xn-2
xn
....+ f ( x )dx + f ( x )dx
+" +"
. . .
xn-4
xn-2
x
b
x0 x2 xn-2 xn
f ( x0 ) + 4 f ( x1 ) + f ( x2 )ń#
f ( x )dx = ( x2 - x0 )Ą# + ...
+"
ó# Ą#
6
Ł# Ś#
a
xi - xi-2 = 2h
f ( x2 )+ 4 f ( x3 )+ f ( x4 )ń# i = 2, 4, ..., n
+( x4 - x2 )Ą# + ...
ó# Ą#
6
Ł# Ś#
f ( xn-4 ) + 4 f ( xn-3 ) + f ( xn-2 )ń#
...+ ( xn-2 - xn-4 )Ą# + ...
ó# Ą#
6
Ł# Ś#
f ( xn-2 )+ 4 f ( xn-1 )+ f ( xn )ń#
+( xn - xn-2 )Ą#
ó# Ą#
6
Ł# Ś#
43
Met.Numer. wykład 5
Metoda Simpsona
b
f ( x0 )+ 4 f ( x1 ) + f ( x2 )ń#
f ( x )dx = 2hĄ# + ...
+"
ó# Ą#
6
a Ł# Ś#
f ( x2 ) + 4 f ( x3 )+ f ( x4 )ń#
+ 2hĄ# + ...
ó# Ą#
6
Ł# Ś#
f ( xn-4 ) + 4 f ( xn-3 ) + f ( xn-2 )ń#
+ 2hĄ# + ...
ó# Ą#
6
Ł# Ś#
f ( xn-2 ) + 4 f ( xn-1 ) + f ( xn )ń#
+ 2hĄ#
ó# Ą#
6
Ł# Ś#
44
Met.Numer. wykład 5
Metoda Simpsona
b
h
f ( x )dx = [f ( x0 ) + 4{f ( x1 )+ f ( x3 )+ ...+ f ( xn-1 )}+ ...]
+"
3
a
...+ 2{f (x2) + f (x4) + ...+ f (xn-2)}+ f (xn )}]
Ą# ń#
n-1 n-2
h
ó#
= f (x0 ) + 4 f (xi ) + 2 f (xi ) + f (xn )Ą#
" "
3 ó# Ą#
i=1 i=2
i=odd i=even
Ł# Ś#
Ą# ń#
n-1 n-2
b - a
ó#
= f (x0 ) + 4 f (xi ) + 2 f (xi ) + f (xn )Ą#
" "
3n ó# Ą#
i=1 i=2
i=odd i=even
Ł# Ś#
45
Met.Numer. wykład 5
15
2013-11-23
Metoda Simpsona  analiza błędu
Wartości przybliżone
przykładu stosując regułę 1/3 Simpsona z wieloma segmentami
n Wartość przybliżona Et |t |
2 11065.72 4.38 0.0396%
4 11061.64 0.30 0.0027%
6 11061.40 0.06 0.0005%
8 11061.35 0.01 0.0001%
10 11061.34 0.00 0.0000%
46
Met.Numer. wykład 5
Metoda Simpsona  analiza błędu
(b - a)5 (4)
Błąd dla jednego segmentu Et = - f (ś), a < ś < b
2880
Błąd w metodzie wielosegmentowej
( x2 - x0 )5 ( 4 ) = - h5 f ( 4 )
x0 < ś1 < x2
( ś1 ),
E1 = - f ( ś1 )
90
2880
( x4 - x2 )5 ( 4 ) h5 ( 4 ) x2 < ś2 < x4
E2 = - f ( ś2 ) = - f ( ś2 ),
2880 90
( x2i - x2( i-1 ) )5 ( 4 ) h5 ( 4 ) x2( < śi < x2i
Ei = - f ( śi ) - f ( śi ), i-1 )
=
2880
90
47
Met.Numer. wykład 5
Metoda Simpsona  analiza błędu
Całkowity błąd
n
n n
h5 2 (4) (b - a)5 2 (4)
2
"
Et = Ei = - i=1 f (śi ) = - "
f (śi )
"
90
i=1 90n5 i=1
n
2
(4)
f (śi )
"
(b - a)5 i=1
= -
90n4 n
n
2
(4)
f (śi ) (b - a)5 (4)
"
(4) Et = - f
i=1
f =
90n4
n
średnia wartość pochodnej
48
Met.Numer. wykład 5
16
2013-11-23
Metoda Gaussa
Całkę metodą kwadratury Gaussa
przedstawia wzór:
b
H" c1 f ( x1 )+ c2 f ( x2 )
I = f ( x)dx
+"
a
stałe współczynniki
Punkty x1 i x2 , w których określamy wartość funkcji podcałkowej nie
są ustalone (jak poprzednio na granicach przedziału czyli a i b), ale
są a priori dowolnie rozmieszczone w przedziale .
49
Met.Numer. wykład 5
Metoda Gaussa
Niewiadome x1, x2, c1, c2 znajduje się zakładając, że
funkcja podcałkowa spełnia warunek:
f ( x ) = a0 + a1x + a2x2 + a3x3.
b b
f ( x )dx = (a0 + a1x + a2x2 + a3x3)dx
+" +"
a a
b
Ą# x2 x3 x4 ń#
= x + a1 + a2 + a3 Ą#
0
ó#a 2 3 4
Ł# Ś#a
�# - a2 ź# ś# b3 - a3 ź# ś# b4 - a4 ź#
b2 ś# �# ś# �# ś#
ś#
= a0(b - a)+ a1ś# ź# + a2 ś# ź# + a3ś# ź#
2 3 4
�# # �# # �# #
50
Met.Numer. wykład 5
Metoda Gaussa
b
H" c1 f ( x1 )+ c2 f ( x2 )
I = f ( x)dx
+"
a
z jednej strony
b
f ( x )dx = c1(a0 + a1x1 + a2 x12 + a3x13)+ c2(a0 + a1x2 + a2 x2 2 + a3x23)
+"
a
= a0(c1 + c2 )+ a1(c1x1 + c2x2 )+ a2(c1x12 + c2 x2 2)+ a3(c1x13 + c2x23)
z drugiej strony
b b
f ( x )dx = (a0 + a1x + a2x2 + a3x3)dx
+" +"
a a
�# - a2 ź# ś# b3 - a3 ź# ś# b4 - a4 ź#
b2 ś# �# ś# �# ś#
ś#
= a0(b - a)+ a1ś# ź# + a2 ś# ź# + a3ś# ź#
2 3 4
�# # �# # �# #
51
Met.Numer. wykład 5
17
2013-11-23
Metoda Gaussa
a0(c1 + c2)+ a1(c1x1 + c2x2)+ a2(c1x12 + c2x22)+ a3(c1x13 + c2x23)
�# - a2 ź# ś# b3 - a3 ź# ś# b4 - a4 ź#
b2 ś# �# ś# �# ś#
ś#
= a0(b - a)+ a1ś# ź# + a2 ś# ź# + a3ś# ź#
2 3 4
�# # �# # �# #
b2 - a2
b - a = c1 + c2
= c1x1 + c2x2
2
b3 - a3
-
= c1x12 + c2x2 2 b4 a4 c1x13 + c2 x2 3
=
3
4
52
Met.Numer. wykład 5
Metoda Gaussa
Rozwiązując układ równań:
b - a = c1 + c2
b2 - a2
= c1x1 + c2 x2
2
b3 - a3
= c1x12 + c2x2 2
3
b4 - a4
= c1x13 + c2 x23
4
otrzymujemy dla metody dwupunktowej:
b
�# - a 1 b + a
ś#�# ś#
b
�# - a 1 b + a
ś#�# ś#
x2 = +
ś# ź#ś# ź#
x1 = +
ś# ź#ś#- ź#
2 3 2
�# #�# #
2 3 2
�# #�# #
b - a
b - a
c2 =
c1 =
2
2
53
Met.Numer. wykład 5
Metoda Gaussa
W dwu-punktowej metodzie Gaussa, całka funkcji f(x) wyraża
się wzorem:
b
f (x)dx H" c1 f (x1)+ c2 f (x2 )
+"
a
�# ś# �#
b - a b - a �# 1 ś# b + a b - a b - a �# 1 ś# b + a ś#
= f ś# + ź# + f ś# + ź#
ś#- ź# ś# ź#
ś# ź# ś# ź#
2 2 2 2 2 2
3 3
�# # �# #
�# # �# #
Uogólniając dla n punktów:
b
f ( x )dx H" c1 f ( x1 ) + c2 f ( x2 ) + . . . . . . . + cn f ( xn )
+"
a
54
Met.Numer. wykład 5
18
2013-11-23
Metoda Gaussa
W n-punktowej metodzie Gaussa, współczynniki ci oraz argumenty
xi są stabelaryzowane dla całek w granicach od -1 do 1
1
n
g( x )dx E" ci g( xi )
"
+"
i=1
-1
n współczynniki argumenty funkcji
2 c1 = 1.000000000 x1 = -0.577350269
c2 = 1.000000000 x2 = 0.577350269
3 c1 = 0.555555556 x1 = -0.774596669
c2 = 0.888888889 x2 = 0.000000000
c3 = 0.555555556 x3 = 0.774596669
4 c1 = 0.347854845 x1 = -0.861136312
c2 = 0.652145155 x2 = -0.339981044
c3 = 0.652145155 x3 = 0.339981044
c4 = 0.347854845 x4 = 0.861136312
55
Met.Numer. wykład 5
Metoda Gaussa
n współczynniki argumenty funkcji
5 c1 = 0.236926885 x1 = -0.906179846
c2 = 0.478628670 x2 = -0.538469310
c3 = 0.568888889 x3 = 0.000000000
c4 = 0.478628670 x4 = 0.538469310
c5 = 0.236926885 x5 = 0.906179846
6 c1 = 0.171324492 x1 = -0.932469514
c2 = 0.360761573 x2 = -0.661209386
c3 = 0.467913935 x3 = -0.2386191860
c4 = 0.467913935 x4 = 0.2386191860
c5 = 0.360761573 x5 = 0.661209386
c6 = 0.171324492 x6 = 0.932469514
56
Met.Numer. wykład 5
Metoda Gaussa
1
g( x )dx
Jeżeli dane są w tablicach dla to jak obliczamy
+"
-1
b
?
f ( x )dx
+"
a
Granice całkowania należy zamienić na
[a, b] [-1,1]
Niech x = mt + c
Wynika stąd, że:
Dla
x = a, -1
t =
b - a
m =
Dla t = 1 2
x = b,
b + a
c =
2
57
Met.Numer. wykład 5
19
2013-11-23
Metoda Gaussa
b - a b + a
Stąd:
x = t +
2 2
b - a
dx = dt
2
Ostatecznie:
b 1
b
�# - a b + a b - a
ś#
f (x)dx = f t + dt
ś# ź#
+" +"
2 2 2
�# #
a -1
58
Met.Numer. wykład 5
Metoda Gaussa
Przykład 5:
Prędkość rakiety w przedziale czasu od t1=8 s do t2=30 s jest
opisana wzorem:
140000
ń#
v(t) = 2000lnĄ# - 9.8t
ó#140000 - 2100t Ą#
Ł# Ś#
a) Przy pomocy dwu-punktowej metody Gaussa znajdz

przemieszczenie rakiety w tym przedziale czasu czyli "x=x(t2)-
x(t1)
b) Wyznacz błąd względny odniesiony do wartości dokładnej (true
relative error)
59
Met.Numer. wykład 5
Metoda Gaussa
Najpierw zmieniamy granice całkowania z [8,30] na [-1,1]
30 1
30 - 8 30 - 8 30 + 8ś#
�#
f ( t )dt = f x +
+" +" ś# ź#dx
2 -1
2 2
8 �# #
1
= 11 f (11x +19)dx
+"
-1
c1 = 1.000000000
Następnie odczytujemy z tablic dla n=2
x1 = -0.577350269
c2 = 1.000000000
x2 = 0.577350269
60
Met.Numer. wykład 5
20
2013-11-23
Metoda Gaussa
Korzystamy ze wzoru kwadratury Gaussa
1
11 f (11x +19)dx H" 11c1 f (11x1 +19)+11c2 f (11x2 +19)
+"
-1
= 11 f (11( -0.5773503) +19)+11 f (11( 0.5773503) +19)
= 11 f (12.64915 ) +11 f ( 25.35085 )
= 11( 296.8317 ) + 11( 708.4811)
= 11058.44 m
61
Met.Numer. wykład 5
Metoda Gaussa
Skorzystaliśmy z tego, że:
Ą# 140000 ń#
f (12.64915 ) = 2000lnó# - 9.8(12.64915 )
Ł#140000 - 2100(12.64915 )Ą#
Ś#
= 296.8317
Ą# 140000 ń#
f ( 25.35085 ) = 2000lnó# - 9.8( 25.35085 )
Ł#140000 - 2100( 25.35085 )Ą#
Ś#
= 708.4811
62
Met.Numer. wykład 5
Metoda Gaussa
b)
Et
Błąd bezwzględny (true error)
Et =11061.34 -11058.44
= 2.9000 m
c)
"t
Błąd względny,
11061.34 -11058.44
"t = �100%
11061.34
= 0.0262%
63
Met.Numer. wykład 5
21