[R] Numerical stability of: 1/(1 - cos(x)) - 2/x^2
Martin Maechler
m@ech|er @end|ng |rom @t@t@m@th@ethz@ch
Wed Aug 16 10:00:18 CEST 2023
>>>>> 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 <ikwsimmo using 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 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://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 using 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.
More information about the R-help
mailing list