ランダムウォークと拡散

• ランダムウォークによる拡散現象について
  • スケーリングについて
  • こちらの記事を正確に
• 下の教科書p.23
• 1次元ランダムウォーク
  • 空間と時間は離散的に考えることにして
    • の整数倍 の格子点上
    • の整数倍 の時刻
  • 移動の確率は
    • 確率でからへ
    • 確率でからへ
• 存在確率について
  • 
• Taylor展開と以下の仮定を用いて
  • 仮定1.
  • 仮定2.
• 
  • 仮定1 単位時間あたりの移動量の期待値が
  • 仮定2 単位時間あたりの移動量の分散が
  • 仮定1と2 から がでる
• 以下３種類の計算方法
  • A（赤）：乱数にしたがってランダムウォーク
  • B（緑）：確率で計算
  • C（青）：係数を求めてから拡散方程式を計算
  • いずれも境界はNeumannで反射壁
    • Cの境界条件の計算は正しくないのでこちらのようにすべき
    • Cのの一階微分の項はで近似してある
    • の一階微分の項がある（）ので一方に流れていく移流拡散方程式

#パラメータと初期分布
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")