ローレンツアトラクタ

計算 
# ローレンツ方程式
p<-10
r<-28
b<-8/3
f<-function(x,y,z){-p*x+p*y}
g<-function(x,y,z){-x*z+r*x-y}
h<-function(x,y,z){x*y-b*z}

# 平衡点を求める
# 初期の値
x<-10
y<-10
z<-10

T<-100
d<-0.01

for(t in 1:T){

A<-matrix(c(
(f(x+d,y,z)-f(x-d,y,z))/(2*d),(f(x,y+d,z)-f(x,y-d,z))/(2*d),(f(x,y,z+d)-f(x,y,z-d))/(2*d),
(g(x+d,y,z)-g(x-d,y,z))/(2*d),(g(x,y+d,z)-g(x,y-d,z))/(2*d),(g(x,y,z+d)-g(x,y,z-d))/(2*d),
(h(x+d,y,z)-h(x-d,y,z))/(2*d),(h(x,y+d,z)-h(x,y-d,z))/(2*d),(h(x,y,z+d)-h(x,y,z-d))/(2*d)
),byrow=T,3,3)

dX<- -solve(A)%*%c(f(x,y,z),g(x,y,z),h(x,y,z))
x<-x+dX[1]
y<-y+dX[2]
z<-z+dX[3]
}

x;y;z
f(x,y,z);g(x,y,z);h(x,y,z)  # 確認

# 軌跡を書いてみる
library(rgl)

x<-0.1;y<-0;z<-0
V<-matrix(c(x,y,z),3,1)
dt<-0.01
T<-100
Nt<-T/dt

for(t in 1:Nt){
dX<-c(f(x,y,z),g(x,y,z),h(x,y,z))*dt
x<-x+dX[1]
y<-y+dX[2]
z<-z+dX[3]
V<-cbind(V,c(x,y,z))
}

plot3d(t(V))

library(rgl)

# 平行移動と拡大縮小で(0,1)におさめる
U<-V
U[1,]<-(V[1,]-min(V[1,]))/(max(V[1,])-min(V[1,]))
U[2,]<-(V[2,]-min(V[2,]))/(max(V[2,])-min(V[2,]))
U[3,]<-(V[3,]-min(V[3,]))/(max(V[3,])-min(V[3,]))
plot3d(t(U),xlim=c(0,1),ylim=c(0,1),zlim=c(0,1))

# 精度を与えてデータを集計する
L<-40
M<-array(0,c(L,L,L))
for(l in 1:dim(U)[2]){
i<-floor(U[1,l]*L)
j<-floor(U[2,l]*L)
k<-floor(U[3,l]*L)
M[i,j,k]<-M[i,j,k]+1
}

plot3d(slice.index(M,1),slice.index(M,2),slice.index(M,3),alpha=ifelse(M==0,0,1))

# このままでは重いので点のあるところのみを選んでプロット
M2<-rbind((which(M > 0)-1) %% L,(which(M > 0)-1) %% L^2,(which(M > 0)-1) %% L^3)+1
plot3d(t(M2))

# 各位置での速度を計算する
Nx<-Ny<-Nz<-40

v1<-seq(from=min(V[1,],to=max(V[1,])),length=Nx)
v2<-seq(from=min(V[2,],to=max(V[2,])),length=Ny)
v3<-seq(from=min(V[3,],to=max(V[3,])),length=Nz)

# 各位置での速度を計算する
F<-array(0,c(Nx,Ny,Nz,3))
for(x in 1:Nx){
for(y in 1:Ny){
for(z in 1:Nz){
lx<-v1[x]
ly<-v2[y]
lz<-v3[z]
F[x,y,z,1]<-f(lx,ly,lz)
F[x,y,z,2]<-g(lx,ly,lz)
F[x,y,z,3]<-h(lx,ly,lz)
}}}

# 速度を色にしてプロットする
Col1<-(F[,,,1]-min(F[,,,1]))/(max(F[,,,1])-min(F[,,,1]))
Col2<-(F[,,,2]-min(F[,,,2]))/(max(F[,,,2])-min(F[,,,2]))
Col3<-(F[,,,3]-min(F[,,,3]))/(max(F[,,,3])-min(F[,,,3]))


plot3d(v1[slice.index(F,1)],v2[slice.index(F,2)],v3[slice.index(F,3)],xlab="x",ylab="y",zlab="z",col=rainbow(101)[Col1*100+1],main="velocity x")
plot3d(v1[slice.index(F,1)],v2[slice.index(F,2)],v3[slice.index(F,3)],xlab="x",ylab="y",zlab="z",col=rainbow(101)[Col2*100+1],main="velocity y")
plot3d(v1[slice.index(F,1)],v2[slice.index(F,2)],v3[slice.index(F,3)],xlab="x",ylab="y",zlab="z",col=rainbow(101)[Col3*100+1],main="velocity z")

# 偏微分係数を計算してプロット
# z方向について
Fz<-F[,,-1,3]-F[,,-Nz,3]

err<-0.000001
if(max(Fz) - min(Fz) < err){
ColFz<-1
}else{
ColFz<-(Fz-min(Fz))/(max(Fz)-min(Fz))
}

plot3d(v1[slice.index(Fz,1)],v2[slice.index(Fz,2)],v3[slice.index(Fz,3)],xlab="x",ylab="y",zlab="z",col=rainbow(101)[ColFz*100+1],main="Fz")

# 時系列データから速度を計算し色にしてプロット
du1<-U[1,-1]-U[1,-length(U[1,])]
du1<-(du1-min(du1))/(max(du1)-min(du1))

du2<-U[2,-1]-U[2,-length(U[2,])]
du2<-(du2-min(du2))/(max(du2)-min(du2))

du3<-U[3,-3]-U[3,-length(U[3,])]
du3<-(du3-min(du3))/(max(du3)-min(du3))

C1<-c(du1[1],du1*100+1)
C2<-c(du2[1],du2*100+1)
C3<-c(du3[1],du3*100+1)

plot3d(U[1,],U[2,],U[3,],col=rainbow(101)[C1])
plot3d(U[1,],U[2,],U[3,],col=rainbow(101)[C2])
plot3d(U[1,],U[2,],U[3,],col=rainbow(101)[C3])