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

      

[戻る]