我在物理模拟和 C++ 编程方面拥有专业经验，尽管我没有天体物理学模拟方面的具体经验。
我正在尝试自己构建一个两体进化系统，但我正在努力获得即使是最简单的稳定轨道（不仅仅是几圈）。
因此，我正在寻找一些建议或理想情况下是一些源代码项目作为参考。
我从随机参数开始，因为我不想强制任何东西到位。我已经按照建议使用了来自 boost 的 Runge-Kutta，但现在我怀疑这是否足够。行星往往会相互碰撞或立即分离。我还尝试将初始参数限制在某个合理范围内。
我得到的最好结果是转了几圈，根本不稳定（在轨迹方面也是如此）。
提前致谢。
我不能声称拥有非常深厚的计算专业知识，所以如果它可能对您有所帮助，我可以发布一些代码，其中包含两个版本的固定步长直接辛算法（leap-frog / Verlet 的类型）。不过我对C++不是很精通，所以我可以提供一个python版本，所以也许你可以把它翻译成C++。从好的方面来说，python 在方法论思想和组织方面更加透明，因此可能是有益的。该代码也不具有实际常量，而是缩放版本。
此外，代码设置为平面中的两个质点，因为您指定要从两体模拟开始。我希望这可以让你开始。
'''
Two body leap-forg algorithm simulation
'''
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
'''
Calculation of the gravitational acceleration
r is a 2 x 2 matrix, each row is the position of one of the two mass-points
'''
def g_accel(r, mass1, mass2):
r_12 = r[1,:] - r[0,:]
g_force = mass1 * mass2 * r_12 / np.linalg.norm( r_12 )**3
return np.array([g_force / mass1, - g_force / mass2])
'''
propagation/numerical integration using leap-from / verlet's type time-reversible symplectic algorithm of convergence rate dt**2
'''
def propagate_system(r1_in, r2_in, v1_in, v2_in, mass1, mass2, n_steps, dt):
r = np.empty((n_steps, 2, 2),dtype=float)
v = np.empty((n_steps, 2, 2),dtype=float)
r[0,0,:] = r1_in
r[0,1,:] = r2_in
v[0,0,:] = v1_in
v[0,1,:] = v2_in
for n in range(n_steps-1):
r[n+1,:,:] = r[n,:,:] + dt * v[n,:,:] / 2
v[n+1,:,:] = v[n,:,:] + dt * g_accel(r[n+1,:], mass1, mass2)
r[n+1,:,:] = r[n+1,:,:] + dt * v[n+1,:,:] / 2
return r, v
'''
another propagation/numerical integration using leap-from / verlet's type time-reversible symplectic algorithm of convergence rate dt**4
'''
def propagate_system_1(r1_in, r2_in, v1_in, v2_in, mass1, mass2, n_steps, dt):
w0 = pow(2, 1/3)
w1 = 1 / (2 - w0)
w0 = - w0*w1
c = np.array([w1, w0+w1 , w0+w1, w1]) / 2
d = np.array([w1, w0, w1])
r = np.empty((n_steps, 2, 2),dtype=float)
v = np.empty((n_steps, 2, 2),dtype=float)
r[0,0,:] = r1_in
r[0,1,:] = r2_in
v[0,0,:] = v1_in
v[0,1,:] = v2_in
for n in range(n_steps-1):
r[n+1,:,:] = r[n,:,:]
v[n+1,:,:] = v[n,:,:]
for i in range(3):
r[n+1,:,:] = r[n+1,:,:] + c[i]*dt * v[n+1,:,:]
v[n+1,:,:] = v[n+1,:,:] + d[i]*dt * g_accel(r[n+1,:,:], mass1, mass2)
r[n+1,:,:] = r[n+1,:,:] + c[3]*dt * v[n+1,:,:]
return r, v
'''
Test simulation:
Initial conditions:
'''
mass1 = 2
mass2 = 1
initial_position1 = np.array([-1,0])
initial_position2 = np.array([1.3,0])
initial_velocity1 = np.array([0, 0.5])
initial_velocity2 = - mass1 * initial_velocity1 / mass2
'''
Integration parameters:
'''
t_step = 0.05
N = 20000
'''
System time-propagation (numerical integration):
'''
r, v = propagate_system_1(initial_position1,
initial_position2,
initial_velocity1,
initial_velocity2,
mass1, mass2,
N, t_step)
'''
reducing the frequency of time-measurements for faster simulation speed
'''
r = r[np.arange(1, N, 10), :, :]
'''
animation plot of time-evolution:
'''
plt.style.use('seaborn-whitegrid')
fig = plt.figure()
ax = plt.axes()
ax.set_aspect('equal')
ax.set_xlim(-12, 12)
ax.set_ylim(-12, 12)
line = np.empty(2, dtype=type(ax.plot(r[0, 0, 0], r[0, 0, 1])))
for point in range(2):
line[point], = ax.plot( r[0, point, 0], r[0, point, 1] )
def animate(i):
'''
update plot
'''
for p in range(2):
line[p].set_xdata(r[0:i, p, 0])
line[p].set_ydata(r[0:i, p, 1])
return line
intervals = 50
frames = N
frames = int(frames)
anim = FuncAnimation(fig, animate, frames=frames, interval=intervals)
plt.show()
```