简介压缩包chp5_ex2.zip是配合丁承先向量式结构力学内容编写的质点法练习二面向正在学习计算结构力学、MATLAB编程的本科生或工程师目的是通过一个可直接运行的源代码示例理解将连续体离散为质量点、按相互作用力建立运动方程并数值求解的完整流程。质点法的核心优势在于无需复杂网格生成对非均匀、复杂几何结构尤其有效因此这一练习也适合作为从有限元思维过渡到无网格/粒子方法的入门参考。包内仅含1个m文件压缩后约2KB文件虽小但结构完整覆盖质点位置与刚度定义、弹性/重力等相互作用力计算、牛顿第二定律组建代数方程组、欧拉或龙格-库塔时间积分以及边界条件处理和结果输出等多个关键环节。已有116人学习浏览对于想快速上手质点法的读者来说该示例能提供从理论到代码的直观映射运行和分析时还可借鉴到适合结构力学场景的矩阵构建与求解思路并能进一步迁移到流体、地震工程等更广领域。1. 拿到 chp5_ex2.zip质点法练习的第二部分要先想清楚什么chp5_ex2.zip 这种命名方式在图形学与物理仿真的课程包里很常见chp5 表示第 5 章ex2 表示第 2 个练习后面跟的“质点法2”说明这不是第一次接触粒子类模拟。压缩包里边通常是源码、初始数据和题目说明你拿到手的第一反应不该是双击运行而是先把“这套代码在算哪种质点法、要你改哪一块”搞清楚。按我的分类质点法是一个宽泛说法既包括游戏引擎里的粒子系统也包括质点弹簧、光滑粒子流体动力学SPH和物质点法MPM它们的共同点是用一堆离散质点代替连续的物体或场然后用牛顿第二定律驱动每个点运动。这篇文章就顺着 chp5_ex2 这个练习包的常见结构把解包、最小实现、积分器选择、碰撞与并行化、以及如何验证结果一条线讲到位。2. 把 chp5_ex2.zip 摊开目录结构、依赖检查与质点法最小模型2.1 先看压缩包里有什么再决定用哪条命令把它解开我的习惯是解压之前先列清单不解压也能知道包里是不是你要的东西。用 unzip 的 -l 选项只打印压缩包内的文件列表不落盘unzip -l chp5_ex2.zip输出会列出每个文件的权限、大小、日期和完整路径。看到 src/、data/、config/ 这类目录结构基本就能判断是一个带工程骨架的练习包而不是几页散装代码。确认无误后解压到独立目录避免把一堆文件直接撒在下载目录里mkdir -p chp5_ex2 unzip chp5_ex2.zip -d chp5_ex2-d 参数指定目标目录名后续编译或运行都在 chp5_ex2 下进行不会污染其他地方。如果这一行报了 error read zip archive 之类的错先跑一遍unzip -t chp5_ex2.zip做完整性测试多半是文件在传输中被截断这种情况重下远比修 zip 省时间。我整理了一份课程练习包常见的布局对照路径典型内容练习时需要动的地方src/质点类、力模型、积分器实现题目要求补全或改写的部分data/初始坐标、边界盒、材质参数换测试场景时在这里改config/JSON、YAML 或 INI 参数每轮实验都会调README.md题目说明、编译命令、运行要求动手之前先读如果你是从 GitHub 上下载的 zip 包目录里通常还会多一层仓库根目录解压后 cd 进去即可其余逻辑和这里一致。把目录结构看过一遍再去看 README 里要求的依赖版本一般练习包的代码不会依赖太新的编译器或 SDK但 Python 类的练习经常要 numpy、scipy 中的固定小版本有版本偏差时先按项目要求的版本建虚拟环境省得后面跑出莫名其妙的报错。提示unzip -t 只校验压缩包内部数据的完整性不对每个文件做业务层面的校验。文本文件解出来是乱码通常是编码问题不是压缩包损坏。2.2 质点法的力学模型不管文件多少核心就 F ma质点法的出发点是拉格朗日视角不追踪网格而是把质量离散到一组点上。每个质点的状态由位置 x 和速度 v 描述运动方程可以写成两个一阶常微分方程组dx/dt vdv/dt F(x, v) / m这里的 F 是所有外力与内力之和。课程练习的第一部分通常只做自由粒子在重力场中的运动到了“质点法2”这个位置一般会加入粒子与粒子之间的相互作用也就是在 F 里补上弹簧力、阻尼力或者压力梯度项。写实现时我不建议把每种力单独写成一段顺次执行的代码而是先定义一个累加器把每个力模型的返回值加进去全部累加完成后再除以质量得到加速度。这样做的好处是新增一种力只需要补一个函数不需要改动积分流程。常见模型的力表达式集中列成一张表方便对照力的类型表达式作用重力F m * g提供外力驱动方向恒定弹簧胡克F -k * (r线性阻尼F -c * v_rel消耗相对动能抑制持续震荡边界碰撞法向速度反向或位置修正防止质点穿出计算域计算顺序也有讲究先求出同类质点与邻近质点的相互作用再叠加重力和边界力最后统一用积分器推进。练习包里常见的错误是把重力在每次迭代里重复乘以 dt或者把阻尼因子写成 dt 的幂次这两种写法都会让结果偏离物理直觉调参时尤其难找原因。2.3 最小可运行示例用纯 Python 把一次自由落体模拟跑通在没有题目代码的情况下我会先用一段最小实现验证环境再往里面加约束。下面的脚本只包含一个质点和重力外加 200 步显式欧拉积分# demo_particle.py import numpy as np class Particle: def __init__(self, x, v, mass1.0): self.x np.array(x, dtypefloat) # 初始位置 [x, y] self.v np.array(v, dtypefloat) # 初始速度 [vx, vy] self.mass mass # 质量默认 1.0 def step(self, dt, gravity9.8): # 只算重力方向朝下 force np.array([0.0, -self.mass * gravity]) acc force / self.mass self.v self.v acc * dt # 先更新速度 self.x self.x self.v * dt # 再更新位置 p Particle([0.0, 10.0], [1.0, 0.0], mass1.0) for _ in range(200): p.step(0.01) print(final x , p.x, final v , p.v)逻辑上这版实现遵循先速度后位置的显式欧拉更新顺序加速度恒定为 9.8 向下速度增量等于加速度乘以 dt位置增量等于更新后的速度乘以 dt。参数上 dt0.01 是迭代步长200 步对应 2 秒的物理时间最终位置的 y 坐标大约在 4.6 附近。想要验证解析解把 dt 改小到 0.002 并跑 1000 步数值结果会向 y 10 - 0.5 * 9.8 * 2² 靠近。这段最小代码跑通了再往 src 里加弹簧约束心里才有底。3. 时间积分与弹簧约束质点法参数调优的三个关键点3.1 显式欧拉不是不能用但你要知道它误差从哪里来显式欧拉的每一步都沿着当前位置的切线走直线切线方向会随着曲率变化而偏移累加几百步之后轨迹和真实解之间的偏差就会肉眼可见。在弹簧或轨道这类强弯曲运动中显式欧拉会让系统看起来“越蹦越高”本质是离散化在向系统内持续注入能量。质点法练习包里几乎都会给你一个积分器替换任务从欧拉改成 Verlet。位置 Verlet 不需要显式保存速度它用上一帧和当前帧的位置推算下一帧位置加速度只在这里出现def euler_update(x, v, accel, dt): v_new v accel * dt x_new x v_new * dt return x_new, v_new def verlet_update(x, x_prev, accel, dt): # x_prev 是上一帧位置accel 是当前受力计算的加速度 x_next 2.0 * x - x_prev accel * (dt * dt) return x_next两种函数对比Verlet 少了一次速度存储代码更短但在碰撞处理里需要恢复速度时反而要多一次差分运算。实际工程里我一般这么选粒子系统用位置 Verlet因为它在约束投影场景下实现简单需要在碰撞瞬间获得准确速度的刚体粒子模拟用 RK2 或 half-step Verlet。half-step Verlet 的做法是先用半个步长更新速度再用新速度更新一整步位置最后再用半个步长更新速度速度的获取代价只增加一次加法。参数上dt 直接决定两类误差截断误差随 dt 变小而变小但浮点舍入误差随迭代次数增加而累积。对弹簧系统先跑一个保守的 dt 1/240也就是每帧最多 4 个子步观察最大位移是否随时间线性增长如果位置曲线在第 2000 步附近出现折返就说明积分器在共振点上发散了。3.2 弹簧约束与阻尼把一堆独立的点做成一块布质点法练习到了第二部分最常见的任务是把自由粒子用弹簧连起来模拟布料或弹性绳。弹簧力方向沿两质点连线大小按当前伸长量与静止长度的差计算再叠加沿连线方向的阻尼速度# spring_force.py import numpy as np def spring_force(a, b, rest_len, k, c): a, b 为两个质点对象需要 x 和 v 属性 rest_len 是弹簧静止长度k 是刚度c 是阻尼系数 r b.x - a.x dist np.linalg.norm(r) if dist 1e-9: # 防止除零 return np.zeros(2), np.zeros(2) direction r / dist # 从 a 指向 b 的单位向量 v_rel b.v - a.v force_mag -k * (dist - rest_len) - c * np.dot(v_rel, direction) fa force_mag * direction # a 受到的力 fb -fa # 作用力与反作用力 return fa, fb参数上 k 是拉伸刚度数值越大抵抗形变越强c 是阻尼系数控制相对运动被削弱的程度。c 相对 k 过小时弹簧会沿着轴向反复震荡表现为布料上出现持续抖动c 相对 k 过大时相邻粒子黏在一起布料像面团。以粒子质量 m1 为单位我通常会从 k100、c0.5 起步再把 k 每档降低一倍看效果。调参时最容易忽略的是 rest_len 与初始网格边长的一致性。网格生成时如果相邻粒子距离是 0.3而 rest_len 被设成 1.0弹簧从第一帧就在做大行程收缩系统会在收敛前产生剧烈的初始振动。先把 rest_len 打印出来和初始网格边长对比再去调 k、c顺序不要反。下面是组合参数表按上面这段代码直接可用参数起步值调大后的效果常见的坑rest_len等于初始网格边长布料被拉伸或压缩与初始距离不匹配会先振动k100更硬约束更强k 过大配大 dt 会爆c0.5振动衰减更快c 过大让运动迟钝dt1/240单步更准太小则总耗时线性上升3.3 时间步长怎么压先看波形再看数值接手一个几百行的质点法练习代码我的调参顺序是这样的先不碰任何物理参数只改 dt把 dt 1/240 和 dt 1/480 各跑一遍对比同一时刻的质点位置。如果两条轨迹几乎重合说明当前 dt 已经足够小可以回退到 1/240 省时间如果差距超过 5%说明约束或碰撞引入了大量高频分量继续细化 dt 才是正路而不是转而加大阻尼去掩盖问题。判断发散不看单帧截图要看趋势。每 100 步打印一次最大速度的模长出现连续蹿升或者出现 inf/NaN立刻把 dt 减半重跑。注意这里说的减半不是把迭代次数减半而是在相同物理时长内加密子步数所以总耗时不变。还有一种现象是代码本身没毛病但报错出在文件读取而不是积分器上例如配置里写了data/particles.csv实际路径是config/../data/particles.csv。这类问题在练习包里出现频率很高解包后先把 data 目录相对路径和 README 对照一遍比我刚才提到的 error read zip archive 更隐蔽因为它是运行时才出现而且报错内容不会提示路径深度不对。注意打印最大速度时不要用自带的 print 直接刷屏每 100 步输出一行就够了否则日志文件会拖慢整体模拟。4. 大规模质点模拟空间哈希、碰撞修正与并行改写边界4.1 空间网格哈希把 O(n²) 的碰撞查询压下来几千个质点彼此两两判断距离一次迭代就要上千万次运算浏览器里的 JS 演示还能勉强跑桌面程序也扛不住上万的 n。最常用、也最容易写对的优化是空间网格把平面按 cell_size 切成正方形单元格每个质点按所在单元格编号存入哈希表查询某质点邻居时只进入周围 3×3二维或 3×3×3三维的相邻单元格做距离判断。# spatial_hash.py import numpy as np def grid_key(pos, cell_size): 计算质点 pos 所在单元格的整数坐标 return (int(np.floor(pos[0] / cell_size)), int(np.floor(pos[1] / cell_size))) def build_map(positions, cell_size): hash: (cx, cy) - [质点索引...] table {} for i, p in enumerate(positions): key grid_key(p, cell_size) table.setdefault(key, []).append(i) return table def neighbor_indices(i, positions, table, cell_size): 返回与第 i 个质点可能接触的质点索引列表 key grid_key(positions[i], cell_size) hits [] for dx in (-1, 0, 1): for dy in (-1, 0, 1): k (key[0] dx, key[1] dy) for j in table.get(k, []): if j ! i: dist np.linalg.norm(positions[j] - positions[i]) if dist cell_size: hits.append(j) return hits逻辑说明build_map 把每个质点分配到单元格neighbor_indices 只查当前格及相邻格内的候选点再做一次精确距离判断来过滤。参数上 cell_size 是难点它要大于质点之间的最大相互作用距离通常是弹性碰撞直径的 1 到 2 倍太小会让同一对相互作用被切到两个格子导致漏检太大又退化成全量扫描。一个稳妥的起步值是粒子直径的 1.5 倍之后按漏检率调整。4.2 接触处理先修位置再看速度质点穿过边界那一帧如果你等速度反向后才去更新位置结果依然是穿透的。Verlet 风格的模拟更推荐做位置修正检测到越界后直接把位置提到边界上再根据修正后的位置重算本帧速度。墙体是最简单的例子# 假设地面在 y ground_y if pos_new[1] ground_y: pos_new[1] ground_y # 位置修正阻止穿透 vel[0] * 0.8 # 沿墙方向的摩擦衰减 vel[1] 0.0 # 法向速度清零这段代码说明的是修正顺序先位置后速度。摩擦系数 0.8 是保守起步值0 表示完全光滑、1 表示完全不打滑法向速度清零等价于恢复系数 0若想做出有弹性的落地改成vel[1] -restitution * vel[1]restitution 取 0.2 到 0.6 之间比较自然。惩罚力写法比位置修正更容易实现但调参要同时对付刚度和阻尼初学者往往调一晚上都得不到像样的反弹。4.3 并行化的三个容易踩的坑质点法本身很容易并行每个粒子求力的输入只是自己的位置、速度以及相邻粒子的位置。但把它写成多线程版本时有三件事比线程数更值得先想清楚。第一更新顺序。同一帧里线程 A 正在写粒子 i 的新位置线程 B 同时在读粒子 j 的旧位置去求力读到的数据就混了帧。常见做法是双缓冲计算时只读旧位置把新位置写进另一块内存整帧结束再交换指针。第二数据竞争。把空间哈希表建好后两个线程可能同时向同一个单元格追加索引加锁会把并行度拖垮。更稳的切法是按格子划分任务每个线程负责若干单元格内所有粒子的力计算线程之间在格子边界上交换一份只读的邻居列表写入阶段完全避开共享结构。第三内存访问模式。位置和速度分开存比把粒子对象 packed 在一起更有利于缓存命中。这里用 C 语言写一段示意// SoA 布局坐标和速度分别用独立 vector 存储 struct Particles { std::vectorfloat x, y, vx, vy; };用 x 和 y 分开的两个数组会让相邻粒子访问时连续读取同一段缓存而 AoS数组里每个元素是一整个 Particle 对象的连续访问会夹杂着无关字段。对 n 上万、迭代几万次的模拟这个布局差异能省下可观的时间。要不要上 SIMD取决于你的编译器是否支持自动向量化先跑基线再逐个打开优化选项对比不要一开始就上手工 intrinsics。5. 验证质点法算对了能量曲线、CSV 导出与轨迹对比5.1 用总能量曲线发现积分器的隐患结构完整的质点法代码长时间跑下来总能量应当在一个窄带内波动如果只盯着某一个质点的位置看很难区分是物理现象还是数值误差。我一般让模拟每 10 步输出一次位置、速度、重力势能和动能脚本汇总成一条能量曲线# energy_check.py def total_energy(pos_list, vel_list, gravity9.8, mass1.0): kin 0.5 * mass * sum(np.dot(v, v) for v in vel_list) pot -mass * gravity * sum(p[1] for p in pos_list) # 以 y0 为势能零点 return kin pot # 假设模拟过程中按帧记录了两个快照数组 for frame_i, (pos, vel) in enumerate(zip(frames_pos, frames_vel)): e total_energy(pos, vel) print(frame_i, e)量级参照以没有外力时的单体自由平移为例总能量应该完全不变出现重力和弹簧时能量曲线允许小幅波动。波动幅度超过总量的 1%先做一件事把 dt 除以 2 重跑一次如果波动一起缩小说明问题是数值解在注入能量而不是模型写错。如果能量曲线形状对但数值上整体偏移优先检查势能零点设定。5.2 导出低开销的 CSV 与轨迹重叠对比zip 包里如果自带 data 输出接口通常给的是二进制序列或按帧写出的 CSV。我推荐的导出列固定为frame,id,x,y,vx,vy一帧一个文件或者合并成一个总文件方便 pandas 直接切分import csv def save_trajectory(history, filenametrajectory.csv): with open(filename, w, newline) as f: writer csv.writer(f) writer.writerow([frame, id, x, y, vx, vy]) for frame_id, particles in enumerate(history): for pid, p in enumerate(particles): writer.writerow([frame_id, pid, p.x[0], p.x[1], p.v[0], p.v[1]])把同一初始条件用 dt 1/240 和 dt 1/480 各跑一遍导出两份 CSV在绘图脚本里按 frame_id 对齐画轨迹线。两条曲线在第 100 帧之前几乎重合之后逐渐分开是正常的但如果分开的时间点出现在第 20 帧之前说明系统对初始扰动过于敏感要么是刚度参数太高要么是积分误差超过预期。最后一个小技巧是给轨迹曲线叠加半透明的速度向量箭头能立刻看出哪一段路径上的速度方向与受力方向不一致通常那就是积分器在该处发散的起点。本文还有配套的精品资源点击获取 质点法练习包全解析：从解压到Verlet积分与碰撞优化
简介压缩包chp5_ex2.zip是配合丁承先向量式结构力学内容编写的质点法练习二面向正在学习计算结构力学、MATLAB编程的本科生或工程师目的是通过一个可直接运行的源代码示例理解将连续体离散为质量点、按相互作用力建立运动方程并数值求解的完整流程。质点法的核心优势在于无需复杂网格生成对非均匀、复杂几何结构尤其有效因此这一练习也适合作为从有限元思维过渡到无网格/粒子方法的入门参考。包内仅含1个m文件压缩后约2KB文件虽小但结构完整覆盖质点位置与刚度定义、弹性/重力等相互作用力计算、牛顿第二定律组建代数方程组、欧拉或龙格-库塔时间积分以及边界条件处理和结果输出等多个关键环节。已有116人学习浏览对于想快速上手质点法的读者来说该示例能提供从理论到代码的直观映射运行和分析时还可借鉴到适合结构力学场景的矩阵构建与求解思路并能进一步迁移到流体、地震工程等更广领域。1. 拿到 chp5_ex2.zip质点法练习的第二部分要先想清楚什么chp5_ex2.zip 这种命名方式在图形学与物理仿真的课程包里很常见chp5 表示第 5 章ex2 表示第 2 个练习后面跟的“质点法2”说明这不是第一次接触粒子类模拟。压缩包里边通常是源码、初始数据和题目说明你拿到手的第一反应不该是双击运行而是先把“这套代码在算哪种质点法、要你改哪一块”搞清楚。按我的分类质点法是一个宽泛说法既包括游戏引擎里的粒子系统也包括质点弹簧、光滑粒子流体动力学SPH和物质点法MPM它们的共同点是用一堆离散质点代替连续的物体或场然后用牛顿第二定律驱动每个点运动。这篇文章就顺着 chp5_ex2 这个练习包的常见结构把解包、最小实现、积分器选择、碰撞与并行化、以及如何验证结果一条线讲到位。2. 把 chp5_ex2.zip 摊开目录结构、依赖检查与质点法最小模型2.1 先看压缩包里有什么再决定用哪条命令把它解开我的习惯是解压之前先列清单不解压也能知道包里是不是你要的东西。用 unzip 的 -l 选项只打印压缩包内的文件列表不落盘unzip -l chp5_ex2.zip输出会列出每个文件的权限、大小、日期和完整路径。看到 src/、data/、config/ 这类目录结构基本就能判断是一个带工程骨架的练习包而不是几页散装代码。确认无误后解压到独立目录避免把一堆文件直接撒在下载目录里mkdir -p chp5_ex2 unzip chp5_ex2.zip -d chp5_ex2-d 参数指定目标目录名后续编译或运行都在 chp5_ex2 下进行不会污染其他地方。如果这一行报了 error read zip archive 之类的错先跑一遍unzip -t chp5_ex2.zip做完整性测试多半是文件在传输中被截断这种情况重下远比修 zip 省时间。我整理了一份课程练习包常见的布局对照路径典型内容练习时需要动的地方src/质点类、力模型、积分器实现题目要求补全或改写的部分data/初始坐标、边界盒、材质参数换测试场景时在这里改config/JSON、YAML 或 INI 参数每轮实验都会调README.md题目说明、编译命令、运行要求动手之前先读如果你是从 GitHub 上下载的 zip 包目录里通常还会多一层仓库根目录解压后 cd 进去即可其余逻辑和这里一致。把目录结构看过一遍再去看 README 里要求的依赖版本一般练习包的代码不会依赖太新的编译器或 SDK但 Python 类的练习经常要 numpy、scipy 中的固定小版本有版本偏差时先按项目要求的版本建虚拟环境省得后面跑出莫名其妙的报错。提示unzip -t 只校验压缩包内部数据的完整性不对每个文件做业务层面的校验。文本文件解出来是乱码通常是编码问题不是压缩包损坏。2.2 质点法的力学模型不管文件多少核心就 F ma质点法的出发点是拉格朗日视角不追踪网格而是把质量离散到一组点上。每个质点的状态由位置 x 和速度 v 描述运动方程可以写成两个一阶常微分方程组dx/dt vdv/dt F(x, v) / m这里的 F 是所有外力与内力之和。课程练习的第一部分通常只做自由粒子在重力场中的运动到了“质点法2”这个位置一般会加入粒子与粒子之间的相互作用也就是在 F 里补上弹簧力、阻尼力或者压力梯度项。写实现时我不建议把每种力单独写成一段顺次执行的代码而是先定义一个累加器把每个力模型的返回值加进去全部累加完成后再除以质量得到加速度。这样做的好处是新增一种力只需要补一个函数不需要改动积分流程。常见模型的力表达式集中列成一张表方便对照力的类型表达式作用重力F m * g提供外力驱动方向恒定弹簧胡克F -k * (r线性阻尼F -c * v_rel消耗相对动能抑制持续震荡边界碰撞法向速度反向或位置修正防止质点穿出计算域计算顺序也有讲究先求出同类质点与邻近质点的相互作用再叠加重力和边界力最后统一用积分器推进。练习包里常见的错误是把重力在每次迭代里重复乘以 dt或者把阻尼因子写成 dt 的幂次这两种写法都会让结果偏离物理直觉调参时尤其难找原因。2.3 最小可运行示例用纯 Python 把一次自由落体模拟跑通在没有题目代码的情况下我会先用一段最小实现验证环境再往里面加约束。下面的脚本只包含一个质点和重力外加 200 步显式欧拉积分# demo_particle.py import numpy as np class Particle: def __init__(self, x, v, mass1.0): self.x np.array(x, dtypefloat) # 初始位置 [x, y] self.v np.array(v, dtypefloat) # 初始速度 [vx, vy] self.mass mass # 质量默认 1.0 def step(self, dt, gravity9.8): # 只算重力方向朝下 force np.array([0.0, -self.mass * gravity]) acc force / self.mass self.v self.v acc * dt # 先更新速度 self.x self.x self.v * dt # 再更新位置 p Particle([0.0, 10.0], [1.0, 0.0], mass1.0) for _ in range(200): p.step(0.01) print(final x , p.x, final v , p.v)逻辑上这版实现遵循先速度后位置的显式欧拉更新顺序加速度恒定为 9.8 向下速度增量等于加速度乘以 dt位置增量等于更新后的速度乘以 dt。参数上 dt0.01 是迭代步长200 步对应 2 秒的物理时间最终位置的 y 坐标大约在 4.6 附近。想要验证解析解把 dt 改小到 0.002 并跑 1000 步数值结果会向 y 10 - 0.5 * 9.8 * 2² 靠近。这段最小代码跑通了再往 src 里加弹簧约束心里才有底。3. 时间积分与弹簧约束质点法参数调优的三个关键点3.1 显式欧拉不是不能用但你要知道它误差从哪里来显式欧拉的每一步都沿着当前位置的切线走直线切线方向会随着曲率变化而偏移累加几百步之后轨迹和真实解之间的偏差就会肉眼可见。在弹簧或轨道这类强弯曲运动中显式欧拉会让系统看起来“越蹦越高”本质是离散化在向系统内持续注入能量。质点法练习包里几乎都会给你一个积分器替换任务从欧拉改成 Verlet。位置 Verlet 不需要显式保存速度它用上一帧和当前帧的位置推算下一帧位置加速度只在这里出现def euler_update(x, v, accel, dt): v_new v accel * dt x_new x v_new * dt return x_new, v_new def verlet_update(x, x_prev, accel, dt): # x_prev 是上一帧位置accel 是当前受力计算的加速度 x_next 2.0 * x - x_prev accel * (dt * dt) return x_next两种函数对比Verlet 少了一次速度存储代码更短但在碰撞处理里需要恢复速度时反而要多一次差分运算。实际工程里我一般这么选粒子系统用位置 Verlet因为它在约束投影场景下实现简单需要在碰撞瞬间获得准确速度的刚体粒子模拟用 RK2 或 half-step Verlet。half-step Verlet 的做法是先用半个步长更新速度再用新速度更新一整步位置最后再用半个步长更新速度速度的获取代价只增加一次加法。参数上dt 直接决定两类误差截断误差随 dt 变小而变小但浮点舍入误差随迭代次数增加而累积。对弹簧系统先跑一个保守的 dt 1/240也就是每帧最多 4 个子步观察最大位移是否随时间线性增长如果位置曲线在第 2000 步附近出现折返就说明积分器在共振点上发散了。3.2 弹簧约束与阻尼把一堆独立的点做成一块布质点法练习到了第二部分最常见的任务是把自由粒子用弹簧连起来模拟布料或弹性绳。弹簧力方向沿两质点连线大小按当前伸长量与静止长度的差计算再叠加沿连线方向的阻尼速度# spring_force.py import numpy as np def spring_force(a, b, rest_len, k, c): a, b 为两个质点对象需要 x 和 v 属性 rest_len 是弹簧静止长度k 是刚度c 是阻尼系数 r b.x - a.x dist np.linalg.norm(r) if dist 1e-9: # 防止除零 return np.zeros(2), np.zeros(2) direction r / dist # 从 a 指向 b 的单位向量 v_rel b.v - a.v force_mag -k * (dist - rest_len) - c * np.dot(v_rel, direction) fa force_mag * direction # a 受到的力 fb -fa # 作用力与反作用力 return fa, fb参数上 k 是拉伸刚度数值越大抵抗形变越强c 是阻尼系数控制相对运动被削弱的程度。c 相对 k 过小时弹簧会沿着轴向反复震荡表现为布料上出现持续抖动c 相对 k 过大时相邻粒子黏在一起布料像面团。以粒子质量 m1 为单位我通常会从 k100、c0.5 起步再把 k 每档降低一倍看效果。调参时最容易忽略的是 rest_len 与初始网格边长的一致性。网格生成时如果相邻粒子距离是 0.3而 rest_len 被设成 1.0弹簧从第一帧就在做大行程收缩系统会在收敛前产生剧烈的初始振动。先把 rest_len 打印出来和初始网格边长对比再去调 k、c顺序不要反。下面是组合参数表按上面这段代码直接可用参数起步值调大后的效果常见的坑rest_len等于初始网格边长布料被拉伸或压缩与初始距离不匹配会先振动k100更硬约束更强k 过大配大 dt 会爆c0.5振动衰减更快c 过大让运动迟钝dt1/240单步更准太小则总耗时线性上升3.3 时间步长怎么压先看波形再看数值接手一个几百行的质点法练习代码我的调参顺序是这样的先不碰任何物理参数只改 dt把 dt 1/240 和 dt 1/480 各跑一遍对比同一时刻的质点位置。如果两条轨迹几乎重合说明当前 dt 已经足够小可以回退到 1/240 省时间如果差距超过 5%说明约束或碰撞引入了大量高频分量继续细化 dt 才是正路而不是转而加大阻尼去掩盖问题。判断发散不看单帧截图要看趋势。每 100 步打印一次最大速度的模长出现连续蹿升或者出现 inf/NaN立刻把 dt 减半重跑。注意这里说的减半不是把迭代次数减半而是在相同物理时长内加密子步数所以总耗时不变。还有一种现象是代码本身没毛病但报错出在文件读取而不是积分器上例如配置里写了data/particles.csv实际路径是config/../data/particles.csv。这类问题在练习包里出现频率很高解包后先把 data 目录相对路径和 README 对照一遍比我刚才提到的 error read zip archive 更隐蔽因为它是运行时才出现而且报错内容不会提示路径深度不对。注意打印最大速度时不要用自带的 print 直接刷屏每 100 步输出一行就够了否则日志文件会拖慢整体模拟。4. 大规模质点模拟空间哈希、碰撞修正与并行改写边界4.1 空间网格哈希把 O(n²) 的碰撞查询压下来几千个质点彼此两两判断距离一次迭代就要上千万次运算浏览器里的 JS 演示还能勉强跑桌面程序也扛不住上万的 n。最常用、也最容易写对的优化是空间网格把平面按 cell_size 切成正方形单元格每个质点按所在单元格编号存入哈希表查询某质点邻居时只进入周围 3×3二维或 3×3×3三维的相邻单元格做距离判断。# spatial_hash.py import numpy as np def grid_key(pos, cell_size): 计算质点 pos 所在单元格的整数坐标 return (int(np.floor(pos[0] / cell_size)), int(np.floor(pos[1] / cell_size))) def build_map(positions, cell_size): hash: (cx, cy) - [质点索引...] table {} for i, p in enumerate(positions): key grid_key(p, cell_size) table.setdefault(key, []).append(i) return table def neighbor_indices(i, positions, table, cell_size): 返回与第 i 个质点可能接触的质点索引列表 key grid_key(positions[i], cell_size) hits [] for dx in (-1, 0, 1): for dy in (-1, 0, 1): k (key[0] dx, key[1] dy) for j in table.get(k, []): if j ! i: dist np.linalg.norm(positions[j] - positions[i]) if dist cell_size: hits.append(j) return hits逻辑说明build_map 把每个质点分配到单元格neighbor_indices 只查当前格及相邻格内的候选点再做一次精确距离判断来过滤。参数上 cell_size 是难点它要大于质点之间的最大相互作用距离通常是弹性碰撞直径的 1 到 2 倍太小会让同一对相互作用被切到两个格子导致漏检太大又退化成全量扫描。一个稳妥的起步值是粒子直径的 1.5 倍之后按漏检率调整。4.2 接触处理先修位置再看速度质点穿过边界那一帧如果你等速度反向后才去更新位置结果依然是穿透的。Verlet 风格的模拟更推荐做位置修正检测到越界后直接把位置提到边界上再根据修正后的位置重算本帧速度。墙体是最简单的例子# 假设地面在 y ground_y if pos_new[1] ground_y: pos_new[1] ground_y # 位置修正阻止穿透 vel[0] * 0.8 # 沿墙方向的摩擦衰减 vel[1] 0.0 # 法向速度清零这段代码说明的是修正顺序先位置后速度。摩擦系数 0.8 是保守起步值0 表示完全光滑、1 表示完全不打滑法向速度清零等价于恢复系数 0若想做出有弹性的落地改成vel[1] -restitution * vel[1]restitution 取 0.2 到 0.6 之间比较自然。惩罚力写法比位置修正更容易实现但调参要同时对付刚度和阻尼初学者往往调一晚上都得不到像样的反弹。4.3 并行化的三个容易踩的坑质点法本身很容易并行每个粒子求力的输入只是自己的位置、速度以及相邻粒子的位置。但把它写成多线程版本时有三件事比线程数更值得先想清楚。第一更新顺序。同一帧里线程 A 正在写粒子 i 的新位置线程 B 同时在读粒子 j 的旧位置去求力读到的数据就混了帧。常见做法是双缓冲计算时只读旧位置把新位置写进另一块内存整帧结束再交换指针。第二数据竞争。把空间哈希表建好后两个线程可能同时向同一个单元格追加索引加锁会把并行度拖垮。更稳的切法是按格子划分任务每个线程负责若干单元格内所有粒子的力计算线程之间在格子边界上交换一份只读的邻居列表写入阶段完全避开共享结构。第三内存访问模式。位置和速度分开存比把粒子对象 packed 在一起更有利于缓存命中。这里用 C 语言写一段示意// SoA 布局坐标和速度分别用独立 vector 存储 struct Particles { std::vectorfloat x, y, vx, vy; };用 x 和 y 分开的两个数组会让相邻粒子访问时连续读取同一段缓存而 AoS数组里每个元素是一整个 Particle 对象的连续访问会夹杂着无关字段。对 n 上万、迭代几万次的模拟这个布局差异能省下可观的时间。要不要上 SIMD取决于你的编译器是否支持自动向量化先跑基线再逐个打开优化选项对比不要一开始就上手工 intrinsics。5. 验证质点法算对了能量曲线、CSV 导出与轨迹对比5.1 用总能量曲线发现积分器的隐患结构完整的质点法代码长时间跑下来总能量应当在一个窄带内波动如果只盯着某一个质点的位置看很难区分是物理现象还是数值误差。我一般让模拟每 10 步输出一次位置、速度、重力势能和动能脚本汇总成一条能量曲线# energy_check.py def total_energy(pos_list, vel_list, gravity9.8, mass1.0): kin 0.5 * mass * sum(np.dot(v, v) for v in vel_list) pot -mass * gravity * sum(p[1] for p in pos_list) # 以 y0 为势能零点 return kin pot # 假设模拟过程中按帧记录了两个快照数组 for frame_i, (pos, vel) in enumerate(zip(frames_pos, frames_vel)): e total_energy(pos, vel) print(frame_i, e)量级参照以没有外力时的单体自由平移为例总能量应该完全不变出现重力和弹簧时能量曲线允许小幅波动。波动幅度超过总量的 1%先做一件事把 dt 除以 2 重跑一次如果波动一起缩小说明问题是数值解在注入能量而不是模型写错。如果能量曲线形状对但数值上整体偏移优先检查势能零点设定。5.2 导出低开销的 CSV 与轨迹重叠对比zip 包里如果自带 data 输出接口通常给的是二进制序列或按帧写出的 CSV。我推荐的导出列固定为frame,id,x,y,vx,vy一帧一个文件或者合并成一个总文件方便 pandas 直接切分import csv def save_trajectory(history, filenametrajectory.csv): with open(filename, w, newline) as f: writer csv.writer(f) writer.writerow([frame, id, x, y, vx, vy]) for frame_id, particles in enumerate(history): for pid, p in enumerate(particles): writer.writerow([frame_id, pid, p.x[0], p.x[1], p.v[0], p.v[1]])把同一初始条件用 dt 1/240 和 dt 1/480 各跑一遍导出两份 CSV在绘图脚本里按 frame_id 对齐画轨迹线。两条曲线在第 100 帧之前几乎重合之后逐渐分开是正常的但如果分开的时间点出现在第 20 帧之前说明系统对初始扰动过于敏感要么是刚度参数太高要么是积分误差超过预期。最后一个小技巧是给轨迹曲线叠加半透明的速度向量箭头能立刻看出哪一段路径上的速度方向与受力方向不一致通常那就是积分器在该处发散的起点。本文还有配套的精品资源点击获取