Hoofdstuk 1. Illustratie 2

Hoofdstuk 1

Numerical Methods College 2 A. Floating-point representatie (Hoofdstuk 1) B. Matlab

Twee belangrijke onderwerpen die moeten leiden tot een beter begrip van de numerieke problematiek: 1. §1.3 Floating-point representatie

A.A.N. Ridder

2. §1.4 Loss of significance

normalsize Department EOR Vrije Universiteit Amsterdam

Deze leiden tot analyse waarom fouten optreden.

Huispagina: http://personal.vu.nl/a.a.n.ridder/numprog/default.htm

7 september 2015

c Ad Ridder (VU)

Numerical Methods – Periode 1 (2015-2016)

1 / 49

Illustratie 1

c Ad Ridder (VU)


2 / 49


4 / 49

Illustratie 2

>> format long >> x = 2.2; >> x = x+0.2

Een MATLAB sessie >> x = 59/19; y = 59 - 19*x

x = 2.400000000000000

y = 0

>> x = x+0.2

>> x = 59/190; y = 59 - 190*x

x = 2.600000000000001

y = -7.1054e-15

>> x = 2.2; x = x+0.4 x = 2.600000000000000

c Ad Ridder (VU)


3 / 49

c Ad Ridder (VU)

Getalrepresentatie

Getalrepresentatie heeft betrekking hoe getallen worden geschreven in de wiskunde en in een computertaal.

§1.3 Floating-Point Representatie

We onderscheiden 1. Klassieke voorstelling voor gehele getallen (type int) 2. Floating-point representatie voor reële getallen (type double) In Matlab zijn alle getallen (ook de integers!) van het type double, dus die behandelen we hier.

c Ad Ridder (VU)


5 / 49

Genormalizeerde Floating-Point Representatie

I

Ook wel wetenschappelijke notatie genoemd;

I

In wiskunde, elke x ∈

x = ± d0 .d1 d2 · · · × 10 = ±


6 / 49

Binaire Systeem

R (x 6= 0) kun je schrijven als: n

c Ad Ridder (VU)

∞ X

dk 10

−k

I

Op vorige slide had x decimale vorm;

I

Net zo goed kan x geschreven worden t.o.v. basis 2;

n

10 ,

x=±

k=0

|

{z

=r

∞ X

bk 2−k

2m ,

k=0

}

|

{z

=q

}

met cijfers met bk ∈ {0, 1}, b0 = 1 (dus 1 ≤ q < 2);

dk ∈ {0, 1, . . . , 9} zodat d0 6= 0;

I

I

r heet genormalizeerde mantisse (1 ≤ r < 10);

I

n∈

Z heet exponent.

c Ad Ridder (VU)


7 / 49

Dus genormalizeerde floating point representatie van reële getallen in het binaire systeem is x = ± (1.b1 b2 · · · )2 × 2m .

c Ad Ridder (VU)


8 / 49

Voorbeeld

Machinegetallen

(1.3625)10 × 10 = (13.625)10 = 1 × 101 + 3 × 100 + 6 × 10−1 + 2 × 10−2 + 5 × 10−3

I

Computers zijn eindige machines;

I

Getallen worden in binaire genormalizeerde floating point representatie opgeslagen met eindig veel (binaire) cijfers in de mantisse q, en eindig veel binaire cijfers voor exponent m;

I

De gerepresenteerde reële getallen heten machinegetallen;

I

Dat betekent bv dat

= 10 + 3 + 6/10 + 2/100 + 5/1000 = 8 + 4 + 1 + 1/2 + (1/10 + 1/50 + 1/200) = 8 + 4 + 1 + 1/2 + (20 + 4 + 1)/200 = 8 + 4 + 1 + 1/2 + 1/8 = 1 × 23 + 1 × 22 + 1 × 20 + 1 × 2−1 + 1 × 2−3 = (1101.101)2 = (1.101101)2 × 23

1 5

= (0.2)10 niet exact grepresenteerd zal kunnen worden;

NB: kan natuurlijk niet altijd (eigenlijk bijna nooit) in een eindige reeks! (0.2)10 = (1.10011001100110011 · · · )2 × 2−3

c Ad Ridder (VU)


9 / 49

Voorbeeld

c Ad Ridder (VU)


