!    2011/06/23
!    euler法
!
!    dx/dt=-wx
!    x(0)=1,w=pi
!    (0<=t<=1.0)


!---------------メインプログラム-------------------
program euler2
  implicit none
  integer :: loop ! ループカウンター
  integer :: in   ! ループ変数
  integer :: hh   ! 刻む回数
  real :: min     ! 時間の始端
  real :: max     ! 時間の終端
  real :: time    ! 変数t
  real :: xx      ! x_nの値
  real :: xn      ! x_(n+1)の値
  real :: error   ! 数値解と解析解との誤差
  real :: pi      ! 円周率

  ! 円周率の計算
  pi=4.0e0*atan(1.0e0)

  ! x(0)=0であると設定
  xx=1.0e0

  !時間の始端と終端を設定
  min=0.0e0
  max=1.0e0

  !刻み幅hhを決める
  write(*,*) '刻む回数を入力 >'
  read(*,*) hh

  ! ループ変数の初期値を設定
  loop=0

  ! timeが1.0になるまでループ, 最初の0.0も含んで表示
  do in=0,hh

     ! 時間tの設定
     time=(max-min)*real(in)/real(hh)

     ! ループ変数を更新
     loop=loop+1

     ! xnを求める
     xn=xx+(max-min)/real(hh)*(-pi*xx)

     ! 誤差の更新
     error=abs(xx-exp(-pi*time))

     ! ループ回数と時間と解と誤差を表示
     write(*,*) loop,time,xx,error

     ! xxの更新
     xx=xn

  end do

end program euler2
