[CR3BP] 주기궤도의 안정성

어떤 \(\bar{\mathbf{x}}(t)\) 가 다음 미분방정식의 해로 주어지는 주기(period)가 \(T\) 인 주기궤도라고 하자.

 

\[ \dot{\bar{\mathbf{x}}}(t)= \mathbf{f}( \bar{\mathbf{x}} (t)) \tag{1} \]

 

\(\bar{\mathbf{x}}(t)\) 에 약간의 섭동 \(\delta \mathbf{x}(t)\) 을 주고 식 (1)에 대입한 후 테일러 시리즈 1차 근사식을 구하면 다음과 같이 된다.

 

\[ \begin{align} & \dot{\bar{\mathbf{x}}} (t)+ \delta \dot{\mathbf{x}}(t) \approx \mathbf{f}( \bar{\mathbf{x}}(t))+ \left. \frac{ \partial \mathbf{f}}{\partial \mathbf{x} } \right\vert_{\bar{\mathbf{x}}(t) } \delta \mathbf{x}(t) \\ \\ \to \ & \delta \dot{\mathbf{x}}(t)=A(t) \delta \mathbf{x}(t) \tag{2} \end{align} \]

 

\(\bar{\mathbf{x}}(t)\) 가 주기함수이므로 미분도 주기함수다. 따라서 \(A(t)\) 도 주기함수다.

 

\[ A(t)=A(t+T) \tag{3} \]

 

 