Standaard Dubbele Precisie Floating-Point

Stel dat een computer machinegetallen heeft waarvoor m

x = ± (1.b1 b2 )2 × 2 ,

I

In Matlab, Java en Ox zijn er veel meer (eindig veel) machinegetallen;

I

Namelijk x = ± (1. b1 b2 · · · b52 )2 × 2c−1023 , |{z} | {z }

met exponent m ∈ {−3, −2, . . . , 4}. I

10 / 49

1 bit

Merk op dat m = c − 3 en dat 3 bits nodig zijn voor c:

52 bits

met c ∈ 0, 1, . . . , (11111111111)2 = (2047)10 ; | {z }

c = (c2 c1 c0 )2 = c2 × 22 + c1 × 21 + c0 × 20 ∈ {0, 1, . . . , 7}.

11 bits

I

Is representatie uniek?

I

Wat zijn alle machinegetallen? Hoeveel?

I I

Kleinste positieve machinegetal: (1.00)2 ×

2−3

Grootste positieve machinegetal: (1.11)2 ×

24

Noteer f voor het fractionele deel van de mantisse;

I

Dus machinegetal x = ±q × 2m met q = 1 + f en

1 ; 8

f = (0.b1 b2 · · · b52 )2 = b1 × 2−1 + b2 × 2−2 + · · · + b52 × 2−52

= 28; =

I

Geen representatie voor 0 (!)

I

Teken de "getallenrechte" (vgl Figure 1.7);

c Ad Ridder (VU)

=

I

52 X

bi 2−k ∈

n

o 0, 2−52 , 2 × 2−52 , . . . , 1 − 2−52 × 2−52

k=1


m ∈ {−1023, −1022, . . . , 1024}

11 / 49

c Ad Ridder (VU)


12 / 49

Gereserveerde Getallen

Getallenrechte

Sommige machinegetallen zijn gereserveerd. I

I

f = (00 · · · 0)2 , c = 0: dus x = (1 + 0) × 2−1023 wordt gebruikt om 0 te representeren;

I

Krijg je als je een positief getal deelt door 0;

I

Of als je een te groot getal invoert; etc;

Idem wordt x = −(1 + 0) × 21024 voor -inf gebruikt;

I

f 6= (00 · · · 0)2 , c = 2047: dus x > (1 + 0) × 21024 wordt gebruikt als NaN (not a number); I

I

(w1 < w2 < · · · < wn );

wn = (1.11 · · · 1)2 × 22046−1023 = 1 + (1 − 2−52 ) 21023 ≈ 1.8 10308 ; I

Dus w1 = −wn ;

I

Het kleinste positieve getal heet realmin en is (1.00 · · · 0)2 × 21−1023 = (1 + 0) 2−1022 ≈ 2.2 10−308 ;

Concludeer dat exponenten c = 0 en c = 2047 niet gebruikt worden voor "gewone" getallen;


13 / 49

Opeenvolgende Machinegetallen

c Ad Ridder (VU)

I

Stel wi = q × 2m met q = (1.b1 b2 · · · b52 );

I

Voor het gemak, stel b52 = 0, dan is het volgende machinegetal wi+1 = (1.b1 b2 · · · b51 1) × 2m = q + 2−52 × 2m m−52

= wi + 2


Wat gebeurt er als je een getal invoert dat geen machinegetal, bv x = 0.2 of x = sqrt(2)?

Afrondafspraak

R

Elke x ∈ [w1 , wn ] ⊂ wordt afgerond naar het dichtstbijzijnde machinegetal, aangegeven door fl(x).

;

I

Dat geeft een afrondfout;

I

Stel x ∈ (wi , wi+1 ), dan voldoet de absolute fout van afronding aan

Het gat tussen twee opeenvolgende getallen die exact gerepresenteerd kunnen worden is 2m−52 ;

|fl(x) − x| ≤

Dat kan zo groot worden als 21022−52 = 2970 ≈ 10292 (!)

I

I


15 / 49

1 1 (wi+1 − wi ) = × 2m−52 = 2m−53 . 2 2

De relatieve fout is maximaal (voor ’t gemak wi > 0) |fl(x) − x| ≤ |x|

c Ad Ridder (VU)

14 / 49

Afronden I

I

R

