Fortran 95中泰勒级数求和与内置sin函数值对比错误排查

我尝试用Fortran 95编写代码，对比泰勒级数求和的sin值与内置sin函数的计算结果，但泰勒级数的输出结果存在错误。以下是我的原代码：

program taylor_sin
    implicit none
    
    integer, parameter :: dp = selected_real_kind(15, 307)
    integer :: i
    real(dp) ::  x_rad , term, sum, fact, y, delta
    
   write(*,*) "Enter the value of x in degree:"
    read(*,*) y
        x_rad = y * (3.14159265358979323846 / 180.0)
    
    ! Initialize the sum and the factorial
    sum = 0.0_dp
    fact = 1.0_dp
    delta = 10.0E-8
    
    ! Calculate the Taylor series expansion
    do 
     i = 1  
        term = ((-1)**i) * (x_rad**(2*i + 1)) / fact
        fact  = fact*(2*i + 1)
        sum = sum + term
       if( (sum - sin(y) ) < delta)exit
        
        
    end do
    
    
    ! Print the result
  write(*,*) " from the Taylor series"
    write(*,*) "sin(", y, ") = ", sum

        ! Print the result
    write(*,*) "intrinsic value for sin x", sin(y)
    
end program taylor_sin

错误分析 #

• 循环变量重置错误：do循环内每次将i=1，导致循环永远只计算i=1的项，无法迭代后续项
• 泰勒级数起始项错误：sin(x)的泰勒级数展开式为 $\sin(x) = \sum_{i=0}^{\infty} \frac{(-1)^i x^{2i+1}}{(2i+1)!}$，原代码从i=1开始，跳过了首项x，且阶乘计算逻辑错误
• 终止条件错误：
  1. 内置sin函数的参数应为弧度值x_rad，而非角度值y
  2. 用sum - sin(y)判断终止不严谨，应该判断当前项的绝对值小于精度阈值delta，避免因累加误差导致提前退出或死循环
• 阶乘计算逻辑错误：原代码每次仅乘以2i+1，无法正确生成(2i+1)!的阶乘值

修正后的代码 #

program taylor_sin
    implicit none
    
    integer, parameter :: dp = selected_real_kind(15, 307)
    integer :: i
    real(dp) :: x_rad, term, sum, y, delta
    real(dp), parameter :: pi = 3.14159265358979323846_dp
    
    write(*,*) "Enter the value of x in degree:"
    read(*,*) y
    x_rad = y * (pi / 180.0_dp)  ! 角度转弧度
    
    ! 初始化：泰勒级数首项为x_rad（i=0时的项）
    sum = x_rad
    term = x_rad
    delta = 1.0E-12_dp  ! 匹配双精度的精度阈值
    i = 1
    
    ! 迭代计算泰勒级数后续项
    do
        ! 利用前一项推导当前项，避免直接计算高次幂和大阶乘（防止溢出）
        term = term * (-1.0_dp) * x_rad**2 / (real(2*i, dp) * real(2*i + 1, dp))
        sum = sum + term
        
        ! 当当前项的绝对值小于精度阈值时退出循环
        if (abs(term) < delta) exit
        i = i + 1
    end do
    
    ! 输出结果，内置sin函数使用弧度参数x_rad
    write(*,*) "From the Taylor series:"
    write(*,*) "sin(", y, ") = ", sum
    write(*,*) "Intrinsic sin value: ", sin(x_rad)
    
end program taylor_sin

修正说明 #

• 循环变量正确迭代：将i=1移到循环外，每次循环后i=i+1，实现逐项迭代
• 泰勒级数首项正确初始化：sum初始化为x_rad（对应i=0的项），后续项通过前一项推导，避免直接计算高次幂和大阶乘，同时防止数值溢出
• 终止条件修正：判断当前项的绝对值小于精度阈值，确保迭代到足够小的项才停止，同时内置sin函数使用转换后的弧度值x_rad
• 精度优化：将delta设为1.0E-12_dp，匹配双精度的精度范围，同时用pi参数定义圆周率，保证精度

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