如何避免运行时错误scipy.optimize.fs文件

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的空空间中找到一个解mat @ y = vec的4d向量y，以及一个向量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也将为任何标量a求解mat @ 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]]