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
integrate
and setting the conversion limit to something like 1e-9 via therel.tol
andabs.tol
arguments will get you an answer that's plenty accurate.plot(data$wave,fx.spline(data$wave),type='l')