速度Verlet法が、2段のオイラー法と同じ結果になってしまう。

 投稿者:SECOND  投稿日:2017年 3月29日(水)12時16分20秒
 
速度Verlet法が、2段のオイラー法と区別できない観測値となり、
どなたか、以下のレポートについて、修正の加筆を御願い。 山中さん助けて!

--------------------------------------------------------
慶大の計算物理学特論

※シンボル名は、位置q() 速度v() 力F() 質量m  時間t

6.1.4 速度Verlet法(もう1つの同等な方法)

              F(q(t)) (⊿t)^2
 q(t+⊿t)=q(t)+v(t)⊿t+──── ───
               m       2

         ⊿t ┌ F(q(t))   F(q(t+⊿t)) ┐
 v(t+⊿t)=v(t)+──│────+────── │
          2 └    m          m       ┘

--------------------------------------------------------
以下の様に、カッコを、括りなおす。

         ┌       F(q(t)) ⊿t ┐
 q(t+⊿t)=q(t)+│v(t)+──── ──│⊿t
         └          m      2 ┘

      ┌       F(q(t)) ⊿t ┐  F(q(t+⊿t))  ⊿t
 v(t+⊿t)=│v(t)+──── ──│+────── ──
      └          m      2 ┘       m         2

--------------------------------------------------------
上の式を、十進BASICプログラム風に直すと・・
DO
    LET v=(v+ F/m *⊿t/2)
    LET q=q+v *⊿t
    LET F= force(q)         ! 更新された q の、F
    LET v=v+ F/m *⊿t/2
LOOP

ループなので、4行目の計算は、Fの更新なしに、1行目へ連続しており、
以下の様にしても、最初の1回だけの違い。( 初期値 F=0 なら、影響なし)
DO
    LET v=v+F/m *⊿t/2+F/m *⊿t/2
    LET q=q+v *⊿t
    LET F= force(q)
LOOP

式としては、以下と同一。但し、BASIC 側の計算の誤差が異なる分、累積で違っていく。
DO
    LET v=v+F/m *⊿t
    LET q=q+v *⊿t
    LET F= force(q)
LOOP

-----------------------------------------------------------------
実際に、これが、同一か否か、既存の投稿プログラムで、測定した。
結果は、やはり一致する。
-----------------------------------------------------------------

#4244 のプログラムの
以下の、現状のパーツ2つを、すげ替えて、比較します。

FUNCTION sysTemp
   LET v2 = 0
   FOR i=1 TO Nmt
      LET v2 = v2 + ABS(vxy(i))^2-vz(i)^2        ! vx^2 + vy^2 + vz^2
      !LET v2 = v2 + ABS(ABS(vxy(i))+vz(i))^2     ! vx^2 + vy^2 + vz^2
   NEXT i
   LET sysTemp = mass*v2/Nmt/3/k_                ! {∑(.5*m*v^2) /N} *2/3 /k
END FUNCTION

SUB moveParticles
   FOR i=1 TO Nmt
      LET vxy(i) = vxy(i)+a_*fxy(i)              !a_ = 0.5*dt/mass
      LET  vz(i) =  vz(i)+ a_*fz(i)
      LET pxy(i) = pxy(i)+vxy(i)*dt
      LET  pz(i) =  pz(i)+vz(i)*dt
   NEXT i
   CALL calcForce
   FOR i=1 TO Nmt
      LET vxy(i) = vxy(i)+a_*fxy(i)              !a_ = 0.5*dt/mass
      LET  vz(i) =  vz(i)+ a_*fz(i)
   NEXT i
   LET sysTime=sysTime+dt
   LET Az=Az+dt*1e11
END SUB

現状の、上式での実行結果
time =    .00 (ps) temp = 150.00 (K)
time =   2.00 (ps) temp = 155.62 (K)
time =   4.00 (ps) temp = 194.59 (K)
time =   6.00 (ps) temp = 180.24 (K)
time =   8.00 (ps) temp = 181.89 (K)
time =  10.00 (ps) temp = 188.48 (K)
time =  12.00 (ps) temp = 192.71 (K)
time =  14.00 (ps) temp = 196.33 (K)
time =  16.00 (ps) temp = 193.79 (K)
time =  18.00 (ps) temp = 187.88 (K)
time =  20.00 (ps) temp = 171.75 (K)
         (
          )
time = 992.00 (ps) temp = 180.61 (K)
time = 994.00 (ps) temp = 194.11 (K)
time = 996.00 (ps) temp = 188.18 (K)
time = 998.00 (ps) temp = 196.90 (K)
time =1000.00 (ps) temp = 205.22 (K)

-----------------------------------------------------
上のパーツ2つを、以下に、取替える。

FUNCTION sysTemp
   LET v2 = 0                                     !v2:速度 絶対値の2乗
   FOR i=1 TO Nmt
   !---------------- 温度表示を、同位置にする為、vxy,vz の1段先を、観測。
      LET w1=vxy(i)+a_*fxy(i)                     !a_ = 0.5*dt/mass
      LET w2=vz(i) +a_*fz(i)
      LET v2 = v2 + ABS( w1)^2 -( w2)^2           !re(vx+⊿)^2 +im(vy+⊿)^2 +im(vz+⊿)^2
      !---------------- 元の、温度表示。(停止)
      ! LET v2 = v2 + ABS(vxy(i))^2 -vz(i)^2      !re(vxy)^2 +im(vxy)^2 +im(vz)^2
   NEXT i
   LET sysTemp = mass*v2/Nmt/3/k_                 ! {∑(.5*mass*v^2) /Nmt} *2/3 /k
