如何修改Fortran实现的Simpson's Rule代码以获取正确积分结果？

嘿，我帮你找出代码里的问题啦！你尝试用Fortran实现辛普森法则计算函数 ( f(x) = x^5 + (x-2)\sin(x) + (x-1) ) 的积分，已知当 ( a=-1 )、( b=1 )、( m=20 ) 或40时正确结果应为-1.398，但当前代码输出离谱的7015869.0，问题出在两个关键地方：

1. 步长h的计算完全错误 #

辛普森法则要求把区间 ([a,b]) 分成m个等距子区间（m必须是偶数），正确的步长计算公式应该是：

h = (b - a) / real(m)

而你代码里写的是 h=(b-a)/2.0，这直接导致后续所有采样点的位置完全偏离原区间，数值计算爆炸，结果自然不对。

2. 循环内的采样点计算逻辑错误 #

原代码的循环逻辑没有正确对应辛普森法则的采样点顺序。辛普森法则的标准公式是：

( \int_a^b f(x)dx \approx \frac{h}{3} \left[ f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + ... + 2f(x_{m-2}) + 4f(x_{m-1}) + f(x_m) \right] )
其中 ( x_i = a + i \times h )，( i ) 从0到m，且m为偶数。

原代码里的点计算基于错误的h，同时循环内的奇偶项处理也有逻辑漏洞，比如重复计算或者遗漏部分点。

修正后的完整代码 #

program simpsons
    implicit none
    real a, b, h, integ, fa, fb, x
    integer i, m

    ! 输入区间边界
    write(*,*) 'enter the lower boundary'
    read(*,*) a
    write(*,*) 'enter the upper boundary'
    read(*,*) b
    do while(a.ge.b)
        write(*,*) 'reenter the lower boundary (must be less than upper)'
        read(*,*) a
        write(*,*) 'reenter the upper boundary'
        read(*,*) b
    enddo

    ! 输入区间数（要求偶数）
    write(*,*) 'enter the number of intervals (must be even)'
    read(*,*) m
    ! 额外判断m是否为偶数，避免后续逻辑出错
    do while(mod(m,2)/=0)
        write(*,*) 'reenter an even number of intervals'
        read(*,*) m
    enddo

    ! 计算正确步长
    h = (b - a) / real(m)

    ! 计算区间端点的函数值
    fa = a**5.0 + (a-2.0)*sin(a) + (a-1.0)
    fb = b**5.0 + (b-2.0)*sin(b) + (b-1.0)

    ! 初始化积分值
    integ = 0.0

    ! 遍历中间采样点，按辛普森法则累加奇偶项
    do i = 1, m-1
        x = a + real(i)*h
        if(mod(i,2)==1) then
            ! 奇数项（i=1,3,...m-1）乘以4
            integ = integ + 4.0 * (x**5.0 + (x-2.0)*sin(x) + (x-1.0))
        else
            ! 偶数项（i=2,4,...m-2）乘以2
            integ = integ + 2.0 * (x**5.0 + (x-2.0)*sin(x) + (x-1.0))
        endif
    enddo

    ! 应用辛普森法则的最终公式
    integ = (fa + fb + integ) * (h / 3.0)

    write(*,*) 'integration = ', integ
end program simpsons

修正说明 #

• 新增了m必须为偶数的判断，因为辛普森法则要求区间数是偶数，避免非法输入导致错误
• 修正了步长h的计算方式
• 重新设计了循环逻辑，直接遍历所有中间点，按奇偶性分别乘以4或2，更符合辛普森法则的标准实现
• 增加了变量x存储当前采样点，让代码更清晰易读

现在输入 ( a=-1 )、( b=1 )、( m=20 )，你会得到接近-1.398的正确结果啦！

内容的提问来源于stack exchange，提问作者Karim Abulazm