初步认识一些动力系统，如太阳系系统，混沌系统等，并研究其数值解法
本文公式较多，在浏览器中将会花较长时间用于渲染公式。
最近看boost::numeric::odeint的文档，发现了一些很有趣的动力系统。
太阳系系统
这个看起来简单，无非就是万有引力的动力学，但是深入探究并不简单。这里我们想考虑太阳、木星、土星、天王星、海王星等天体的运动以及他们的相互作用。
NOTE
一个简单的物理史：二体问题是可积的，可以理论分析得到解析解，而三体问题则不存在解析解。在19世纪末，庞加莱发现了这个问题，并成为混沌理论的始祖。
理想的天体运动都是保守场，因此可以列出哈密顿方程。首先是动量的表达式（当然这里的动量的量纲是速度的量纲，但是物理含义是动量）
\begin{equation}
\frac{dq_i}{dt}=p_i
\end{equation}
\begin{equation}
\frac{dp_i}{dt}=\frac{1}{m_i} \sum_{i\neq j}F_{ij}
\end{equation}
那么哈密顿函数的表达式为
\begin{equation}
H=\sum_i \frac{p_i^2}{2m_i}+\sum_j V(q_i,q_j)
\end{equation}
$V(q_i,q_j)$是$i,j$相互的势能。哈密顿方程的动力学方程为
在boost库中求解这个动力系统，有专门的辛算法：Runge-Kutta-Nystroem算法，在boost库中是一个叫symplectic_rkn_sb3a_mclachlan的步进器。事实上，对于很多实际生活中的动力系统，只要它有惯性，不耗散，就能够写出哈密顿函数，从而写出正则方程组，用这个算法求解。三体问题就是混沌的了，因此在这个世纪太阳系问题的求解中，必须要选择符合实际的初始条件，否则得不到想要的结果。
而对于一个实际的航天器，考虑的动力学方程通常如下：
这里某几个项的含义是（不一一赘述了）：
- $\nabla \phi_{sj}^o$ 中心天体非圆带来的扰动项
- 其它天体的引力扰动项
- $\frac{\dot m_s}{m}\frac{d\mathbf{r}}{dt}$发动机推力
- 大气阻力
- $\frac{P_{SR}C_RA_{sun}}{m_s}\hat{\mathbf{r}}_{sun}$是太阳辐射光压
- $\frac{\mu}{c^2r^3}((4\frac{\mu}{r}-v^2)\mathbf{r}+4(\mathbf{r\cdot v})\mathbf{v})$为史瓦西解（Schwarzsehild solution)
- $2(\mathbf{\Omega}\times\mathbf{v})$是测地线精度(广义相对论造成的空间扭曲形成的误差)
- $2\frac{\mu}{c^2r^3}(\frac{3}{r^2}(\mathbf{r\times v})(\mathbf{r \cdot J})+(\mathbf{v\times J}))$是兰斯-蒂林效应(Lens-Thirring Precession)
广义相对论对受力的修正作用就是式(\ref{eqs_force})的最后三项，即[1] 如果在J2000坐标系下，那么
- 是一个近似
- $c$是光速
- $\mathbf{r}$是在J2000坐标系下的位置矢量
- $\mathbf{v}$是在J2000坐标系下的速度矢量
- $\mathbf{r}_{B/S}$是中心天体相对于太阳的位置矢量
- $\mathbf{v}_{B/S}$是中心天体相对于太阳的速度矢量
- $\mathbf{J}$是中心天体的单位质量角动量，有$\mathbf{J}=\mathbf{R}^{I/F}_B[0\quad 0\quad \frac{2}{5}R_B^2\omega_B]^T$，而$\mathbf{R}^{I/F}_B$是中心天体的体坐标系到惯性系的旋转矩阵，$R_B$是中心天体的平均赤道半径，$\omega_B$是中心天体的自旋速度。
式(\ref{eqs_force})数值求解通常的方法是广义的龙格库塔法，见我之前的讨论。
混沌系统与李雅普诺夫指数
这个系统我是没怎么看懂。。。设$x$是自变量，$\delta x$是扰动量。扰动量满足线性微分方程，但是是时间相关的：
\begin{equation}
\frac{d\delta x}{dt}=J(x)\delta x
\end{equation}
李雅普诺夫指数被定义为扰动随指数形式增长的那个对数。如果存在一个李雅普诺夫指数大于0那么扰动会发散，形成混沌。如果所有李雅普诺夫指数都小于0，那么扰动将收敛至一个点。
为了计算李雅普诺夫指数，需要首先求解扰动动力学方程，每过k步后将结果正交化，李雅普诺夫指数可以用取对数的方法得到，再多次计算取平均值。boost::numeric::odeint的test有相关代码。这个东西我是真的没搞明白是什么……
刚性方程
在数学领域中，刚性方程是指一个微分方程，其数值分析的解只有在时间间隔很小时才会稳定，只要时间间隔略大，其解就会不稳定。特点是特征值的实部均为负数并且实部最大的的特征值的实部与实部最小的的特征值的实部的比值远大于1。化学反应的动力学方程可能会是这样的方程，如果各个子反应的速率相差非常大。 刚性系统需要用Rosenbrock method来求解。boost::numeric::odeint提供了这个算法，然而坑的是必须和boost::ublas一起用。。。odeint的文档提供了一个对比，如果用Rosenbrock算法，只需要71步迭代，而龙格库塔算法，准确说是Dormand-Prince 5算法需要1531次迭代。
Stuart-Landau振子
微分方程也可以是复数的形式，如下面这个例子，Stuart-Landau oscillator: \begin{equation} \frac{d\Psi}{dt}=(1=i\eta)\Psi+(1+i\alpha)|\Psi|^2\Psi \end{equation} 这个系统表示了什么，我也不知道。从数值求解来说，这个系统虽然数域是复数域，不过依然可以用4阶龙格库塔算法求解。
网格动力系统
一个突出的例子是Fermi-Pasta-Ulam系统，这是一个非线性的哈密顿系统，哈密顿函数为
\begin{equation}
H=\sum_i \frac{p_i^2}{2}+\frac{1}{2}(q_{i+1}-q_i)^2+\frac{\beta}{4}(q_{i+1}-q_i)^4
\end{equation}
和太阳系动力系统一样，这个系统也适合用辛求解算法，但是更为简单，因为$\frac{dq_i}{dt}=p_i$.
有纪念意义的是，FPU系统的求解是世界上最早的计算机上的数值实验。在1953年，洛斯阿拉莫斯，在世界上最早的计算机之一（MANIAC I）进行了这项研究，并触发了一个数学物理学的新领域。
振子集合
N-N耦合相位振子（N all-to all coupled phase oscillator）定义为 \begin{equation} \frac{d\phi_k}{dt}=\omega_k+\frac{\epsilon}{N}\sin(\phi_i-\phi_k) \end{equation} 各个振子的自然频率满足某种分布，$\epsilon$是耦合强度。如果我们取$\omega_i$的分布为洛伦兹分布（其实就是柯西分布，物理学上一般叫做洛伦兹分布），那么我们会观察到一个有趣的现象：当耦合强度超过某个临界值后，振子间会发生相位转移，伺候振子会以共同的频率运动。这种转移被称为Kuramoto transition。这个动力学系统本身的数值求解直接用4阶龙格库塔就行。
参考文献
[1]Huang C, Ries J C, Tapley B D, et al. Relativistic effects for near-earth satellite orbit determination[J]. Celestial Mechanics & Dynamical Astronomy, 1990, 48(2):167-185.