handout - Vrije Universiteit Amsterdam

Onderwerpen
Numerieke Methoden I
College 5: Numerieke Integratie
(Hoofdstuk 5)
1. Introductie integratie
2. Toepassing in Bayesiaanse statistiek
3. Matlab functies
A.A.N. Ridder
4. Kwadratuurregels
Department EOR
Vrije Universiteit Amsterdam
Huispagina: http://personal.vu.nl/a.a.n.ridder/numprog/default.htm
I
Rechthoek
I
§5.1 Trapezium
I
§5.3 Simpson
I
§5.4 Gauss
30 september 2014
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
1 / 77
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
2 / 77
Probleemschets
R, continu.
I
Gegeven f : [a, b] →
I
Bekend uit de Analyse: dan is f Riemann integreerbaar.
Probleem
Bereken de bepaalde integraal
Toepassing
b
Z
I(f ) =
f (x) dx
a
a
c Ad Ridder (VU)
b
Numerieke Methoden I – Najaar 2014
3 / 77
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
4 / 77
Bayesiaanse Statistiek I
I
Bayesiaanse Statistiek II
Zij X een s.v. met de discrete logaritmische reeksverdeling:
θ
P(X = x) = −x log(1
,
− θ)
x
p(x|θ) =
I
Veronderstel θ is onbekend;
I
Veronderstel data (onafhankelijke realisaties van X)
I
y = (1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 4, 5, 5, 7)
P
Q
Dus aantal data is n = 19, de som is
yi = 41 en het product is yi = 16800;
I
Probleem: schat θ.
I
Klassiek: maximum likelihood:
voor x = 1, 2, . . ., met parameter θ ∈ (0, 1).
`(θ) = log
n
Y
p(yi |θ)
i=1
Etc.
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
5 / 77
Bayesiaanse Statistiek III
Numerieke Methoden I – Najaar 2014
Bayesiaans: we beschouwen θ als een realisatie van een stochastische variabele
Θ (op (0,1));
I
Als we geen of imperfecte informatie hebben, geven we Θ een voor de hand
liggende dichtheid;
I
De data y bezorgen ons informatie;
I
De voorwaardelijke simultane kansdichtheid dat de data zouden optreden als
gegeven is Θ = θ
f (y|θ) =
π(θ) = 6θ(1 − θ),
n
Y
P(X = yi |Θ = θ) =
i=1
Die heet prior, genoteerd π; bv een Beta(2, 2):
=
0<θ<1
n
Y
i=1
I
c Ad Ridder (VU)
6 / 77
Bayesiaanse Statistiek IV
I
I
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
7 / 77
n
Y
p(yi |θ)
i=1
P
θyi
(−1)n θ yi
−θ41
= Q
=
n
−yi log(1 − θ)
( yi )(log(1 − θ))
16800(log(1 − θ))19
Die heet de likelihood; ook wel L(θ);
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
8 / 77
Bayesiaanse Statistiek V
I
Bayesiaanse Statistiek VI
Omgekeerd, de voorwaardelijke kansdichtheid van Θ als de data y gegeven zijn, is
(pas de regel van Bayes toe):
π(θ|y) = fΘ (θ|X1 = y1 , . . . , Xn = yn ) = π(θ)f (y|θ)/c(y),
I
Die heet de posterior kansdichtheid;
I
Met normeringsconstante
I
E[Θk |y] =
0 < θ < 1,
I
In het bijzonder h(θ) ≡ 1 geeft de normeringsconstante: c(y) =
R
Daarna posteriormomenten via θk p∗ (θ) dθ.
π(θ)f (y|θ) dθ.
De ongenormaliseerde posterior (geen kansverdeling!) is
∗
p (θ) = π(θ)f (y|θ) ∝ π(θ|y)
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
I
9 / 77
Bayesiaanse Statistiek VII
I
θk π(θ|y) dθ
0
Wegens proportionaliteit is het voldoende de ongenormaliseerde posterior te
gebruiken voor diverse functies h
Z
h(θ)p∗ (θ) dθ
0
I
1
Z
I
1
Z
c(y) =
In Bayesiaanse statistiek is belangrijk verwachtingen te kunnen berekenen van de
posterior verdeling; bv de posteriormomenten
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
R
p∗ (θ) dθ;
10 / 77
Bayesiaanse Statistiek VIII
Bij de gegeven data y en prior π(θ):
p∗ (θ) = 6θ(1 − θ) ×
I
Integreer deze functie over (0, 1). Geeft c(y) = 2.1648 × 10−13 .
I
Dus de posterior dichtheid van Θ gegeven de data y (en prior π(θ)) is
π(θ|y) =
c Ad Ridder (VU)
I
−θ41
−θ42 (1 − θ)
=
.
16800(log(1 − θ))19
2800(log(1 − θ))19
−0.1650θ42 (1 − θ)1010
(log(1 − θ))19
Als schatting voor θ zouden we de verwachting kunnen nemen:
E[Θ|y] =
1
Z
θπ(θ|y) dθ = 0.7176
0
I
MLE (klassieke statistiek) zou schatting 0.7488 geven;
I
Uniforme prior π(θ) = 1 op (0,1), geeft 0.7338.
(0 < θ < 1).
Numerieke Methoden I – Najaar 2014
11 / 77
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
12 / 77
Bayesiaanse Statistiek IX
MATLAB
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
13 / 77
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
14 / 77
Beschikbare Functies
I
trapz: trapeziumregel voor univariate functies
I
quad: regel van Simpson voor univariate functies
I
quadl: Gauss-Lobatto kwadratuur voor univariate functies
I
dblquad: dubbele integraal over een rechthoek
I
Enz
§5.1 Kwadratuur
Zie Pratap §5.4.
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
15 / 77
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
16 / 77
Oppervlakte
Numerieke Integratie Methoden
I
De integraal kan ge¨ınterpreteerd worden als de oppervlakte tussen de grafiek van
f en de x-as (tussen a en b).
I
Numerieke methoden zijn gebaseerd op het benaderen van deze oppervlakte
door “eenvoudige” (meetkundige) figuren: kwadratuurregels.
I
I
Numerieke integratiemethoden woren ook wel kwadratuurregels genoemd;
I
Reflecterend aan klassieke technieken van de kwadratuur van een meetkundig
object.
A. Rechthoeksregel (of Middelpuntregel);
B. Trapeziumregel;
Historie: kwadratuur is het construeren van een vierkant met dezelfde oppervlakte
als een gegeven meetkundige figuur, bv driehoek, cirkel, ...
C. Simpson’s regel;
D. Gauss kwadratuur.
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
17 / 77
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
18 / 77
Middelpuntregel I
Middelpuntregel
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
19 / 77
I
Zie opgaven 6 en 6 van Exercises 5.1;
I
Verdeel [a, b] in n gelijke deelintervallen van lengte h = (b − a)/n;
I
Het i-de deelinterval is [xi , xi+1 ] met xi = a + ih (i = 0, . . . , n; x0 = a);
I
Noem mi = (xi + xi+1 )/2 het middelpunt van het i-de deelinterval;
I
Beschouw de rechthoek Ri met basis [xi , xi+1 ] en hoogte f (mi ), de functiewaarde
in het middelpunt (zie volgende slide);
I
Bepaal de oppervlakte van Ri .
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
20 / 77
Middelpuntregel II
Middelpuntregel III
functie f (x)
Oppervlakte Ri
= (xi+1 − xi )f (mi )
rechthoek Ri
I
De benadering van de integraal is de som van de oppervlaktes van de
rechthoeken:
I(f ) ≈ Mn (f ) =
n−1
X
(xi+1 − xi ) f (mi ) = h
i=0
n−1 X
2i + 1
f a+h
2
i=0
n−1 X
2i + 1
=h
f a + (b − a)
2n
i=0
xi
xi+1
mi
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
21 / 77
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
22 / 77
Middelpuntregel IV
§5.1 Trapeziumregel
x0
c Ad Ridder (VU)
x1
x2
Numerieke Methoden I – Najaar 2014
x3
23 / 77
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
24 / 77
Uniforme Opsplitsing
Trapeziumregel II
functie f (x)
xi+1
Z
I
`i (x) dx
Verdeel [a, b] in n gelijke deelintervallen van lengte h = (b − a)/n;
xi
I
Het i-de deelinterval is [xi , xi+1 ] met xi = a + ih (i = 0, . . . , n; x0 = a);
I
Bepaal de lineaire interpolatie tussen de grafiekpunten (xi , f (xi )) en (xi+1 , f (xi+1 ));
I
Dat geeft een functie `i (x) (rechte lijn);
I
Bepaal de oppervlakte tussen de x-as en de grafiek van `i (x) op [xi , xi+1 ].
I
Merk op dat de figuur een trapezium is.
I
Zie Figure 5.1 en volgende slide.
lineair
polynoom `i (x)
=
xi+1
xi
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
25 / 77
Trapeziumregel III
I
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
26 / 77
Trapeziumregel IV
De benadering van de integraal is de som van de oppervlaktes van de trapezia:
I(f ) ≈ Tn (f ) =
n−1
X
(xi+1 − xi )
i=0
1
f (xi ) + f (xi+1 )
2
h
=
f (x0 ) + 2f (x1 ) + 2f (x2 ) + · · · + 2f (xn−1 ) + f (xn )
2
n−1
X
h
f (x0 ) + f (xn ) + h
f (xi )
=
2
i=1
=
I
1
xi+1 − xi
2
× f (xi ) + f (xi+1 )
(1)
n−1
X
h
f (a) + 2
f (a + ih) + f (b)
2
i=1
x0
Zie Figure 5.2 en volgende slide.
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
27 / 77
c Ad Ridder (VU)
x1
x2
Numerieke Methoden I – Najaar 2014
x3
28 / 77
Integratiefout van Trapeziumregel
Bewijs Integratiefout bij e´ e´ n Trapezium I
R
(1). Schrijf T1 (f ) als een integraal: T1 (f ) = ab `(x) dx.
Hierin is `(·) de lijn door (a, f (a)) en (b, f (b)) (lineaire interpolatie).
Theorem 1
(2). Pas stelling over interpolatiefout toe (zie boek op pag. 181):
voor alle x ∈ [a, b] bestaat een ξ(x) ∈ [a, b] zodat
Veronderstel dat f tweemaal continu differentieerbaar is op [a, b]. Splits [a, b] op in n
gelijke intervallen van lengte h = (b − a)/n. Dan is er een ζ ∈ [a, b] zodat
f (x) − `(x) = 12 (x − a)(x − b)f 00 (ξ(x)).
1
1
I(f ) − Tn (f ) = − nh3 f 00 (ζ) = − (b − a)h2 f 00 (ζ) = O(h2 )
12
12
(3). Dus
Bewijs.
Z
I(f ) − T1 (f ) =
´ trapezium op [a, b].
Beschouw eerst n = 1; dwz e´ en
b
b
Z
f (x) − `(x) dx =
a
a
1
(x
2
− a)(x − b)f 00 (ξ(x)) dx.
(4). Pas middelwaardestelling voor integralen toe (zie volgende slide).
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
29 / 77
Middelwaardestelling voor Integralen
c Ad Ridder (VU)
30 / 77
Bewijs Integratiefout bij e´ e´ n Trapezium II
Stelling
(4). Pas middelwaardestelling voor integralen toe met
Veronderstel twee functies g1 en g2 op [a, b] met de eigenschappen
g1 (x) = f 00 (ξ(x)),
(i). g1 is continu op [a, b];
Z b
Z b
1
(x − a)(x − b) dx
I(f ) − T1 (f ) = g1 (η)
g2 (x) dx = f 00 (ξ(η))
2
a
a
−1
1
(b − a)3 = − h3 f 00 (ζ)
= f 00 (ζ)
12
12
(iii). g2 wisselt niet van teken op [a, b].
Dan is er een η ∈ [a, b] zodat
Z
b
Z
g1 (x)g2 (x) dx = g1 (η)
a
g2 (x) = 21 (x − a)(x − b).
Dan krijgen we
(ii). g2 is Riemann integreerbaar op [a, b];
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
b
g2 (x) dx.
(5)
a
Numerieke Methoden I – Najaar 2014
31 / 77
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
32 / 77
Bewijs Integratiefout bij n Trapezia
Toepassing van Integratiefout I
I
Beschouw nu weer de trapeziumregel bij een opsplitsing van [a, b] in n
deelintervallen van lengte h = (b − a)/n;
I
´ trapezium:
Pas op elk deelinterval de integratiefout toe van e´ en
I(f ) − Tn (f ) =
n−1 Z x
X
i+1
xi
i=0
f (x) dx − T1 (f |[xi ,xi+1 ] ) = −
n−1
X
i=0
1 3 00
h f (ζi ),
12
(7)
met ζi ∈ [xi , xi+1 ], i = 0, . . . , n − 1;
I
Er geldt
min f 00 (ζi ) ≤
i
X
f 00 (ζi )/n ≤ max f 00 (ζi )
I
Pas de tussenwaardestelling van continue
functies toe: er bestaat een
P
00
ζ ∈ [ζ0 , ξn−1 ] ⊂ [a, b] zodat f 00 (ζ) = n−1
i=0 f (ζi )/n;
I
Conclusie:
1
1 x
Numerieke integratie van I =
I
Vind het aantal trapezia n en intervalbreedte h zodat gegarandeerd
|I − Tn | < 10−8 is;
I
Relatie: nh = 5;
I
De integrand f (x) = 1/x heeft f 00 (x) = 2/x3 ;
I
Merk op dat f 00 (x) positief is en monotoon dalend is op [1, 6];
I
Dus voor alle x ∈ [1, 6] is
i
i
R6
I
dx volgens trapeziumregel Tn met n trapezia;
|f 00 (x)| ≤ f 00 (1) = 2;
n−1
I(f ) − Tn (f ) = −
c Ad Ridder (VU)
1
1 3 X 00
f (ζi ) = − nh3 f 00 (ζ)
h
12 i=0
12
Numerieke Methoden I – Najaar 2014
33 / 77
Toepassing van Integratiefout II
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
34 / 77
Als Aantal Deelintervallen Tweemacht is I
Bijzonder geval: n = 2k een tweemacht. Bv k = 3.
I
I
Dus
1
1 53
|I − Tn | = nh3 f 00 (ξ) ≤
.
12
6 n2
Dus n2 > (125/6) × 108 geeft n ≥ 45644. Bij n = 45644 is
h = 5/n = 0.000109544511501033.
1
2
2
2
2
2
2
2
1
x0
x1
x2
x3
x4
x5
x6
x7
x8
h
h
2h
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
35 / 77
c Ad Ridder (VU)
h
h
2h
h
h
2h
Numerieke Methoden I – Najaar 2014
h
h
2h
36 / 77
Als Aantal Deelintervallen Tweemacht is II
I
Zie Theorem 2;
I
Recursieve trapezium formule:
§5.3 Simpson’s Regel
h
Tn (f ) =
f (x0 ) + 2f (x2 ) + 2f (x4 ) + · · · + 2f (xn−2 ) + f (xn )
2
+ h f (x1 ) + f (x3 ) + · · · + f (xn−1 )
1
1
1
= Tn/2 (f ) + h f (x1 ) + f (x3 ) + · · · + f (xn−1 ) = Tn/2 (f ) + Mn/2 (f )
2
2
2
I
Geeft Tn/2 (f ) + Mn/2 (f ) = 2Tn (f );
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
37 / 77
Uniforme Opsplitsing
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
38 / 77
Simpson’s Regel II
functie f (x)
I
Verdeel [a, b] in n gelijke deelintervallen van lengte h = (b − a)/n;
I
Aanname: n is even!
I
Het i-de deelinterval is [xi , xi+1 ] met xi = a + ih (i = 0, . . . , n; x0 = a);
I
Beschouw achtereenvolgens het i − 1-ste en het i-de interval tezamen voor
i = 1, 3, 5, . . .;
I
Bepaal de kwadratische interpolatie tussen de grafiekpunten (xi−1 , f (xi−1 )),
(xi , f (xi )) en (xi+1 , f (xi+1 ));
I
Dat geeft een functie (parabool) qi (x);
I
Bepaal de oppervlakte tussen de x-as en de grafiek van qi (x) op [xi−1 , xi+1 ];
I
Zie Figure 5.5 en volgende slide;
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
xi+1
Z
qi (x) dx
kwadratisch
polynoom qi (x)
xi−1
=
xi−1
39 / 77
c Ad Ridder (VU)
xi
1
(xi+1 − xi−1 )
6
× f (xi−1 ) + 4f (xi ) + f (xi+1 )
xi+1
Numerieke Methoden I – Najaar 2014
40 / 77
Simpson’s Regel III
Simpson’s Regel IV
De benadering van de integraal is de som van deze oppervlaktes:
I(f ) =
n/2 Z
X
x2i
f (x) dx
x2i−2
i=1
≈ Sn (f ) =
n/2
X
1
(x2i − x2i−2 ) f (x2i−2 ) + 4f (x2i−1 ) + f (x2i )
6
i=1
h
f (x0 ) + 4f (x1 ) + 2f (x2 ) + 4f (x3 ) + · · · + 2f (xn−2 ) + 4f (xn−1 ) + f (xn )
3
n/2
n/2−1
X
X
h
=
f (x0 ) + 4
f (x2i−1 ) + 2
f (x2i ) + f (xn )
3
i=1
i=1
=
=
n/2
n/2−1
X
X
h
f (a) + 4
f (a + (2i − 1)h) + 2
f (a + 2ih) + f (b)
3
i=1
i=1
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
x0
41 / 77
Simpson’s Regel V
x1
x2
c Ad Ridder (VU)
x3
x4
x5
x6
Numerieke Methoden I – Najaar 2014
42 / 77
Integratiefout bij Simpson’s Regel
De samengestelde Simpson’s formule op slide 41 kan worden herschreven:
Stelling
h
f (x0 ) + 2f (x2 ) + 2f (x4 ) + · · · + 2f (xn−2 ) + f (xn )
3
4 + h f (x1 ) + f (x3 ) + · · · + f (xn−1 )
3
2
1
= Tn/2 (f ) + Mn/2 (f )
3
3
Sn (f ) =
Veronderstel dat f viermaal continu differentieerbaar is op [a, b]. Splits [a, b] op in n
gelijke intervallen van lengte h = (b − a)/n. Dan is er een ξ ∈ [a, b] zodat
I(f ) − Sn (f ) = −
1
1
nh5 f (4) (ξ) = −
(b − a)h4 f (4) (ξ) = O(h4 )
180
180
(8)
Bewijs: Aanaloog als voor de trapeziumregelfout (slides 29–33); iets complexer.
Dat geeft de relatie tussen trapezium, middelpunt en Simpson’s regel.
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
43 / 77
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
44 / 77
Grote Oh Symbool
Als Aantal Deelintervallen Tweemacht is
I
Fouten worden meestal slechts aangegeven in termen van orde van grootte;
I
Trapeziumregel:
I
I
per deelinterval van lengte h: I(f ) − T1 (f ) = O(h3 );
I
n deelintervallen van lengte h = (b − a)/n: I(f ) − Tn (f ) = nO(h3 ) = O(h2 );
I
Herinner recursieve trapezium formule Tn/2 (f ) + Mn/2 (f ) = 2Tn (f );
I
Gebruik in Simpson’s regel:
1
2
Tn/2 (f ) + Mn/2 (f )
3
3
1
2
=
Tn/2 (f ) + Mn/2 (f ) − Tn/2 (f )
3
3
4
1
= Tn (f ) − Tn/2 (f )
3
3
1
= Tn +
Tn (f ) − Tn/2 (f )
3
Sn (f ) =
Simpson’s regel:
I
per deelinterval van lengte h: I(f ) − S1 (f ) = O(h5 );
I
n deelintervallen van lengte h = (b − a)/n: I(f ) − Sn f ) = nO(h5 ) = O(h4 );
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
45 / 77
Adaptieve Recursieve Integratie I
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
46 / 77
Adaptieve Recursieve Integratie II
Volgens de help faciliteit implementeert Matlab’s quad een adaptieve recursieve
Simpson’s regel.
(a). Recursief:
meer intervallen
c Ad Ridder (VU)
minder intervallen
Numerieke Methoden I – Najaar 2014
47 / 77
I
We hebben gezien dat Tn/2 , Tn en Sn de integraal I(f ) benaderen waarbij
Sn = Tn + (Tn (f ) − Tn/2 (f ))/3;
I
Term (Tn − Tn/2 )/3 kan worden gebruikt als foutschatting;
I
Dit gaat ook zo voor Sn/2 , Sn , en de kwadratuurregel Qn die ontstaat door 4-de graads
polynomen te gebruiken;
I
Er geldt Qn = Sn + (Sn − Sn/2 )/15;
I
Term (Sn − Sn/2 )/15 kan worden gebruikt als foutschatting;
I
Verdubbel het aantal deelintervallen en bereken Simpson’s regel net zo lang totdat
|(Sn − Sn/2 )/15| < .
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
48 / 77
Adaptieve Recursieve Integratie III
Adaptieve Recursieve Integratie IV
(b). Adaptief:
I
¨
Merk op: voor de volgende iteratie (in de recursie) wordt elk deelinterval in tweeen
gedeeld;
I
Merk op: per deelinterval hoeven slechts S1 en S2 berekend te worden;
I
Voer deze deling niet meer uit voor deelintervallen waarvoor al |(S2 − S1 )/15| klein
genoeg is;
I
Namelijk |S2 − S1 | < × lengte deelinterval/(b − a).
Wetende hoe het algoritme precies werkt, kan een voorbeeld geconstrueerd worden
om quad om de tuin te leiden;
Moet nog bij worden gezegd dat quad begint door eerst [a, b] in drie intervallen splitst
en dan op die intervallen S1 , S2 berekent enz.
Dit voorbeeld is zo geconstrueerd dat op het eerste startinterval [a, a + h]:
S1 (a, a + h) = S1 (a, a + h/2) + S1 (a + h/2, a + h) = S2 (a, a + h)
Idem op de overige startintervallen. Inderdaad: quad geeft waarde 158.67 terwijl
trapz en gaudl allebei 99.97 geven.
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
49 / 77
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
50 / 77
Beschrijving
§5.4 Gauss Regels
Rb
I
Herinner: I(f ) =
I
Kies n + 1 knooppunten in [a, b]: a ≤ x0 < x1 < · · · < xn ≤ b.
a
f (x) dx;
R.
I
Kies n + 1 gewichten A0 , A1 , . . . , An ∈
I
De benadering van I(f ) is de bijbehorende kwadratuur:
I(f ) ≈ Gn (f ) =
n
X
Ai f (xi ).
(1)
i=0
I
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
51 / 77
Interpretatie: benader oppervlak onder de grafiek van f door oppervlak van een
aantal rechthoeken.
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
52 / 77
Gauss Kwadratuurregels II
Voorbeeld I
De knoopunten en gewichten worden op een speciale manier geconstrueerd en
gebruikt.
Constructie
I
n = 2, a = 0, b = 7;
De benadering is exact (dwz Gn (f ) = I(f )) voor alle polynomen f op [a, b] t/m graad
2n + 1.
I
De knooppunten zijn (x0 , x1 , x2 ) = (0.7889, 3.5000, 6.2111) en de gewichten zijn
(A0 , A1 , A2 ) = (1.9444, 3.1111, 1.9444);
I
Pas toe op f (x) = 30 + 64x − 56x2 + 14x3 − x4 ; (vierde graadspolynoom);
R
Eenvoudig: I(f ) = 07 f (x) dx = 417.433333;
P
Ook eenvoudig: G3 (f ) = 2i=0 Ai f (xi ) = 417.433333.
I
Eigenschap
De knooppunten en gewichten hangen alleen af van n, a, b.
I
I
Dwz, exact.
Gebruik
Dezelfde knooppunten en gewichten worden voor alle integreerbare functies op [a, b]
gebruikt.
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
53 / 77
Voorbeeld I in Figuren
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
54 / 77
Voorbeeld II
A2
A0
I
Nog steeds n = 2, a = 0, b = 7;
I
Pas nu toe op f (x) = x4 e−x +
I
I
A1
f (x0 )
0
I
7
0
Een vierde graads polynoom f (x) op
[0, 7].
c Ad Ridder (VU)
x0
x1
x2
1
x+0.5
− 12 ;
R
Mbv trapeziumregel (106 trapezia): I(f ) = 07 f (x) dx ≈ 19.056252;
P
Mbv 3-punts Gausskwadratuur: G3 (f ) = 3i=1 wi f (xi ) = 19.324277.
Relatieve fout van 1.4%.
7
Drie
rechthoeken
met
breedtes
A0 , A1 , A2 en hoogtes f (x0 ), f (x1 ), f (x2 ).
Numerieke Methoden I – Najaar 2014
55 / 77
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
56 / 77
Voorbeeld II in Figuren
Gauss Kwadratuurregels III
A1
Uitwerking van de eis van slide 53.
A2
I
2n + 2 onbekende variabelen x0 , x1 , . . . , xn , A0 , . . . , An ;
I
Construeer 2n + 2 vergelijkingen in die variabelen;
I
Vergelijking k is Gn (f ) = I(f ) voor f = xk ;
I
Hierin is
Gn (f ) =
n
X
I(f ) =
I
Een integreerbare f (x) op [0, 7].
c Ad Ridder (VU)
7
n
X
0
x0
x1
x2
7
i=0
Ai xik =
bk+1 − ak+1
k+1
(1)
Drie rechthoeken met breedtes A0 , A1 en
A2 , en met hoogtes f (x0 ), f (x1 ), f (x2 ).
Numerieke Methoden I – Najaar 2014
57 / 77
Gauss Kwadratuurregels IV
I
bk+1 − ak+1
.
k+1
De vergelijkingen zijn: voor k = 0, 1, . . . , 2n + 1
f (x0 )
0
xk dx =
a
i=0
A0
b
Z
Ai xik ;
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
58 / 77
Terug naar Voorbeelden I en II
I
n = 2, a = 0, b = 7;
I
Stelsel van 6 vergelijkingen:
R
Wegens lineariteit van is dan voor elk polynoom
f (x) = p0 + p1 x + · · · + p2n+1 x2n+1 van graad ≤ 2n + 1 (ga na):
Gn (f ) = I(f ).
I
Stelsel (1) is een stelsel van 2n + 2 niet-lineaire vergelijkingen in de 2n + 2
onbekenden x0 , . . . , xn en A0 , . . . , An ;
I
In principe oplosbaar.
I
k=0:
A0 + A1 + A2 = 7
k=1:
A0 x0 + A1 x1 + A2 x2 = 72 /2
k=2:
A0 x02 + A1 x12 + A2 x22 = 73 /3
k=3:
A0 x03 + A1 x13 + A2 x23 = 74 /4
k=4:
A0 x04 + A1 x14 + A2 x24 = 75 /5
k=5:
A0 x05 + A1 x15 + A2 x25 = 76 /6
Oplossing
x = (0.7889, 3.5000, 6.2111) A = (1.9444, 3.1111, 1.9444)
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
59 / 77
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
60 / 77
Gauss Kwadratuurregels V
I
Niet-lineair stelsel is lastig;
I
Alternatieve procedure;
Gauss Kwadratuurregels VI
1. Bepaal (of probeer) eerst knooppunten x0 , . . . , xn ;
I
Stelsel (2) is een stelsel van n + 1 lineaire vergelijkingen in de n + 1 onbekenden
A0 , . . . , An ;
I
In matrix/vector notatie XA = y met
2. Bepaal dan gewichten zo dat voor k = 0, 1, . . . , n
n
X
k+1
k
Ai xi =
i=0
I
b
−a
k+1
1
1 ···
1
x
x
x
·
·
·
x
 0
n
1
1

 2
x12 x22 · · · xn2 
x
X =  03

x0 x13 x23 · · · xn3 


.....................
x0n x1n x2n · · · xnn
Dan is voor elk polynoom f (x) = p0 + p1 x + · · · + pn xn van graad ≤ n (ga na):
Gn (f ) = I(f )
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
61 / 77
Gauss Kwadratuurregels VII
c Ad Ridder (VU)




y=




b−a

(b2 − a2 )/2


(b3 − a3 )/3


(b4 − a4 )/4


..

.
n+1
n+1
(b
−a
)/(n + 1)
Numerieke Methoden I – Najaar 2014
62 / 77
Voorbeelden I en II revisited
Matlab programma, bij gegeven integratie-interval [a, b], graad n + 1, en knooppunten x
(als kolomvector).
I
n = 2 Gauss kwadratuur proberen met x0 = 1, x1 = 3.5, x2 = 6;
I
Stelsel van 3 vergelijkingen:
function A = gaussweights(n,a,b,x)
X = zeros(n+1,n+1);
r = zeros(n+1,1);
A0 + A1 + A2 = 7
A0 x0 + A1 x1 + A2 x2 = 72 /2
A0 x02 + A1 x12 + A2 x22 = 73 /3
X(1,:) = ones(1,n+1);
y(1) = b-a;
for i=2:n+1
X(i,:) = X(i-1,:) .* x’;
y(i) = (bˆi - aˆi)/i;
end
I
Oplossing A = (2.2867, 2.4267, 2.2867);
I
Voorbeeld I:
I
A = X\y;
end
c Ad Ridder (VU)

1
k+1
Numerieke Methoden I – Najaar 2014
63 / 77
I
Benadering G3 (f ) =
I
Relatieve fout 7.5 %.
P2
Ai f (xi ) = 448.875000; exact I(f ) = 417.433333;
P3
Ai f (xi ) = 18.166383; exact I(f ) = 19.056252;
i=0
Voorbeeld II:
I
Benadering G3 (f ) =
I
Relatieve fout 4.7 %.
c Ad Ridder (VU)
i=0
Numerieke Methoden I – Najaar 2014
64 / 77
Overeenkomst met Interpolatie
I
I
Gauss Kwadratuurregels VIII
Gegeven de knooppunten x0 , . . . , xn ;
De gewichten Ai voldoen aan
b
Z
`i (x) dx;
Ai =
a
n
Y
x − xj
`i (x) =
x
− xj
j=0 i
I
In plaats van voor elk interval [a, b] en elke n een niet-lineair stelsel op te lossen,
I
Doen we dit alleen voor het interval [−1, 1] (wel voor elke n),
I
Op een speciale manier.
I
De knooppunten oplossing x en gewichten oplossing A kunnen worden
opgeslagen in een tabel;
I
Namelijk n = 1, 2, . . . en bij elke n twee vectoren van lengte n + 1.
I
Zie Table 5.1;
I
De knooppunten en gewichten bij [a, b] volgen door translatie en substitutie, zie
slide 69.
j6=i
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
65 / 77
Gauss Kwadratuurregels IX
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
66 / 77
Legendre Polynomen
Die worden gedefinieerd door q0 (x) = 1, q1 (x) = x,
qn (x) =
2n − 1
n−1
xqn−1 (x) −
qn−2 (x),
n
n
n = 2, 3, . . . .
Afspraken
(A). Beschouw alleen interval [−1, 1]; dwz a = −1 en b = 1;
(B). Knooppunten x0 , . . . , xn zijn de nulpunten van het n + 1-ste graads Legendre
polynoom;
(C). Gewichten A0 , . . . , An zijn oplossing van lineair stelsel (2).
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
67 / 77
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
68 / 77
Gauss Kwadratuurregels X
Gauss Kwadratuurregels XI
I
Van afspraak naar algemeen;
I
Gegeven f : [a, b] →
I
Bereken knooppunten −1 ≤ x0 < x1 < · · · < xn ≤ 1 en gewichten A0 , . . . , An
volgens de afspraken;
I
R;
1. Zeer goede foutanalyse;
2. Eenvoudig aan te passen dat de knooppunten nulpunten zijn van andere
orthogonale polynomen; zie Theorem 2.
Gebruik de substitutieregel van integratie
Z
b
f (x) dx =
a
I
Voordelen:
Concludeer
Gn (f ) =
c Ad Ridder (VU)
b−a
2
Z
1
f
−1
b + a + x(b − a)
2
dx
3. Eenvoudig uit te breiden naar meerdimensionale ingtegralen.
n
X
b + a + xi (b − a)
b−a
Ai f
2
2
i=0
Numerieke Methoden I – Najaar 2014
69 / 77
Gauss Kwadratuurregels XII
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
Matlab Voorbeeld I
Bereken:
3
Z
I
70 / 77
De functie quadl in Matlab is gebaseerd op Gauss kwadratuur;
cx−1 sin(dx−1 ) dx.
1
I
Waarin het eerste en het laatste knooppunt de eindpunten vormen van het
integratie-interval):
x0 = −1, xn = 1
I
Aantal n + 1 wordt zo gekozen dat de fout ‘klein genoeg’ is (weer via recursief
adaptief);
60
I
Resterende 2n onbekenden x1 , . . . , xn−1 en A0 , . . . , An zorgen dat de benadering
exact is voor polynomen tot en met graad 2n − 1;
30
I
Namelijk oplossing van
(−1)k A0 +
n−1
X
i=1
I
Plot bij c = 100 en d = 10.
0
Ai xik + An =
1 − (−1)k+1
2
(k = 0, . . . , 2n − 1);
−30
Heet Gauss-Lobatto.
c Ad Ridder (VU)
−60
1.0
Numerieke Methoden I – Najaar 2014
71 / 77
c Ad Ridder (VU)
1.5
2.0
2.5
Numerieke Methoden I – Najaar 2014
3.0
72 / 77
Matlab Voorbeeld II
Matlab Voorbeeld III
% trapeziumregel voor 1,2,4,...1024 trapezia
function trapeziumregel(a,b,c,d)
n = 1;
for k=1:10
x = linspace(a,b,n+1);
y = integrand(x,c,d);
q = trapz(x,y); disp(q);
n = 2*n;
end
end
function intvb
main;
function main
a = 1; b = 3; c = 100; d = 10;
trapeziumregel(a,b,c,d);
simpsonregel(a,b,c,d);
gausskwadratuur(a,b,c,d);
end
function simpsonregel(a,b,c,d)
y = @(x)integrand(x,c,d);
q = quad(y,a,b); disp(q);
end
function y = integrand(x,c,d)
y = c./x.*sin(d./x);
end
function gausskwadratuur(a,b,c,d)
y = @(x)integrand(x,c,d);
q = quadl(y,a,b); disp(q);
end
end % EOF
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
73 / 77
Meerdimensionale Integratie I
I
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
74 / 77
Meerdimensionale Integratie II
Numeriek mogelijk voor integratie op een ‘rechthoekige’ ruimte
[a1 , b1 ] × [a2 , b2 ], × · · ·
I
Bijvoorbeeld twee-dimensionaal:
Z
b2
Voorbeeld:
Z
Z
b1
f (x1 , x2 ) dx1 dx2
a2
0
a1
I
Middelpuntregel, Simpson’s regel, Gauss kwadratuur zijn eenvoudig naar hogere
dimensies uit te breiden;
I
Bijvoorbeeld Gauss kwadratuur in 2-D:
n1
n2
X
X
Ai1 Ai2 f xi1 , xi2
1
1
Z
sin2 π(x − y) dx dy
0
q = dblquad(@(x,y)(sin(pi*(x-y))).ˆ2,0,1,0,1)
i1 =0 i2 =0
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
75 / 77
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
76 / 77
Meerdimensionale Integratie III
I
Slechte prestatie bij veel dimensies;
I
¨ modellen;
Zoals in Bayesiaanse statistiek; financiele
I
Moeilijk toepasbaar bij niet-rechthoekige integratiegevieden;
I
Alternatief: simulatie (Monte Carlo Integratie).
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
77 / 77

Download Report