我有一个椭球的一般公式：
A*x**2 + C*y**2 + D*x + E*y + B*x*y + F + G*z**2 = 0
其中A，B，C，D，E，F，G是恒定因子。
如何在matplotlib中将该方程绘制为3D图？ （最好使用线框。）
我看到了这个example，但是它是参数形式，我不确定如何在此代码中放置z坐标。有没有一种方法可以保留一般形式而不用参数形式来绘制它？
我开始将其放在这样的某种代码中：
from mpl_toolkits import mplot3d
%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
def f(x, y):
return ((A*x**2 + C*y**2 + D*x + E*y + B*x*y + F))
def f(z):
return G*z**2
x = np.linspace(-2200, 1850, 30)
y = np.linspace(-100, 60, 30)
z = np.linspace(-100, 60, 30)
X, Y, Z = np.meshgrid(x, y, z)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(X, Y, Z, rstride=10, cstride=10)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z');
我收到此错误：
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-1-95b1296ae6a4> in <module>()
18 fig = plt.figure()
19 ax = fig.add_subplot(111, projection='3d')
---> 20 ax.plot_wireframe(X, Y, Z, rstride=10, cstride=10)
21 ax.set_xlabel('x')
22 ax.set_ylabel('y')
C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\lib\site-packages\mpl_toolkits\mplot3d\axes3d.py in plot_wireframe(self, X, Y, Z, *args, **kwargs)
1847 had_data = self.has_data()
1848 if Z.ndim != 2:
-> 1849 raise ValueError("Argument Z must be 2-dimensional.")
1850 # FIXME: Support masked arrays
1851 X, Y, Z = np.broadcast_arrays(X, Y, Z)
ValueError: Argument Z must be 2-dimensional.
最佳答案
附带说明，但您拥有的并不是3D椭球的最通用方程。您的方程式可以改写为
A*x**2 + C*y**2 + D*x + E*y + B*x*y = - G*z**2 - F,
这意味着实际上对于每个
z值，您都会获得不同级别的2d椭圆，并且这些切片相对于z = 0平面对称。这说明您的椭球体不是一般性的，并有助于检查结果以确保我们得到的结果有意义。假设我们得到一个通用点
r0 = [x0, y0, z0]，r0 @ M @ r0 + b0 @ r0 + c0 == 0
哪里
M = [ A B/2 0
B/2 C 0
0 0 G],
b0 = [D, E, 0],
c0 = F
其中
@ stands for matrix-vector or vector-vector product。您可以使用函数和plot its isosurface，但这不是次优的：您将需要函数的网格近似，这对于获得足够的分辨率非常昂贵，并且您必须明智地选择采样域。
相反，您可以对数据执行principal axis transformation来概括您自己链接的parametric plot of a canonical ellipsoid。
第一步是将
M对角化为M = V @ D @ V.T，其中D是diagonal。由于它是一个实对称矩阵，因此总是可能的，并且V是orthogonal。那我们有r0 @ V @ D @ V.T @ r0 + b0 @ r0 + c0 == 0
我们可以将其重组为
(V.T @ r0) @ D @ (V.T @ r0) + b0 @ V @ (V.T @ r0) + c0 == 0
这激发了辅助坐标
r1 = V.T @ r0和矢量b1 = b0 @ V的定义，为此r1 @ D @ r1 + b1 @ r1 + c0 == 0.
由于
D是一个对称矩阵，其特征值d1, d2, d3在对角线中，因此上式是d1 * x1**2 + d2 * x2**2 + d3 * x3**3 + b11 * x1 + b12 * x2 + b13 * x3 + c0 == 0
其中
r1 = [x1, x2, x3]和b1 = [b11, b12, b13]。剩下的就是从
r1切换到r2，以便我们删除线性项：d1 * (x1 + b11/(2*d1))**2 + d2 * (x2 + b12/(2*d2))**2 + d3 * (x3 + b13/(2*d3))**2 - b11**2/(4*d1) - b12**2/(4*d2) - b13**2/(4*d3) + c0 == 0
所以我们定义
r2 = [x2, y2, z2]
x2 = x1 + b11/(2*d1)
y2 = y1 + b12/(2*d2)
z2 = z1 + b13/(2*d3)
c2 = b11**2/(4*d1) b12**2/(4*d2) b13**2/(4*d3) - c0.
对于这些，我们终于有了
d1 * x2**2 + d2 * y2**2 + d3 * z2**2 == c2,
d1/c2 * x2**2 + d2/c2 * y2**2 + d3/c2 * z2**2 == 1
这是二阶曲面的规范形式。为了使它有意义地对应于椭球，我们必须确保
d1，d2，d3和c2都严格为正。如果可以保证，则标准形式的半长轴为sqrt(c2/d1)，sqrt(c2/d2)和sqrt(c2/d3)。所以这是我们的工作：
确保参数对应于椭球
为极角和方位角生成theta和phi网格
计算转换后的坐标
[x2, y2, z2]将它们移回（按
r2 - r1）以获得[x1, y1, z1]将坐标转换回
V以获得r0，即我们感兴趣的实际[x, y, z]坐标。这是我的实现方式：
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def get_transforms(A, B, C, D, E, F, G):
""" Get transformation matrix and shift for a 3d ellipsoid
Assume A*x**2 + C*y**2 + D*x + E*y + B*x*y + F + G*z**2 = 0,
use principal axis transformation and verify that the inputs
correspond to an ellipsoid.
Returns: (d, V, s) tuple of arrays
d: shape (3,) of semi-major axes in the canonical form
(X/d1)**2 + (Y/d2)**2 + (Z/d3)**2 = 1
V: shape (3,3) of the eigensystem
s: shape (3,) shift from the linear terms
"""
# construct original matrix
M = np.array([[A, B/2, 0],
[B/2, C, 0],
[0, 0, G]])
# construct original linear coefficient vector
b0 = np.array([D, E, 0])
# constant term
c0 = F
# compute eigensystem
D, V = np.linalg.eig(M)
if (D <= 0).any():
raise ValueError("Parameter matrix is not positive definite!")
# transform the shift
b1 = b0 @ V
# compute the final shift vector
s = b1 / (2 * D)
# compute the final constant term, also has to be positive
c2 = (b1**2 / (4 * D)).sum() - c0
if c2 <= 0:
print(b1, D, c0, c2)
raise ValueError("Constant in the canonical form is not positive!")
# compute the semi-major axes
d = np.sqrt(c2 / D)
return d, V, s
def get_ellipsoid_coordinates(A, B, C, D, E, F, G, n_theta=20, n_phi=40):
"""Compute coordinates of an ellipsoid on an ellipsoidal grid
Returns: x, y, z arrays of shape (n_theta, n_phi)
"""
# get canonical grid
theta,phi = np.mgrid[0:np.pi:n_theta*1j, 0:2*np.pi:n_phi*1j]
r2 = np.array([np.sin(theta) * np.cos(phi),
np.sin(theta) * np.sin(phi),
np.cos(theta)]) # shape (3, n_theta, n_phi)
# get transformation data
d, V, s = get_transforms(A, B, C, D, E, F, G) # could be *args I guess
# shift and transform back the coordinates
r1 = d[:,None,None]*r2 - s[:,None,None] # broadcast along first of three axes
r0 = (V @ r1.reshape(3, -1)).reshape(r1.shape) # shape (3, n_theta, n_phi)
return r0 # unpackable to x, y, z of shape (n_theta, n_phi)
这是一个带有椭球的例子，并证明它可以工作：
A,B,C,D,E,F,G = args = 2, -1, 2, 3, -4, -3, 4
x,y,z = get_ellipsoid_coordinates(*args)
print(np.allclose(A*x**2 + C*y**2 + D*x + E*y + B*x*y + F + G*z**2, 0)) # True
从这里开始的实际绘图是微不足道的。使用3d scaling hack from this answer保留相等的轴：
# create 3d axes
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# plot the data
ax.plot_wireframe(x, y, z)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
# scaling hack
bbox_min = np.min([x, y, z])
bbox_max = np.max([x, y, z])
ax.auto_scale_xyz([bbox_min, bbox_max], [bbox_min, bbox_max], [bbox_min, bbox_max])
plt.show()
结果如下所示：
绕其旋转很明显，该表面确实相对于
z = 0平面是反射对称的，这从方程式中可以明显看出。您可以将函数的
n_theta和n_phi关键字参数更改为生成具有不同网格的网格。有趣的是，您可以获取单位球面上的任何分散点，并将其插入函数r2中的get_ellipsoid_coordinates定义中（只要此数组的第一维大小为3），然后输出坐标将具有相同的形状，但是它们将转换为实际的椭球。您也可以使用其他库来可视化表面，例如mayavi，您可以在其中绘制我们刚刚计算的表面，也可以在其中内置compare it with an isosurface。