Het grootste getal heet realmax en is

Krijg je bij 0/0, inf-inf, etc;

c Ad Ridder (VU)

I

De machinegetallen vormen een eindige verzameling W = {w1 , w2 , . . . , wn } ⊂

f = (00 · · · 0)2 , c = 2047: dus x = (1 + 0) × 21024 wordt gebruikt om inf te representeren;

I

I

I

Waarin u =

c Ad Ridder (VU)

2−53

1 (wi+1 2

wi

− wi )

≤

2m−53 = 2−53 = u, 2m

de eenheid van de afrondfout.


16 / 49

Machine Epsilon

Machine Epsilon (vervolg)

I

Het midden tussen wj en wj+1 is het getal h = 1 + 2−53 (geen machinegetal!);

I

Tel > 0 bij wj = 1 op: y = 1 + :

Definitie

I

Als 1 < y < h is y geen machinegetal en wordt in Matlab afgerond naar dichtstbijzijnde machinegetal, en dat is wj = 1;

I

Maw: als 0 < < 2−53 dan is fl(1 + ) = 1;

I

Als h ≤ y < wj+1 is y geen machinegetal en wordt in Matlab afgerond naar dichtstbijzijnde machinegetal, en dat is wj+1 6= 1;

I

Maw: als 2−53 ≤ < 2−52 dan is fl(1 + ) 6= 1;

m is het kleinste getal zodat fl(1 + m ) 6= 1; Berekening van m : I

Zij wj ∈ W; wj = 1; dus wj = (1 + 0) × 21023−1023 ;

I

Dwz in exponent m = 0 en mantisse q = 1;

I

Het volgende machinegetal is

Conclusie

Machine epsilon is m = 2−53 .

wj+1 = 1 + 2−52 × 20 = 1 + 2−52 ;

Gevolg In Matlab zijn machine epsilon en eenheid van afrondfout aan elkaar gelijk.

c Ad Ridder (VU)


17 / 49

c Ad Ridder (VU)


18 / 49

Nauwkeurigheid en Precisie

I

We hebben gezien voor machine epsilon m en de eenheid van afronding u:

§1.4 Loss of Significance

m = u = 2−53 ≈ 1.11 · 10−16 . I

Dus de floating-point representatie van doubles is in 15 decimalen nauwkeurig, dwz heeft 15 juiste (decimale) cijfers.

c Ad Ridder (VU)


19 / 49

c Ad Ridder (VU)


20 / 49

Effect van Floating Point Aritmetiek

Klassieke LoS

I

Bij verschil nemen van bijna evengrote getallen.

I

Voorbeeld.

I

Omdat getallen worden afgerond naar machinegetallen;

I

Dat gebeurt niet alleen als je getallen invoert;

I

Maar ook (en eigenlijk de meeste keren) bij elke berekening of operatie met die getallen;

I

Na een berekening kan de relatieve fout groter zijn dan de afrondeenheid;

I

Dat betekent dat het aantal correcte decimalen afneemt;

>> x = (1+0.6*eps)*(2^50); >> y = (1+0.4*eps)*(2^50); >> s = x-y

I

Bij elke volgende berekening kan dat propageren;

s =

I

NB: eps = 2−52 is het verschil in de mantisses van twee opeenvolgende machinegetallen.

0.2500

c Ad Ridder (VU)


21 / 49

Berekeningen

I

c Ad Ridder (VU)


22 / 49

Reductie van Loss of Significance

Exact zou moeten zijn 1 + 0.6 · 2−52 250 − 1 + 0.4 · 2−52 250 = 0.2 · 2−52 · 250 = 0.2 · 2−2 = 0.05.

I

Het berekende resultaat ontstaat doordat x naar boven wordt afgerond, en y naar beneden;

I

Matlab berekent

I

Soms mogelijk.

I

Voorbeeld 1: voor x ≈ 0 is sin(x) ≈ x.

I

Dus Matlab geeft x − sin(x) = 0 als x dichtbij 0 is.

I

Remedie: ontwikkel sin(x) in Taylorreeks, geeft x − sin(x) =

I

x3 x5 − + ··· . 3! 5!

Voor (heel) kleine x is x3 /3! een goede benadering.

