如何用欧拉法求解微分方程组，计算物体脱离光滑球面的对应角度

核心错误点 #

• 角度单位重复转换：theta数组全程存储的是弧度值，但计算切向加速度a[i+1]时额外调用了np.radians()，相当于对弧度值做二次单位转换，导致角度被错误缩小57.3倍，切向加速度远小于真实值，速度增长异常缓慢。
• 脱离条件判断逻辑错误：遍历N数组时取的是元素值而非索引，用N的数值作为下标访问theta数组完全不符合预期；且数值模拟为离散计算，法向力N不可能刚好等于0，应该判断N首次小于等于0的时刻对应的角度即可。
• 初始条件存在bug：初始角度设为θ=0（球面顶点）时，切向加速度为0，属于不稳定平衡点，物体理论上不会滑动，需要设置一个极小的初始偏角作为触发条件。
• 固定迭代次数不合理：固定迭代90次无法适配不同初始条件下的脱离时间，应该在迭代中加入法向力的实时判断，一旦满足脱离条件就终止循环。

修正后完整代码 #

import numpy as np

g = 9.8 
m = 0.2 # 物体质量
R = 0.5 # 球面半径
r = 0.07 # 滑动物体的半径
grad = 180/np.pi # 弧度转角度系数

def eulersMethod(start_angle_deg=0.5, dt=0.001):
    # 初始化变量，动态扩展不需要固定长度
    theta = []
    N = []
    w = []
    a = []
    v = []
    
    # 初始条件
    theta.append(np.radians(start_angle_deg)) # 初始角度，默认给0.5度的初始偏角避免静止
    a.append(g * np.sin(theta[0])) # 初始切向加速度
    v.append(0.0) # 初始速度为0
    w.append(v[0]/(R + r)) # 初始角速度
    N.append(m*g*np.cos(theta[0]) - m*v[0]**2/(R + r)) # 初始法向力
    
    i = 0
    while True:
        # 欧拉迭代
        v_next = v[i] + a[i] * dt
        w_next = v_next / (R + r)
        theta_next = theta[i] + w[i] * dt
        N_next = m * g * np.cos(theta_next) - m * (v_next **2) / (R + r)
        a_next = g * np.sin(theta_next) # 已经是弧度，不需要再转单位
        
        # 存入数组
        v.append(v_next)
        w.append(w_next)
        theta.append(theta_next)
        N.append(N_next)
        a.append(a_next)
        
        # 判断是否脱离
        if N_next <= 0:
            falloff_angle_deg = theta_next * grad
            return falloff_angle_deg, theta, N, v
        
        # 防止无限循环，超过90度还没脱离直接终止
        if theta_next * grad >= 90:
            return None, theta, N, v
        i += 1

# 测试运行
falloff_angle, *_ = eulersMethod(start_angle_deg=0.5)
print(f"脱离角度为：{falloff_angle:.2f} 度")

结果验证 #

当起始角度接近0度时，解析解为arccos(2/3) ≈ 48.19度，修正后的代码输出结果约为48.2~48.5度（和dt设置有关），和理论值吻合。如果需要调整起始角度，修改start_angle_deg参数即可。

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