在 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

最简单的方法是以参数形式重写，这样您就可以得到表达式 x = A cos(t)； y = B sin(t)。然后通过分配 t = [0, 2 * pi] 并计算相应的 x 和 y.

创建一个完整的椭圆

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

Answer 2

简介

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

每个值的 x, y 对

这里要注意的是您的椭圆已旋转（x*y 耦合项表明了这一点）。要分离方程，您必须首先变换方程以去除耦合项。这与找到一组与椭圆旋转相同角度的轴，沿这些轴进行参数化，然后将结果旋转回来是一样的。查看 this forum post 了解总体概览。

数学

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

$\theta = \frac{1}{2} tan^{-1}(\frac{b}{-1 - a})$

你的等式然后转换为

$a'x'^2+b'y'^2+c'x+d'y'+e'=0$

哪里

$a'=-cos^2(\theta)+b*sin(\theta)cos(\theta)+a*sin^2(\theta)$

$b'=-sin^2(\theta)-b*sin(\theta)cos(\theta)+a*cos^2(\theta)$

$c'=c*cos(\theta)+d*sin(\theta)$

$d'=-c*sin(\theta)+d*cos(\theta)$

$e'=e$

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

$\frac{(x'-h)^2}{a''^2}+\frac{(y'-k)^2}{b''^2}=c''$

哪里

$a''=\frac{1}{\sqrt{a'}}$

$b''=\frac{1}{\sqrt{b'}}$

$c''=\frac{c'^2}{4a'}+\frac{d'^2}{4b'}-e'$

$h=-\frac{c'}{2a'}$

$k=-\frac{d'}{2b'}$

既然你知道 $\theta$ ，你现在可以参数化 $x'$ 和 $y'$ 的方程：

$x'=h+a''\sqrt{c''}sin(t)$

$y'=k+b''\sqrt{c''}cos(t)$

然后您将使用公式

旋转回正常的 x-y space

$x=x'cos(\theta)-y'sin(\theta)$

和

$y=x'sin(\theta)+y'cos(\theta)$

代码

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

结果如下图：

备注

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

Answer 3

y 作为 x 的函数可以求解该问题

要注意的是每个有效的 x 都有 2 个 y 值，并且在 x 的运行ge 之外没有（或虚构的）y 解，椭圆跨越

下面是 3.5 代码，sympy 1.0 应该没问题，但是打印，list comps 可能无法向后兼容 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 到 "feel out" 有效 x 运行ge 或更多数学运算

喜欢下面的try/except：（编辑为"close"绘图重新评论）

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

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

Drawing elliptical orbit in Python (using numpy, matplotlib)

python

matplotlib

orbit

简介

数学

代码

备注