R: Integrieren: Max. Anzahl der Unterteilungen erreicht, Rundungsfehler

Question

Ich bin auf ein interessantes aber eher nerviges Problem gestoßen.

Ich versuche eine Funktion zu integrieren, die aus einem Datensatz berechnet wurde. Die Daten finden Sie hier: Link to sample.txt.

Ich fange an, eine Zeile an meine Daten anzupassen. Dies kann linear mit approxfun oder nichtlinear mit splinefun erfolgen. In meinem Beispiel unten benutze ich Letzteres. Nun, wenn ich versuche, die angepasste Funktion zu integrieren, laufe ich in den Fehler

• maximum number of subdivisions reached

aber wenn ich die Unterteilung erhöhen ich

• roundoff error

Anhand der Werte in meinem Beispielcode können Sie das für diese spezifischen Daten sehen Stellen Sie den Schwellenwert auf 754-> 755.

Mein Kollege hat kein Problem, diesen Datensatz in Matlab zu integrieren. Gibt es eine Möglichkeit, meine Daten zu integrieren, um sie zu integrieren? Gibt es eine andere Methode zur numerischen 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

Es gibt viele Integrationsroutinen in R, und Sie können einige von ihnen von ‚RSiteSearch'ing oder mit dem‚sos‘Paket.

Zum Beispiel Paket pracma hat mehrere Implementierungen, beispielsweise

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 

Bitte nicht, daß diese reine R-Skripte und daher langsamer als beispielsweise die kompilierte integrate Routine mit einem solchen oszillierenden Funktion.