R help - Jan 2011 - Integrate and subdivisions limit

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

> 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