如何使用Python及Sympy库获取非线性方程组的全部解

问题复现 #

尝试用Sympy求解如下二元非线性方程组时，仅得到5组解，与Maxima计算得到的10组总解数不符：

from sympy import *
x,y = symbols('x,y')
rea1 = (1.0*10**(-4)*(x+2)*(4*y+3*x+1)**3) - 1.3 * ((-2*y-x+1)*(-y-x+1)*(2*y+2*x+6)**2)
rea2 = (1.0*10**(-4)*(y+1)*(4*y+3*x+1)**4) - 2.99 * ((-2*y-x+1)**2*(-y-x+1)*(2*y+2*x+6)**2)
solu = solve([rea1,rea2],[x,y])
sol = nsolve([rea1,rea2],[x,y],[-0.1,1])
print(solu)
print(sol)

运行输出：

• solve()返回4组解析解：[(-11.0, 8.0), (-3.0, 2.0), (-2.0, -1.0), (5.0, -4.0)]
• 固定初始点的nsolve()返回1组数值解：Matrix([[1.32358494772278], [-0.324195048403443]])

Maxima计算得到的全部10组解如下：

[[x=-11,y=8],
[x=-3,y=2],
[x=5,y=-4],
[x=-2,y=-1],
[x=1.3236,y=-0.3242],
[x=-2.0091,y=-0.98836],
[x=-3.8143,y=0.84582],
[x=3.004,y=-1.0016],
[x=-4.0297,y=0.9959],
[x=-8.4744,y=9.4724]]

原因说明 #

Sympy的solve()默认仅返回具有解析闭式表达的有理根，不会自动计算无简单解析形式的数值根；而nsolve()是基于牛顿迭代的局部求解方法，单次调用仅能收敛到初始点所在收敛域内的1个根，无法自动遍历所有解。

可行实现方案 #

方案1：Groebner基消元后求全部多项式根 #

该方程组本质是二元高次多项式方程组，可通过Groebner基消元将二元方程组转化为单变量多项式，求出单变量多项式的所有根后回代，校验残差过滤增根即可得到全部解，参考代码如下：

from sympy import symbols, Poly, groebner, N
x, y = symbols('x y')
rea1 = (1e-4*(x+2)*(4*y+3*x+1)**3) - 1.3 * ((-2*y-x+1)*(-y-x+1)*(2*y+2*x+6)**2)
rea2 = (1e-4*(y+1)*(4*y+3*x+1)**4) - 2.99 * ((-2*y-x+1)**2*(-y-x+1)*(2*y+2*x+6)**2)
# 计算格罗布纳基完成消元
gb = groebner([rea1, rea2], [x, y])
# 取消元后仅含y的单变量多项式，求所有数值根
poly_y = Poly(gb[-1], y)
y_candidates = poly_y.nroots(n=8)
all_solutions = []
for y_val in y_candidates:
    # 代入y值，得到关于x的多项式，求x的所有根
    poly_x = Poly(rea1.subs(y, y_val), x)
    x_candidates = poly_x.nroots(n=8)
    for x_val in x_candidates:
        # 校验残差，过滤增根
        residual1 = abs(float(rea1.subs({x:x_val, y:y_val})))
        residual2 = abs(float(rea2.subs({x:x_val, y:y_val})))
        if residual1 < 1e-6 and residual2 < 1e-6:
            all_solutions.append( (round(x_val,5), round(y_val,5)) )
# 去重
unique_solutions = []
for sol in all_solutions:
    is_dup = False
    for exist in unique_solutions:
        if abs(sol[0]-exist[0]) < 1e-4 and abs(sol[1]-exist[1]) < 1e-4:
            is_dup = True
            break
    if not is_dup:
        unique_solutions.append(sol)
print(unique_solutions)

运行后即可输出全部10组实根。

方案2：多初始点网格搜索 #

如果不想处理多项式消元逻辑，可以根据已知解的范围划定搜索区间（本例中x范围约为[-12,6]，y范围约为[-2,10]），在区间内生成大量均匀分布的初始点，逐个调用nsolve迭代，收集所有收敛结果，经残差校验、去重后得到全部解。
注意：该方法需要足够密的初始点才能避免漏根，计算效率低于消元法。

方案3：同伦延拓法求解 #

针对多项式方程组全局求全根的场景，同伦延拓法是理论上可保证找到所有根的专业数值方法，求解效率远高于多初始点牛顿迭代，可直接调用对应实现完成求解。

内容的提问来源于stack exchange，提问作者Alby Stalks