>>>>> Iris Simmons >>>>> on Wed, 16 Aug 2023 02:57:48 -0400 writes:
> You might also be able to rewrite > log(1 - cos(x)) > as > log1p(-cos(x)) > On Wed, Aug 16, 2023, 02:51 Iris Simmons <ikwsi...@gmail.com> wrote: >> You could rewrite >> >> 1 - cos(x) >> >> as >> >> 2 * sin(x/2)^2 >> >> and that might give you more precision? Thank you, Iris, yes, both suggestions are excellent *and* necessary. Note that this is part "numerical analysis 101" and is called "cancellation": Direct evaluation of 1 - cos(x) for small 'x' *cannot* ever be numerically accurate and suffers from cancellation. log(1+x) for small x is slightly more subtle than pure cancellation, but exactly the same reason we introduced log1p() {and expm1(); then even cospi(), sinpi() etc} into R even before it was widely adopted in C math libraries. Martin >> 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://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. >>> >> > [[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. ______________________________________________ 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.