如何在Fortran中嵌套调用积分函数求解二重积分？

核心思路 #

二重积分$\int_{x=a}^{x=b} \int_{y=c}^{y=d} h(x,y) dy dx$可拆分为两步：先固定x，对y做一维积分得到关于x的函数，再对x做一维积分。利用已有的INTEGRAL2函数，通过嵌套调用即可实现：

1. 定义中间函数，接收x作为参数，内部调用INTEGRAL2完成对y的积分
2. 外层调用INTEGRAL2，将中间函数作为被积函数，完成对x的积分

关键：用Module共享信息 #

Fortran的外部函数无法直接携带额外参数，因此必须通过Module共享目标函数h(x,y)的定义，以及y方向的积分上下限、子区间数，避免修改原INTEGRAL2函数。

完整代码示例 #

! 定义模块，共享h(x,y)、y方向积分参数
module integral_utils
    implicit none
    double precision :: y_low, y_high  ! y的积分上下限
    integer :: n_y                     ! y方向的子区间数
contains
    ! 定义目标被积函数h(x,y)，替换成你的实际表达式即可
    double precision function h(x, y)
        double precision, intent(in) :: x, y
        ! 示例表达式：h(x,y) = x*y + sin(x+y)
        h = x*y + sin(x + y)
    end function h

    ! 中间函数：固定x，对y积分h(x,y)
    double precision function integral_y(x)
        double precision, intent(in) :: x
        integral_y = INTEGRAL2(func_y, y_low, y_high, n_y)
    contains
        ! 内部函数：将h(x,y)转为仅关于y的函数（x已固定）
        double precision function func_y(y)
            double precision, intent(in) :: y
            func_y = h(x, y)
        end function func_y
    end function integral_y
end module integral_utils

! 原INTEGRAL2函数，保持完全不变
!----------------------------------------------------------------------------
double precision function INTEGRAL2(FUNC,a,b,N)
!----------------------------------------------------------------------------
! ***  numerical integration (Simpson-rule) with equidistant spacing      ***
!----------------------------------------------------------------------------
  implicit none
  double precision,external :: FUNC                    ! the function to be integrated
  double precision,intent(in) :: a,b                   ! boundary values
  integer,intent(in) :: N                    ! number of sub-intervals
  double precision :: dx,x1,x2,xm,f1,f2,fm,int         ! local variables
  integer :: i
  dx  = (b-a)/DBLE(N)                        ! x subinterval
  x1  = a                                    ! left   
  f1  = FUNC(a)
  int = 0.d0
  do i=1,N
    x2  = a+DBLE(i)*dx                       ! right
    xm  = 0.5d0*(x1+x2)                      ! midpoint
    f2  = FUNC(x2)
    fm  = FUNC(xm)
    int = int + (f1+4.d0*fm+f2)/6.d0*(x2-x1) ! Simpson rule
    x1  = x2
    f1  = f2                                 ! save for next subinterval
  enddo
  INTEGRAL2 = int
end function INTEGRAL2

! 主程序：调用嵌套积分求解二重积分
program double_integral_demo
    use integral_utils
    implicit none
    double precision :: x_low, x_high, result
    integer :: n_x

    ! 设置积分区间与子区间数（建议n_x、n_y取偶数，符合辛普森法则要求）
    x_low = 0.0d0    ! x的积分下限
    x_high = 1.0d0   ! x的积分上限
    y_low = 0.0d0    ! y的积分下限
    y_high = 2.0d0   ! y的积分上限
    n_x = 100        ! x方向子区间数
    n_y = 100        ! y方向子区间数

    ! 外层调用INTEGRAL2，对x积分，被积函数为integral_y
    result = INTEGRAL2(integral_y, x_low, x_high, n_x)

    ! 输出结果
    write(*,*) "二重积分结果：", result
end program double_integral_demo

代码说明 #

1. Moduleintegral_utils：
   • 存储y方向积分的参数，让中间函数integral_y可直接访问
   • 定义目标函数h(x,y)，替换为你实际需要计算的函数表达式即可
   • 中间函数integral_y通过内部函数func_y将二元函数转为一元函数，再调用INTEGRAL2完成y方向积分
2. 原INTEGRAL2函数：完全保留，无需任何修改
3. 主程序：设置积分区间与子区间数，外层调用INTEGRAL2完成x方向积分，最终得到二重积分结果

注意事项 #

• 辛普森法则要求子区间数N为偶数，建议确保n_x和n_y取偶数，否则会影响计算精度
• 若y的积分上下限是x的函数（如$y=c(x)$到$y=d(x)$），只需在integral_y函数内动态设置y_low和y_high即可，模块依然适用

内容的提问来源于stack exchange，提问作者Cameron Duff