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

[Iris Simmons]

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?
>
> 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]]