在 Python 中绘制椭圆轨道（使用 numpy、matplotlib）

Question

我想知道如何使用方程 ay² + bxy + cx + dy + e = x² 绘制椭圆轨道？

我首先确定了 a、b、c、d、e 常数，现在我假设通过给出 x 值我将获得 y，这将给我想要的图形，但我无法通过使用 matplotlib 来做到这一点。

如果您能帮助我，我将不胜感激！

编辑：我在这里添加了代码。

from numpy import linalg
from numpy import linspace
import numpy as np
from numpy import meshgrid
import random
import matplotlib.pyplot as plt
from scipy import optimize

x = [1.02, 0.95, 0.87, 0.77, 0.67, 0.56, 0.44, 0.30, 0.16, 0.01]
y = [0.39, 0.32, 0.27, 0.22, 0.18, 0.15, 0.13, 0.12, 0.12, 0.15]

my_list = [] #It is the main list.
b = [0] * len(x) # That is the list that contains the results that are given as x^2 from the equation.

def fxn():  # That is the function that solves the given equation to find each parameter.
    global my_list
    global b
    for z in range(len(x)):
        w = [0] * 5
        w[0] = y[z] ** 2
        w[1] = x[z] * y[z]
        w[2] = x[z]
        w[3] = y[z]
        w[4] = 1
        my_list.append(w)
        b[z] = x[z] ** 2

    t = linalg.lstsq(my_list, b)[0]
    print 'List of list representation is', my_list
    print 'x^2, the result of the given equation is', b
    print '\nThe list that contains the parameters is', t

fxn()
t = linalg.lstsq(my_list, b)[0]

print '\nThe constant a is', t[0]
print 'The constant b is', t[1]
print 'The constant c is', t[2]
print 'The constant d is', t[3]
print 'The constant e is', t[4]

编辑：这里是常量值：

a = -4.10267300566
b = 1.10642410023
c = 0.39735696603
d = 3.05101004127
e = -0.370426134994

Answer 1

简介

最简单的方法是参数化这个方程。正如@Escualo 建议的那样，您可以引入一个变量t 并参数化x 和y。参数化意味着根据t，将您的方程分成两个单独的方程x 和y。因此，对于t 的某些值，您将拥有x = f(t) 和y = g(t)。然后，您可以绘制每个 t 值的 x, y 对。

这里的问题是你的椭圆是旋转的（x*y 耦合项就是一个指示）。要分离方程，您必须首先变换方程以消除耦合项。这与找到一组与椭圆旋转相同角度的轴相同，沿这些轴进行参数化，然后将结果旋转回来。查看this forum post 了解总体概况。

数学

您首先需要找到椭圆轴相对于 x-y 坐标平面的旋转角度。

然后你的方程转换为

在哪里

要找到椭圆的（几乎）标准形式，您可以完成 和 部分的正方形并稍微重新排列方程：

在哪里

既然您知道，您现在可以参数化 和 的方程：

然后您将使用公式旋转回正常的 x-y 空间

和

代码

将 x 和 y 数组传递给 plt.plot 的代码现在相对简单：

def computeEllipse(a, b, c, d, e):
    """
    Returns x-y arrays for ellipse coordinates.
    Equation is of the form a*y**2 + b*x*y + c*x + d*y + e = x**2
    """
    # Convert x**2 coordinate to +1
    a = -a
    b = -b
    c = -c
    d = -d
    e = -e
    # Rotation angle
    theta = 0.5 * math.atan(b / (1 - a))
    # Rotated equation
    sin = math.sin(theta)
    cos = math.cos(theta)
    aa = cos**2 + b * sin * cos + a * sin**2
    bb = sin**2 - b * cos * sin + a * cos**2
    cc = c * cos + d * sin
    dd = -c * sin + d * cos
    ee = e
    # Standard Form
    axMaj = 1 / math.sqrt(aa)
    axMin = 1 / math.sqrt(bb)
    scale = math.sqrt(cc**2 / (4 * aa) + dd**2 / (4 * bb) - ee)
    h = -cc / (2 * aa)
    k = -dd / (2 * bb)
    # Parametrized Equation
    t = np.linspace(0, 2 * math.pi, 1000)
    xx = h + axMaj * scale * np.sin(t)
    yy = k + axMin * scale * np.cos(t)
    # Un-rotated coordinates
    x = xx * cos - yy * sin
    y = xx * sin + yy * cos

    return x, y

实际使用代码：

from matplotlib import pyplot as plt

a = -4.10267300566
b = 1.10642410023
c = 0.39735696603
d = 3.05101004127
e = -0.370426134994

lines = plt.plot(*computeEllipse(a, b, c, d, e))

在椭圆上绘制原始点：

x = [1.02, 0.95, 0.87, 0.77, 0.67, 0.56, 0.44, 0.30, 0.16, 0.01]
y = [0.39, 0.32, 0.27, 0.22, 0.18, 0.15, 0.13, 0.12, 0.12, 0.15]
ax = lines[0].axes
ax.plot(x, y, 'r.')

结果如下图：

注意

