Я работаю над проектированием траектории космического корабля и внедряю алгоритм для нахождения периодов c орбит, которые необходимо обнаруживать, когда траектория пересекает заданную c плоскость. Я использую scipy.integrate solve_ivp и его events функциональность.
В идеале я хотел бы подсчитать, когда инициируемые события останавливают интеграцию только после того, как сработало определенное количество событий, но поскольку я не смог найти подобную опцию в документации solve_ivp , Я просто использую цикл for для объединения этих дуг интеграции и подсчета количества событий. Однако моя проблема в том, что solve_ivp застревает, когда начальное условие вызывает событие, и event.terminal = True.
Ниже приведен пример кода:
import numpy as np
from scipy.integrate import solve_ivp
def eq_motion(t,Y):
Ydot = np.zeros((6,1))
r = (Y[0,0]**2 + Y[1,0]**2 + Y[2,0]**2)**0.5
Ydot[0] = Y[3,0]
Ydot[1] = Y[4,0]
Ydot[2] = Y[5,0]
Ydot[3] = 2*Y[4,0] + 3*Y[0,0] - Y[0,0]/r**3
Ydot[4] = -2*Y[3,0] - Y[1,0]/r**3
Ydot[5] = -Y[2,0] - Y[2,0]/r**3
return Ydot
def event_xz_plane(t,Y):
return Y[1]
event_xz_plane.terminal = True
event_xz_plane.direction = 0
# initial conditions
Y0 = np.array([-0.006489114696252, 0, 0.107462057974196 + 0.0006, 0, 4.165257687202607,0])
t_span = np.array([0, 1.892845635529040*2])
n_plane_intersections = 4
for i in np.arange(n_plane_intersections):
integrate = solve_ivp(fun=eq_motion, t_span=t_span, y0=Y0, method = 'LSODA',
vectorized = True, rtol = 1e-13, atol = 1e-14,
events = (event_xz_plane),dense_output=True)
print('Number of integrations before event: {}'.format(integrate.y.shape[1]-1))
print(integrate.y_events[0])
print('\n')
# attempting to chain the different integration arcs
Y0 = integrate.y[:,-1]
Выбранные начальные условия близки к периодам c орбиты и дают следующую траекторию где вы можете видеть, что есть 4 пересечения с плоскостью xz, включая начальное условие, которое находится на этой плоскости.
Тем не менее, вышеприведенный скрипт выводит следующее:
Number of integrations before event: 1
[[-0.00648911 0. 0.10806206 0. 4.16525769 0. ]]
Number of integrations before event: 1
[[-0.00648911 0. 0.10806206 0. 4.16525769 0. ]]
Number of integrations before event: 1
[[-0.00648911 0. 0.10806206 0. 4.16525769 0. ]]
Number of integrations before event: 1
[[-0.00648911 0. 0.10806206 0. 4.16525769 0. ]]
где вы можете видеть, что это в основном застревает в начальном состоянии. Это происходит, конечно, потому что начальное условие находится в терминальном условии, однако, кажется, нет никакой возможности «пройти» первую итерацию интегрирования.
Я попытался написать функцию события как
def event_xz_plane(t,Y):
return Y[1] + 1e3*(t==0)
, чтобы она игнорировала начальное состояние, но затем просто «зависала» на следующем пересечении плоскости, и в любом случае это не не похоже на хорошее решение.
Единственное другое решение, о котором я подумал, это установить terminate=False для функции события и проверить количество запущенных событий и вектор состояния в этих случаях. Тем не менее, это не идеально для моей реализации, так как это зависит от установленного времени t_span, которое неизвестно в алгоритме. То есть, если мне нужно получить и n количество пересечений плоскостей, мне придется итеративно корректировать t_span до тех пор, пока я не достигну этого числа, вместо того, чтобы просто позволить интегрированию работать до тех пор, пока не сработают n пересечения плоскостей.
Есть ли способ преодолеть это?