移動項 拡散項 境界 繊維方向

• 拡散方程式のシミュレーション
• 昨日のものとほとんど同じにみえるが
  • 一部修正
  • 移動項の中身を修正
  • 境界条件の計算方法を違う方法でかくほうがいいのだろうか
  • パラメータ（cやnu）にdxやdyの値をかけて補正がきくようにしたい
• ムービ−はDirichlet条件で値０にして計算している
• 移動について
• 人の移動を考えると
  • 道を両方向に移動する
  • 移動の平均が移動項であらわされる
  • 移動の分散は拡散項であらわされる
• 道の走る方向に拡散しやすいと考えて
  • ベクトル場を道の走向とみて、拡散係数nuにかける
  • 斜め方向の移動がないので、ベクトル場による斜め方向の影響がうまくあらわせない

X<-20    #空間
Y<-20
T<-50    #時間


dx<-1
dy<-1
dt<-0.5  #時間間隔

Nx<-X/dx   #格子の数
Ny<-Y/dy
Nt<-T/dt   #計算する回数

phi<-pi/20
Ox<-Nx/2
Oy<-Ny/2
V<-array(0,c(Nx,Ny,2))

f<-function(x,y,Ox,Oy,phi){
#r<-sqrt((x-Ox)^2+(y-Oy)^2)
r<-1
th<-atan2(y-Oy,x-Ox)
v<-c(-r*sin(th+phi),r*cos(th+phi))
v
}

for(i in 1:Nx){
for(j in 1:Ny){
V[i,j,]<-f(i,j,Ox,Oy,phi)
}
}

U<-array(0,c(Nx,Ny,Nt))
#U[5:10,5:10,1]<-1     #初期分布
U[floor(Nx/6):floor(Nx/3),floor(Ny/6):floor(Ny/3),1]<-1


cx<-0.1
cy<-0.1
c1<-cx*dx/dt
c2<-cy*dy/dt

nux<-0.05
nuy<-0.05
nu1<-nux*dx^2/dt
nu2<-nuy*dy^2/dt
 
#nux<-0.4
#nuy<-0.4

for(t in 1:(Nt-1)){
#for(i in 1:Nx){
#for(j in 1:Ny){

for(i in 2:(Nx-1)){
for(j in 2:(Ny-1)){

kx<-length(max(1,i-1):min(Nx,i+1))
ky<-length(max(1,j-1):min(Ny,j+1))

dux<-nu1/2*(sum(U[max(1,i-1):min(Nx,i+1),j,t])-kx*U[i,j,t])/(dx^2)*dt+c1*(max(0,V[max(i-1,1),j,1]*(i != 1))*U[max(i-1,1),j,t]-min(0,V[min(Nx,i+1),j,1]*(i != Nx))*U[min(i+1,Nx),j,t]-(max(0,V[i,j,1]*(i != Nx))-min(0,V[i,j,1]*(i != 1)))*U[i,j,t])/dx*dt

duy<-nu2/2*(sum(U[i,max(1,j-1):min(Ny,j+1),t])-ky*U[i,j,t])/(dy^2)*dt+c2*(max(0,V[i,max(j-1,1),2]*(j != 1))*U[i,max(j-1,1),t]-min(0,V[i,min(j+1,Ny),2]*(j != Ny))*U[i,min(j+1,Ny),t]-(max(0,V[i,j,2]*(i != Ny))-min(0,V[i,j,2]*(i != 1)))*U[i,j,t])/dy*dt

du<-dux+duy

U[i,j,t+1]<-U[i,j,t]+du
}
}
}
min(U)

for(t in 1:Nt){
persp(U[,,t],col=6,theta=20,phi=60,zlim=(c(0,1)))
#image(U[,,t],col=topo.colors(50))
print(sum(U[,,t]))
}

X<-20    #空間
Y<-20
T<-500    #時間


dx<-1
dy<-1
dt<-1  #時間間隔

Nx<-X/dx   #格子の数
Ny<-Y/dy
Nt<-T/dt   #計算する回数

phi<--pi/4
Ox<-Nx/2
Oy<-Ny/2
V<-array(0,c(Nx,Ny,2))

f<-function(x,y,Ox,Oy,phi){
#r<-sqrt((x-Ox)^2+(y-Oy)^2)
r<-1
th<-atan2(y-Oy,x-Ox)*0
v<-c(-r*sin(th+phi),r*cos(th+phi))
v
}

for(i in 1:Nx){
for(j in 1:Ny){
V[i,j,]<-f(i,j,Ox,Oy,phi)
}
}

U<-array(0,c(Nx,Ny,Nt))
#U[5:10,5:10,1]<-1     #初期分布
U[floor(Nx/6):floor(Nx/3),floor(Ny/6):floor(Ny/3),1]<-1


cx<-0.1
cy<-0.1
c1<-cx*dx/dt
c2<-cy*dy/dt

nux<-0.1
nuy<-0.1
nu1<-nux*dx^2/dt
nu2<-nuy*dy^2/dt
 
#nux<-0.4
#nuy<-0.4

for(t in 1:(Nt-1)){
for(i in 1:Nx){
for(j in 1:Ny){

kx<-length(max(1,i-1):min(Nx,i+1))
ky<-length(max(1,j-1):min(Ny,j+1))

dux<-abs(V[i,j,1])*nu1/2*(sum(U[max(1,i-1):min(Nx,i+1),j,t])-kx*U[i,j,t])/(dx^2)*dt
duy<-abs(V[i,j,2])*nu2/2*(sum(U[i,max(1,j-1):min(Ny,j+1),t])-ky*U[i,j,t])/(dy^2)*dt
du<-dux+duy

U[i,j,t+1]<-U[i,j,t]+du
}
}
}
min(U)

for(t in 1:Nt){
persp(U[,,t],col=6,theta=20,phi=60,zlim=(c(0,1)))
image(U[,,t],col=topo.colors(50))
print(sum(U[,,t]))
}