２つの周期の積

• 書きかけ

n<-4

t<-seq(from=0,to=1,length.out=1000)*2*pi*3



nt<-10000
dt<-0.001
p1<-3
p2<-2
a1<-1
a2<-0.4
delta<-runif(1)
delta<-1
k<-1
x<-matrix(0,nt,n)

X<-x
X[,1]<-a1*cos(p1*t)*(a2*cos(p2*(t+delta))+k)
X[,2]<-a1*sin(p1*t)*(a2*cos(p2*(t+delta))+k)
X[,3]<-a2*cos(p2*(t+delta))
X[,4]<-a2*sin(p2*(t+delta))


x<-matrix(0,nt,n)
x[1,]<-X[1,]
for(i in 2:nt){
	v<-rep(0,n)
	v[1]<--p1*x[i-1,2]-p2*x[i-1,1]*(x[i-1,4])/sqrt(x[i-1,1]^2+x[i-1,2]^2)
	v[2]<-p1*x[i-1,1]-p2*x[i-1,2]*(x[i-1,4])/sqrt(x[i-1,1]^2+x[i-1,2]^2)
	v[3]<--p2*x[i-1,4]
	v[4]<-p2*x[i-1,3]
	
	v<-v*dt
	x[i,]<-x[i-1,]+v
}

plot3d(x[,1],x[,2],x[,4])
plot(as.data.frame(x))

plot3d(X[,1],X[,2],X[,4])
plot3d(x[,1],x[,2],x[,4])
plot(as.data.frame(x))
matplot(x,type="l")

• 条件を変えて、複数本の軌道を描くソースはこちら

Nrep<-10
xssum<-NULL
col<-c()
for(rep in 1:Nrep){

n<-4

t<-seq(from=0,to=1,length.out=1000)*2*pi*3



nt<-1000
dt<-0.01
p1<-3
p2<-2
p2<-runif(1)*5
a1<-1
#a2<-runif(1)
a2<-0.4
delta<-runif(1)
delta<-1

k<-1
#k<-runif(1)*10
x<-matrix(0,nt,n)

X<-x
X[,1]<-a1*cos(p1*t)*(a2*cos(p2*(t+delta))+k)
X[,2]<-a1*sin(p1*t)*(a2*cos(p2*(t+delta))+k)
X[,3]<-a2*cos(p2*(t+delta))
X[,4]<-a2*sin(p2*(t+delta))


x<-matrix(0,nt,n)
x[1,]<-X[1,]
for(i in 2:nt){
	v<-rep(0,n)
	v[1]<--p1*x[i-1,2]-p2*x[i-1,1]*(x[i-1,4])/sqrt(x[i-1,1]^2+x[i-1,2]^2)
	v[2]<-p1*x[i-1,1]-p2*x[i-1,2]*(x[i-1,4])/sqrt(x[i-1,1]^2+x[i-1,2]^2)
	v[3]<--p2*x[i-1,4]
	v[4]<-p2*x[i-1,3]
	
	v<-v*dt
	x[i,]<-x[i-1,]+v
}

plot3d(x[,1],x[,2],x[,4])
plot(as.data.frame(x))

plot3d(X[,1],X[,2],X[,4])
plot3d(x[,1],x[,2],x[,4])
plot(as.data.frame(x))
matplot(x,type="l")
xssum<-rbind(xssum,x)
col<-c(col,rep(rep,nt))
}

plot3d(xssum[,1],xssum[,2],xssum[,4],col=col)