寻找复杂解的Sympy解，即使real=True?

我的意思是通过sympy.solve来解一个符号方程。这花了很长时间，所以我尝试了一些变化，直到我找到了一个可能的原因。它寻找复数解，而我只对实数解感兴趣。

下面的代码…

import sympy as sp
import numpy as np
x, y, z = sp.symbols('x, y, z', real=True, positive=True)
# Option X1: works
expr3 = 0.6*x**(9/8)*y**(7/5)*z**2.2
expr4 = 0.9*x**(1/2)*y**(4/5)*z**1.2
s3 = sp.solve(expr3 - expr4, x)
print('Exponents of x: 9/8, 1/2. Solution =', s3)
# Option X3: works
expr3 = 0.6*x**(11/9)*y**(7/5)
expr4 = 0.9*x**(1/2)*y**(4/5)
s3 = sp.solve(expr3 - expr4, x, check=True, minimal=True, rational=True, force=True)
print('Exponents of x: 11/9, 1/2. Solution =', s3)
# Option X2: takes "forever"
expr3 = 0.6*x**(11/9)*y**(7/5)
expr4 = 0.9*x**(1/2)*y**(4/5)
s3 = sp.solve(expr3 - expr4, x)
print('Exponents of x: 11/9, 1/2. Solution =', s3)

…产生以下输出

Exponents of x: 9/8, 1/2. Solution = [1.91313675093869/(y**(24/25)*z**(8/5)), 0.351080874270527*(-y**(-0.12)*z**(-0.2) - 0.726542528005361*I*y**(-0.12)*z**(-0.2))**8, 0.351080874270527*(-y**(-0.12)*z**(-0.2) + 0.726542528005361*I*y**(-0.12)*z**(-0.2))**8, 1.28055023108529*(0.324919696232906*y**(-0.12)*z**(-0.2) - I*y**(-0.12)*z**(-0.2))**8, 1.28055023108529*(0.324919696232906*y**(-0.12)*z**(-0.2) + I*y**(-0.12)*z**(-0.2))**8]
Exponents of x: 11/9, 1/2. Solution = [3*2**(8/13)*3**(5/13)/(4*y**(54/65)), (-2**(12/13)*3**(1/13)*cos(pi/13)/(2*y**(3/65)) - 2**(12/13)*3**(1/13)*I*sin(pi/13)/(2*y**(3/65)))**18, (-2**(12/13)*3**(1/13)*cos(pi/13)/(2*y**(3/65)) + 2**(12/13)*3**(1/13)*I*sin(pi/13)/(2*y**(3/65)))**18, (2**(12/13)*3**(1/13)*cos(2*pi/13)/(2*y**(3/65)) - 2**(12/13)*3**(1/13)*I*sin(2*pi/13)/(2*y**(3/65)))**18, (2**(12/13)*3**(1/13)*cos(2*pi/13)/(2*y**(3/65)) + 2**(12/13)*3**(1/13)*I*sin(2*pi/13)/(2*y**(3/65)))**18, (-2**(12/13)*3**(1/13)*cos(3*pi/13)/(2*y**(3/65)) - 2**(12/13)*3**(1/13)*I*sin(3*pi/13)/(2*y**(3/65)))**18, (-2**(12/13)*3**(1/13)*cos(3*pi/13)/(2*y**(3/65)) + 2**(12/13)*3**(1/13)*I*sin(3*pi/13)/(2*y**(3/65)))**18, (2**(12/13)*3**(1/13)*cos(4*pi/13)/(2*y**(3/65)) - 2**(12/13)*3**(1/13)*I*sin(4*pi/13)/(2*y**(3/65)))**18, (2**(12/13)*3**(1/13)*cos(4*pi/13)/(2*y**(3/65)) + 2**(12/13)*3**(1/13)*I*sin(4*pi/13)/(2*y**(3/65)))**18, (-2**(12/13)*3**(1/13)*cos(5*pi/13)/(2*y**(3/65)) - 2**(12/13)*3**(1/13)*I*sin(5*pi/13)/(2*y**(3/65)))**18, (-2**(12/13)*3**(1/13)*cos(5*pi/13)/(2*y**(3/65)) + 2**(12/13)*3**(1/13)*I*sin(5*pi/13)/(2*y**(3/65)))**18, (2**(12/13)*3**(1/13)*cos(6*pi/13)/(2*y**(3/65)) - 2**(12/13)*3**(1/13)*I*sin(6*pi/13)/(2*y**(3/65)))**18, (2**(12/13)*3**(1/13)*cos(6*pi/13)/(2*y**(3/65)) + 2**(12/13)*3**(1/13)*I*sin(6*pi/13)/(2*y**(3/65)))**18]

来自选项X1和X3，然后它被困在选项X2。我必须手动中断计算。
X3是对solve的选项进行修修补补以获得X2的解的许多尝试之一。从它的结果，我可以理解为什么X2需要这么长时间。

