Octave Ungenauigkeit

[Jan K. (jan_k)]

Hi Leute,

bin gerade nach studenlangem Suchen nach dem Fehler meiner 2d DFT 
Implementation auf eine numerische Ungenauigkeit gestoßen, die alles 
kaputt gemacht hat.

Könnte von euch vielleicht mal jemand die Ausgaben von 
"exp(-i*2*pi*1*1./2)" in eurem Octave und bitte auch in Matlab (hab' ich 
nämlich nicht hier) posten? Ich bekomme sowas raus:

>>> exp(-i*2*pi*1*1./2)
ans = -1.00000000000000e+00 - 1.22460635382238e-16i

wie kann ich denn das am besten beheben?

[Stefan H. (Firma: dm2sh) (stefan_helmert)]

Welche Genauigkeit brauchst du denn mindestens?

[Marek N. (Gast)]

Gerne, Matlab:

  To get started, select MATLAB Help or Demos from the Help menu.

version

ans =

7.3.0.267 (R2006b)

>> exp(-i*2*pi*1*1./2)

ans =

  -1.0000 - 0.0000i

>>

[Purzel H. (hacky)]

Mein Maple ergibt exp(-i*pi). Was soll's denn sein ?

Allenfalls den Turbo reinwerfen, sodass noch ein paar mehr Stellen 
rauskommen ?

[Marek N. (Gast)]

Und Octave:

octave:3> version
ans = 3.0.0
octave:4> exp(-i*2*pi*1*1./2)
ans = -1.0000e+000 - 1.2246e-016i

[Jan K. (jan_k)]

Witzig witzig,

>>> version
ans = 3.4.3
>>> exp(-i*pi)
ans = -1.00000000000000e+00 - 1.22460635382238e-16i

Es soll natürlich -1 rauskommen, meine benötigte Genauigkeit ist 
allerallerhöchstens E-6. Ich ging dennoch davon aus, dass Octave das 
berechnen kann - abgesehen von der Genauigkeit. Also nicht rein 
numerisch - hab ich mich offensichtlich vertan :D

Danke fürs Feedback

[Marek N. (Gast)]

Scilab macht auch Murks:

-->ver
 ans  =
 
!Scilab Version:             5.3.3.1311149487                                         !
!                                                                                     !
!Operating System:           Windows XP 5.1                                           !
!                                                                                     !
!Java version :              1.6.0_18                                                 !
!                                                                                     !
!Java runtime information :  Java(TM) SE Runtime Environment (build 1.6.0_18-b07)     !
!                                                                                     !
!Java vm information :       Java HotSpot(TM) Client VM (build 16.0-b13, mixed mode)  !
!                                                                                     !
!Vendor specification:       Sun Microsystems Inc.                                    !
 
-->exp(-%i*2*%pi*1*1./2)
 ans  =
 
  - 1. - 1.225D-16i

[Florian R. (Firma: TU Wien) (frist)]

...und was bekommt ihr so für exp(pi)-pi raus?

SCNR
Flo

[Jan K. (jan_k)]

?? Wofür war das jetzt?

Aber mal ernsthaft, ich kann doch wohl von einem Mathematik Programm 
ausgehen, dass es sowas rechnen kann :

octave:16> cos(pi)
ans = -1
octave:17> sin(pi)
ans =  1.22460635382238e-16
octave:18> cos(pi/2)
ans =  6.12303176911189e-17

Oder etwa nicht?

[Florian R. (Firma: TU Wien) (frist)]

Hallo Jan

> ?? Wofür war das jetzt?

Ich musste nur an das hier denken: http://xkcd.com/217/


Ansonsten hab ich nicht den Eindruck als würde irgendwer falsch rechnen.

In allen geposteten Beispielen wurde scheinbar mit ca. 16 (dezimalen) 
Stellen gerechnet, das wären dann Zahlen vom Type double. Da ist es doch 
kein wunder, dass nicht genau 0, sondern eben fast null, ca. 10^-16 
auskommt. Auch das Pi, mit dem die Rechunung als Eingabe begann, war ja 
nur auf ca. +/- 10^-16 genau angegeben, als double eben. Wenn man's 
genau wissen wollte warum und wieso müsste man sich eh die binäre 
Darstellung der Zahlen anschauen, z.B. hier: 
http://babbage.cs.qc.edu/IEEE-754/

