python執行3d圖時會凍結

Question

從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() 

Answer 1

您正在運行到以下幾個原因數值問題： 首先，你不應該問ODE求解器，爲864萬點的返回數據。其次，您的參數和初始條件包含大量數字，您可能通過引入適當的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() 

結果是