!    台形則
!    f(x)=e^(x) の 0<=x<=1 における積分
!    分割数が1, 2, 4, 8, 16, 32の場合
!    2013/08/01 村橋究理基 新規作成
!    2013/08/08 村橋究理基 更新 積分範囲を99 - 100に変更

program daikei2
implicit none

integer :: n     ! ループ変数
integer :: set   ! 台形の分割数
integer :: setmax ! 分割数の最大値
real(8) :: max   ! 積分区間の上端
real(8) :: min   ! 積分区間の下端
real(8) :: fmax  ! 被積分関数の引数maxの値
real(8) :: fmin  ! 被積分関数の引数minの値
real(8) :: div   ! 分割した一区間の幅
real(8) :: sum   ! 部分和の保存する変数
real(8) :: f     ! 関数の値を計算する変数
real(8) :: res   ! 数値計算の結果を入れる変数
real(8) :: error ! 数値解と解析解の差から得られる誤差

! 積分区間の上端と下端を設定
min = 99.0d0
max = 100.0d0

! 積分区間の上端と下端における関数fの値を計算
fmin = exp(min)
fmax = exp(max)

! 分割数を最初は1に設定
set = 1
setmax = 512

write(*,*) '分割数, 数値解, 誤差'

! 分割数がsetmaxになるまでループを回す
do while(set <= setmax)

   ! 一区間の幅を計算
   div = abs(max - min) / dble(set)

   ! 和を0で初期化
   sum = 0.0d0

   do n = 1, set - 1
      f = exp(min + dble(n)*div)

      sum = sum + f
      !write(*,*) n
   enddo

   ! 数値解の計算
   res = div*(fmin/2.0d0 + sum + fmax/2.0d0)

   ! 誤差の計算
   error = res - (exp(max) - exp(min))

   ! 結果の表示
   !write(*,*) set, error
   write(*, *) set, '&', error, '&'!, error*4.0d0, '\\ \hline'

   ! 分割数を二倍に増やす
   set = set * 2

enddo


endprogram

