sample 2
c 運動方程式(a=-x)をEuler法で解く
c -----------------------------------------------------------------------------
implicit none
integer j,jmax
parameter(jmax=600)
real x0,v0,xh,vh,ah
real x(0:jmax),v(0:jmax),a(0:jmax)
real t(0:jmax),dt,tmax
c 初期条件
call getinit(x0,v0,dt,tmax)
x(0) = x0; v(0) = v0;
t(0) = 0.0
c 運動方程式を解く
do j=0, jmax-1
a(j)= -sin(x(j))
xh = x(j) + v(j)*(dt/2)
vh = v(j) + a(j)*(dt/2)
ah =-sin(xh)
x(j+1) = x(j) + vh*dt
v(j+1) = v(j) + ah*dt
t(j+1) = t(j) + dt
if( t(j) .gt. tmax ) exit
enddo
c 結果を出力する
do j=0,jmax-1
write(*,*) t(j),x(j)
if( t(j) .gt. tmax ) exit
enddo
stop
end
subroutine getinit(xinit,vinit,dt,tmax)
write(0,*) 'initial condition: x,v ='
read(*,*) xinit,vinit
write(0,*) 'dt, tmax='
read(*,*) dt,tmax
return
end
[戻る]