fl fl(x) − fl(y) = fl 1 + 2−52 250 − 1 + 0 250 = fl 2−2 = 2−2 = 0.25; I

De relatieve fout is (0.25 − 0.05)/0.05 = 4 = 400%.

c Ad Ridder (VU)


23 / 49

c Ad Ridder (VU)


24 / 49

Reductie van Loss of Significance

I

Propogatie van Fout

Voorbeeld 2: abc-formule van nulpunt van parabool: √ −b ± b2 − 4ac x1,2 = . 2a

I

Beschouw x1 (met de +) en stel b2 4ac. √ Dan · ≈ b en fl(x1 ) = 0.

I

Remedie:

I

I

Integraalvoorbeeld;

I

Probleem: bereken voor n = 0, 1, . . .

−b − −b −

√

I

√

Met een beetje analyse kun je afleiden dat In ≥ 0 en

!

1

b2 − 4ac

I0 =

2c √ = . −b − b2 − 4ac

ex−1 dx = 1 − e−1 = 0.63212 · · · .

0 I

En voor n ≥ 1 pas PI toe. Geeft recursie 1 Z In = xn ex−1 − 0

c Ad Ridder (VU)


25 / 49

Matlabprogramma

c Ad Ridder (VU)

1

nxn−1 ex−1 dx = 1 − nIn−1 .

0


26 / 49

Resultaten

Matlab geeft

function integraalvb

I( 0) = 0.632120559 I( 1) = 0.367879441 I( 2) = 0.264241118 I( 3) = 0.207276647 -------------------

start; function start n = input(’geef n: ’); y = integraal(n); fprintf(’I(%2d) = %.9f\n’,n,y); end


I(16) I(17) I(18) I(19) I(20)

= = = = =

0.055459302 0.057191871 -0.029453671 1.559619744 -30.192394886

Grote afwijkingen, belachelijke uitkomsten. Reden: afrondfouten in I0 als In−1 een fout bevat, wordt die bij berekening van In versterkt door de vermenigvuldiging met n.

function y = integraal(n) y = 1-exp(-1); for k=1:n y = 1 - k*y; end; end end

c Ad Ridder (VU)

In ↓ 0.

Eenvoudig: Z

b2 − 4ac

xn ex−1 dx.

0

I

√ −b + b2 − 4ac x1 = 2a

1

Z In =

In dit geval is er een remedie: omdat limn→∞ In = 0 zet (bv) I100 = 0 en werk de recursie terug: 1 In−1 = (1 − In ). n

27 / 49

c Ad Ridder (VU)


28 / 49

Matrices: Invoeren

Wiskundige notatie: 

16 5  A= 9 4

MATLAB

3 10 6 15

2 11 7 14

 13 8  12 1

In Matlab voer je A als volgt in: A = [16 3 2 13; 5 10 11 8; 9 6 7 12; 4 15 14 1] I

Kolommen worden gescheiden door een spatie of een komma.

I

Rijen worden gescheiden door een puntkomma.

I

Het geheel omsluiten met vierkante haken.

I

Dus kan ook A=[[16;5;9;4], [3;10;6;15], [2;11;7;14], [13;8;12;1]]

c Ad Ridder (VU)


29 / 49

Speciale Vectoren

I

I

I

rijvector x die begint met first, optelt met 1, en eindigt op of voor last.

x = first:increment:last I

rijvector x die begint met first, optelt met increment, en eindigt voor of op last.

I

Nummering begint bij 1!

I

ij-de element aij krijg je door A(i,j).

I

Grotere delen door gebruik te maken dubbele punt.

x = linspace(first,last,n) I

rijvector x die met begint met first, en eindigt met last en n elementen heeft op gelijke afstand.

I

Voorbeelden: x = 1:8; x = 0.6:0.3:2.0; x = linspace(10,20,101);

c Ad Ridder (VU)


30 / 49

Matrices: Elementoperaties

x = first:last I

c Ad Ridder (VU)


31 / 49

I

Bv, 3-de rij door A(3,:).

I

3 × 3 submatrix van rijen 2,3,4 en kolommen 1,2,3 door A(2:4,1:3).

Merk op dat je vectoren kunt beschouwen als matrices, nl met één kolom of met één rij.

c Ad Ridder (VU)


32 / 49

Matrices: Invoegen en Verwijderen

