如何用纯数学二次方程求直线与椭圆的交点（Python实现）

需求说明 #

需要用Python计算椭圆与直线的交点，仅通过纯数学二次方程推导实现，已知条件包括直线的两个端点、椭圆的中心坐标及长/短半轴长度。要求不得使用numpy、scipy等第三方库，也不能采用三角学方法。此前参考的二次方程推导方案计算出的交点存在错误，既不在椭圆上也不在直线上，现需修正实现。

基础公式 #

椭圆标准方程：

((x - center_x)² / rad_x²) + ((y - center_y)² / rad_y²) = 1

直线斜截式方程（非垂直x轴时）：

y = m*x + c

原推导的二次方程系数（对应椭圆方程中a=rad_x，b=rad_y，椭圆中心(h,k)=(center_x, center_y)）：

A = a²m² + b²
B = 2a²m(c - k) - 2b²h
C = b²h² + a²(c - k)² - a²b²

原代码问题分析 #

原代码存在两处关键错误：

1. 二次项系数C计算错误：将(c - center_y)²误写为(c - center_y)*2，导致整个方程偏离正确解。
2. 垂直x轴的直线处理错误：当直线垂直x轴时（x1=x2），直接将斜率m设为0，这会导致直线方程完全错误，应单独处理这种情况。

修正后的实现代码 #

from math import sqrt

def get_intersection_points(line_point_a, line_point_b, center_x, center_y, rad_x, rad_y):
    x1, y1 = line_point_a
    x2, y2 = line_point_b
    intersections = []
    
    # 处理直线垂直于x轴的特殊情况
    if abs(x1 - x2) < 1e-9:
        x = x1
        # 代入椭圆方程求解y
        term = 1 - ((x - center_x)**2) / (rad_x**2)
        if term < -1e-9:
            return intersections  # 无交点
        term = max(term, 0)  # 避免浮点误差导致的负数开方
        y_offset = rad_y * sqrt(term)
        y3 = center_y + y_offset
        y4 = center_y - y_offset
        # 检查点是否在直线线段上（可选，若需要仅线段交点则保留）
        min_y = min(y1, y2)
        max_y = max(y1, y2)
        if min_y - 1e-9 <= y3 <= max_y + 1e-9:
            intersections.append([x, y3])
        if abs(y3 - y4) > 1e-9 and (min_y - 1e-9 <= y4 <= max_y + 1e-9):
            intersections.append([x, y4])
        return intersections
    
    # 计算直线的斜率m和截距c
    m = (y2 - y1) / (x2 - x1)
    c = y1 - m * x1
    
    # 计算二次方程的系数A、B、C
    a_sq = rad_x ** 2
    b_sq = rad_y ** 2
    k = center_y
    h = center_x
    A = a_sq * (m ** 2) + b_sq
    B = 2 * a_sq * m * (c - k) - 2 * b_sq * h
    C = b_sq * (h ** 2) + a_sq * ((c - k) ** 2) - a_sq * b_sq
    
    # 计算判别式，判断交点数量
    discriminant = B ** 2 - 4 * A * C
    if discriminant < -1e-9:
        return intersections  # 无实交点
    discriminant = max(discriminant, 0)  # 处理浮点误差
    sqrt_d = sqrt(discriminant)
    
    # 求解x
    x3 = (-B - sqrt_d) / (2 * A)
    x4 = (-B + sqrt_d) / (2 * A)
    
    # 求解对应的y
    y3 = m * x3 + c
    y4 = m * x4 + c
    
    # 可选：检查点是否在直线线段范围内（若只需要线段与椭圆交点则启用）
    min_x = min(x1, x2)
    max_x = max(x1, x2)
    min_y = min(y1, y2)
    max_y = max(y1, y2)
    
    def is_on_segment(px, py):
        return (min_x - 1e-9 <= px <= max_x + 1e-9) and (min_y - 1e-9 <= py <= max_y + 1e-9)
    
    if is_on_segment(x3, y3):
        intersections.append([round(x3, 6), round(y3, 6)])
    if abs(x3 - x4) > 1e-9 and is_on_segment(x4, y4):
        intersections.append([round(x4, 6), round(y4, 6)])
    
    return intersections

说明 #

1. 浮点误差处理：加入了1e-9的容差，避免因浮点计算精度问题导致的错误判断。
2. 线段交点可选：代码中包含了检查交点是否在直线线段上的逻辑，若需要求解直线（无限延伸）与椭圆的交点，可移除该部分判断。
3. 特殊情况覆盖：单独处理了垂直x轴的直线，避免斜率不存在导致的计算错误。

内容的提问来源于stack exchange，提问作者Lion In A Box