使用 scipy.ivp() 实现自由落体球弹跳

Question

我想使用 scipy_ivp() 求解常微分方程，球在 x 方向上以 1 的初始速度落下，在 y 方向上以 0 的初始速度落下。重力加速度为 g = 9.82，当球撞击地面时，其速度应该改变符号并乘以 0.9。但是，使用 events 参数，我发现它无法正常工作。这是我的代码和结果：

from scipy.integrate import solve_ivp
def f_derivs_ivp(t, vars, g = 9.82):
    
    dxdt = vars[2]
    dydt = vars[3]
    
    dvxdt = 0
    dvydt = -g
    
    return dxdt, dydt, dvxdt, dvydt

def bounce(t, y, g = 9.82):
    if y[1] <= 0:
        y[3] = -0.9*y[3]
    #print(y[1])
    return y[1]

#bounce.terminal = True
#bounce.direction = -1

vars0 = np.array([0, 10, 1, 0])

sol = solve_ivp(f_derivs_ivp, [0, 7], vars0, max_step=0.01, events = bounce)


plt.plot(sol.y[0], sol.y[1], "ko")

print(sol.y_events)
print(sol.t_events)

使用 scipy.ivp() 之外的另一种方法，结果应如下所示：

我做错了什么，事件参数如何工作？另请注意，如果我在弹跳函数中写入

return y[1] - 10

，则不会发生任何变化。例如，如果我在弹跳函数的 if 语句中写入

y[3] = 10

，它会上下弹跳，但不是应有的样子。

Answer 1

scipy ODE 求解器不具备用户可定义操作的事件操作机制。事件函数只是轨迹上的标量值函数，其零点（如果设置了交叉方向）是事件。这意味着事件函数用于搜索事件（首先检测符号变化，然后使用布伦特方法进行细化）。每个时间步将被调用多次，无论时间步是否包含事件。

内置操作是“注册”（默认）和“终止”。

您必须立即中断积分，并以反射状态作为初始条件重新启动。分别绘制各个部分或使用串联以获得一个大的解决方案数组。

Answer 2

您确实可以使用solve_ivp的事件处理程序来检测球击中地面时的影响，方法是在y<=0.

时触发模拟终止

终止后，保护最后一个时间步的所有状态，降低速度并反转 y 方向（在您的情况下为 -0.9）。然后将数组作为下一次运行的初始条件传递。循环直到到达定义的模拟结束时间。要完全退出模拟循环，请添加名为“end_of_simulation”的第二个事件检测。如果倾向于 - t < 0, the simulation gets terminated completely.

注意，时间向量以及所有状态向量需要在每个循环中堆叠在一起。

此外，一旦球速度低于某个阈值（例如 0.1m/s），请关闭速度和重力。否则，球就会“掉”到地面上。如果有人有更好的主意来处理这种情况，请告诉我。

总而言之，它看起来像这样：

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


def falling_obj_rhs(t, Y, *args):
    
    g, _= args
    
    dxdt = Y[2,0]
    dydt = Y[3,0]
    
    dvxdt = 0
    dvydt = g
    
    return dxdt, dydt, dvxdt, dvydt

def hit_ground(t, Y, *args): 
    return Y[1]


def end_of_simulation(t, Y, *args):
    _, tend = args
    
    return tend - t
    

def main():

    g       = -9.81
    t_start = 0
    t_end   = 7
    
    # the function y is approaching the x axis from above --> negative slope
    hit_ground.direction       = -1  
    hit_ground.terminal        = True
    end_of_simulation.terminal = True
    
    #       x, y, vx, vy
    ic   = [0, 10, 1, 0]
    tsim = np.array([])
    ypos = np.array([])
    
    while t_start < t_end:
        
        sol = solve_ivp( falling_obj_rhs, 
                         [t_start, t_end], 
                         ic, 
                         method     = 'LSODA',
                         max_step   = 0.01, 
                         events     = [hit_ground, end_of_simulation], 
                         vectorized = True,
                         rtol       = 1e-6, 
                         atol       = 1e-9,
                         first_step = 1e-6,
                         args       = (g, t_end) )
        
        
        # store solution from last time step as initial condition
        ic = sol.y[:,-1]
        # reduce initial condition of y-velocity (energy loss du to collistion)
        # and change sign
        ic[-1] = ic[-1]*-0.9
        # turn off gravity and velocity once the ball is on the ground
        if sol.y[:,-1][-1] < 0.05:
            ic[-1] = 0
            g      = 0
        
        t_start = sol.t[-1]
        tsim = np.concatenate([tsim, sol.t[:]], axis=0)
        ypos = np.concatenate([ypos, sol.y[1]], axis=0)
    
    plt.plot(tsim, ypos, label='obj. trajectory')
    plt.xlabel('t')
    plt.ylabel('y')
    plt.grid()
    plt.legend()


if __name__ == '__main__':
    main()