[R] Problem with integrate(function(x) x^3 / sin(x), -pi/2, pi/2)

[Andrew Simmons]

`subdivisions` is the maximum number of subintervals. Looking here

https://github.com/wch/r-source/blob/79298c499218846d14500255efd622b5021c10ec/src/appl/integrate.c#L1275

I'm not surprised that changing `subdivisions` has no effect on the
outcome. The integration method from {pracma} might work, maybe it
never calculates the function at 0, maybe it's using some alternate
method. Or maybe it did calculate f(0) and got NaN and just assumed
that the limit converged. Maybe that's a fair assumption, and maybe
it's not.

However:

integrate(function(x) ifelse(x == 0, 0, x^3/sin(x), -pi/2, pi/2)

works perfectly fine and is a better function definition anyway.
`integrate` in the {stats} package is unlikely to change, so use the
alternate function definition or use {pracma}.

On Sat, Jan 7, 2023 at 10:41 PM Leonard Mada <leo.mada using syonic.eu> wrote:
>
> Dear Andrew,
>
>
> The limit when x approaches 0 is 0. I think most integration algorithms handle this situation.
>
>
> This actually works, and it "includes" formally the point x = 0:
>
> integrate(function(x) x^3 / sin(x), -pi/2, pi/2 + 1E-10)
>
>
> Trying to "bypass" the 0 using subdivisions unfortunately did not work:
>
> integrate(function(x) x^3 / sin(x), -pi/2, pi/2, subdivisions=4097) # or 4096
>
>
> Sincerely,
>
>
> Leonard
>
>
> On 1/8/2023 5:32 AM, Andrew Simmons wrote:
>
> You're dividing 0 by 0, giving you NaN, perhaps you should try
>
> function(x) ifelse(x == 0, 0, x^3/sin(x))
>
> On Sat, Jan 7, 2023, 22:24 Leonard Mada via R-help <r-help using r-project.org> wrote:
>>
>> Dear List-Members,
>>
>> I encounter a problem while trying to integrate the following function:
>>
>> integrate(function(x) x^3 / sin(x), -pi/2, pi/2)
>> # Error in integrate(function(x) x^3/sin(x), -pi/2, pi/2) :
>> #  non-finite function value
>>
>> # the value should be finite:
>> curve(x^3 / sin(x), -pi/2, pi/2)
>> integrate(function(x) x^3 / sin(x), -pi/2, 0)
>> integrate(function(x) x^3 / sin(x), 0, pi/2)
>> # but this does NOT work:
>> integrate(function(x) x^3 / sin(x), -pi/2, pi/2, subdivisions=4096)
>> integrate(function(x) x^3 / sin(x), -pi/2, pi/2, subdivisions=4097)
>> # works:
>> integrate(function(x) x^3 / sin(x), -pi/2, pi/2 + 1E-10)
>>
>>
>> # Note: works directly with other packages
>>
>> pracma::integral(function(x) x^3 / sin(x), -pi/2, pi/2 )
>> # 3.385985
>>
>>
>> I hope that integrate() gets improved in base R as well.
>>
>>
>> Sincerely,
>>
>>
>> Leonard
>>
>> ______________________________________________
>> 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.