[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