Python数值求解非线性方程组异常：返回初始值问题排查

我尝试数值求解一个非线性方程组，先定义了函数F：

def F(x, na=1, nb=2, L=6, A=0.09, V=0.54, dh=0.03, k=10**-2, dB=0.005, fB=38, ur=0.39, rhoW=998, rhoG=1.2, g=9.81, pu=10**5, myW = 10**-3):
    # 变量提取
    uWa = x[0]
    uWb = x[1]
    ta = x[2]
    tb = x[3]
    ea = x[4]
    eb = x[5]
    rhoa = x[6]
    rhob = x[7]
    ua = x[8]
    ub = x[9]
    pa0 = x[10]
    pb0 = x[11]
    Rea = x[12]
    Reb = x[13]
    lambdaa = x[14]
    lambdab = x[15]

    # 构建非线性方程组（形式为f=0）
    f = np.zeros(16)
    f[0] = ta - L/(ur+uWa)
    f[1] = tb - L/(ur-uWb)
    f[2] = ea - na*ta*fB*np.pi*dB**3/(6*V)
    f[3] = eb - nb*tb*fB*np.pi*dB**3/(6*V)
    f[4] = rhoa - rhoW + ea*(rhoW-rhoG)
    f[5] = rhob - rhoW + eb*(rhoW-rhoG)
    f[6] = rhoa*ua - ea*rhoG*(ur+uWa) - (1-ea)*rhoW*uWa
    f[7] = rhob*ub + eb*rhoG*(ur-uWb) - (1-eb)*rhoW*uWb
    f[8] = (-rhoa*ua + rhob*ub + rhoG*(na+nb)*fB*np.pi*dB**3/(6*A**2))*A
    f[9] = pb0 - pa0 - rhoa*rhob*(ua**2-ub**2)/(rhoa+rhob)
    f[10] = pa0 - pu - rhoa*((ua**2*lambdaa*dh*np.pi)/(8*A)+g)*L    
    f[11] = pb0 - pu + rhob*((ub**2*lambdab*dh*np.pi)/(8*A)+g)*L
    f[12] = Rea - rhoa*ua*dh/myW
    f[13] = Reb - rhob*ub*dh/myW  
    f[14] = lambdaa - fixpunktiteration(colebrook, 0.05, args=(Rea,k,dh)) 
    f[15] = lambdab - fixpunktiteration(colebrook, 0.05, args=(Reb,k,dh))

    return f

为了计算f[14]和f[15]，我定义了两个辅助函数：

def colebrook(lam,args):
    Re, k, d = args
    return ( 1/(4*np.log10( 2.51/(Re*lam**0.5) + (k/d)/3.71 )**2) )

def fixpunktiteration(f,xalt,args,kmax=500,tol=1.e-7):
    for k in range(1,kmax):
        xneu = f(xalt,args)
        xdiff = xneu - xalt
        if abs(xdiff/xneu) < tol:
            break
        xalt = xneu
    else:
        xneu = None
    return xneu

之后我定义了所有参数和初始值，用scipy的fsolve求解：

xinit = np.array([uWa_0, uWb_0, ta_0, tb_0, ea_0, eb_0, rhoa_0, rhob_0, ua_0, ub_0, pa0_0, pb0_0, Rea_0, Reb_0, lambdaa_0, lambdab_0]) 
loesung = fsolve(func=F,x0=xinit)

程序运行没有报错，但返回的结果和初始值完全一致。如果把f[14]替换成f[14]=64/Rea，程序会返回和初始值不同的结果，但这个结果没有物理意义。我怀疑问题出在f[14]、f[15]或者colebrook函数里，但找不到具体位置。

问题分析与解决思路 #

1. 不动点迭代的核心问题

你的fixpunktiteration函数在fsolve迭代流程中存在两个致命缺陷：

• 固定初始迭代值：每次调用fixpunktiteration时，都用固定的0.05作为初始值，而非基于当前lambdaa/lambdab的迭代值。这意味着不管fsolve传入的lambdaa是什么，你都从0.05重新计算，导致f[14] = lambdaa - 固定值，相当于给lambdaa加了硬约束，fsolve会直接判定初始值满足约束，停止迭代。
• 迭代失败返回None：如果不动点迭代未收敛（进入else分支），返回None会导致f[14]/f[15]出现NaN，破坏fsolve的数值稳定性。

2. 修正方案

• 传递当前lambda值作为迭代初始值：把fsolve传入的lambdaa/lambdab作为fixpunktiteration的初始值，而非固定的0.05。这样fixpunktiteration的结果会跟随fsolve的迭代变化，让方程组约束生效。
• 处理迭代失败的情况：增加迭代失败后的 fallback 逻辑，比如返回最后一次迭代值，避免返回None。

修正后的代码片段：

# 修正F函数中的f[14]和f[15]
f[14] = lambdaa - fixpunktiteration(colebrook, lambdaa, args=(Rea,k,dh)) 
f[15] = lambdab - fixpunktiteration(colebrook, lambdab, args=(Reb,k,dh))

# 修正不动点迭代函数
def fixpunktiteration(f,xalt,args,kmax=500,tol=1.e-7):
    for k in range(kmax):
        xneu = f(xalt,args)
        if np.isnan(xneu):
            break  # 避免NaN扩散
        xdiff = xneu - xalt
        if abs(xdiff/xneu) < tol:
            return xneu
        xalt = xneu
    # 迭代未收敛时返回最后一次值并抛出警告
    import warnings
    warnings.warn("不动点迭代未收敛")
    return xalt

3. 额外排查点

• 检查初始值合理性：确保xinit中的Rea_0、Reb_0为正数，否则colebrook函数中的对数计算会出错（产生NaN）。
• 启用fsolve调试参数：开启full_output=True和verbose=1，查看迭代过程中的残差变化，确认是哪个方程导致迭代停滞：

  loesung, infodict, ier, mesg = fsolve(func=F,x0=xinit, full_output=True, verbose=1)
  print(mesg)
  print(infodict['fvec'])  # 查看每个方程的残差
  

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