我一直在嘗試解決兩體問題，將來應該將其用于更多行星，但它不起作用，我應該得到的情節是圓形的，但我收到了兩體的直線系統。有誰知道我該如何解決這個問題并獲得正確的情節？
這是我使用的代碼：
import numpy as np
import matplotlib.pyplot as plt
aEarth = 1 #semi-major axis (AU)
eEarth = 0.0167 #eccentricity (no units)
aMercury = 0.387098 #semi-major axis (AU)
eMercury = 0.205635 #eccentricity (no units)
Msun = 1 #Mass of Sun (Solar Mass)
Mearth = 3.0024584*10**(-6) #Mass of Earth (Solar Mass)
Mmercury = 1.65956463*10**(-7) #Mass of Mercury (Solar Mass)
Mes = (Msun Mearth)
Mms = (Msun Mmercury)
G = 1 #Gravitational Constant
apocentreEarth = (aEarth*(1 eEarth))
apocentreMercury = (aMercury*(1 eMercury))
vapocentreEarth = (np.sqrt((G*(Mearth Msun)/aEarth)*((1-eEarth)/(1 eEarth))))
vapocentreMercury = (np.sqrt((G*(Mmercury Msun)/aMercury)*(1-eMercury/1 eMercury)))
xEarth = apocentreEarth
yEarth = 0.0
zEarth = 0.0
vxEarth = 0.0
vyEarth =(vapocentreEarth)
vzEarth = 0.0
xMercury = apocentreMercury
yMercury = 0.0
zMercury = 0.0
vxMercury = 0.0
vyMercury =(vapocentreMercury)
vzMercury = 0.0
class CelBody(object):
# Constants of nature
def __init__(self, id, name, x0, v0, mass, color, lw):
# Name of the body (string)
self.id = id
self.name = name
# Mass of the body (solar mass)
self.M = mass
# Initial position of the body (au)
self.x0 = np.asarray(x0, dtype=float)
# Position (au). Set to initial value.
self.x = self.x0.copy()
# Initial velocity of the body (au/s)
self.v0 = np.asarray(v0, dtype=float)
# Velocity (au/s). Set to initial value.
self.v = self.v0.copy()
self.a = np.zeros([3], dtype=float)
self.color = color
self.lw = lw
# All Celestial Objects
t = 0
dt = 0.01
Bodies = [
CelBody(0, 'Sun', [0, 0, 0], [0, 0, 0], Msun, 'yellow', 10),
CelBody(1, 'Earth', [xEarth, yEarth, zEarth], [vxEarth, vyEarth, vzEarth], Mearth, 'blue', 3),
CelBody(2, 'Mercury', [xMercury, yMercury, zMercury], [ vxMercury, vyMercury, vzMercury], Mmercury, 'pink', 3),
]
paths = [ [ b.x[:2].copy() ] for b in Bodies]
# loop over ten astronomical years
v = 0
for t in range(0,1000):
# compute forces/accelerations
for body in Bodies:
body.a *= 0
for other in Bodies:
# no force on itself
if (body == other): continue # jump to next loop
rx = body.x - other.x
r3 = (np.sqrt(rx[0]**2 rx[1]**2 rx[2]**2))**3
body.a = -G*other.M*rx/r3
for n, planet in enumerate(Bodies):
# use the Forward Euler algorithm
planet.a = -G*other.M*rx/r3
planet.v = planet.v planet.a*dt
planet.x = planet.x planet.v*dt
paths[n].append( planet.x[:2].copy() )
#print("s x:Ss v:Ss"%(planet.name,planet.x, planet.v))
plt.figure(figsize=(8,8))
for n, planet in enumerate(Bodies):
px, py=np.array(paths[n]).T;
plt.plot(px, py, color=planet.color, lw=planet.lw)
plt.show()
uj5u.com熱心網友回復：
這幾行有問題：
planet.v = planet.v planet.a*dt
planet.x = planet.x planet.v*dt
它們相當于
planet.v = 2*planet.v planet.a*dt
planet.x = 2*planet.x planet.v*dt
這不是你想要的。
要么不使用 =，要么將這些行更改為：
planet.v = planet.a*dt
planet.x = planet.v*dt
仍然存在一個問題：第一行更改planet.v，然后第二行使用new planet.v進行更新planet.x，這不是顯式歐拉積分應該如何進行的。（代碼中有一條關于使用辛歐拉的注釋，但對此答案的注釋說您打算使用前向（即顯式）歐拉。）
對于這個系統，一個簡單的解決方法是切換陳述句：
planet.x = planet.v*dt
planet.v = planet.a*dt
可能還有其他問題。如果您需要更多幫助，請盡可能簡化代碼以創建最小的可重現示例。就像現在一樣，看起來宣告了很多不相關的變數，計算歐拉方法的兩個不同的地方，你分配的兩個地方，dt當你說你正在解決兩體問題時定義的三個體，等等。
編輯更新后，還有幾個錯誤：
在開始的回圈
for body in Bodies:中，您在其中計算body.a每個主體，結果應該被累積，所以更新的行body.a應該是body.a = -G*other.M*rx/r3（注意更改為
=）。在您應用歐拉方法步驟的第二個內部回圈 (
for n, planet in enumerate(Bodies):) 中，洗掉該行planet.a = -G*other.M*rx/r3您已經在上一個回圈中計算了加速度并將它們存盤在
planet.a.
如果我僅對問題中的代碼進行這兩項更改，則該圖將顯示您期望的圓形軌道。
轉載請註明出處，本文鏈接：https://www.uj5u.com/caozuo/441558.html