Monday, March 28, 2011

How to plot curvature in R

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:

Visitors