移動項 境界 ベクトル場
- 拡散方程式のシミュレーション
- 昨日の記事と同じベクトル場を用意する
- 境界からの流出入をなくす
- ソースのなかで (i != 1)や(i != Nx)の論理値 をかけることで境界では流出入が0になる
- 総量の保存について
- 初期条件によることであるが t=0 で36、 t=300 で 34.2 になった
- 端で流出している
- 今回は真ん中まで集まってくる
- やはりまだスケーリングを気にしないといけない
- 移動量の計算を3回繰り返しているの部分があるので、計算を工夫したほうがいいかもしれない
X<-30 #空間 Y<-30 T<-100 #時間 dx<-1 dy<-1 dt<-0.5 #時間間隔 Nx<-X/dx #格子の数 Ny<-Y/dy Nt<-T/dt #計算する回数 phi<-pi/8 Ox<-X/2 Oy<-Y/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/3 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<-array(c(runif(Nx*Ny)/10,rep(0,Nx*Ny*(Nt-1))),c(Nx,Ny,Nt)) U[5:10,5:10,1]<-1 #初期分布 c<-0.5 #c1<-0.1 #c2<-0.1 #nux<-(1-c1^2)*dx^2/dt #nuy<-(1-c2^2)*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<-nux/2*(sum(U[max(1,i-1):min(Nx,i+1),j,t])-kx*U[i,j,t])/(dx^2)*dt+c*(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]-abs(V[i,j,1])*U[i,j,t])*dt duy<-nuy/2*(sum(U[i,max(1,j-1):min(Ny,j+1),t])-ky*U[i,j,t])/(dy^2)*dt+c*(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]-abs(V[i,j,2])*U[i,j,t])*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,0.3))) #image(U[,,t],col=topo.colors(50)) print(sum(U[,,t])) }