R: Integration: maximum number of units achieved, rounding error

Question

I ran into an interesting but rather unpleasant problem.

I am trying to integrate a function that was calculated from a dataset. Data can be found here: Link to sample.txt .

I start by setting a string to my data. this can be done linear with approxfun or nonlinear with splinefun . In my example below, I use the latter. Now when I try to integrate the installed function, I encounter an error

• maximum number of subdivisions reached

but when I increase the unit, I get

• roundoff error

From the values ​​in my code example, you can see that for this particular dataset the threshold is 754-> 755.

My colleague has no problem integrating this dataset into Matlab. Is there a way to manipulate my data for integration? Is there another method for numerical integration in R?

 data<-read.table('sample.txt',sep=',') colnames(data)<-c('wave','trans') plot(data$wave,data$trans,type='l') trans<- -1 * log(data$trans) plot(data$wave,trans,type='l') fx.spline<-splinefun(data$wave,trans) #Try either Fx.spline<-integrate(fx.spline,min(data$wave),max(data$wave)) #Above: Number of subdivision reached Fx.spline<-integrate(fx.spline,min(data$wave),max(data$wave),subdivisions=754) #Above: Number of subdivision reached Fx.spline<-integrate(fx.spline,min(data$wave),max(data$wave),subdivisions=755) #Above: Roundoff error 

Answer 1

There are many integration routines in R, and you can find some of them using RSiteSearch or the sos package.

For example, the pracma package has several implementations, for example

 quad(fx.spline,min(data$wave),max(data$wave)) # adaptive Simpson # [1] 2.170449 # 2.5 sec quadgk(fx.spline,min(data$wave),max(data$wave)) # adaptive Gauss-Kronrod # [1] 2.170449 # 0.9 sec quadl(fx.spline,min(data$wave),max(data$wave)) # adaptive Lobatto # [1] 2.170449 # 0.8 sec 

Please do not consider these to be pure R-scripts and therefore slower than, for example, a compiled integrate routine with such an oscillating function.