在Python中绘制轨道轨迹

Question

我如何在Python中设置三体问题？如何定义函数来解决ODEs？

三个方程式为： 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公里，火星距离太阳227.9 * 10 ** 6公里。

#!/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

- dustin

---

你具体在问什么？如果你有更具体的问题，人们可能更愿意提供帮助。 - askewchan

@askewchan 我不知道如何将6个一阶ODE输入为一个函数，因为integrate.odeint只能接受一阶方程。 - dustin

你可以使用 integrate.ode 来解决高阶方程。 - askewchan

@askewchan 如果你将它们作为一阶函数编写，是可以的。但我不知道如何在此之后定义该函数。 - dustin

请参考 http://www.scipy.org/Cookbook/CoupledSpringMassSystem 获取另一个示例。 - Warren Weckesser

顺便提一下，不要写成-149.6 * 10 ** 6，直接写成-149.6e6即可 :) - astrojuanlu

Answer 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

- askewchan

---

@dustin 你最好单独提出一个问题来询问。 - askewchan

我遇到了一个值错误，提示有太多的值需要解包（期望是6个）。 - dustin

@Juanlu001，那是因为我一开始漏掉了.T。 - askewchan