7

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
8
  • Consider posting the fitted function. I do not seem to be getting parameter estimates, although maybe that is not how splinefun works, or I am doing something wrong. Commented May 4, 2012 at 17:58
  • Your data looks spectacularly "flat," so maybe calling integrate and setting the conversion limit to something like 1e-9 via the rel.tol and abs.tol arguments will get you an answer that's plenty accurate. Commented May 4, 2012 at 18:01
  • @MarkMiller There is a good tutorial here: casoilresource.lawr.ucdavis.edu/drupal/node/896. To plot the function you can write plot(data$wave,fx.spline(data$wave),type='l')
    – Martin H
    Commented May 4, 2012 at 18:16
  • @CarlWitthoft I changed the tolerances to the lowest possible value but it does not change anything. The question is, why is Matlab doing it and R not :). I am aware that this curve is flat. It is a spectra of hydrogen chloride. Sadly I can't change chemistry
    – Martin H
    Commented May 4, 2012 at 18:25
  • I got this to work: Fx.spline <- integrate(fx.spline, 2550, 3100, subdivisions = 754) but so far not for the full range of x values. Commented May 4, 2012 at 18:50

1 Answer 1

7

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.

1
  • Good answer and a match to Matlab results. I am still learning R, but finding things in the documentations or searching for the correct terms is sometimes not that straight forward
    – Martin H
    Commented May 4, 2012 at 19:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.