如何在 solve_ivp (Python) 中均匀分配时间？

Question

我想针对给定的时间步求解一个函数，但即使在将 t_eval 设置为均匀分布的时刻数组之后，该函数也会对彼此大小不同的时间步执行计算（我知道这可能有用，但是为了我的目的，它必须是偶数）。 此代码：

t0 = 0
tf = 10
dt = 0.001

#def projectile_motion(t, y):
    #here is the function evaluated by solve_ivp
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))

导致函数在 1E-8 秒这样的时间步被评估，尽管 t_eval 明确表示每 10^-3 秒执行一次。我该怎么做才能让 solve_ivp 均匀地选择时间步长，一次一个 dt？

编辑： 整个代码：

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]
m_total = 27 # rocket + fuel mass [kg]
m_fuel = 4 # fuel mass [kg]
v_e =800 # velocity of gas [m/s]
v0 =1000 # muzle velocity # [m/s]
thetadeg =25 # 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 = 100
dt = 0.001

if v_e > 0 and m_fuel>0:
    fuel_duration = m_fuel/(v_e*capacity)
    m_fuel_duration = m_fuel/fuel_duration
else:
    m_fuel = 0
    m_fuel_duration = 0
    fuel_duration = 0

check = 0
check_hit = 0
x_of_fuel_out = None
x_hit = None

def projectile_motion(t, y):
    global check, x_of_fuel_out, x_hit, check_hit
    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 projectile_motion.m_fuel>0:
        projectile_motion.t_fuel -= dt
        projectile_motion.m_fuel -= dt*m_fuel_duration
        projectile_motion.m_total -= dt*m_fuel_duration
        momentum_fuel = v_e * dt*m_fuel_duration
        F_thrust = momentum_fuel* np.array([np.cos(theta),np.sin(theta)])/dt
    else:
        if check ==0:
            check = 1
            x_of_fuel_out = x
        if t>fuel_duration:
            projectile_motion.fuel_cutoff = t
            projectile_motion.fuel = 0
        F_thrust = np.array([0, 0])
    if y<0 and check_hit==0:
        check_hit = 1
        x_hit = x
    F_total = F_gravity + F_drag + F_thrust
    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}")
    t+=dt
    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))
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], c="#000000", label="start")
ax.set_aspect('equal', 'box')
ax.set_ylim([0, y_max*1.05])

ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')
ax.set_title('Projectile Trajectory')

if x_of_fuel_out != None:
    print(f"Fuel ran out at x = {x_of_fuel_out}m")
    plt.axvline(x=x_of_fuel_out, linestyle='--')
else:
    print("There was enough fuel.")

if x_hit != None:
    print(f"Range: {x_hit} m")
    ax.set_xlim([0,x_hit])
print(f"max height: {y_max} m")

plt.show()