ローレンツアトラクタ - ryamadaのコンピュータ・数学メモ

ローレンツアトラクタも描いておこう(アトラクタ)
$x,y,z$ の時間微分をわかりやすいように式変形すると以下のようになる
$\frac{dx}{dt}=-10x+10y=-10(x-y)$
$\frac{dy}{dt}=28x-y-xz=(x-y)-(x-0)(z-27)$
- $\frac{dz}{dt}=-\frac{8}{3}(z-27) +(xy-(6\sqrt{2})^2)$
$\frac{dz}{dt}=-\frac{8}{3}(z-27) +(x\pm 6\sqrt{2})(y\mp 6\sqrt{2}) \pm 6\sqrt{2}(x-y)$
$x=y$ , $z=27$ , $x=0$ , $xy=72$ という面が変数の時間微分の正負に影響を与えることがわかる
また、自明な固定点 $x=y=z=0$ のほかに $x=y=\pm6\sqrt{2},z=27$ という固定点があることもわかる
この固定点の周囲では $x,y,z$ のそれぞれの増減が $(+,+,+),(+,+,-),(+,-,+),...,(-,-,-)$ の８パターン(３次元市松模様)を取っていることを確かめるソースもつける
色は黒が $(\frac{dx}{dt},\frac{dy}{dt},\frac{dz}{dt})$ について(負、負、負)、赤が（正、負、負）、緑が（負、正、負）、青が（正、正、負）、水色が（負、負、正）、紫が（正、負、正）、黄色が（負、正、正）、灰色が（正、正、正）

$z=27$ 平面では $\frac{dx}{dt}=-10(x-y)$ なので、 $x=y$ 直線でx軸方向動きはない
$\frac{dy}{dt}=28(x-y)-x(z-27)=x-y$ なので、y軸方向も $x=y$ で動きがない
$\frac{dz}{dt}=xy-72$ なので、xyに関して双曲線の内側か外側かで向きが変わる。
つまり、 $y=x$ 直線上では、双曲線の外側ではz軸正に内側ではz軸負の動きをする。
これがローレンツアトラクタの円板に見える部分（z軸方向のみの動きでx,y軸への動きがないため平らになる
固定点からずれると、動きのない面が傾くので、軌道を図示したときに見える円盤は斜めになっている
zがより大きくなると $x=y$ 面上では相変わらずx軸方向の動きはないが、y軸の動きは負となり、ある高さでzの動きがなくなるので、そこではy軸負方向のみの動きになる
同様にzを小さくするとy軸正方向のみの動きになる点もある
似たように考えてx軸正のみ、x軸負のみの動きの点も固定点の周りにあり、それが固定点のごく近傍では平面上に並び、固定点から遠くなると傾きが大きくなってきたのがアトラクタの絵
この「平らな面」が２個の非原点固定点のまわりで異なるのだが、その折り合いをつけた部分が、２つの円盤形の移行部分

library(MCMCpack)
xyz1<-c(6*sqrt(2),6*sqrt(2),27) # 原点ではない固定点を中心に
xyz1<-rep(0,3) # 原点
#xyz1<-runif(3)*5 # 適当な固定点
d<-20 # 調べる範囲を決めるパラメタ

#pmpmpm<-expand.grid(rep(list(c(-1,1)),3))
Niter<-100000

loc<-out<-matrix(0,Niter,3)
for(i in 1:Niter){
	loc[i,]<-xyz1+rnorm(3)*d
	#loc[i,]<-xyz1+rdirichlet(1,rep(1,3))*d
	#loc[i,]<-xyz1+runif(3)*d
	tmp<-loc[i,]
	out[i,]<-sign(c(-10*(tmp[1]-tmp[2]),(tmp[1]-tmp[2])-tmp[1]*(tmp[3]-27),-8/3*(tmp[3]-27)+tmp[1]*tmp[2]-72))
}
out<-(out+1)/2
col<-t(t(out)*c(1,2,4))
col2<-apply(col,1,sum)+1
plot3d(loc[,1],loc[,2],loc[,3],col=col2)

#X <- par3d("userMatrix") 
xlim<-ylim<-zlim<-c(min(loc),max(loc))
for(i in 1:8){
	open3d()
	coltobeplotted<-i
	#filename<-paste("out",i,sep="")
plot3d(loc[which(col2==coltobeplotted),1],loc[which(col2==coltobeplotted),2],loc[which(col2==coltobeplotted),3],col=coltobeplotted,xlim=xlim,ylim=ylim,zlim=zlim)
	#movie3d( par3dinterp( userMatrix=list(X,rotate3d(X, pi/2, 1, 0, 0),rotate3d(X, pi/2, 0, 1, 0) )), duration=5 ,movie=filename,dir=".")

}

