Fortran中一维数组的积分

Question

任务是编写计算多项式函数积分的代码。我所附图像中显示的函数id。我写了代码并编译了，答案就出来了。但是，它与解析解完全不同。代码：

     program rectangularApproximation
     
     write(*,*) "Input values of a ,b and eps"
     read(*,*) a,b,eps



1    continue
     n=1000
     h=(b-a)/n
     s=0.0
     do i=1,n
         x=a+h*i
         s=s+f(x)*h
     enddo

     sprev=s
     n=10*n

     h=(b-a)/n

     s=0.0

     do i=1,n
         x=a+h*i
         s=s+f(x)*h
     enddo

     snext=s

     if (abs(sprev-snext)<eps) then
         write(*,*) snext,n
         stop
     end if
         goto 1

     write(*,*) s
     end
       real function f(x)

            implicit none
            real, intent(in)    :: x
            integer             :: i
            real, dimension(8) :: numbers
            numbers = (/1,3,1,4,2,3,0,1 /)


            do i = 1,8
                f = f + numbers(i) * x**(numbers(i))
            end do
        end function

运行代码得到的结果为588189248（区间（a,b）为（1,2），我选择epsillon=0.001）解析解如下：

解析解的答案是 169.256 。我的代码可能出了什么问题？

Answer 1

你的多项式代码是错误的：

            do i = 1,8
                f = f + numbers(i) * x**(numbers(i))
            end do

应该是这样的

            do i = 1,8
                f = f + numbers(i) * x**i
            end do