从plotting orbital trajectories,我们有以下代码。该值已更改为可以工作的已知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()
最佳答案
您遇到数字问题的原因有几个:
首先,您不应该要求ODE求解器返回8640000点的数据。
其次,您的参数和初始条件包含大量数字,可以通过引入适当的non-dimensional quantities来消除它们。
设置好之后,下面的代码将产生明智的输出:
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()
结果是
关于python - 执行3D图时python卡住,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/16051523/