sample 1
C 振り子の方程式をEuler法で解く
C-----------------------------------------------------------------------------
implicit none
integer j,jmax
parameter(jmax=500)
real x0,v0
real x(0:jmax),v(0:jmax),a(0:jmax)
real t(0:jmax),dt,tmax
real xh,vh,ah
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)) !運動方程式:a=-sinx
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
[戻る]