【メモ帳】Romberg 積分

計算結果

   1.64493103803760        1.64493061750701        1.64492430954817
   1.64492523592403        1.64490909117507        1.64486343605578
 pi^2/6
   1.64493406684823
続行するには何かキーを押してください . . .

プログラム

module m_mod
    implicit none
    integer, parameter :: kd = kind(0.0d0)
    real(kd), parameter :: pi = 4 * atan(1.0_kd)
contains
    pure real(kd) function f(x)
        real(kd), intent(in) :: x
        f = acos(1.0_kd / (1.0_kd + 2.0_kd * cos(x)))
    end function f
    
    pure function romberg(func, x0, x1, n) result(s) 
        procedure(f) :: func
        real(kd) :: s(6)
        real(kd), intent(in) :: x0, x1
        integer , intent(in) :: n
        real(kd) :: h, y(0:n), s0, s1, s2, t1, t2, r2
        integer :: i
        h = (x1 - x0) / n
        forall(i = 0:n) y(i) = func(x0 + i * h)
        s2 =     h * ( y(0) / 2 + sum(y(1:n-1  )) + y(n) / 2 )
        s1 = 2 * h * ( y(0) / 2 + sum(y(2:n-1:2)) + y(n) / 2 )
        s0 = 4 * h * ( y(0) / 2 + sum(y(4:n-1:4)) + y(n) / 2 )
        t2 = (4 * s2 - s1) / 3
        t1 = (4 * s1 - s0) / 3
        r2 = (16 * t2 - t1) / 15
        s = [r2,t2,t1,s2,s1,s0]
    end function romberg
end module m_mod

program test
    use m_mod
    implicit none
    integer, parameter :: n = 2**11
    real(kd) :: x0, x1
    integer :: i
    x0 = 0.0_kd 
    x1 = pi / 2.0_kd 
    print *, romberg(f, x0, x1, n)
    print *, 'pi^2/6'
    print *, pi**2 / 6
end program test