!    2011/06/23
!    euler法(倍精度)
!
!    2011/06/23 福士 邦彦 新規作成
!    2012/08/19 増田 和孝 最終更新
!    dx/dt=-wx
!    x(0)=1,w=pi
!    (0<=t<=1.0)
!    解析解はx = exp(-wt)


!---------------メインプログラム-------------------
program euler2
  implicit none
  integer :: in   ! ループ変数
  integer :: t_step   ! 刻む回数
  real(8) :: t_min     ! 時間の始端
  real(8) :: t_max     ! 時間の終端
  real(8) :: time    ! 変数t
  real(8) :: x_before      ! x_nの値
  real(8) :: x_current      ! x_(n+1)の値
  real(8) :: error   ! 数値解と解析解との誤差
  real(8) :: dt      ! 刻み幅
  real(8) :: pi      ! 円周率
  real(8) :: x_init  ! xの初期値

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

  ! 初期値の設定
  x_init = 1.0d0

  !時間の始端と終端を設定
  t_min=0.0d0
  t_max=1.0d0

  ! 初期化
  x_before=x_init
  x_current = x_before
  time = t_min
  error = 0.0d0
 
  !ステップ数t_stepを決める
  write(*,*) 'ステップ数を入力 >'
  read(*,*) t_step

  !ステップ数0以下の時のエラー処理
  if (t_step <= 0) then
   write(*,*) 'ステップ数は正の数でなければいけません'
   stop
  endif

  open(10, file = 'eular.data')

  ! 初期値を出力
  write(*,*) 'time ', 'x_current ', 'error'
  write(*,*) time, x_current, error
  write(10,*) time, x_current

  ! 刻み幅を計算する
  dt = abs(t_max - t_min)/real(t_step)

  ! timeが1.0になるまで表示
  do in=1,t_step

     ! 時間tの設定
     time = dt*real(in)

     ! x_currentを求める
     x_current=x_before + dt*(-pi*x_before)

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

     ! ループ回数と時間と解と誤差を表示
     write(*,*) time,x_current ,error
     write(10,*) time, x_current
   
  ! x_beforeの更新
     x_before=x_current

  end do
close(10)
end program euler2
