ランダムウォークと拡散
- ランダムウォークによる拡散現象について
- スケーリングについて
- こちらの記事を正確に
- 下の教科書p.23
- 1次元ランダムウォーク
- 空間と時間は離散的に考えることにして
- の整数倍 の格子点上
- の整数倍 の時刻
- 移動の確率は
- 確率でからへ
- 確率でからへ
- 空間と時間は離散的に考えることにして
- 存在確率について
- Taylor展開と以下の仮定を用いて
- 仮定1.
- 仮定2.
-
- 仮定1 単位時間あたりの移動量の期待値が
- 仮定2 単位時間あたりの移動量の分散が
- 仮定1と2 から がでる
- 以下3種類の計算方法
#パラメータと初期分布 p<-1/2-0.1 #右へ移動する確率 q<-1-p #左へ移動する確率 X<-100 #空間 T<-10 #時間 M<-100 #点の数の最大値 dx<-1 #格子幅 dt<-0.1 #時間間隔 Nx<-X/dx #格子の数 Nt<-T/dt #計算する回数 c<-(p-q)*dx/dt #移動量の平均 仮定1 nu<-(1-(p-q)^2)*dx^2/dt #移動量の分散 仮定2 A<-matrix(0,Nx,Nt) s<-sample(2:Nx,1) A[s:(s-1),1]<-M #初期分布 C<-B<-A
#Aの計算 #実際にランダムウォークする for(j in 1:(Nt-1)){ for(i in 1:Nx){ counter<-A[i,j] while(counter>0){ counter<-counter-1 r<-runif(1) if(r < p){A[min(i+1,Nx),j+1]<-A[min(i+1,Nx),j+1]+1 }else if(r >= p){A[max(i-1,1),j+1]<-A[max(i-1,1),j+1]+1 } } } }
#Bの計算 #確率を計算していく for(j in 2:Nt){ for(i in 2:(Nx-1)){ B[i,j]<-B[i-1,j-1]*p+B[i+1,j-1]*q B[1,j]<-B[1,j-1]*q+B[2,j-1]*q #反射壁 B[Nx,j]<-B[Nx,j-1]*p+B[Nx-1,j-1]*p #反射壁 } }
#Cの計算 #拡散方程式を計算 for(j in 1:(Nt-1)){ for(i in 2:(Nx-1)){ dc<-(nu/2*(C[i-1,j]+C[i+1,j]-2*C[i,j])/(dx^2)-c*(C[i+1,j]-C[i-1,j])/(2*dx))*dt C[i,j+1]<-C[i,j]+dc } C[1,j+1]<-C[2,j+1] C[Nx,j+1]<-C[Nx-1,j+1] }
#結果のプロット library(rgl) plot3d(row(A),col(A),A,col=2,xlab="position",ylab="time",zlab="count",main="A") open3d() plot3d(row(B),col(B),B,col=3,xlab="position",ylab="time",zlab="count",main="B") open3d() plot3d(row(C),col(C),C,col=4,xlab="position",ylab="time",zlab="count",main="C")