拡散の係数が関数の場合

• 拡散の係数kの値を自由に与える
• 昨日の記事の係数(k1,..,k4)を各点で与える行列を作る

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

k1[1:ceiling(Nx/2),]<-k2[1:ceiling(Nx/2),]<-k3[1:ceiling(Nx/2),]<-k4[1:ceiling(Nx/2),]<-0

for (t in 1:Nt){
	dUx<-U[2:Nx,]-U[1:(Nx-1),]
	dUy<-U[,2:Ny]-U[,1:(Ny-1)]
	dUxy<-U[2:Nx,2:Ny]-U[1:Nx-1,1:Ny-1]
	dUyx<-U[1:Nx-1,2:Ny]-U[2:Nx,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
	
	U[1:Nx-1,1:Ny-1]<-U[1:Nx-1,1:Ny-1]+k3*dUxy
	U[2:Nx,2:Ny]<-U[2:Nx,2:Ny]-k3*dUxy
	U[1:Nx-1,2:Ny]<-U[1:Nx-1,2:Ny]-k4*dUyx
	U[2:Nx,1:Ny-1]<-U[2:Nx,1:Ny-1]+k4*dUyx
	
	persp(U,col=6,theta=30,phi=30,zlim=c(0,1))
#	image(U,col=topo.colors(100))
#	print(sum(U))
}