Somewhat nasty function computation

[Herman Rubin]

Are there any good ways of handling the computation of
functions where for part of the region the "usual" formula
gives rise to large roundoff errors? The USER who is adept
at numerical analysis can go through the effort of writing
a procedure to do so, but is this not contrary to the
spirit of computer algebra systems?


--
This address is for information only. I do not claim that these views
are those of the Statistics Department or of Purdue University.
Herman Rubin, Department of Statistics, Purdue University
hru...@stat.purdue.edu Phone: (765)494-6054 FAX: (765)494-0558

[Richard Fateman]

On 3/21/2012 11:32 AM, Herman Rubin wrote:
> Are there any good ways of handling the computation of
> functions where for part of the region the "usual" formula
> gives rise to large roundoff errors? The USER who is adept
> at numerical analysis can go through the effort of writing
> a procedure to do so, but is this not contrary to the
> spirit of computer algebra systems?
>
>

google for kahan automated error
for a collection of links.

Short answer: it is hard to do right, and the people who claim to
do it (Mathematica, significance arithmetic) don't do it right.

You can use interval arithmetic, or you can just heuristically
try higher and higher precision and hope that when 2 answers are
the same, you have enough precision.

Or you can find someone who is adept.

Or find a library written by adept people.

RJF

[Herman Rubin]

I was suggesting something simpler than the google
search above. What I seem to run into are situations
where there are alternate methods of computing the
function over parts of the range.

Here are two examples. One is

2*csch(x)^2*(x*coth(x)-1)

and the other is

(1/x - 1/(2*sinh(x/2)))/x.

The second might be reasonably efficient for x sufficiently
far from 0, which might not need to be too far to be useful.
For x near 0, either the power series or a procedure to
calculate sinh(z)-z to precision accuracy would be a possibiliy.

For the first, I would not use the stated expression for
computation anywhere. For x large, express the function
in exp(-2*x), or even in z = 1/(exp(2*x) - 1), whence the
function to be calculated becomes z*(z+1)*(x*(1+2*z)-1).
but this is not good near 0. Here I would calculate
coth(x) - 1/x, either by series or by the excellent
continued fraction

x/3/(1+x^2/(3*5)/(1+x^2/(5*7)/(1+ ...

and use the identity csch(x)^2 = coth(x)^2-1.

You users of CAS, tell me how you would instruct the system,
if you can. And more, tell me how you would go about designing
the system so that the user would not need to supply that much
information. If we are going to automate it, that needs to be
done; I do not see that much automation at this time.