Как мога да направя 3D диаграма в matplotlib на елипсоид, дефиниран от квадратно уравнение?

Странична бележка, но това, което имате, не е най-общото уравнение за 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

където @ означава матрица-вектор или вектор-вектор продукт.

Можете да вземете функцията си и да начертайте нейната изоповърхност, но това би било неоптимално: ще ви трябва мрежово приближение за вашата функция, което е много скъпо да се направи с достатъчна разделителна способност, и ще трябва да изберете домейна за това вземане на проби разумно.

Вместо това можете да извършите трансформация на главната ос върху вашите данни, за да обобщите параметричен график на каноничен елипсоид, който вие сами сте свързали.

Първата стъпка е да диагонализирате M като M = V @ D @ V.T, където D е диагонал. Тъй като това е реална симетрична матрица, това винаги е възможно и V е ортогонална. Тогава имаме

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).

И така, ето какво правим:

1. уверете се, че параметрите съответстват на елипсоид
2. генерирайте тета и фи мрежа за полярни и азимутални ъгли
3. изчислете трансформираните координати [x2, y2, z2]
4. преместете ги назад (с r2 - r1), за да получите [x1, y1, z1]
5. трансформирайте координатите обратно с 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 мащабиране от този отговор за запазване на равни оси:

# 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, където можете или да начертаете повърхността, която току-що изчислихме, или да я сравните с isosurface, който е вграден там.