1. No category

10.4 Richardsonova ekstrapolacija
Z razpolavljanjem h pridemo pri sestavljenih formulah do toˇcnejˇsih rezultatov. Iz pribliˇzkov pri razliˇcnih h lahko ocenimo napako in ugotovimo, ˇce je potrebno ˇse nadaljnje
razpolavljanje.
Naj bo Sh(f ) Simpsonova formula s korakom h in Rh(f ) napaka Simpsonove formule pri
koraku h. Vemo, da velja
Zb
f (x)dx = Sh(f ) + Rh(f ) = Sh/2(f ) + Rh/2(f ).
I(f ) =
a
Za napaki velja
−(b − a)h4 (4)
Rh(f ) =
f (ξ1),
180
−(b − a)h4 (4)
Rh/2(f ) =
f (ξ2).
16 · 180
ˇ predpostavimo, da je f (4)(ξ1) ≈ f (4)(ξ2), dobimo
Ce
Rh(f ) ≈ 16Rh/2(f ).
Bor Plestenjak - NM 2 2009
Richardsonova ekstrapolacija
Iz
Rh/2(f )
=
I(f ) − Sh/2(f )
=
Sh(f ) + Rh(f ) − Sh/2(f )
≈
Sh(f ) + 16Rh/2(f ) − Sh/2(f )
dobimo
Rh/2(f ) ≈
Sh/2(f ) − Sh(f )
15
.
To formulo lahko uporabimo za oceno napake Simpsonove formule.
Poleg ocene lahko iz Sh/2(f ) in Sh(f ) z ekstrapolacijo dobimo ˇse boljˇsi pribliˇzek, saj je
I(f ) = Sh/2(f ) + Rh/2(f ) ≈
16Sh/2(f ) − Sh(f )
Podobno lahko naredimo pri trapezni in sredinski formuli.
Bor Plestenjak - NM 2 2009
15
.
Richardsonova ekstrapolacija - zgled
Z1
Raˇcunamo
I =
x2
e dx = 1.462651745907182.
0
S Simpsonovim sestavljenim pravilom dobimo:
S(0.1) = 1.462681400099797
in
S(0.05) = 1.462653624886297.
Ocena za napako je
(S(0.05) − S(0.1))/15 = −1.85 · 10
prava ocena za napako pa je
I − S(0, 0.05) = −1.87 · 10
Vidimo, da je v tem primeru ocena izredno dobra.
Bor Plestenjak - NM 2 2009
−6
.
−6
,
10.5 Adaptivne metode
Pri adaptivnih metodah sproti ocenjujemo napako in ˇce je potrebno dodatno delimo
podintervale.
Metoda se prilagaja funkciji, kar pomeni, da pri pohlevnem obnaˇsanju uporablja veˇcji h, pri
divjem obnaˇsanju pa manjˇsi h.
Bor Plestenjak - NM 2 2009
Adaptivna metoda na osnovi Simpsonove metode
function [Q,etol] = quadtx(F,a,b,tol)
c = (a + b)/2;
fa = F(a); fc = F(c); fb=F(b);
[Q,etol] = quadtxstep(F, a, b, tol, fa, fc, fb);
function [Q,etol] = quadtxstep(F,a,b,tol,fa,fc,fb)
h = b - a; c = (a + b)/2;
fd = F((a+c)/2); fe = F((c+b)/2);
Q1 = h/6 * (fa + 4*fc + fb);
Q2 = h/12 * (fa + 4*fd + 2*fc + 4*fe + fb);
etol=abs(Q2-Q1);
if etol <= tol
Q = Q2 + (Q2 - Q1)/15;
else
[Qa,etol1] = quadtxstep(F, a, c, tol/2, fa, fd, fc);
[Qb,etol2] = quadtxstep(F, c, b, tol-etol1, fc, fe, fb);
Q = Qa + Qb;
etol=etol1+etol2;
end
Bor Plestenjak - NM 2 2009
10.6 Rombergova metoda
Izrek 1. Za neskonˇcnokrat zvezno odvedljivo funkcijo f velja Euler-Maclaurinova sumacijska formula
Zb
f (x)dx = Th(f ) −
I(f ) =
a
∞
X
k=1
B2k 2k
(2k−1)
(2k−1)
h
f
(b) − f
(a) .
(2k)!
Bistveno je, da ima napako obliko
Zb
2
4
6
f (x)dx = Th(f ) + a1h + a2h + a3h + · · · ,
a
koeficienti ai pa so neodvisni od h.
Bor Plestenjak - NM 2 2009
Euler-Maclaurinova sumacijska formula
Zb
f (x)dx = Th(f ) −
I(f ) =
a
∞
X
k=1
B2k 2k
(2k−1)
(2k−1)
f
(b) − f
(a) .
h
(2k)!
Bk so Bernoullijeva ˇstevila, ki so doloˇcena z razvojem
x
=
ex − 1
∞
X
Bk
k=0
k!
k
x ,
|x| < 2π.
Vsa Bernoullijeva ˇstevila so racionalna, vsa liha ˇstevila razen B1 pa so enaka 0. Nekaj prvih
Bernoullijevih ˇstevil je
B0 = 1,
1
B1 = − ,
2
B2 =
1
,
6
B4 = −
1
,
30
B6 =
1
,
42
....
Od tod sledi
h2 0
h4
0
000
000
I(f ) = Th(f ) −
(f (b) − f (a)) +
(f (b) − f (a)) + · · · ,
12
720
vrsta praviloma ni konvergentna, je pa asimptotska, kar pomeni, da velja, ko gre h → 0.
Bor Plestenjak - NM 2 2009
Rombergova ekstrapolacija
2
4
6
I(f )
=
Th (f ) + a1,0 h + a2,0 h + a3,0 h + · · ·
I(f )
=
Th/2 (f ) + a1,0
I(f )
=
Th/4 (f ) + a1,0
h 2
2
h 2
4
+ a2,0
+ a2,0
h 4
2
h 2
4
+ a3,0
+ a3,0
h 6
2
h 6
4
+ ···
+ ···
ˇ enaˇcbo za h/2 pomnoˇzimo s 4 in odˇstejemo od enaˇcbe za h, se znebimo ˇclena h2 in dobimo toˇcnejˇsi
Ce
pribliˇzek:
I(f )
I(f )
(1)
4
6
=
Th/2 (f ) + a2,1 h + a3,1 h + · · ·
=
(1)
Th/4 (f ) + a2,1
h 4
2
+ a3,1
h 6
2
+ ···,
kjer sta
(1)
Th/2 (f ) =
4Th/2 (f ) − Th (f )
Postopek sedaj nadaljujemo v
3
,
(1)
Th/4 (f ) =
(2)
4Th/4 (f ) − Th/2 (f )
6
3
8
I(f ) = Th/4 (f ) + a3,2 h + a4,2 h + · · · ,
kjer je
(1)
(2)
Th/4 (f ) =
Bor Plestenjak - NM 2 2009
(1)
16Th/4 (f ) − Th/2 (f )
15
.
.
Rombergova shema
V sploˇsnem postopku tvorimo shemo
O(h2)
(0)
Th (f )
(0)
Th/2(f )
napaka
O(h4)
(1)
Th/4(f )
(0)
Th/8(f )
Th/8(f )
O(h8)
Th/2(f )
(0)
Th/4(f )
O(h6)
(1)
Th/4(f )
(2)
(1)
Th/8(f )
(3)
(4)
Th/8(f )
kjer je sploˇsna formula
T
(j)
h/2k
Bor Plestenjak - NM 2 2009
(j−1)
(f )
h/2k
4j T
(f ) =
(j−1)
(f )
h/2k−1
−T
4j − 1
.
···
Zgled uporabe Rombergove metode
Z Rombergovo metodo bomo izraˇcunali
in naredimo dve razpolavljanji.
R 2.2 ln xdx = 0.5346062. Zaˇcnemo s h = 0.6
1
Dobimo:
(0)
=
Th/2
(0)
=
(0)
=
Th
Th/4
1
1
0.6( ln 1.0 + ln 1.6 + ln 2.2) = 0.5185394
2
2
1 (0)
T + 0.3(ln 1.3 + ln 1.9) = 0.5305351
2 h
1 (0)
T
+ 0.15(ln 1.15 + ln 1.45 + ln 1.75 + ln 2.05) = 0.5335847
2 h/2
Pomembno je, da Th/2k vedno raˇcunamo kot
Th/2k =
1
h
Th/2k−1 + k (y1 + y3 + · · · + y2k −1).
2
2
Tako vsako funkcijsko vrednost izraˇcunamo enkrat in imamo z Rombergom zanemarljivo
(0)
dodatnega dela v primerjavi z raˇcunanjem T k , rezultat pa je lahko mnogo natanˇcnejˇsi.
h/2
Bor Plestenjak - NM 2 2009
Rezultat po uporabi Rombergove ekstrapolacije
Za primerjavo, toˇcen rezultat je
R 2.2 ln xdx = 0.5346062.
1
Z Rombergovo ekstrapolacijo dobimo
(0)
(1)
Th/2
=
(0)
4Th/2 − Th
(0)
(1)
Th/4
=
(0)
4Th/4 − Th/2
3
(1)
(2)
Th/4
Bor Plestenjak - NM 2 2009
=
= 0.5345337
3
= 0.5346013
(1)
16Th/4 − Th/2
15
= 0.5346058.
10.7 Gaussove kvadraturne formule
Zb
Integral
f (x)ρ(x)dx,
a
kjer je ρ nenegativna uteˇz, aproksimiramo s kvadraturno formulo
Zb
f (x)ρ(x)dx =
a
n
X
(n)
(n)
Ai f (xi ) + R(f ).
i=0
(n)
Ai
R
b
=
L
Koeficienti so doloˇceni z vozli, saj velja
a
poljubni izbiri vozlov toˇcna za polinome stopnje vsaj n.
n,i (x)ρ(x)dx,
formula pa je pri
S primerno izbiro vozlov lahko doseˇzemo, da bo formula toˇcna za polinome stopnje vsaj
2n + 1, v ozadju pa so ortogonalni polinomi.
Bor Plestenjak - NM 2 2009
Ortogonalni polinomi
Za funkcije lahko na [a, b] definiramo skalarni produkt kot
Zb
hf, gi: =
f (x)g(x)ρ(x)dx.
a
Funkciji f in g sta ortogonalni, ˇce je hf, gi = 0.
Iz standardne baze polinomov
2
1, x, x , . . .
z ortogonalizacijo dobimo ortonormirano bazo
P0(x), P1(x), P2(x), . . . ,
kjer je Pi polinom stopnje i in velja hPi, Pk i = δik .
Bor Plestenjak - NM 2 2009
Toˇ
cnost do stopnje vsaj 2n + 1
(n)
Naj bo Pn+1 = kn+1(x−x0 ) · · · (x−x(n)
n ) tak normiran polinom, da je hPn+1 , qi = 0
za vsak polinom q stopnje kveˇcjemu n. Pri Gaussovi kvadraturni formuli za vozle izberemo
(n)
niˇcle x0 , . . . , x(n)
n , torej je
(n)
(n)
ω(x) = (x − x0 ) · · · (x − xn ).
Poljuben polinom f stopnje 2n + 1 lahko zapiˇsemo kot f (x) = q(x)ω(x) + r(x), kjer
sta q, r polinoma stopnje kveˇcjemu n. To pomeni:
Zb
Zb
f (x)ρ(x)dx
=
a
q(x)ω(x)ρ(x)dx +
a
=
Zb
0+
X
X
r(x)ρ(x)dx
a
n
(n)
(n)
Ai r(xi )
i=0
n
=
(n)
(n)
Ai f (xi ).
i=0
Pokazali smo, da je pravilo toˇcno za vse polinome stopnje 2n + 1 ali manj.
Bor Plestenjak - NM 2 2009
Napaka Gaussovih kvadraturnih formul
(n)
Iz teorije ortogonalnih polinomov sledi, da so vse niˇcle xi
(a, b).
enostavne, realne in leˇzijo na
Koeficiente Ai lahko izraˇcunamo z integriranjem Lagrangevih koeficientov, ˇse bolje pa je,
ˇce uporabimo Darboux-Cristoffelove formule:
(n)
Ai
=
Pn
1
2 (x(n) )
P
j=0 j
i
,
i = 0, . . . , n.
Lema 1.
Uteˇzi Gaussovih kvadraturnih pravil so pozitivne.
Izrek 2.
Za f ∈ C (2n+2)[a, b] velja
Zb
f (x)ρ(x) =
a
Bor Plestenjak - NM 2 2009
n
X
i=0
(n)
(n)
Ai f (xi )
f (2n+2)(ξ)
+
.
2
(2n + 2)!kn+1
Znani ortogonalni polinomi
Najpogostejˇsi ortogonalni polinomi so:
[−1, 1],
ρ(x) = 1,
[−1, 1],
ρ(x) = (1 − x2)− 2 ,
[−1, 1],
[−1, 1],
ρ(x) = (1 − x2) 2 ,
ρ(x) = (1 − x)α(1 + x)β ,
[−1, 1],
[0, ∞),
ρ(x) = (1 − x2)σ− 2 ,
ρ(x) = xσ e−x,
(−∞, ∞),
ρ(x) = e−x ,
(Legendre)
ˇ sev 1. vrste)
(Cebiˇ
ˇ sev 2. vrste)
(Cebiˇ
1
1
1
2
α, β > −1,
(Jacobi)
σ > 12 ,
σ > −1
(Gegenbauer)
(Laguerre)
(Hermite).
Za te ortogonalne polinome imamo tabelirane niˇcle in Gaussove kvadraturne formule.
Bor Plestenjak - NM 2 2009
Gauss–Legendrove kvadraturne formule
Gauss–Legendrovi kvadraturni formuli na dveh in treh toˇckah sta
Z1
f (x)dx = f
−1
Z1
f (x)dx =
−1
5
f
9
r1!
−
3
+f
r3!
−
5
r1!
3
1 (4)
+
f (ξ),
135
5
8
+ f (0) + f
9
9
r3!
5
+
1
(6)
f (ξ).
15750
Za primerjavo, pri trapeznem in Simpsonovem pravilu dobimo
Z1
2 (2)
f (x)dx = f (−1) + f (1) − f (ξ),
3
−1
Z1
1
4
1
1 (4)
f (x)dx = f (−1) + f (0) + f (1) +
f (ξ).
3
3
3
90
−1
Bor Plestenjak - NM 2 2009
10.8 Izlimitirani integrali
Zb
Pol v krajiˇsˇcu:
I =
a
g(x)
dx,
p
(x − a)
0 < p < 1,
kjer je g zvezna na [a, b].
Zb
Po definiciji je
I = lim
→0
g(x)
dx,
p
a+ (x − a)
a je ta konvergenca lahko prepoˇcasna, da bi bila uporabna za praktiˇcno raˇcunanje. Zato
iˇsˇcemo boljˇse moˇznosti.
Bor Plestenjak - NM 2 2009
Pol v krajiˇsˇ
cu 1
Zb
Raˇcunamo
I =
a
g(x)
dx,
(x − a)p
0 < p < 1,
kjer je g zvezna na [a, b].
Varianta 1: I razdelimo na I = I1 + I2, kjer sta
Z a+
I1 =
a
g(x)
dx,
p
(x − a)
I2 =
Zb
g(x)
dx.
p
a+ (x − a)
I2 izraˇcunamo s standardnimi metodami, pri raˇcunanju I1 pa g razvijemo v Taylorjevo
vrsto okoli a:
(x − a)2 00
g (a) + · · ·
g(x) = g(a) + (x − a)g (x) +
2
0
in dobimo
0
I1 = Bor Plestenjak - NM 2 2009
1−p
2 00
g(a)
g (a)
g (a)
+
+
+ ···
1−p
1!(2 − p)
2!(3 − p)
!
.
Pol v krajiˇsˇ
cu 2
Zb
Raˇcunamo
I =
a
g(x)
dx,
(x − a)p
0 < p < 1,
kjer je g zvezna na [a, b].
Varianta 2: g razvijemo v Taylorjevo vrsto okoli a:
g(x) = Ps(x) + ostanek
I zapiˇsemo kot I = I1 + I2, kjer sta
Zb
I1 =
a
Ps(x)
dx,
(x − a)p
I2 =
Z b g(x) − Ps(x)
a
(x − a)p
dx.
I1 lahko izraˇcunamo eksplicitno, pri I2 lahko uporabimo standardne metode.
Bor Plestenjak - NM 2 2009
Pol v krajiˇsˇ
cu 3
Zb
Raˇcunamo
I =
a
g(x)
dx,
(x − a)p
0 < p < 1,
kjer je g zvezna na [a, b].
Varianta 3: Uporabimo substitucijo, npr.
m
(x − a) = t ,
m−1
dx = pt
Z (b−a)1/p
Dobimo
I =p
p
k−1
g(a + t )t
0
in lahko uporabimo standardne metode.
Bor Plestenjak - NM 2 2009
dt,
k
,
m=
1−p
dt
k ∈ N.
Pol v krajiˇsˇ
cu 4
Varianta 4: pomagamo si lahko tudi z Gaussovimi kvadraturnimi formulami. Npr., ˇce
imamo
Z1
I =
potem je uteˇz √ 1
1−x2
g(x)
dx,
√
2
1−x
−1
ˇ seve kvadraturne formule.
in uporabimo Gauss-Cebiˇ
Primer:
Z1
g(x)
π
dx =
√
3
1 − x2
−1
Bor Plestenjak - NM 2 2009
√
g
−
3
2
!
√
+ g(0) + g
3
2
!!
+ Cf
(6)
(ξ).
Neskonˇ
cen interval
Z∞
Raˇcunamo I =
f (x)dx. I razdelimo na I = I1 + I2, kjer sta
a
Zb
I1 =
Z∞
f (x)dx,
I2 =
a
f (x)dx.
b
I1 izraˇcunamo normalno s standardnimi metodami. Variante za I2 so
• I2 zanemarimo, ˇce poznamo kakˇsno dobro oceno za |I2|.
• Naredimo substitucijo u = 1/x in dobimo integral po konˇcnem intervalu
Z0
I2 = −
g(1/u)u
−2
du.
1/b
Tretja varianta je, da za I uporabimo ustrezno Gaussovo kvadraturno formulo.
Bor Plestenjak - NM 2 2009