轨迹模型没有给出正确的结果

Question

我在求解非线性常微分方程组时遇到了困难。他们似乎没有给我正确的结果。你知道吗

下面是我需要解的方程：

我想得到的结果是：

以下是我的输出：

“”“

def model(t, z):
    #Variables
    r, zeta, phi, V, Gamma, psi = z

    print(t)

    #Simplify
    rho = planet.density(r)
    gr = planet.g(r)
    w = planet.angular_rotation_rate

    #Postition
    drdt = V*np.sin(Gamma)
    dzetadt = (V*np.cos(Gamma)*np.cos(psi))/(r*np.cos(phi))
    dphidt = (V*np.cos(Gamma)*np.sin(psi))/r

    #Forces
    dVdt = -1*(rho*V**2)/(2*beta) - gr*np.sin(Gamma) + r*w**2*np.cos(phi)*(np.sin(Gamma)*np.cos(phi) - np.cos(Gamma)*np.sin(phi)*np.sin(psi))
    dGammadt = (-1*gr*np.cos(Gamma) + V**2*np.cos(Gamma)/(r) + 2*w*V*np.cos(phi)*np.cos(psi) + r*w**2*np.cos(phi)*(np.cos(Gamma)*np.cos(phi) + np.sin(Gamma)*np.sin(phi)*np.sin(psi)))/(V)
    dpsidt = (-1*(V**2*np.cos(Gamma)*np.cos(psi)*np.tan(phi))/r + 2*w*V*(np.tan(Gamma)*np.cos(phi)*np.sin(psi) - np.sin(phi)) - (r*w**2*np.sin(phi)*np.cos(phi)*np.cos(psi))/(np.cos(Gamma)))/V

    #Pack into dzdt
    dzdt = [drdt, dzetadt, dphidt, dVdt, dGammadt, dpsidt]
    return dzdt

def hit_ground(t, y): 
    return y[0]
hit_ground.terminal = True
hit_ground.direction = -1

#Solv    print("Solving trajectory for beta = {}.".format(beta))
sol = solve_ivp(model, t_span = [start, end], y0 = z0, events = hit_ground, max_step = dt)

“”“

据我所知，术语V**2*np.cos(Gamma)/r会引起一些麻烦，因为当我将其乘以0时，结果看起来更像正确的结果。你知道吗

我的代码还有其他部分（比如初始条件和后期处理），但我不久前在这里发布了一篇文章，结果发现这不是提问的正确方式。如果你还需要什么，请告诉我。你知道吗

谢谢你。你知道吗