426 lines
16 KiB
Python
426 lines
16 KiB
Python
"""
|
|
随机初始条件示例 - 探索不同的三体系统
|
|
"""
|
|
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
import sys
|
|
import os
|
|
|
|
# 添加父目录到路径
|
|
sys.path.append(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
|
|
|
|
from three_body_problem import ThreeBodySolver, ThreeBodyConfig, ThreeBodyVisualizer
|
|
|
|
|
|
def run_random_example(seed: int = 42, total_time: float = 20.0):
|
|
"""
|
|
运行随机初始条件示例
|
|
|
|
参数:
|
|
seed: 随机种子
|
|
total_time: 总模拟时间（年）
|
|
"""
|
|
np.random.seed(seed)
|
|
|
|
print("=" * 60)
|
|
print(f"随机初始条件示例 (种子: {seed})")
|
|
print("=" * 60)
|
|
|
|
# 创建随机配置
|
|
particles = ThreeBodyConfig.create_random_config(
|
|
masses=None, # 使用随机质量
|
|
position_range=2.0, # 位置范围 ±2 AU
|
|
velocity_scale=3.0 # 速度缩放因子
|
|
)
|
|
|
|
# 打印配置摘要
|
|
ThreeBodyConfig.print_config_summary(particles)
|
|
|
|
# 创建求解器
|
|
dt = 0.001 # 时间步长（年）
|
|
solver = ThreeBodySolver(particles, dt=dt)
|
|
|
|
print(f"\n开始模拟，总时间: {total_time}年，时间步长: {dt}年")
|
|
|
|
solver.simulate(total_time=total_time, progress_interval=2000)
|
|
|
|
# 计算守恒误差
|
|
momentum_error, angular_momentum_error, energy_error = solver.get_conservation_errors()
|
|
print(f"\n守恒定律误差:")
|
|
print(f" 动量误差: {momentum_error:.6e}")
|
|
print(f" 角动量误差: {angular_momentum_error:.6e}")
|
|
print(f" 能量相对误差: {energy_error:.6e}")
|
|
|
|
# 分析系统行为
|
|
analyze_system_behavior(solver)
|
|
|
|
# 可视化
|
|
print("\n生成可视化图形...")
|
|
visualize_random_system(solver, seed)
|
|
|
|
return solver
|
|
|
|
|
|
def analyze_system_behavior(solver: ThreeBodySolver):
|
|
"""分析三体系统的行为"""
|
|
print("\n" + "-" * 40)
|
|
print("系统行为分析")
|
|
print("-" * 40)
|
|
|
|
trajectories = solver.get_trajectories()
|
|
|
|
# 计算每个质点的运动范围
|
|
for i, (traj, particle) in enumerate(zip(trajectories, solver.particles)):
|
|
pos_range = np.ptp(traj, axis=0) # 位置范围 (max - min)
|
|
avg_speed = np.mean(np.linalg.norm(np.gradient(traj, solver.dt, axis=0), axis=1))
|
|
|
|
print(f"\n质点 {i+1} ({particle.name}):")
|
|
print(f" 质量: {particle.mass:.4f} M_sun")
|
|
print(f" 位置范围: X={pos_range[0]:.3f}, Y={pos_range[1]:.3f}, Z={pos_range[2]:.3f} AU")
|
|
print(f" 平均速度: {avg_speed:.3f} AU/年")
|
|
|
|
# 计算质点之间的最小距离
|
|
min_distances = []
|
|
for i in range(3):
|
|
for j in range(i+1, 3):
|
|
traj_i = trajectories[i]
|
|
traj_j = trajectories[j]
|
|
distances = np.linalg.norm(traj_i - traj_j, axis=1)
|
|
min_dist = np.min(distances)
|
|
min_distances.append((i, j, min_dist))
|
|
|
|
print(f"\n质点间最小距离:")
|
|
for i, j, min_dist in min_distances:
|
|
print(f" 质点{i+1}-质点{j+1}: {min_dist:.4f} AU")
|
|
|
|
# 检查是否有碰撞或近距离接近
|
|
collision_threshold = 0.1 # AU
|
|
close_encounters = [(i, j, d) for i, j, d in min_distances if d < collision_threshold]
|
|
|
|
if close_encounters:
|
|
print(f"\n警告: 检测到近距离接近!")
|
|
for i, j, d in close_encounters:
|
|
print(f" 质点{i+1}和质点{j+1}的最小距离: {d:.4f} AU < {collision_threshold} AU")
|
|
else:
|
|
print(f"\n系统稳定: 所有质点间距离都大于 {collision_threshold} AU")
|
|
|
|
# 计算系统质心运动
|
|
com_trajectory = []
|
|
for t in range(len(trajectories[0])):
|
|
com = np.zeros(3)
|
|
total_mass = 0.0
|
|
for i, traj in enumerate(trajectories):
|
|
com += solver.particles[i].mass * traj[t]
|
|
total_mass += solver.particles[i].mass
|
|
com_trajectory.append(com / total_mass)
|
|
|
|
com_trajectory = np.array(com_trajectory)
|
|
com_range = np.ptp(com_trajectory, axis=0)
|
|
print(f"\n系统质心运动范围: X={com_range[0]:.4f}, Y={com_range[1]:.4f}, Z={com_range[2]:.4f} AU")
|
|
|
|
|
|
def visualize_random_system(solver: ThreeBodySolver, seed: int):
|
|
"""可视化随机三体系统"""
|
|
trajectories = solver.get_trajectories()
|
|
|
|
# 创建图形
|
|
fig = plt.figure(figsize=(16, 12))
|
|
|
|
# 1. 3D轨迹图
|
|
ax1 = plt.subplot(2, 3, 1, projection='3d')
|
|
|
|
colors = ['red', 'green', 'blue']
|
|
for i, (traj, particle) in enumerate(zip(trajectories, solver.particles)):
|
|
color = particle.color if particle.color else colors[i % len(colors)]
|
|
label = particle.name if particle.name else f"质点 {i+1}"
|
|
|
|
ax1.plot(traj[:, 0], traj[:, 1], traj[:, 2],
|
|
color=color, alpha=0.7, linewidth=1.5, label=label)
|
|
ax1.scatter(traj[-1, 0], traj[-1, 1], traj[-1, 2],
|
|
color=color, s=100, edgecolors='black', linewidth=1.5)
|
|
|
|
# 绘制质心轨迹
|
|
com_trajectory = []
|
|
for t in range(len(trajectories[0])):
|
|
com = np.zeros(3)
|
|
total_mass = 0.0
|
|
for i, traj in enumerate(trajectories):
|
|
com += solver.particles[i].mass * traj[t]
|
|
total_mass += solver.particles[i].mass
|
|
com_trajectory.append(com / total_mass)
|
|
|
|
com_trajectory = np.array(com_trajectory)
|
|
ax1.plot(com_trajectory[:, 0], com_trajectory[:, 1], com_trajectory[:, 2],
|
|
'k--', alpha=0.5, linewidth=1, label='质心轨迹')
|
|
ax1.scatter(com_trajectory[-1, 0], com_trajectory[-1, 1], com_trajectory[-1, 2],
|
|
color='black', marker='x', s=200, label='质心', linewidth=2)
|
|
|
|
ax1.set_xlabel('X (AU)', fontsize=12)
|
|
ax1.set_ylabel('Y (AU)', fontsize=12)
|
|
ax1.set_zlabel('Z (AU)', fontsize=12)
|
|
ax1.set_title('3D轨迹图', fontsize=14, fontweight='bold')
|
|
ax1.legend(fontsize=10)
|
|
ax1.grid(True, alpha=0.3)
|
|
|
|
# 2. XY平面投影
|
|
ax2 = plt.subplot(2, 3, 2)
|
|
for i, (traj, particle) in enumerate(zip(trajectories, solver.particles)):
|
|
color = particle.color if particle.color else colors[i % len(colors)]
|
|
label = particle.name if particle.name else f"质点 {i+1}"
|
|
ax2.plot(traj[:, 0], traj[:, 1], color=color, alpha=0.7, linewidth=1.5, label=label)
|
|
ax2.scatter(traj[-1, 0], traj[-1, 1], color=color, s=100, edgecolors='black', linewidth=1.5)
|
|
|
|
ax2.plot(com_trajectory[:, 0], com_trajectory[:, 1], 'k--', alpha=0.5, linewidth=1)
|
|
ax2.scatter(com_trajectory[-1, 0], com_trajectory[-1, 1],
|
|
color='black', marker='x', s=200, linewidth=2)
|
|
|
|
ax2.set_xlabel('X (AU)', fontsize=12)
|
|
ax2.set_ylabel('Y (AU)', fontsize=12)
|
|
ax2.set_title('XY平面投影', fontsize=14, fontweight='bold')
|
|
ax2.legend(fontsize=10)
|
|
ax2.grid(True, alpha=0.3)
|
|
ax2.set_aspect('equal', adjustable='box')
|
|
|
|
# 3. 距离随时间变化
|
|
ax3 = plt.subplot(2, 3, 3)
|
|
|
|
time_points = np.arange(len(trajectories[0])) * solver.dt
|
|
|
|
# 计算所有质点对之间的距离
|
|
distances = []
|
|
labels = []
|
|
for i in range(3):
|
|
for j in range(i+1, 3):
|
|
dist = np.linalg.norm(trajectories[i] - trajectories[j], axis=1)
|
|
distances.append(dist)
|
|
labels.append(f"质点{i+1}-质点{j+1}")
|
|
|
|
for dist, label in zip(distances, labels):
|
|
ax3.plot(time_points, dist, linewidth=1.5, alpha=0.8, label=label)
|
|
|
|
ax3.set_xlabel('时间 (年)', fontsize=12)
|
|
ax3.set_ylabel('距离 (AU)', fontsize=12)
|
|
ax3.set_title('质点间距离变化', fontsize=14, fontweight='bold')
|
|
ax3.legend(fontsize=10)
|
|
ax3.grid(True, alpha=0.3)
|
|
|
|
# 4. 速度大小随时间变化
|
|
ax4 = plt.subplot(2, 3, 4)
|
|
|
|
for i, (traj, particle) in enumerate(zip(trajectories, solver.particles)):
|
|
color = particle.color if particle.color else colors[i % len(colors)]
|
|
label = particle.name if particle.name else f"质点 {i+1}"
|
|
|
|
# 计算速度大小（使用位置差分）
|
|
if len(traj) > 1:
|
|
velocities = np.gradient(traj, solver.dt, axis=0)
|
|
speed = np.linalg.norm(velocities, axis=1)
|
|
ax4.plot(time_points, speed, color=color, linewidth=1.5, alpha=0.8, label=label)
|
|
|
|
ax4.set_xlabel('时间 (年)', fontsize=12)
|
|
ax4.set_ylabel('速度大小 (AU/年)', fontsize=12)
|
|
ax4.set_title('质点速度变化', fontsize=14, fontweight='bold')
|
|
ax4.legend(fontsize=10)
|
|
ax4.grid(True, alpha=0.3)
|
|
|
|
# 5. 能量分布饼图
|
|
ax5 = plt.subplot(2, 3, 5)
|
|
|
|
# 计算每个质点的动能和势能
|
|
kinetic_energies = []
|
|
potential_energies = []
|
|
|
|
for i, particle in enumerate(solver.particles):
|
|
# 动能
|
|
v_squared = np.sum(particle.velocity**2)
|
|
kinetic_energy = 0.5 * particle.mass * v_squared
|
|
kinetic_energies.append(kinetic_energy)
|
|
|
|
# 势能（与其他质点的相互作用）
|
|
potential_energy = 0.0
|
|
for j, other in enumerate(solver.particles):
|
|
if i != j:
|
|
r_vec = other.position - particle.position
|
|
r = np.linalg.norm(r_vec)
|
|
if r > 1e-10:
|
|
potential_energy -= ThreeBodySolver.G * particle.mass * other.mass / r
|
|
potential_energies.append(potential_energy)
|
|
|
|
# 只考虑势能的一半（每对质点计算了两次）
|
|
potential_energies = [pe/2 for pe in potential_energies]
|
|
|
|
labels = [f"质点{i+1}" for i in range(3)]
|
|
colors_pie = ['lightcoral', 'lightgreen', 'lightblue']
|
|
|
|
ax5.pie(kinetic_energies, labels=labels, autopct='%1.1f%%',
|
|
colors=colors_pie, startangle=90)
|
|
ax5.set_title('动能分布', fontsize=14, fontweight='bold')
|
|
|
|
# 6. 相空间图（所有质点）
|
|
ax6 = plt.subplot(2, 3, 6)
|
|
|
|
for i, (traj, particle) in enumerate(zip(trajectories, solver.particles)):
|
|
color = particle.color if particle.color else colors[i % len(colors)]
|
|
label = particle.name if particle.name else f"质点 {i+1}"
|
|
|
|
if len(traj) > 1:
|
|
velocities = np.gradient(traj, solver.dt, axis=0)
|
|
x_positions = traj[:, 0]
|
|
x_velocities = velocities[:, 0]
|
|
|
|
# 使用颜色表示时间
|
|
scatter = ax6.scatter(x_positions, x_velocities, c=time_points,
|
|
cmap='viridis', alpha=0.6, s=10, label=label)
|
|
|
|
plt.colorbar(scatter, ax=ax6, label='时间 (年)')
|
|
ax6.set_xlabel('X 位置 (AU)', fontsize=12)
|
|
ax6.set_ylabel('X 速度 (AU/年)', fontsize=12)
|
|
ax6.set_title('相空间图 (X维度)', fontsize=14, fontweight='bold')
|
|
ax6.legend(fontsize=10)
|
|
ax6.grid(True, alpha=0.3)
|
|
|
|
plt.suptitle(f'随机三体系统 (种子: {seed})', fontsize=16, fontweight='bold')
|
|
plt.tight_layout()
|
|
|
|
# 保存图形
|
|
output_file = f"random_system_seed_{seed}.png"
|
|
plt.savefig(output_file, dpi=300, bbox_inches='tight')
|
|
print(f"\n图形已保存到: {output_file}")
|
|
|
|
# 显示图形
|
|
plt.show()
|
|
|
|
|
|
def run_multiple_random_simulations(n_simulations: int = 5, total_time: float = 10.0):
|
|
"""运行多个随机模拟并比较结果"""
|
|
print("=" * 60)
|
|
print(f"运行 {n_simulations} 个随机三体系统模拟")
|
|
print("=" * 60)
|
|
|
|
results = []
|
|
|
|
for sim_idx in range(n_simulations):
|
|
seed = 100 + sim_idx # 不同的随机种子
|
|
np.random.seed(seed)
|
|
|
|
print(f"\n模拟 {sim_idx+1}/{n_simulations} (种子: {seed})")
|
|
|
|
# 创建随机配置
|
|
particles = ThreeBodyConfig.create_random_config(
|
|
masses=None,
|
|
position_range=2.0,
|
|
velocity_scale=2.0 + np.random.random() * 2.0 # 随机速度缩放
|
|
)
|
|
|
|
# 创建求解器
|
|
solver = ThreeBodySolver(particles, dt=0.001)
|
|
solver.simulate(total_time=total_time, progress_interval=5000)
|
|
|
|
# 分析结果
|
|
trajectories = solver.get_trajectories()
|
|
|
|
# 计算系统特性
|
|
final_distances = []
|
|
for i in range(3):
|
|
for j in range(i+1, 3):
|
|
final_dist = np.linalg.norm(trajectories[i][-1] - trajectories[j][-1])
|
|
final_distances.append(final_dist)
|
|
|
|
avg_final_distance = np.mean(final_distances)
|
|
std_final_distance = np.std(final_distances)
|
|
|
|
# 计算质心移动距离
|
|
initial_com = solver.get_center_of_mass()
|
|
# 需要重新计算初始质心
|
|
total_mass = sum(p.mass for p in particles)
|
|
initial_com = np.zeros(3)
|
|
for p in particles:
|
|
initial_com += p.mass * p.position
|
|
initial_com /= total_mass
|
|
|
|
final_com = solver.get_center_of_mass()
|
|
com_movement = np.linalg.norm(final_com - initial_com)
|
|
|
|
results.append({
|
|
'seed': seed,
|
|
'avg_final_distance': avg_final_distance,
|
|
'std_final_distance': std_final_distance,
|
|
'com_movement': com_movement,
|
|
'energy_error': solver.get_conservation_errors()[2]
|
|
})
|
|
|
|
print(f" 平均最终距离: {avg_final_distance:.3f} AU")
|
|
print(f" 质心移动: {com_movement:.3f} AU")
|
|
print(f" 能量相对误差: {solver.get_conservation_errors()[2]:.6e}")
|
|
|
|
# 绘制比较图
|
|
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
|
|
|
|
seeds = [r['seed'] for r in results]
|
|
avg_distances = [r['avg_final_distance'] for r in results]
|
|
com_movements = [r['com_movement'] for r in results]
|
|
energy_errors = [r['energy_error'] for r in results]
|
|
|
|
# 平均最终距离
|
|
axes[0, 0].bar(range(n_simulations), avg_distances, color='skyblue', edgecolor='black')
|
|
axes[0, 0].set_xlabel('模拟编号', fontsize=12)
|
|
axes[0, 0].set_ylabel('平均最终距离 (AU)', fontsize=12)
|
|
axes[0, 0].set_title('质点间平均距离', fontsize=14, fontweight='bold')
|
|
axes[0, 0].set_xticks(range(n_simulations))
|
|
axes[0, 0].set_xticklabels([f"#{i+1}" for i in range(n_simulations)])
|
|
axes[0, 0].grid(True, alpha=0.3, axis='y')
|
|
|
|
# 质心移动
|
|
axes[0, 1].bar(range(n_simulations), com_movements, color='lightgreen', edgecolor='black')
|
|
axes[0, 1].set_xlabel('模拟编号', fontsize=12)
|
|
axes[0, 1].set_ylabel('质心移动距离 (AU)', fontsize=12)
|
|
axes[0, 1].set_title('系统质心移动', fontsize=14, fontweight='bold')
|
|
axes[0, 1].set_xticks(range(n_simulations))
|
|
axes[0, 1].set_xticklabels([f"#{i+1}" for i in range(n_simulations)])
|
|
axes[0, 1].grid(True, alpha=0.3, axis='y')
|
|
|
|
# 能量误差
|
|
axes[1, 0].bar(range(n_simulations), energy_errors, color='lightcoral', edgecolor='black')
|
|
axes[1, 0].set_xlabel('模拟编号', fontsize=12)
|
|
axes[1, 0].set_ylabel('能量相对误差', fontsize=12)
|
|
axes[1, 0].set_title('数值积分误差', fontsize=14, fontweight='bold')
|
|
axes[1, 0].set_xticks(range(n_simulations))
|
|
axes[1, 0].set_xticklabels([f"#{i+1}" for i in range(n_simulations)])
|
|
axes[1, 0].set_yscale('log')
|
|
axes[1, 0].grid(True, alpha=0.3, axis='y')
|
|
|
|
# 散点图：质心移动 vs 平均距离
|
|
axes[1, 1].scatter(avg_distances, com_movements, s=100, c=range(n_simulations),
|
|
cmap='viridis', edgecolors='black', alpha=0.8)
|
|
axes[1, 1].set_xlabel('平均最终距离 (AU)', fontsize=12)
|
|
axes[1, 1].set_ylabel('质心移动距离 (AU)', fontsize=12)
|
|
axes[1, 1].set_title('系统稳定性关系', fontsize=14, fontweight='bold')
|
|
|
|
# 添加标签
|
|
for i, (x, y) in enumerate(zip(avg_distances, com_movements)):
|
|
axes[1, 1].annotate(f"#{i+1}", (x, y), textcoords="offset points",
|
|
xytext=(0, 10), ha='center', fontsize=9)
|
|
|
|
axes[1, 1].grid(True, alpha=0.3)
|
|
|
|
plt.suptitle(f'{n_simulations}个随机三体系统模拟比较', fontsize=16, fontweight='bold')
|
|
plt.tight_layout()
|
|
|
|
output_file = "multiple_random_simulations.png"
|
|
plt.savefig(output_file, dpi=300, bbox_inches='tight')
|
|
print(f"\n比较图形已保存到: {output_file}")
|
|
plt.show()
|
|
|
|
return results
|
|
|
|
|
|
if __name__ == "__main__":
|
|
# 运行单个随机示例
|
|
print("运行单个随机三体系统示例...")
|
|
solver = run_random_example(seed=42, total_time=15.0)
|
|
|
|
# 运行多个随机模拟（可选）
|
|
# print("\n" + "="*60)
|
|
# print("运行多个随机模拟比较...")
|
|
# results = run_multiple_random_simulations(n_simulations=5, total_time=5.0) |