如何使用Python Sympy求解两条函数曲线的最近点对

核心原理 #

两条曲线上距离最近的点对连线，必然是两条曲线的公法线段：即连线与两点处的切线同时垂直。这个几何条件对任意次多项式曲线都成立，据此可以列方程组求解，不需要针对二次、三次多项式写特殊逻辑。
该条件和直接最小化两点距离平方的优化条件完全等价：对距离平方函数求偏导令其为0，得到的方程组和几何法推导的结果没有差异。

实现步骤 #

• 设第一条曲线表达式为y = f(x)，取曲线上任意点的横坐标为x1，对应点坐标为(x1, f(x1))；第二条曲线表达式为y = g(x)，取曲线上任意点的横坐标为x2，对应点坐标为(x2, g(x2))
• 计算两条曲线的导函数f'(x)、g'(x)，分别代入x1、x2得到两点处的切线斜率
• 列公法线约束方程，为了避免分母为0的问题，交叉相乘整理为整式方程：
  1. (g(x2)-f(x1)) * f'(x1) = -(x2 - x1)
  2. (g(x2)-f(x1)) * g'(x2) = -(x2 - x1)
• 求解上述二元方程组，筛选所有实数解，计算每组解对应的两点距离，排序后即可得到最近点对，存在多组等距解时会全部返回
• 注意：5次及以上多项式不存在通用根式解析解，此时SymPy无法返回闭式结果，自动切换数值求解即可

通用可运行代码 #

import sympy as sp
from sympy import N

def closest_poly_pair(f, g, x=sp.symbols('x')):
    # 定义两个独立参数
    x1, x2 = sp.symbols('x1 x2', real=True)
    # 两点坐标
    y1 = f.subs(x, x1)
    y2 = g.subs(x, x2)
    # 两点处切线斜率
    dy1 = sp.diff(f, x).subs(x, x1)
    dy2 = sp.diff(g, x).subs(x, x2)
    # 列公法线整式方程
    eq1 = sp.Eq((y2 - y1)*dy1, -(x2 - x1))
    eq2 = sp.Eq((y2 - y1)*dy2, -(x2 - x1))
    
    # 优先尝试解析求解
    sol_list = []
    try:
        sol_list = sp.solve([eq1, eq2], [x1, x2], dict=True)
    except Exception:
        pass
    
    # 无解析解时切换数值求解，扫描多初始点找实根
    if not sol_list:
        print("检测到高次多项式无解析解，自动切换数值求解")
        # 扫描区间可根据实际需求调整
        for init_x1 in range(-5, 6):
            for init_x2 in range(-5, 6):
                try:
                    sol = sp.nsolve([eq1, eq2], [x1, x2], [init_x1, init_x2], maxsteps=100)
                    # 去重保留实根
                    sol_dict = {x1: N(sol[0]), x2: N(sol[1])}
                    is_dup = False
                    for s in sol_list:
                        if abs(N(s[x1]-sol_dict[x1])) < 1e-6 and abs(N(s[x2]-sol_dict[x2])) <1e-6:
                            is_dup = True
                            break
                    if not is_dup:
                        sol_list.append(sol_dict)
                except Exception:
                    continue
    
    # 过滤复数解，计算距离排序
    res = []
    for sol in sol_list:
        px1, py1 = N(sol[x1]), N(y1.subs(sol))
        px2, py2 = N(sol[x2]), N(y2.subs(sol))
        if not (px1.is_real and py1.is_real and px2.is_real and py2.is_real):
            continue
        dist = sp.sqrt((px2-px1)**2 + (py2-py1)**2)
        res.append( (float(dist), (float(px1), float(py1)), (float(px2), float(py2))) )
    # 按距离从小到大排序
    res.sort(key=lambda x:x[0])
    return res

# 测试示例1：y=x² 和 y=-(x-3)²
x = sp.symbols('x')
f1 = x**2
g1 = -(x-3)**2
print("示例1结果：", closest_poly_pair(f1, g1, x))

# 测试示例2：y=x⁴+x²+x 和 y=x³+3x²
f2 = x**4 + x**2 + x
g2 = x**3 + 3*x**2
print("示例2结果：", closest_poly_pair(f2, g2, x))

上述代码对任意次多项式组合生效，返回结果按点对距离从小到大排序，支持返回多组有效解。

兜底方案 #

如果SymPy数值求解的效率或根覆盖度不满足需求，可以直接用通用数值优化方案：

• 将距离平方(x2-x1)**2 + (g(x2)-f(x1))**2作为二元目标函数
• 调用scipy.optimize中的全局优化接口（如差分进化算法）搜索全局最小值，再用局部优化接口精修结果，该方案对任意连续可导函数都生效，不局限于多项式场景。

内容的提问来源于stack exchange，提问作者Eisverygoodletter