加速 SymPy 方程求解器
Speed up SymPy equation solver
我正在尝试使用以下 python 代码（当然使用 SymPy）来求解一组方程式：
def Solve(kp1, kp2):
a, b, d, e, f = S('a b d e f'.split())
equations = [
Eq(a+b, 2.6),
Eq(2*a + b + d + 2*f, 7),
Eq(d + e, 2),
Eq(a*e,kp2*b*d),
Eq( ((b * (f**0.5))/a)*((a+b+d+e+f+13.16)**-0.5), kp1)
]
return solve(equations)
代码成功求解了方程，但大约需要 35 秒。 Solve() 函数在另一个文件的迭代块（大约 2000 次迭代）中使用，所以速度对我来说真的很重要。
有什么方法可以加快求解器的速度吗？如果不能，您能否推荐另一种使用 python?
求解方程组的方法
你只需要解一次方程。之后，您将得到以下形式的方程式：
- a = f_1(kp1, kp2)
- b = f_2(kp1, kp2)
- ...
因此您可以根据 kp1 和 kp2 简单地计算 a, ..., e。例如求解第一个、第三个和第四个方程得到：
- b: -a + 2.6
- e: 2.0*kp2*(5.0*a - 13.0)/(5.0*a*kp2 - 5.0*a - 13.0*kp2),
- d: 10.0*a/(-5.0*a*kp2 + 5.0*a + 13.0*kp2)
在我的电脑上求解所有五个方程太慢，但如果它给你一个表达式，你只需插入（替换）kp1 和 kp2 的值，你不必再次求解方程。对于替换，请查看 sympy documentation。
所以你的循环应该是这样的：
solutions = sympy.solve(eqs, exclude=[kp1, kp2])
for data_kp1, data_kp2 in data:
for key, eq in solutions:
solution = eq.subs([(kp1, data_kp1), (kp2, data_kp2)])
final_solutions.append("{key}={solution}".format(key=key, solution=solution))
求解 a、b、d、e 的第 4 个线性方程。这会生成两个解决方案。将这些中的每一个代入第 5 个方程（以这种形式：Eq(b**2*f, d**2*kp1**2*(a + b + 2*d + f + 13.16))) 给出两个方程（e51 和 e52）。这两个方程在 f 中是非线性的，如果你对它们使用 unrad，你将得到 2 第六f 中的阶多项式——这可能是您得不到解决方案的原因，因为通常只有四次方程是完全可解的。将这两个方程称为 e51u 和 e52u。
如果您对实根感兴趣，可以使用 real_roots 为这些多项式 提供根，然后 将所需的 kp1 和 kp2 值代入只留下 f 作为未知数。例如，
>>> ans = solve(equations[:-1],a,b,d,e)
>>> l=equations[-1] # modified as suggested
>>> l=l.lhs-l.rhs
>>> l
b**2*f - d**2*kp1**2*(a + b + 2*d + f + 13.16)
>>> e51=l.subs(ans[0]); e51u=unrad(e51)[0]
>>> e52=l.subs(ans[1]); e52u=unrad(e52)[0]
>>> import random
>>> for r1,r2 in [[random.random() for i in range(2)] for j in range(3)]:
... print [i.n(2) for i in real_roots(e51u.subs(dict(kp1=r1,kp2=r2)))]
... print [i.n(2) for i in real_roots(e52u.subs(dict(kp1=r1,kp2=r2)))]
... print '^_r1,r2=',r1,r2
...
[1.7, 2.9, 3.0, 8.2]
[1.7, 2.9, 3.0, 8.2]
^_r1,r2= 0.937748743197 0.134640776315
[1.3, 2.3, 4.7, 7.4]
[1.3, 2.3, 4.7, 7.4]
^_r1,r2= 0.490002815309 0.324553144174
[1.1, 2.1]
[1.1, 2.1]
^_r1,r2= 0.308803300429 0.595356213169
看来e51u和e52u都给出了相同的解，所以你可能只需要使用其中一个。你应该检查原始方程式中的答案，看看哪些是真实的解决方案：
>>> r1,r2
(0.30880330042869408, 0.59535621316941589)
>>> [e51.subs(dict(kp1=r1,kp2=r2,f=i)).n(2) for i in real_roots(e51u.subs(dict(kp1=r1,kp2=r2)))]
[1.0e-12, 13.]
所以在这里你看到只有第一个解（从上面看到是1.1）实际上是一个解； 2.1 是一个虚假的解决方案。
您的方程组可以根据 kp1 和 kp2 作为符号参数转换为多项式系统。在这种情况下，将其转换为多项式系统是有利的，以便它可以部分求解，为参数的特定值的数值求解做准备。
她的技术涉及计算 Groebner 基数，您不希望在方程式中有任何浮点数，因此让我们首先确保使用精确的有理数（例如，使用 sqrt 而不是 **0.5 和 Rational('2.6') 而不是 2.6):
In [20]: kp1, kp2 = symbols('kp1, kp2')
In [21]: a, b, d, e, f = symbols('a, b, d, e, f')
In [22]: equations = [
...: Eq(a+b, Rational('2.6')),
...: Eq(2*a + b + d + 2*f, 7),
...: Eq(d + e, 2),
...: Eq(a*e,kp2*b*d),
...: Eq( ((b * sqrt(f))/a)/sqrt((a+b+d+e+f+Rational('13.16'))), kp1)
...: ]
In [23]: for eq in equations:
...: pprint(eq)
...:
a + b = 13/5
2⋅a + b + d + 2⋅f = 7
d + e = 2
a⋅e = b⋅d⋅kp₂
b⋅√f
─────────────────────────────── = kp₁
_________________________
╱ 329
a⋅ ╱ a + b + d + e + f + ───
╲╱ 25
最终方程式可以简单地转换为多项式：将分母相乘并对两边进行平方，例如：
In [24]: den = denom(equations[-1].lhs)
In [25]: neweq = Eq((equations[-1].lhs*den)**2, (equations[-1].rhs*den)**2)
In [26]: neweq
Out[26]:
2 2 2 ⎛ 329⎞
b ⋅f = a ⋅kp₁ ⋅⎜a + b + d + e + f + ───⎟
⎝ 25⎠
In [27]: equations[-1] = neweq
In [28]: for eq in equations:
...: pprint(eq)
...:
a + b = 13/5
2⋅a + b + d + 2⋅f = 7
d + e = 2
a⋅e = b⋅d⋅kp₂
2 2 2 ⎛ 329⎞
b ⋅f = a ⋅kp₁ ⋅⎜a + b + d + e + f + ───⎟
⎝ 25⎠
既然我们有了一个多项式系统，我们就可以计算它的 Groebner 基了。我们将 select a, b, d, e, f 作为 Groebner 基础的生成器，并留下 kp1 和 kp2 作为系数环的一部分：
In [29]: basis = list(groebner(equations, [a, b, d, e, f]))
现在这是 abdef 的简化方程组，具有依赖于 kp1 和 kp2 的复杂外观系数。我将以其中一个等式为例：
In [30]: basis[-1]
Out[30]:
⎛ 4 2 2 2 2 ⎞ ⎛
4 ⎛ 4 2 2 2 2 ⎞ 3 ⎜778⋅kp₁ ⋅kp₂ 399⋅kp₁ ⋅kp₂ 384⋅kp₁ 1⎟ 2 ⎜102481⋅
f ⋅⎝kp₁ ⋅kp₂ - kp₁ ⋅kp₂ - kp₁ + 1⎠ + f ⋅⎜───────────── - ───────────── - ──────── + ─⎟ + f ⋅⎜───────
⎝ 25 25 25 5⎠ ⎝ 6
4 2 2 2 2 2 ⎞ ⎛ 4 2 2 2
kp₁ ⋅kp₂ 15579⋅kp₁ ⋅kp₂ 13⋅kp₁ ⋅kp₂ 5148⋅kp₁ 1 ⎟ ⎜ 3799752⋅kp₁ ⋅kp₂ 8991⋅kp₁ ⋅kp₂
───────── + ─────────────── - ─────────── + ───────── + ───⎟ + f⋅⎜- ───────────────── - ──────────────
25 500 5 125 100⎠ ⎝ 3125 625
2 2⎞ 4 2
5772⋅kp₁ ⋅kp₂ 15984⋅kp₁ ⎟ 23853456⋅kp₁ ⋅kp₂
- ───────────── - ──────────⎟ + ──────────────────
125 625 ⎠ 15625
尽管这些方程式可能看起来很复杂，但它们具有特定的结构，这意味着一旦您为 kp1 和 kp2 代入值，最终方程式仅取决于 f，而所有其他方程式根据 f 线性给出符号。这意味着您可以替换您的值，然后使用 sympy 的 real_roots 或 numpys np.roots 来获得 f 的精确根或近似根，将它们代入您的其他方程并求解一个简单的系统剩余未知数的线性方程组。这些是在最终循环中需要完成的唯一步骤，我们可以使用 lambdify 使它们快速发生。那么最后的结果就是
from sympy import *
# First solve the system as much as possible in terms of symbolic kp1, kp2
kp1, kp2 = symbols('kp1, kp2')
a, b, d, e, f = symbols('a, b, d, e, f')
# Don't use any floats (yet)
equations = [
Eq(a+b, Rational('2.6')),
Eq(2*a + b + d + 2*f, 7),
Eq(d + e, 2),
Eq(a*e,kp2*b*d),
Eq( ((b * sqrt(f))/a)/sqrt((a+b+d+e+f+Rational('13.16'))), kp1)
]
# Convert to full polynomial system:
den = denom(equations[-1].lhs)
neweq = Eq((equations[-1].lhs*den)**2, (equations[-1].rhs*den)**2)
equations2 = equations[:-1] + [neweq]
# Compute the Groebner basis and split the linear equations for a, b, d, e
# from the polynomial equation for f
basis = list(groebner(equations2, [a, b, d, e, f]))
lineqs, eq_f = basis[:-1], basis[-1]
# Solve for a, b, d, e in terms of f, kp1 and kp2
[sol_abde] = linsolve(lineqs, [a, b, d, e])
# Use lambdify to be able to efficiently evaluate the solutions for a, b, d, e
sol_abde_lam = lambdify((f, kp1, kp2), sol_abde)
# Use lambdify to efficiently substitute kp1 and kp2 into the coefficients for eq_f
eq_f_lam = lambdify((kp1, kp2), Poly(eq_f, f).all_coeffs())
def solve_k(kp1val, kp2val):
eq_f_coeffs = eq_f_lam(kp1val, kp2val)
# Note that sympy's real_roots function is more accurate and does not need
# a threshold. It is however slower than np.roots:
#
# p_f = Poly(eq_f_coeffs, f)
# f_vals = real_roots(p_f) # exact (RootOf)
# f_vals = [r.n() for r in real_roots(p_f)] # approximate
#
f_vals = np.roots(eq_f_coeffs)
f_vals = f_vals.real[abs(f_vals.imag) < 1e-10] # arbitrary threshold
sols = []
for f_i in f_vals:
abde_vals = sol_abde_lam(f_i, kp1val, kp2val)
sol = {f: f_i}
sol.update(zip([a,b,d,e], abde_vals))
sols.append(sol)
return sols
现在您可以在大约 1 毫秒内为 kp1 和 kp2 的特定值生成解决方案：
In [40]: %time solve_k(2.1, 3.8)
CPU times: user 1.59 ms, sys: 0 ns, total: 1.59 ms
Wall time: 1.51 ms
Out[40]:
[{a: 49.77432869943, b: -47.174328699430006, d: -0.7687860254920669, e: 2.768786025492067, f: -22.30277
1336968966}, {a: 3.4616794447024692, b: -0.8616794447024684, d: 36.964491584342106, e: -34.964491584342
11, f: -18.01308551452229}, {a: -0.5231222050642267, b: 3.1231222050642273, d: -0.0922228459726158, e:
2.092222845972616, f: 2.507672525518422}, {a: 0.34146680496125464, b: 2.258533195038746, d: 0.076528664
56901455, e: 1.9234713354309858, f: 1.9910022652348665}]
注意事项：上述方法在大多数情况下都适用于像这样的多项式系统，尽管在某些情况下，如果不额外使用分裂变量，Groebner 基将无法完全分离。基本上只是引入一个新变量，比如 t 和一个随机方程，比如 t = a + 2*b + 3*d + 4*e + 5*f。然后将 t 作为 Groebner 基的最后一个符号，并计算 t 而不是 f 和 np.roots。同样在某些情况下，计算 Groebner 基础可能会非常慢，但请记住，它只需要计算 ever 一次（如果需要，您可以将结果保存在文件中）。
另一个警告是平方变换为多项式系统会引入不满足原始根方程的虚假解。代入原系统即可查看。在这种情况下，只有最后一个等式很重要：
In [41]: sols = solve_k(2.1, 3.8)
In [42]: equations[-1].subs(sols[0])
Out[42]: -2.10000000000001 = kp₁
In [43]: equations[-1].subs(sols[1])
Out[43]: -2.10000000000006 = kp₁
In [44]: equations[-1].subs(sols[2])
Out[44]: -2.09999999999995 = kp₁
In [45]: equations[-1].subs(sols[3])
Out[45]: 2.09999999999993 = kp₁
所以我们看到返回的 4 个解决方案中只有最后一个是正确的，因为它（大约）满足 kp1 = 2.1，这是传入的值。这意味着需要一点 post-processing获得您真正想要的解决方案。
我正在尝试使用以下 python 代码（当然使用 SymPy）来求解一组方程式：
def Solve(kp1, kp2):
a, b, d, e, f = S('a b d e f'.split())
equations = [
Eq(a+b, 2.6),
Eq(2*a + b + d + 2*f, 7),
Eq(d + e, 2),
Eq(a*e,kp2*b*d),
Eq( ((b * (f**0.5))/a)*((a+b+d+e+f+13.16)**-0.5), kp1)
]
return solve(equations)
代码成功求解了方程，但大约需要 35 秒。 Solve() 函数在另一个文件的迭代块（大约 2000 次迭代）中使用，所以速度对我来说真的很重要。
有什么方法可以加快求解器的速度吗？如果不能，您能否推荐另一种使用 python?
求解方程组的方法你只需要解一次方程。之后，您将得到以下形式的方程式：
- a = f_1(kp1, kp2)
- b = f_2(kp1, kp2)
- ...
因此您可以根据 kp1 和 kp2 简单地计算 a, ..., e。例如求解第一个、第三个和第四个方程得到：
- b: -a + 2.6
- e: 2.0*kp2*(5.0*a - 13.0)/(5.0*a*kp2 - 5.0*a - 13.0*kp2),
- d: 10.0*a/(-5.0*a*kp2 + 5.0*a + 13.0*kp2)
在我的电脑上求解所有五个方程太慢，但如果它给你一个表达式，你只需插入（替换）kp1 和 kp2 的值，你不必再次求解方程。对于替换，请查看 sympy documentation。
所以你的循环应该是这样的：
solutions = sympy.solve(eqs, exclude=[kp1, kp2])
for data_kp1, data_kp2 in data:
for key, eq in solutions:
solution = eq.subs([(kp1, data_kp1), (kp2, data_kp2)])
final_solutions.append("{key}={solution}".format(key=key, solution=solution))
求解 a、b、d、e 的第 4 个线性方程。这会生成两个解决方案。将这些中的每一个代入第 5 个方程（以这种形式：Eq(b**2*f, d**2*kp1**2*(a + b + 2*d + f + 13.16))) 给出两个方程（e51 和 e52）。这两个方程在 f 中是非线性的，如果你对它们使用 unrad，你将得到 2 第六f 中的阶多项式——这可能是您得不到解决方案的原因，因为通常只有四次方程是完全可解的。将这两个方程称为 e51u 和 e52u。
如果您对实根感兴趣，可以使用 real_roots 为这些多项式 提供根，然后 将所需的 kp1 和 kp2 值代入只留下 f 作为未知数。例如，
>>> ans = solve(equations[:-1],a,b,d,e)
>>> l=equations[-1] # modified as suggested
>>> l=l.lhs-l.rhs
>>> l
b**2*f - d**2*kp1**2*(a + b + 2*d + f + 13.16)
>>> e51=l.subs(ans[0]); e51u=unrad(e51)[0]
>>> e52=l.subs(ans[1]); e52u=unrad(e52)[0]
>>> import random
>>> for r1,r2 in [[random.random() for i in range(2)] for j in range(3)]:
... print [i.n(2) for i in real_roots(e51u.subs(dict(kp1=r1,kp2=r2)))]
... print [i.n(2) for i in real_roots(e52u.subs(dict(kp1=r1,kp2=r2)))]
... print '^_r1,r2=',r1,r2
...
[1.7, 2.9, 3.0, 8.2]
[1.7, 2.9, 3.0, 8.2]
^_r1,r2= 0.937748743197 0.134640776315
[1.3, 2.3, 4.7, 7.4]
[1.3, 2.3, 4.7, 7.4]
^_r1,r2= 0.490002815309 0.324553144174
[1.1, 2.1]
[1.1, 2.1]
^_r1,r2= 0.308803300429 0.595356213169
看来e51u和e52u都给出了相同的解，所以你可能只需要使用其中一个。你应该检查原始方程式中的答案，看看哪些是真实的解决方案：
>>> r1,r2
(0.30880330042869408, 0.59535621316941589)
>>> [e51.subs(dict(kp1=r1,kp2=r2,f=i)).n(2) for i in real_roots(e51u.subs(dict(kp1=r1,kp2=r2)))]
[1.0e-12, 13.]
所以在这里你看到只有第一个解（从上面看到是1.1）实际上是一个解； 2.1 是一个虚假的解决方案。
您的方程组可以根据 kp1 和 kp2 作为符号参数转换为多项式系统。在这种情况下，将其转换为多项式系统是有利的，以便它可以部分求解，为参数的特定值的数值求解做准备。
她的技术涉及计算 Groebner 基数，您不希望在方程式中有任何浮点数，因此让我们首先确保使用精确的有理数（例如，使用 sqrt 而不是 **0.5 和 Rational('2.6') 而不是 2.6):
In [20]: kp1, kp2 = symbols('kp1, kp2')
In [21]: a, b, d, e, f = symbols('a, b, d, e, f')
In [22]: equations = [
...: Eq(a+b, Rational('2.6')),
...: Eq(2*a + b + d + 2*f, 7),
...: Eq(d + e, 2),
...: Eq(a*e,kp2*b*d),
...: Eq( ((b * sqrt(f))/a)/sqrt((a+b+d+e+f+Rational('13.16'))), kp1)
...: ]
In [23]: for eq in equations:
...: pprint(eq)
...:
a + b = 13/5
2⋅a + b + d + 2⋅f = 7
d + e = 2
a⋅e = b⋅d⋅kp₂
b⋅√f
─────────────────────────────── = kp₁
_________________________
╱ 329
a⋅ ╱ a + b + d + e + f + ───
╲╱ 25
最终方程式可以简单地转换为多项式：将分母相乘并对两边进行平方，例如：
In [24]: den = denom(equations[-1].lhs)
In [25]: neweq = Eq((equations[-1].lhs*den)**2, (equations[-1].rhs*den)**2)
In [26]: neweq
Out[26]:
2 2 2 ⎛ 329⎞
b ⋅f = a ⋅kp₁ ⋅⎜a + b + d + e + f + ───⎟
⎝ 25⎠
In [27]: equations[-1] = neweq
In [28]: for eq in equations:
...: pprint(eq)
...:
a + b = 13/5
2⋅a + b + d + 2⋅f = 7
d + e = 2
a⋅e = b⋅d⋅kp₂
2 2 2 ⎛ 329⎞
b ⋅f = a ⋅kp₁ ⋅⎜a + b + d + e + f + ───⎟
⎝ 25⎠
既然我们有了一个多项式系统，我们就可以计算它的 Groebner 基了。我们将 select a, b, d, e, f 作为 Groebner 基础的生成器，并留下 kp1 和 kp2 作为系数环的一部分：
In [29]: basis = list(groebner(equations, [a, b, d, e, f]))
现在这是 abdef 的简化方程组，具有依赖于 kp1 和 kp2 的复杂外观系数。我将以其中一个等式为例：
In [30]: basis[-1]
Out[30]:
⎛ 4 2 2 2 2 ⎞ ⎛
4 ⎛ 4 2 2 2 2 ⎞ 3 ⎜778⋅kp₁ ⋅kp₂ 399⋅kp₁ ⋅kp₂ 384⋅kp₁ 1⎟ 2 ⎜102481⋅
f ⋅⎝kp₁ ⋅kp₂ - kp₁ ⋅kp₂ - kp₁ + 1⎠ + f ⋅⎜───────────── - ───────────── - ──────── + ─⎟ + f ⋅⎜───────
⎝ 25 25 25 5⎠ ⎝ 6
4 2 2 2 2 2 ⎞ ⎛ 4 2 2 2
kp₁ ⋅kp₂ 15579⋅kp₁ ⋅kp₂ 13⋅kp₁ ⋅kp₂ 5148⋅kp₁ 1 ⎟ ⎜ 3799752⋅kp₁ ⋅kp₂ 8991⋅kp₁ ⋅kp₂
───────── + ─────────────── - ─────────── + ───────── + ───⎟ + f⋅⎜- ───────────────── - ──────────────
25 500 5 125 100⎠ ⎝ 3125 625
2 2⎞ 4 2
5772⋅kp₁ ⋅kp₂ 15984⋅kp₁ ⎟ 23853456⋅kp₁ ⋅kp₂
- ───────────── - ──────────⎟ + ──────────────────
125 625 ⎠ 15625
尽管这些方程式可能看起来很复杂，但它们具有特定的结构，这意味着一旦您为 kp1 和 kp2 代入值，最终方程式仅取决于 f，而所有其他方程式根据 f 线性给出符号。这意味着您可以替换您的值，然后使用 sympy 的 real_roots 或 numpys np.roots 来获得 f 的精确根或近似根，将它们代入您的其他方程并求解一个简单的系统剩余未知数的线性方程组。这些是在最终循环中需要完成的唯一步骤，我们可以使用 lambdify 使它们快速发生。那么最后的结果就是
from sympy import *
# First solve the system as much as possible in terms of symbolic kp1, kp2
kp1, kp2 = symbols('kp1, kp2')
a, b, d, e, f = symbols('a, b, d, e, f')
# Don't use any floats (yet)
equations = [
Eq(a+b, Rational('2.6')),
Eq(2*a + b + d + 2*f, 7),
Eq(d + e, 2),
Eq(a*e,kp2*b*d),
Eq( ((b * sqrt(f))/a)/sqrt((a+b+d+e+f+Rational('13.16'))), kp1)
]
# Convert to full polynomial system:
den = denom(equations[-1].lhs)
neweq = Eq((equations[-1].lhs*den)**2, (equations[-1].rhs*den)**2)
equations2 = equations[:-1] + [neweq]
# Compute the Groebner basis and split the linear equations for a, b, d, e
# from the polynomial equation for f
basis = list(groebner(equations2, [a, b, d, e, f]))
lineqs, eq_f = basis[:-1], basis[-1]
# Solve for a, b, d, e in terms of f, kp1 and kp2
[sol_abde] = linsolve(lineqs, [a, b, d, e])
# Use lambdify to be able to efficiently evaluate the solutions for a, b, d, e
sol_abde_lam = lambdify((f, kp1, kp2), sol_abde)
# Use lambdify to efficiently substitute kp1 and kp2 into the coefficients for eq_f
eq_f_lam = lambdify((kp1, kp2), Poly(eq_f, f).all_coeffs())
def solve_k(kp1val, kp2val):
eq_f_coeffs = eq_f_lam(kp1val, kp2val)
# Note that sympy's real_roots function is more accurate and does not need
# a threshold. It is however slower than np.roots:
#
# p_f = Poly(eq_f_coeffs, f)
# f_vals = real_roots(p_f) # exact (RootOf)
# f_vals = [r.n() for r in real_roots(p_f)] # approximate
#
f_vals = np.roots(eq_f_coeffs)
f_vals = f_vals.real[abs(f_vals.imag) < 1e-10] # arbitrary threshold
sols = []
for f_i in f_vals:
abde_vals = sol_abde_lam(f_i, kp1val, kp2val)
sol = {f: f_i}
sol.update(zip([a,b,d,e], abde_vals))
sols.append(sol)
return sols
现在您可以在大约 1 毫秒内为 kp1 和 kp2 的特定值生成解决方案：
In [40]: %time solve_k(2.1, 3.8)
CPU times: user 1.59 ms, sys: 0 ns, total: 1.59 ms
Wall time: 1.51 ms
Out[40]:
[{a: 49.77432869943, b: -47.174328699430006, d: -0.7687860254920669, e: 2.768786025492067, f: -22.30277
1336968966}, {a: 3.4616794447024692, b: -0.8616794447024684, d: 36.964491584342106, e: -34.964491584342
11, f: -18.01308551452229}, {a: -0.5231222050642267, b: 3.1231222050642273, d: -0.0922228459726158, e:
2.092222845972616, f: 2.507672525518422}, {a: 0.34146680496125464, b: 2.258533195038746, d: 0.076528664
56901455, e: 1.9234713354309858, f: 1.9910022652348665}]
注意事项：上述方法在大多数情况下都适用于像这样的多项式系统，尽管在某些情况下，如果不额外使用分裂变量，Groebner 基将无法完全分离。基本上只是引入一个新变量，比如 t 和一个随机方程，比如 t = a + 2*b + 3*d + 4*e + 5*f。然后将 t 作为 Groebner 基的最后一个符号，并计算 t 而不是 f 和 np.roots。同样在某些情况下，计算 Groebner 基础可能会非常慢，但请记住，它只需要计算 ever 一次（如果需要，您可以将结果保存在文件中）。
另一个警告是平方变换为多项式系统会引入不满足原始根方程的虚假解。代入原系统即可查看。在这种情况下，只有最后一个等式很重要：
In [41]: sols = solve_k(2.1, 3.8)
In [42]: equations[-1].subs(sols[0])
Out[42]: -2.10000000000001 = kp₁
In [43]: equations[-1].subs(sols[1])
Out[43]: -2.10000000000006 = kp₁
In [44]: equations[-1].subs(sols[2])
Out[44]: -2.09999999999995 = kp₁
In [45]: equations[-1].subs(sols[3])
Out[45]: 2.09999999999993 = kp₁
所以我们看到返回的 4 个解决方案中只有最后一个是正确的，因为它（大约）满足 kp1 = 2.1，这是传入的值。这意味着需要一点 post-processing获得您真正想要的解决方案。