如何求解非线性方程组

Question

我有一组可以用矩阵表示的方程组，我还要求对于解中的所有变量x_I（x_I）³-x_I=0。例如

A = [0 1 0 0]
    [0 0 1 0]
    [1 0 0 1]

我也有，不是所有的变量都是0。在

这意味着

x₂=0
x₃=0
x₁+x₄=0
（x₁）³-x₁=0
（x₂）³-x₂=0
（x₃）³-x₃=0
（x₄）³-x₄=0

一个简单的解决方案是x₁=1和x₄=-1。在

你怎么能解决像这样的方程组的小例子？最好是我想要一个至少可以从python调用的解决方案。在

我目前解决这个问题的方法是用-1,0,1的值尝试所有3ⁿ不同的向量。在

^{pr2}$

编辑

这应该是对@alko答案的评论，但它太长了。让我用一个例子来说明这个方法。在

A = np.matrix([[0, 1, 1, 1, 1, 0, 1], [1, 0, 1, 1, 1, 1, 0], [0, 1, 0, 1, 1, 1, 1], [1, 0, 1, 0, 1, 1, 1], [0, 1, 0, 1, 0, 1, 1]])
p,l,u=scipy.linalg.lu(A)
print u
[[ 1.  0.  1.  1.  1.  1.  0.]
 [ 0.  1.  1.  1.  1.  0.  1.]
 [ 0.  0. -1.  0.  0.  1.  0.]
 [ 0.  0.  0. -1.  0.  0.  1.]
 [ 0.  0.  0.  0. -1.  0.  0.]]

我不清楚下一步会是什么。。。？在

Answer 1

这种问题称为约束规划。有一些python库可以解决这个问题。例如，下面的代码使用or-tools

from constraint_solver import pywrapcp as cs
import numpy as np

A = np.array(
    [[0, 1, 0, 0],
     [0, 0, 1, 0],
     [1, 0, 0, 1]], np.bool)

#A = np.array( [[0, 1, 1, 1, 1, 0, 1],
#               [1, 0, 1, 1, 1, 1, 0], 
#               [0, 1, 0, 1, 1, 1, 1], 
#               [1, 0, 1, 0, 1, 1, 1], 
#               [0, 1, 0, 1, 0, 1, 1]], np.bool)

values = [-1, 0, 1]
solver = cs.Solver("name")
X = np.array([solver.IntVar(values) for i in range(A.shape[1])])
Xl = X.tolist()

for row in A:
    x = X[row].tolist()
    solver.Add(solver.Sum(x) == 0)

db = solver.Phase(Xl, 
                  solver.INT_VAR_DEFAULT, 
                  solver.INT_VALUE_DEFAULT)

solver.NewSearch(db)
while solver.NextSolution():
    solution = [x.Value() for x in Xl]
    print solution

输出：

^{pr2}$

Answer 2

您可以将A分解为^{} decomposition：

import numpy as np
import scipy.linalg
A = np.matrix([[0, 1, 1, 1, 1, 0, 1],
               [1, 0, 1, 1, 1, 1, 0], 
               [0, 1, 0, 1, 1, 1, 1], 
               [1, 0, 1, 0, 1, 1, 1], 
               [0, 1, 0, 1, 0, 1, 1]])
p, l, u = scipy.linalg.lu(A)

为了简单起见，我们假设A具有最大秩，也就是说，它的某些子类是可逆的。如果不是真的，您可以使用如下相同的LU分解将矩阵A缩减为更小的A'，A'满足此属性并在等式中等价（最后一行u将为零，删除它们并继续）。你的方程Ax=0等于(plu)x=0，asp和{}是可逆矩阵，ux=0。u，依次是上三角形

^{pr2}$

您只对last（两）列感兴趣，后者根据last variable value定义所有变量。探测所有x4个可能的非零值，1和-1您可以检查您的方程组是否有解决方案：

from itertools import product
idx = u.shape[0] - u.shape[1] 
m, v = u[:, :idx], u[:,idx:]
for spare in product({1, 0, -1}, repeat=-idx):
    if any(spare): # to exclude all zeros
        sol = scipy.linalg.solve(m, -np.dot(v, spare))
        if not set(sol).difference({0,1,-1}):
            print 'solution: ', sol, spare
# solution:  [-1.  1.  1. -1. -0.] (1, -1)
# solution:  [ 1. -1. -1.  1. -0.] (-1, 1)

Answer 3

既然你加上了SymPy标签，我会指出SymPy可以很容易地解决这个问题

In [6]: x1, x2, x3, x4 = symbols('x1:5')

In [7]: solve([x2, x3, x1 + x4, x1**3 - x1, x2**3 - x2, x3**3 - x3, x4**3 - x4], [x1, x2, x3, x4], dict=True)
Out[7]: [{x₁: -1, x₂: 0, x₃: 0, x₄: 1}, {x₁: 0, x₂: 0, x₃: 0, x₄: 0}, {x₁: 1, x₂: 0, x₃: 0, x₄: -1}]

如果你的解不是整数，事情会变得更棘手，因为根式的解可能不存在，或者至少SymPy可能找不到它们。如果是这样，并且您只对数值解感兴趣，那么您应该使用numpy或scipy之类的数值库，但是由于您的解都是-1、0或1，所以这不是问题所在。在

编辑：

如果你有矩阵，说：

^{pr2}$

然后把它转换成一个系统很容易。只需将其乘以包含符号的向量（为了方便起见，我在这里切换到基于0的索引）：

In [13]: syms = symbols("x:4")

In [14]: s = Matrix(syms)

In [15]: constraints = [xi**3 - xi for xi in syms]

In [16]: A*s
Out[16]:
⎡  x₁   ⎤
⎢       ⎥
⎢  x₂   ⎥
⎢       ⎥
⎣x₀ + x₃⎦

In [17]: solve(list(A*s) + constraints, syms, dict=True)
Out[17]: [{x₀: -1, x₁: 0, x₂: 0, x₃: 1}, {x₀: 0, x₁: 0, x₂: 0, x₃: 0}, {x₀: 1, x₁: 0, x₂: 0, x₃: -1}]

以下是您的大型系统的解决方案：

In [35]: A = np.matrix([[0, 1, 1, 1, 1, 0, 1], [1, 0, 1, 1, 1, 1, 0], [0, 1, 0, 1, 1, 1, 1], [1, 0, 1, 0, 1, 1, 1], [0, 1, 0, 1, 0, 1, 1]])

In [36]: A = Matrix(A).applyfunc(int)

In [37]: syms = symbols("x:7")

In [38]: s = Matrix(syms)

In [39]: constraints = [xi**3 - xi for xi in syms]

In [40]: solve(list(A*s) + constraints, syms, dict=True)
Out[40]:
[{x₀: -1, x₁: 1, x₂: 1, x₃: -1, x₄: 0, x₅: 1, x₆: -1}, {x₀: 0, x₁: 0, x₂: 0, x₃: 0, x₄: 0, x₅: 0, x₆: 0}, {x₀: 1, x₁: -1, x₂: -1, x₃: 1, x₄: 0, x₅: -1, x₆:
1}]

两个注意事项：

• 你不需要在陆上拿到。SymPy的解算将为您解决这个问题（如果您的计算时间太长，您可以开始担心这些问题）。

• In [36]将矩阵项转换为int（默认情况下它们是浮点型）。这是不必要的，但一般来说，当你知道精确数是精确数时，SymPy会做得更好，特别是因为你知道你的解无论如何都是整数。如果您从一开始就使用SymPyMatrix，您就不必担心这个问题。