|
|
速度Ｖｅｒｌｅｔ法が、２段のオイラー法と区別できない観測値となり、
どなたか、以下のレポートについて、修正の加筆を御願い。 山中さん助けて！
--------------------------------------------------------
慶大の計算物理学特論
※シンボル名は、位置q() 速度v() 力F() 質量m 時間t
6.1.4 速度Ｖｅｒｌｅｔ法（もう１つの同等な方法）
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
ループなので、４行目の計算は、Ｆの更新なしに、１行目へ連続しており、
以下の様にしても、最初の１回だけの違い。( 初期値 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 のプログラムの
以下の、現状のパーツ２つを、すげ替えて、比較します。
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)
-----------------------------------------------------
上のパーツ２つを、以下に、取替える。
FUNCTION sysTemp
LET v2 = 0 !v2：速度 絶対値の２乗
FOR i=1 TO Nmt
!---------------- 温度表示を、同位置にする為、vxy,vz の１段先を、観測。
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 !左記は、式として、上の２行と同一ですが、
! 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)
---------------------------------------------------------------
このプログラムの雛型になっている２次元においても、同様。
#4233 の、
現状を、下のパーツ２つに、取り替えて、比較。
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 の１段先を、観測。
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
|
|