マスを仕切ることについて

拡散 シミュレーション

KABIRA2011-02-25

  • こちらなど、隣り合うマスにどれだけ移動するかを考えている
  • 平均や分散で移動量を計算しようと思っているが
  • そもそも平均を使って移動量を計算することで誤差がどれくらい生じるだろうか
    • これが不可避なのかどうかはまた別問題で、この方法自体のもつ誤差がどれくらいかということ
    • むしろ平均をとることではなく、マスに区切っていることによる誤差がみられる
  • 以下のソースでは等速直線運動することを考える
    • A[x,t]
    • 最初A[2,1]にのみ分布している
    • 速度vで等速直線運動させる
      • 時間 s の間に 位置 x から x + v*s へ移動するが
      • A[x,t] から A[j,t+1]とA[j+1,t+1]の２マスに分けられてしまう
      • ２マスを平均して移動量が x + v*s になるようにそれぞれのマスへの移動量は
      • j=[x+vs] として
        • A[x,t]*(j+1-(x+v*s))
        • A[x,t]*(x+v*s-j)
      • とする
  • Aでは徐々にピークが低くなっているということは裾がのびていくようだ
    • おそらく各時間の位置を計算していくときに、正確な位置情報が失われている
    • 精度が1マス分の精度に切り捨てられるために誤差が生まれるということだろう
  • 一方、Bは常に初期位置の情報を持っているので情報が失われない（分散しない）
X<-60
T<-30
s<-1   #時間ステップ
N<-1
v<-1.3   #速度
B<-A<-matrix(0,X,T)
B[2,1]<-A[2,1]<-N
for(t in 1:(T-1)){
for(x in 1:X){
j<-floor(x + v*s)
if(j  <= X){A[j,t+1]<-A[j,t+1]+A[x,t]*(j+1-(x+v*s))}
if(j + 1 <= X){A[j+1,t+1]<-A[j+1,t+1]+A[x,t]*(x+v*s-j)}
}
}

persp(A,col=rainbow(T*X),theta=30,phi=30,xlab="X",ylab="T")
#image(A,main="A")

for(t in 1:(T-1)){
j<-floor(2 + v*t)
if(j <= X){B[j,t+1]<-N*(j+1-(2+v*t))}
if(j+1 <= X){B[j+1,t+1]<-N*(2+v*t-j)}
}
#persp(B,col=rainbow(T*X),theta=30,phi=30,xlab="X",ylab="T")
#image(B,main="B")

#ここからは平均と分散を計算してみる
mA<-mB<-mA2<-mB2<-rep(0,length=T)
VA<-VB<-rep(0,length=T)
for(t in 1:T){
for(x in 1:X){
mA[t]<-mA[t]+A[x,t]*x    #これがAの場合の移動の平均
mB[t]<-mB[t]+B[x,t]*x
mA2[t]<-mA2[t]+A[x,t]*x^2
mB2[t]<-mB2[t]+B[x,t]*x^2
}
}

VA<-mA2-mA^2    #これがAの場合の分散の経時変化
VB<-mB2-mB^2

par(mfrow=c(2,2))
plot(mA,main="mean A",ylim=c(0,max(mA)),col=2)
plot(mB,main="mean B",ylim=c(0,max(mA)),col=2)
plot(VA,main="Var A",ylim=c(0,max(VA)),col=3)
plot(VB,main="Var B",ylim=c(0,max(VA)),col=3)

SA<-SB<-rep(0,length=T)
for(t in 1:T){
for(x in 1:X){
if(A[x,t]>0){SA[t]<-SA[t]-log(A[x,t])*A[x,t]}
if(B[x,t]>0){SB[t]<-SB[t]-log(B[x,t])*B[x,t]}
}
}

plot(SA,ylim=c(0,max(SA)),col=4)
par(new="TRUE")
plot(SB,ylim=c(0,max(SA)),col=6)