解决方案列表中的第一个元素是我正在寻找的，可能稍后会简化。我的问题如下:

Q1:如果我设置real=True, positive=True，为什么solve寻找并返回复解?
Q2:即使不理解Q1，有没有办法告诉solve不要寻找复杂的解决方案?(注意:我也希望positive=True为x，而不是x=0)
Q3:我可以使用solve的选项，在合理的时间内进行计算，即使我得到复杂的解决方案。
但是我怎样才能自动选择真正的解呢?
我会选择元素0，但我不确定它总是正确的。

与:

忽略sympy中的虚根
使用solve

默认情况下，solve将表达式中的浮点数转换为有理数，然后用于求解某些多项式方程。当计算需要很长时间时，我们可以使用关键字参数。对于第三种情况:

expr3 = 0.6*x**(11/9)*y**(7/5)
expr4 = 0.9*x**(1/2)*y**(4/5)
s3 = sp.solve(expr3 - expr4, x, check=False)
print('Exponents of x: 11/9, 1/2. Solution =', s3)

计算是立即的。另一个选项是要求solve避免将浮点数转换为有理数:

s3 = sp.solve(expr3 - expr4, x, rational=False)
print('Exponents of x: 11/9, 1/2. Solution =', s3)

同样，计算非常快。

现在到另一个问题。解决这个问题的一种方法可能是告诉solve只寻找特定的区域。下面的命令需要几秒钟来计算:

sp.solve([expr3 - expr4, x >= 0, y >= 0], x)
# [(1.75314834631487/y**(54/65),)]

但是，请注意省略了平凡解x=0。

让我们试一试你的第一个案子:

expr3 = 0.6*x**(9/8)*y**(7/5)*z**2.2
expr4 = 0.9*x**(1/2)*y**(4/5)*z**1.2
sp.solve([expr3 - expr4, x >= 0, y >= 0, z >= 0], x)
# [(1.91313675093869/(y**(24/25)*z**(8/5)),)]

如果使用理数，最后一个跑得更快:

In [45]: # Option X1: works
...: expr3 = Rational(3,5)*x**(Rational(9,8))*y**(Rational(7,5))
...: expr4 = Rational(9,10)*x**(Rational(1,2))*y**(Rational(4,5))
...: s3 = sp.solve(expr3 - expr4, x)
In [46]: print(s3)
[0, 3*2**(2/5)*3**(3/5)/(4*y**(24/25)), (-2**(4/5)*3**(1/5)/(8*y**(3/25)) + 2**(4/5)*3**(1/5)*sqrt(5)/(8*y**(3/25)) - 2**(4/5)*3**(1/5)*I*sqrt(sqrt(5)/8 + 5/8)/(2*y**(3/25)))**8, (-2**(4/5)*3**(1/5)/(8*y**(3/25)) + 2**(4/5)*3**(1/5)*sqrt(5)/(8*y**(3/25)) + 2**(4/5)*3**(1/5)*I*sqrt(sqrt(5)/8 + 5/8)/(2*y**(3/25)))**8, (-2**(4/5)*3**(1/5)*sqrt(5)/(8*y**(3/25)) - 2**(4/5)*3**(1/5)/(8*y**(3/25)) - 2**(4/5)*3**(1/5)*I*sqrt(5/8 - sqrt(5)/8)/(2*y**(3/25)))**8, (-2**(4/5)*3**(1/5)*sqrt(5)/(8*y**(3/25)) - 2**(4/5)*3**(1/5)/(8*y**(3/25)) + 2**(4/5)*3**(1/5)*I*sqrt(5/8 - sqrt(5)/8)/(2*y**(3/25)))**8]

我不知道Q1和Q2的答案。但是我设法获得了代码，这给了我一个真正的正解，一个人可以通过手动计算expr3和expr4的形式与OP相同(expr3 = number * product(xi**bi, i, 1, n),n变量和指数bi通常是有理的)。到目前为止，代码对所有用例都有效，包括那些没有通过发布的任何内容正确解决的用例。而且效果很快。

# Option X6: works
expr3 = 0.6*x**(4/9)*y**(1/3)
expr4 = 0.9*x**(2/5)*y**(2/5)
s3_allsoln = sp.solve([expr3 - expr4, x>0, y>0], x, check=False, minimal=True, rational=True, force=True, dict=True)
s3_real = [s3[x] for s3 in s3_allsoln if s3[x].is_real]
s3_realpos = [s3 for s3 in s3_real if s3.is_positive]
if(len(s3_realpos) != 1):
print('Warning, found', len(s3_realpos), 'real positive solutions')
s3 = s3_realpos[0]
print('Exponents of x: 4/9, 2/5. Solution =', s3)

输出:

Exponents of x: 4/9, 2/5. Solution = 31381059609*sqrt(6)*y**(3/2)/8388608

我还需要把数字因子化简成一个十进制数。但这可能要容易得多。

相关内容

最新更新

热门标签：