#python #events #scipy #ode #montecarlo
Вопрос:
Мне нужно смоделировать пенальти, когда мяч летит во все стороны. Для этого я использую solve_ivp и прекращаю интеграцию, когда мяч пересекает заднюю линию. Когда это произойдет, я хочу использовать значения, найденные при остановке интеграции, чтобы увидеть, находится ли мяч в этой точке в пределах размеров цели (и считать его целью). Тем не менее, решение.y не дает требуемой точности, которую я хочу. Желаемые значения находятся в пределах решения.y_events, однако, я, похоже, не могу их получить. Я также не могу найти информацию об этом в Интернете. Мой код:
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 def is_it_goal(vector): if vector.status == 1: if (vector.y[1][len(vector.y[1])-1] gt; 0) and (-3.36 lt; vector.y[0][len(vector.y[1])-1] lt; 3.36) and (vector.y[2][len(vector.y[1])-1] lt; 2.44): print("GOAAAAAAAAAAAAL!") elif (vector.y[1][len(vector.y[1])-1] lt; 0) and (-3.36 lt; vector.y[0][len(vector.y[1])-1] lt; 3.36) and (vector.y[2][len(vector.y[1])-1] lt; 2.44): print("Own goal?! Why?") else: print("Awwwwh") 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)) is_it_goal(solution) Здесь я хочу изменить vector.y[1][len(vector.y[1])-1] на что-то вроде vector.y_events…
Комментарии:
1. Вы можете сократить
vector.y[1][len(vector.y[1])-1]доvector.y[1][-1]. Вам нужно оценитьsolution.status, со смыслом,solution.messageчтобы выяснить, было ли событие..t_events[k]является списком времени срабатывания для события№k,.y_events[k]содержит список соответствующих точек (несовместимых с форматированиемsolution.y).2. .y_events состоит только из двух входных данных. К сожалению, это моя проблема. т. е.: y.события = [массив([], dtype=float64), массив([[ 88.53280031, -100. , 59.34496498, 48.46126737, -54.73820686, 32.48436969]])], y_events[0] = [], и, y_events[1] = [[ 88.53280031 -100. 59.34496498 48.46126737 -54.73820686 32.48436969]]. Из y_events[1] я хотел бы получить первое и третье значения, чтобы проверить, находятся ли они внутри цели.
3. Если
len(sol.t_events[k]gt;0), то вы можете оценитьsol.y_events[k][0,[0,2]], чтобы получить первый и третий компоненты.4. Не знаю, что ты сделал, но это работает! Большое спасибо!!!