原点周囲はどうなっているだろうか。ローレンツアトラクタ軌道のよじれは原点周辺で起きている
原点周辺のベクトル正負を８パターン別にプロットすると単純な８分割になっていないことがわかる

空間のいろいろなところからすべてこの２重渦巻に取り込まれていく様子を、いろいろな点からスタートさせた多数の軌道によって示したのが以下の図
すべての始点は黒で、それ以降は軌道ごとに色を変えてある。
$z=0$ 平面では、２方向に向かって吹き別れていく様子が見える

Nrep<-400
xssum<-c()
col<-c()
for(rep in 1:Nrep){
	Niter<-40
	#col<-c(col,rep(rep,Niter))
	col<-c(col,c(1,rep(rep,Niter-1)))
dt<-0.002
xs<-matrix(0,Niter,3)

xs[1,]<-runif(3,min=-20,max=20)
xs[1,3]<-abs(xs[1,3])*2
#xs[1,]<-c(2*sqrt(2)*3,2*sqrt(2)*3,27)
vs<-runif(3)*50
vs<-rep(0,3)
for(i in 2:Niter){
	xs[i,1]<-xs[i-1,1]+(-10*xs[i-1,1]+10*xs[i-1,2]+vs[1])*dt
	xs[i,2]<-xs[i-1,2]+(28*xs[i-1,1]-xs[i-1,2]-xs[i-1,1]*xs[i-1,3]+vs[2])*dt
	xs[i,3]<-xs[i-1,3]+(-8/3*xs[i-1,3]+xs[i-1,1]*xs[i-1,2]+vs[3])*dt
}
#tmpx1<-xs[,1]/sqrt(2)-xs[,2]/sqrt(2)
#tmpy1<-xs[,1]/sqrt(2)+xs[,2]/sqrt(2)
#xs[,1]<-tmpx1
#xs[,2]<-tmpy1
xssum<-rbind(xssum,xs)
}

library(rgl)
open3d()
xlim<-ylim<-zlim<-c(min(xssum),max(xssum))
#plot3d(xssum[,1],xssum[,2],xssum[,3],col=rainbow(1000),cex=0.1)
#plot3d(xssum[,1],xssum[,2],xssum[,3],col=gray(col/Nrep),cex=0.1,xlim=xlim,ylim=ylim,zlim=zlim)
plot3d(xssum[,1],xssum[,2],xssum[,3],col=col,xlim=xlim,ylim=ylim,zlim=zlim,cex=0.1)

この軌道の通り方は、x,y,zの時間微分の式 $\frac{dx}{dt}=0,\frac{dy}{dt}=0,\frac{dz}{dt}=0$ の面の重なり具合が決める
この面は一つはxy平面に垂直な面、一つは $y=(28-z)x$ が表す曲面で、これは、zの値が与えられたときに $y=Cx$ という直線だから、原点を通る直線がzの値に応じて変化しながら作る曲面である。また、この曲面は $z=27$ のときにxy平面に垂直な面と直線で交わる。最後の面は双曲線が移行しながら作る面

library(lattice)
library(rgl)

maxz<-100

n<-100
x<-seq(from=-30,to=30,length.out=n)
y<-seq(from=-30,to=30,length.out=n)
xy2<-expand.grid(x,y)

z3<-xy2[,1]*xy2[,2]*3/8
tobeplotted3<-which(abs(z3)<100)

n<-400

x<-seq(from=-30,to=30,length.out=n)
z<-seq(from=-maxz,to=maxz,length.out=n)
xz<-expand.grid(x,z)

y2<-28*xz[,1]-xz[,1]*xz[,2]

tobeplotted2<-which(abs(y2)<100)

n<-50
y<-seq(from=-30,to=30,length.out=n)
z<-seq(from=-maxz,to=maxz,length.out=n)
yz<-expand.grid(y,z)
x1<-yz[,1]

xx<-c(x1,xz[tobeplotted2,1],xy2[tobeplotted3,1])
yy<-c(yz[,1],y2[tobeplotted2],xy2[tobeplotted3,2])
zz<-c(yz[,2],xz[tobeplotted2,2],z3[tobeplotted3])


gr<-c(rep(1,length(x1)),rep(2,length(y2[tobeplotted2])),rep(3,length(z3[tobeplotted3])))

