シミュレーション物理4.pptVIP

シミュレーション物理6 運動方程式の方法: 惑星の軌道 出席のメール（件名に学生番号と氏名）に，中点法をサブルーチンを使って書いたプログラムを添付 今回の授業の目的 4次のRunge-Kutta法を用いて，惑星の軌道のシミュレーションを行う rk4.f90 4次のルンゲークッタ法 Dt4まで正しい。 高精度の数値計算に向いている。 これ以上次数を上げると，表式が複雑になりかえって計算時間がかかる。 自由落下のプログラムをrk4.f90を使うよう改良 program freefall !------------------------- ! This is a program to investigate the free fall ! 2005/5/12 Written by T. Ohtsuki ! Runge-Kutta method, subroutine used ! derivs is used as a subprogram, using contains. ! This makes parameters in the procedure derivs global variables !------------------------- implicit none ! Always begin with this statement integer,parameter::double=selected_real_kind(14) real(kind=double), parameter::zero=0.0_double,one=1.0_double,half=0.5_double real(kind=double),parameter::pi=3.1415926535897932 integer,parameter::nd=2 ! dimension of the equation real(kind=double),dimension(nd)::y,dydx,yout !位置，速度はベクトルで定義 real(kind=double)::deltat,t,tmax,analytic integer::i,nmax real(kind=double),parameter::g=9.8_double,gamma=1.0_double external:: derivsfreefall !---------main program--------- deltat=0.05_double t=zero y(1)=zero ! initial velocity，初期値も当然変わる y(2)=zero ! initial position tmax=5._double nmax=int((tmax-deltat/2._double)/deltat)+1 do i=1,nmax call derivsfreefall(nd,t,y,dydx) ! 微分係数を求める call rk4(nd,y,dydx,t,deltat,yout,derivsfreefall) !微分係数を使ってΔtだけ時間発展 y=yout end do tmax=deltat*nmax analytic=-g*tmax/gamma-g/gamma**2*exp(-gamma*tmax)+g/gamma**2 !ここら辺もいらない print *,tmax,y(2),analytic end program freefall 微分を求めるサブルーチン, derivsfreefall.f90として保存 !------------------- Subroutines --------------------- SUBROUTINE derivsfreefall(nd,x,y,dydx) !------------------------- ! subroutine derivs ! nd : dimension of y ! x : time (usually time independent) ! dydx: first derivative w.r.t. time, not a partial derivative !------------------------- IMPLICIT NONE integer,parameter::double=selected_real_kind(14) integer::nd REAL(kind=double), INTENT(IN) :: x REAL(kind=double), DIMENSION(nd), INTENT(IN) :: y