In Mathematik 8 sie die Sache übrigens so aus:

In[39]:= $MachinePrecision
Out[39]= 15.9546

In[40]:= NumberForm[Evaluate[Exp[-I*2*Pi*1*1./2]], 20]
Out[40]=-1. -(1.2246467991473532x10^1-16)i

In[41]:= $MachineEpsilon 
Out[41]= 2.22045*10^-16

In[42]:= Chop[Exp[-I*2*Pi*1*1./2]]
Out[42]= -1.

Wie man sieht, mit 64bit floating points gerechnet, das gibt ca. 16 
Dezimalstellen. $MachineEpsilon ist übrigens die kleinste (in der 
gewählten Genauigkeit) darstellbare Differenz zu 1. Der Fehler ist 
kleiner, also alles OK, find ich.

Für das Problem, dass bei numerischer Berechnung z.B. ein Imaginärteil 
nicht wirklich weg fällt gibt's übrigens zumindest in Mathematica auch 
eine Lösung, siehe letzten beiden Zeilen.


Die oben gepostete Matlab Ausgabe ist übrigens wohl irreführend, weil 
die Zahlen gerundet ausgegeben wurden. Ich weiß grad nicht wie man die 
Ausgabe aller Stellen erzwingt, vielleicht mit "format long".

Grüße
Flo

[Florian R. (Firma: TU Wien) (frist)]

... noch ein kleiner Nachtrag:

Der Fehler mit dem Pi`15.9546 (also Pi als 64bit float) behaftet ist ist 
übrigens ca. 1.2246467991473531772e-16. Sprich das ganze rumgerechne hat 
den Eingangsfehler nicht größer gemacht - was will man sich mehr 
wünschen.

Grüße
Flo

[Jan K. (jan_k)]

Hey, danke für deine Ausführungen. Den Comic kannte ich noch nicht ;)

Bedeutet das denn nun, dass in den Matlab Ausgaben oben automatisch 
gerundet wird, sobald das Ergebnis kleiner ist als die Genauigkeit und 
bei octave eben nicht, oder wie kann man das deuten?

Ich meine, es ist ja schön, dass numerisch richtig gerechnet wird, ich 
als Benutzer kann aber mit den Zahlen, solange sie ungerundet bleiben 
nix anfangen.

Und ich habe gelernt, dass floor() bei sowas ganz böse ist, da wird aus 
-x*E-16 ganz schnell -1 :D

[Marek N. (Gast)]

Viel schlimmer ist noch, dass exp(pi) ungefähr 23 ist! Ich vermute da 
die Illuminaten dahinter!

[j. c. (jesuschristus)]

Matlabimaginärteil:
-1.2246e-016i

Genau das gleiche. Wenn sowas aber zu einem Fehler führt, liegt das 
Problem woanders.

[Jan K. (jan_k)]

Wie gesagt, der eigentliche Fehler lag darin, dass ich zur besseren 
Übersicht in der Konsole die Matrix Elemente (falsch) gerundet habe, 
ohne mir vorher das echte Ergebnis angeguckt zu haben.

[Timm T. (Gast)]

Wolfram Alpha sagt: -1

http://www.wolframalpha.com/input/?i=exp%28-i*2*pi*1*1.%2F2%29

[j. c. (jesuschristus)]

Dann überleg Dir mal den Unterschied zwischen Octave und Wolfram Alpha.

[Daniel -. (root)]

noch ein Tipp:

Python
>>> "%.50f" % (0.1*3)
'0.30000000000000004440892098500626161694526672363281'

# perl -e 'printf("%.50f\n", 0.1*3)'
0.30000000000000004440892098500626161694526672363281

es gibt eben die numerische Mathematiksoftware und die symbolisch 
evaluirende wie Maxima