lbl-anp / becquerel

Becquerel is a Python package for analyzing nuclear spectroscopic measurements.
Other
43 stars 16 forks source link

Floating-point precision in exponential decays #65

Open bplimley opened 7 years ago

bplimley commented 7 years ago

IsotopeQuantity.decays_from and IsotopeQuantity.from_decays can suffer in precision if halflife >> time interval. This is because the calculations involve a 1 - np.exp(-decay_const * time) kind of expression, and if decay_const * time is very small, the exponent is very close to 1, and the difference loses precision because of these two floating point quantities. See np.spacing.

Better precision could be achieved in some cases by a Taylor series approximation:

1 - exp(-lambda * t)
~= 1 - (1 - lambda * t)
~= lambda * t

The code could use np.spacing to check whether the series approximation or floating point has better precision.

This is relevant for long-lived isotopes like K-40. (I'm getting ~0.1% loss of precision for an hour's interval with K-40, which was causing tests to fail.)

jvavrek commented 2 years ago

https://numpy.org/doc/stable/reference/generated/numpy.expm1.html