Andrew Simmons
2023-Jan-08 04:02 UTC
[R] Problem with integrate(function(x) x^3 / sin(x), -pi/2, pi/2)
`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 at 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 at
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 at 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.