handout - Vrije Universiteit Amsterdam

Onderwerpen
Numerieke Methoden I
College 4: Interpolatie
(Hoofdstuk 4)
1. Introductie en voorbeeld van interpolatie
2. Matlab functies
A.A.N. Ridder
3. §4.1 Polynomiale interpolatie
Department EOR
Vrije Universiteit Amsterdam
I
Lagrange
I
Vandermonde
I
Newton
4. §4.2 Fouten analyse van interpolatie
Huispagina: http://personal.vu.nl/a.a.n.ridder/numprog/default.htm
23 september 2014
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
c Ad Ridder (VU)
1 / 41
Numerieke Methoden I – Najaar 2014
2 / 41
Data
R2 (i = 0, 1, . . . , n);
I
Gegeven zijn n + 1 punten (xi , yi ) ∈
I
Verkregen via simulatie of metingen;
Introductie
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
3 / 41
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
4 / 41
Probleem
Interpolatie
Interpolerende modellen zijn polynomen en splines;
I
Gevraagd een model op te stellen voor die data;
I
M.a.w. een functie f :
“representeert”;
I
Vele modelleringen zijn mogelijk:
R → R die deze punten “beschrijft” of “reproduceert” of
I
Specifieke klasse van functies (bv polynomen; splines; ...)?
I
Continue functies?
I
Differentieerbare functies (een, twee, drie, ... keer)?
I
Interpolerende functies (f (xi ) = yi voor alle i) of hoeft dat niet?
polynomiale interpolatie alleen door
eerste 10 punten(!)
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
5 / 41
Regressie
c Ad Ridder (VU)
cubic spline
Numerieke Methoden I – Najaar 2014
6 / 41
Smoothing
Fittende (benaderde) modellen zijn regressie en smoothing.
smoothing met splines
lineaire regressie
c Ad Ridder (VU)
smoothing met splines
smoothing met splines
kwadratische regressie
Numerieke Methoden I – Najaar 2014
7 / 41
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
8 / 41
Spline regressie
Overzicht
regressie met splines
I
Vandaag: polynomiale interpolatie;
I
Later in Numerieke Methoden I: splines;
I
Regressie en smoothing in andere colleges van (bedrijfs)econometrie;
regressie met splines
Deze functies interpoleren een deelverzameling van de data!
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
9 / 41
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
10 / 41
Beschikbare Functies
Vele, bv
Matlab
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
11 / 41
I
polyfit: polynomiale interpolatie
I
interp1: splines (lineair en cubic)
I
spline: cubic splines
I
Zie Pratap §5.2;
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
12 / 41
Voorbeeld
I
Programma
xdata = [-4,-1,3,5,8,10];
ydata = [2,4,-9,1,12,0];
n = length(xdata)-1; % graad van polynoom
p = polyfit(xdata,ydata,n);
disp(p)
x = linspace(min(x)-1,max(x)+1,301); % x-as voor plotten
y = polyval(p,x);
plot(xdata,ydata,’dr’,x,y,’k’,’LineWidth’,2,’MarkerSize’,10);
grid on
set(gca,’GridLineStyle’,’:’);
legend(’data’, ’polynoom’,’Location’,’SouthEast’);
title(’Polynomiale interpolatie’,’FontSize’,16,’FontWeight’,’bold’);
Data 6 punten (n = 5):
xi
yi
−4
2
−1
4
3
−9
5
1
I
De x-coordinaten heten knopen of knooppunten;
I
Gevraagd 5-de graads interpolerend polynoom
8
12
10
0
p(x) = p0 + p1 x + p2 x2 + p3 x3 + p4 x4 + p5 x5 ,
zodat
p(xi ) = yi
c Ad Ridder (VU)
(i = 0, 1, . . . , 5).
Numerieke Methoden I – Najaar 2014
13 / 41
Uitvoer
I
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
14 / 41
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
16 / 41
Plot
p = polyfit(xdata,ydata,n) met n = 5 geeft
0.002, −0.055, 0.320, 0.989, −6.488, −3.099.
I
Dit zijn de 6 coefficienten van een 5-de graadspolynoom gerepresenteerd door de
rijvector p.
I
De coefficient van de hoogste macht x5 staat voorop.
I
Dus het polynoom is
p(x) = 0.002x5 − 0.055x4 + 0.320x3 + 0.989x2 − 6.488x − 3.099.
I
Plot, zie volgende slide.
I
Zie ook Pratap §5.2.2.
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
15 / 41
Intermezzo
Horner’s Algoritme
I
y = polyval(p,x)
I
Zie pp 8-11 in het boek;
I
met p de vector van de 6 coefficienten van het 5-de graadpolynoom,
I
Efficient algoritme om polynomen te evalueren;
I
en met x een vector met de 301 punten die de ‘x-as’ representeren,
I
Zo is Matlab’s polyval geprogrammeerd;
I
berekent voor elke x(i) (i = 1, . . . , 301) de bijbehorende y(i) volgens
5
4
p(x) = 0.002x5 − 0.055x4 + 0.320x3 + 0.989x2 − 6.488x − 3.099
!
= −3.099 + x − 6.488 + x 0.989 + x 0.320 + x − 0.055 + 0.002x
3
y(i) = p(1) × x(i) + p(2) × x(i) + p(3) × x(i)
+ p(4) × x(i)2 + p(5) × x(i) + p(6).
I
Het resultaat staat in de vector y van lengte 301.
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
I
17 / 41
Complexiteit (aantal rekenoperaties) is 2n voor n-de graads polynoomevaluatie.
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
18 / 41
Probleem
§4.1 Polynomiale Interpolatie
I
I
R2 (i = 0, 1, . . . , n), met xi 6= xj voor alle i, j;
Vind een n-de graad polynoom p : R → R die deze punten verbindt; i.e.,
Gegeven n + 1 punten (xi , yi ) ∈
p(x) = p0 + p1 x + p2 x2 + · + pn xn ;
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
19 / 41
c Ad Ridder (VU)
zodat p(xi ) = yi ∀ i;
Numerieke Methoden I – Najaar 2014
20 / 41
Discussie
Existentie
Stelling (Theorem 1, blz 157)
Als alle xi 6= xj dan bestaat er een uniek n-de graad interpolerend polynoom.
I
Motivering?
I
In welke situaties doe je deze interpolatie?
I
Geschikt voor veel data?
I
Bestaat er altijd zo een polynoom?
I
Zijn er meer oplossingen of is een gevonden polynoom uniek?
I
Welke methoden en algorititmes zijn beschikbaar?
I
Wanneer nog zelf programmeren ipv Matlab functies?
Bewijs van existentie:
I
Definieer voor i = 0, 1, . . . , n de kardinale polynomen
`i (x) =
I
Y x − xj
.
x − xj
j6=i i
(2)
Dan is de Lagrange vorm van interpolerend polynoom:
p(x) =
n
X
(1)
yi `i (x);
i=0
I
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
21 / 41
Bewijs
Bewering: p(x) is een oplossing.
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
22 / 41
Existentie en Uniciteit
I
Gezocht coefficienten p0 , p1 , . . . , pn van polynoom p zodat voor i = 0, . . . , n:
p(x0 ) = y0 ⇔ p0 + p1 x0 + p2 x02 + · · · + pn x0n = y0
p(x1 ) = y1 ⇔ p0 + p1 x1 + p2 x12 + · · · + pn x1n = y1
I
Elke `i (x) is een n-de graad polynoom;
I
Som van n-de graad polynomen is een n-de graad polynoom;
I
Substitueer x = xk in (1):
(
1
`i (xk ) =
0
..
.
..
.
p(xn ) = yn ⇔ p0 + p1 xn + p2 xn2 + · · · + pn xnn = yn
voor i = k;
voor i =
6 k;
I
⇒ p(xk ) = yk ;
Dit is een stelsel lineaire vergelijkingen Vp = y, met onbekende p en gegeven
matrix V en gegeven rechterkant y:

