如何避免使用 scipy.optimize.fsolve 的运行时错误

Question

问题

我正在尝试使用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))

Answer 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 @ y = vecy4D 向量以及mat零空间中的向量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也将解决任何标量amat @ z = vec的事实。这使我们能够找到一些z令人满意的z[1]**2 = z[2].对于这样的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]]