两个三次表达式之间的解析交集
Analytic intersection between two cubic expressions
我正在求解 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
问题
你对解的假设是真实的，加上 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)]
而不是您可能期望的浮点数。该解决方案是精确且完全准确的。但是，在我看来，它不适合进一步计算。这就是我在每个 sympy 结果上使用 evalf() 的原因。如果在此示例中对所有结果使用数值评估，则输出变为：
[-1.41421356237310, 1.41421356237310]
您可能会问为什么它不适合进一步计算？记住你的第一个代码。找到的第一个 root sympy 是
-6.32436145176552 + 0.e-19*I
呵呵，虚部为零，不错。但是，如果您打印 sym.im(x) == 0，则输出为 False。计算机和语句 '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]
但是 请注意，这与我的第一种方法并没有太大的不同，如果有人也会为实数部分实现 "cut off"。然而不同之处在于：你对这强加的准确性一无所知。
如果您从不使用任何其他多项式，那可能没问题。以下代码应该可以说明问题：
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() 后仍然精确，但它们被切掉了，因为它们太小了。这就是为什么我会一直建议使用最简单、最通用的解决方案。之后，查看您的多项式并了解所需的数值精度。
让我们从第一个结果开始：
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]
您的解决方案意味着处理 complex python 的内置函数，它将上述结果转换为更好的元组结果，其中 "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)]
由于type(complex_evalf_result)returnscomplex，现在正好可以使用complex_evalf_result.real策略获取x中每个x的实部=20=]。这是你的策略，我同意。
应用 evalf 和 complex 函数后，您现在可以实现 is_close 函数方法，我发现它非常有趣：
"If the abs difference between the real part and the complex part are less than 10E-10, then discard the complex part."
这通常适用于复数部分小于 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 函数只是为了丢弃复杂的部分。
如果有什么我误解了的地方，或者你不太同意的地方，请告诉我 - 我很乐意讨论 :) 你的所有帮助都很棒，非常感谢！
要求多项式的根，请使用专用方法 roots 而不是通用方法 solve。
sol = sym.roots(H_I(P) - H_II(P) , P)
这个returns根作为具有多重性的字典，
{-6.32436145176552: 1, 1.79012202335501: 1, 34.4111917095165: 1}
获取根列表通常更方便（如果有多个根，将重复）：
sol = sym.roots(H_I(P) - H_II(P) , P, multiple=True)
returns [-6.32436145176552, 1.79012202335501, 34.4111917095165]
如果这个方程不是多项式，我建议使用 SciPy 的求解器，例如 fsolve 而不是 SymPy。 SymPy 不是寻找充满浮点系数的方程的数值解的正确工具。它旨在进行符号数学，而符号数学和浮点数不能很好地混合。
我正在求解 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
问题
你对解的假设是真实的，加上 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)]
而不是您可能期望的浮点数。该解决方案是精确且完全准确的。但是，在我看来，它不适合进一步计算。这就是我在每个 sympy 结果上使用 evalf() 的原因。如果在此示例中对所有结果使用数值评估，则输出变为：
[-1.41421356237310, 1.41421356237310]
您可能会问为什么它不适合进一步计算？记住你的第一个代码。找到的第一个 root sympy 是
-6.32436145176552 + 0.e-19*I
呵呵，虚部为零，不错。但是，如果您打印 sym.im(x) == 0，则输出为 False。计算机和语句 '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]
但是 请注意，这与我的第一种方法并没有太大的不同，如果有人也会为实数部分实现 "cut off"。然而不同之处在于：你对这强加的准确性一无所知。 如果您从不使用任何其他多项式，那可能没问题。以下代码应该可以说明问题：
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() 后仍然精确，但它们被切掉了，因为它们太小了。这就是为什么我会一直建议使用最简单、最通用的解决方案。之后，查看您的多项式并了解所需的数值精度。
让我们从第一个结果开始：
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]
您的解决方案意味着处理 complex python 的内置函数，它将上述结果转换为更好的元组结果，其中 "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)]
由于type(complex_evalf_result)returnscomplex，现在正好可以使用complex_evalf_result.real策略获取x中每个x的实部=20=]。这是你的策略，我同意。
应用 evalf 和 complex 函数后，您现在可以实现 is_close 函数方法，我发现它非常有趣：
"If the abs difference between the real part and the complex part are less than 10E-10, then discard the complex part."
这通常适用于复数部分小于 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 函数只是为了丢弃复杂的部分。
如果有什么我误解了的地方，或者你不太同意的地方，请告诉我 - 我很乐意讨论 :) 你的所有帮助都很棒，非常感谢！
要求多项式的根，请使用专用方法 roots 而不是通用方法 solve。
sol = sym.roots(H_I(P) - H_II(P) , P)
这个returns根作为具有多重性的字典，
{-6.32436145176552: 1, 1.79012202335501: 1, 34.4111917095165: 1}
获取根列表通常更方便（如果有多个根，将重复）：
sol = sym.roots(H_I(P) - H_II(P) , P, multiple=True)
returns [-6.32436145176552, 1.79012202335501, 34.4111917095165]
如果这个方程不是多项式，我建议使用 SciPy 的求解器，例如 fsolve 而不是 SymPy。 SymPy 不是寻找充满浮点系数的方程的数值解的正确工具。它旨在进行符号数学，而符号数学和浮点数不能很好地混合。