简介一份哈尔滨工业大学版《理论力学》教学课件面向工科院校学生与工程技术人员适合在学完静力学、运动学后系统学习动力学核心内容。资源聚焦第十章“质点动力学”以牛顿三定律为基础详细梳理动力学基本定律、国际单位制与工程单位制并重点讲解质点运动微分方程的矢量形式、直角坐标形式和自然坐标形式进而归纳质点动力学的两类基本问题已知运动求力、已知力求运动。课件结合曲柄连杆机构经典例题演示如何选取研究对象、建立运动方程并利用微分积分求解帮助读者掌握动力学分析的一般方法。全包内含1个演示文稿文件共1010KB结构紧凑、要点突出可用于课堂同步学习、课后复习或考研一轮梳理。该资源已有5626人学习下载是深入理解力与运动关系、夯实理论力学基础的实用资料。1. 质点动力学从静力平衡到运动变化的那一步大部分工程师接触力学是从静力学开始的一个物体放在那里受力平衡求出约束反力就算完事。到了运动学阶段我们开始描述位置、速度和加速度但完全不问“是什么力让它变成这样”。真正把这两者接起来的是动力学。而动力学里最先遇到的就是质点动力学——把物体抽象成一个有质量但没有尺寸的点研究它受的力与运动变化之间的关系。哈工大版的《理论力学》第十章就是围绕这个核心展开的三个基本定律、运动微分方程、以及两类典型问题。这里有个容易被忽略但很关键的点第一类问题是“已知运动求力”只需要微分第二类问题是“已知力求运动”需要积分。前者是正问题后者是反问题。实际工程里大多数情况是已知激励求响应也就是第二类但第一类在传感器数据处理、逆动力学分析里同样常见。本文从定律和单位制讲起把三类坐标系下的微分方程写清楚然后给出曲柄连杆、圆锥摆、弹簧振子三个例题的Python数值解法最后补一点关于数值积分稳定性和精度验证的实操经验。2. 三大定律与单位制先搞清楚Fma里每个符号的物理边界2.1 牛顿第一定律的真正用途定义惯性参考系牛顿第一定律惯性定律描述的是质点不受力或合力为零时它对惯性参考系保持静止或匀速直线运动。很多初学者把它当作第二定律的推论实际上它在逻辑上更靠前——它定义了“惯性参考系”这个概念。没有惯性参考系第二定律中的加速度就没有可对比的基准。地面通常被近似为惯性系但在涉及地球自转的精密计算里科氏力就必须被考虑进来这时地面不再是严格的惯性系。2.2 牛顿第二定律的矢量性与瞬时性第二定律的表达式是 F ma但要注意两点第一这是矢量方程F 和 a 方向相同第二它是瞬时关系某一时刻的力决定该时刻的加速度而不是速度。常见的错误是把 ma 和 mv 混在一起或者把平均力代入瞬时方程。在数值计算中我们通常把二阶微分方程降阶成一阶方程组用状态向量 [x, v] 来表示这样更容易交给求解器处理。2.3 单位制的转换SI 与工程单位制之间的桥国际单位制SI中质量kg、长度m、时间s是基本单位力是导出单位牛顿N1 N 1 kg·m/s²。工程单位制里则反过来力kgf、长度m、时间s是基本单位质量是导出单位。1 kgf 约等于 9.80665 N这个数值来源于标准重力加速度。换算公式如下单位制基本量力的单位质量单位换算关系国际单位制 (SI)质量 kg长度 m时间 sN (kg·m/s²)kg1 N 1 kg·m/s²工程单位制力 kgf长度 m时间 skgfkgf·s²/m1 kgf 9.80665 N换算实例—1 kgf ≈ 9.8 N1 工程质量单位 ≈ 9.8 kg—实际工程计算中如果把 kgf 直接当 N 用误差是 9.8 倍这在动力学仿真里会直接导致错误结果。我一般建议所有计算统一走 SI 制只在最终输出或图纸标注时再转换。用代码落地的话可以写一个小工具函数def kgf_to_newton(kgf, g9.80665): 工程单位制力到SI制力的转换 return kgf * g def emu_to_kg(emu): 工程单位质量(kgf·s^2/m)到kg的转换 return emu * 9.80665这段代码的逻辑很简单kgf 乘以标准重力加速度就是牛顿工程单位质量同理。注意这里的 g 只是单位换算用的常量不是当地重力加速度除非题目特别说明否则不要混用。3. 运动微分方程的三种坐标系形式选对坐标系才能少算三步3.1 矢量形式的通式质点动力学基本方程写成矢量形式是m * d²r/dt² ΣF这个式子最简单也最抽象。实际解题时必须把它投影到具体坐标系中。选择坐标系的原则是让受力或运动表达起来最简洁。例如直线运动用直角坐标圆周运动用自然坐标有心力问题用极坐标。3.2 直角坐标形式适合力是位置或时间的显函数直角坐标形式把矢量方程拆成三个标量方程m * d²x/dt² ΣFx m * d²y/dt² ΣFy m * d²z/dt² ΣFz每个方向独立求解互不耦合前提是力没有耦合项。这种形式适合抛体运动、滑块运动等场景。比如曲柄连杆机构里滑块 B 沿 x 轴运动只需要写 x 方向的方程即可。3.3 自然坐标形式切向与法向分开圆周运动利器自然坐标轴系沿着轨迹的切线τ、主法线n和副法线b方向。运动微分方程写为切向: m * dv/dt ΣFτ 法向: m * v²/ρ ΣFn 副法线: ΣFb 0这里 ρ 是轨迹的曲率半径。切向方程决定速度大小的变化法向方程决定方向的变化。对于匀速圆周运动v 是常量切向加速度为零法向方程直接给出向心力。圆锥摆就是典型的自然坐标应用场景。3.4 坐标系选择的工程判断问题特征首选坐标系原因直线轨道、平动直角坐标每个方向独立好解圆周运动、已知轨迹自然坐标切向/法向解耦有心力、极坐标问题极坐标角动量守恒更明显三维空间受任意力直角坐标数值积分通用性强易于程序化在实际编程时我通常先用直角坐标把状态方程写出来再根据约束条件做降维。自然坐标更适合手算推导但数值求解时需要额外判断曲率半径反而麻烦。所以例题里我用直角坐标做数值解用自然坐标做解析校验。4. 两类问题的数值解法从微积分到离散迭代4.1 第一类问题运动已知求力第一类问题的计算量不大核心是求导。已知 x(t)对它求二阶导数得到 a(t)再乘以质量得到力。工程上如果运动数据来自传感器采样直接微分会把噪声放大通常要先做平滑或滤波。常见做法是滑移窗口多项式拟合然后对拟合函数求导。下面是一个简单的例子import numpy as np from numpy.polynomial import polynomial as P def acceleration_from_position(t, x, window5, order2): 从离散位移信号估计加速度 方法局部多项式拟合后求二阶导 a np.zeros_like(x) half window // 2 for i in range(len(t)): lo max(0, i - half) hi min(len(t), i half 1) t_loc t[lo:hi] x_loc x[lo:hi] # 归一化时间避免数值病态 t0 t_loc.mean() t_norm (t_loc - t0) / (t_loc.std() 1e-12) coefs np.polynomial.polynomial.polyfit(t_norm, x_loc, order) # 对归一化时间求二阶导再换算到原始时间尺度 c2 coefs[2] * 2 a[i] c2 / (t_loc.std() 1e-12) ** 2 return a这段代码用局部窗口拟合二次多项式再解析求二阶导。参数window控制平滑程度窗口越大越平滑但延迟越大order一般取 2 或 3取高了容易过拟合噪声。注意时间归一化操作这是防止矩阵条件数过大的关键。4.2 第二类问题力求运动积分求解第二类问题需要解微分方程。对于线性系统可以用解析法对于非线性系统只能用数值积分。最常用的是 Runge-Kutta 法SciPy 的solve_ivp封装了 RK45 等算法。把二阶方程降阶为一阶方程组令 y1 x, y2 v dx/dt dy1/dt y2 dy2/dt F(t, y1, y2) / m下面以弹簧振子为例物块质量 m弹簧刚度 k初始拉长 a释放后求运动规律。微分方程为m * d²x/dt² -k * x边界条件x(0) a, v(0) 0。理论解是 x(t) a * cos(ω t)其中 ω sqrt(k/m)。from scipy.integrate import solve_ivp import matplotlib.pyplot as plt def spring_state(t, y, m, k): x, v y dxdt v dvdt -k * x / m return [dxdt, dvdt] m 1.0 # kg k 10.0 # N/m a 0.1 # m初始拉长量 t_span (0, 5) t_eval np.linspace(0, 5, 500) sol solve_ivp(spring_state, t_span, [a, 0], t_evalt_eval, args(m, k), rtol1e-8, atol1e-10) omega np.sqrt(k / m) x_analytic a * np.cos(omega * t_eval) x_numeric sol.y[0] error np.max(np.abs(x_numeric - x_analytic)) print(数值解与解析解最大误差:, error)solve_ivp的关键参数是rtol和atol它们控制相对误差和绝对误差容限。这里取了 1e-8 和 1e-10能得到机器精度级别的结果。如果误差不达标优先放宽atol而不是rtol因为位移量级可能很小。4.3 两类问题的对偶关系第一类是微分求导第二类是积分求解微分方程两者互为逆过程。数值上微分是病态问题积分是稳定问题。所以当你面对的是“由测量位移推受力”的逆动力学问题时结果往往对噪声极其敏感需要引入正则化。这是工程中经常被忽视的边界条件。5. 哈工大版第十章三个例题的完整复现5.1 例10-1曲柄连杆机构求连杆力滑块 B 的运动方程近似为x(t) r*(1 - λ²/4) r*cos(ωt) r*λ/2*cos(2ωt)其中 r 是曲柄半径ω 是曲柄角速度λ r/l 是曲柄与连杆长度之比。滑块质量 m忽略摩擦和连杆质量求某时刻连杆 AB 对滑块的力。注意连杆是二力杆力沿 AB 方向。为了求力需要对 x(t) 求二阶导得到 a_x然后根据牛顿第二定律求水平力再几何投影到杆方向。这里直接用符号微分是更稳的做法import sympy as sp t, r, l, m, omega, lam sp.symbols(t r l m omega lambda, positiveTrue) x r*(1 - lam**2/4) r*sp.cos(omega*t) r*lam/2*sp.cos(2*omega*t) v sp.diff(x, t) a sp.diff(v, t) print(加速度表达式:, sp.simplify(a))加速度表达式: -omega**2*r*(cos(omega*t) 2*lam*cos(2*omega*t))在 φ 0 时cos(0)1cos(0)1得到 a_x -ω² r (1 2λ)。代入数值 r0.1 ml0.5 mω10 rad/sm2 kg则 λ0.2a_x -1000.1(10.4) -14 m/s²水平力 F_x m*a_x -28 N。负号表示力与 x 正向相反。当 φ π/2 时cos(π/2)0cos(π)-1a_x -ω² r (0 - 2λ) 2 λ ω² r方向变为正。这就是为什么曲柄连杆机构中连杆有时受拉、有时受压。5.2 例10-2圆锥摆的速度与张力圆锥摆小球质量 m0.1 kg绳长 l0.3 m与铅直线夹角 θ60°在水平面内做匀速圆周运动。用自然坐标方程法向F sinθ m v² / ρ其中 ρ l sinθ切向F cosθ - mg 0由切向方程得 F mg / cosθ 0.1*9.8/0.5 1.96 N。再由法向方程得 v² F sinθ * ρ / m F sinθ * l sinθ / m F l sin²θ / m。代入数值import math m 0.1 g 9.8 l 0.3 theta math.radians(60) F m * g / math.cos(theta) rho l * math.sin(theta) v math.sqrt(F * math.sin(theta) * rho / m) print(f张力 F {F:.3f} N, 速度 v {v:.3f} m/s)张力 F 1.960 N, 速度 v 2.100 m/s这个例子说明自然坐标下匀速圆周运动的切向方程退化为恒等式法向方程直接给出向心力来源。如果用力学单位制结果会差 9.8 倍这是最常见的低级错误。5.3 例10-3弹簧振子的数值解与能量验证前面已经用 solve_ivp 解过弹簧振子。现在验证数值解的长期稳定性。理论角频率 ωsqrt(k/m)周期 T2π/ω。取 k10, m1, a0.1数值解 5 秒内误差保持在 1e-10 量级。再多算一步总能量 E 1/2 m v² 1/2 k x²应该守恒v_num sol.y[1] E 0.5 * m * v_num**2 0.5 * k * sol.y[0]**2 E0 0.5 * k * a**2 print(能量最大相对误差:, np.max(np.abs(E - E0) / E0))能量最大相对误差: 2.2e-10对于无阻尼系统能量守恒是检验数值积分器性能的重要指标。如果误差随时间线性增长说明步长过大或方法不稳定。RK45 自适应步长在大多数情况下足够用但更长的仿真时间建议改用更高阶方法或 symplectic 积分器比如scipy.integrate.solve_ivp中的Radau不过它并不保能量。工程上如果跑轨道仿真建议用 split-operator 或 leapfrog 方法避免能量漂移。6. 进阶技巧用自然坐标验证数值解用能量判据检查积分器数值解算完之后怎么确认没有出错至少三重验证解析特解对照、物理不变量守恒、极限情况退化。这三件事都做一遍比单纯盯误差曲线更可靠。第一招退化为已知特解。把曲柄连杆的 λ 设成 0此时 l 无穷大滑块做简谐运动 x r r cos(ωt)加速度 a -ω² r cos(ωt)与简谐振动一致。用上面的代码把 lam0 跑一遍看数值结果是否匹配。如果匹配说明代码框架没问题。第二招自然坐标解析验证。圆锥摆的数值解可以按直角坐标写但手算解析解用自然坐标更方便。我做完数值解后会特意算一下切向加速度是否为零把速度分解到轨迹切线方向求导看是否为 0。如果非零说明稳态假设不成立或者角度选错了这往往是坐标系选择用力过猛导致的。第三招能量或动量守恒检验。对于保守系统能量应该守恒。上面弹簧振子的能量误差已经到 1e-10 量级说明积分器没问题。如果跑的是有阻尼或阻尼可变系统则改用动量守恒或者阻尼功耗散率对照理论公式。关于时间步长我的经验是对于快速振荡问题周期 T初始步长建议小于 T/100对于慢变问题的长时间仿真用自适应步长时最好限制最大步长防止在加速度变化快的区间计算次数不足。如果你发现数值解发散第一步不是减小步长而是检查状态方程是否写错比如符号颠倒、质量用错或者单位制没统一。第二步才调整积分容差。最后补充一个实用技巧把上面的代码封装成一个particle_solver类支持任意外力函数这样你可以复用于不同的动力学问题不必每道例题重新写一遍积分逻辑。统一的接口可以大幅减少第二类问题反复调试带来的时间成本。本文还有配套的精品资源点击获取 理论力学质点动力学：三类典型问题的Python数值解法
简介一份哈尔滨工业大学版《理论力学》教学课件面向工科院校学生与工程技术人员适合在学完静力学、运动学后系统学习动力学核心内容。资源聚焦第十章“质点动力学”以牛顿三定律为基础详细梳理动力学基本定律、国际单位制与工程单位制并重点讲解质点运动微分方程的矢量形式、直角坐标形式和自然坐标形式进而归纳质点动力学的两类基本问题已知运动求力、已知力求运动。课件结合曲柄连杆机构经典例题演示如何选取研究对象、建立运动方程并利用微分积分求解帮助读者掌握动力学分析的一般方法。全包内含1个演示文稿文件共1010KB结构紧凑、要点突出可用于课堂同步学习、课后复习或考研一轮梳理。该资源已有5626人学习下载是深入理解力与运动关系、夯实理论力学基础的实用资料。1. 质点动力学从静力平衡到运动变化的那一步大部分工程师接触力学是从静力学开始的一个物体放在那里受力平衡求出约束反力就算完事。到了运动学阶段我们开始描述位置、速度和加速度但完全不问“是什么力让它变成这样”。真正把这两者接起来的是动力学。而动力学里最先遇到的就是质点动力学——把物体抽象成一个有质量但没有尺寸的点研究它受的力与运动变化之间的关系。哈工大版的《理论力学》第十章就是围绕这个核心展开的三个基本定律、运动微分方程、以及两类典型问题。这里有个容易被忽略但很关键的点第一类问题是“已知运动求力”只需要微分第二类问题是“已知力求运动”需要积分。前者是正问题后者是反问题。实际工程里大多数情况是已知激励求响应也就是第二类但第一类在传感器数据处理、逆动力学分析里同样常见。本文从定律和单位制讲起把三类坐标系下的微分方程写清楚然后给出曲柄连杆、圆锥摆、弹簧振子三个例题的Python数值解法最后补一点关于数值积分稳定性和精度验证的实操经验。2. 三大定律与单位制先搞清楚Fma里每个符号的物理边界2.1 牛顿第一定律的真正用途定义惯性参考系牛顿第一定律惯性定律描述的是质点不受力或合力为零时它对惯性参考系保持静止或匀速直线运动。很多初学者把它当作第二定律的推论实际上它在逻辑上更靠前——它定义了“惯性参考系”这个概念。没有惯性参考系第二定律中的加速度就没有可对比的基准。地面通常被近似为惯性系但在涉及地球自转的精密计算里科氏力就必须被考虑进来这时地面不再是严格的惯性系。2.2 牛顿第二定律的矢量性与瞬时性第二定律的表达式是 F ma但要注意两点第一这是矢量方程F 和 a 方向相同第二它是瞬时关系某一时刻的力决定该时刻的加速度而不是速度。常见的错误是把 ma 和 mv 混在一起或者把平均力代入瞬时方程。在数值计算中我们通常把二阶微分方程降阶成一阶方程组用状态向量 [x, v] 来表示这样更容易交给求解器处理。2.3 单位制的转换SI 与工程单位制之间的桥国际单位制SI中质量kg、长度m、时间s是基本单位力是导出单位牛顿N1 N 1 kg·m/s²。工程单位制里则反过来力kgf、长度m、时间s是基本单位质量是导出单位。1 kgf 约等于 9.80665 N这个数值来源于标准重力加速度。换算公式如下单位制基本量力的单位质量单位换算关系国际单位制 (SI)质量 kg长度 m时间 sN (kg·m/s²)kg1 N 1 kg·m/s²工程单位制力 kgf长度 m时间 skgfkgf·s²/m1 kgf 9.80665 N换算实例—1 kgf ≈ 9.8 N1 工程质量单位 ≈ 9.8 kg—实际工程计算中如果把 kgf 直接当 N 用误差是 9.8 倍这在动力学仿真里会直接导致错误结果。我一般建议所有计算统一走 SI 制只在最终输出或图纸标注时再转换。用代码落地的话可以写一个小工具函数def kgf_to_newton(kgf, g9.80665): 工程单位制力到SI制力的转换 return kgf * g def emu_to_kg(emu): 工程单位质量(kgf·s^2/m)到kg的转换 return emu * 9.80665这段代码的逻辑很简单kgf 乘以标准重力加速度就是牛顿工程单位质量同理。注意这里的 g 只是单位换算用的常量不是当地重力加速度除非题目特别说明否则不要混用。3. 运动微分方程的三种坐标系形式选对坐标系才能少算三步3.1 矢量形式的通式质点动力学基本方程写成矢量形式是m * d²r/dt² ΣF这个式子最简单也最抽象。实际解题时必须把它投影到具体坐标系中。选择坐标系的原则是让受力或运动表达起来最简洁。例如直线运动用直角坐标圆周运动用自然坐标有心力问题用极坐标。3.2 直角坐标形式适合力是位置或时间的显函数直角坐标形式把矢量方程拆成三个标量方程m * d²x/dt² ΣFx m * d²y/dt² ΣFy m * d²z/dt² ΣFz每个方向独立求解互不耦合前提是力没有耦合项。这种形式适合抛体运动、滑块运动等场景。比如曲柄连杆机构里滑块 B 沿 x 轴运动只需要写 x 方向的方程即可。3.3 自然坐标形式切向与法向分开圆周运动利器自然坐标轴系沿着轨迹的切线τ、主法线n和副法线b方向。运动微分方程写为切向: m * dv/dt ΣFτ 法向: m * v²/ρ ΣFn 副法线: ΣFb 0这里 ρ 是轨迹的曲率半径。切向方程决定速度大小的变化法向方程决定方向的变化。对于匀速圆周运动v 是常量切向加速度为零法向方程直接给出向心力。圆锥摆就是典型的自然坐标应用场景。3.4 坐标系选择的工程判断问题特征首选坐标系原因直线轨道、平动直角坐标每个方向独立好解圆周运动、已知轨迹自然坐标切向/法向解耦有心力、极坐标问题极坐标角动量守恒更明显三维空间受任意力直角坐标数值积分通用性强易于程序化在实际编程时我通常先用直角坐标把状态方程写出来再根据约束条件做降维。自然坐标更适合手算推导但数值求解时需要额外判断曲率半径反而麻烦。所以例题里我用直角坐标做数值解用自然坐标做解析校验。4. 两类问题的数值解法从微积分到离散迭代4.1 第一类问题运动已知求力第一类问题的计算量不大核心是求导。已知 x(t)对它求二阶导数得到 a(t)再乘以质量得到力。工程上如果运动数据来自传感器采样直接微分会把噪声放大通常要先做平滑或滤波。常见做法是滑移窗口多项式拟合然后对拟合函数求导。下面是一个简单的例子import numpy as np from numpy.polynomial import polynomial as P def acceleration_from_position(t, x, window5, order2): 从离散位移信号估计加速度 方法局部多项式拟合后求二阶导 a np.zeros_like(x) half window // 2 for i in range(len(t)): lo max(0, i - half) hi min(len(t), i half 1) t_loc t[lo:hi] x_loc x[lo:hi] # 归一化时间避免数值病态 t0 t_loc.mean() t_norm (t_loc - t0) / (t_loc.std() 1e-12) coefs np.polynomial.polynomial.polyfit(t_norm, x_loc, order) # 对归一化时间求二阶导再换算到原始时间尺度 c2 coefs[2] * 2 a[i] c2 / (t_loc.std() 1e-12) ** 2 return a这段代码用局部窗口拟合二次多项式再解析求二阶导。参数window控制平滑程度窗口越大越平滑但延迟越大order一般取 2 或 3取高了容易过拟合噪声。注意时间归一化操作这是防止矩阵条件数过大的关键。4.2 第二类问题力求运动积分求解第二类问题需要解微分方程。对于线性系统可以用解析法对于非线性系统只能用数值积分。最常用的是 Runge-Kutta 法SciPy 的solve_ivp封装了 RK45 等算法。把二阶方程降阶为一阶方程组令 y1 x, y2 v dx/dt dy1/dt y2 dy2/dt F(t, y1, y2) / m下面以弹簧振子为例物块质量 m弹簧刚度 k初始拉长 a释放后求运动规律。微分方程为m * d²x/dt² -k * x边界条件x(0) a, v(0) 0。理论解是 x(t) a * cos(ω t)其中 ω sqrt(k/m)。from scipy.integrate import solve_ivp import matplotlib.pyplot as plt def spring_state(t, y, m, k): x, v y dxdt v dvdt -k * x / m return [dxdt, dvdt] m 1.0 # kg k 10.0 # N/m a 0.1 # m初始拉长量 t_span (0, 5) t_eval np.linspace(0, 5, 500) sol solve_ivp(spring_state, t_span, [a, 0], t_evalt_eval, args(m, k), rtol1e-8, atol1e-10) omega np.sqrt(k / m) x_analytic a * np.cos(omega * t_eval) x_numeric sol.y[0] error np.max(np.abs(x_numeric - x_analytic)) print(数值解与解析解最大误差:, error)solve_ivp的关键参数是rtol和atol它们控制相对误差和绝对误差容限。这里取了 1e-8 和 1e-10能得到机器精度级别的结果。如果误差不达标优先放宽atol而不是rtol因为位移量级可能很小。4.3 两类问题的对偶关系第一类是微分求导第二类是积分求解微分方程两者互为逆过程。数值上微分是病态问题积分是稳定问题。所以当你面对的是“由测量位移推受力”的逆动力学问题时结果往往对噪声极其敏感需要引入正则化。这是工程中经常被忽视的边界条件。5. 哈工大版第十章三个例题的完整复现5.1 例10-1曲柄连杆机构求连杆力滑块 B 的运动方程近似为x(t) r*(1 - λ²/4) r*cos(ωt) r*λ/2*cos(2ωt)其中 r 是曲柄半径ω 是曲柄角速度λ r/l 是曲柄与连杆长度之比。滑块质量 m忽略摩擦和连杆质量求某时刻连杆 AB 对滑块的力。注意连杆是二力杆力沿 AB 方向。为了求力需要对 x(t) 求二阶导得到 a_x然后根据牛顿第二定律求水平力再几何投影到杆方向。这里直接用符号微分是更稳的做法import sympy as sp t, r, l, m, omega, lam sp.symbols(t r l m omega lambda, positiveTrue) x r*(1 - lam**2/4) r*sp.cos(omega*t) r*lam/2*sp.cos(2*omega*t) v sp.diff(x, t) a sp.diff(v, t) print(加速度表达式:, sp.simplify(a))加速度表达式: -omega**2*r*(cos(omega*t) 2*lam*cos(2*omega*t))在 φ 0 时cos(0)1cos(0)1得到 a_x -ω² r (1 2λ)。代入数值 r0.1 ml0.5 mω10 rad/sm2 kg则 λ0.2a_x -1000.1(10.4) -14 m/s²水平力 F_x m*a_x -28 N。负号表示力与 x 正向相反。当 φ π/2 时cos(π/2)0cos(π)-1a_x -ω² r (0 - 2λ) 2 λ ω² r方向变为正。这就是为什么曲柄连杆机构中连杆有时受拉、有时受压。5.2 例10-2圆锥摆的速度与张力圆锥摆小球质量 m0.1 kg绳长 l0.3 m与铅直线夹角 θ60°在水平面内做匀速圆周运动。用自然坐标方程法向F sinθ m v² / ρ其中 ρ l sinθ切向F cosθ - mg 0由切向方程得 F mg / cosθ 0.1*9.8/0.5 1.96 N。再由法向方程得 v² F sinθ * ρ / m F sinθ * l sinθ / m F l sin²θ / m。代入数值import math m 0.1 g 9.8 l 0.3 theta math.radians(60) F m * g / math.cos(theta) rho l * math.sin(theta) v math.sqrt(F * math.sin(theta) * rho / m) print(f张力 F {F:.3f} N, 速度 v {v:.3f} m/s)张力 F 1.960 N, 速度 v 2.100 m/s这个例子说明自然坐标下匀速圆周运动的切向方程退化为恒等式法向方程直接给出向心力来源。如果用力学单位制结果会差 9.8 倍这是最常见的低级错误。5.3 例10-3弹簧振子的数值解与能量验证前面已经用 solve_ivp 解过弹簧振子。现在验证数值解的长期稳定性。理论角频率 ωsqrt(k/m)周期 T2π/ω。取 k10, m1, a0.1数值解 5 秒内误差保持在 1e-10 量级。再多算一步总能量 E 1/2 m v² 1/2 k x²应该守恒v_num sol.y[1] E 0.5 * m * v_num**2 0.5 * k * sol.y[0]**2 E0 0.5 * k * a**2 print(能量最大相对误差:, np.max(np.abs(E - E0) / E0))能量最大相对误差: 2.2e-10对于无阻尼系统能量守恒是检验数值积分器性能的重要指标。如果误差随时间线性增长说明步长过大或方法不稳定。RK45 自适应步长在大多数情况下足够用但更长的仿真时间建议改用更高阶方法或 symplectic 积分器比如scipy.integrate.solve_ivp中的Radau不过它并不保能量。工程上如果跑轨道仿真建议用 split-operator 或 leapfrog 方法避免能量漂移。6. 进阶技巧用自然坐标验证数值解用能量判据检查积分器数值解算完之后怎么确认没有出错至少三重验证解析特解对照、物理不变量守恒、极限情况退化。这三件事都做一遍比单纯盯误差曲线更可靠。第一招退化为已知特解。把曲柄连杆的 λ 设成 0此时 l 无穷大滑块做简谐运动 x r r cos(ωt)加速度 a -ω² r cos(ωt)与简谐振动一致。用上面的代码把 lam0 跑一遍看数值结果是否匹配。如果匹配说明代码框架没问题。第二招自然坐标解析验证。圆锥摆的数值解可以按直角坐标写但手算解析解用自然坐标更方便。我做完数值解后会特意算一下切向加速度是否为零把速度分解到轨迹切线方向求导看是否为 0。如果非零说明稳态假设不成立或者角度选错了这往往是坐标系选择用力过猛导致的。第三招能量或动量守恒检验。对于保守系统能量应该守恒。上面弹簧振子的能量误差已经到 1e-10 量级说明积分器没问题。如果跑的是有阻尼或阻尼可变系统则改用动量守恒或者阻尼功耗散率对照理论公式。关于时间步长我的经验是对于快速振荡问题周期 T初始步长建议小于 T/100对于慢变问题的长时间仿真用自适应步长时最好限制最大步长防止在加速度变化快的区间计算次数不足。如果你发现数值解发散第一步不是减小步长而是检查状态方程是否写错比如符号颠倒、质量用错或者单位制没统一。第二步才调整积分容差。最后补充一个实用技巧把上面的代码封装成一个particle_solver类支持任意外力函数这样你可以复用于不同的动力学问题不必每道例题重新写一遍积分逻辑。统一的接口可以大幅减少第二类问题反复调试带来的时间成本。本文还有配套的精品资源点击获取