Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

R: Integrate: Max number of subdivisions reached, roundoff error

Tags:

function

r

I came across an interesting but rather annoying problem.

I am trying to integrate a function which has been calculated from a dataset. The data can be found here: Link to sample.txt.

I start by fitting a line to my data. this can be done linear with approxfun or non-linear with splinefun. In my example below I use the latter. Now, when I try to integrate the fitted function I run into the error

  • maximum number of subdivisions reached

but when I increase the subdivision I get

  • roundoff error

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

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

enter image description here

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
like image 502
Martin H Avatar asked May 04 '12 16:05

Martin H


1 Answers

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

For example, package pracma has several implementations, for instance

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 not that these are pure R scripts and therefore slower than, e.g., the compiled integrate routine with such an oscillating function.

like image 155
Hans W. Avatar answered Sep 20 '22 02:09

Hans W.