0

问题

我正在尝试使用 数值求解代数方程的非线性系统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))
4

1 回答 1

1

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

从你原来的方程组开始

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

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

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令人满意的东西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]]
于 2019-08-15T22:28:55.550 回答