I

Matrices: Enige Functies

Voeg een nieuwe 3-de kolom (van lengte 4) toe aan A, de oude 3-de en 4-de kolommen schuiven op:

I

Transponeren door apostrophe, A’ of A(3,:)’, etc.

I

Kolomsommen door sum(A) of sum(A,1).

I

Rijsommen door sum(A,2).

I

diag(A) geeft de kolomvector van de diagonaalelementen van A.

I

Maximaal element per kolom max(A) of max(A,[],1).

I

Maximaal element per rij max(A,[],2).

I

Maximaal element in de matrix max(max(A)).

A = [A(:,1:2), [1;2;-1;-2], A(:,3:4)] I

Verwijder kolom 4 (kolom 5 wordt dan de 4-de): A(:,4)=[]

c Ad Ridder (VU)


33 / 49

Enige Operaties met Matrices

c Ad Ridder (VU)


34 / 49

Matrices: Elementen Zoeken

I

Welke (i, j)-elementen in A zijn gelijk aan 1? [i,j] = find(A==1); disp([i,j]’) I

Stel A en B zijn 3 × 4-matrices.

Toelichting: find geeft het resultaat in 2 kolomvectoren.

I

Elementsgewijze optelling: A+B.

I

Elementsgewijze vermenigvuldiging: A.*B.

f = find(isprime(A)); p = A(f); disp([f,p]’)

I

Matrixvermenigvuldiging (lineaire algebra) A*B geeft foutmelding (waarom?).

Toelichting:

I

Matrixvermenigvuldiging C=A*B’ geeft een 3 × 3-matrix C (waarom?).

I I

I

Welke elementen zijn priem? En welke priemgetallen zijn die?

I

Elementen van een matrix kunnen rij/kolom genoteerd worden A(i,j).

I

Maar ook met één index A(k).

I

Dan wordt A beschouwd als een lange kolomvector gevormd uit de kolommen van A.

I

Dus A(2) is hetzelfde als A(2,1); en A(5) is hetzelfde als A(1,2) als A een matrix is met 4 rijen.

I

De find opdracht geeft een kolomvector f met de indices k waarvoor A(k) een priemgetal is.

I

A(f) zijn precies al die elementen, in een kolomvector.

Elementsgewijs machtheffen: A.^3. Machtheffen A^3 geeft foutmelding; C^3 is wel correct (waarom?).

c Ad Ridder (VU)


35 / 49

c Ad Ridder (VU)


36 / 49

Speciale Matrices

Plotten van een Grafiek 2

(x) .

I

Voorbeeld: f (x) = esin

I

Gevraagd de grafiek van f op [−2π, 2π] in een figuur te plotten.

I

Methode:

I

A = zeros(5,7) is een matrix van 5 rijen en 7 kolommen met 0-en.

1. Genereer voldoende veel x-waarden in [−2π, 2π] op gelijke afstand en inclusief de randen.

I

A = ones(5,7) idem met 1-en.

2. Bv door x = linspace(-2*pi,2*pi,200)

I

A = rand(5,7) idem met random uniform (0,1) getallen.

3. Bereken in deze x-waarden de bijbehorende functiewaarden: hulp = sin(x); y = exp(hulp.^2)

I

A = randn(5,7) idem met random normaal (0,1) getallen.

4. Je hebt nu 200 (xi , yi ) punten in de

R. 2

5. Verbind punten (x1 , y1 ) en (x2 , y2 ) met een rechte lijn. 6. Idem (x2 , y2 ) en (x3 , y3 ); etc tot en met (x199 , y199 ) en (x200 , y200 ). 7. En geef al die rechte lijn(tjes) weer in een grafiek (lijkt een differentieerbare kromme): plot(x,y)

Let op dat x en y twee vectoren moeten zijn met dezelfde orientatie en dezelfde lengte.

c Ad Ridder (VU)


37 / 49

Matlabprogramma

c Ad Ridder (VU)


38 / 49

Uitvoer

Je kunt aan de plot functie meer argumenten meegeven voor kleur, dikte, enz.

De grafiek wordt weergegeven in een figuur dat op je computerscherm komt.

