python - 两个三次表达式之间的解析交集

Question

我正在求解 2 条三次曲线的解析交集，其参数在以下代码中的两个单独函数中定义。

通过绘制曲线，很容易看出有一个交点:

放大版:

但是，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

Answer 1

问题

是您对解决方案的真实假设与 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)

因为你问:我使用复合体是一种品味问题。不需要，这取决于您的进一步目的。不过，这样做没有任何限制。出于一般性原因，我将此 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)]

而不是您可能期望的 float 。该解决方案是精确且完全准确的。但是，在我看来，它不适合进一步计算。这就是我在每个 sympy 结果上使用 evalf() 的原因。如果在此示例中对所有结果使用数值评估，则输出将变为:

[-1.41421356237310, 1.41421356237310]

为什么它不适合您可能会问的进一步计算？记住你的第一个代码。找到的第一个 root sympy 是

-6.32436145176552 + 0.e-19*I

呵呵，虚部为零，不错。但是，如果您打印 sym.im(x) == 0，则输出为 False。计算机和“精确”语句是敏感的组合。那里要小心。

解决方案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() 后仍然精确，但它们会因为太小而被切掉。这就是为什么我会一直建议使用最简单、最通用的解决方案。之后，查看您的多项式并了解所需的数值精度。