台形公式をリチャードソン補外でシンプソン法に。

ローレンツ型の関数を‐1から1の範囲で、いくつかの分割数に対して台形公式で数値積分します。出てきた積分値をリチャードソン補外をして改善します。その結果は理屈上シンプソン法によって求めたものと一致します。

実行結果

まずもっとも細かい分割数での座標値を配列生成子で生成し、被積分関数に入れて値を配列として求めます。ここで被積分関数を要素型で宣言することで、引数と関数値のDO LOOPの交換をします。
『y = f( [(x0 + i * h, i = 0, n)] )』


円周率πが求まるはずです。
『integer, parameter :: kd = kind(0.0e0)』の精度を変えることで、単精度や倍精度で計算できます。

ソース・プログラム

program Integrate
  implicit none
  integer, parameter :: kd = kind(0.0e0) 
  integer, parameter :: nn = 5, n = 2**nn
  real(kd) :: y(0:n), y1(nn), y2(nn - 1), hh(nn), h, x0, x1
  integer :: i, k(nn)
  x0 = -1.0_kd
  x1 =  1.0_kd
  h = (x1 - x0) / real(n, kd)
  y = f( [(x0 + i * h, i = 0, n)] )
  k = [(n / 2**i, i = 1, nn)]
  
  forall(i = 1:nn) & ! trapezoid rule for 2**i intervals
  y1(i) = h * k(i) * ( y(0) / 2 + sum(y(k(i):n - 1:k(i))) + y(n) / 2 )
  y2 = ( 4 * y1(2:nn) - y1(1:nn - 1) ) / 3 ! Richardson extrapolation
  
  print *, 2, real(k(1) * h), y1(1)
  do i = 2, nn
   print *, 2**i, real(k(i) * h), y1(i), y2(i - 1) ! trapezoidal & Simpson's rule
  end do
  stop
 
contains
   elemental real(kd) function f(x)
     real(kd), intent(in) :: x
     f = 2.0_kd / (1.0_kd + x**2)
     return
   end function f
end program Integrate