两个三次表达式之间的分析交点

5
我正在解决两条三次曲线的分析交点问题,这两条曲线的参数在下面代码的两个不同函数中定义。
通过绘制曲线,很容易看出它们有一个交点:

enter image description here

放大版本:

enter image description here

然而,sym.solve 没有找到交点,即当要求 print 'sol_ H_I(P) - H_II(P) =', sol 时,没有返回结果:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
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

fig = plt.figure()

# Linspace for plotting the curves:
P_lin = np.linspace(-5.0, 12.5, 10000)

# Plotting the curves:
p1, = plt.plot(P_lin, H_I(P_lin), color='black' )
p2, = plt.plot(P_lin, H_II(P_lin), color='blue' )

# Labels:
fontP = FontProperties()
fontP.set_size('15')
plt.legend((p1, p2), ("Curve 1", "Curve 2"), prop=fontP)
plt.ticklabel_format(useOffset=False)

plt.savefig('2_curves.pdf', bbox_inches='tight')

plt.show()
plt.close()

# Solving the intersection:
P = sym.symbols('P', real=True)

sol = sym.solve(H_I(P) - H_II(P) , P)
print 'sol_ H_I(P) - H_II(P)  =', sol
3个回答

5

问题

你对解决方案的假设与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()后,根数仍然完好无损且准确,但由于它们太小而被截断。这就是为什么我建议始终使用最简单、最通用的解决方案。之后,请查看您的多项式并了解所需的数值精度。


1
非常感谢您的回答。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*Ievalf()0.e-19*I上做了什么,使其转换为1.0842021724855e-19*I - DavidC.
1
@DavidC。很高兴能帮助到您。我在我的答案中添加了更多信息,作为编辑 :) - Banana
非常感谢您更新的答案:)。我已经在思考这个问题,并考虑了一些替代方案。很抱歉让您等待了几天,但我需要时间来思考:) 请看下面是我的建议。我非常乐意讨论。如果我误解了您的某些方法,请见谅。 - DavidC.

0

让我们从第一个结果开始:

 sol = sym.solve(H_I(P) - H_II(P) , P)
 print 'sol_ H_I(P) - H_II(P)  =', sol

这段代码会输出以下内容:

[-6.32436145176552 + 0.e-19*I, 1.79012202335501 + 0.e-19*I, 34.4111917095165 - 0.e-20*I]

我完全同意你使用 evalf() 的策略,因为它可以提高精度:

 evalf_result = [x.evalf() for x in sol]
 print '[x.evalf() for x in sol] = ', evalf_result

这将产生:

[-6.32436145176552 + 1.0842021724855e-19*I, 1.79012202335501 + 1.0842021724855e-19*I, 34.4111917095165 - 1.35525271560688e-20*I]

您的解决方案意味着使用 Python 的内置 complex 函数,将上述结果转换为更好的元组式结果,其中 "I" 符号被 "j" 取代:

 complex_evalf_result = complex(x.evalf()) for x in sol
 print 'complex(x.evalf()) for x in sol = ', complex_evalf_result

这将得到以下结果:

[(-6.324361451765517+1.0842021724855044e-19j), (1.7901220233550066+1.0842021724855044e-19j), (34.41119170951654-1.3552527156068805e-20j)]

由于 complex_evalf_result 的类型为 complex,因此您现在可以使用 complex_evalf_result.real 策略来获取 complex_evalf_result 中每个 x 的实部。这一策略是您一直在采用的,我也同意。

一旦应用了 evalfcomplex 函数,您现在可以实现 is_close 函数方法,这是我觉得非常有趣的方法:

"如果实部和虚部之间的绝对差小于 10E-10,则舍弃虚部。"

一般来说,这通常适用于复数部分小于 10E-10 的情况。例如,对于

[(-6.324361451765517+1.0842021724855044e-19j), (1.7901220233550066+1.0842021724855044e-19j), (34.41119170951654-1.3552527156068805e-20j)]

发生了这样的事情:

abs(-6.324361451765517+1.0842021724855044e-19j - (-6.324361451765517)) = 1.0842021724855044e-19

这总是小于10E-10。

您的函数本质上是丢弃了复数部分(如果此函数还有其他应用程序并且我有点傻无法理解它,请原谅我)。

那么,为什么不使用这个更简单的解决方案呢?

 import numpy as np
 import sympy as sym
 from sympy.functions import re

 # Crude intersection:
 sol = sym.solve(H_I(P) - H_II(P) , P)
 print 'sol_ H_I(P) - H_II(P)  =', sol

 print """

 """
 # Use of evalf to obtain better precision:
 evalf_result = [x.evalf() for x in sol]
 print '[x.evalf() for x in sol] = ', evalf_result

 print """

 """
 # Now, let's grab the real part of the evalf_result:
 real_roots = []
 for x in evalf_result:
   each_real_root = re(x)
   real_roots.append(each_real_root)

 print 'real_roots = ', real_roots

这将直接打印:

real_roots = [-6.32436145176552, 1.79012202335501, 34.4111917095165]

通过遵循这种策略,发现:

1)不需要调用Python的complex内置策略。一旦evalf函数完成了工作,实部可以通过re(x)简单地提取。

2)不需要通过is_close函数传递我们的交点结果,只需丢弃复数部分即可。

请告诉我是否有什么我误解的地方,或者您不太同意的地方 - 我非常乐意讨论 :) 您所有的帮助都非常棒,非常感谢!


请注意,Stack Overflow 不是一个论坛。回答区域的帖子应该直接回答问题,而不是讨论其他答案。如果您想与用户 Banana 开展讨论,请找到另一个场所,比如 chat - user6655984
我认为你有一些误解,所以我又再次更新了我的答案(完全)。大多数问题的简单答案都是:通用性。我写的答案不仅适用于你提到的那2个多项式,而是基本上适用于所有多项式。我建议你再次阅读我的更新版本,但我同意@CrazyIvan的观点,如果还有任何未解决的问题,我很乐意与你交流。 - Banana

0

要找出一个多项式的根,请使用专用方法roots,而不是通用的solve

sol = sym.roots(H_I(P) - H_II(P) , P)

这将返回一个带有重复次数的字典形式的根,{-6.32436145176552: 1, 1.79012202335501: 1, 34.4111917095165: 1}

通常更方便的是获取根的列表(如果有多个根,则会重复):

sol = sym.roots(H_I(P) - H_II(P) , P, multiple=True) 

返回[-6.32436145176552, 1.79012202335501, 34.4111917095165]


如果这个方程不是多项式,我建议使用SciPy的求解器,例如fsolve,而不是SymPy。 SymPy不是找到充满浮点系数的方程的数值解的正确工具。它被设计用于符号数学,符号数学和浮点数不太搭配。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接