求两相交圆与已知点间的公切圆弧绘制方法

问题定义 #

已知：

• 圆C₁：圆心坐标$(x_1,y_1)$，半径$r_1$
• 圆C₂：圆心坐标$(x_2,y_2)$，半径$r_2$
• 定点P：坐标$(x_p,y_p)$
  目标：构造一段圆弧，满足：

1. 与C₁、C₂均相切（以图中外切为例，可扩展到内切/混合相切）
2. 经过点P
3. 圆弧圆心为O$(x,y)$，半径为R

数学推导 #

设目标圆弧的圆心为O$(x,y)$，半径为R，根据条件列出方程：

1. 与C₁外切：$\sqrt{(x-x_1)^2 + (y-y_1)^2} = R + r_1$
2. 与C₂外切：$\sqrt{(x-x_2)^2 + (y-y_2)^2} = R + r_2$
3. 经过P点：$\sqrt{(x-x_p)^2 + (y-y_p)^2} = R$

将三个方程两边平方消去根号：
$$
\begin{cases}
(x-x_1)^2 + (y-y_1)^2 = (R + r_1)^2 \quad (1) \
(x-x_2)^2 + (y-y_2)^2 = (R + r_2)^2 \quad (2) \
(x-x_p)^2 + (y-y_p)^2 = R^2 \quad (3)
\end{cases}
$$

消元简化 #

用方程(1)-(3)、(2)-(3)消去$R^2$，得到两个线性方程：
从(1)-(3)展开整理：
$$2(x_p - x_1)x + 2(y_p - y_1)y + 2Rr_1 = (x_1^2 + y_1^2) - (x_p^2 + y_p^2) - r_1^2 \quad (1a)$$
从(2)-(3)展开整理：
$$2(x_p - x_2)x + 2(y_p - y_2)y + 2Rr_2 = (x_2^2 + y_2^2) - (x_p^2 + y_p^2) - r_2^2 \quad (2a)$$

从(1a)解出R的表达式：
$$R = \frac{(x_1^2 + y_1^2 - x_p^2 - y_p^2 - r_1^2) - 2(x_p - x_1)x - 2(y_p - y_1)y}{2r_1} \quad (R1)$$

将(R1)代入(2a)，得到一个关于x、y的线性方程；再结合方程(3)，就得到一个二元非线性方程组，可通过数值方法求解。

编程实现思路 #

1. 数值解法（牛顿迭代法） #

因为方程组是非线性的，适合用牛顿迭代法求解，步骤如下：

• 初始值选取：用C₁与C₂的中点M$(\frac{x_1+x_2}{2}, \frac{y_1+y_2}{2})$作为初始圆心，计算初始R为$\sqrt{(M_x - x_p)^2 + (M_y - y_p)^2}$
• 迭代优化：
  1. 定义残差函数：
     $$
     \begin{cases}
     f_1(x,y,R) = (x-x_1)^2 + (y-y_1)^2 - (R + r_1)^2 \
     f_2(x,y,R) = (x-x_2)^2 + (y-y_2)^2 - (R + r_2)^2 \
     f_3(x,y,R) = (x-x_p)^2 + (y-y_p)^2 - R^2
     \end{cases}
     $$
  2. 计算雅可比矩阵J，包含每个残差函数对x、y、R的偏导数
  3. 解线性方程组$J \cdot \Delta = -f$，得到修正量$\Delta(x,y,R)$
  4. 更新$x = x + \Delta x$，$y = y + \Delta y$，$R = R + \Delta R$
  5. 重复步骤2-4，直到残差的绝对值小于设定的精度阈值（如$1e-6$）

2. 解的有效性验证 #

迭代结束后，验证以下条件是否满足（允许微小误差）：

• $| \sqrt{(x-x_1)^2 + (y-y_1)^2} - (R + r_1) | < \epsilon$
• $| \sqrt{(x-x_2)^2 + (y-y_2)^2} - (R + r_2) | < \epsilon$
• $| \sqrt{(x-x_p)^2 + (y-y_p)^2} - R | < \epsilon$

如果外切情况无解，尝试内切（将方程中的$R + r_1$替换为$|R - r_1|$，$R + r_2$替换为$|R - r_2|$）或混合相切的组合（共4种相切组合，逐一尝试）。

3. 代码示例伪代码 #

import numpy as np

def solve_tangent_arc(C1, C2, P, epsilon=1e-6, max_iter=100):
    x1, y1, r1 = C1
    x2, y2, r2 = C2
    px, py = P
    
    # 初始值：中点M作为圆心
    x = (x1 + x2) / 2.0
    y = (y1 + y2) / 2.0
    R = np.sqrt((x - px)**2 + (y - py)**2)
    
    for _ in range(max_iter):
        # 计算残差
        f1 = (x - x1)**2 + (y - y1)**2 - (R + r1)**2
        f2 = (x - x2)**2 + (y - y2)**2 - (R + r2)**2
        f3 = (x - px)**2 + (y - py)**2 - R**2
        
        # 检查残差是否足够小
        if abs(f1) < epsilon and abs(f2) < epsilon and abs(f3) < epsilon:
            break
        
        # 计算雅可比矩阵
        J = np.array([
            [2*(x - x1), 2*(y - y1), -2*(R + r1)],
            [2*(x - x2), 2*(y - y2), -2*(R + r2)],
            [2*(x - px), 2*(y - py), -2*R]
        ])
        
        # 解线性方程组 J * delta = -f
        f = np.array([-f1, -f2, -f3])
        delta = np.linalg.solve(J, f)
        
        # 更新变量
        x += delta[0]
        y += delta[1]
        R += delta[2]
    
    return (x, y, R)

# 使用示例
C1 = (0, 0, 2)  # (x1,y1,r1)
C2 = (5, 0, 3)  # (x2,y2,r2)
P = (2, 3)      # (px,py)
center, radius = solve_tangent_arc(C1, C2, P)[:2], solve_tangent_arc(C1, C2, P)[2]

注意事项 #

• 如果迭代不收敛，可能是初始值选取得不好，可尝试其他初始点（比如P点到C1、C2的垂线交点）
• 需处理无解的情况：当P点位置无法构造满足条件的圆弧时，需返回错误提示
• 相切组合有4种（外切+外切、外切+内切、内切+外切、内切+内切），需逐一尝试找到有效解

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