含空气阻力自由下落ODE数值解与解析解差异过大问题排查

核心错误：解析解与模型不匹配 #

你编写的数值求解代码针对的是平方空气阻力模型，运动方程为：

m*a = m*g -k*v**2  # 其中a=vdot

对应的ODE定义也正确：

def f(v, g, k, m):
    return g - k/m * v**2

但你用来对比的解析解却是线性空气阻力（阻力与速度一次方成正比）下的公式：

v_ana = m*g*(1.-np.exp(-k/m*t))/k

这两种阻力模型的物理规律完全不同，因此计算结果必然出现巨大差异。

图表显示：数值解收敛到较低的终端速度，而误用的解析解持续上升到一个高得多的速度，二者差异显著。

平方阻力下的正确解析解 #

对于平方阻力的自由下落问题，速度的解析解为：
$$v(t) = v_t \tanh\left(\frac{g t}{v_t}\right)$$
其中终端速度 $v_t = \sqrt{\frac{mg}{k}}$，代入你的参数（$g=9.81, k=0.1, m=1$），可得终端速度 $v_t=\sqrt{9.81/0.1}\approx9.90\ \text{m/s}$，你的数值解最终会收敛到这个值，求解逻辑是正确的。

数值求解的小优化建议 #

• 位移更新逻辑：你当前用 $x_{i+1}=x_i - v_{i+1} \cdot h$，若要更贴合显式欧拉的前向更新逻辑，也可以改为 $x_{i+1}=x_i - v_i \cdot h$，步长较小时二者差异极小；
• 时间范围优化：平方阻力下速度会在几秒内就接近终端速度，无需设置到1000秒，可缩短时间范围提升计算效率。

内容的提问来源于stack exchange，提问作者Axel