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?
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
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.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With