轨迹图上的奇怪异常

Question

我正在编写一个程序来绘制火箭的轨迹。当我将参数设置为特定值时，我遇到了一个奇怪的结果 - 火箭在燃料结束前不久加速得更快。你知道是什么原因造成的吗？

那是vy vs time的图，黑线是fuel用完的一个时间点。 我尝试深入研究 solve_ivp 以查看在那里作用的力有何变化，但时间步长太大，无法从中理解任何意义。 这是代码：

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

Cd = 0.47 #coefficient of drag
r = 1 #radius of the ball [meters]
A = np.pi * r**2 #cross-section area of the ball [m^2]
g = 9.81 # [m/s^2]
rho = 1.2 # air density [kg/m^3]
rocket_mass = 37 # [kg]
m_fuel = 10 # fuel mass [kg]
m_total = rocket_mass-m_fuel # rocket + fuel mass [kg]
v_e =1000 # velocity of gas [m/s]
v0 =50 # muzle velocity # [m/s]
thetadeg = 30 # angle between the muzzle velocity and the ground [degrees]
capacity =0.01 #amount of fuel (in kilograms) that will exit the rocket during one second when the gas speed is 1 m/s [kg/m]

theta = np.deg2rad(thetadeg)

x0 = 0
y0 = 0
vx0 = v0 * np.cos(theta)
vy0 = v0 * np.sin(theta)

t0 = 0
tf = 10
dt = 0.00001

if v_e > 0 and m_fuel>0:
    fuel_duration = m_fuel/(v_e*capacity)
    dm_over_dt = m_fuel/fuel_duration
else:
    m_fuel = 0
    dm_over_dt = 0
    fuel_duration = 0
    v_e = 0
thr = []
times = []
check = 0
check_hit = 0
x_of_fuel_out = None
x_hit = None
t_impact = None

def mass_of_fuel_from_time(t):
    if t<fuel_duration:
        return m_fuel-dm_over_dt*t
    else:
        return 0
    
def projectile_motion(t, y):
    global check, x_of_fuel_out, x_hit, check_hit, m_fuel,fuel_duration,dm_over_dt, rocket_mass, t_impact

    if 'm_total' not in projectile_motion.__dict__:
        projectile_motion.m_total = m_total

    x, y, vx, vy = y
    v = np.sqrt(vx**2 + vy**2)
    rho_air = rho
    F_gravity = np.array([0, -projectile_motion.m_total * g])
    F_drag = -0.5 * rho_air * v**2 * Cd * A * np.array([vx, vy])/v 

    if 'm_fuel' not in projectile_motion.__dict__:
        projectile_motion.m_fuel = m_fuel
        projectile_motion.t_fuel = fuel_duration
        projectile_motion.fuel = 1

    if t<fuel_duration:
        projectile_motion.m_fuel  = mass_of_fuel_from_time(t)
        projectile_motion.m_total = projectile_motion.m_fuel+rocket_mass
        F_thrust = dm_over_dt*v_e*np.array([np.cos(theta),np.sin(theta)])
    else:
        if check ==0:
            check = 1
            x_of_fuel_out = x
        if t>fuel_duration:
            projectile_motion.fuel_cutoff = t
            projectile_motion.fuel = 0
            projectile_motion.m_fuel  = 0
            projectile_motion.m_total = rocket_mass
        F_thrust = np.array([0, 0])
    if y<0 and check_hit==0:
        check_hit = 1
        x_hit = x
        t_impact = t
    F_total = F_gravity + F_drag + F_thrust
    thr.append(F_thrust)
    times.append(t)
    a = F_total / projectile_motion.m_total
    projectile_motion.t_prev = t
    print(f"t: {t}      x: {round(x,2)},     v: {round(vx,2)}, {round(vy,2)} ,     thrust: {round(F_thrust[0],2)}, {round(F_thrust[1],2)},     drag: {round(F_drag[0],2)}, {round(F_drag[1],2)}     fuel left: {projectile_motion.m_fuel}")
    return [vx, vy, a[0], a[1]]

projectile_motion.t_prev = t0

sol = solve_ivp(projectile_motion, (t0, tf), [x0, y0, vx0, vy0], t_eval = np.linspace(t0, tf, int((tf-t0)/dt)+1))
print(f"IMPACT AT {t_impact} seconds")
vx_values = sol.y[2]
vy_values = sol.y[3]

# Plotting vx against time
fig, ax = plt.subplots()
ax.plot(sol.t, vx_values)
ax.set_xlabel('Time (s)')
ax.set_ylabel('vx (m/s)')
ax.set_xlim([0, t_impact])
ax.set_title('vx for: v0 = {} m/s, v_e = {} m/s, m_rocket = {} kg, m_fuel = {} kg, theta = {} deg, fuel end at {} s'.format(v0,v_e,rocket_mass,m_fuel,thetadeg,fuel_duration))
if x_of_fuel_out != None:
    print(f"Fuel ran out at x = {x_of_fuel_out} m")
    plt.axvline(x=fuel_duration,color="#000000", linestyle='--')
plt.show()

# Plotting vy against time
fig, ax = plt.subplots()
ax.plot(sol.t, vy_values)
ax.set_xlabel('Time (s)')
ax.set_ylabel('vy (m/s)')
ax.set_xlim([0, t_impact])
ax.set_title('vy for: v0 = {} m/s, v_e = {} m/s, m_rocket = {} kg, m_fuel = {} kg, theta = {} deg, fuel and at {} s'.format(v0,v_e,rocket_mass,m_fuel,thetadeg,fuel_duration))
if x_of_fuel_out != None:
    plt.axvline(x=fuel_duration,color="#000000", linestyle='--')
plt.show()

y_max = max(sol.y[1])

fig, ax = plt.subplots()
ax.plot(sol.y[0], sol.y[1], c='red')
ax.plot(sol.y[0][0], sol.y[1][0], color="#000000", label="start")
ax.set_aspect('equal', 'box')
ax.set_xlim([0, x_hit*1.05])
ax.set_ylim([0,y_max*1.05])
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')
ax.set_title('v0 = {} m/s, v_e = {} m/s, m_rocket = {} kg, m_fuel = {} kg, theta = {} deg'.format(v0,v_e,rocket_mass,m_fuel,thetadeg))
if x_of_fuel_out != None:
    plt.axvline(x=x_of_fuel_out, color="#000000",linestyle='--')
else:
    print("There was enough fuel.")
if x_hit != None:
    print(f"Range: {x_hit} m")
print(f"max height: {y_max} m")
print(f"fuel lasted {fuel_duration} s")
plt.show()
print(sol)