如何使 sympy 在求解高阶方程时只求解 REAL 结果？

Question

• 嗨，那里。当我使用 sympy 求解高阶非线性方程时，我发现获取所有结果太慢了。我想这是因为 sympy 也解决了complex的结果，事实是我只需要REAL的结果。那么，请问如何使 sympy 只求解 REAL 域中的结果，还是应该使用其他模块来求解这些方程？

• 这是代码。

from sympy import *
from sympy.solvers.solveset import nonlinsolve


h_0 = 170
b = 120
Es = 205
h = 200
M = 19760400
fc = 83
fy = 585


def cal(para):
    Ec_f, As_f, ft_f = para
    vp_c, vp_s, x_n = symbols('vp_c, vp_s, x_n', real=True)
    f1 = vp_c / vp_s - x_n / (h_0 - x_n)
    f2 = 1 / 2 * Ec_f * vp_c * b * x_n - As_f * Es * vp_s - ft_f * b * (h - x_n)
    f3 = M - As_f * Es * vp_s * (h_0 - x_n / 3) - ft_f * b * (h - x_n) * (h / 2 + x_n / 6)
    re = nonlinsolve([f1, f2, f3], [vp_c, vp_s, x_n])
    return re

para = [42.4, 226, 3]
result = cal(para)

• 上面的代码在我的 i5-5500U 笔记本电脑上大约需要 20 秒。因为我需要很多循环，所以加快 sympy 的速度很重要。

• 这些是第一个循环的结果。

    {(-0.469378305833599, 1.02614290514871, -143.317861965125), (1.00656308634739, 2.01548035330581, 56.6225231688573), ((-183.969547582998 - 351.920528939668*I)*(0.00608328163162871 - 0.00630805832877277*I + 4.48116677633078e-8*(-183.969547582998 - 351.920528939668*I)**2)/(2.44289841241567 + 1.15340253748558e-5*(-183.969547582998 - 351.920528939668*I)**2 + 2.76016101129151*I), -0.429603644119307*(-1.36072460310392 + 0.690040252822878*I)*(1.03415787737688 - 1.07236991589137*I + 7.61798351976232e-6*(-183.969547582998 - 351.920528939668*I)**2), -183.969547582998 - 351.920528939668*I), ((-183.969547582998 + 351.920528939668*I)*(0.00608328163162871 + 4.48116677633078e-8*(-183.969547582998 + 351.920528939668*I)**2 + 0.00630805832877277*I)/(2.44289841241567 - 2.76016101129151*I + 1.15340253748558e-5*(-183.969547582998 + 351.920528939668*I)**2), -0.429603644119307*(-1.36072460310392 - 0.690040252822878*I)*(1.03415787737688 + 7.61798351976232e-6*(-183.969547582998 + 351.920528939668*I)**2 + 1.07236991589137*I), -183.969547582998 + 351.920528939668*I)}

Answer 1

加快求解速度的最佳方法可能是更改问题的表示形式：

1. 使每个方程都为多项式（适用于 f1），
2. 使 Rational 实例中的每个数字
3. 导出一个Groebner Basis Representation

我为要演示中间结果，请参阅 https:// nbviewer.org/github/cknoll/demo-material/blob/main/symbolic/groebner_basis_demo.ipynb

这是代码（没有中间结果）：

from sympy import *
from sympy.solvers.solveset import nonlinsolve

import sympy as sp

def clean_numbers(expr):
    
    numbers = expr.atoms(sp.Number)
    rplmts = [(n, sp.Rational(n)) for n in numbers]
    return expr.subs(rplmts)


h_0 = 170
b = 120
Es = 205
h = 200
M = 19760400
fc = 83
fy = 585


Ec_f, As_f, ft_f = (42.4, 226, 3)
vp_c, vp_s, x_n = xx = symbols('vp_c, vp_s, x_n', real=True)


f1 = vp_c * (h_0 - x_n) - x_n * vp_s # ← converted to a polynomial equation


f2 = 1 / 2 * Ec_f * vp_c * b * x_n - As_f * Es * vp_s - ft_f * b * (h - x_n)
f3 = M - As_f * Es * vp_s * (h_0 - x_n / 3) - ft_f * b * (h - x_n) * (h / 2 + x_n / 6)

ff = clean_numbers(Matrix([f1, f2, f3]).expand())

# reformulate the polynomial equations in terms of a groebner basis
# this set of equations has the same set of solutions but is described by different (often easier) equations
# note that this description depends on the (lexical) order of the variables
gb = groebner(ff, xx, order="lex")

# convert to matrix for easier access
ff2 = Matrix(gb.args[0]) ##:

# solve the last equation (only depends on one variable)

sol2 = solve(ff2[-1], xx[-1]) ##

sol = []
for s in sol2:
    
    tmp_eqns = ff2[:2, :].subs(xx[-1], s.evalf())
    tmp_sol = solve(tmp_eqns, xx[:2]) # → dict
    tmp_sol[xx[-1]] =  s.evalf()
    sol.append(tmp_sol)

print(sol)

→

[{vp_c: 1.00656308634739, vp_s: 2.01548035330581, x_n: 56.6225231688573},
{vp_c: -0.469378305833599, vp_s: 1.02614290514871, x_n: -143.317861965125}]

大约需要1s。

Answer 2

如果您可以将系统简化为需要求解的单个方程，则可以使用real_roots，这应该非常快。

我假设您的循环将提供不同的 para_guessed 值。现在我们将这些称为 x、y、z，并根据这些变量得到方程组。我将对 cal 中对“nonlinsolve”的调用替换为 return ([f1, f2, f3], [vp_c, vp_s, x_n])

>>> from sympy.abc import x, y, z
>>> eqs, v = cal((x,y,z))

最后两个方程很容易求解前两个变量：

>>> v2 = solve(eqs[-2:], v[:2], dict=True)

现在将它们代入第一个方程——只有一个解，所以我们只需要处理一个可能的方程：

>>> xeq = eqs[0].xreplace(v2[0])

为了做好准备，让我们现在得到分子：

>>> top = xeq.as_numer_denom()[0]

现在...插入您的值 - 或您感兴趣的任何其他值 - 并获得该单个变量的解决方案：

>>> reps = dict(zip((x,y,z),(42.4, 226, 3)))
>>> e = xeq.subs(reps)
>>> xns = real_roots(e)

Backsubstitute 以获得其他变量：

>>> for i in xns:
...     s = {i:j.n() for i,j in Dict(v2[0]).xreplace(reps).xreplace({v[-1]: i}).items()}
...     s.update({v[-1]: i.n()})
...     print(s)

给出

{vp_s: 1.02614290514871, vp_c: -0.469378305833599, x_n: -143.317861965125}
{vp_s: 2.01548035330581, vp_c: 1.00656308634739, x_n: 56.6225231688573}

所以回顾一下：您使用 CAS 来减少您的系统 象征性地为一个单一方程，在代入已知值后可以很快求解出实根。然后用它来确定其他方程的值，这些方程已经很容易用符号求解以获得完整的数值解。没有虚部......而且很快。