在 python 中绘制轨道轨迹答案

【问题标题】：plotting orbital trajectories in python在 python 中绘制轨道轨迹
【发布时间】：2013-04-09 14:04:57
【问题描述】：

如何在 python 中设置三体问题？如何定义求解 ODE 的函数？

这三个方程是
x'' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * x、
y'' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * y和
z'' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * z。

写成 6 一阶我们有

x' = x2,

y' = y2,

z' = z2,

x2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * x,

y2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * y，和

z2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * z

我还想添加地球轨道和火星的路径，我们可以假设它是圆形的。地球距离太阳149.6 * 10 ** 6 km，距离火星227.9 * 10 ** 6 km。

#!/usr/bin/env python                                                             
#  This program solves the 3 Body Problem numerically and plots the trajectories      

import pylab
import numpy as np
import scipy.integrate as integrate
import matplotlib.pyplot as plt
from numpy import linspace

mu = 132712000000  #gravitational parameter
r0 = [-149.6 * 10 ** 6, 0.0, 0.0]
v0 = [29.0, -5.0, 0.0]
dt = np.linspace(0.0, 86400 * 700, 5000)  # time is seconds

【问题讨论】：

您特别想问什么？如果您有更具体的问题，人们可能会更有可能提供帮助。
@askewchan 我不知道如何输入 6 个一阶 ODE 作为函数，因为 integrate.odeint 只会采用一阶方程。
您可以使用integrate.ode 求解高阶方程。
@askewchan 是的，如果您将它们写为第一顺序。但是我不知道之后如何定义函数。
顺便说一句，不要写-149.6 * 10 ** 6，而是写-149.6e6 :)

标签： python numpy matplotlib scipy differential-equations

【解决方案1】：

正如您所展示的，您可以将其写为六个一阶颂歌的系统：

x' = x2
y' = y2
z' = z2
x2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * x
y2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * y
z2' = -mu / np.sqrt(x ** 2 + y ** 2 + z ** 2) * z

您可以将其保存为矢量：

u = (x, y, z, x2, y2, z2)

并因此创建一个返回其导数的函数：

def deriv(u, t):
    n = -mu / np.sqrt(u[0]**2 + u[1]**2 + u[2]**2)
    return [u[3],      # u[0]' = u[3]
            u[4],      # u[1]' = u[4]
            u[5],      # u[2]' = u[5]
            u[0] * n,  # u[3]' = u[0] * n
            u[1] * n,  # u[4]' = u[1] * n
            u[2] * n]  # u[5]' = u[2] * n

给定一个初始状态u0 = (x0, y0, z0, x20, y20, z20)，以及一个时间变量t，这可以像这样输入scipy.integrate.odeint：

u = odeint(deriv, u0, t)

u 将是上面的列表。或者你可以从一开始就解压u，忽略x2、y2和z2的值（你必须先用.T转置输出）

x, y, z, _, _, _ = odeint(deriv, u0, t).T

【讨论】：

@dustin 您可能应该尝试作为一个单独的问题提出。
我收到值错误太多值无法解包（预计 6 个）
@Juanlu001 那是因为我一开始忽略了.T。