【问题标题】：how to avoid runtime errors with scipy.optimize.fsolve如何使用 scipy.optimize.fsolve 避免运行时错误
【发布时间】：2019-08-16 06:34:39
【问题描述】：

问题

我正在尝试使用scipy.optimize.fsolve 对代数方程的非线性系统进行数值求解。

我解决了系统的几个不同参数值（下面的k1, k2, k3）。对于参数fsolve 的某些值，会找到正确的解决方案，而对于其他参数值，则会出现以下警告

RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last five Jacobian evaluations.
  warnings.warn(msg, RuntimeWarning)

从这些情况下的结果来看，很明显出现了问题，因为h(result) 不是零向量。毫无疑问，解决方案确实存在，并且与找到正确解决方案的情况在性质上没有什么不同。

在这些情况下通常推荐什么？是初始条件的问题吗？

示例

下面我将展示我如何求解方程组的想法：

import numpy as np
from scipy.optimize import fsolve

# Parameters for the system of equations
k1, k2, k3 = 2., 4.5, 0.1

# Function for evaluating the residual of the system
def h(x):
    x1, x2, x3=x
    eqn1 = k1*x1 + k2*x2**2 - x3
    eqn2 = x1 - 2.*x2 + k3*x3 
    eqn3 = x1 + x2 + x3 - 1.
    return eqn1, eqn2, eqn3

# An initial guess
x0 = np.array([1., 0.5, 1.])

# Get the solution of the system of equations
result = fsolve(h, x0=x0, xtol=1e-5)

# Now, verify that the solution is correct
print(h(result)) # Should be near (0,0,0)

这有时很好用，但对于 k1, k2, k3 的某些值，它会引发上面讨论的 RuntimeWarning 并返回错误的结果

# Bad parameters
k1, k2, k3 = 2., 4.5, -1.

# Bad result
result = fsolve(h, x0=x0, xtol=1e-5)

# Verification shows the residual is not near (0,0,0)
print(h(result))

【问题讨论】：

您在这里忘记了操作员：eqn1=x1 k1+k2*x2**2-x3
这里有语法错误：eqn3=x1+x2+x3=1.
如果开始条件不好，通常数值方法会给你不好/没有结果。因此，通常您会通过选择更好的初始值来获得更好的结果。自动实现该目标的一种方法是在循环中使用不同的初始参数调用 fsolve 两次，并检查哪些结果使用从 0 的较小导数求解您的方程。
@3sm1r 你能分享一个k1,k2,k3 的例子，你遇到了你提到的RuntimeWarning 吗？ k1,k2,k3=np.array([2.,4.5,0.1]) 对我来说运行良好。
@Moormanly Ty 为您提供帮助。我已根据您的建议编辑了问题。为了清楚起见，我只想强调，问题不在于解决方案不存在。我知道它确实可以，但算法没有找到它。

标签： python scipy scipy-optimize

【解决方案1】：

我不知道如何使用fsolve 摆脱您的RuntimeWarning。如果您不介意使用一些代数来求解方程组，那么您可以这样做。

从你原来的方程组开始

def h(x):
    x1, x2, x3=x
    eqn1 = k1*x1 + k2*x2**2 - x3
    eqn2 = x1 - 2.*x2 + k3*x3 
    eqn3 = x1 + x2 + x3 - 1.
    return eqn1, eqn2, eqn3

矩阵化

mat = np.array([
    [ k1, 0, k2, -1  ],
    [ 1, -2, 0,   k3 ],
    [ 1,  1, 0,   1  ]
], dtype=float)

vec = np.array([0, 0, 1], dtype=float)

def h(x):
    x1, x2, x3=x
    return mat @ (x1, x2, x2**2, x3) - vec

在mat的零空间中找到一个4d向量y求解mat @ y = vec以及一个向量ns

y = np.linalg.lstsq(mat, vec, rcond=None)[0]

from scipy.linalg import null_space
ns = null_space(mat).flatten()

现在使用z = y + a * ns 也将解决mat @ z = vec 的任何标量a 的事实。这让我们可以找到一些满足z[1]**2 = z[2] 的z。对于这样的z，我们有x = z[[0,1,3]] 作为原始方程组的解。

将z = y + a * ns 插入约束z[1]**2 = z[2]，我们需要(y[1] + a * ns[1])**2 = y[2] + a * ns[2]。求解a 的二次方程得到

a = np.roots([ns[1]**2, 2 * y[1] * ns[1] - ns[2], y[1]**2 - y[2]])

由于您声称原始方程组有解，a 应该具有真实值。

# This assertion will fail if (and only if) the original system has no solutions.
assert np.isreal(a).all()

最后，我们可以得到原方程组的两个解

soln1 = (y + a[0] * ns)[[0,1,3]]
soln2 = (y + a[1] * ns)[[0,1,3]]

【讨论】：