sympy解决vs.solveset vs.nsolve

如何解决sympy解决vs.solveset vs.nsolve

我正在尝试求解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,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,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())

返回：

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

有趣的是，先替换值，然后使用solveset（）或solve（）得出真实结果：

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

{0.000182122477993494}

相反，除非开始点包含solution（！）的前7个有效数字，否则nsolve不能通过该方程式进行计算：

nsolve(expr.subs(vdict),0.000182122 )

ValueError：在给定的公差范围内找不到根。 （9562985778.9619347103> 2.16840434497100886801e-19）

不应该那么难，这是情节：

我的问题：

1. 为什么nsolve在这里这么没用？
2. 如何使用solveset返回的解来计算任何数值解？
3. 如果先求解然后替换值，为什么不能从solve获得真实的解决方案？

解决方法

您的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。 （请考虑提供不同的度量单位。例如，代替米/年，请考虑公里/小时）

2

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

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

3

我相信soln1有3个解决方案，而不是您所说的4个。否则，您的循环将打印4行而不是3行。从技术上讲，所有这些行都是复杂的（Python中的float就是这种情况）。但是，您提供的第三个解决方案的虚构分量很小。为了消除这些挑剔的东西，有一个名为chop的参数应该可以帮助您

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

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

,

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

您替换的值都是浮点数，并且这些值使多项式处于不良状态。这意味着您替换为表达式的形式会影响返回结果的准确性。这就是为什么在解决方案中将值代入解决方案并不一定提供与您在调用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]

因为多项式是病态的，并且在根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⋅ⅈ]

,

这里没有解决以下基本问题：如何有效地计算上述方程的根？

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

from sympy import pi,Poly,roots,solve,solveset,nsolve,nroots,symbols
(n_go,P_l,T,gamma_w,P_g,R_mol) = symbols(
    'n_go,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,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吗？