[R] Integrate and subdivisions limit

Hans W Borchers hwborchers at googlemail.com
Wed Jan 12 20:02:22 CET 2011


> Dear all,
> I have some issues with integrate in R thus I would like to request
> your help. I am trying to calculate the integral of f(x)*g(x).
> The f(x) is a step function while g(x) is a polynomial.
> If f(x) (step function) changes its value only few times (5 or 6 'steps')
> everything is calulated ok(verified results in scrap paper) but if f(x)
> takes like 800 different values I receive the error
> 
>     "Error in integrate number of subdivisions reached"
> 
> I did some checks on the internet and I found that I can increase the
> number of subdivisions (as this is a parameter in integrate().
> Thus I raised it from 100 to 1000 (and then to 10000).
> 
> A. Does this makes the error produced higher or does it only stress
>    the computer?
> B. When the number was raised to 10.000 I started getting the error message
>    "roundoff error was detected"
> 
> What do you think I should do to solve that?
> I would like to thank u in advance for your help
> 
> Best Regards
> Alex

There's obviously a more numerically stable approach. If g(x) is a polynomial
you do know its polynomial antiderivative. Take that and sum up all intervals
where the step function is constant.

Example:  g(x) = 1 constant, integrate x^2 over 0..1 in 10000 subdivisions.
          The antiderivative is 1/3 * x^3, thus

    g <- function(x) 1
    x <- seq(0, 1, len=10001)
    sum((x[2:10001]^3 - x[1:10000]^3)*g(2:10001))/3  #=> 0.3333333

The antiderivative of a polynomial a_1 x^n + ... + a_{n+1} given as vector
c(a1,...) can also be determined automatically, no manual work is necessary.

Hans Werner



More information about the R-help mailing list