模拟两个小球的弹性碰撞，展示动量守恒和能量守恒的过程。

碰撞实验
模拟两个小球的弹性碰撞，展示动量守恒和能量守恒的过程。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib.patches import Circle

# 设置图片清晰度
plt.rcParams['figure.dpi'] = 300

# 设置中文字体
plt.rcParams["font.family"] = ["SimHei", "WenQuanYi Micro Hei", "Heiti TC"]

# 调整字体大小
plt.rcParams.update({
    'font.size': 3,  # 基础字体大小
    'axes.titlesize': 3,  # 标题字体大小
    'axes.labelsize': 3,  # 坐标轴标签字体大小
    'xtick.labelsize': 3,  # x轴刻度字体大小
    'ytick.labelsize': 3,  # y轴刻度字体大小
    'legend.fontsize': 3,  # 图例字体大小
})

# 创建画布和子图
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(8, 12))
fig.subplots_adjust(hspace=0.4)

# 碰撞模拟区域
ax1.set_xlim(0, 10)
ax1.set_ylim(0, 5)
ax1.set_aspect('equal')
ax1.set_title('小球弹性碰撞模拟')
ax1.set_xlabel('X 位置')
ax1.set_ylabel('Y 位置')
ax1.grid(True)

# 动量守恒图
ax2.set_xlim(0, None)
ax2.set_ylim(0, None)
ax2.set_title('动量变化（有外力）')
ax2.set_xlabel('时间')
ax2.set_ylabel('动量')
ax2.grid(True)

# 能量守恒图
ax3.set_xlim(0, None)
ax3.set_ylim(0, None)
ax3.set_title('能量变化（有外力）')
ax3.set_xlabel('时间')
ax3.set_ylabel('能量')
ax3.grid(True)

# 物理参数
g = 0.5  # 重力加速度
mu = 0.1  # 摩擦系数

# 小球参数
m1 = 1.0  # 小球1质量
m2 = 1.0  # 小球2质量
r1 = 0.5  # 小球1半径
r2 = 0.5  # 小球2半径

# 初始位置和速度
x1, y1 = 2.0, 2.5
x2, y2 = 8.0, 2.5
vx1, vy1 = 1.0, 0.0
vx2, vy2 = -1.0, 0.0

# 创建小球
ball1 = Circle((x1, y1), r1, color='red')
ball2 = Circle((x2, y2), r2, color='blue')
ax1.add_patch(ball1)
ax1.add_patch(ball2)

# 初始化动量和能量图表
time_data = []
p1_data, p2_data, ptot_data = [], [], []
ke1_data, ke2_data, ketot_data = [], [], []
pe1_data, pe2_data, petot_data = [], [], []  # 势能数据
me1_data, me2_data, metot_data = [], [], []  # 机械能数据
p1_line, = ax2.plot([], [], 'r-', label='小球1动量')
p2_line, = ax2.plot([], [], 'b-', label='小球2动量')
ptot_line, = ax2.plot([], [], 'g-', label='总动量')
ke1_line, = ax3.plot([], [], 'r-', label='小球1动能')
ke2_line, = ax3.plot([], [], 'b-', label='小球2动能')
pe1_line, = ax3.plot([], [], 'r--', label='小球1势能')
pe2_line, = ax3.plot([], [], 'b--', label='小球2势能')
me1_line, = ax3.plot([], [], 'r-.', label='小球1机械能')
me2_line, = ax3.plot([], [], 'b-.', label='小球2机械能')
metot_line, = ax3.plot([], [], 'g-.', label='总机械能')
ax2.legend()
ax3.legend()


# 弹性碰撞处理函数
def handle_collision():
    global x1, y1, x2, y2, vx1, vy1, vx2, vy2

    # 计算两球之间的距离
    dx = x2 - x1
    dy = y2 - y1
    d = np.sqrt(dx ** 2 + dy ** 2)

    # 如果两球相撞
    if d < r1 + r2:
        # 计算碰撞法线方向
        nx = dx / d
        ny = dy / d

        # 计算相对速度
        dvx = vx2 - vx1
        dvy = vy2 - vy1

        # 计算相对速度在法线方向的分量
        vn = dvx * nx + dvy * ny

        # 如果两球正在靠近，进行碰撞处理
        if vn < 0:
            # 计算碰撞后的速度变化
            j = -(1 + 1.0) * vn / (1 / m1 + 1 / m2)  # 弹性系数为1（完全弹性碰撞）
            jx = j * nx
            jy = j * ny

            # 更新速度
            vx1_new = vx1 - jx / m1
            vy1_new = vy1 - jy / m1
            vx2_new = vx2 + jx / m2
            vy2_new = vy2 + jy / m2

            # 分离两球，防止粘连
            overlap = 0.5 * (r1 + r2 - d + 0.01)
            x1 -= overlap * nx
            y1 -= overlap * ny
            x2 += overlap * nx
            y2 += overlap * ny

            return vx1_new, vy1_new, vx2_new, vy2_new

    return vx1, vy1, vx2, vy2


