为什么 SymPy 会计算出错误的平面交点？

Question

我有一个奇怪的问题，SymPy 中的平面相交适用于简单的示例，但对于具有更复杂坐标的示例却失败了。我 post 一个有效的简单示例，一个失败的示例。如 Povray 图像所示，我有三个平面穿过 a polyhedron and are perpendicular to the line through the respective vertex and the center. I would like to calculate the point where these planes intersect, but SymPy gives wrong results for the lines in which pairs of planes intersect. In the images the correct intersections can be seen as short lines (created with CSG 交点的顶点）。平行于它们的长线是SymPy计算出来的。

我是不是做错了什么，或者这是 SymPy 中的错误？

更多图片在这里：http://paste.watchduck.net/1712/sympy_planes/
有谁知道如何在页面上放置许多图像，而不被 post 排除在问题之外？ ("Your post appears to contain code that is not properly formatted as code.")

有效

代码：

from sympy import Point3D, Plane


pointR = Point3D(1/2, 0, 1/2)
pointG = Point3D(1, 0, 0)

planeR = Plane(pointR, pointR)
planeG = Plane(pointG, pointG)

print('\n######## Intersection of the planes:')
lineRG = planeR.intersection(planeG)[0]  # yellow
print(lineRG)

print('\n######## Intersection of plane and contained line returns the line:')
lineRG_again = planeR.intersection(lineRG)[0]
print(lineRG_again.equals(lineRG))

输出：

######## Intersection of the planes:
Line3D(Point3D(1, 0, 0), Point3D(1, 1/2, 0))

######## Intersection of plane and contained line returns the line:
True

失败

代码：

from sympy import sqrt, Point3D, Plane

pointR = Point3D(-1, 1 + sqrt(2), -2*sqrt(2) - 1)
pointG = Point3D(-sqrt(2) - 1, 1, -2*sqrt(2) - 1)
pointB = Point3D(-1, -sqrt(2) - 1, -2*sqrt(2) - 1)

planeR = Plane(pointR, pointR)
planeG = Plane(pointG, pointG)
planeB = Plane(pointB, pointB)

print('\n######## Intersections of the planes:')

lineRG = planeR.intersection(planeG)[0]  # yellow
lineRB = planeR.intersection(planeB)[0]  # magenta
lineGB = planeG.intersection(planeB)[0]  # cyan

print(lineRG)
print(lineRB)
print(lineGB)

print('\n######## Lines RG (yellow) and GB (cyan) intersect:')
print(lineRG.intersection(lineGB))
print('\n######## But line RB (magenta) is skew to both of them:')
print(lineRB.intersection(lineRG))
print(lineRB.intersection(lineGB))

print('\n######## Intersection of plane and contained line fails:')
lineRG_again = planeR.intersection(lineRG)

输出：

######## Intersections of the planes:
Line3D(Point3D(-1, 1, 0), Point3D((1 + sqrt(2))*(-2*sqrt(2) - 1) + 2*sqrt(2), -2*sqrt(2) + (-2*sqrt(2) - 1)*(-sqrt(2) - 1), -1 - (1 + sqrt(2))*(-sqrt(2) - 1)))
Line3D(Point3D(-1, 0, 0), Point3D((1 + sqrt(2))*(-2*sqrt(2) - 1) - (-2*sqrt(2) - 1)*(-sqrt(2) - 1) - 1, 0, 2 + 2*sqrt(2)))
Line3D(Point3D(-1, 1, 0), Point3D(-(-2*sqrt(2) - 1)*(-sqrt(2) - 1) - 2*sqrt(2) - 2, -(-2*sqrt(2) - 1)*(-sqrt(2) - 1) + 2 + 2*sqrt(2), 1 + (-sqrt(2) - 1)**2))

######## Lines RG (yellow) and GB (cyan) intersect:
[Point3D(-1, 1, 0)]

######## But line RB (magenta) is skew to both of them:
[]
[]

######## Intersection of plane and contained line fails:
Traceback (most recent call last):
  File "planes2.py", line 47, in <module>
    lineRG_again = planeR.intersection(lineRG)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/geometry/plane.py", line 390, in intersection
    p = a.subs(t, c[0])
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/core/basic.py", line 916, in subs
    rv = rv._subs(old, new, **kwargs)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/core/cache.py", line 93, in wrapper
    retval = cfunc(*args, **kwargs)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/core/basic.py", line 1030, in _subs
    rv = fallback(self, old, new)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/core/basic.py", line 1007, in fallback
    rv = self.func(*args)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/geometry/point.py", line 1104, in __new__
    args = Point(*args, **kwargs)
  File "/home/tilman/Code/024_polyhedron/env/lib/python3.5/site-packages/sympy/geometry/point.py", line 159, in __new__
    raise ValueError('Imaginary coordinates are not permitted.')
ValueError: Imaginary coordinates are not permitted.

图片：

编辑：适用于 SymPy 1.1.2

安装 SymPy 开发版后 (pip install git+https://github.com/sympy/sympy.git) 我得到了正确的结果：

######## Intersections of pairs of planes:
Line3D(Point3D(-7 + sqrt(2)/2, -sqrt(2)/2 + 7, 0), Point3D((1 + sqrt(2))*(-2*sqrt(2) - 1) - 6 + 5*sqrt(2)/2, -5*sqrt(2)/2 + 6 + (-2*sqrt(2) - 1)*(-sqrt(2) - 1), -1 - (1 + sqrt(2))*(-sqrt(2) - 1)))
Line3D(Point3D(-13 - 6*sqrt(2), 0, 0), Point3D(-13 + (1 + sqrt(2))*(-2*sqrt(2) - 1) - (-2*sqrt(2) - 1)*(-sqrt(2) - 1) - 6*sqrt(2), 0, 2 + 2*sqrt(2)))
Line3D(Point3D(-13/2 - 3*sqrt(2), -7*sqrt(2)/2 + 1/2, 0), Point3D(-(-2*sqrt(2) - 1)*(-sqrt(2) - 1) - 15/2 - 5*sqrt(2), -(-2*sqrt(2) - 1)*(-sqrt(2) - 1) - 3*sqrt(2)/2 + 3/2, 1 + (-sqrt(2) - 1)**2))

######## Intersection of all planes:
[Point3D(0, 0, -20*sqrt(2)/7 - 11/7)]

Answer 1

在 SymPy 1.1.1 及更早版本中，当法向量包含根时，intersection returns 会出现错误结果。这是一个更简单的例子：

p1 = Plane((1, 0, 0), (sqrt(2), 0, 0))
p2 = Plane((1, 1, 1), (1, 1, 1))
line = p1.intersection(p2)[0]      # this line is wrong
print(line.arbitrary_point())

这个 returns 线的参数方程 Point3D(3, -sqrt(2)*t, sqrt(2)*t) 这是错误的，因为第一个平面有方程 sqrt(2)*(x-1) = 0，即 x=1 .

您仍然可以找到与

相交的正确方程

solve([p1.equation(), p2.equation()])

但这不会那么容易用于绘图。

好消息

该错误（在 linsolve 方法中）已在当前开发版本 1.1.2.dev 中修复。从 GitHub repo 获取。

早期版本的解决方法

用浮点数替换部首：

pointR = Point3D(-1, N(1 + sqrt(2)), N(-2*sqrt(2) - 1))
pointG = Point3D(N(-sqrt(2) - 1), 1, N(-2*sqrt(2) - 1))
pointB = Point3D(-1, N(-sqrt(2) - 1), N(-2*sqrt(2) - 1))

这不会使一切都变得完美，但应该会减少错误的影响，并且您可能能够为您的图表获得合理的交叉点。

Answer 2

这个错误在正在开发的版本上没有完全解决。特别是关于非简单坐标的线和圆之间的交点。您可以改用 shapely：

[

为什么 SymPy 会计算出错误的平面交点？

Why does SymPy calculate wrong intersections of planes?

python

symbolic-math

sympy

有效

失败

编辑：适用于 SymPy 1.1.2

好消息

早期版本的解决方法