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

[戻る]