Scipy fsolve无法收敛至正确值，请求代码问题排查

现象描述 #

• 尝试用Scipy的fsolve求解微分方程$\frac{dh}{dt}$的初始值（t=0时的h），迭代值始终呈现如下模式：

1
1.0
1.0
1.0000000149011612
101.0
nan
1.0000000149011612
101.0
nan
nan
nan
nan
nan
nan
nan

• 无论调整P、T、kf或kp参数，迭代模式完全一致。
• 运行时触发警告：

RuntimeWarning: invalid value encountered in sqrt
  return 2*L*np.sqrt(R**2-(R-h)**2)

• 核心疑问：该异常是否由代码本身问题导致？

代码实现 #

# -*- coding: utf-8 -*-
"""
Created on Mon Apr  3 16:46:15 2023

@author: houle
"""

import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.optimize import fsolve

class parametre_antoine():
    A = 8.13484
    B = 1662.48
    C = 238.131
    mmhg_atm = 760
    
prm_antoine = parametre_antoine() 

rho = 0.8 #kg/L
Tin = 110 #C   
R = 1
kp = -200
kf = 300
M = 46.068/1000 #kg/mol
L = 5
M = 46.068 #g/mol
Vtot = np.pi*R**2*L

theta = [rho,R,kp,kf,M,L,prm_antoine,Vtot]

Fin = 20000
u = [Fin]

P = 3.3
T = 300

def Aire(R,L,h):
    return 2*L*np.sqrt(R**2-(R-h)**2)

def dv_dh(R,L,h):
    dVdh = R/np.sqrt(1-(h/R-1)**2)+L*np.sqrt((2*R-h)*h)-(L*(2*R-2*h)*(R-h))/(2*np.sqrt((2*R-h)*h))
    return dVdh

def dh_dt(h,theta,u,T,P):
    [rho,R,kp,kf,M,L,prm_antoine,Vtot] = theta
    [Fin] = u
    dhdt =(Fin- kp*(psat(prm_antoine,T)-P)*Aire(R,L,h)-kf*h**0.5)/rho/ dv_dh(R,L,h)
    print(h)
    return dhdt

def psat(prm,T):
    #T en Celcius
    #Valide entre -114 et 243C
    #Retourne la tension de vapeur en atm
    p_mmhg = 10**(prm.A-prm.B/(T+prm.C))
    p_atm = p_mmhg / prm.mmhg_atm
    return p_atm

x0 = [1]
u0 = [Fin]
x0 = fsolve(dh_dt, x0, args=(theta, u0,T,P))

问题分析与结论 #

该异常完全由代码本身的问题导致，核心原因及修复建议如下：

1. 函数定义域未约束，导致无效计算 #

• Aire函数的合法性问题：
  函数中np.sqrt(R**2-(R-h)**2)展开后为np.sqrt(2Rh - h²)，根号内的值必须≥0，即$0 ≤ h ≤ 2R$（当前R=1，h需在0~2之间）。但迭代中出现h=101.0，远超出定义域，导致根号内为负数，触发invalid value in sqrt警告，最终生成nan。
• dv_dh函数的同样问题：
  函数中多处用到np.sqrt((2*R-h)*h)，同样要求h在0~2R范围内，h=101时直接产生nan，导致dh_dt返回无效值，迭代彻底失控。

2. 代码中的其他隐患 #

• 变量重复赋值：代码中两次定义M，最终M被覆盖为46.068（g/mol），若后续逻辑需要kg/mol单位，会引发计算错误。
• dv_dh表达式冗余：(2*R-2*h)可简化为2*(R-h)，当前写法增加计算量且易出错。

3. 修复建议 #

• 添加定义域检查：在Aire和dv_dh中加入h的范围判断，超出范围时返回合理值（如0，或根据物理意义返回边界值）：

  def Aire(R,L,h):
      if h < 0 or h > 2*R:
          return 0
      return 2*L*np.sqrt(R**2-(R-h)**2)
  

• 限制fsolve搜索范围：改用scipy.optimize.root并传入bounds参数，约束h在0~2R之间，避免迭代跳出合法区间。
• 修正变量重复赋值：删除重复的M定义，统一单位。
• 简化dv_dh计算：优化表达式减少出错概率：

  def dv_dh(R,L,h):
      if h <= 0 or h >= 2*R:
          return 0
      sqrt_term = np.sqrt((2*R - h)*h)
      term1 = R / np.sqrt(1 - (h/R - 1)**2)
      term2 = L * sqrt_term
      term3 = (L * (R-h)**2) / sqrt_term
      return term1 + term2 - term3
  

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