问题
你对解决方案的假设与Sympy对于数值不确定性的错误判断相结合,导致了问题。
如果您删除该赋值语句,将得到以下代码:
import sympy as sym
def H_I(P):
return (-941.254840173) + (0.014460465765)*P + (-9.41529726451e-05)*P**2 + (1.23485231253e-06)*P**3
def H_II(P):
return (-941.254313412) + (0.014234188877)*P + (-0.00013455013645)*P**2 + (2.58697027372e-06)*P**3
P = sym.symbols('P')
sol = sym.solve(H_I(P) - H_II(P) , P)
sol = [x.evalf() for x in sol]
print(sol)
输出结果为:
[-6.32436145176552 + 1.0842021724855e-19*I, 1.79012202335501 +
1.0842021724855e-19*I, 34.4111917095165 - 1.35525271560688e-20*I]
您可以通过 sym.re(x)
访问解的实部。
解决方法
如果您有特定的数字精度,我认为收集实际结果最简单的方法类似于以下代码:
def is_close(a,b,tol):
if abs(a-b)<tol: return True
else: return False
P = sym.symbols('P')
sol = sym.solve(H_I(P) - H_II(P) , P)
sol = [complex(x.evalf()) for x in sol]
real_solutions = []
for x in sol:
if is_close(x,x.real,10**(-10)): real_solutions.append(x.real)
print(real_solutions)
因为你问了: 我使用complex作为口味的问题。根据你进一步的目的,没有必要这样做。但是没有限制这样做。我将is_close()编写为一个通用函数。您可能希望在其他多项式中使用此代码或在不同的上下文中使用此函数,那么为什么不以聪明和可读的方式编写代码呢?然而,最初的目的是告诉我变量x及其实部re(x)是否在某个数值精度上相同,即虚部可以忽略不计。您还应该检查微小的实部,但我省略了它们。
编辑
小的虚部通常是求解过程中出现的复数减法残留物。sympy被视为精确处理他们。evalf()给出了精确解的数值计算或近似值。这与更好的精度无关。例如:
import sympy as sym
def square(P):
return P**2-2
P = sym.symbols('P')
sol2 = sym.solve(square(P),P)
print(sol2)
这段代码输出为:
[-sqrt(2), sqrt(2)]
而不是像你可能期望的那样输出浮点数。这个解决方案是准确无误的,但我认为它不适合进一步的计算。这就是为什么我在每个Sympy结果上使用evalf()的原因。如果你在这个例子中对所有的结果进行数字化评估,输出就会变成:
[-1.41421356237310, 1.41421356237310]
你可能会问为什么它不适合进一步的计算?还记得你的第一个代码吗?Sympy找到的第一个根是:
-6.32436145176552 + 0.e-19*I
哎呀,虚部是零,很好。然而,如果你打印 sym.im(x) == 0,输出就是错误的。电脑和语句"exact"是敏感的组合。要小心。
解决方案2
如果你只想摆脱微小的虚部而不是真正施加显式的数字精度,你可以在数字化评估中使用关键字.evalf(chop=True)。这实际上忽略了不必要的小数并且在你的原始代码中仅仅截去了虚部。考虑到你甚至满意于忽略任何虚部,正如你在你的回答中所说,这可能是你最好的解决方案。为了完整起见,下面是相应的代码。
P = sym.symbols('P')
sol = sym.solve(H_I(P) - H_II(P) , P)
sol = [x.evalf(chop=True) for x in sol]
但是需要注意的是,这并不与我的第一种方法非常不同,如果也针对实部进行“截断”,那么两种方法的差异在于:您不知道这会带来多少精度问题。如果您从未使用其他多项式,则可能没有问题。下面的代码应该可以解决问题:
def H_0(P):
return P**2 - 10**(-40)
P = sym.symbols('P')
sol = sym.solve(H_0(P) , P)
sol_full = [x.evalf() for x in sol]
sol_chop = [x.evalf(chop=True) for x in sol]
print(sol_full)
print(sol_chop)
即使您使用evalf()后,根数仍然完好无损且准确,但由于它们太小而被截断。这就是为什么我建议始终使用最简单、最通用的解决方案。之后,请查看您的多项式并了解所需的数值精度。
sym.solve(H_I(P) - H_II(P) , P)
返回了-6.32436145176552 + 0.e-19*I
和另外两个根。然而,[x.evalf() for x in sol]
返回了-6.32436145176552 + 1.0842021724855e-19*I
。evalf()
在0.e-19*I
上做了什么,使其转换为1.0842021724855e-19*I
? - DavidC.