如何在Python中求解具有未知成员乘法的非线性方程组？

Question

我有以下非线性方程组：

​

​

这与Cholesky的分解非常相似，但不幸的是，并不是这样。

我已经找到了非常熟悉的解决方案：用Python (scipy.optimize.fsolve)求解非线性方程组，但我不知道如何动态地设置系统的所有方程。

如何在Numpy、Scipy或Python中的任何其他包中解决这个系统？

Answer 1

SymPy是你要找的包。这里是一个动态设置方程的函数，如您的系统中所描述的：

import sympy as sp

def not_choleskys_decomposition(N):
    # Initialisation of list of variables.
    c = []
    K = []
    for n in range(N + 1):
        c.append(sp.Symbol("c_{}".format(n), real=True))
        K.append(sp.Symbol("K_{}".format(n), real=True))

    # Setup your N+1 equations.
    equations = []
    for n in range(N + 1):
        num_c = N + 1 - n
        c_terms = [c_i * c_j for c_i, c_j in zip(c[:num_c], c[n:][:num_c])]
        lhs = sp.Add(*c_terms)
        equations.append(sp.Eq(lhs, K[n]))
    return equations, c, K

打印用于验证，可以使用sympy.latex函数完成，该函数直接将数学转换为LaTeX或sympy.pprint函数。

>>> test, _, _ = not_choleskys_decomposition(4)
>>> for t in test:
...    sp.pprint(t)

c₀⋅c₀ + c₁⋅c₁ + c₂⋅c₂ + c₃⋅c₃ + c₄⋅c₄ = K₀
c₀⋅c₁ + c₁⋅c₂ + c₂⋅c₃ + c₃⋅c₄ = K₁
c₀⋅c₂ + c₁⋅c₃ + c₂⋅c₄ = K₂
c₀⋅c₃ + c₁⋅c₄ = K₃
c₀⋅c₄ = K₄

可以使用属于每个对象的.subs方法来替换变量：

>>> equations, c_sym, K_sym = not_choleskys_decomposition(4)
>>> dict_replace = {c_sym[0]: 1.0}
>>> new_equation = equations[0].subs(dict_replace)
>>> sympy.pprint(new_equation)

c₁⋅c₂ + 1.0⋅c₁ + c₂⋅c₃ + c₃⋅c₄ = K₁

，这里是一个示例函数，它取代了，您的c或K列表。

def substitute_sym(equations, sym_list, val_list):
    # Confirm all dimensions are consistent.
    try:
        assert len(equations) == len(sym_list) == len(val_list)
    except AssertionError:
        raise IndexError("Inconsistent dimensions.")

    # Replace c symbols with values in equations.
    substituted_equations = []
    for eq in equations:
        _eq = eq.subs(dict(zip(sym_list, val_list)))
        substituted_equations.append(_eq)

    return substituted_equations

在c上进行测试

>>> equations, c_sym, K_sym = not_choleskys_decomposition(4)
>>> test_2 = substitute_sym(test_1, c_sym, [1] * 5)
>>> for t in test_2:
...    sp.pprint(t)

5 = K₀
4 = K₁
3 = K₂
2 = K₃
1 = K₄

在K上进行测试

>>> equations, c_sym, K_sym = not_choleskys_decomposition(4)
>>> K_vals = [1, 1.4, 0.2, 0.5, 3.0]
>>> test_3 = substitute_sym(test_1, K_sym, K_vals)
>>> for t in test_3:
...    sp.pprint(t)

c₀⋅c₀ + c₁⋅c₁ + c₂⋅c₂ + c₃⋅c₃ + c₄⋅c₄ = 1
c₀⋅c₁ + c₁⋅c₂ + c₂⋅c₃ + c₃⋅c₄ = 1.4
c₀⋅c₂ + c₁⋅c₃ + c₂⋅c₄ = 0.2
c₀⋅c₃ + c₁⋅c₄ = 0.5
c₀⋅c₄ = 3.0

可以用sympy.solvers.solveset.nonlinsolve(system, *symbols) 记录在这里。来求解你的方程组

给定c的K的求解

