Python scipy牛顿法求解原油密度迭代发散报错求助

核心错误是用错了scipy.optimize.newton的接口规则：

• 这个函数的作用是求解f(x) = 0的根，要求传入的目标函数返回值是「当前输入x对应的f(x)取值」，迭代过程会自动调整x让这个返回值趋近于0。
• 你现在写的目标函数返回的是「根据输入rho_p0计算出的更新后rho_p0值」，相当于让牛顿法去找让这个更新值等于0的解，完全偏离了要计算的不动点目标，第二次迭代就得到了负数rho_a，后续数值直接爆炸发散。

另外写在目标函数里的try-except RuntimeError是无效代码：迭代不收敛的异常是newton函数在外部抛出的，根本不会进入目标函数内部的异常捕获逻辑。

把目标函数改成残差形式即可：返回「根据输入rho_p0计算出的新rho_p0」和「输入rho_p0」的差值，当差值为0时，就满足需要的不动点关系，也就是正确的密度值。

修正后的可运行代码：

from scipy.optimize import newton

def mccain_hill_rhop0(rho_p0, Rs, Sg, So):
    a0 = -49.8930
    a1 = 85.0149
    a2 = -3.70373
    a3 = 0.0479818
    a4 = 2.98914
    a5 = -0.0356888
    rho_a = a0 + a1*Sg + a2*Sg*rho_p0 + a3*Sg*rho_p0**2 + a4*rho_p0 + a5*rho_p0**2
    rho_p0_calc = (Rs*Sg + 4600*So) / (73.71 + Rs*Sg/rho_a)
    return rho_p0_calc - rho_p0

Rs = 1000
Sg = 0.55
So = 0.8
rho_p0_init = 52.8 - 0.01*Rs
rho_p0 = newton(mccain_hill_rhop0, rho_p0_init, args=(Rs, Sg, So), tol=1e-3, maxiter=100)
print(rho_p0)

运行后可以正常收敛，得到结果约为41.51，不会再出现数值溢出的报错。

额外提示：如果只是要做不动点迭代，完全没必要调用牛顿法，自己写个while循环判断两次迭代的差值小于容差就可以终止。用scipy的求根函数必须先把方程整理成f(x)=0的标准形式。

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