来自 scipy.optimize 的 fsolve 失败答案

【问题标题】：fsolve from scipy.optimize fails来自 scipy.optimize 的 fsolve 失败
【发布时间】：2021-12-17 08:33:36
【问题描述】：

我正在尝试计算向量值函数的根，但我从 fsolve 得到的解决方案是完全错误的。我确实收到了警告：175: RuntimeWarning: 迭代没有取得良好的进展，这是通过最近十次迭代的改进来衡量的。 warnings.warn(msg, RuntimeWarning)

下面是工作示例：


import numpy as np 
from scipy.optimize import fsolve


def AMOC(amoc_state,
         gamma= 1/0.34,
         theta = 1,
         mu = 7.5,
         sigma = 0.85):
    T = amoc_state[0]
    S = amoc_state[1]
    dT = -gamma * (T-theta) - T * (1+ mu*np.abs(T-S)) 
    dS = sigma-S*(1+mu*np.abs(T-S))
    return (dT, dS) 

test = fsolve(AMOC, (0.8,0.3), 2.3,xtol = 1e-7, factor = 0.1, maxfev = 1000)

这给了我test = array([0.57372733, 0.47660354])，这显然不是 AMOC 函数的根，因为

AMOC((0.57372733, 0.47660354), gamma = 2.3) 返回 (-0.01121948095040759, 0.026224895476309684)

非常感谢任何帮助！

【问题讨论】：

仅供参考：不同的初始猜测给出了很好的结果：test = fsolve(AMOC, (0.5, 1.0), 2.3, xtol=1e-7, factor=0.1, maxfev=1000)，AMOC(test, 2.3) 返回(1.454392162258955e-14, 1.1657341758564144e-14)。
一个可能的困难来源是您使用np.abs(T-S)。此函数在 T == S 处不可微分。 fsolve 是 MINPACK 函数 HYBRD 的包装器，它使用数值计算的雅可比行列式。在T == S 的点附近，雅可比将不正确。
除了Warren的评论，值得一提的是abs是不可微分的，所以这里要小心。
试试stackoverflow.com/questions/69759905/…中使用的方法
谢谢你们。这一切都说得通。我尝试了几个初始条件，但都没有接近解决方案。

标签： python scipy solver

【解决方案1】：

正如 cmets 中已经提到的，abs 函数是不可微分的，您的函数 AMOC 也是不可微分的。这一点很重要，因为fsolves 底层算法使用导数信息来求解F(x) = 0，因此假设函数F 是可微分的。

但是，作为一种解决方法，您可以将 abs(x) 替换为平滑近似 sqrt(x**2 + 1.0e-4)：

def AMOC(amoc_state, gamma= 1/0.34, theta = 1, mu = 7.5, sigma = 0.85):
    T, S = amoc_state
    smooth_abs = np.sqrt((T-S)**2 + 1.0e-4)
    dT = -gamma * (T-theta) - T * (1+ mu*smooth_abs) 
    dS = sigma-S*(1+mu*smooth_abs)
    return np.array((dT, dS))

接下来，您可以尝试 Warren 已经提到的另一个初步猜测。否则，您可以通过scipy.optimize.minimize 将问题解决为等效最小化问题。根据我的经验，求解器更加强大：

from scipy.optimize import minimize

res = minimize(lambda x: np.sum(AMOC(x, 2.3)**2), x0=(0.8,0.3))

产生解决方案array([2.25836418e-01, 1.38576534e-06])。

【讨论】：