3D拡散のシミュレーション

• こちらを見ていて
  • 高次元で状態空間などを考えてみようと思った
  • 高次元をグラフにするのは難しいが、３次元+時間 なら今までも少し書いているのでできるはずと思って関係ないが書いてみた
  • ３次元空間×半径の大きさ×時間
  • Rのslice.index()で元の情報を取り出せる
• 拡散のシミュレーション
  • このときは2D
• 今回は3Dでのランダムウォーク
  • こちらのAと同じような計算方法
• Rの関数について
  • arrayの扱いの中で
  • slice.index(A,2) はarrayに対する row(),col()の一般化
  • arrayのdim()の中で知りたい空間の要素番号を返してくれる
  • rglのプロットについてはこちら

X <- 10
Y <- 10
Z <- 10
T <- 30

A<-array(0,c(X,Y,Z,T))
A[2:3,2:3,2:3,1]<-2

for(t in 1:(T-1)){
for(x in 1:X){
for(y in 1:Y){
for(z in 1:Z){
counter <- A[x,y,z,t]

while(counter >0){
counter <- counter -1
dx<-sample(-1:1,1);dy<-sample(-1:1,1);dz<-sample(-1:1,1)
newx<-sort(c(1,x+dx,X))[2]      #反射壁
newy<-sort(c(1,y+dy,Y))[2]
newz<-sort(c(1,z+dz,Z))[2]
A[newx,newy,newz,t+1]<-A[newx,newy,newz,t+1]+1
}

}
}
}
}

library(rgl)
for(t in 1:T){
open3d()
spheres3d(slice.index(A[,,,t],1),slice.index(A[,,,t],2),slice.index(A[,,,t],3),radius=A[,,,t]/4,col=6)
rgl.viewpoint(theta=30,phi=30)
rgl.bbox(color="#ffffff")
#filename<-paste(100+t,".png",sep="")
#rgl.snapshot(filename)
}