sympy 解决与 solveset 与 nsolve

Question

我正在尝试为 r 求解以下等式：

from sympy import pi, S, solve, solveset, nsolve, symbols

(n_go, P_l, T, gamma_w, P_g, r, R_mol) = symbols(
    'n_go, P_l, T, gamma_w, P_g, r, R_mol', real=True)

expr = -P_g + P_l - 3*R_mol*T*n_go/(4*r**3*pi) + 2*gamma_w/r
soln = solveset(expr, r, domain=S.Reals)
soln1 = solve(expr, r)

soln 是 Complement(Intersection(FiniteSet(...))) 的形式，我真的不知道该怎么办。 soln1 是 3 个表达式的列表，其中两个是复数。事实上，如果我用符号替换值并计算 soln1 的解，所有的都是复数：

vdict = {n_go: 1e-09, P_l: 101325, T: 300, gamma_w: 0.07168596252716256, P_g: 3534.48011713030, R_mol: 8.31451457896800}
for result in soln1:
    print(result.subs(vdict).n())

returns:

-9.17942953565355e-5 + 0.000158143657514283*I
-9.17942953565355e-5 - 0.000158143657514283*I
0.000182122477993494 + 1.23259516440783e-32*I

有趣的是，先代入值，然后使用 solveset() 或 solve() 给出真实结果：

solveset(expr.subs(vdict), r, domain=S.Reals).n()

{0.000182122477993494}

相反，nsolve 使用此等式会失败，除非起点包含解的前 7 位有效数字 (!):

nsolve(expr.subs(vdict), r,0.000182122 )

ValueError: Could not find root within given tolerance. (9562985778.9619347103 > 2.16840434497100886801e-19)

应该不难，剧情如下：

我的问题：

为什么 nsolve 在这里没用？
如何使用从 solveset 返回的解来计算任何数值解？
为什么先求解再代入值，solve得不到真正的解？

Answer 1

你的 expr 本质上是一个三次方程。

在解决之前或之后应用 subs 应该不会有任何实质性的改变。

`soln`

soln 的形式为 Complement(Intersection(FiniteSet(<3 cubic solutions>), Reals), FiniteSet(0)) 即实数域上不包括 0 的三次解。

下面应该给你一个简单的 FiniteSet 返回，但是 evalf 似乎没有很好地实现集合。

print(soln.subs(vdict).evalf())

希望很快能有所作为。

1

之所以nsolve没有用，是因为图形几乎是渐近垂直的。根据您的图表，梯度大约为 1.0e8。我不认为 nsolve 对于如此陡峭的图表有用。

绘制您替换的表达式我们得到：缩小我们得到：

这是一个非常疯狂的函数，我怀疑 nsolve 使用的 epsilon 在这种情况下非常有用。要解决此问题，您可以在替换时提供更接近 1 的更合理的数字。（考虑提供不同的测量单位。例如，而不是 meters/year 考虑 km/hour）

2

很难告诉您一般如何处理 solveset 的输出，因为每种类型的集合都需要以不同的方式处理。它在数学上也不合理，因为 soln.args[0].args[0].args[0] 应该给出第一个三次解，但它忘记了这必须是实数且非零。

您可以使用 args 或 preorder_traversal 或其他东西来导航树。阅读各种集合的文档也应该有所帮助。 solve 和 solveset 需要“交互”使用，因为有很多可能的输出，有很多方法可以理解它。

3

我相信 soln1 有 3 个解决方案，而不是您所说的 4 个。否则，您的循环将打印 4 行而不是 3 行。所有这些在技术上都很复杂（Python 中浮点数的性质）。但是，您提供的第三个解决方案的虚部非常小。要删除这些挑剔的东西，有一个名为 chop 的参数应该有所帮助：

for result in soln1:
    print(result.subs(vdict).n(chop=True))

其中一个结果是 0.000182122477993494，它看起来像您的根。

Answer 2

Maelstrom 的回答很好，但我想补充几点。

您替换的值都是浮点数，这些值的多项式是 ill-conditioned。这意味着您代入的表达式的形式会影响返回结果的准确性。这就是为什么将值代入 solve 的解决方案不一定给出与调用 solve 之前代入的值完全相同的原因之一。

此外，在您替换符号之前，求解不可能知道三个根中的哪一个是真实的。这就是为什么您从 solve(expr, r) 中得到三个解而从 solve(expr.subs(vdict), r) 中只得到一个解的原因。替换后为实数的第三个解与替换后 solve 返回的相同（忽略微小的虚部）：

In [7]: soln1[2].subs(vdict).n()                                                                                                                              
Out[7]: 0.000182122477993494 + 1.23259516440783e-32⋅ⅈ

In [8]: solve(expr.subs(vdict), r)                                                                                                                            
Out[8]: [0.000182122477993494]

因为多项式是ill-conditioned并且在根处有很大的梯度nsolve很难找到这个根。但是，如果给定足够窄的间隔，nsolve 可以找到根：

In [9]: nsolve(expr.subs(vdict), r, [0.0001821, 0.0001823])                                                                                                   
Out[9]: 0.000182122477993494

因为这本质上是一个多项式，所以最好的办法实际上是将其转换为多项式并使用 nroots。最快的方法是使用 as_numer_denom，尽管在这种情况下会在零处引入伪根：

In [26]: Poly(expr.subs(vdict).as_numer_denom()[0], r).nroots()                                                                                               
Out[26]: [0, 0.000182122477993494, -9.17942953565356e-5 - 0.000158143657514284⋅ⅈ, -9.17942953565356e-5 + 0.000158143657514284⋅ⅈ]

Answer 3

下面是对基本问题的回答：如何有效地计算上述方程的根？

根据@OscarBenjamin 的建议，我们可以通过使用 Poly 和 roots 而不是 nroots 来做得更好更快。下面，sympy 立即计算出 P_g 的 100 个不同值的方程根，同时保持其他一切不变：

from sympy import pi, Poly, roots, solve, solveset, nsolve, nroots, symbols
(n_go, P_l, T, gamma_w, P_g, r, R_mol) = symbols(
    'n_go, P_l, T, gamma_w, P_g, r, R_mol', real=True)
vdict = {pi:pi.n(), n_go:1e-09, P_l:101325, T:300, gamma_w:0.0717, R_mol: 8.31451457896800}
expr = -P_g + P_l - 3*R_mol*T*n_go/(4*r**3*pi) + 2*gamma_w/r
expr_poly = Poly(expr.as_numer_denom()[0], n_go, P_l, T, gamma_w, P_g, r, R_mol, domain='RR[pi]')
result = [roots(expr_poly.subs(vdict).subs(P_g, val)).keys() for val in range(4000,4100)]

剩下的就是检查解决方案是否满足我们的条件（积极的，真实的）。谢谢@OscarBenjamin！

PS: 我是否应该扩展上面的主题以包括 nroots 和 roots？

sympy 解决与 solveset 与 nsolve

sympy solve vs. solveset vs. nsolve

solver

sympy

`soln`

1

2

3