Chronology Current Month Current Thread Current Date [Year List] [Month List (current year)] [Date Index] [Thread Index] [Thread Prev] [Thread Next] [Date Prev] [Date Next]

# Re: [Phys-L] numerical methods

The textbook formula

(exp(-a*t1) - exp(-a*t2))/a

would run into significant precision subtraction error when |a*(t2 - t1)| gets small compared to 1, and also it divides by zero when a is zero. However I don't see any way out of the fundamental precision subtraction error that occurs when |t2 - t1| << 1 (& the original integral has a subtraction error in a single Riemann dt interval when the upper and lower integration limits are infinitesimally separated), for which the possibly errant dt is multiplied by an ostensibly finite integrand. So I think we have to live with any global precision subtraction error in the t2 - t1 value. Given that being the case, I suppose, if I was writing a fairly general routine to evaluate the above formula for a fairly unrestricted real domain of (a, t1, t2) and needed to avoid the precision subtraction error between the exponentials when they are nearly equal and avoid the division by zero when a is zero I would do something like the following:

Set threshhold = some appropriate positive number tiny compared to 1, but not too close to the smallest nonzero floating point number for the precision used.
Let dt = t2 - t1
if abs(x) < threshold
Set compsinc = 1st few initial terms of Taylor series for sinh(x)/x
else
Set compsinc = sinh(x)/x
endif
Set return value = dt* compsinc/exp(ad2*(t2 + t1))
end

If there was any possibility of values of the arguments of the intrinsic exp() and sinh() functions getting too big & cause an overflow, or the possibility a greatly negative argument of exp() causing an underflow I would put in further checks to handle those possible exceptional cases. Also, the above outline assumes that I wouldn't mind having the computer do whatever work was involved in evaluating both the intrinsic functions exp() *and* sinh(). If I was trying really hard to be miserly with the clock cycles I might try something else that didn't require the evaluation of 2 different transcendental functions. But the above outline treats both endpoints t1 & t2 on a similar footing without unduly biasing one end of the interval relative to the other.

Dave Bowman

JSD wrote:

Hi Folks --

This is relevant to measuring radioactive decays ... and to many other things. Very very many.

Consider the definite integral of exp(-a t) dt from t1 to t2.