如何访问scipy.solve_ivp返回的y_events以获取点球模拟终止时刻精确值

解决方案 #

核心说明 #

你调用solve_ivp时传入了两个终端事件(goal, own_goal)，返回的solution.y_events是长度为2的列表：

• 索引0对应goal事件的触发结果，存储该事件触发时刻的状态值，形状为(触发次数, 状态变量总数)
• 索引1对应own_goal事件的触发结果，结构同上
  由于两个事件都设置了terminal=True，每次积分终止时只会有一个事件被触发，对应位置的数组长度为1，另一个为空，我们直接取触发事件对应的状态值即可，该值是scipy做根求解得到的事件触发时刻的精确值，精度远高于积分步长的最后一个点。

修改点 #

只需修改is_it_goal函数的取值逻辑即可，替换原有取积分最后一步状态的写法：

def is_it_goal(solution):
    if solution.status == 1:
        # 读取触发事件对应的精确状态值
        if len(solution.y_events[0]) > 0:
            # 触发对方球门底线事件 y=11
            hit_x, hit_y, hit_z = solution.y_events[0][0][:3]
            if -3.36 < hit_x < 3.36 and 0 < hit_z < 2.44:
                print("GOAAAAAAAAAAAAL!")
            else:
                print("Awwwwh")
        elif len(solution.y_events[1]) > 0:
            # 触发己方球门底线事件 y=-100
            hit_x, hit_y, hit_z = solution.y_events[1][0][:3]
            if -3.36 < hit_x < 3.36 and 0 < hit_z < 2.44:
                print("Own goal?! Why?")
            else:
                print("Awwwwh")
        # 若加入了落地事件可额外补充落地判断逻辑
    else:
        print("Not even close, lol")

修改后完整可运行代码 #

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

#### Constants
# Number of simulations
number_of_penalty_shots = 10

# Angle of the shots
theta = np.random.uniform(0, 2.0*np.pi, number_of_penalty_shots)
phi = np.random.uniform(0, np.pi, number_of_penalty_shots)

# Velocity of the ball
v_magnitude = 80

### Starting Position Ball (defined as the penalty stip)
pos_x = 0.0
pos_y = 0.0
pos_z = 0.0

in_position = np.array([pos_x, pos_y, pos_z]) # Inital position in m

def homo_magnetic_field(t, vector):
    vx = vector[3] # dx/dt = vx
    vy = vector[4] # dy/dt = vy
    vz = vector[5] # dz/dt = vz   
        
    # 开启阻力和重力的话取消下方注释即可
    # ax = -0.05*vector[3] # dvx/dt = ax
    # ay = -0.05*vector[4] # dvy/dy = ay
    # az = -0.05*vector[5] - constants.g #dvz/dt = az
    
    ax = 0
    ay = 0
    az = 0
        
    dvectordt = (vx,vy,vz,ax,ay,az)
    return(dvectordt)

def goal(t, vector):
    return vector[1] - 11

def own_goal(t,vector):
    return vector[1] + 100

def ground(t,vector):
    return vector[2]

goal.terminal=True
own_goal.terminal=True
ground.terminal=True # 如果需要判断是否落地可以把ground也加入events列表

def is_it_goal(solution):
    if solution.status == 1:
        # 读取触发事件对应的精确状态值
        if len(solution.y_events[0]) > 0:
            # 触发对方球门底线事件 y=11
            hit_x, hit_y, hit_z = solution.y_events[0][0][:3]
            if -3.36 < hit_x < 3.36 and 0 < hit_z < 2.44:
                print("GOAAAAAAAAAAAAL!")
            else:
                print("Awwwwh")
        elif len(solution.y_events[1]) > 0:
            # 触发己方球门底线事件 y=-100
            hit_x, hit_y, hit_z = solution.y_events[1][0][:3]
            if -3.36 < hit_x < 3.36 and 0 < hit_z < 2.44:
                print("Own goal?! Why?")
            else:
                print("Awwwwh")
        elif len(solution.y_events[2]) > 0:
            # 触发落地事件
            print("Awwwwh, ball hit the ground")
    else:
        print("Not even close, lol")

# Integrating
time_range = (0.0, 10**2)

for i in range(number_of_penalty_shots):
    v_x = v_magnitude*np.sin(phi[i])*np.cos(theta[i])
    v_y = v_magnitude*np.sin(phi[i])*np.sin(theta[i])
    v_z = v_magnitude*np.cos(phi[i])
    in_velocity = np.array([v_x, v_y, v_z])
    initial_point = np.array([in_position, in_velocity])
    start_point = initial_point.reshape(6,)
    
    solution = integrate.solve_ivp(homo_magnetic_field , time_range, start_point,events=(goal, own_goal, ground))
    is_it_goal(solution)

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