使用Odeint求解耦合微分方程结果异常，求代码错误排查

我为物理项目求解标枪的耦合微分方程，编写了Python代码能运行但结果完全错误，已检查过所用系数，求排查代码问题。参考文档为标枪运动相关论文的第5/6页。

原代码 #

#résolution vectorielle des 2 équations couplées p5 et p6 par la méthode de Runge Kutta  ordre 4 (?) (odeint)

#importation des bibliothèques
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Mise de l'équation sous forme dY/dt = f(Y,t) (cf classeur TIPE)
# Y = [[u],[v]] , u = dr / dt et v = ds /dt
# odeint prend en argument (la fonction f, les conditions initiales, )

#? Calcul des coefficients
epsilon = 30e-3 #rayon du javelot en m
gamma = 2 #Cx amimensionné
rho = 1.204 #masse volumique de l'air en kg/m^3
L = 2700e-3 #longueur javelot en m
m = 800e-3 #masse javelot en kg
delta = 10*(np.pi/180) #angle d'attaque en rad
alpha = np.pi/4 #angle de lancé en rad
g = 9.81

betaj = 8*gamma*rho*epsilon*L  #8gamma*rho*epsilon*L  de dim M.L^-1
gammaj = 2*gamma*np.pi*rho*np.square(epsilon) #en M.L^-1

a = -(gammaj + betaj*np.abs(np.sin(delta)))/m
b = -g + (betaj/m) * np.sin(delta)*np.sin(alpha)



# définition de f(Y,t)
def SecondMembre(Y,t) : #u vitesse selon ux, v vitesse selon uy
    u = Y[0]
    v = Y[1]
    dudt = a*u**2 + a*v**2
    dvdt = b*u**2 + b*v**2
    return([dudt, dvdt])

#tableau de temps en lesquels la fonction fournit la solution Y(t)
    #n est le nombre de temps dans la liste, ti le temps initial et tf le temps final en seconde
ti=0
tf=30
n=200
t=np.linspace(ti,tf,n)  
 
#conditions initiales
u0 = 7.07
v0 = 7.07
CondInit=np.array([u0,v0])

#utilisation de odeint
Ys = odeint(SecondMembre, CondInit, t)  #renvoie une matrice avec u et v  pour les 200 pas de temps
u=Ys[:,0]
v=Ys[:,1]  #! c'est une fonction numpy et pas un tableau python classique

print(u)

代码核心错误分析 #

1. 微分方程实现完全偏离物理模型
   参考文档中的标枪运动方程，空气阻力和升力产生的加速度分量与速度分量和速度模长的乘积成正比，而非你代码中直接使用的速度平方和。

   • 正确的速度模长应为 V = np.sqrt(u² + v²)
   • 加速度分量需结合速度分量与模长计算，而非a*(u²+v²)这类错误形式。

2. 系数定义混淆
   你定义的b错误地将重力项与升力项合并，文档中这两项是独立的，且升力项需要分别作用在水平和竖直速度分量上，需要拆分出独立的系数处理交叉项。

3. 积分逻辑不完整
   当前代码仅积分了速度，未计算位移；且积分时间设置为30秒，远超过标枪落地时间，会产生无效计算。

修正后的代码 #

# 标枪耦合运动微分方程求解（Runge-Kutta 4阶，odeint实现）
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 参数定义
epsilon = 30e-3       # 标枪半径(m)
gamma = 2             # 无量纲阻力系数
rho = 1.204           # 空气密度(kg/m³)
L = 2700e-3           # 标枪长度(m)
m = 800e-3            # 标枪质量(kg)
delta = 10*(np.pi/180)# 攻角(rad)
alpha = np.pi/4       # 发射角(rad)
g = 9.81              # 重力加速度(m/s²)

# 计算核心系数
betaj = 8 * gamma * rho * epsilon * L
gammaj = 2 * gamma * np.pi * rho * np.square(epsilon)

# 拆分阻力、升力相关系数
coeff_drag = -(gammaj + betaj * np.abs(np.sin(delta))) / m
coeff_lift_u = (betaj * np.sin(delta) * np.cos(alpha)) / m
coeff_lift_v = (betaj * np.sin(delta) * np.sin(alpha)) / m

# 定义微分方程组：Y = [水平速度u, 竖直速度v, 水平位移x, 竖直位移y]
def SecondMembre(Y, t):
    u, v, x, y = Y
    V = np.sqrt(u**2 + v**2)
    # 防护速度为0的极端情况
    if V < 1e-6:
        V = 1e-6
    # 参考文档中的加速度方程
    dudt = coeff_drag * u * V + coeff_lift_u * v * V
    dvdt = -g + coeff_drag * v * V + coeff_lift_v * u * V
    dxdt = u
    dydt = v
    return [dudt, dvdt, dxdt, dydt]

# 时间数组：缩短至标枪合理飞行时间
ti = 0
tf = 5
n = 200
t = np.linspace(ti, tf, n)

# 初始条件：初速度10m/s（45度分解为u0=v0≈7.07），初始位置原点
u0 = 7.07
v0 = 7.07
x0 = 0
y0 = 0
CondInit = np.array([u0, v0, x0, y0])

# 求解微分方程
Ys = odeint(SecondMembre, CondInit, t)
u = Ys[:, 0]
v = Ys[:, 1]
x = Ys[:, 2]
y = Ys[:, 3]

# 绘制运动轨迹
plt.figure(figsize=(10, 6))
plt.plot(x, y)
plt.xlabel('水平位移(m)')
plt.ylabel('竖直位移(m)')
plt.title('标枪运动轨迹')
plt.grid(True)
plt.show()

# 打印前10个时间步的速度数据
print("水平速度u（前10步）:", u[:10])
print("竖直速度v（前10步）:", v[:10])

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