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@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@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]] ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