plot3d(xx,yy,zz,col=gr)

Niter<-100000
dt<-0.001
xs<-matrix(0,Niter,3)

xs[1,]<-runif(3)

for(i in 2:Niter){
	xs[i,1]<-xs[i-1,1]+(-10*xs[i-1,1]+10*xs[i-1,2])*dt
	xs[i,2]<-xs[i-1,2]+(28*xs[i-1,1]-xs[i-1,2]-xs[i-1,1]*xs[i-1,3])*dt
	xs[i,3]<-xs[i-1,3]+(-8/3*xs[i-1,3]+xs[i-1,1]*xs[i-1,2])*dt
}

library(rgl)

plot3d(xs[,1],xs[,2],xs[,3],col=rainbow(1000),cex=0.1)

係数が、-10,8/3,28など風変りなので、上記の式変形を参考に次のような係数にしてソースを書き直してみる
$\frac{dx}{dt}=-10x+10y=-Ks(1)\times (x-Ks(11)\times y)$
$\frac{dy}{dt}=28x-y-xz=Ks(2)(x-Ks(11)y)-Ks(4)(x-Ks(3))(z-Ks(5))$
$\frac{dz}{dt}=-Ks(6)\times (z-Ks(4)) +Ks(10)\times xy-Ks(4)\times Ks(6)$

Niter<-10000
dt<-0.01
xs<-matrix(0,Niter,3)

xs[1,]<-runif(3)*3
#xs[1,]<-c(2*sqrt(2)*3,2*sqrt(2)*3,27)
vs<-runif(3)
vs<-rep(0,3)
Ks<-runif(11)
Ks[4]<-Ks[4]*27
Ks[6]<-Ks[6]*8/3
Ks[11]<-Ks[11]*2
Ks[2]<-Ks[2]*2
Ks[1]<-Ks[1]*10
Ks[3]<-Ks[3]*2
Ks[5]<-Ks[5]*2
#Ks[11]<-1
#Ks[2]<-1
#Ks[1]<-10
#Ks[5]<-1
#Ks[3]<-0
#Ks[4]<-27
#Ks[6]<-8/3
Ks[10]<-1
Ks[7]<-6*sqrt(2)
Ks[9]<-6*sqrt(2)
for(i in 2:Niter){
	xs[i,1]<-xs[i-1,1]+(-Ks[1]*(xs[i-1,1]-Ks[11]*xs[i-1,2]))*dt
	xs[i,2]<-xs[i-1,2]+(Ks[2]*(xs[i-1,1]-Ks[11]*xs[i-1,2])-Ks[5]*(xs[i-1,1]-Ks[3])*(xs[i-1,3]-Ks[4]))*dt
	xs[i,3]<-xs[i-1,3]+(-Ks[6]*(xs[i-1,3]-Ks[4])+Ks[10]*(xs[i-1,1]*xs[i-1,2])-Ks[6]*Ks[4])*dt
	#xs[i,1]<-xs[i-1,1]+(-10*xs[i-1,1]+10*xs[i-1,2]+vs[1])*dt
	#xs[i,2]<-xs[i-1,2]+(28*xs[i-1,1]-xs[i-1,2]-xs[i-1,1]*xs[i-1,3]+vs[2])*dt
	#xs[i,3]<-xs[i-1,3]+(-8/3*xs[i-1,3]+xs[i-1,1]*xs[i-1,2]+vs[3])*dt
}
#tmpx1<-xs[,1]/sqrt(2)-xs[,2]/sqrt(2)
#tmpy1<-xs[,1]/sqrt(2)+xs[,2]/sqrt(2)
#xs[,1]<-tmpx1
#xs[,2]<-tmpy1
library(rgl)
open3d()
plot3d(xs[,1],xs[,2],xs[,3],col=rainbow(1000),cex=0.1)

軌跡とx,y,zの微分が0になる面とを重ねて描いてみよう