1 x0 x02 · · · x0n
2
n
1 x 1 x
· · · x1 
1

V=
. . . . . . . . . . . . . . . . . . . .
2
n
1 xn xn · · · xn

c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
23 / 41
c Ad Ridder (VU)
 
y0
y 1 
 
y=.
 .. 
Numerieke Methoden I – Najaar 2014
yn
24 / 41
Existentie en Uniciteit (vervolg)
Voorbeeld
I
Dit is een bekende matrix, een zogeheten Vandermonde matrix;
I
Als alle xi 6= xj , dan is de matrix V niet-singulier; i.e.,
xi
yi
1
1

1
V=
1
1
1

det(V) 6= 0
I
Stelsel is oplosbaar met oplossing p = V −1 y;
I
Dat wil zeggen dat een oplossing bestaat;
I
En dat de oplossing uniek is.
I
´ interpolerend polynoom.
Dus er is precies e´ en
−4
2
−1
4
3
−9
−4
−1
3
5
8
10
16
1
9
25
64
100
−64
−1
27
125
512
1000
5
1
8
12
256
1
81
625
4096
10000
10
0

−1024
−1 

243 

3125 
32768 
100000
Stelsel Vp = y oplossen in Matlab via p = V \ y geeft waarschuwing dat de matrix V
slecht geconditioneerd is.
Dat wil zeggen dat een kleine wijziging in de data een heel ander oplossing geeft.
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
25 / 41
Illustratie
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
26 / 41
Numerieke Methoden I – Najaar 2014
28 / 41
Uniciteit
Oplossing van p = V \ y is
−3.0993 − 6.4882 0.9885 0.3203 − 0.0552 0.0020
Verander V(5, 2) van 8 in 10. Oplossing is
I
Alternatief bewijs van uniciteit;
I
Zie boek blz 157.
−1.1430 − 5.2819 0.1636 0.2910 − 0.0128 − 0.0013
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
27 / 41
c Ad Ridder (VU)
Newton’s Vorm
Newton’s gedeelde differenties
Definieer f [xi ] = yi , en voor 0 ≤ i < j ≤ n:
f [xi , . . . , xj ] =
I
Derde constructie (na Lagrange en Vandermonde);
I
Het interpolerende polynoom wordt geconstrueerd door
p(x) =
n
X
ai
i=0
I
f [xi+1 , . . . , xj ] − f [xi , . . . , xj−1 ]
.
xj − xi
(13)
De tabel van gedeelde differenties (voor n = 3):
i−1
Y
(x − xi ).
(3)
xi
x0
f [·]
f [x0 ]
x1
f [x1 ]
f [·, ·]
f [·, ·, ·]
f [·, ·, ·, ·]
f [x0 , x1 ]
j=0
¨
¨
Waarin de coeffici
enten
ai bepaald worden als de divided differences (gedeelde
differenties) van de gegevensvector y (zie volgende slide).
f [x0 , x1 , x2 ]
f [x1 , x2 ]
x2
f [x2 ]
x3
f [x3 ]
f [x0 , x1 , x2 , x3 ]
f [x1 , x2 , x3 ]
f [x2 , x3 ]
De constanten ak in Newton’s interpolatieformule staan in de bovenste diagonale rij:
ak = f [x0 , . . . , xk ] (k = 0, . . . , n).
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
29 / 41
Voorbeeld
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
30 / 41
Bewijs Newton’s interpolatie
−1
2
xi
yi
xi
−1
f [·]
2
f [·, ·]
−4−2
1−(−1)
1
5
5
10
f [·, ·, ·]
5−(−3)
3−(−1)
=2
2−5
5−1
− 34
=5
6
10−6
5−3
3
6
Lemma
Zij p(x) een k-de graads polynoom die (xi , yi )ki=0 interpoleert. En Zij q(x) een k-de
graads polynoom die (xi , yi )k+1
i=1 interpoleert. Dan is
f [·, ·, ·, ·]
= −3
−4
6−(−4)
3−1
3
1
−4
=
r(x) =
− 11
24
een k + 1-ste graads polynoom dat (xi , yi )k+1
i=0 interpoleert.
=2
Bewijs Lemma: substitutie.
10
Stelling (Theorem 2; blz 163
Bovenste rij geeft de constanten ak :
ak = f [x0 , . . . , xk ] in Newton’s interpolatieformule.
p(x) = 2 − 3(x + 1) + 2(x + 1)(x − 1) −
c Ad Ridder (VU)
(x − x0 )q(x) − (x − xk+1 )p(x)
xk+1 − x0
11
(x + 1)(x − 1)(x − 3).
24
Numerieke Methoden I – Najaar 2014
Bewijs Stelling: inductie naar k en het Lemma gebruiken.
31 / 41
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
32 / 41
Methoden en Algoritmen
Discussie
Waarom verschillende algoritmes?
We hebben gezien:
1. Lagrange vorm;
Ze produceren allemaal hetzelfde interpolerende polynoom p bij gegeven data
(xi , yi )ni=0 .
2. Vandermonde matrix;
Elk algoritme heeft voordeel/nadeel t.o.v. de andere algoritmes.
3. Newton’s vorm;
Bv:
4. Er zijn nog vele andere methoden, zie pp 169–172; en de website
http://mathworld.wolfram.com/topics/Interpolation.html
I
Vandermonde matrix is slecht geconditioneerd;
I
De polyfit van Matlab geeft de coefficienten van het interpolerend polynoom;
NB: Matlab’s polyfit is gebaseerd op 2 (Vandermonde matrix), maar lost
bijbehorende stelsel geavanceerd en efficient op.
I
Lagrange geeft direct de polynoomwaarde in elke x (maar niet de coefficienten);
I
Als een datapunt (xn+1 , yn+1 ) wordt toegevoegd, vereist de Newton’s formule veel
minder rekenwerk dan de Lagrange formule om het nieuwe polynoom te bepalen;
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
33 / 41
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
34 / 41
Probleem
Gegeven
§4.2 Fouten in
Polynomiale Interpolatie
R;
I
Een functie f : [a, b] →
I
knoopunten a ≤ x0 < x1 < · · · < xn ≤ b;
I
yi = f (xi ), i = 0, . . . , n;
I
p(·) het interpolerend polynoom door (xi , yi )ni=0 ;
Polynoom p(·) kan beschouwd worden als een benadering van de functie f (·) op het
interval [a, b];
Duidelijk is dat f (xi ) − p(xi ) = 0 in alle knopen xi ;
Gevraagd de fout te bepalen in elk punt x ∈ [a, b], dwz
f (x) − p(x).
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
35 / 41
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
36 / 41
Stelling
Bewijs
I
Kies een x ∈ [a, b];
I
Definieer de functies w(t) =
Stelling (Theorem 1; blz 181)
− xi ) en
n
Y
1
f (n+1) (ξ)
(x − xi ).
(n + 1)!
i=0
I
φ(t) heeft n + 2 nulpunten: x0 , . . . , xn , x.
Pas nu de stelling van Rolle n + 1 keer toe (op φ en zijn afgeleiden)!
I
φ0 heeft n + 1 nulpunten;
(2)
Numerieke Methoden I – Najaar 2014
37 / 41
Intermezzo: Uit de Analyse
f (x) − p(x)
w(t).
w(x)
I
I
c Ad Ridder (VU)
i=0 (t
φ(t) = f (t) − p(t) −
Veronderstel dat f (·) een n + 1 keer continu differentieerbare functie is, en dat p(·) het
n-de graads polynoom is dat de punten (xi , f (xi ))ni=0 interpoleert, waarbij
a ≤ x0 < x1 < · · · < xn ≤ b. Dan bestaat voor elke x ∈ [a, b] een punt ξ ∈ (a, b) zodat
f (x) − p(x) =
Qn
I
φ00 heeft n nulpunten;
I
···
I
φ(n+1) heeft 1 nulpunt, zeg ξ;
Werk uit φ(n+1) (ξ) = 0.
Merk op dat p(n+1) (ξ) = 0 en w(n+1) (ξ) = (n + 1)!.
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
38 / 41
Illustratie
f (x) = (x2 + 1)−1 op [−5, 5] en 21 punten xi ∈ [−5, 5] (i = 0, 1, . . . , 20).
Stelling van Rolle
Als een functie f voldoet aan de voorwaarden:
1. f is continu op het gesloten interval [a, b];
2. f is differentieerbaar op het open interval (a, b);
3. f (a) = f (b);
dan bestaat er een getal c in het open interval (a, b), waarin de afgeleide van f gelijk is
aan 0 (f 0 (c) = 0).
Punten op gelijke afstand:
xi = −5 + 12 i.
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
39 / 41
c Ad Ridder (VU)
Chebyshev punten:
xi = 5 cos(iπ/20).
Numerieke Methoden I – Najaar 2014
40 / 41
Toelichting
I
Bekend fenomeen van polynomiale interpolatie is het Runge’s effect: het
interpolerend polynoom p geeft steeds grotere afwijkingen van de functie f als het
aantal punten n toeneemt en de knopen op gelijke afstand liggen;
I
Het effect is veel minder bij een slimmere keuze van knopen x0 < x1 < · · · < xn ;
I
Bv zogeheten Chebyshev abscissen; dwz nulpunten van Chebyshev polynomen;
c Ad Ridder (VU)
Numerieke Methoden I – Najaar 2014
41 / 41

Download Report