MATLAB牛顿摆仿真：事件检测与碰撞建模实践

简介牛顿摆仿真项目以Matlab实现面向物理仿真入门者及有一定Matlab基础的开发者可用于理解动量守恒、碰撞动力学以及数值仿真思路。压缩包内共1个文件为单个m脚本大小仅2KB结构轻量几乎不依赖额外工具箱适合快速运行与二次修改。源码经过实测校正可直接运行并通过物体碰撞、摆动动画及参数调整等环节清晰展示完整建模流程能为课程设计或兴趣实验提供参考。代码组织简洁注释与变量命名便于新手逐行阅读开发者也可以据此扩展为多球牛顿摆、加入能量损耗或阻尼效果进一步提升对物理模型与Matlab绘图机制的理解。项目规模虽小却涵盖了从参数设定、碰撞检测到动画输出的核心步骤很适合作为物理仿真学习的起点目前已有189人学习下载适合希望用较少代码体验物理仿真与动画展示的Matlab用户。1. 从“碰撞球”到事件驱动newton_pendulum 仿真的第一道门槛牛顿摆的仿真难点不在“算碰撞”而在碰撞瞬间的速度不连续。用固定步长直接递推球会在 1 毫秒内穿透一个微米级的 overlap能量漂移得让你怀疑自己在写随机游走。这道题的正确打开方式是把连续运动摆的慢动力学和瞬时碰撞离散事件分开处理用 ode45 的事件检测去“卡”住碰撞发生的精确时刻碰撞那一瞬再用恢复系数公式把速度改掉。这篇博文用 MATLAB 从零搭一个 newton_pendulum 仿真覆盖事件检测、碰撞更新、动画和能量验证适合正在做机器人碰撞、机械冲击或物理仿真的人也适合刚学 MATLAB 数值求解但想找一个完整案例的读者。跑通之后你会发现牛顿摆“一端进、另一端出”的直觉在数值上并不便宜。2. 碰撞建模与恢复系数动量守恒公式和 MATLAB 里的三个关键量在写代码之前先把牛顿摆简化成什么样想清楚。我一般会把每个球抽象成带质量的刚体球运动限制在水平一维方向碰撞只在相邻球之间发生且碰撞过程瞬时完成。这样处理丢掉了绳子的弯矩、球的形变和滚动摩擦但保留了牛顿摆的核心物理动量沿链传递、能量被恢复系数消耗。三个关键量分别是质量m、恢复系数e和碰撞阈值2R。前两个决定碰撞后的速度分配第三个决定事件函数什么时候“扣动扳机”。下面先定模型再给公式最后解释为什么事件检测比小步长递推更适合这类问题。2.1 先把牛顿摆拆成“链式碰撞”而不是“整体摆动”真正的牛顿摆每个球都是一个单摆球心位置应该是(x, y) (L sinθ, -L cosθ)广义坐标是摆角。如果直接用二维刚体摆建模碰撞检测要处理圆弧轨迹上的接触点复杂度会立刻上去。工程上常见的做法是假设摆长L远大于球的横向位移用水平位移x近似弧长于是每根绳子给球的恢复力写成a -(g / L) * sin(x / L); % 单摆线性项-g/L * x 是更粗糙的近似当初始摆幅小于 15° 时sin(x/L) ≈ x/L的误差不到 1%如果你把第一个球拉高到 30° 以上线性化会高估恢复力仿真出来的摆动周期会偏短这时就要保留sin项。代码里我直接用sin版本代价可以忽略物理上更稳。碰撞只发生在相邻两球之间这是“链式”二字的含义。第 1 个球撞第 2 个第 2 个撞第 3 个能量沿链传递除非有异常的多次反弹非相邻球不会直接接触。模型里我把第 25 个球的初始位置设为紧贴第 1 个球拉开一个水平位移让它在重力分量作用下回去撞链。2.2 恢复系数与动量守恒公式先写对速度更新一维对心碰撞有两个约束动量守恒一定成立机械能是否守恒取决于恢复系数e。设碰撞前两球速度为v1, v2碰撞后为v1, v2恢复系数定义为分离速度与接近速度之比e (v2 - v1) / (v1 - v2)联立动量守恒方程解出碰撞后的速度更新公式% 两球碰撞后的速度m1 m2 为质量e 为恢复系数 v1_new ((m1 - e * m2) * v1 (1 e) * m2 * v2) / (m1 m2); v2_new ((1 e) * m1 * v1 (m2 - e * m1) * v2) / (m1 m2);当e 1时公式退化为v1_new v2, v2_new v1也就是等质量理想弹性碰撞的速度交换这正是牛顿摆“中间球不动、两端球交换”的理论基础。当e 0时两球碰撞后速度相同完全非弹性能量损失最大。真实牛顿摆的不锈钢球e大概在 0.930.98 之间能量损失并不小这也是为什么它最终会停下来。参数表如下我后面的代码都按这组值跑参数符号建议取值说明球数N5经典配置奇数球更容易观察动量传递球质量m0.1 kg等质量时公式最漂亮非等质量也可跑球半径R0.02 m决定碰撞阈值2R摆长L0.3 m远大于位移满足水平近似条件恢复系数e0.981 为理想弹性0 为完全非弹性重力加速度g9.81 m/s²默认即可2.3 事件检测为什么比固定步长更适合碰撞瞬间如果你用dt 1e-4的固定步长递推碰撞发生的那一步两球中心距会从2R 0.001直接跳到2R - 0.002穿模之后要么把球强行推开导致能量突变要么让碰撞处理逻辑判断错方向。缩小步长只是把问题往后推计算量却涨了一个量级。MATLAB 的odeset(Events, ...)机制本质上是在积分过程中监视某个标量函数的零点。对牛顿摆我监视“所有相邻球间隙的最小值”value min(x(2:N) - x(1:N-1)) - 2R只要value 0说明没有任何一对球接触当它从正变负穿过零求解器就会用二分法把碰撞时刻卡出来而不是让球继续穿模。这样做的额外好处是如果多个相邻对同时接近min只挑出最早发生碰撞的那一对天然规避了“多事件同时触发”的歧义。事件函数里三个返回值的含义用一张表说清楚返回值作用本例设置value要监视的标量函数穿过 0 触发事件最小间隙减2Risterminal是否为终止事件1 表示停止积分恒为 1碰撞必须停下来处理direction只捕获方向-1表示只在递减时触发设为-1忽略分离过程direction -1在实际中很重要。球碰撞后分开时value也会从负变正穿过零如果不限定方向求解器会在碰撞后的第一个积分步又触发一次停止白白浪费计算。3. 用 ode45 Events 搭出可复现的 newton_pendulum 主循环理论部分够用了现在落到代码。整个仿真的主循环只有一个思路从当前时刻积分到下一个动画采样点或第一次碰撞时刻碰撞发生就更新速度然后继续积分。把这段跑通后面加动画、加能量统计都只是往这个骨架里填。3.1 主程序骨架状态向量、初始摆角和数值积分选项状态向量我用一个2N维行向量表示前N个是球的水平位移后N个是水平速度。这样的好处是 ode45 的y0不需要做矩阵变换事件函数里也能用切片直接取位移。% newton_pendulum_demo.m % 状态 z [x1..xN, v1..vN] N 5; m 0.1 * ones(1, N); % 每个球质量 0.1 kg R 0.02; % 球半径 L 0.3; % 摆长 g 9.81; % 重力加速度 e 0.98; % 恢复系数 % 初始位置第 2~5 个球紧贴第 1 个球向右拉开 0.2 m x0 zeros(1, N); for k 2:N x0(k) (k - 1) * (2*R 1e-6); % 微间隙避免事件函数在 t0 误触发 end x0(1) 0.2; v0 zeros(1, N); z [x0, v0]; t 0; T 5; % 仿真 5 秒 dt_plot 0.02; % 动画/记录间隔 t_next dt_plot; state_log []; t_log []; % 后面画图和能量统计用 while t T opts odeset(Events, (t,z) collide_events(t,z,N,R), ... RelTol, 1e-9, AbsTol, 1e-12); [t_seg, z_seg] ode45((t,z) pendulum_ode(t,z,N,g,L), [t t_next], z, opts); % 记录当前段的离散点4 倍降采样以控制内存 idx 1:max(1, floor(end/4)):end; state_log [state_log; z_seg(idx, :)]; t_log [t_log; t_seg(idx)]; z z_seg(end, :); t t_seg(end); % 如果这段里触发了碰撞事件函数会在 te 时刻停住 if t t_next z handle_collision(z, N, R, e, m); end t_next t_next dt_plot; end代码里有两个容易忽略的细节。x0(k) (k-1)*(2R 1e-6)给每个接触面留一个微米级缝隙是为了避免value在初始时刻正好为零否则 ode45 的事件检测可能第一时间把它当作“已经穿过零点”而误触。5.4行的floor(end/4)降采样是防止 5 秒仿真存出几十万行状态画图时卡到无法拖动。3.2 事件函数与碰撞检测的容差设置事件函数只需要 4 行function [value, isterminal, direction] collide_events(~, z, N, R) x z(1:N); value min(x(2:N) - x(1:N-1)) - 2*R; % 最小相邻间隙 isterminal 1; % 碰到就停 direction -1; % 只捕获接近过程 end注意我监视的是最小间隙不是每一对都监视。如果同时有第 1/2 对和第 3/4 对都在接近min会返回最先接触的那一对对应的间隙求解器只会停在最早碰撞时刻。碰撞处理函数里再统一扫描所有可能接触的对这样事件函数不用关心“谁先撞”只需回答“什么时候撞”。RelTol和AbsTol这里给到1e-9和1e-12原因有两个碰撞时刻的定位精度直接影响后续能量统计默认的RelTol1e-3会让碰撞时刻的误差到毫秒量级在 5 秒仿真里累计出来的相位移相当明显另外碰撞后速度的跳变会让积分器在事件点附近反复试探容差太松会更快触发“检测到奇点”的假错误。把 AbsTol 设成1e-12就够了别写成1e-100MATLAB 里浮点数当然能表示它但这个量级会逼迫步长收缩到不可理喻的程度积分器直接报废。3.3 碰撞后速度更新把非连续跳变合回求解流程事件函数负责“卡住”碰撞时刻真正的物理更新在handle_collision里function z handle_collision(z, N, R, e, m) x z(1:N); v z(N1:2*N); for k 1:N-1 dv v(k) - v(k1); % 接近速度正数表示正在接近 gap x(k1) - x(k); % 当前中心距 if gap 2*R dv 0 % 位置修正把 overlap 按质量比例分摊回两个球 overlap 2*R - gap; x(k) x(k) - overlap * m(k1) / (m(k) m(k1)); x(k1) x(k1) overlap * m(k) / (m(k) m(k1)); % 速度更新一维对心碰撞公式 v1 v(k); v2 v(k1); v(k) ((m(k) - e*m(k1))*v1 (1e)*m(k1)*v2) / (m(k)m(k1)); v(k1) ((1e)*m(k)*v1 (m(k1) - e*m(k))*v2) / (m(k)m(k1)); end end z [x, v]; end位置修正这一步很多人会漏。理论上事件检测触发时两球已经轻微 overlap虽然只有微米量级但不推开的话下一轮积分里事件函数会立刻再次触发造成死循环。修正方式是按质量比分摊 overlap质量大的球被推得更少符合质心不变的物理直觉。提示handle_collision里的循环是顺序扫描相邻对。如果五个球被挤压成一团这种逐对处理严格来说有方向偏差对牛顿摆演示场景这个偏差肉眼不可见但如果你拿这套代码做科研建议换成多体冲量矩阵求解或对碰撞对做 23 轮 Jacobi 迭代。这也是很多物理引擎里速度级 LCP 松弛的雏形。4. 动画驱动与能量曲线让 MATLAB 把仿真“画”出来状态算出来只是第一半牛顿摆的价值一半在“看”。MATLAB 里把球和绳子画出来只需要plot和rectangle两条核心指令但要让动画不卡、导出不糊还有几个参数要调。4.1 绘制牛顿摆绳子、球与摆线的三条绘图指令每根绳子的悬挂点固定在横梁上球心位置由状态决定。注意悬挂点的横坐标和球静止时的横坐标一致而不是都在原点function draw_pendulum(t, z, N, R, L) x z(1:N); x_pin (0:N-1) * (2*R 1e-6); % 悬挂点横坐标 y -sqrt(max(L^2 - (x - x_pin).^2, 0)); % 摆线末端高度 cla; hold on; for k 1:N % 画绳子从悬挂点到球心 plot([x_pin(k) x(k)], [0 y(k)], k-, LineWidth, 1.2); % 画球用 rectangle 加圆角 rectangle(Position, [x(k)-R, y(k)-R, 2*R, 2*R], ... Curvature, [1 1], FaceColor, [0.3 0.7 1.0], ... EdgeColor, k); end axis equal; xlim([-L L]); ylim([-L 0.05]); title(sprintf(newton_pendulum t %.2f s, t)); drawnow limitrate; % 限制刷新率动画不卡 enddrawnow limitrate是动画的关键。它允许 MATLAB 丢帧保证动画节奏跟实时时间接近而不是为了渲染每一帧把仿真拖慢 10 倍。如果是离线导出视频则可以不要这个逐帧全渲染。rectangle加Curvature, [1 1]是 MATLAB 里画圆的惯用做法比viscircles更适合动画场景因为它可以直接接受FaceColor填充不需要额外拼接图像对象。球的数量不多时这种方式足够不用考虑scatter或patch的性能优化。4.2 能量曲线与守恒量计算画图前先把物理量算对动画只是视觉能量曲线才是验证物理正确性的工具。仿真过程中动能和势能要分别累积function [Ek, Ep] compute_energy(z, N, m, g, R, L) x z(1:N); v z(N1:2*N); x_pin (0:N-1) * (2*R 1e-6); y -sqrt(max(L^2 - (x - x_pin).^2, 0)); Ek 0.5 * sum(m .* v.^2); Ep sum(m .* g .* y); % 以横梁为 0 势能面 end注意势能这里用的是y也就是绳子的垂直分量高度。如果你在状态向量里只存了水平位移千万别直接用-x当高度那会把摆的弧线运动错误地当成竖直下落。用-sqrt(L^2 - x^2)才能把摆长约束考虑进来。有碰撞时总能量是阶梯状下降的每一步下降量可以用公式预估ΔE (1 - e^2) / 2 * μ * (v_rel)^2其中μ m1*m2 / (m1 m2)是约化质量v_rel v1 - v2是碰撞前接近速度。我习惯把这个式子写在能量图旁边做对照如果仿真里能量下降量和公式差超过 2%基本可以断定事件检测少抓了一次碰撞。4.3 导出视频与状态数据仿真结果怎么沉淀动画看完要能带走。VideoWriter 导出 MP4 的标准流程如下vp VideoWriter(newton_pendulum.mp4, MPEG-4); vp.FrameRate 60; open(vp); for i 1:size(state_log, 1) draw_pendulum(t_log(i), state_log(i, :), N, R, L); writeVideo(vp, getframe(gcf)); end close(vp); % 同时把状态存成 CSV方便后续分析 writematrix([t_log, state_log], pendulum_states.csv);FrameRate可以按需调整60 fps 播放 5 秒仿真会得到 300 帧文件体积和流畅度比较平衡。存 CSV 是因为后续你可能不只想在 MATLAB 里看把 CSV 读进 Python 或别的工具做 FFT 分析、统计碰撞间隔都比在 MATLAB 里维护一个大struct方便。5. 参数扫描与结果验证把恢复系数、球数、初始摆角调到物理自洽仿真能跑、动画能看之后真正的工程问题是这组参数下的结果是否符合物理直觉把e从 1.0 调到 0.9末端球的速度会掉多少把球数从 5 改成 7动量传递还会那么干净吗这一章用表格和验证脚本回答这些问题。5.1 参数影响对照表与物理直觉我做了三组仿真每组只改一个参数结果如下参数改动现象物理解释e 1.0第 1 球停住第 5 球几乎以原速弹出中间球不动理想弹性碰撞速度完全交换e 0.98第 5 球弹出速度约为入射的 96%第 1 球有微小回弹能量损失导致分离速度略低球的微小反向速度来自残余动量e 0.8第 5 球速度降到约 60%中间球也明显运动每次碰撞损失约 20% 相对动能链式传递迅速衰减N 5 → 7末端球弹出速度不变但延迟约 2 倍碰撞间隔等质量等弹性时动量传递速度与球数无关初始摆幅从 0.2 m 加到 0.35 m碰撞间隔变短周期有明显非线性大角度下sin(x/L)偏离线性恢复力不再对称这个表格说明一个容易被忽视的点恢复系数对动量守恒没有影响但对能量传递的衰减非常敏感。e 0.98看着只差 2%经过五次碰撞后末端球的动能只剩理论值的0.98^5 ≈ 0.90也就是 90%。这也是为什么真实牛顿摆不能无限摆下去每一次碰撞都在消耗能量。5.2 守恒量验证脚本动量与能量的数值误差怎么看判断仿真是否可信我只看两个指标动量守恒残差和能量守恒曲线。动量在理想条件下应该严格守恒因为内力不改变总动量能量的总下降量则要和恢复系数公式对得上。% 读回状态验证守恒量 data readmatrix(pendulum_states.csv); t data(:, 1); z data(:, 2:end); N 5; m 0.1 * ones(1, N); p_total zeros(size(t)); E_total zeros(size(t)); for i 1:length(t) v z(i, N1:2*N); p_total(i) sum(m .* v); [~, Ep] compute_energy(z(i, :), N, m, g, R, L); Ek 0.5 * sum(m .* v.^2); E_total(i) Ek Ep; end fprintf(动量残差: %.3e\n, max(abs(p_total - p_total(1)))); plot(t, E_total, LineWidth, 1.5);动量残差在RelTol1e-9下通常能到1e-14量级如果残差到了1e-10以上先查事件检测有没有漏碰撞能量曲线则要注意它是不是只减不增出现“碰撞后能量反而变大”基本说明handle_collision里的接近速度方向判断写反了。5.3 数据后处理把碰撞时刻“卡”得更准的两个小技巧仿真数据拿到手后有两个高频需求分析末端球的速度频谱以及检查碰撞时刻序列。末端球速度做 FFT 可以直接用 MATLAB 自带函数但横轴时间点可能多到看不清这时要稀疏化刻度fs 1 / mean(diff(t)); v5 z(:, end); % 最后一个球的速度 V5 fft(v5 - mean(v5)); f_axis (0:length(V5)-1) * fs / length(V5); plot(f_axis(1:floor(end/2)), abs(V5(1:floor(end/2)))); xlim([0 10]); ylabel(幅值); xlabel(频率 (Hz)); xticks(0:0.5:10); % 刻度稀疏化避免糊成一团谱线里你会看到两个特征频率一个是摆的固有频率sqrt(g/L)/(2π) ≈ 0.91 Hz另一个是碰撞产生的宽带脉冲尖峰频率和碰撞间隔成正比。用findpeaks把这个尖峰位置找出来就能推算出每次碰撞的精确时刻反推到事件检测里验证一下——如果找出的碰撞次数小于你代码里实际处理的次数说明有碰撞没被卡住。最后一个实用技巧把handle_collision里的碰撞次数用一个全局计数器累积仿真结束打印出来。这个数字和findpeaks检出的碰撞时刻数对得上你的事件检测就是准的对不上优先检查direction和value初始符号这两个是事件检测里最常埋坑的地方。本文还有配套的精品资源点击获取