Nrep<-5
ct<-c(6*sqrt(2),6*sqrt(2),27)
xssum<-c()
col<-c()
for(rep in 1:Nrep){
	Niter<-1000
	#col<-c(col,rep(rep,Niter))
	col<-c(col,c(1,rep(rep,Niter-1)))
dt<-0.004
xs<-matrix(0,Niter,3)

xs[1,]<-runif(3,min=-10,max=10)+ct
#xs[1,3]<-abs(xs[1,3])*2

#xs[1,]<-c(2*sqrt(2)*3,2*sqrt(2)*3,27)
vs<-runif(3)*50
vs<-rep(0,3)
for(i in 2:Niter){
	xs[i,1]<-xs[i-1,1]+(-10*xs[i-1,1]+10*xs[i-1,2]+vs[1])*dt
	xs[i,2]<-xs[i-1,2]+(28*xs[i-1,1]-xs[i-1,2]-xs[i-1,1]*xs[i-1,3]+vs[2])*dt
	xs[i,3]<-xs[i-1,3]+(-8/3*xs[i-1,3]+xs[i-1,1]*xs[i-1,2]+vs[3])*dt
}
#tmpx1<-xs[,1]/sqrt(2)-xs[,2]/sqrt(2)
#tmpy1<-xs[,1]/sqrt(2)+xs[,2]/sqrt(2)
#xs[,1]<-tmpx1
#xs[,2]<-tmpy1
xssum<-rbind(xssum,xs)
}



library(rgl)
open3d()
xlim<-ylim<-zlim<-c(min(xssum),max(xssum))
#plot3d(xssum[,1],xssum[,2],xssum[,3],col=rainbow(1000),cex=0.1)
#plot3d(xssum[,1],xssum[,2],xssum[,3],col=gray(col/Nrep),cex=0.1,xlim=xlim,ylim=ylim,zlim=zlim)
plot3d(xssum[,1],xssum[,2],xssum[,3],col=col,xlim=xlim,ylim=ylim,zlim=zlim,cex=0.1)

library(lattice)
library(rgl)

maxz<-100

n<-100
x<-seq(from=-30,to=30,length.out=n)
y<-seq(from=-30,to=30,length.out=n)
xy2<-expand.grid(x,y)

z3<-xy2[,1]*xy2[,2]*3/8
tobeplotted3<-which(abs(z3)<100)

n<-400

x<-seq(from=-30,to=30,length.out=n)
z<-seq(from=-maxz,to=maxz,length.out=n)
xz<-expand.grid(x,z)

y2<-28*xz[,1]-xz[,1]*xz[,2]

tobeplotted2<-which(abs(y2)<100)

n<-50
y<-seq(from=-30,to=30,length.out=n)
z<-seq(from=-maxz,to=maxz,length.out=n)
yz<-expand.grid(y,z)
x1<-yz[,1]

xx<-c(x1,xz[tobeplotted2,1],xy2[tobeplotted3,1])
yy<-c(yz[,1],y2[tobeplotted2],xy2[tobeplotted3,2])
zz<-c(yz[,2],xz[tobeplotted2,2],z3[tobeplotted3])


gr<-c(rep(1,length(x1)),rep(2,length(y2[tobeplotted2])),rep(3,length(z3[tobeplotted3])))



plot3d(x1,yz[,1],yz[,2],col=1)
open3d()
plot3d(xz[tobeplotted2,1],y2[tobeplotted2],xz[tobeplotted2,2],col=2)
open3d()
plot3d(xy2[tobeplotted3,1],xy2[tobeplotted3,2],z3[tobeplotted3],col=3)
open3d()
plot3d(xx,yy,zz,col=gr)

xsurf<-cbind(x1,yz[,1],yz[,2])
ysurf<-cbind(xz[tobeplotted2,1],y2[tobeplotted2],xz[tobeplotted2,2])
zsurf<-cbind(xy2[tobeplotted3,1],xy2[tobeplotted3,2],z3[tobeplotted3])
Ax<-rbind(xssum,xsurf)
plot3d(Ax[,1],Ax[,2],Ax[,3],col=c(rep(4,length(xssum[,1])),rep(1,length(Ax[,1])-length(xssum[,1]))))
open3d()
Bx<-rbind(xssum,ysurf)
plot3d(Bx[,1],Bx[,2],Bx[,3],col=c(rep(4,length(xssum[,1])),rep(2,length(Bx[,1])-length(xssum[,1]))))
open3d()
Cx<-rbind(xssum,zsurf)
plot3d(Cx[,1],Cx[,2],Cx[,3],col=c(rep(4,length(xssum[,1])),rep(3,length(Bx[,1])-length(xssum[,1]))))
open3d()

Dx<-rbind(xssum,xsurf,ysurf,zsurf)
plot3d(Dx[,1],Dx[,2],Dx[,3],col=c(rep(4,length(xssum[,1])),rep(1,length(Ax[,1])-length(xssum[,1])),rep(2,length(Bx[,1])-length(xssum[,1])),rep(3,length(Bx[,1])-length(xssum[,1]))))