изображение - 1000 слов
Прямые ошибки в вашем коде:
Вы вычисляетефорсировать в неправильном направлении, это должно быть rx = b[n].x-b[x].x и т. д., или вам нужно удалить знак минус через несколько строк.
При вычислении в отдельных координатах возникают ошибки копирования-вставки, как в
x = int(Bodies[n].x) + int(Bodies[n].vx) * dt
y = int(Bodies[n].y) + int(Bodies[n].vx) * dt
z = int(Bodies[n].z) + int(Bodies[n].vz) * dt
, где в y координате вы все еще используете vx.Промежуточное округление до целочисленных значений не имеет смысла, оно лишь несколько снижает точность.
Я изменил ваш код, чтобы использовать массивы в качестве векторов, отделив вычисление ускорения от Эйлераобновления, убрал бессмысленное округление до целочисленных значений во время численного моделирования, удалил неиспользуемые переменные и поля, удалил промежуточные переменные для вычислений силы / ускорения, чтобы напрямую обновить поле ускорения, изменил цикл, чтобы использовать время, чтобы замечать, когда год(или 10) прошло (ваш код повторяется более 100 лет с шагом 0,1 дня, это было задумано?), ... и добавил Венеру к телам и добавил код для создания изображений, результат см. выше.
Это спирали характерно для метода Эйлера.Вы можете легко улучшить эту модель, изменив обновление Эйлера на симплектическое обновление Эйлера, что означает сначала обновить скорость и вычислить положение с новой скоростью.Это дает, при всем остальном, изображение
day = 60*60*24
# Constants
G = 6.67408e-11
au = 1.496e11
class CelBody(object):
# Constants of nature
# Universal constant of gravitation
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 (kg)
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 Bodies
t = 0
dt = 0.1*day
Bodies = [
CelBody(0, 'Sun', [0, 0, 0], [0, 0, 0], 1.989e30, 'yellow', 10),
CelBody(1, 'Earth', [-1*au, 0, 0], [0, 29783, 0], 5.9742e24, 'blue', 3),
CelBody(2, 'Venus', [0, 0.723 * au, 0], [ 35020, 0, 0], 4.8685e24, 'red', 2),
]
paths = [ [ b.x[:2].copy() ] for b in Bodies]
# loop over ten astronomical years
v = 0
while t < 10*365.242*day:
# 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 = sum(rx**2)**1.5
body.a += -G*other.M*rx/r3
for n, planet in enumerate(Bodies):
# use the symplectic Euler method for better conservation of the constants of motion
planet.v += planet.a*dt
planet.x += planet.v*dt
paths[n].append( planet.x[:2].copy() )
#print("%10s x:%53s v:%53s"%(planet.name,planet.x, planet.v))
if t > v:
print("t=%f"%t)
for b in Bodies: print("%10s %s"%(b.name,b.x))
v += 30.5*day
t += dt
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()