两个立方表达式之间的解析交集

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

问题

你对解决方案的假设是真实的,并且同意对数值不确定性的错误判断.如果您取消作业,最终会得到以下代码:

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)是否相同,直到某个数值精度,即虚部可忽略不计.您还应检查我忽略的可忽略不计的实部.

编辑

小的虚部通常是在求解过程中出现的复数上的减法残差.被视为准确,同情不会抹去它们.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)]

而不是您可能预期的浮点数.解决方案准确且完全准确.但是,在我看来,它不适合进一步计算.这就是我在每个同情结果上使用evalf()的原因.如果对此示例中的所有结果使用数值计算,则输出变为:

[-1.41421356237310,1.41421356237310]

为什么它不适合您可能会问的进一步计算？记住你的第一个代码.发现的第一个根发现是

-6.32436145176552 + 0.e-19*我

嗯,想象中的部分是零,很好.但是,如果打印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()之后仍然准确,它们会被切断,因为它们太小了.这就是为什么我会一直建议使用最简单,最通用的解决方案.之后,查看多项式并了解所需的数值精度.