拡散とNeumann条件

• こちらのつづき
  • 拡散反応系までもっていきたいが...
• Neumann条件になっているか不安だったので
• 全体の量が保存されているかを見てみた
  • 下の枠の5行が総量のプロット
  • 初期値は33
  • 36くらいまで増加している
• 保存されていない...?
  • 立ち上がりまでに少しかかるところをみると
  • 境界からの流入が0になっていないのだろうか
• 追記
  • 保存されるようにしたのがこちら

N<-30
M<-40
T<-200

d<-1    #キョリの次元
dt<-1   #時間の次元
D<-0.12  #拡散係数

U<-array(0,c(N,M,T))

U[12:10,20:10,1]<-1    #初期条件

#拡散
for(t in 1:(T-1)){

for(i in 2:(N-1)){
for(j in 2:(M-1)){
du<-D*(sum(U[(i-1):(i+1),(j-1):(j+1),t])-9*U[i,j,t])/(d^2)*dt        #uの変化量duを計算
U[i,j,t+1]<-U[i,j,t]+du
}
}
U[,1,t+1]<-U[,2,t+1]  #Neumann条件
U[1,,t+1]<-U[2,,t+1]
U[N,,t+1]<-U[N-1,,t]
U[,M,t+1]<-U[,M-1,t]


}
min(U)

for (t in 1:T){
persp(U[,,t],col=3,theta=20,phi=30,zlim=c(0,1))
#image(U[,,t],col=topo.colors(10))
}

v<-c()
for(t in 1:T){
v<-c(v,sum(U[,,t]))
}
plot(v)