拡散の方向

• こちらでは拡散する方向はxとy方向に限られている
• ナナメ方向にも移動する項を考えたい
• 拡散のシンプルなシミュレーションから考え直して、斜め方向を考慮することにする
• まずはそのシンプルな拡散

Nx<-30
Ny<-30
Nt<-30
U<-tempU<-matrix(0,Nx,Ny)
U[5:10,5:10]<-1
k1<-k2<-0.2

for (t in 1:Nt){
for (x in 1:Nx){
for (y in 1:Ny){
	du<-0
	if(x != 1){du <- du + k1*(U[x-1,y]-U[x,y])}
	if(x != Nx){du <- du + k1*(U[x+1,y]-U[x,y])}
	if(y != 1){du <- du + k2*(U[x,y-1]-U[x,y])}
	if(y != Ny){du <- du + k2*(U[x,y+1]-U[x,y])}
	tempU[x,y] <- U[x,y] + du
}
}
U<-tempU
persp(U,col=3,theta=30,phi=30,zlim=c(0,1))
}

• 今の場合構造が一様であるのでXY平面を行列として扱って計算することもできてこちらのほうが速い

Nx<-30
Ny<-30
Nt<-30
U<-tempU<-matrix(0,Nx,Ny)
U[5:10,5:10]<-1
k1<-k2<-0.2

for (t in 1:Nt){
	dUx<-U[2:Nx,]-U[1:(Nx-1),]
	dUy<-U[,2:Ny]-U[,1:(Ny-1)]
	
	U[1:(Nx-1),]<-U[1:(Nx-1),]+k1*dUx
	U[2:Nx,]<-U[2:Nx,]-k1*dUx
	U[,1:(Ny-1)]<-U[,1:(Ny-1)]+k2*dUy
	U[,2:Ny]<-U[,2:Ny]-k2*dUy
	
	persp(U,col="lawngreen",theta=30,phi=30,zlim=c(0,1))
}