한편 식 (2)의 기본행렬 \(\Psi(t)\) 는 다음과 같다 (https://pasus.tistory.com/274).

 

\[ \frac{d}{dt} \Psi (t)=A(t) \Psi (t), \ \ \ \ \ \Psi (t_0 )=I \tag{4} \]

 

여기서 \(A(t)\) 가 주기함수이므로 다음과 같은 상수 비특이행렬(non-singular matrix) \(M\) 이 존재한다.

 

\[ \Psi (t+T)= \Psi (t) M \tag{5} \]

 

여기서 \(M\) 을 모노드로미 행렬(monodromy matrix)이라고 하며 다음과 같이 주어진다.

 

\[ M= \Psi^{-1} (t) \Psi (t+T)= \Phi (t+T,t) \tag{6} \]

 

여기서 \(\Phi (t, \tau )\) 는 상태천이행렬(state transition matrix)로서 다음 미분방정식을 만족한다.

 

\[ \frac{d}{dt} \Phi (t, \tau)=A(t) \Phi (t, \tau ), \ \ \ \ \ \Phi(\tau, \tau)=I \tag{7} \]

 

또한 두 싯점간의 상태변수를 매핑해 준다.

 

\[ \delta \mathbf{x}(t)= \Phi (t, \tau ) \delta \mathbf{x}( \tau ) \tag{8} \]

 

 

 

식 (6)과 (8)에 의하면,

 

\[ \begin{align} \delta \mathbf{x}(t_0+kT) &= \Phi (t_0+kT, t_0+(k-1)T) \delta \mathbf{x}(t_0+(k-1)T) \tag{9} \\ \\ &=M \delta \mathbf{x}(t_0+(k-1)T) \\ \\ &= M \Phi (t_0+(k-1)T, t_0+(k-2)T) \delta \mathbf{x}(t_0+(k-2)T) \\ \\ &=M^2 \delta \mathbf{x}(t_0+(k-2)T) \\ \\ & \ \ \ \ \ \cdots \\ \\ &=M^k \delta \mathbf{x} (t_0 ) \end{align} \]

 

이므로 주기궤도 \(\bar{\mathbf{x}}(t)\) 의 안정성(stability)은 모노드로미 행렬 \(M\) 의 고유값(eigenvalue)에 달려있다. 모노드로미 행렬의 고유값 \(\mu_i\) 을 Floquet multipliers 라고 한다. 고유값이 \(|\mu_j | \lt 1\) 이면 궤도의 섭동은 \(0\) 으로 수렴할 것이고 \( |\mu_j | \gt 1\) 이면 발산할 것이다. 참고로 Floquet exponents \(\gamma_j\) 는 다음과 같이 Floquet multipliers의 로그값으로 정의한다.

 

\[ \gamma_j= \frac{1}{T} \log \mu_j \tag{10} \]

 

이제 라그랑지 포인트(Langrage point)에서의 주기함수의 모노드로미 행렬의 고유값을 구해보도록 하자. 라그랑지 포인트에서의 주기궤도의 선형화 미분방정식은 식 (2)로 주어지며 그 때의 상태변수는 다음과 같다 (https://pasus.tistory.com/280).

 

\[ \delta \mathbf{x}(t)=[ \delta x \ \ \delta y \ \ \delta z \ \ \delta v_x \ \ \delta v_y \ \ \delta v_z ]^T \tag{11} \]

 

한편 원궤도제한 삼체문제(CR3BP)의 운동 방정식은 해밀톤(Hamilton) 방정식으로도 유도할 수 있다 (https://pasus.tistory.com/158). 일반화 좌표 \(\mathbf{q}\) 와 일반화 운동량 \(\mathbf{p}\) 를 다음과 같이 정의하고,

 

\[ \mathbf{q}= \begin{bmatrix} q_1 \\ q_2 \\ q_3 \end{bmatrix} = \begin{bmatrix} x \\ y \\ z \end{bmatrix} , \ \ \ \ \mathbf{p}= \begin{bmatrix} p_1 \\ p_2 \\ p_3 \end{bmatrix} = \begin{bmatrix} \dot{x}-y \\ \dot{y}+x \\ \dot{z} \end{bmatrix} \tag{12} \]

 

해밀토니안 함수(Hamiltonian function)를 다음과 같이 정의한다면,

 

\[ H= \frac{1}{2} \left( (p_1+y)^2+(p_2-x)^2+p_3^2 \right) + U_{eff} \tag{13} \]

 

해밀톤 표준 방정식에 의해서 CR3BP의 운동 방정식을 얻을 수 있다 (https://pasus.tistory.com/147).

 

\[ \begin{align} \dot{\mathbf{q}} &= \frac{ \partial H}{ \partial \mathbf{p} } = \begin{bmatrix} p_1+y \\ p_2-x \\ p_3 \end{bmatrix} \tag{14} \\ \\ \dot{\mathbf{p}} &= - \frac{\partial H}{ \partial \mathbf{q} } = \begin{bmatrix} p_2-x- \bar{U}_x \\ -p_1-y-\bar{U}_y \\ - \bar{U}_z \end{bmatrix} \end{align} \]

 

식 (14)를 상태변수 \(\mathbf{z}= \begin{bmatrix} \mathbf{q} \\ \mathbf{p} \end{bmatrix} \) 의 방정식으로 표현하면 다음과 같다.

 

\[ \begin{align} \dot{\mathbf{z}} &= \begin{bmatrix} \dot{\mathbf{q}} \\ \dot{\mathbf{p}} \end{bmatrix} = \begin{bmatrix} 0 & I \\ -I & 0 \end{bmatrix} \begin{bmatrix} \frac{\partial H}{ \partial \mathbf{q} } \\ \frac{\partial H}{ \partial \mathbf{p} } \end{bmatrix} \tag{15} \\ \\ &= \begin{bmatrix} 0 & I \\ -I & 0 \end{bmatrix} \frac{\partial H}{ \partial \mathbf{z} } \end{align} \]

 

위 식으로부터 기준 주기궤도 \(\bar{\mathbf{z}}\) 에서의 선형화 미분방정식을 다음과 같이 얻을 수 있다.

 

\[ \delta \dot{\mathbf{z}} = J \frac{\partial }{ \partial \mathbf{z}} \left( \frac{\partial H}{\partial \mathbf{z}} \right) \delta \mathbf{z}= J S(\bar{\mathbf{z}} ) \delta \mathbf{z} \tag{16} \]

 

여기서

 

\[ J= \begin{bmatrix} 0 & I \\ -I & 0 \end{bmatrix} , \ \ \ \ \ S(\bar{\mathbf{z}} )=\frac{\partial }{ \partial \mathbf{z}} \left( \frac{\partial H}{\partial \mathbf{z}} \right) \]

 

이다. 선형 미분방정식 (16)의 상태천이행렬 \(\tilde{\Phi} (t, \tau )\) 는 다음과 같이 구할 수 있다.

 

\[ \frac{d}{dt} \tilde{\Phi} (t, \tau)= JS(\bar{\mathbf{z}} ) \tilde{\Phi}(t, \tau ) \tag{17} \]

 

그리고 시간 \(t\) 의 상태변수 \(\delta \mathbf{z}(t)\) 는 다음과 같이 계산할 수 있다.

 

\[ \delta \mathbf{z}(t) =\tilde{\Phi} (t, t_0 ) \delta \mathbf{z}(t_0) \tag{18} \]

 

 

 

이제 상태천이행렬 \(\tilde{\Phi} (t, t_0)\) 가 심플렉틱 행렬(symplectic matrix)이라는 것을 증명하려고 한다 (https://pasus.tistory.com/275). 먼저 행렬 \(U(t)\) 를 다음과 같이 정의한다.

 

\[ U(t)= \tilde{\Phi}^T (t, t_0)J \tilde{\Phi} (t, t_0), \ \ \ \ \ U(t_0 )=J \tag{19} \]

 

위 식을 미분하면 식 (17)에 의해서 다음과 같이 된다.

 

\[ \begin{align} \dot{U}(t) &= \dot{\tilde{\Phi}}^T (t, t_0 )J \tilde{\Phi} (t, t_0 ) + \tilde{\Phi}^T (t, t_0 )J \dot{\tilde{\Phi}} (t, t_0 ) \tag{20} \\ \\ &= \tilde{\Phi}^T (t, t_0 )S( \bar{\mathbf{z}} ) J^T J \tilde{\Phi}(t, t_0 )+ \tilde{\Phi}^T (t, t_0 )JJS( \bar{\mathbf{z}}) \tilde{\Phi} (t, t_0 ) \\ \\ &= \tilde{\Phi}^T (t, t_0 )S( \bar{\mathbf{z}}) \tilde{\Phi} (t, t_0 )- \tilde{\Phi}^T (t, t_0 )S(\bar{\mathbf{z}}) \tilde{\Phi} (t, t_0 ) \\ \\ &=0 \end{align} \]

 

따라서

 

\[ U(t)= \tilde{\Phi}^T (t, t_0 )J \tilde{\Phi}(t, t_0 )= J \tag{21} \]

 

이기 때문에 \(\tilde{\Phi}(t, t_0 )\) 는 심플렉틱 행렬이다. 심플렉틱 행렬의 특성 (https://pasus.tistory.com/275)에 의해서 \(\tilde{\Phi}(t, t_0 )\) 의 행렬식(determinant)은 항상 \(1\) 이고, \(\lambda\) 가 \(\tilde{\Phi}(t, t_0 )\) 의 고유값이면, \(\lambda^{-1}, \ \bar{\lambda}, \ \bar{\lambda}^{-1}\) 도 같은 다중도(multiplicity)를 갖는 고유값이 된다.

이제 시스템 (2)의 상태천이행렬을 구하기 위하여 식 (18)을 이용하기로 한다. 상태변수 \(\delta \mathbf{z}\) 와 \(\delta \mathbf{x}\) 의 좌표변환 행렬 \(Q\) 는 식 (11)과 (12)에 의해서 다음과 같이 구할 수 있다.

 

\[ \delta \mathbf{z}=Q \delta \mathbf{x}, \ \ \ \ \ Q= \begin{bmatrix} I & 0 \\ K & I \end{bmatrix}, \ \ \ \ \ K=\begin{bmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \tag{22} \]

 

식 (22)를 식 (18)에 대입하면 다음과 같이 된다.

 

\[ \delta \mathbf{x}(t) =Q^{-1} \tilde{\Phi} (t, t_0 ) Q \delta \mathbf{x} (t_0) \tag{23} \]

 

위 식을 식 (8)과 비교하면 시스템 (2)의 상태천이행렬 \(\Phi (t, t_0 )\) 를 다음과 같이 계산할 수 있다.

 

\[ \Phi (t, t_0) =Q^{-1} \tilde{\Phi} (t, t_0 ) Q \tag{24} \]

 

위 식에 의하면 \(\Phi (t, t_0 )\) 와 \(\tilde{\Phi} (t, t_0 )\) 는 서로 상사(similarity)관계에 있으므로 고유값은 서로 같다. 그리고 행렬식도 \(\det \Phi (t, t_0 )= \det \tilde{\Phi}(t, t_0 ) =1\) 이다.

식 (24)를 (21)에 대입하면,

 

\[ \begin{align} & \Phi^T (t, t_0 ) Q^T J Q \Phi (t, t_0 )= Q^T JQ \tag{25} \\ \\ & Q^T JQ= \begin{bmatrix} 2K & I \\ -I & 0 \end{bmatrix} \end{align} \]

 

이 되어서 \(\Phi (t, t_0 )\) 는 \(Q^T JQ\) 에 대한 심플렉틱 행렬이 된다.

시스템 (2)의 모노드로미 행렬은 식 (6)에 의하면 다음과 같다.

 

\[ M= \Phi (t_0+T, t_0) \tag{26} \]

 

따라서 모노드로미 행렬은 행렬식이 항상 \(1\) 이며 \(\mu\) 가 \(M \) 의 고유값이면, \( \mu^{-1}, \ \bar{\mu}, \ \bar{\mu}^{-1}\) 도 같은 다중도를 갖는 고유값이 된다.

한편 식 (2)에서 A(t)가 주기함수이므로 식 (8)에 의하면

 

\[ \begin{align} \delta \mathbf{x}(t_0+T) = \delta \mathbf{x} (t_0 ) &= \Phi (t_0 +T, t_0 ) \delta \mathbf{x}(t_0 ) \tag{27} \\ \\ &=M \delta \mathbf{x} (t_0 ) \end{align} \]

 

이 되므로 \(M\) 은 고유값 \(\mu_1=1\) 을 가지며 그에 해당하는 고유벡터는 \(\delta \mathbf{x}(t_0 )\) 이다. \( M \in \mathbb{R}^{6 \times 6}\) 이므로 \(M\) 의 고유값은 중복된 값을 포함 \(6\) 개가 있으며 행렬식이 \(1\) 이므로 고유값을 모두 곱하면 \(1\) 이다. 따라서 각 고유값은 다음과 같이 주어진다.

 

\[ \begin{align} & \mu_1= \mu_2=1 \tag{28} \\ \\ & \mu_3 \mu_4=1, \ \ \ \ \ \mu_5 \mu_6=1 \end{align} \]

 

주기궤도의 안정성은 \(\mu_3, \ \mu_4, \ \mu_5, \ \mu_6\) 에 의해서 결정된다. 고유값은 서로 역수의 관계에 있으므로 \(4\) 개의 고유값 중 실수(real number)가 있다면 한 고유값은 \(1\) 보다 크거나 \(-1\) 보다 작을 것이기 때문에 궤도는 불안정해진다. 모든 고유값이 단위 원(unit circle)상에 있는 결레복소수라면 궤도는 안정하다.

 

 

모노드로미 행렬 \(M\) 은 주기궤도를 따라 상태천이행렬 \(\Phi (t_0+t, t_0 )\) 를 시간 \(t_0\) 에서 \(t_0+T\) 까지 수치적분하여 얻을 수 있다. 행렬이 얻어지면 고유값을 수치적으로 계산할 수 있다.

한편 식 (1)을 시간 미분하면 다음과 같다.

 

\[ \frac{ d \dot{\bar{\mathbf{x}}} (t) }{dt} = \frac{ \partial \mathbf{f}( \bar{\mathbf{x}}(t)) }{\partial \bar{\mathbf{x}}(t) } \dot{\bar{\mathbf{x}}} (t)= A(t) \dot{\bar{\mathbf{x}}} (t) \tag{29} \]

 

식 (29)와 (2)에 의하면

 

\[ \delta \mathbf{x}(t)= \dot{\bar{\mathbf{x}}}(t)= \mathbf{f}( \bar{\mathbf{x}}(t)) \tag{30} \]

 

이다. 따라서 \(M\) 이 고유값 \(\mu_1=1\) 을 가질 때의 고유벡터는 \(\delta \mathbf{x}(t_0 )\) 는 \(\delta \mathbf{x}(t_0 )= \dot{\bar{\mathbf{x}}}(t_0)\) 로서 주기궤도의 접선 방향이다.