function grafiekplotten start; function start n = input(’geef n: ’); x = linspace(-2*pi,2*pi,n); y = funf(x); plot(x,y,’k’,’LineWidth’,2); end function y = funf(x) hulp = sin(x); y = exp(hulp.^2); % puntsgewijs want x is een vector end end

c Ad Ridder (VU)


39 / 49

c Ad Ridder (VU)


40 / 49

Figuur Aanpassen

Plotopties Ook in je programma kun je opties opgeven. Hier is een code voor drie grafieken in verschillende kleuren, met toevoeging van titel, legenda, x en y-as benaming. Zie help plot.

Het figuurscherm heeft menu’s om de figuur te werken.

function grafieken start; function start x = linspace(0,10,301); y1 = sin(x); y2 = 3*exp(-x); y3 = sqrt(0.5*x); plot(x,y1,’k’,x,y2,’r’,x,y3,’b’); title(’drie grafieken’); legend(’fun1’, ’fun2’, ’fun3’); xlabel(’x’); ylabel(’y’); end end

c Ad Ridder (VU)


41 / 49

Uitvoer

c Ad Ridder (VU)


42 / 49

Meer Tips voor Opdracht 1 Stel p1 , . . . , pn zijn getallen 0 < pi < 1. Bereken dan log Voorbeeld:

Qn

i=1 pi

door

Pn

i=1

log pi .

p = 2.^(-(1:n)); x = log(prod(p)); y = sum(log(p)); Geeft n 10 20 30 40 50 60 70 80 90 100

c Ad Ridder (VU)


43 / 49

c Ad Ridder (VU)

Q log -38.12 -145.56 -322.31 -568.38 -Inf -Inf -Inf -Inf -Inf -Inf

P

log -38.12 -145.56 -322.31 -568.38 -883.76 -1268.46 -1722.47 -2245.80 -2838.44 -3500.39


44 / 49

Ook in Plot

Meer Tips voor Opdracht 1

Log likelihood functie log

Q

i

f (xi |λ, µ) versus

P

i

log f (xi |λ, µ) als functie van λ. Gegeven een functie van twee variabelen θ, x: f (θ, x) = xθ e−x

(θ, x > 0).

Bij een gegeven x-waarde, bv x = 4, wordt de functie een functie van alleen θ. Gevraagd de grafiek ervan te plotten op [0, 5]. x = 4; theta = linspace(0,5,200); L = x.^theta * exp(-x); plot(theta,L); Merk op: x.^theta is de vector

c Ad Ridder (VU)


45 / 49


c Ad Ridder (VU)

4θ1 , 4θ2 , . . . , 4θ200 .


46 / 49

Meer Tips voor Opdracht 1 Oplossing 1: theta = linspace(0,5,200); L = zeros(1,200); for i=1:200 L(i) = prod(x.^theta(i) .* exp(-x)); end plot(theta,L);

Stel x1 , . . . , x10 zijn 10 waarden van x (data), en je beschouwt het product als functie van θ: 10 Y L(θ) = xiθ e−xi (θ > 0). i=1

Schrijf een Matlab script om de grafiek van L(θ) te plotten op [0, 5]. θi Want: x.^theta(i) is de vector x1θi , x2θi , . . . , x10 , en exp(-x) is de vector −x −x −x e 1 , e 2 , . . . , e 10 . En dus is x.^theta(i) .* exp(-x) de vector

theta = linspace(0,5,200); L = prod(x.^theta .* exp(-x)); plot(theta,L);

Geeft foutmelding; want wat is

θi −x10 x1θi e−x1 , x2θi e−x2 , . . . , x10 e .

Hiervan het prod van nemen geeft (x1 , . . . , x10 )(θ1 ,...,θ200 ) ? Li =

10 Y

xjθi e−xj .

j=1

En dit afzonderlijk voor i = 1, 2, . . . , 200 (zie de for loop).

c Ad Ridder (VU)


47 / 49

c Ad Ridder (VU)


48 / 49


Oplossing 2: Voor elke θ is 10 Y j=1

xjθ e−xj =

10 Y θ P10 xj e− j=1 xj . j=1

Dus xp = prod(x); xs = sum(x); theta = linspace(0,5,200); L = xp.^theta * exp(-xs); plot(theta,L);

c Ad Ridder (VU)


49 / 49

Hoofdstuk 1. Illustratie 2

Recommend Documents