Voorbeeld van een instabiel algoritme
.Definieer de integraal
En :=
Z 1
0
xnex−1dx voor n = 0, 1, 2, 3, · · ·.
De waarde van E0 is eenvoudig uit te rekenen, E0 =
Z 1
0
ex−1dx = ex−1
1 0
= 1− e−1 = 0.63212055882856 . Voor andere waarden van n leiden we via parti¨ele integratie de volgende recursie af:
En = Z 1
0
xnex−1dx = xnex−1
1
0 − n Z 1
0
xn−1ex−1dx = 1 − nEn−1. De voorwaartse recursie,
E0 := 1− e−1 , En := 1− nEn−1 (n = 1, 2, · · ·),
is instabiel, zoals blijkt uit de negatieve waarde voor n = 18 in de tabel hieronder. De reden is, dat een fout ε in Ek−1 een fout kε in Ek veroorzaakt. De fout in E18 is dus ongeveer 18! ≈ 1016 maal die in E0. De achterwaartse recursie,
Em = willekeurig , En−1 = (1− En)/n (n = m, m − 1, · · ·),
is stabiel ongeacht de startwaarde voor Em, als m maar voldoend groot wordt gekozen. Dit zien we in kolom 4, waar E18 = 0 gekozen is; hier neemt de fout bij iedere iteratie verder af en bij E5 is deze in de afrondfout verdwenen.
1 14/10/1993
n voorwaarts achterwaarts achterwaarts verschil tussen vanaf k = 0 vanaf n = 50 vanaf n = 18 kolommen 3 en 4 0 0.63212055882856 0.63212055882856 0.63212055882856 0.00000000000000 1 0.36787944117144 0.36787944117144 0.36787944117144 0.00000000000000 2 0.26424111765712 0.26424111765712 0.26424111765712 0.00000000000000 3 0.20727664702865 0.20727664702865 0.20727664702865 -0.00000000000000 4 0.17089341188538 0.17089341188538 0.17089341188538 0.00000000000000 5 0.14553294057308 0.14553294057308 0.14553294057308 -0.00000000000000 6 0.12680235656152 0.12680235656153 0.12680235656152 0.00000000000001 7 0.11238350406936 0.11238350406930 0.11238350406934 -0.00000000000004 8 0.10093196744509 0.10093196744559 0.10093196744528 0.00000000000032 9 0.09161229299417 0.09161229298966 0.09161229299250 -0.00000000000284 10 0.08387707005829 0.08387707010339 0.08387707007499 0.00000000002841 11 0.07735222935878 0.07735222886266 0.07735222917515 -0.00000000031248 12 0.07177324769464 0.07177325364803 0.07177324989825 0.00000000374978 13 0.06694777996972 0.06694770257562 0.06694775132275 -0.00000004874714 14 0.06273108042387 0.06273216394138 0.06273148148148 0.00000068245990 15 0.05903379364190 0.05901754087930 0.05902777777778 -0.00001023689848 16 0.05545930172957 0.05571934593124 0.05555555555556 0.00016379037568 17 0.05719187059731 0.05277111916899 0.05555555555556 -0.00278443638656 18 -0.02945367075154 0.05011985495809 0 0.05011985495809
2 14/10/1993
Voorbeeld 2, berekening van de variantie.
De variantie van een stel metingen kan berekend worden met twee Mathematisch equivalente formules. Gegeven n metingen {x1, x2, · · · , xn} van een grootheid X, dan is zijn gemiddelde g en variantie S2 gegeven door
g :=
Xn k=1
xk , Sn2 :=
Xn k=1
(xk − g)2 = Xn k=1
x2k − ng2 .
De tweede formule is potentieel numeriek instabiel (als Sn2 << g2).
Experiment. (met Matlab op een Mac)
>> rand(’normal’), n=10000;
>> x = rand(n,1) + 1e8 * ones(n,1);
>> g = sum(x)/n;
>> sig1 = x’*x - n*g*g;
>> sig2 = (x - g*ones(x))’ * (x - g*ones(x));
>> [sig2, sig1, sig2-sig1]
ans = 1.0027e+04 -1.1469e+05 1.2471e+05
3 14/10/1993