END FUNCTION

SUB moveParticles
   FOR i=1 TO Nmt
      LET vxy(i) = vxy(i)+a_*fxy(i) +a_*fxy(i)    !a_ = 0.5*dt/mass
      LET  vz(i) =  vz(i)+ a_*fz(i) + a_*fz(i)
      LET pxy(i) = pxy(i)+vxy(i)*dt
      LET  pz(i) =  pz(i)+vz(i)*dt
   NEXT i
   CALL calcForce
   LET sysTime=sysTime+dt
   LET Az=Az+dt*1e11
END SUB

!  LET vxy(i) = vxy(i)+fxy(i)/mass*dt   !左記は、式として、上の2行と同一ですが、
!  LET  vz(i) =  vz(i)+ fz(i)/mass*dt   ! 計算誤差の違いで、time=34.00(ps) 以降ズレる

上式での実行結果
time =    .00 (ps) temp = 150.00 (K)
time =   2.00 (ps) temp = 155.62 (K)
time =   4.00 (ps) temp = 194.59 (K)
time =   6.00 (ps) temp = 180.24 (K)
time =   8.00 (ps) temp = 181.89 (K)
time =  10.00 (ps) temp = 188.48 (K)
time =  12.00 (ps) temp = 192.71 (K)
time =  14.00 (ps) temp = 196.33 (K)
time =  16.00 (ps) temp = 193.79 (K)
time =  18.00 (ps) temp = 187.88 (K)
time =  20.00 (ps) temp = 171.75 (K)
         (
          )
time = 992.00 (ps) temp = 180.61 (K)
time = 994.00 (ps) temp = 194.11 (K)
time = 996.00 (ps) temp = 188.18 (K)
time = 998.00 (ps) temp = 196.90 (K)
time =1000.00 (ps) temp = 205.22 (K)


---------------------------------------------------------------
このプログラムの雛型になっている2次元においても、同様。
#4233 の、
現状を、下のパーツ2つに、取り替えて、比較。

EXTERNAL SUB setInitialCondition
   DECLARE EXTERNAL SUB ajustVelocity
   ! set particles
   RANDOMIZE 5       ←再現させるため引数を、追加する。

    (
     )

EXTERNAL SUB moveParticles ! velocity Verlet method
   DECLARE EXTERNAL SUB calcForce
   LET a = 0.5*dt/mass
   FOR i=1 TO nMolec
      LET vx(i) = vx(i)+a*ffx(i)+a*ffx(i)
      LET vy(i) = vy(i)+a*ffy(i)+a*ffy(i)
      LET xx(i) = xx(i)+vx(i)*dt
      LET yy(i) = yy(i)+vy(i)*dt
   NEXT i
   CALL calcForce
   LET sysTime=sysTime+dt
END SUB

    (
     )

EXTERNAL FUNCTION systemTemprature
   LET kB = 1.38e-23 ! Boltzman's constant (J/K)
   LET ek= 0.0       !kinetic energy (J)
   FOR i=1 TO nMolec
   !---------------- 温度表示を、同位置にする為、vx,vy の1段先を、観測。
      LET a = 0.5*dt/mass
      LET w1=vx(i) +a*ffx(i)
      LET w2=vy(i) +a*ffy(i)
      LET ek = ek + 0.5*mass*( w1^2 + w2^2)
      !---------------- 元の、温度表示。
      ! LET ek = ek + 0.5*mass*(vx(i)^2+vy(i)^2)
   NEXT i
   LET systemTemprature = ek/(nMolec*kB)
END FUNCTION
 

Re: 速度Verlet法が、2段のオイラー法と同じ結果になってしまう。

 投稿者:mike  投稿日:2017年 3月30日(木)10時30分46秒
  SECONDさんへのお返事です。SECONDさん、興味深い事実を提供していただき、ありがとうございます。

フォローしてみたのですが、述べられていることは正しいと思います。

velocity-Verlet法では本来
  t=0、x(t=0)、v(t=0) を与える(初期化)

  f = f(x) !****
  do
  v(t+dt/2) = v(t) + (f(t)/m)*dt/2
  x(t+dt) = x(t) + v(t+dt/2)*dt
  f(f+dt) = f(x(t+dt))
  v(t+dt) = v(t+dt/2) + (f(t+dt)/m)*dt/2
  loop
とする必要がありますが、雛形プログラムでは初期の配置がランダムなため、
最初のループでf=0でも影響が少ないだろうと考えて!****の行を省略しました。
このため、オイラー法のような形のプログラムと同じ結果を与えたと考えます。

オイラー法のような形のにする利点は計算がO(N)少なくてすむことですが、
時間がx,fとvでずれているため、分子のと速度から色々な統計量(たとえば全エネルギーΣi(ΣUj(xi)+0.5mvi^2)
を計算するため、時間のズレを補正する(O(N)の計算量が必要)ため、実際にはメリットが少ないと思えます。

慶大の計算物理学特論の
能勢先生は分子動力学ではNose-Hoover法など、とても有名な先生でした。

 

戻る