ich habe eine Beispieldatei, bei der die etwas weiter unten stehenden 3 Spalten von Bedeutung sind: 1. Zeit, 3. eine vorgegebene Sinuskurve (Sollkurve), und 2. eine Kurve, die die Sinuskurve nachzuahmen versucht (Istkurve). Ich möchte ab einem Startpunkt, den ich festgelegt habe, die Fläche unter der Kurve aus der 2. Spalte berechnen, konkret für eine Wellenlänge. Dazu nutze ich den auc-Befehl. Ich habe das auch für die "1. Wellenlänge" der Kurve aus Spalte 3 gemacht, da klappt alles. Bei Spalte 2 jedoch erscheint der Fehler "Error in integrate(myfunction, lower = from, upper = to) : maximum number of subdivisions reached".
Ich habe bereits gelesen, dass man die number of subdivisions erhöhen kann, allerdings bekomme ich es nicht hin, das auf dieses Beispiel zu übertragen. Über eine Hilfe wäre ich sehr dankbar. Den Beispieldatensatz_1 habe ich hoch geladen. (die Routine ist sicherlich kürzer bzw. einfacher möglich, bitte nicht wundern).
Vielen Dank im Voraus!
Dies ist der aktuelle Code:
Code: Alles auswählen
library(MESS)
require(MESS)
a<-list.files("C:/Users/schlg/Desktop/AUC_Test",
full.names=TRUE)
Name<-list.files("C:/Users/schlg/Desktop/AUC_Test",
full.names=FALSE)
for (i in 1:length(a)){
options(warn=1)
wo=scan(a[i], sep="\n", what=" ")
ab=grep("X_Value", wo)
länge=length(wo)
segment=länge-ab
b<-read.table(a[i], skip=ab, nrows=segment, dec=",", fill=TRUE)
Name1=strsplit(Name[i], "_")[[1]][2]
xWerte=b[[1]]
yIst=b[[2]]
ySoll=b[[3]]
trigger1soll=(which(ySoll>12.5)[1])-1 #determie start of sine
length=length(xWerte)-trigger1soll
lengthh=as.numeric(length)
Sollkurve=ySoll[trigger1soll:(trigger1soll+lengthh)]
Istkurve=yIst[trigger1soll:(trigger1soll+lengthh)]
xWerte_neu=xWerte[trigger1soll:(trigger1soll+lengthh)]
#horizontal graph to detect intersections with sine
x=rep(12.5, each=(length(Sollkurve)))
#determining intersection Points
above=Sollkurve>x
intersect.points<-which(diff(above)!=0)
Sollkurve.slopes<-Sollkurve[intersect.points+1]-
Sollkurve[intersect.points]
x.slopes<-x[intersect.points+1]-x[intersect.points]
x.points<-intersect.points + ((x[intersect.points] -
Sollkurve[intersect.points]) / (Sollkurve.slopes-x.slopes))
y.points<-Sollkurve[intersect.points] +
(Sollkurve.slopes*(x.points-intersect.points))
#Integral Istkurve 1. Wellenlänge
xx=x.points[c(TRUE, FALSE)]#Remove every other element to
#obtain full wave length (red circles)
yy=y.points[c(TRUE, FALSE)]#only necessary for plot
xx1=xx[1]
xx2=xx[2]
xWerte_1=xWerte_neu[xx1:xx2]
yIst_1=Istkurve[xx1:xx2]
#head(xWerte_1)#zeigt erste Daten an
#tail(xWerte_1)#zeigt letzte Daten an
Int_Ist_Total1=auc(xWerte_1,
yIst_1,
type = 'spline')
#Plot Gesamtdaten
plot(yIst[trigger1soll:(trigger1soll+lengthh)],
type="l", col="red", main=
"Kurvenverläufe", ylab='',
xlab="Time", xlim=c(0, 9000))
lines(Sollkurve, type="l", col="blue")
lines(x, type="l")
points(x.points, y.points, col="green")
points(xx, yy, col="red")
legend(6000, 2, c("Istkurve", "Sollkurve"),
fill=c("red", "blue"), cex=0.85)
abline(v=xx)
axis(side=1, at=seq(0, 8000, by = 1000))
}