1

通过绘制轨道轨迹，我们有以下代码。这些值已更改为可以工作的已知 IC。

如果这段代码是正确的（虽然不可能），它会生成 在此处输入图像描述

运行此代码只会冻结我的计算机或输出绝对错误的图。有人可以帮我找到解决此问题的方法吗？

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from numpy import linspace
from mpl_toolkits.mplot3d import Axes3D

mu = 398600
# r0 = [-149.6 * 10 ** 6, 0.0, 0.0]  #  Initial position
# v0 = [29.9652, -5.04769, 0.0]      #  Initial velocity
u0 = [-4069.503, 2861.786, 4483.608, -5.114, -5.691, -5]


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

dt = np.arange(0.0, 24 * 3600, .01)   #  Time to run code in seconds'
u = odeint(deriv, u0, dt)
x, y, z, x2, y2, z2 = z.T

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x, y, z)
plt.show()
4

1 回答 1

4

您遇到数值问题有几个原因：首先，您不应该要求 ODE 求解器返回 8640000 点的数据。其次，您的参数和初始条件包含大量数字，您可以通过引入适当的无量纲量来摆脱它们。

设置后，下面的代码会产生合理的输出：

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from numpy import linspace
from mpl_toolkits.mplot3d import Axes3D

u0 = [0, 0, 1, 0, -1, 0]
mu = .1


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

times = np.linspace(0.0, 200, 100)
u = odeint(deriv, u0, times)
x, y, z, x2, y2, z2 = u.T

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x, y, z)
plt.show()

结果是结果

于 2013-04-17T08:14:14.070 回答