# 更新函数，用于动画每一帧的绘制
def update(frame):
    global x1, y1, x2, y2, vx1, vy1, vx2, vy2

    # 时间步长
    dt = 0.05

    # 应用重力
    vy1 -= g * dt
    vy2 -= g * dt

    # 应用摩擦力（与速度方向相反）
    speed1 = np.sqrt(vx1 ** 2 + vy1 ** 2)
    if speed1 > 0:
        friction_force1 = -mu * g * m1
        acc_x1 = friction_force1 * vx1 / speed1 / m1  # 避免与子图ax1冲突
        acc_y1 = friction_force1 * vy1 / speed1 / m1
        vx1 += acc_x1 * dt
        vy1 += acc_y1 * dt

    speed2 = np.sqrt(vx2 ** 2 + vy2 ** 2)
    if speed2 > 0:
        friction_force2 = -mu * g * m2
        acc_x2 = friction_force2 * vx2 / speed2 / m2  # 避免与子图ax2冲突
        acc_y2 = friction_force2 * vy2 / speed2 / m2
        vx2 += acc_x2 * dt
        vy2 += acc_y2 * dt

    # 更新位置
    x1 += vx1 * dt
    y1 += vy1 * dt
    x2 += vx2 * dt
    y2 += vy2 * dt

    # 边界条件（添加能量损失）
    if x1 - r1 < 0:
        x1 = r1
        vx1 = -vx1 * 0.9  # 能量损失
    elif x1 + r1 > 10:
        x1 = 10 - r1
        vx1 = -vx1 * 0.9  # 能量损失

    if y1 - r1 < 0:
        y1 = r1
        vy1 = -vy1 * 0.9  # 能量损失
    elif y1 + r1 > 5:
        y1 = 5 - r1
        vy1 = -vy1 * 0.9  # 能量损失

    if x2 - r2 < 0:
        x2 = r2
        vx2 = -vx2 * 0.9  # 能量损失
    elif x2 + r2 > 10:
        x2 = 10 - r2
        vx2 = -vx2 * 0.9  # 能量损失

    if y2 - r2 < 0:
        y2 = r2
        vy2 = -vy2 * 0.9  # 能量损失
    elif y2 + r2 > 5:
        y2 = 5 - r2
        vy2 = -vy2 * 0.9  # 能量损失

    # 处理碰撞
    vx1, vy1, vx2, vy2 = handle_collision()

    # 更新小球位置
    ball1.center = (x1, y1)
    ball2.center = (x2, y2)

    # 计算动量和能量
    p1 = m1 * np.sqrt(vx1 ** 2 + vy1 ** 2)
    p2 = m2 * np.sqrt(vx2 ** 2 + vy2 ** 2)
    ptot = p1 + p2
    ke1 = 0.5 * m1 * (vx1 ** 2 + vy1 ** 2)
    ke2 = 0.5 * m2 * (vx2 ** 2 + vy2 ** 2)
    pe1 = m1 * g * (y1 - r1)  # 势能参考点为地面
    pe2 = m2 * g * (y2 - r2)
    me1 = ke1 + pe1  # 机械能 = 动能 + 势能
    me2 = ke2 + pe2
    ketot = ke1 + ke2
    petot = pe1 + pe2
    metot = me1 + me2

    # 更新数据
    t = frame * dt
    time_data.append(t)
    p1_data.append(p1)
    p2_data.append(p2)
    ptot_data.append(ptot)
    ke1_data.append(ke1)
    ke2_data.append(ke2)
    pe1_data.append(pe1)
    pe2_data.append(pe2)
    me1_data.append(me1)
    me2_data.append(me2)
    metot_data.append(metot)

    # 更新图表
    p1_line.set_data(time_data, p1_data)
    p2_line.set_data(time_data, p2_data)
    ptot_line.set_data(time_data, ptot_data)
    ke1_line.set_data(time_data, ke1_data)
    ke2_line.set_data(time_data, ke2_data)
    pe1_line.set_data(time_data, pe1_data)
    pe2_line.set_data(time_data, pe2_data)
    me1_line.set_data(time_data, me1_data)
    me2_line.set_data(time_data, me2_data)
    metot_line.set_data(time_data, metot_data)

    # 自动调整y轴范围
    if frame > 5:
        ax2.set_ylim(0, max(p1_data + p2_data + ptot_data) * 1.1)
        ax3.set_ylim(0, max(ke1_data + ke2_data + pe1_data + pe2_data + me1_data + me2_data + metot_data) * 1.1)
        ax2.set_xlim(0, t)
        ax3.set_xlim(0, t)

    return ball1, ball2, p1_line, p2_line, ptot_line, ke1_line, ke2_line, pe1_line, pe2_line, me1_line, me2_line, metot_line


# 创建动画
ani = FuncAnimation(fig, update, frames=range(300), interval=30, blit=True)

# 显示动画
plt.tight_layout()
plt.show()

# 如果需要保存动画，取消下面一行的注释
# ani.save('elastic_collision.gif', writer='pillow', fps=20)