Fortran中用割线法求解多输入非线性方程遇除零错误求助

直接触发除零错误的原因 #

你的secant子程序存在变量屏蔽问题：形参已经声明了mu和u，但局部变量里又重复声明了REAL*8:: x_next, dx, derivate, mu, u，导致传入的mu(i)和int_energy(i)被未初始化的局部变量覆盖。未定义的mu和u会让函数f计算出异常值，进而导致derivate（导数近似值）为0，后续执行x_next = x0-f(x0, mu, u)/derivate时触发除零错误。

修复方法：删除局部变量声明中的mu和u，修改为：

REAL*8:: x_next, dx, derivate

核心逻辑矛盾：mu与T的依赖关系 #

你提到mu是温度T的非线性函数，但当前f函数将mu作为独立参数传入，这意味着迭代求解T时，mu被当成了常量，完全违背了你的实际需求。这会导致f(T)退化为线性函数，不仅无法正确求解，还可能引发计算异常。

修复方法：让f函数根据输入的T，从已有的mu数组中通过插值获取对应值（假设你有存储温度点的数组T_mu(:)和对应mu值的数组mu_vals(:)）：

REAL*8 FUNCTION f(T, u)
  IMPLICIT NONE
  REAL*8, INTENT(in):: T, u
  REAL*8, PARAMETER:: mp = 1.67e-24, gammaa = 1.666666666d0, k_boltz = 1.380649e-16
  REAL*8:: mu
  INTEGER:: idx
  REAL*8:: frac

  ! 线性插值获取当前T对应的mu（假设T_mu是升序排列）
  DO idx = 1, SIZE(T_mu)-1
    IF (T >= T_mu(idx) .AND. T <= T_mu(idx+1)) EXIT
  END DO
  ! 处理T超出数组范围的边界情况
  IF (idx == SIZE(T_mu)) THEN
    mu = mu_vals(idx)
  ELSE
    frac = (T - T_mu(idx)) / (T_mu(idx+1) - T_mu(idx))
    mu = mu_vals(idx) + frac*(mu_vals(idx+1)-mu_vals(idx))
  END IF

  f = T - ((gammaa - 1) * ((u * mu * mp) / k_boltz))
END FUNCTION f

修正割线法实现 #

你当前的代码不是标准割线法，而是固定导数的简化牛顿法，对于非线性函数（因mu依赖T，f(T)是非线性的），固定导数会导致迭代不收敛或异常。标准割线法需要两个初始点，每次迭代用两点的函数值更新斜率：

SUBROUTINE secant(f, x0, x1, n_max, accuracy, u, res)
  IMPLICIT NONE
  REAL*8, EXTERNAL:: f
  REAL*8, INTENT(in):: x0, x1, accuracy, u
  INTEGER, INTENT(in):: n_max
  REAL*8, INTENT(out):: res(2)
  REAL*8:: x_prev, x_curr, x_next, f_prev, f_curr, slope
  INTEGER:: counter

  x_prev = x0
  x_curr = x1
  f_prev = f(x_prev, u)
  f_curr = f(x_curr, u)

  ! 检查初始点是否满足精度
  IF (ABS(f_curr) < accuracy) THEN
    res(1) = x_curr
    res(2) = f_curr
    RETURN
  END IF

  counter = 0
  DO
    counter = counter + 1
    IF (counter > n_max) THEN
      PRINT*, "Too many iterations with no result."
      res = [0.0d0, 0.0d0]
      RETURN
    END IF

    ! 计算割线斜率，避免除零
    slope = (f_curr - f_prev) / (x_curr - x_prev)
    IF (ABS(slope) < 1e-15) THEN
      PRINT*, "Slope is zero, adjust initial points."
      res = [0.0d0, 0.0d0]
      RETURN
    END IF

    ! 迭代计算下一个点
    x_next = x_curr - f_curr / slope
    f_curr = f(x_next, u)

    ! 检查收敛条件
    IF (ABS(x_next - x_curr) < accuracy .OR. ABS(f_curr) < accuracy) THEN
      res(1) = x_next
      res(2) = f_curr
      RETURN
    END IF

    ! 更新迭代变量
    x_prev = x_curr
    x_curr = x_next
    f_prev = f_curr
  END DO
END SUBROUTINE secant

调用代码调整 #

需要传入两个初始点，并增加最大迭代次数（原5次太少，非线性迭代通常需要50-100次）：

DO i = 1, n
  CALL secant(f, 1.d0, 2.d0, 50, 1e-5d0, int_energy(i), temperature(i))
ENDDO

额外注意事项 #

• 确保T_mu数组是升序排列，插值逻辑才能正常工作；
• 处理T超出T_mu范围的边界情况，避免数组越界；
• 初始点尽量选择在解的附近，提升割线法的收敛概率。

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