4.2 运动微分方程解析与数值求解

4.2 运动微分方程解析与数值求解

本节摘要：受力分析五步走完，交到你手里的是一个微分方程——未知的是"位置随时间怎么变"这条函数曲线。本节讲两条求解路线：阻力与弹簧这类少数幸运情形有解析解，要什么给什么；多数真实情形只能数值积分，欧拉法教思想，龙格-库塔法干活，步长与精度的取舍本节算给你看。

解析是特例，数值是通例

上一节的五步流程把力逐根标清、投影方程立好，账目右端出现 ma，而 a 是速度对时间的导数——方程从此带上了微分号，解的对象从"数"升级为"函数"。这一步升级决定了我方武器的分工：少数情形（线性阻力、小角度摆、弹簧振子）方程能用手笔积出来，解是干净的公式，参数一改答案随手重算；多数情形（平方阻力、非线性弹簧、摆角一大）积分积不动，就得让机器一小步一小步往前挪。两条路线不是竞争关系，而是上下游关系——解析解给数值解当裁判，数值解给解析解够不着的病例动刀。本节承接上一节的方程产出，向下为第六章的振动仿真与第八章的计算力学铺路。

路线一：解析积分，以线性阻力落体为例

拿一台标准病例演示：质量一千克的雨滴下落，空气阻力与速率成正比，比例系数取每秒零点四千克。五步流程照样走完：隔离雨滴，标重力竖直向下、阻力竖直向上，方程写成质量乘加速度等于重力减阻力项。分离变量积分，配上初速为零的初始条件，速度随时间的公式与收尾速度一并到手。

# 解析路线：线性阻力落体 m*dv/dt = m*g - b*v，初速为零
m, g, b = 1.0, 9.8, 0.4
v_terminal = m * g / b                 # 收尾速度：阻力扛平重力
tau = m / b                            # 时间常数
t = 2.0
v = v_terminal * (1 - pow(2.718281828, -t / tau))
import math
v_exact = v_terminal * (1 - math.exp(-t / tau))
print("收尾速度 v_t =", v_terminal, "m/s")     # 24.5
print("时间常数 tau =", tau, "s")              # 2.5
print("t=2s 时 v =", round(v_exact, 4), "m/s") # 13.4914
# 物理读数：一个时间常数后达到收尾速度的 63.2%
print("一个 tau 后占比 =", 1 - math.exp(-1))   # 0.6321205588
# 检查 5 个 tau：工程上视为"已经收尾"
print("5 tau 后占比 =", 1 - math.exp(-5))      # 0.993262053
​

公式给出的读数很有分量：雨滴不用无限久才"不加速"，约五个时间常数（这里十二秒半）后速度已达收尾值的九成九以上。小雾滴的 b 比值更大，收尾速度更低——这就是雾悬浮而雨落地的原因，一条解析式把两种天气现象串在了一根线上。

路线二：数值积分，欧拉与龙格-库塔同台

换成平方阻力（阻力正比于速率平方），分离变量积不动了，数值方法登场。欧拉法是最朴素的思路：用当前时刻的斜率往前挪一步。它便宜但误差随步长线性累积；龙格-库塔四阶法在每个步长内采样四次斜率做加权平均，误差随步长的四次方缩小——步长砍半，欧拉误差减半，龙格-库塔误差缩到约十六分之一。

# 数值路线：平方阻力落体 dv/dt = g - (c/m)*v*|v|，初速 20 m/s 向下
import math
m, g, c = 80.0, 9.8, 0.25          # 跳伞员 80 kg，阻力系数 0.25 kg/m
vt = math.sqrt(m * g / c)          # 收尾速度
f = lambda t, v: g - (c / m) * v * abs(v)

def euler(dt, n):
    t, v = 0.0, 20.0
    for _ in range(n):
        v += f(t, v) * dt
        t += dt
    return t, v

def rk4(dt, n):
    t, v = 0.0, 20.0
    for _ in range(n):
        k1 = f(t, v)
        k2 = f(t + dt/2, v + dt/2*k1)
        k3 = f(t + dt/2, v + dt/2*k2)
        k4 = f(t + dt, v + dt*k3)
        v += dt * (k1 + 2*k2 + 2*k3 + k4) / 6
        t += dt
    return t, v

print("收尾速度 =", round(vt, 4), "m/s")          # 56.0
for dt in (0.1, 0.05):
    te, ve = euler(dt, int(2.0/dt))
    tr, vr = rk4(dt, int(2.0/dt))
    print(f"dt={dt}: 欧拉 v={ve:.6f}  龙格库塔 v={vr:.6f}")
# dt=0.1: 欧拉 v=34.777404  龙格库塔 v=34.671768
# dt=0.05: 欧拉 v=34.724388  龙格库塔 v=34.671768
​

同一时刻同一状态，龙格-库塔在两种步长下给出六位小数一致的答案，欧拉法却在第一位小数之后就分了家（差约零点一），步长砍半也只把差距追回一半——这就是"误差阶数"的直观含义。工程仿真几乎清一色用高阶法，欧拉法留作理解数值积分的教学起点。

⚠️ 常见坑：数值解不是"算出来就对了"。步长太大时龙格-库塔也会给出看似平滑、实则畸变的曲线；换一个步长重跑一遍，两次结果在小数位上咬合，才敢签字。第八章把这套检验固化成极限检验清单。

再进一步：刚性问题与隐式思想

还有一类方程会让显式方法吃暗亏：各项时间尺度悬殊巨大的"刚性"问题。典型如强阻尼环节，衰减的时间常数比运动过程短几个量级，显式法的步长被最快的那把尺子卡死——步长稍大，数值解不但不准，还会指数式放大直接爆炸。隐式方法换了一个思路：斜率不取当前点，改取未知终点处的斜率，每步多解一个小方程，换来"想走多大步就走多大步"的稳定性。结构动力学与化学反应动力学的商用软件默认隐式，原因就在这里。

# 刚性问题演示：dv/dt = -20*v，显式欧拉 vs 隐式欧拉，dt=0.12
v0, lam, dt, n = 1.0, 20.0, 0.12, 5
v_ex, v_im = v0, v0
for _ in range(n):
    v_ex = v_ex * (1 - lam * dt)          # 显式：因子 1-lam*dt = -1.4，绝对值大于 1，振荡发散
    v_im = v_im / (1 + lam * dt)          # 隐式：因子 1/(1+lam*dt)，恒小于 1，稳定
print("显式五步 =", [round(v0 * (1-lam*dt)**i, 4) for i in range(6)])
# 显式五步 = [1.0, -1.4, 1.96, -2.744, 3.8416, -5.3782]  —— 振荡放大，炸了
print("隐式五步 =", [round(v0 / (1+lam*dt)**i, 4) for i in range(6)])
# 隐式五步 = [1.0, 0.2941, 0.0865, 0.0254, 0.0075, 0.0022]  —— 平滑衰减，物理正确
​

同一道方程，两套步法，一套发散一套守规——刚性问题的第一课就是认出它，再换武器。

要点复盘：

• 五步流程的产出是微分方程：解的对象从数升级为函数，求解路线按方程性质分流；
• 解析解是裁判：线性阻力落体的收尾速度与时间常数给出"多快不加速"的定量答案；
• 误差阶数决定取舍：步长砍半，欧拉误差减半、龙格-库塔缩到约十六分之一；
• 数值解要自检：变步长重跑对比，是确认曲线可信的最低成本手段。

方程立得住、也解得动了。下一节讨论一个更省力的问题：有些题根本不必知道全过程，始末状态两条账目就能定案——动量与动能两本近路账本登场。