라그랑지 행성 방정식 (Lagrange Planetary Equation)

이체문제는 우주에는 두 개의 질점만 존재하며, 중력이 두 질점 사이에 작용하는 유일한 힘이라는 가정을 기반으로 한다. 이체문제에서 이 힘을 제외한 모든 힘을 섭동력(perturbation force)이라고 한다. 두 질점 운동의 일반적인 섭동력에는 비구형 중심체, 대기 항력, 추진 추력, 태양 복사 압력, 제3의 질점에 의한 중력 등이 있다. 섭동력은 이체문제의 케플러 궤도에 교란을 가하여 정상적인 궤도에서 벗어나는 현상을 초래한다.

파라미터 변분법(VOP, variation of parameters)은 섭동력에 의해 교란된 동적 시스템의 풀이에 적합한 방법이다. 이 방법은 교란되지 않은 시스템 해(solution)의 상수(constant) 파라미터를 시변(time-varying) 파라미터로 일반화할 수 있다면 교란되지 않은 시스템의 해를 사용하여 교란된 시스템의 해를 나타낼 수 있다는 것을 기반으로 한다 (https://pasus.tistory.com/249).

교란되지 않은 시스템은 다음과 같이 이체문제 운동방정식이다.

 

\[ \ddot{\mathbf{r}} + \frac{\mu}{r^3} \mathbf{r}=0 \tag{1} \]

 

위 식의 해는 다음과 같다.

 

\[ \mathbf{r}(t)= \mathbf{f}( \mathbf{c},t) \tag{2} \]

 

여기서 \(\mathbf{c}\) 는 초기값의 조합으로 이루어진 상수벡터로서 여기서는 6개의 궤도요소(orbital elements)를 나타낸다. 위 식은 교란되지 않은 시스템의 해이므로 미분하면 다음과 같다.

 

\[ \dot{\mathbf{r}} = \mathbf{v}= \frac{\partial \mathbf{f}}{\partial t} \tag{3} \]

 

여기서 상미분 대신 편미분을 사용한 이유는 벡터 \(\mathbf{c}\) 가 상수벡터임을 강조하기 위해서다. 위 식을 한번 더 미분하고 식 (1)에 대입하면 다음과 같다.

 

\[ \frac{\partial^2 \mathbf{f} }{ \partial t^2 }+ \frac{ \mu}{ \lVert \mathbf{f} \rVert^3} \mathbf{f}=0 \tag{4} \]

 

임의의 섭동 가속도로 인해 교란된 이체문제 운동방정식은 다음과 같다.

 

\[ \ddot{\mathbf{r}}+ \frac{\mu }{r^3} \mathbf{r} = \mathbf{a}_p \tag{5} \]

 

여기서 \(\mathbf{a}_p\) 는 섭동 가속도이다. VOP에 의하면 식 (2)에서 상수 \(\mathbf{c}\) 를 시변 파라미터로 간주한다면 식 (2)도 식 (5)의 해가 될 수 있다. 즉, 식 (5)의 해는 다음과 같이 쓸 수 있다.

 

\[ \mathbf{r}= \mathbf{f}( \mathbf{c}(t), t) \tag{6} \]

 

위 식을 시간 미분하면 다음과 같다.

 

\[ \dot{\mathbf{r}}= \mathbf{v}= \frac{\partial \mathbf{f}}{\partial t} + \frac{\partial \mathbf{f} }{\partial \mathbf{c}} \dot{\mathbf{c }} \tag{7} \]

 

이제 수식을 보다 간단히 하기 위해서, 또는 \(\mathbf{v}\) 가 교란되지 않은 경우의 \(\mathbf{v}\) 와 동일한 수학적 형태를 갖게 만들기 위해서 다음과 같은 제약조건을 부과하기로 한다.

 

\[ \frac{\partial \mathbf{f} }{\partial \mathbf{c}} \dot{\mathbf{c }} = 0 \tag{8} \]

 

그러면 식 (7)은 식 (3)과 동일해 진다.

 

\[ \dot{\mathbf{r}} = \mathbf{v} = \frac{ \partial \mathbf{f}}{\partial t } \tag{9} \]

 

식 (8)을 라그랑지 제약조건(Lagrange’s constraint) 또는 오스큘레이션 제약조건(osculation constraint)이라고 한다. 이는 물리적으로 교란된 운동과 교란되지 않은 운동의 속도 벡터가 동일해야 한다는 것을 의미한다.

 

 

식 (9)를 미분하면 다음과 같다.

 

\[ \ddot{\mathbf{r}} = \frac{\partial}{\partial t} \left( \frac{ \partial \mathbf{f}}{\partial t} \right)+ \frac{\partial}{\partial \mathbf{c}} \left( \frac{\partial \mathbf{f}}{\partial t} \right) \dot{\mathbf{c}} \tag{10} \]

 

식 (6)과 (10)을 식 (5)에 대입하면,

 

\[ \frac{\partial^2 \mathbf{f} }{\partial t^2 }+ \frac{\partial }{\partial \mathbf{c}} \left( \frac{\partial \mathbf{f} }{\partial t} \right) \dot{\mathbf{c}}+ \frac{\mu }{ \lVert \mathbf{f} \rVert^3} \mathbf{f}= \mathbf{a}_p \tag{11} \]

 

이 되는데 식 (4)에 의해서 위 식은 다음과 같이 된다.

 

\[ \frac{\partial }{\partial \mathbf{c}} \left( \frac{\partial \mathbf{f} }{\partial t} \right) \dot{\mathbf{c}} = \mathbf{a}_p \tag{12} \]

 

식 (8)과 (12)는 시변 파라미터 \(\mathbf{c}(t)\) 가 만족해야할 제약조건이다. 이를 한 개의 식으로 표현하면 다음과 같다.

 

\[ \begin{bmatrix} \frac{\partial \mathbf{f} }{\partial \mathbf{c}} \\ \frac{ \partial }{\partial \mathbf{c}} \left( \frac{\partial \mathbf{f}}{\partial t} \right) \end{bmatrix} \dot{\mathbf{c}} = \mathcal{L} \dot{\mathbf{c}} =\begin{bmatrix} 0 \\ \mathbf{a}_p \end{bmatrix} \tag{13} \]

 

여기서 \( \mathcal{L}=\begin{bmatrix} \frac{\partial \mathbf{f} }{\partial \mathbf{c}} \\ \frac{ \partial }{\partial \mathbf{c}} \left( \frac{\partial \mathbf{f}}{\partial t} \right) \end{bmatrix} \) 을 라그랑지 행렬이라고 한다. 위 식에 의하면 \(\dot{\mathbf{c}} \) 은 다음과 같이 계산할 수 있다.

 

\[ \dot{\mathbf{c}}= \mathcal{L}^{-1} \begin{bmatrix} 0 \\ \mathbf{a}_p \end{bmatrix} \tag{14} \]

 

 

 

라그랑지는 섭동력이 보존력(conservative force)이라는 가정 하에 라그랑지 행렬의 역행렬을 계산하는 편리한 방법을 개발했다. 보존력은 위치만의 함수인 어떤 스칼라 함수의 그래디언트로 표현될 수 있는 힘을 의미한다 (https://pasus.tistory.com/153). 그러면 식 (5)의 섭동 가속도는 다음과 같이 쓸 수 있다.

 

\[ \mathbf{a}_p= \frac{dR(\mathbf{r})}{d \mathbf{r}} \tag{15} \]

 

여기서 \(R(\mathbf{r})\) 은 단위 질량당 섭동력의 포텐셜 에너지이다. 그러면 식 (8)과 (12)는 다음과 같이 쓸 수 있다.

 

\[ \begin{align} & \frac{\partial \mathbf{r}}{\partial \mathbf{c} } \dot{\mathbf{c}}=0 \tag{16} \\ \\ & \frac{\partial \mathbf{v} }{\partial \mathbf{c}} \dot{\mathbf{c}} = \frac{dR}{d \mathbf{r}} \tag{17} \end{align} \]

 

식 (16)의 양변에 \( \left( \frac{\partial \mathbf{v}}{\partial \mathbf{c}} \right)^T \) 을 곱하고, 식 (17)의 양변에 \( \left( \frac{\partial \mathbf{r}}{\partial \mathbf{c} } \right)^T\) 을 곱한 후 빼면 다음과 같은 식이 된다.

 

\[ \begin{align} \left( \frac{\partial \mathbf{r}}{\partial \mathbf{c}} \right)^T\frac{ \partial \mathbf{v}}{\partial \mathbf{c}} \dot{\mathbf{c}}- \left( \frac{\partial \mathbf{v}}{\partial \mathbf{c}} \right)^T \frac{\partial \mathbf{r}}{\partial \mathbf{c}} \dot{\mathbf{c}} = \left( \frac{\partial \mathbf{r}}{\partial \mathbf{c}} \right)^T \frac{dR}{d\mathbf{c}} \tag{18} \end{align} \]

 

또는

 

\[ \begin{align} L \dot{\mathbf{c}} = \left( \frac{\partial \mathbf{r}}{\partial \mathbf{c}} \right)^T \frac{dR}{d\mathbf{c}} \tag{19} \end{align} \]

 

이다. 여기서 \(L\) 을 새로운 라그랑지 행렬 또는 라그랑지 계수 행렬이라고 한다.

 

\[ \begin{align} L = \left( \frac{\partial \mathbf{r}}{\partial \mathbf{c}} \right)^T\frac{ \partial \mathbf{v}}{\partial \mathbf{c}} - \left( \frac{\partial \mathbf{v}}{\partial \mathbf{c}} \right)^T \frac{\partial \mathbf{r}}{\partial \mathbf{c}} \ \tag{20} \end{align} \]

 

그리고 \(L\) 의 각 구성요소 \(L_{ij}\) 를 라그랑지 브래킷(Lagrangian brackets)이라고 한다.

 

\[ \begin{align} L_{ij} = \left( \frac{\partial \mathbf{r}}{\partial c_i} \right)^T\frac{ \partial \mathbf{v}}{\partial c_j} - \left( \frac{\partial \mathbf{v}}{\partial c_i} \right)^T \frac{\partial \mathbf{r}}{\partial c_j} = [ c_i, \ c_j ] \tag{21} \end{align} \]

 

라그랑지 브래킷은 다음 세 가지 특성을 갖는다.

 

\[ \begin{align} & L_{ij}=-L_{ji} \tag{22} \\ \\ & L_{ii}= 0 \\ \\ & \frac{\partial L_{ij}}{\partial t}=0 \end{align} \]

 

위에서 첫번째와 두번째 특성은 라그랑지 행렬이 빗대칭(skew-symmetric) 행렬, 즉 \(L=-L^T\) 임을 나타내는데 이는 식 (20)의 정의로 부터 쉽게 알 수 있다. 세번째 특성은 라그랑지 행렬이 명시적인 시간 함수가 아니라는 것을 의미하는 데, 이는 다음과 같이 증명할 수 있다.

 

\[ \begin{align} \frac{\partial L_{ij}}{\partial t} &= \left(\frac{ \partial }{\partial t} \frac{\partial \mathbf{r}}{\partial c_i } \right)^T \frac{ \partial \mathbf{v}}{\partial c_j }-\left(\frac{ \partial }{\partial t} \frac{\partial \mathbf{v}}{\partial c_i } \right)^T \frac{ \partial \mathbf{r}}{\partial c_j } \tag{23} \\ \\ & \ \ \ \ \ + \left( \frac{\partial \mathbf{r} }{\partial c_i } \right)^T \frac{\partial }{\partial t} \frac{\partial \mathbf{v} }{\partial c_j }-\left( \frac{\partial \mathbf{v} }{\partial c_i } \right)^T \frac{\partial }{\partial t} \frac{\partial \mathbf{r} }{\partial c_j } \\ \\ &= \left( \frac{\partial \mathbf{v}}{\partial c_i } \right)^T \frac{ \partial \mathbf{v}}{\partial c_j }-\left( \frac{ \partial }{\partial t} \frac{\partial \mathbf{v}}{\partial c_i } \right)^T \frac{ \partial \mathbf{r}}{\partial c_j } \\ \\ & \ \ \ \ \ + \left( \frac{\partial \mathbf{r} }{\partial c_i } \right)^T \frac{\partial }{\partial t} \frac{\partial \mathbf{v} }{\partial c_j }-\left( \frac{\partial \mathbf{v} }{\partial c_i } \right)^T \frac{\partial \mathbf{v} }{\partial c_j } \\ \\ &=-\left( \frac{ \partial }{\partial t} \frac{\partial \mathbf{v}}{\partial c_i } \right)^T \frac{ \partial \mathbf{r}}{\partial c_j } + \left( \frac{\partial \mathbf{r} }{\partial c_i } \right)^T \frac{\partial }{\partial t} \frac{\partial \mathbf{v} }{\partial c_j } \end{align} \]

 

여기서 잠시 이체문제의 운동 방정식 유도 과정을 보면,

 

\[ \begin{align} \frac{\partial \mathbf{v}}{\partial t} =- \frac{dV}{d\mathbf{r}}, \ \ \ V(\mathbf{r})= - \frac{\mu}{r} \tag{24} \end{align} \]

 

이므로 식 (23)에 대입하면,

 

\[ \begin{align} \frac{\partial L_{ij}}{\partial t} &= \left( \frac{\partial }{\partial c_i} \frac{\partial V}{\partial \mathbf{r} } \right)^T \frac{ \partial \mathbf{r}}{\partial c_j } - \left( \frac{\partial}{\partial c_j } \frac{\partial V }{\partial \mathbf{r} } \right)^T \frac{\partial \mathbf{r} }{\partial c_i } \tag{25} \\ \\ &= \frac{\partial^2 V}{\partial c_i \partial c_j }-\frac{\partial^2 V}{\partial c_j \partial c_i } \\ \\ &=0 \end{align} \]

 

이 된다. 따라서 \(L\) 은 시간에 명시적으로 의존하지 않으므로 \(L\) 을 계산할 때 궤도 상의 위치는 중요하지 않다. 계산을 쉽게 하는 위치라면 궤도 상의 어느 위치든 관계 없다.

 

 

이제 벡터 \(\mathbf{c}\) 를 다음과 같이 6개의 궤도요소(orbital elements)라고 하자 (https://pasus.tistory.com/286).

 

\[ \begin{align} \mathbf{c} = \begin{bmatrix} a & e & i & \Omega & \omega & M_0 \end{bmatrix}^T \tag{26} \end{align} \]

 

여기서 \(M_0\) 는 평균비행각(mean anomaly)의 초기값이다. 섭동력이 없다면 식 (26)의 벡터 \(\mathbf{c}\) 는 상수벡터다.

Battin의 책 'An Introduction to the Mathematics and Methods of Astrodynamics'에 의하면 식 (21)의 라그랑지 브래킷이 다음과 같이 계산되어 있다.

 

\[ \begin{align} & [\omega, \ \Omega]= \left( \frac{\partial \mathbf{r} }{\partial \omega} \right)^T \frac{\partial \mathbf{v}}{\partial \Omega}-\left( \frac{\partial \mathbf{v} }{\partial \omega} \right)^T \frac{\partial \mathbf{r}}{\partial \Omega}=0, \ \ \ [\omega, \ i]=0 \tag{27} \\ \\ & [a,\ i]=0, \ \ \ [e,\ i]=0, \ \ \ [e, \ a]=0 \\ \\ &[ M_0, \ \Omega ]=0, \ \ \ [M_0,\ i]=0, \ \ \ [M_0, \ \omega ]=0, \ \ \ [M_0, \ e]=0 \\ \\ & [i, \ \Omega ]=nab \sin i, \ \ \ [a, \ \omega]=-\frac{nb}{2}, \ \ \ [e, \ \omega]= \frac{na^3 e}{b} \\ \\ &[M_0, \ a]= \frac{na}{2}, \ \ \ [a, \ \Omega]=- \frac{nb}{2} \cos i, \ \ \ [e, \ \Omega ]= \frac{na^3 e}{b} \cos i \end{align} \]

 

여기서 \(n\) 은 평균 각속도, \(b=a \sqrt{1-e^2}\) 이다. 식 (27)을 이용하여 식 (19)의 \(L\) 을 만들고 역행렬을 구한 다음, \(\mathbf{c}\) 의 시간변화율을 계산하면 다음과 같다.

 

\[ \begin{align} & \frac{da}{dt} = \frac{2}{na} \frac{\partial R}{\partial M_0 } \tag{28} \\ \\ & \frac{de}{dt} =-\frac{b}{na^3 e} \frac{\partial R}{\partial \omega}+ \frac{b^2}{na^4 e} \frac{ \partial R}{\partial M_0 } \\ \\ & \frac{di}{dt}=-\frac{1}{nab \sin i} \frac{\partial R}{\partial \Omega} + \frac{\cos i}{nab \sin i} \frac{\partial R}{\partial \omega} \\ \\ & \frac{d\Omega}{dt}= \frac{1}{nab \sin i} \frac{\partial R}{\partial i} \\ \\ & \frac{d\omega}{dt}=- \frac{\cos i}{nab \sin i} \frac{\partial R}{\partial i} + \frac{b}{na^3 e} \frac{\partial R}{\partial e} \\ \\ & \frac{dM_0}{dt}=-\frac{2}{na} \frac{\partial R}{\partial a}- \frac{b^2}{na^4 e} \frac{\partial R}{\partial e} \end{align} \]

 

식 (28)을 라그랑지 행성 방정식(Lagrange planetary equation)이라고 한다. 위 식은 경사각이 \(i=0^0\) 일 때, 이심율이 \(e=0\) 일 때 특이값을 갖는다.