为什么在 Fortran 上进行黎曼和近似时中点法则比辛普森法则更准确

Question

大家。

我只是在使用中点法则和辛普森法则从 [1, 2] 计算 x^2 的积分。而且我发现，在子区间数相同的情况下，中点法则近似似乎比辛普森法则近似更准确，这真的很奇怪。

中点法则近似源码为：

program midpoint
implicit none                 ! Turn off implicit typing
Integer, parameter :: n=100   ! Number of subintervals
integer :: i                  ! Loop index
real :: xlow=1.0, xhi=2.0     ! Bounds of integral
real :: dx                    ! Variable to hold width of subinterval
real :: sum                   ! Variable to hold sum
real :: xi                    ! Variable to hold location of ith subinterval
real :: fi                    ! Variable to value of function at ith subinterval
dx = (xhi-xlow)/(1.0*n)       ! Calculate with of subinterval
sum = 0.0                     ! Initialize sum
xi = xlow+0.5*dx              ! Initialize value of xi
do i = 1,n,1                  ! Initiate loop
 ! xi = xlow+(0.5+1.0*i)*dx
 write(*,*) "i,xi ",i,xi      ! Print intermidiate result
 fi = xi**2                   ! Evaluate function at ith point
 sum = sum+fi*dx              ! Accumulate sum
 xi = xi+dx                   ! Increment location of ith point
end do                        ! Terminate loop
write(*,*) "sum =",sum
stop                          ! Stop execution of the program
end program midpoint

对应的执行是：

  ......                   .....      ..................  
 i,xi          100   1.99499905    
 sum =   2.33332348

辛普森法则近似源码为：

program simpson
implicit none                 ! Turn off implicit typing
integer, parameter :: n=100   ! Number of subintervals
integer :: i=0                ! Loop index
real :: xlow=1.0, xhi=2.0     ! Bounds of integral
real :: h                     ! Variable to hold width of subinterval
real :: sum                   ! Variable to hold sum
real :: xi                    ! Variable to hold location of ith subinterval
real :: fi                    ! Variable to value of function at ith subinterval
real :: Psimp                 ! Variable of simpson polynomial of xi interval
h = (xhi-xlow)/(1.0*n)        ! Calculate width of subinterval
sum = 0.0                     ! Initialize sum
do while (xi<=xhi-h)          ! Initiate loop
 xi = xlow+i*2.0*h            ! Increment of xi
 i=i+1
 write(*,*) "i,xi ",i,xi      ! Print intermidiate result
 Psimp=xi**2+4.0*(xi+h)**2+(xi+2.0*h)**2
                              ! Evaluate function at ith point
 sum = sum+(h/3.0)*Psimp      ! Accumulate sum
end do                        ! Terminate loop
write(*,*) "sum =",sum
end program simpson

对应的执行是：

 ........                  ......    ...................  
 i,xi          101   2.00000000    
 sum =   2.37353396

为了获得与中点结果相同的数字精度，我必须将 Simpson 程序中的子区间数设置为 100000，这是中点程序的 1000 倍（我最初将两个数字子区间设置为 100 )

我检查了 Simpson 程序中的代码，但找不到问题所在。

如果我没记错的话，辛普森法则应该比中点法则收敛得更快。

Answer 1

Craig Burley 曾经说过，WHILE 循环看起来就像一旦违反循环的前提，循环就会立即退出。这里循环的前提在 x=xhi 时被违反，但循环不会在那个点中断，只有当整个 nother 迭代完成并且可以在循环顶部应用测试时。您可以更一致地使用 Fortran 习语将循环转换为计数 DO 循环，例如

DO i = 0, n/2-1

然后注释掉

  i=i+1

行。或者干脆在修改xi:

后立即测试循环前提

xi = xlow+i*2.0*h ! Increment of xi
if(xi>xhi-h) exit ! Test loop premise

根据辛普森规则，任何一种方法都会导致次数不高于 3 的多项式的预期结果。

为什么在 Fortran 上进行黎曼和近似时中点法则比辛普森法则更准确

Why midpoint rule turns out more accurate than Simpson's rule when doing riemann sum approximation on Fortran

fortran

gfortran