SymPy无法求解两平面相交方程组的问题咨询

问题重现 #

用户尝试用SymPy求解两个平面的交线方程组，两次调用solve()均返回空列表，代码如下：

第一次尝试代码：

import sympy as sm
x,y,z = sm.symbols('x y z ')
p1 = sm.Plane(sm.Point3D(0,0,0), (1,3,5))
p1Eq = p1.equation()
display(p1Eq)
p2 = sm.Plane(sm.Point3D(1,3,3), (-1,2,2))
p2Eq = p2.equation()
display(p2Eq)
sm.solve( (p1Eq, p2Eq), (x,y,z))

第二次尝试代码：

eq1 = sm.Eq(p1Eq,0)
eq2 = sm.Eq(p2Eq,0)
x,y,z = sm.symbols('x y z ')
sm.solve( (eq1, eq2), (x,y,z))

原因分析 #

两个平面相交的结果是一条直线，包含无穷多组连续解，而SymPy的solve()函数默认只求解方程组的有限离散解（孤立点）。由于直线无法用有限个点表示，solve()无法返回符合预期的结果，因此返回空列表。

解决方法 #

方法1：直接使用Plane类的intersection()方法 #

SymPy的Plane类提供了专门计算平面相交的方法，可直接得到交线对象：

import sympy as sm

x,y,z = sm.symbols('x y z ')
p1 = sm.Plane(sm.Point3D(0,0,0), (1,3,5))
p2 = sm.Plane(sm.Point3D(1,3,3), (-1,2,2))

# 获取交线（返回列表，取第一个元素即交线对象）
intersection_line = p1.intersection(p2)[0]
print(intersection_line)

# 输出交线的参数方程形式
print(intersection_line.arbitrary_point())

方法2：用参数表示无穷解 #

如果要通过solve()处理，可引入一个自由参数来表示无穷多解：

import sympy as sm

x,y,z,t = sm.symbols('x y z t')
p1_eq = sm.Plane(sm.Point3D(0,0,0), (1,3,5)).equation()
p2_eq = sm.Plane(sm.Point3D(1,3,3), (-1,2,2)).equation()

# 以z为自由变量，求解x、y的表达式
sol = sm.solve( (p1_eq, p2_eq), (x,y), dict=True )[0]
# 将z替换为参数t，得到参数化的解
param_sol = {x: sol[x].subs(z, t), y: sol[y].subs(z, t), z: t}
print(param_sol)

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