构建Romberg积分表

Question

我正在尝试编写一个 Fortran 程序来生成一个 Romberg 积分表。在 R.L.Burden 和 J.D.Faires 第 9 版的Numerical Analysis 一书中第 4.5 章中有一个算法。到目前为止，我已经写了这个

implicit none
integer,parameter::n=4
real::a,b,f,r(n,n),h,sum1
integer::i,k,j,m,l
open(1,file='out.txt')
a=0.
b=1.
h=b-a
r(1,1)=.5*h*(f(a)+f(b))
write(1,*)r(1,1)
do i=2,n
    sum1=0.
    do k=1,2**(i-2)
        sum1=sum1+f(a+(k-.5)*h)
    enddo
    r(2,1)=.5*(r(1,1)+h*sum1)
    do j=2,i
        r(2,j)=r(2,j-1)+(r(2,j-1)-r(1,j-1))/(4**(j-1)-1)
        write(1,*)((r(m,l),m=2,2),l=1,i)
    enddo
    h=h/2.
    do j=1,i
        r(1,j)=r(2,j)
    enddo
enddo

end

real function f(x)
implicit none
real,intent(in)::x

f=1/(1+x**2)

end function

这个程序给出以下输出：

  0.750000000    
  0.774999976      0.783333302    
  0.782794118      0.785392165       3.56011134E-22
  0.782794118      0.785392165      0.785529435    
  0.784747124      0.785398126      0.785529435       7.30006976E+28
  0.784747124      0.785398126      0.785398543       7.30006976E+28
  0.784747124      0.785398126      0.785398543      0.785396457    

但它应该给出这个：

0.7500000000 
0.7750000000 0.7833333333 
0.7827941176 0.7853921567 0.7855294120 
0.7847471236 0.7853981253 0.7853985227 0.7853964451 
0.7852354030 0.7853981627 0.7853981647 0.7853981590 0.7853981659 

以上是由 Maple 编写的程序完成的。 Maple中的程序是

>   romberg := proc(f::algebraic, a, b, N,print_table) 
 local R,h,k,row,col; 
 R := array(0..N,0..N); 

 # Compute column 0, Trapezoid Rule approximations of 
 #                   1,2,4,8,..2^N subintervals 
 h := evalf(b - a); 
 R[0,0] := evalf(h/2 * (f(a)+f(b))); 
 for row from 1 to N do; 
   h := h/2; 
   R[row,0] := evalf(0.5*R[row-1,0] + 
                     sum(h*f(a+(2*k-1)*h),k=1..2^(row-1))); 
   # Compute [row,1]:[row,row], via Richardson extrapolation 
   for col from 1 to row do; 
     R[row,col] := ((4^col)*R[row,col-1] - R[row-1,col-1]) / 
                       (4^col - 1); 
   end do; 
 end do; 

 # Display results if requested 
 if (print_table) then 
       for row from 0 to N do; 
     for col from 0 to row do; 
       printf("%12.10f ",R[row,col]); 
     end do; 
     printf("\n"); 
   end do; 
 end if; 

 return(R[N,N]); 

end proc:
f:=x->1/(1+x^2);
val:=romberg(f,0,1,4,true)

那么现在如何使用 Fortran 程序来获得与 Maple 程序相同的结果？

Answer 1

您的程序太长，仅通过阅读代码就无法调试（并且超过 5 年我没有写过一行 f95），但是从您的结果表中可以明显看出，在 F90 中表示的是实变量作为单精度（32 位）IEEE 754 浮点数，而 maple 使用更高的精度。

事实上，对于 32 位精度，0.775 变成（大约）0.77499998，而对于 64 位，它是 0.77500000000000002。 完全错误的数字（例如 3.56011134E-22）可能是由于下溢/溢出或取消造成的。

您可以使用合适的编译器选项强制 64 位精度，或者通过指定适当的真实类型，请参阅 http://fortranwiki.org/fortran/show/Real+precision

integer, parameter :: dp = kind(1.d0)
real(kind=dp) :: a

Answer 2

maple 程序和 fortran 源代码之间存在许多差异。

• maple 程序的结果数组的维度是从 0 到 n，而 Fortran 程序运行的维度是从 1 到 n。

• 由于固定的列索引，Fortran 源从不定义（计算值）r(3:,:)。

鉴于这些差异，结果不同也就不足为奇了。

在考虑到浮点运算的常见变数之后，将 maple 源代码简单、相对直接地翻译成 F2008 会得到相同的结果。

module romberg_module
  implicit none

  integer, parameter :: rk = kind(1.0d0)

  abstract interface
    function f_interface(x)
      import :: rk
      implicit none
      real(rk), intent(in) :: x
      real(rk) :: f_interface
    end function f_interface
  end interface
contains
  function romberg(f, a, b, n) result(r)
    procedure(f_interface) :: f
    real(rk), intent(in) :: a
    real(rk), intent(in) :: b
    integer, intent(in) :: n
    real(rk) :: r(0:n,0:n)    ! function result.

    real(rk) :: h
    integer :: row
    integer :: col
    integer :: k

    h = b - a
    r(0,0) = h / 2 * (f(a) + f(b))
    do row = 1, n
      h = h / 2
      r(row, 0) = 0.5_rk * r(row-1, 0)  &
          + sum(h * [(f(a + (2 * k - 1) * h), k = 1, 2**(row-1))])
      do col = 1, row
        r(row, col) = (4**col * r(row, col-1) - r(row-1, col-1))  &
            / (4**col - 1)
      end do
    end do
  end function romberg

  subroutine print_table(unit, r)
    integer, intent(in) :: unit
    real(rk), intent(in) :: r(0:,0:)
    integer :: row
    do row = 0, ubound(r,1)
      write (unit, "(*(F13.10,1X))") r(row, :row)
    end do
  end subroutine print_table
end module romberg_module

program print_romberg_table
  use, intrinsic :: iso_fortran_env, only: output_unit
  use romberg_module
  implicit none
  real(rk), allocatable :: r(:,:)
  r = romberg(f, 0.0_rk, 1.0_rk, 4)
  call print_table(output_unit, r)
contains
  function f(x)
    real(rk), intent(in) :: x
    real(rk) :: f
    f = 1.0_rk / (1.0_rk + x**2)
  end function f
end program print_romberg_table