tmp<-read.table("tmp.txt")
x<-tmp[,1]
y<-tmp[,2]
f<-y
# find 1st derivative using forward difference
k<-length(x)
y1<-c()
for(i in c(1:k-1)) {
t<-(y[i+1]-y[i])/(x[i+1]-x[i])
y1<-c(y1, t)
}
y1<-c(y1,t)
f1<-y1
# find 2nd derivative using backward difference of 1st derivative
y<-f1
y1<-c()
for(i in c(2:k)) {
t<-(y[i]-y[i-1])/(x[i]-x[i-1])
y1<-c(y1, t)
}
y1<-c(y1[1],y1)
f2<-y1
# find curvature and plot
kappa<-f2/(1+f1*f1)^1.5
plot(x,kappa)
-------------------------------------------------------------------------
We can also use the backward difference to approximate the 1st derivative:
tmp<-read.table("tmp.txt")
x<-tmp[,1]
y<-tmp[,2]
f<-y
# find 1st derivative using backward difference
k<-length(x)
y1<-c()
for(i in c(2:k)) {
t<-(y[i]-y[i-1])/(x[i]-x[i-1])
y1<-c(y1, t)
}
y1<-c(y1[1],y1)
f1<-y1
# find 2nd derivative using forward difference of 1st derivative
y<-f1
y1<-c()
for(i in c(1:k-1)) {
t<-(y[i+1]-y[i])/(x[i+1]-x[i])
y1<-c(y1, t)
}
y1<-c(y1,t)
f2<-y1
# find curvature and plot
kappa<-f2/(1+f1*f1)^1.5
plot(x,kappa)
-------------------------------------------------------------------------
Finally, we can use central difference to approximate both the 1st and 2nd derivative. Here we assume the mesh grid is uniform.
tmp<-read.table("tmp.txt")
x<-tmp[,1]
y<-tmp[,2]
f<-y
# find 1st derivative using central difference
k<-length(x)
y1<-c()
for(i in c(2:k-1)) {
t<-(y[i+1]-y[i-1])/(x[i+1]-x[i-1])
y1<-c(y1, t)
}
y1<-c(y1,t)
y1<-c(y1[1],y1)
f1<-y1
# find 2nd derivative using central difference
y<-f1
y1<-c()
for(i in c(2:k-1)) {
t<-(y[i+1]-2*y[i]+y[i-1])/(x[i+1]-x[i])^2
y1<-c(y1, t)
}
y1<-c(y1,t)
y1<-c(y1[1],y1)
f2<-y1
# find curvature and plot
kappa<-f2/(1+f1*f1)^1.5
plot(x,kappa)
No comments:
Post a Comment