一様な電磁場での荷電粒子の挙動のプログラミング

質問

一様な電磁場での電子・陽子の動きをfortran77でプログラミングしようと思ったのですが、なかなかうまくいきません。数値設定も全て自分でやっているのであまり自信はないのですが・・・
何かおかしなところや直した方がいいところがあれば助言の方よろしくお願いします。。。

implicit none

real*8 me,re(3),ve(3),E(3),B(3),qe,Te,dte
real*8 mp,rp(3),vp(3),Tp,dtp,time,qp

c 初期値の設定
me=9.1d-31 ! 電子の質量
mp=1.67d-27 ! 陽子の質量
qe=-1.6d-19 ! 電子の電荷
qp=1.6d-19 ! 陽子の電荷
re(1)=0.0
re(2)=0.0 ! 電子の位置
re(3)=0.0
ve(1)=1.0d4
ve(2)=0.0 ! 電子の速度
ve(3)=0.0
rp(1)=0.0
rp(2)=0.0 ! 陽子の位置
rp(3)=0.0
vp(1)=1.0d4
vp(2)=0.0 ! 陽子の速度
vp(3)=0.0
E(1)=0.0
E(2)=1.0d-3 ! 電場
E(3)=0.0
B(1)=0.0
B(2)=0.0 ! 磁場
B(3)=1.0d-5

dte=-1.0d-5
dtp=1.0d-4

open(unit=10,file='densi.dat',status='unknown')

do time=0.0,0.1,dte
write(10,*) time,re,ve
call lorentz(re,ve,E,B,qe,me,dte,Te)
enddo

close(10)

open(unit=20,file='yousi.dat',status='unknown')

do time=0.0,0.1,dtp
write(20,*) time,rp,vp
call lorentz(rp,vp,E,B,qp,mp,dtp,Tp)
enddo

close(20)

stop
end

c ローレンツの計算
subroutine lorentz(r,v,E,B,q,m,dt,T)
real*8 r(3),v(3),E(3),B(3),q,m,dt,T

v(1)=v(1)+q*(E(1)+(v(2)*B(3)-v(3)*B(2)))*dt/m
v(2)=v(2)+q*(E(2)+(v(3)*B(1)-v(1)*B(3)))*dt/m
v(3)=v(3)+q*(E(3)+(v(1)*B(2)-v(2)*B(1)))*dt/m

r(1)=r(1)+v(1)*dt
r(2)=r(2)+v(2)*dt
r(3)=r(3)+v(3)*dt

T=(m*((v(1)**2)+(v(2)**2)+(v(3)**2)))/2

return
end

回答No.1

foobar

微分方程式を数値的に解くとき、時間の刻み幅が大きすぎると計算がうまくいかないことがままあります。
刻み幅が影響しているかどうかは
・dtの大きさを小さくして（例えば1/10とか）計算してみて、結果が大きく変わらなければOK
というような方法でチェックできるかと思います。