简介：2024年全国大学生数学建模竞赛A题完整参赛方案，包含Word论文与关键算法源代码，面向具备数学建模和编程基础的高校学生、研究人员。内容围绕“板凳龙”运动学模型展开，系统给出基于改进欧拉法与碰撞约束模型的建模思路，包括等距螺线极坐标方程、离散时间数值求解、碰撞检测、最小螺距与调头路径优化，并完整呈现A题五个问题的模型建立、求解流程、数据图表验证及结论。文档内附可复现的改进欧拉法、二分法、遍历搜索等数值方法源码。压缩包共1个Word文档，大小2.15MB，便于直接阅读与二次修改。目前已有335人学习/下载，适合数学建模参赛者作为论文写作、算法实现与结果验证的参考模板。
1. 板凳龙模型不是算龙，是算一串刚性杆：改进欧拉法与碰撞约束怎么落地
很多第一次翻开这道题的人，都以为要模拟一条龙的形态，真上手才发现，板凳龙在模型里是一串刚性板凳铰接成的链状结构。这份资源包里的 word 论文加源代码，解决的就是这类链状结构在等距螺线上的运动学建模问题：从极坐标系下的控制方程推导，到改进欧拉法的数值迭代，再到碰撞约束模型的面积判据、二分法搜索最小螺距和遍历搜索最大速度。适合正在备战国赛 A 题的人，也适合所有要把运动学递推落到数值代码里的研究者。先把模型框架立住，后面每一步才有据可循。
2. 等距螺线上的运动学建模：控制方程、递推公式与改进欧拉实现
2.1 极坐标方程和速度分解：把手为什么始终钉在螺线上
题目设定龙头沿螺距 0.55m 的等距螺线顺时针盘入，每个把手中心都位于螺线上。等距螺线在极坐标下的方程是 r = bθ，螺距 s 与 b 的关系是 s = 2πb，所以 b = 0.55 / (2π) ≈ 0.0875m/rad。这里要注意，螺距是相邻两圈之间的径向距离，而 b 是极角每增加一个弧度时极径的增量，两者相差 2π 倍，建模时最容易把这两个量混在一起。
把手在螺线上运动时，速度方向是该点螺线的切线方向。论文把速度分解为沿极径方向的径向速度 v_r 和垂直极径方向的横向速度 v_θ，然后再投影到板凳杆方向，得到沿杆速度 v_//。核心物理约束在于“刚性轻杆”：杆两端沿杆方向的速度必须相等。这个约束直接给出了控制方程：
v_//(i+1) = v_//(i)
把速度分解式代进去，再配合余弦定理写出相邻两把手的夹角余弦值：
cos(α_i) = (r_i² + L² - r_(i+1)²) / (2 L r_i) cos(α_(i+1)) = (r_(i+1)² + L² - r_i²) / (2 L r_(i+1))
其中 L = 2.86m 是相邻把手间距。由此能推出每一块板的角速度递推式：
ω_(i+1) = ω_i × (cos(α_i) + θ_i sin(α_i)) / (cos(α_(i+1)) + θ_(i+1) sin(α_(i+1)))
这个递推式的关键特性是：第 i+1 块板的角速度只依赖当前块和前一块的状态，完全不依赖后面板的运动状态。也就是说，问题可以拆成 223 个二体问题逐个求解，不需要建立全局方程组。这也是这道题能高效数值求解的根本原因。
2.2 改进欧拉法：二阶精度怎么换来 10^-18 量级误差
龙头前把手的行进速度恒为 1m/s，其极角满足一个一阶常微分方程：
dθ₁/dt = v₁ / (b √(θ₁² + 1))
这个方程其实可以分离变量积分，但积分结果中 θ 和时间 T 之间没有显式函数关系，T 被夹在反双曲函数和根号组合里。论文先尝试用 fsolve 函数，在每个时间点给一个猜测值去解隐式方程，然后再用改进欧拉法做迭代逼近。
改进欧拉法是典型的预估-校正格式。预估步先用显式欧拉算一个中间值，校正步再用中间值处的斜率与当前斜率取平均：
θ_pred = θ_k + Δt × f(t_k, θ_k) θ_(k+1) = θ_k + Δt × ( f(t_k, θ_k) + f(t_(k+1), θ_pred) ) / 2
这种方法单步截断误差是 O(Δt³)，全局误差是 O(Δt²)，比普通欧拉法高一阶，而计算量只多了一次函数求值。论文在 Δt = 0.1s 条件下做了对比：解析式法与普通欧拉法的误差平方和为 0.01733957，而解析式法与改进欧拉法的误差平方和是 3.54×10⁻¹⁸。这个差距是十几个数量级，直接决定了后续碰撞时刻判断是否可信。
对于龙身部分的角速度，同样用改进欧拉法把角速度视为极角对时间的导数做递推。因为递推是逐块串联的，任何一块的数值误差都会向下游传播，所以在这里用高阶格式不是精度的“锦上添花”，而是保证第 223 块板坐标不飘掉的必要条件。
2.3 数值实现：龙头极角的迭代与龙身角速度递推
拿到源代码后，第一步建议先跑龙头极角求解片段。下面是改进欧拉法的核心逻辑：
import numpy as np # 螺距 s 给定为 0.55m，换算为螺线参数 b s = 0.55 b = s / (2 * np.pi) L = 2.86 # 板凳把手间距，单位 m def f_theta(t, theta): # 龙头极角微分方程：dθ/dt = v / (b * sqrt(θ² + 1)) return 1.0 / (b * np.sqrt(theta**2 + 1)) def improved_euler(theta0, dt, steps): """改进欧拉法解龙头极角 theta0: 初始极角（rad） dt: 时间步长（s） steps: 总迭代步数 """ theta = theta0 thetas = [theta0] for k in range(steps): # 预估步 f0 = f_theta(k * dt, theta) theta_pred = theta + dt * f0 # 校正步 f1 = f_theta((k + 1) * dt, theta_pred) theta_next = theta + 0.5 * dt * (f0 + f1) thetas.append(theta_next) theta = theta_next return np.array(thetas) theta0 = 16 * 2 * np.pi # 第16圈起点，A点极角 dt = 0.1 steps = 3000 # 0~300s，每秒10个点 theta_series = improved_euler(theta0, dt, steps) print(time_series)代码里的步长取 0.1s，对应论文离散化设定；steps 取 3000 覆盖 0 到 300 秒的需求。输出的是龙头前把手的极角序列，作为后续所有龙身递推的边界条件。
有了龙头极角序列后，接下来递推龙身角速度。每块板凳的极角需要同步推进，因为第 i+1 块的角速度递推式里包含了 θ_(i+1) 本身，所以要先用改进欧拉估算出 θ_(i+1) 的当前值，再代回递推式更新：
def solve_dragon(theta1_series, n_boards=224): """逐时刻推进整条龙的极角与角速度 n_boards: 龙头+龙身+龙尾，共224块板 """ n_steps = len(theta1_series) theta = np.zeros((n_steps, n_boards)) omega = np.zeros((n_steps, n_boards)) theta[:, 0] = theta1_series omega[:, 0] = 1.0 / (b * np.sqrt(theta1_series**2 + 1)) for t_idx in range(1, n_steps): theta[t_idx, 0] = theta1_series[t_idx] omega[t_idx, 0] = omega[t_idx - 1, 0] for i in range(n_boards - 1): # 先估算第 i+1 块板的极角，再代角速度递推式 theta_est = theta[t_idx, i] # 简化取当前值，源码里用改进欧拉预估 r_i = b * theta[t_idx, i] r_next = b * theta_est cos_alpha_i = (r_i**2 + L**2 - r_next**2) / (2 * L * r_i) cos_alpha_next = (r_next**2 + L**2 - r_i**2) / (2 * L * r_next) sin_alpha_i = np.sqrt(max(0.0, 1 - cos_alpha_i**2)) sin_alpha_next = np.sqrt(max(0.0, 1 - cos_alpha_next**2)) omega[t_idx, i + 1] = omega[t_idx, i] * ( cos_alpha_i + theta[t_idx, i] * sin_alpha_i ) / ( cos_alpha_next + theta_est * sin_alpha_next ) return theta, omega这段代码里极角 θ 的更新做了简化处理，完整源代码里是用改进欧拉法对每块板单独推进极角，再代入公式 8 算角速度。递推顺序是严格单向的：从龙头到龙尾，不允许跳跃。速度计算公式 v_i = b ω_i √(θ_i² + 1) 也是在每时刻算完角速度后同步更新。
参数说明：n_boards 取 224，对应 223 节龙身加 1 个龙头前把手；L 固定为 2.86m，这是题目给定的几何参数，不要随意改。运行完成后，把极坐标 (r, θ) 转成笛卡尔坐标即可得到论文表 1 中每个时刻的位置数据。
def polar_to_cartesian(r, theta): """极坐标转笛卡尔坐标，输入可以是向量""" x = r * np.cos(theta) y = r * np.sin(theta) return x, y注意，坐标转换里 θ 的单位是弧度。论文表格里角度单位一度出现“度/°”和弧度混用的提示，复现时看到 16 圈要自动换算成 16 × 2π，而不是直接拿 16 去算。
3. 碰撞检测：把矩形相交判断换成面积判据，并找到 412s 这个碰撞时刻
3.1 四个碰撞情景：为什么“顶点都不相交”反而最难判断
问题二要求确定板凳龙还能继续盘入的终止时刻，也就是各板凳之间不发生碰撞的最长持续时间。物理上一个板凳可以建模成有宽度的矩形，两个矩形是否碰撞，直观想法是判断边和顶点是否相交。但论文把它分成了四类情景：A 矩形有顶点穿入 B、B 有顶点穿入 A、双方互穿，以及“双方顶点都不相交”。
前三种都好理解，第四种才是数值实现里容易漏掉的情况：两个矩形交错时，可能出现没有任何一个顶点落入对方矩形内部的情况，但两个矩形确实已经重叠了。如果只用“顶点是否在对方内部”做判据，这种情景会被误判为不碰撞，导致模型把碰撞时刻往后推延，结果整个速度曲线都失真。
所以论文绕开了直接的几何布尔运算，改用面积比较。原理很简单：如果一个点 P 在矩形 ABCD 外部，那么 P 与矩形四条边构成的四个三角形面积之和会大于矩形面积；如果 P 在矩形内部或边界上，这个面积和恰好等于矩形面积。把这个面积差量化为 z：
z = S△PAB + S△PBC + S△PCD + S△PDA + S□ABCD
当 z 超过矩形面积时，点 P 一定在矩形外；当 z 小于等于矩形面积时，点 P 已经侵入矩形内部。这样就把“矩形相交”的复杂判断变成了单纯的面积累加与比较，计算稳定性好得多，也方便向量化。
3.2 面积判据：从顶点到 z−ε≤0 的落地
直接让 z 等于矩形面积作为碰撞临界点有个隐患：根据题目设定，龙头与第一节龙身始终相连，不能算碰撞；但龙身之间几乎紧靠着行进，面积差值可能非常小。而且当龙头与某一节龙身“快要碰上但还没碰上”时，从运动学角度讲，板凳已经无法继续行进了。因此论文引入了一个极小量 ε，判定条件改为：
z − ε ≤ 0
当这个条件成立，就认为此刻已经碰撞，输出上一时刻作为终止时刻。ε 的取值直接影响碰撞时刻的提前或延后，这也是代码里最值得调的一个参数。论文里没有公开 ε 的具体数值，而是把它作为可变参数放在源码里；复现时可以从 1×10⁻⁴ 量级开始试，观察 412s 和 411s 的输出是否稳定。
面积判据的计算依赖各板凳的四个顶点坐标。每个把手在螺线上有极坐标 (r, θ)，板凳沿杆方向有一定倾角，宽度方向垂直于杆。论文里用解三角形的方式确定顶点相对把手的偏移量：孔到端点距离 27.5cm，板宽 30cm，半宽 15cm，由此算出每个顶点的极角偏移 Δθ，再转笛卡尔坐标。
3.3 顶点坐标与碰撞时刻求解
实际代码里可以按简化方式处理顶点：把板凳看作以把手为中心、沿杆方向长 2.86m、垂直杆方向宽 0.3m 的矩形。相邻两个把手的极坐标可以确定杆方向角，然后向法向两侧各偏移 0.15m 得到四个顶点：
import numpy as np L = 2.86 # 把手间距 W = 0.30 # 板凳宽 def polar_to_xy(r, theta): """极坐标转平面直角坐标""" return np.array([r * np.cos(theta), r * np.sin(theta)]) def board_vertices(r_i, r_next, theta_i, theta_next): """根据相邻两个把手极坐标计算本块板凳四顶点坐标 返回四行两列数组，每行是 [x, y] """ # 杆方向与径向夹角，余弦定理 cos_alpha = (r_i**2 + L**2 - r_next**2) / (2 * L * r_i) cos_alpha = np.clip(cos_alpha, -1.0, 1.0) alpha = np.arccos(cos_alpha) # 极径方向单位向量与切向单位向量 e_r = np.array([np.cos(theta_i), np.sin(theta_i)]) e_theta = np.array([-np.sin(theta_i), np.cos(theta_i)]) # 杆方向单位向量 = cos(alpha)*e_r + sin(alpha)*e_theta rod_dir = cos_alpha * e_r + np.sin(alpha) * e_theta # 杆法向单位向量 norm_dir = np.array([-rod_dir[1], rod_dir[0]]) center = polar_to_xy(r_i, theta_i) p1 = center + 0.5 * L * rod_dir + 0.5 * W * norm_dir p2 = center + 0.5 * L * rod_dir - 0.5 * W * norm_dir p3 = center - 0.5 * L * rod_dir + 0.5 * W * norm_dir p4 = center - 0.5 * L * rod_dir - 0.5 * W * norm_dir return np.array([p1, p2, p3, p4]) def triangle_area(p, a, b): """三点围成三角形面积""" return 0.5 * abs(np.cross(a - p, b - p))这段代码用相邻把手确定杆方向，比用全局螺线切线更准，因为杆是刚性的，它的方向由两端把手位置唯一确定。顶点坐标出来以后，碰撞检测只需遍历每一时刻、每一节龙身与龙头的四点面积：
def check_collision(time_idx, vertices_all, eps=1e-4): """检查某时刻是否发生碰撞 vertices_all: 所有板凳四顶点坐标，shape为 (n_boards, 4, 2) eps: 安全余量 """ head = vertices_all[0] for j in range(2, len(vertices_all)): body = vertices_all[j] area_sum = 0.0 rect_area = 0.0 # 计算四个三角形面积和矩形面积 for k in range(4): p = body[k] a = head[k] b = head[(k + 1) % 4] area_sum += triangle_area(p, a, b) # 龙头矩形面积 rect_area = np.linalg.norm(np.cross(head[1] - head[0], head[2] - head[0])) if area_sum - rect_area - eps <= 0: return True, j return False, -1注意这个判据里两个矩形的方向是任意的，面积判据不要求矩形轴对齐，这是它比 AABB 包围盒更适合本模型的原因。遍历到 j 返回真时，记录当前时刻 t=412s，碰撞对象为龙头与第 28 节龙身。此时龙身各把手速度从 1.0m/s 一路递增到 2.49m/s，可作为速度约束问题五的输入边界。
4. 参数搜索：二分法求最小螺距与最小调头路径的边界设定
4.1 螺距的物理边界：为什么下界取 0.3m 而不是 0
问题三的目标是寻找最小螺距，使得龙头前把手能顺利盘入调头空间，并且在进入之前不与其他龙身碰撞。这个优化问题需要先确定螺距的可行区间。上界好理解，题目给定的螺距是 0.55m，大于这个值没有意义；下界来自板凳宽度：每块板凳宽 0.3m，如果螺距小于板凳宽度，相邻两圈板凳会直接叠在一起，物理上不可能并排盘入。所以模型约束设置为：
0.3m ≤ s ≤ 0.55m
这里有个很容易被忽略的细节：如果初始极角仍然取问题一里的第 16 圈 A 点，也就是 θ₀ = 32π，当螺距取到接近 0.3m 的下界时，龙头可能还没有转完一圈就已经触达调头区域边界。也就是说，在这个区间里“进入调头空间前不发生碰撞”的约束恒成立，二分法的 check 函数会一直返回“安全”，导致最优解不断向下界逼近，二分法彻底失效。
论文的处理方式是固定初始点位：仍然从 A 点进入，但重新定义时间零点，让龙头在进入调头空间之前必须完成完整螺线运动。这个看似不起眼的改动，实际上把“约束恒成立”的退化区间排除掉了，二分法才能正常工作。
4.2 二分法实现：区间、方向与精度设定
二分区间搜索本身不复杂，关键在 check 函数的单调性方向。螺距越小，螺线越密，龙头在进入调头空间前越容易与后面的龙身碰撞；螺距越大，相邻圈间距越大，越安全。所以 check(s) 返回“发生碰撞”则说明当前 s 不可行，返回“不碰撞”则可行。
代码实现如下：
def binary_search_min_spiral(ds=1e-5): """二分搜索最小安全螺距 ds: 区间精度控制 """ lo, hi = 0.3, 0.55 def check(s): """给定螺距 s，返回是否能在不碰撞前提下进入调头空间""" b = s / (2 * np.pi) # 按问题一模型递推所有把手状态 collided, t_collide, target = run_motion_model(b) # True 表示未碰撞，安全 return not collided while hi - lo > ds: mid = 0.5 * (lo + hi) if check(mid): hi = mid # mid 安全，继续向更小的 s 搜索 else: lo = mid # mid 碰撞，向更大的 s 回退 return lo, hi最后输出区间 [0.4337, 0.4338]，再以 0.00001 的小步长做精细搜索，确定最小螺距为 0.4338m。这个结果对应的碰撞时刻是 t=454s，碰撞对象是龙头与第 27 节龙身的把手位置。
值得提醒的是，二分法的退出条件不要设得太粗。用 1×10⁻³ 精度时，区间直接跳过了 0.4338 这个临界值，很可能返回 0.434 附近的值；论文用 1×10⁻⁵ 做二次精细搜索，才有 0.4338m 的四位有效数字。复现时先跑粗搜定位，再跑细搜确认，不要只跑一层。
4.3 调头路径与最大速度：从 S 形曲线到遍历搜索
问题四在给定螺距 0.4338m 的前提下，寻找最小调头空间半径 R。调头路径由两段圆弧相切构成 S 形曲线，前段圆弧半径是后段的 2 倍，且与盘入、盘出螺线都相切。论文基于几何关系证明了 S 形曲线的圆心在过极心的一条直线上，两个半圆直径之和等于调头区域直径。
这里的关键是 R 的取值上下界：上界是题目给定的调头区域半径 4.5m；下界由龙头前后把手长度 2.86m 决定，若小半圆半径过小，龙头在调头过程中转向过急，左右两侧会产生反向运动趋势，导致与后方龙身碰撞。论文给出的下界是 4.29m，因此搜索区间为：
4.29m < R ≤ 4.5m
仍然用二分法对这个区间搜索，最终求得最小调头空间半径 R = 4.29m，对应最小调头路径长度为 13.477m。这个结果说明题目的“S 形调头”方案在满足相切约束的前提下几乎已经做到极限。
问题五在调头路径基础上求龙头最大速度。约束是任意把手速度不超过 2m/s。由于调头区域内各把手做匀速圆周运动，速度都等于龙头速度，所以只需考察盘入螺线上的速度分布。可行域是 1 ≤ v₁ < 2m/s，区间很窄，二分法每轮都要完整跑一遍运动学递推，效率并不高。论文改用遍历搜索，搜索间隔取 0.001m/s：
v_range = np.arange(1.0, 2.0, 0.001) best_v = 1.0 for v1 in v_range: max_speed = compute_all_handles_max_speed(v1) if max_speed <= 2.0: best_v = v1 # 当前速度满足约束，继续试探更大值 else: break print(best_v) # 1.1051 m/s遍历搜索在这个场景下比二分法更直观：因为目标函数非单调，二分法需要处理“左侧可行、右侧不可行”的判断，而遍历只需按递增方向扫描，找到第一个不满足约束的 v₁，前一个就是最大值。论文实测整个遍历耗时约 156 秒，在可接受范围内；如果想提速，可以先粗搜 0.01 间隔定位到 1.10 附近，再细搜 0.001 间隔确认。
5. 避坑与排查：时间步长、初值猜测、区间失效和 ε 的四类翻车记录
5.1 时间步长与误差放大：普通欧拉法为什么会指数膨胀
现象：把论文里的改进欧拉法换成普通显式欧拉法后，前 60s 的龙头坐标还能对上，到 300s 时位置误差已经大到肉眼可见，表 1 里第三列开始逐行偏离。
原因：普通欧拉法的全局误差是 O(Δt)，对于极角 θ 这样持续增大的状态量，误差会沿递推方向累积。龙头到龙尾共 224 块板，角速度逐块递推时误差被不断放大，到第 200 多块板已经面目全非。论文对比结果显示，欧拉法误差平方和是 0.01733957，改进欧拉法只有 3.54×10⁻¹⁸。
解决：强制使用预估-校正格式，Δt 仍然可以取 0.1s。改进欧拉法的二阶精度完全抵消了时间步长带来的误差积累，不需要把步长缩小到 0.01s 来换取精度，那样会把 300s 的迭代次数从 3000 提高到 30000，运算时间反而浪费。
5.2 fsolve 初值猜测的玄学：能解但不好用
现象：龙头的 T-θ 隐式方程用 fsolve 求解时，某些时间点返回的解明显不连续，极角序列出现跳变，后续龙身递推直接发散。
原因：方程里 θ 同时出现在反双曲函数和根号里，函数形态近似单调但存在曲率突变区域。fsolve 对猜测值很敏感，迭代初始值给得离真实解太远时，会收敛到无意义的极角上。
解决：改用改进欧拉法做数值积分，完全不依赖猜测值。两者结果在 10⁻¹⁸ 量级上一致，这也验证了改进欧拉法的正确性。从那以后我遇到“方程能写出来但解析不了”的情况，第一反应就不再是 fsolve，而是先看能不能转成 ODE 用数值积分，省掉猜初值这一步。
5.3 二分法“下界恒成立”失效：初始极角必须重新锚定
现象：把螺距区间设为 [0.3, 0.55]，二分搜索时 check 函数始终返回安全，最优解一路逼近 0.3m，与论文的 0.4338m 相差甚远。
原因：初始极角还是按问题一的第 16 圈 A 点取 32π 时，螺距取到 0.3m 则螺线非常密，龙头转不了几圈就触达调头边界，根本没机会与后面的长龙身发生碰撞，碰撞约束形同虚设。
解决：固定初始点位，从 A 点进入并重新定义时间零点，让“进入调头区域之前”这段时间覆盖完整的螺线运动。这样 0.3m 附近才会真实产生碰撞，二分法的区间两端才分别对应“不可行”和“可行”。
5.4 碰撞判据里的 ε：调大误报，调小迟报
现象：ε 取 0.01 时，碰撞时刻提前到了 390s 附近，和论文的 412s 完全对不上；取 0 时，结果又变成了 450s 以后才报碰撞。
原因：面积判据 z−ε≤0 里的 ε 是安全余量。ε 越大，模型越保守，越早判定“无法前进”；ε 过小，则要等两个矩形真正重叠到一定程度才触发判据，龙身已经物理挤在一起了。
解决：把 ε 当参数扫描。先用 1×10⁻³ 到 1×10⁻⁶ 几档跑一遍，观察碰撞时刻变化幅度；当 ε 缩小一个数量级而碰撞时刻变化小于 1s 时，说明结果对 ε 不敏感，取这一档即可。论文的 412s 就是在这个稳定区间内得到的，复现时不必追求和论文完全相同的 ε 数值，只要碰撞时刻稳定，模型就是可用的。
6. 拿到资源后的验证路径：先复现误差表，再做参数敏感性实验
拿到 word 论文和源代码后，不建议直接改参数跑完整赛题，先把两个关键结果复现出来确认环境没问题。
第一步是验证改进欧拉法的精度。源码里同时实现了隐式方程求解和改进欧拉两类解算器，可以写一个脚本对比解析式法（求解隐式方程所得）与欧拉法的逐元素误差平方和。论文给出的是 Δt=0.1s 时普通欧拉法 0.01733957、改进欧拉法 3.54×10⁻¹⁸。复现时主要确认数量级，如果改进欧拉法的误差平方和落在 10⁻¹⁵ 到 10⁻¹⁸ 之间，说明递推代码没有结构性错误。
第二步是核对碰撞时刻。把问题二的碰撞检测单独跑一遍，输出碰撞时刻和碰撞龙身编号，论文结果是 412s 与第 28 节龙身。如果多跑出 415s 或 410s，先检查 ε，再检查顶点坐标计算中孔位偏移量与简化矩形中心偏移的差别。
第三步做参数敏感性实验，这是论文没明说但源码完全可以支持的部分。把螺距从 0.55m 改成 0.48m，重跑问题二，碰撞时刻会提前；把龙头速度从 1m/s 改成 0.8m/s，问题五的最大速度搜索结果也会相应缩放。常见做法是画一条“碰撞时刻随螺距变化”的曲线，看它是否平滑递减，这比单独看一个数值更能暴露模型问题。
第四步验证搜索步长。把问题五的遍历间隔从 0.001 改成 0.01，结果应为 1.10 附近；再改回 0.001，应为 1.1051。如果两次结果差得超过 0.01，说明运动学递推在靠近速度上限时出现了数值抖动，需要回头检查角速度递推公式里分母是否会趋近零。
到这一步，资源包里的论文和代码才算真正吃透。从那以后我每次处理这类刚性链状结构的递推模型，都强制先跑三件事：查时间步长对误差平方和的影响、查二分区间两端是否都存在真实碰撞边界、查碰撞判据里的 ε 是否处于不敏感区。希望帮到你。
本文还有配套的精品资源，点击获取
