Ich habe eine numerische Funktion, f(x, y)
die eine doppelte Gleitkommazahl zurückgibt, die eine Formel implementiert, und ich möchte überprüfen, ob sie gegen analytische Ausdrücke für alle Kombinationen der Parameter korrekt ist x
und y
mich interessiert. Was ist der richtige Weg, um die berechneten und zu vergleichen? analytische Gleitkommazahlen?
Nehmen wir an, die beiden Zahlen sind a
und b
. Bisher habe ich sichergestellt, dass sowohl absolute ( abs(a-b) < eps
) als auch relative ( abs(a-b)/max(abs(a), abs(b)) < eps
) Fehler kleiner als eps sind. Auf diese Weise werden numerische Ungenauigkeiten erfasst, selbst wenn die Zahlen etwa 1e-20 betragen.
Heute habe ich jedoch ein Problem entdeckt, der numerische Wert a
und der analytische Wert b
waren:
In [47]: a
Out[47]: 5.9781943146790832e-322
In [48]: b
Out[48]: 6.0276008792632078e-322
In [50]: abs(a-b)
Out[50]: 4.9406564584124654e-324
In [52]: abs(a-b) / max(a, b)
Out[52]: 0.0081967213114754103
Der absolute Fehler [50] ist also (offensichtlich) klein, aber der relative Fehler [52] ist groß. Also dachte ich, dass ich einen Fehler in meinem Programm habe. Durch das Debuggen wurde mir klar, dass diese Zahlen denormal sind . Als solches habe ich die folgende Routine geschrieben, um den richtigen relativen Vergleich durchzuführen:
real(dp) elemental function rel_error(a, b) result(r)
real(dp), intent(in) :: a, b
real(dp) :: m, d
d = abs(a-b)
m = max(abs(a), abs(b))
if (d < tiny(1._dp)) then
r = 0
else
r = d / m
end if
end function
Wo tiny(1._dp)
gibt 2.22507385850720138E-308 auf meinem Computer zurück. Jetzt funktioniert alles und ich bekomme einfach 0 als relativen Fehler und alles ist in Ordnung. Insbesondere ist der obige relative Fehler [52] falsch, er wird einfach durch eine unzureichende Genauigkeit der denormalen Zahlen verursacht. Ist meine Implementierung der rel_error
Funktion korrekt? Sollte ich nur überprüfen, ob abs(a-b)
es weniger als winzig (= denormal) ist, und 0 zurückgeben? Oder sollte ich eine andere Kombination überprüfen, wie
max(abs(a), abs(b))
?
Ich möchte nur wissen, was der "richtige" Weg ist.
quelle
exp(log_gamma(m+0.5_dp) - (m+0.5_dp)*log(t)) / 2
für m = 234, t = 2000. Es geht schnell auf Null, wenn ich zunehmem
. Ich möchte nur sicherstellen, dass meine numerische Routine "korrekte" Zahlen (die Rückgabe von Null ist auch vollkommen in Ordnung) auf mindestens 12 signifikante Stellen zurückgibt. Wenn die Berechnung also eine Denormalzahl zurückgibt, ist sie einfach Null, und es sollte kein Problem geben. Daher muss nur die Vergleichsroutine robust sein.Donald Knuth hat einen Vorschlag für einen Gleitkomma-Vergleichsalgorithmus in Band 2 "Seminumerische Algorithmen" von "The Art of Computer Programming". Es wurde in C von Th implementiert. Belding (siehe fcmp-Paket ) und ist in der GSL verfügbar .
quelle
abs(a-b)/max(a, b) < eps
wir es tunabs(a-b)/2**exponent(max(a, b)) < eps
, anstatt es zu tun , was die Mantissemax(a, b)
so ziemlich einfach fallen lässt , also ist der Unterschied meiner Meinung nach vernachlässigbar.Optimal gerundete denormalisierte Zahlen können tatsächlich einen hohen relativen Fehler aufweisen. (Das auf Null zu setzen, während es immer noch als relativer Fehler bezeichnet wird, ist irreführend.)
Nahe Null ist die Berechnung relativer Fehler jedoch bedeutungslos.
Daher sollten Sie, noch bevor Sie denormalisierte Zahlen erreichen, wahrscheinlich auf absolute Genauigkeit umschalten (nämlich die, die Sie in diesem Fall garantieren möchten).
Dann wissen Benutzer Ihres Codes genau, wie genau sie wirklich sind.
quelle
abs(a-b) < tiny(1._dp)
wie ich es oben mache.tiny(1._dp)=2.22507385850720138E-308
(ich habe in meinem vorherigen Kommentar einen Fehler gemacht, es ist 2e-308, nicht 1e-320). Das ist also mein absoluter Fehler. Dann muss ich den relativen Fehler vergleichen. Ich verstehe deinen Standpunkt, ich denke du hast recht. Vielen Dank!