!       台形則
!       f(x)=1/(1+x)の0<=x<=1での数値積分
!       伊理テキストP21の表4-1の作成
!       作成者    福士 邦彦
!       2011/05/12 新規作成
!       2011/06/08 最終更新


!----------------メインプログラム-----------------------------------

program daikei
  implicit none
  integer :: count     ! countはループ変数
  integer :: bunkatsu  ! bunkatsuは台形の分割数
  real(8) :: max       ! maxは積分の上端
  real(8) :: min       ! minは積分の下端
  real(8) :: fmax      ! fmaxはf(max)の値
  real(8) :: fmin      ! fminはf(min)の値
  real(8) :: floop     ! floopはループ中のそれぞれの値
  real(8) :: fsum      ! fsumはループ内の積算値
  real(8) :: menseki    ! mensekiは最終的な値,つまり面積,教科書のI_Nに相当する
  real(8) :: error     ! errorは真の値と計算値のずれ,教科書のI_N-Iに相当する
  
  !積分の下端を0に設定
  min=0.0d0

  !下端でのfの値,つまりf(0)の値を計算
  fmin=1.0d0/(1.0d0+min)

  !積分の上端を1に設定
  max=1.0d0

  !上端でのfの値,つまりf(1)の値を計算
  fmax=1.0d0/(1.0d0+max)

  !分割数の初期値を1に設定する
  bunkatsu=1

  !分割数が1,2,4,8,16,32,それぞれの時の台形の面積を計算するループ
  do while (bunkatsu<32.1)

     !fsumの初期値を0にする
     fsum=0.0d0

     !台形則のシグマ計算をするループ
     do count=1,bunkatsu-1

        ! f(a-n/N(b-a))の値の計算
        floop=1.0d0/(1.0d0+(min+dble(count)*(max-min)/dble(bunkatsu)))

        ! f(a-n/N(b-a))に対するシグマの部分の計算
        fsum=fsum+floop

     enddo

     !面積の値を計算
     menseki=(max-min)*(fmin/2+fsum+fmax/2)/dble(bunkatsu)

     !解析解と数値解との誤差を計算
     error=menseki-0.6931471805599453

     !分割数,数値解,誤差の順で表示
     write(*,*) bunkatsu,menseki,error

     !分割数を2倍にする
     bunkatsu=bunkatsu*2

  end do


end program daikei
