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 ad2 = a/2
Let dt = t2 - t1
Let x = ad2*dt
if abs(x) < threshold
Set compsinc = 1st few initial terms of Taylor series for sinh(x)/x
Set compsinc = sinh(x)/x
Set return value = dt* compsinc/exp(ad2*(t2 + t1))

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.
Here is the "textbook" answer:

Now explain why that is *not* the formula that I use in my software.

What is the smart way of doing it?

I'll post my answer tomorrow, if nobody nails it before then.

The usual jsd puzzle rules apply. Everything I've said is true to the best of my knowledge. Probably even helpful. No word games. There are undoubtedly more than one right answer, but I expect everybody will agree my answer is as good as any, and vastly better than the "textbook" answer.