>>> from sympy.solvers.solveset import nonlinsolve
>>> test_solve, c_sym, K_sym = not_choleskys_decomposition(1)
>>> test_replaced = substitute_sym(test_solve, K_sym, [sp.Rational(0.5)] * 2)
>>> solved = nonlinsolve(test_replaced, c_sym)
>>> sp.pprint(solved)

⎧⎛ ⎛                    2⎞                           ⎞  ⎛ ⎛                   
⎪⎜ ⎜       ⎛  √6   √2⋅ⅈ⎞ ⎟ ⎛  √6   √2⋅ⅈ⎞    √6   √2⋅ⅈ⎟  ⎜ ⎜       ⎛  √6   √2⋅ⅈ
⎨⎜-⎜-1 + 2⋅⎜- ── - ────⎟ ⎟⋅⎜- ── - ────⎟, - ── - ────⎟, ⎜-⎜-1 + 2⋅⎜- ── + ────
⎪⎝ ⎝       ⎝  4     4  ⎠ ⎠ ⎝  4     4  ⎠    4     4  ⎠  ⎝ ⎝       ⎝  4     4  
⎩                                                                             

 2⎞                           ⎞  ⎛ ⎛                  2⎞                      
⎞ ⎟ ⎛  √6   √2⋅ⅈ⎞    √6   √2⋅ⅈ⎟  ⎜ ⎜       ⎛√6   √2⋅ⅈ⎞ ⎟ ⎛√6   √2⋅ⅈ⎞  √6   √2⋅
⎟ ⎟⋅⎜- ── + ────⎟, - ── + ────⎟, ⎜-⎜-1 + 2⋅⎜── - ────⎟ ⎟⋅⎜── - ────⎟, ── - ───
⎠ ⎠ ⎝  4     4  ⎠    4     4  ⎠  ⎝ ⎝       ⎝4     4  ⎠ ⎠ ⎝4     4  ⎠  4     4 


 ⎞  ⎛ ⎛                  2⎞                       ⎞⎫
ⅈ⎟  ⎜ ⎜       ⎛√6   √2⋅ⅈ⎞ ⎟ ⎛√6   √2⋅ⅈ⎞  √6   √2⋅ⅈ⎟⎪
─⎟, ⎜-⎜-1 + 2⋅⎜── + ────⎟ ⎟⋅⎜── + ────⎟, ── + ────⎟⎬
 ⎠  ⎝ ⎝       ⎝4     4  ⎠ ⎠ ⎝4     4  ⎠  4     4  ⎠⎪
                                                   ⎭

如果您选择使用SymPy，我希望这有助于您开始使用它。我很少使用它的非线性求解器，所以我不能保证任何超过例子的东西。争论也是相关的，而解决问题的人也是情绪化的。您会注意到，对于每个索引，我都将K替换为sp.Rational(0.5)。这样做是为了避免在为某些解决程序定义float类型时抛出错误。祝好运。

编辑：

还请注意，您不需要在SymPy中使用求解器。我很少这么做。我使用中的符号数学包进行LaTeX和方程操作。

编辑：

使它与scipy.optimize.fsolve一起工作的

from sympy.utilities.lambdify import lambdify
from scipy.optimize import fsolve


def lambdify_equations_for_solving_c(equations, c_sym, K_sym, K_val):
    f = []
    equations_subbed = substitute_sym(equations, K_sym, K_val)
    for eq in equations_subbed:
        # Note that I reformulate equation into an expression here
        # for determining the roots.
        f.append(lambdify(c_sym, eq.args[0] - eq.args[1]))
    return f


N = 4
equations, c_sym, K_sym = not_choleskys_decomposition(N)
K_vals = [0.1, 0, 0.1, 0, 0]
f = lambdify_equations_for_solving_c(equations, c_sym, K_sym, K_vals)

def fsolve_friendly(p):
    return [_f(*p) for _f in f]

c_sol = fsolve(fsolve_friendly, x0=(0, 0.1, -1.1, 0, 0))

fsolve返回非收敛结果的警告。

/home/ggarrett/anaconda3/envs/sigh/lib/python3.7/site-packages/scipy/optimize/minpack.py:162: RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last five Jacobian evaluations.
  warnings.warn(msg, RuntimeWarning)

这是一个完整的答案，说明如何将您的方程作为python可调用函数来求解。解决程序本身就是另一个挑战。