C语言完美弹性碰撞运动学模拟：y位置无法回归初始点

问题背景 #

学习C语言期间，以运动学为项目实现完美弹性碰撞模拟，但代码运行后，物体的y位置（d_1）无法回到初始的100，最高只能到97.29左右，且误差随时间累积。原代码及输出问题如下：

原代码 #

#include <stdio.h>
#include <stdbool.h>
#include <string.h>

int main(void)
{
    double v_1;
    double v_0;
    double Δv;
    double Δt;
    double t_0;
    double t_1;
    double d_0;
    double d_1;
    double Δd;
    double a;

    int elapsed_time = 0;
    t_0 = 0;
    t_1 = 0;
    v_0 = 0;
    v_1 = 0;
    a = -9;
    d_0 = 100; // starting displacement
    d_1 = 100; // current displacement, or the y position

    while (elapsed_time < 1000){

        Δt = t_1 - t_0;
        Δd = d_1 - d_0;
        Δv = v_1 - v_0;

        v_1 = v_0 + a * Δt;
        d_1 = d_0 + v_1 * Δt;

        if (d_1 <= 0)
        {
            d_1 = 0;
            v_1 = -v_1;
        }
        d_0 = d_1;
        v_0 = v_1;
        t_0 = t_1;
        t_1 += 0.1;
        elapsed_time += 1;
        printf("d_1 = %lf, v_1 = %lf, v_0 = %lf, a * Δt = %lf\n", d_1, v_1, v_0, a * Δt);

    }

    return 0;
} 

关键输出问题 #

物体下落碰撞地面后，反弹的最高位置仅为97.29，而非初始的100，且后续反弹高度持续降低。

问题成因分析 #

1. 积分方法的顺序错误
   代码中先更新速度，再用新速度计算位移，这偏离了标准显式欧拉积分的逻辑，会引入一阶近似误差。对于无能量损失的弹性碰撞这类保守系统，误差会不断累积，导致能量“泄漏”。

2. 碰撞处理的近似误差
   当检测到d_1 <= 0时直接将位置设为0并反转速度，此时物体已经在时间步长Δt内穿透了地面，这一步截断会损失部分动能，直接导致反弹高度下降。

3. 浮点数精度的累积影响
   浮点数舍入误差存在，但并非主要问题，核心误差来自积分和碰撞处理的近似逻辑。

修复方案 #

方案1：修正积分顺序 + 精确碰撞处理 #

核心调整：

• 改用标准显式欧拉积分：先用当前速度计算位移，再更新速度。
• 碰撞时计算精确的碰撞时间，修正位置和速度，避免能量损失。

修复后的代码

#include <stdio.h>

int main(void)
{
    double v = 0.0;
    double d = 100.0;
    const double a = -9.0;
    const double dt = 0.1;
    int elapsed_time = 0;

    while (elapsed_time < 1000) {
        // 标准显式欧拉：先算位移，再算速度
        double d_new = d + v * dt;
        double v_new = v + a * dt;

        // 精确处理碰撞
        if (d_new <= 0.0) {
            // 计算到达地面的精确时间比例
            double t_collide = (-d) / v;
            // 剩余时间
            double t_remaining = dt - t_collide;

            // 碰撞瞬间速度反转，再计算剩余时间的运动
            double v_collide = v + a * t_collide;
            v_new = -v_collide + a * t_remaining;
            d_new = 0.0 + (-v_collide) * t_remaining;
        }

        // 更新状态
        d = d_new;
        v = v_new;
        elapsed_time += 1;
        printf("d = %.6lf, v = %.6lf\n", d, v);
    }

    return 0;
}

方案2：改用Verlet积分（推荐） #

Verlet积分专为保守系统设计，能更好地保持能量守恒，误差更小，且无需显式存储速度，可通过前后位置推导。

Verlet版本代码

#include <stdio.h>

int main(void)
{
    double d_prev = 100.0; // 上一时刻位置
    double d_curr = 100.0; // 当前位置
    const double a = -9.0;
    const double dt = 0.1;
    int elapsed_time = 0;

    while (elapsed_time < 1000) {
        // Verlet积分核心公式
        double d_next = 2.0 * d_curr - d_prev + a * dt * dt;

        // 碰撞处理：利用对称性修正穿透问题
        if (d_next <= 0.0) {
            double penetration = -d_next;
            d_next = penetration;
            d_prev = d_curr + penetration;
        } else {
            d_prev = d_curr;
        }

        d_curr = d_next;
        elapsed_time += 1;
        printf("d = %.6lf\n", d_curr);
    }

    return 0;
}

修复效果 #

使用上述修复后的代码，物体反弹后会精确回到初始位置100，且后续反弹高度保持稳定，完全符合弹性碰撞的能量守恒特性。

内容的提问来源于stack exchange，提问作者Curulian