使用fortran编写一段变步长的辛普森求积法的程序

Fortran是一种古老的高级编程语言，常用于数值计算，尤其是科学计算领域。下面是一个简单的例子，展示了如何使用Fortran编写一个变步长辛普森求积法（Simpson's Rule with variable step size）的程序。辛普森法则是一种用于近似一维函数曲线下的面积的数值积分方法，对于光滑函数效果良好。

program simpsons_rule
implicit none

real :: a, b, h, f1, f2, integral, step_size
integer :: i, n

! 输入积分区间和步长大小
print *, "Enter the interval (a, b):"
read(*), a, b
print *, "Enter initial step size:"
read(*), h

n = int((b - a) / h) ! 计算步数
integral = (h / 3) * (f(a) + 4*f((a+b)/2) + f(b)) ! 初始化积分值

do i = 1, n
    f1 = f(a)          ! 左端点函数值
    f2 = f(a+h)         ! 右邻节点函数值
    h = h / 2           ! 更新步长

    ! 辛普森求积公式
    integral = integral + (h / 3) * (2*f1 + 4*f2)
    
    a = a + h            ! 移动到下一个节点
end do

! 输出结果
print *, "Estimated integral using Simpson's rule is:", integral
print *, "Total steps taken:", n

contains
function f(x)
    real :: x, f
    f = x**2 * sin(x)   ! 示例函数，你可以替换为你需要积分的实际函数
end function f
end program simpsons_rule