请记住，我链接到的论坛帖子使用的符号与您使用的不同。他们的方程是Ax² + Bxy + Cy² + Dx + Ey + F = 0。这比您的 ay² + bxy - x² + cx + dy + e = 0 形式更标准一些。所有的数学都是根据你的符号来计算的。

Answer 2

问题可以通过 y 作为 x 的函数来解决

问题是每个有效 x 都有 2 个 y 值，并且在椭圆跨越的 x 范围之外没有（或假想的）y 解

下面是 3.5 代码，sympy 1.0 应该没问题，但是打印，列表组合可能无法向后兼容 2.x

from numpy import linalg
from numpy import linspace
import numpy as np
from numpy import meshgrid
import random
import matplotlib.pyplot as plt
from scipy import optimize
from sympy import *

xs = [1.02, 0.95, 0.87, 0.77, 0.67, 0.56, 0.44, 0.30, 0.16, 0.01]
ys = [0.39, 0.32, 0.27, 0.22, 0.18, 0.15, 0.13, 0.12, 0.12, 0.15]

b = [i ** 2 for i in xs] # That is the list that contains the results that are given as x^2 from the equation.

def fxn(x, y):  # That is the function that solves the given equation to find each parameter.
    my_list = [] #It is the main list.
    for z in range(len(x)):
        w = [0] * 5
        w[0] = y[z] ** 2
        w[1] = x[z] * y[z]
        w[2] = x[z]
        w[3] = y[z]
        w[4] = 1
        my_list.append(w)
    return my_list

t = linalg.lstsq(fxn(xs, ys), b)


def ysolv(coeffs):
    x,y,a,b,c,d,e = symbols('x y a b c d e')
    ellipse = a*y**2 + b*x*y + c*x + d*y + e - x**2
    y_sols = solve(ellipse, y)
    print(*y_sols, sep='\n')

    num_coefs = [(a, f) for a, f in (zip([a,b,c,d,e], coeffs))]
    y_solsf0 = y_sols[0].subs(num_coefs)
    y_solsf1 = y_sols[1].subs(num_coefs)

    f0 = lambdify([x], y_solsf0)
    f1 = lambdify([x], y_solsf1)
    return f0, f1

f0, f1 = ysolv(t[0])

y0 = [f0(x) for x in xs]
y1 = [f1(x) for x in xs]

plt.scatter(xs, ys)
plt.scatter(xs, y0, s=100, color = 'red', marker='+')
plt.scatter(xs, y1, s=100, color = 'green', marker='+')
plt.show()  

当上面在 Spyder 中运行时：

runfile('C:/Users/john/mypy/mySE_answers/ellipse.py', wdir='C:/Users/john/mypy/mySE_answers')
(-b*x - d + sqrt(-4*a*c*x - 4*a*e + 4*a*x**2 + b**2*x**2 + 2*b*d*x + d**2))/(2*a)
-(b*x + d + sqrt(-4*a*c*x - 4*a*e + 4*a*x**2 + b**2*x**2 + 2*b*d*x + d**2))/(2*a)

 为 y 值生成的函数并非在任何地方都有效：

f0(0.1), f1(0.1)
Out[5]: (0.12952825130864626, 0.6411040771593166)

f0(2)
Traceback (most recent call last):

  File "<ipython-input-6-9ce260237dcd>", line 1, in <module>
    f0(2)

  File "<string>", line 1, in <lambda>

ValueError: math domain error


In [7]:

域错误需要 try/execpt 来“检测”有效的 x 范围或更多数学运算

喜欢下面的尝试/除了：（编辑为“关闭”绘图重新评论）

def feeloutXrange(f, midx, endx):
    fxs = []
    x = midx
    while True:
        try: f(x)
        except:
            break
        fxs.append(x)
        x += (endx - midx)/100
    return fxs

midx = (min(xs) + max(xs))/2    

xpos = feeloutXrange(f0, midx, max(xs))
xnegs = feeloutXrange(f0, midx, min(xs))
xs_ellipse = xnegs[::-1] + xpos[1:]

y0s = [f0(x) for x in xs_ellipse]
y1s = [f1(x) for x in xs_ellipse]

ys_ellipse = y0s + y1s[::-1] + [y0s[0]] # add y start point to end to close drawing

xs_ellipse = xs_ellipse + xs_ellipse[::-1] + [xs_ellipse[0]] # added x start point


plt.scatter(xs, ys)
plt.scatter(xs, y0, s=100, color = 'red', marker='+')
plt.scatter(xs, y1, s=100, color = 'green', marker='+')
plt.plot(xs_ellipse, ys_ellipse)
plt.show()

编辑：在椭圆点列表的末尾添加重复的起点以关闭绘图图

ys_ellipse = y0s + y1s[::-1] + [y0s[0]] # add y start point to end to close drawing

xs_ellipse = xs_ellipse + xs_ellipse[::-1] + [xs_ellipse[0]] # added x start point

Answer 3

最简单的方法是以参数形式重写，这样你就可以得到表达式x = A cos(t)； y = B sin(t)。然后通过分配t = [0, 2 * pi] 并计算对应的x 和y 来创建一个完整的椭圆。

阅读this article，其中解释了如何从一般二次形式转变为参数形式。