Roundoff error - double precision is not enough

coccoinomane

Hi everybody!

I kindly request your help. I have to compute functions like

[tex]\frac{ \sin (r x) - r x \cos (r x)}{r^3}[/tex]
(primitive function of x sin(rx) )

or

[tex]\frac{ -r x (120 - 20 r^2 x^2 + r^4 x^4) cos(r x) +
5 (24 - 12 r^2 x^2 + r^4 x^4) \sin(r x)}{r^7}[/tex]
(primitive function of x^5 sin(rx) )

when both r and x varies.

The problem with these functions is that the sin and cos factors are very similar to each other when x approaches zero. This is a big issue: double precision is not enough to compute the difference because of roundoff errors or simply because 14-15 digits are not enough to distinguish the two factors.

I kind of solved the problem for the first function. In fact, I could express it up to a factor as the first order spherical bessel function:

[tex]j1(x) = \frac{\sin(x)/x - \cos(x)}{x}[/tex],

which is well computed in the GNU Scientific Library.

I need to calculate those functions to solve the integral I discussed in http://www.physicsforums.com/showthr...08#post2028408 and that uart helped me to solve.

Thank you very much for any suggestion,

Guido

coccoinomane

I forgot to mention that I need to compute the functions for a very large array of "r" values from within a C++ program. Thus, pasting the result from Mathematica is not helpful :)

Cheers,

Guido

mathman

The problem you have is essentially in the numerators of your expressions when rx is small. I suggest that you (on paper) use the power series for sin and cos, where the terms in rx cancel for your expressions. Programming in the cubic term and possibly the fifth order (depending on how precise you want it) should be sufficient for small rx.

NoDoubts

or use any of the high precision c++ libraries...

coccoinomane

Hi mathman, thank you for your answer. You are right indeed.
I tried to expand the sine & cosine as my first approach, but I made a stupid mistake in the calculation of the coefficients and I got wrong results. I re-did everything and now I get results as precise as 10^-6. Thank you!

@NoDoubts
Could you please point me to some of these libraries?
However, I am not sure it would help. Unless they have computed the same function I need, i.e.

[tex]
\int_{0}^{k_0} k^{n} \frac{sin{kr}}{kr} dk
[/tex]

(and GSL doesn't have), I am afraid the precision problem would pop up again. However precise can the library be, it boils down to a difference between very small numbers --> loss of precision.

Cheers,

Guido