[R] Numerical stability of: 1/(1 - cos(x)) - 2/x^2

[Leonard Mada]

Dear Iris,
Dear Martin,

Thank you very much for your replies. I add a few comments.

1.) Correct formula
The formula in the Subject Title was correct. A small glitch swept into 
the last formula:
- 1/(cos(x) - 1) - 2/x^2
or
1/(1 - cos(x)) - 2/x^2 # as in the subject title;

2.) log1p
Actually, the log-part behaves much better. And when it fails, it fails 
completely (which is easy to spot!).

x = 1E-6
log(x) -log(1 - cos(x))/2
# 0.3465291

x = 1E-8
log(x) -log(1 - cos(x))/2
# Inf
log(x) - log1p(- cos(x))/2
# Inf => fails as well!
# although using only log1p(cos(x)) seems to do the trick;
log1p(cos(x)); log(2)/2;

3.) 1/(1 - cos(x)) - 2/x^2
It is possible to convert the formula to one which is numerically more 
stable. It is also possible to compute it manually, but it involves much 
more work and is also error prone:

(x^2 - 2 + 2*cos(x)) / (x^2 * (1 - cos(x)))
And applying L'Hospital:
(2*x - 2*sin(x)) / (2*x * (1 - cos(x)) + x^2*sin(x))
# and a 2nd & 3rd & 4th time
1/6

The big problem was that I did not expect it to fail for x = 1E-4. I 
thought it is more robust and works maybe until 1E-5.
x = 1E-5
2/x^2 - 2E+10
# -3.814697e-06

This is the reason why I believe that there is room for improvement.

Sincerely,

Leonard


On 8/16/2023 9:51 AM, Iris Simmons wrote:
> You could rewrite
>
> 1 - cos(x)
>
> as
>
> 2 * sin(x/2)^2
>
> and that might give you more precision?
>
> On Wed, Aug 16, 2023, 01:50 Leonard Mada via R-help 
> <r-help using r-project.org> wrote:
>
>     Dear R-Users,
>
>     I tried to compute the following limit:
>     x = 1E-3;
>     (-log(1 - cos(x)) - 1/(cos(x)-1)) / 2 - 1/(x^2) + log(x)
>     # 0.4299226
>     log(2)/2 + 1/12
>     # 0.4299069
>
>     However, the result diverges as x decreases:
>     x = 1E-4
>     (-log(1 - cos(x)) - 1/(cos(x)-1)) / 2 - 1/(x^2) + log(x)
>     # 0.9543207
>     # correct: 0.4299069
>
>     I expected the precision to remain good with x = 1E-4 or x = 1E-5.
>
>     This part blows up - probably some significant loss of precision of
>     cos(x) when x -> 0:
>     1/(cos(x) - 1) - 2/x^2
>
>     Maybe there is some room for improvement.
>
>     Sincerely,
>
>     Leonard
>     ==========
>     The limit was part of the integral:
>     up = pi/5;
>     integrate(function(x) 1 / sin(x)^3 - 1/x^3 - 1/(2*x), 0, up)
>     (log( (1 - cos(up)) / (1 + cos(up)) ) +
>          + 1/(cos(up) - 1) + 1/(cos(up) + 1) + 2*log(2) - 1/3) / 4 +
>          + (1/(2*up^2) - log(up)/2);
>
>     # see:
>     https://github.com/discoleo/R/blob/master/Math/Integrals.Trig.Fractions.Poly.R
>
>     ______________________________________________
>     R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
>     https://eu01.z.antigena.com/l/boSAjwics773HEe0HFHDZf3m1AU7fmWr4bglOgXO3sMyE9zLHAAMytf-SnATeHdnKJyeFbBsM6nXG-uPpd0NTW30ooAzNgYV5uhwnlhwxr4i8i21qKJUC~IrUTz2X1a5ioWqOWtHPlgzUrOid926sUOri-_H8XkLDcodDRWb
>
>     PLEASE do read the posting guide
>     https://eu01.z.antigena.com/l/AUS87vWM-isc3qtDXhJTp4jyQv7tuxdolKFlpY6mWcDOjbSlNzcD~~GORwHJFcX866fJF~qQmKc9R6LV9upRYcB4CBlSnLN0U_X8fIqLyhOIiPzDjYTVLEgiilZrKGuUqfW72b_50MVi~TaTlnE_R7fz8zXoZWGrKmGA
>
>     and provide commented, minimal, self-contained, reproducible code.
>
	[[alternative HTML version deleted]]