编辑代码

program main
    implicit none
    
    integer, parameter :: n = 1000  ! 离散点数量
    real(kind=8) :: x_0, y_0, z_0, v, c, b
    real(kind=8) :: a, b, h, integral
    integer :: i
    real(kind=8), dimension(n+1) :: u, t, integrand
    
    ! 定义积分区间
    a = -1.0d0
    b = 1.0d0
    
    ! 计算步长
    h = (b - a) / real(n)
    
    ! 计算函数在每个点的值
    do i = 0, n
        u(i+1) = a + real(i) * h
        t(i+1) = a + real(i) * h
        integrand(i+1) = my_function(u(i+1), t(i+1)) ! 替换 my_function 为实际函数
    end do
    
    ! 执行复合梯形法积分
    integral = 0.5d0 * (integrand(1) + integrand(n+1))
    do i = 1, n-1
        integral = integral + integrand(i+1)
    end do
    integral = integral * h
    
    ! 输出结果
    print *, "Integral =", integral
    
contains

    ! 定义函数
    real(kind=8) function my_function(u, t)
        real(kind=8), intent(in) :: u, t
        
        ! 插入您的函数表达式
        my_function = (z_0*(u-x_0)*(x-u))/((u-x_0)**2+(u-y_0)**2+z_0**2)**(5.0/2.0) * &
                      (-(3*(x-u)**2)/(sqrt((x-u)**2+(y-t)**2+z**2))**2 + &
                      ((1-2*v))/(sqrt((x-u)**2+(y-t)**2+z**2)+z)**2 * &
                      ((sqrt((x-u)**2+(y-t)**2+z**2)+z)**2-3*(sqrt((x-u)**2+(y-t)**2+z**2))**2+ &
                      ((x-u)**2*(3*sqrt((x-u)**2+(y-t)**2+z**2)+z))/((sqrt((x-u)**2+(y-t)**2+z**2)+z)))
    end function my_function
    
end program main