The Quaternion kinematical differential equation (QKDE) plays a key role in navigation, control and guidance systems. Although explicit symplectic geometric algorithms (ESGA) for this problem are available, there is a lack of a unified way for constructing high order symplectic difference schemes with configurable order parameter. We present even order explicit symplectic geometric algorithms to solve the QKDE with diagonal Padè approximation and Cayley transform. The maximum absolute error for solving the QKDE is $\mathcal{O}(\tau^{2\ell})$ where  $\tau$ is the time step and $\ell$ is the order parameter. The linear time complexity and constant space complexity of computation  as well as the simple algorithmic structure show that our algorithms are appropriate for realtime applications in aeronautics, astronautics, robotics, visual-inertial odemetry and so on. The performance of the proposed algorithms are verified and validated by mathematical analysis and numerical simulation.

Keywords: Time-varying system; Guidance navigation and control; Quaternion kinematical differential equation (QKDE); Numerical algorithms;
Symplectic geometric algorithm; Padè approximation; Cayley transform

§ INTRODUCTION

Quaternions have been widely utilized
in aerospace[1, 2, 3, 4, 5] aeronautics[6, 7, 8, 9, 10, 11, 12, 13], robotics  [14, 15, 16, 17, 18, 19, 20], computer vision[21] and computer graphics[22]. A comprehensive review of fundamental results about linear quaternion differential equations is presented by Kou and Xia [23].  The key advantages for quaternions stimulating many researches and applications lie in two facts:

* Geometrical singularity can be avoided when the romation matrices are parameterized with quaternions instead of Euler angles.

* Multiplications and additions are enough for  computing quaterions, which implies that quaternions are appropriate for realtime applications because of low computational complexity.

For the navigation, control and guidance  problems, the key issue is to find stable and precise numeric solution to the Robinson's quaternion kinematical differential equation (QKDE) [24]

\begin{equation} \label{eq-QKDE}
\left\{
\begin{split}
&\fracode{\quat{q}}{t} = \frac{1}{2} \mat{\Omega}(\vec{\omega}(t)) \quat{q}, \quad  t\in [t_0, t_f] \\
&\quat{q}(t_0) = \quat{q}_0
\end{split}
\right.
\end{equation}

where $\quat{q}  = \trsp{[e_0, e_1, e_2, e_3]}$
is the vector representation of the quaternion $\quat{q}$,$[t_0, t_f]$ is the time span, $\quat{q}_0$ is the initial state, $\vec{\omega} =\trsp{[\omega_1(t), \omega_2(t), \omega_3(t)]}\in \ES{R}{3}{1}$
is the angular velocity[In references of aerospace and aeronautics, the notation $\vec{\omega} = \trsp{[p,q,r]}$ is used. In this paper, we use $p_i$ and $q_i$ to denote the $i$-th generalized momentum and displacement.], and $\mat{\Omega} = \mat{\Omega}(\vec{\omega}(t))$ is specified by $\vec{\omega}$ such that

\begin{equation} \label{eq-A}
\mat{\Omega}
=-\trsp{\mat{\Omega}} = \begin{bmatrix}
0 & -\omega_1 & -\omega_2  & -\omega_3 \\
\omega_1 & 0 & \omega_3 & -\omega_2 \\
\omega_2 & -\omega_3 & 0 & \omega_1 \\
\omega_3 & \omega_2 & -\omega_1 & 0
\end{bmatrix}.
\end{equation}

If $\vec{\omega}(t)$ is time independent, then (<ref>)  can be called autonomous QKDE (A-QKDE) or linear time-invariant QKDE (LTI-QKDE) since it is an autonomous system in the sense of classic mechanics and control theory and also it is a linear time-invariant (LTI) system in the sense of signals analysis.  On the contrary, if $\vec{\omega}(t)$ is time dependent,
(<ref>) can be named with non-autonomous QKDE (NA-QKDE) or linear time-varying QKDE (LTV-QKDE) in general.
Although the QKDE (<ref>) is a linear ordinary differential equation (ODE), it is challenging to find its solution for two reasons:

* The time-varying vector $\vec{\omega}(t)$ leads to a time-dependent matrix $\mat{\Omega}(\vec{\omega}(t))$. However, there is no general method for finding the analytical solution (AS) to a LTV system.

* The QKDE is critically (neutrally) stable and the numerical solution is sensitive to the computational errors [1] because the eigenvalues of $\mat{\Omega}$ include two pure imaginary  numbers
$\pm \mj \norm{\vec{\omega}} = \pm \mj \sqrt{\omega_1^2 + \omega_2^2 + \omega_3^2}$ where $ \mj = \sqrt{-1}$.

it is necessary to find a stable and precise integration method in the sense of long term time for solving  (<ref>).

The second approach is built on  symplectic difference scheme (SDS) which avoids accumulative computational errors automatically. It should be pointed out that a TFDS may  be  symplectic. Wang [25] compared the Runge-Kutta scheme with SDS. The Gauss-Legendre (G-L) difference scheme [26, 27], an accurate and implicit Runge-Kutta method, is an implicit SDS. However, its computational complexity is $\mathcal{O}(n^2)$ [28] which means that it is  not appropriate for real-time applications.

Zhang et al [29] proposed   two explicit symplectic geometric algorithms (ESGA) — the ESGA-I and ESGA-II — for the LTI-QKDE (A-QKDE) and LTV-QKDE (NA-QKDE) respectively,   which are not only precise in the sense of long term time, but also with linear time computational complexity and constant space computational complexity. However, there are some disadvantages for these two algorithms:

* they  are not constructed with a unified way and can not be generalized for higher order cases;

* the order of the SDS's accuracy is not configurable, which limits the applicaitons for moving objects with high velocity such as airplane, missile, torpedo, space shuttle and so on.

* the discrete time sequence corresponding to the interval $[t_0,t_f]$ is $\set{t_0, t_0+\frac{\tau}{2}, t_0 +\frac{3\tau}{2}, \cdots, t_0+\frac{(2k+1)\tau}{2}, \cdots}$ instead of $\set{t_0, t_0+\tau, \cdots, t_0+k\tau, \cdots}$ which implies that we have to treat the time intervals  carefully;

The main purpose of this paper is to propose even order explict symplectic geometric algorithms (EoESGA) for solving the QKDE with Padé approximation based on the theorems proved. The flow chart of solving the LTV-QKDE with the EoESGA via Padé approximation and Cayley transform is given in Figure <ref>. There are six key steps for solving the LTV-QKDE: initialization, sampling the angular velocity, calculating the key auxiliary parameters, constructing the system matrix, calculating the symplectic transition matrix and calculating the quaternions iteratively.
The EoESGA is parameterized by a positive integer $\ell$, which is the most imporpant feature and differes from all of the available algorithms for solving QKDE.

Flow chart of the EoESGA for solvingthe the LTV-QKDE with accuracy parameter $\ell$.

The contents of this paper are organized as follows:
The preliminaries are given in Section <ref>. Section <ref> deals with the SDS for the QKDE. In Section <ref> we design the SGA and analyze the computational complexity. Section <ref> deals with performance evaluations. Finally, Section <ref> presents the conclusions.

§ PRELIMINARIES

§.§ Symplectic Transition Matrix

Let $\tau$ be the time step and $\vec{z}[k] = \vec{z}(t_0+k\tau)$ be the sample value at the discrete time $t_k = t_0+k \tau$.
With the
initial condition

\begin{equation}
\vec{z}[0]= \trsp{[p_1(t_0), \cdots, p_N(t_0), q_1(t_0), \cdots, q_N(t_0)]},
\end{equation}

the symplectic difference scheme (SDS) for the Hamiltonian canonical equation, see Appendix <ref> and (<ref>), is given by

\begin{equation}
\vec{z}[k+1] = \mathcal{G}_\tau (\vec{z}[k]), \quad k = 0, 1, 2, \cdots
\end{equation}

where $ \mathcal{G}_\tau:   \ES{R}{2N}{1}\to \ES{R}{2N}{1}$
is the transition operator, also named as the transition mapping, and

\begin{equation}
\mat{G}^{k+1}_{k}(\tau) = \fracpde{\vec{z}[k+1]}{\vec{z}[k]} = \fracpde{(\mathcal{G}_\tau \vec{z}[k])}{\vec{z}[k]}\in \ES{R}{2N}{2N},
\end{equation}

is the transition matrix (Jacobi matrix) of $\mathcal{G}_\tau$
at $\vec{z}[k]$
such that $\mat{G}^{k+1}_{k}(\tau)$  is symplectic, viz.,

\begin{equation}\label{SpMatDef}
\trsp{\mat{G}^{k+1}_{k}(\tau)} \cdot \mat{J} \cdot \mat{G}^{k+1}_{k}(\tau) = \mat{J}.
\end{equation}

The key problem of designing SGA is to determine the $\mathcal{G}_\tau$ or the corresponding  Jacobi matrix $\mat{G}^{k+1}_{k}(\tau)$. In this paper, we  have $\mathcal{G}_\tau  = \mat{G}^{k+1}_{k}(\tau)$ since (<ref>) is a linear ODE.
For the SDS of LTI-QKDE, the notation $\mat{G}^{k+1}_{k}(\tau)$ will be replaced by  $\mat{G}(\tau)$ since it does not depend on the time index $k$. In this paper, we will also take the notations $\mat{G}^{k+1}_{k}(\ell,\tau)$ and $\mat{G}(\ell, \tau)$ where $\ell$ is the order parameter for accuracy.

§.§ Padé Approximation and Symplectic Difference Scheme

For non-negative integers $\ell$ and $m$, the $(\ell+m)$-th order Padé approximation to $\exp(x)$ is defined by

\begin{equation}
\exp(x) =  \frac{N_{\ell m}(x)}{D_{\ell m}(x)}  + \mathcal{O}(\abs{x}^{\ell +m+1})
\end{equation}

\begin{align}
N_{\ell m}(x) &= \sum^m_{k=0} \frac{(\ell +m -k)! m!}{(\ell +m)! k!(m-k)!}x^k,\\
D_{\ell m}(x) &= \sum^m_{k=0} \frac{(\ell +m -k)! \ell!}{(\ell +m)! k!(\ell-k)!}(-x)^k.
\end{align}

\begin{equation}
\eta^{\ell}_r = \frac{\ell-r}{(2\ell-r)(r+1)}, \quad r\in \set{0,1,\cdots,\ell-1},
\end{equation}

then for $\ell = m$ we have the diagonal Padè approximation

\begin{equation}
g_{\ell}(x) =  \frac{N_{\ell \ell}(x)}{D_{\ell \ell}(x)}= \frac{P_\ell(x)}{P_\ell(-x)}
\end{equation}

\begin{equation}
\begin{split}
&= \sum_{k=0}\frac{(2\ell-k)!}{(2\ell)!}\cdot\frac{\ell!}{(\ell-k)!}\cdot\frac{x^k}{k!} \\
%&= 1 + \sum^{\ell}_{k=1}\frac{(2\ell-k)!}{(2\ell)!}\cdot\frac{\ell!}{(\ell-k)!}\cdot\frac{x^k}{k!} \\
%&= 1 + \sum^{\ell}_{k=1}
%\frac{\ell(\ell-1)\cdots(\ell-k+1)}{2\ell(2\ell-1)\cdots (2\ell-k+1)}\cdot \frac{x^{k}}{k!}\\
%&= 1+ \sum^{\ell}_{k=1}x^k\prod^{k-1}_{r=0} \frac{\ell-r}{(2\ell-r)(r+1)}  \\
&= 1+ \sum^{\ell-1}_{k=0}x^{k+1}\prod^{k}_{r=0} \frac{\ell-r}{(2\ell-r)(r+1)}  \\
&= 1+ \sum^{\ell-1}_{k=0}x^{k+1}\prod^k_{r=0}\eta^{\ell}_r.
%&= 1 + x\eta^{\ell}_0 + x^2 \eta^{\ell}_0\eta^{\ell}_1 + x^3 \eta^{\ell}_0\eta^{\ell}_1\eta^{\ell}_2 \\
%&\qquad + \cdots + x\eta^{\ell}_0\eta^{\ell}_1\cdots\eta^{\ell}{\ell-1}.
\end{split}
\end{equation}

For illustration, we have

\begin{align}
P_1(x)  &=   1+\frac{x}{2}, \\
P_2(x)  &=   1 + \frac{x}{2} + \frac{x^2}{12}, \\
P_3(x)  &=   1+ \frac{x}{2} + \frac{x^2}{10} + \frac{x^3}{120},\\
P_4(x)  &=   1+ \frac{x}{2} + \frac{3x^2}{28} + \frac{x^3}{84} + \frac{x^4}{1680}.
\end{align}

It should be pointed out that the constant $\frac{\ell!}{(2\ell)!}$ is kept in the expression of $P_\ell(x)$ to avoid the overflow problem in numerical computations since $(2\ell)!$ may be a large number for $\ell > 5$.

The diagonal Padé approximation $g_\ell(x)$ for $\exp(x)$ can be used to construct SDS. Actually, we have the following theorem [30].

For the linear Hamiltonian system

\begin{equation} \label{eq-LinearHS}
\fracode{\vec{z}}{t} = \mat{F} \vec{z}, \quad \mat{F} = \inv{\mat{J}} \mat{C} , \quad \trsp{\mat{C}} = \mat{C},
\end{equation}

where $\mat{F} = \inv{\mat{J}} \mat{C}$ is infinitesimal symplectic, i.e.,

\begin{equation}
\mat{J}\mat{F} + \trsp{\mat{F}}\mat{J} = \mat{O}.
\end{equation}

\begin{equation}
\begin{split}
&= P_\ell(\tau \mat{F}) \inv{[P_\ell(-\tau \mat{F})]},
\end{split}
\end{equation}

then the difference schemes

\begin{equation}
\vec{z}[k+1] = g_{\ell}(\tau \mat{F})\vec{z}[k], \quad \ell = 1, 2, \cdots
\end{equation}

are symplectic of $2\ell$-th order accuracy for step size $\tau$.

§.§ Cayley Transform and Euler's Formula

In [29], the authors proved two lemmas about Cayley transform and Euler's formula.

Suppose matrix $\mat{M}\in \ES{R}{n}{n}$ is skew-symmetric, i.e., $\trsp{\mat{M}} = -\mat{M}$, and there exists
a positive constant $\gamma$ such that $\mat{M}^2 = -\gamma^2 \mat{I}$. Let $\hat{\mat{M}} = {\gamma}^{-1}\mat{M}$, then $\hat{\mat{M}}^2 = -\mat{I}$.
For any $y\in \mathbb{R}$, we have

\begin{equation}
\me^{y\mat{M}} = \cos(y\gamma)\mat{I}
+ \sin(y\gamma)\hat{\mat{M}}.
\end{equation}

Suppose matrix $\mat{M}\in \ES{R}{n}{n}$ is skew-symmetric, i.e., $\trsp{\mat{M}} = -\mat{M}$, and there exists
a positive constant $\gamma$ such that $\mat{M}^2 = -\gamma^2 \mat{I}$. Let $\hat{\mat{M}} = {\gamma}^{-1}\mat{M}$, then $\hat{\mat{M}}^2 = -\mat{I}$.

\begin{equation}
\phi(\lambda) = \frac{1+\lambda}{1-\lambda}=(1+\lambda)\inv{(1-\lambda)}
\end{equation}
be the
Cayley transformation, then for any $y\in \mathbb{R}$ the Cayley-Euler formula

\begin{equation}
\begin{split}
\phi(y\mat{M})
&= \frac{\mat{I}+y\mat{M} }{\mat{I}-y\mat{M}}
\frac{1}{1+\alpha}[(1-\alpha)\mat{I} +2y\mat{M}]\\
&=\cos \delta(y, \gamma) \mat{I} + \sin \delta(y, \gamma) \hat{\mat{M}}\\
&= \me^{\delta(y, \gamma)\mat{M}}
\end{split}
\end{equation}

holds, in which $\delta = \delta(y, \gamma) = 2\cdot \arctan (y\gamma)$ and $\alpha = y^2\gamma^2$. Furthermore, $\phi(y\mat{M})$ is an orthogonal matrix.

We remark that the primitive definiton of Cayley transform is $\varphi(\lambda) = \frac{1-\lambda}{1+\lambda} = \phi(-\lambda)$ for $\lambda \in \mathbb{C}$ [31]. In this paper, we just need the real Cayley transform and the mappings $\varphi$ and $\phi$ are equivalent.

§ SYMPLECTIC DIFFERENCE SCHEMES FOR QKDE

§.§ Relation of Padè Approximation and Cayley Transform

Generally, the Padè approximation involves the ratio of  polynomials with the same order. It is interesting that the expression of Padè approximation can be simplified into the Cayley transform for solving the QKDE.
We now prove some interesting lemmas and theorems for constructing the SDS for LTI-QKDE.

If there exists some positive constant $c$ such that $x^2 = -c$, then each Padè approximation $g_\ell(x)$ can be represented by the  Cayley transform

\begin{equation}
g_\ell(x) = \frac{1 + \beta(\ell,c)x}{1-\beta(\ell,c)x} = \phi[\beta(\ell,c) x],
\end{equation}

where $\beta(\ell,c)$
is a constant determined by $\ell$ and $c$.

Proof: Let
$\mfloor{x} = m$ in which $\exists m \in \mathbb{Z}$ such that $m\le x < m +1$. Let

\begin{equation}
s_1 = \mfloor{\frac{\ell-1}{2}}, \quad s_2 = \mfloor{\frac{\ell}{2}},
\end{equation}

we immediately have

\begin{equation}
= s_1 + 1 -\mod(2) =
\left\{
\begin{split}
s_1, & \qquad \hbox{for odd $\ell$};\\
s_1 + 1, & \qquad \hbox{for  even $\ell$}.
\end{split}
\right.
\end{equation}

When $x^2 = -c$  we have

\begin{equation}
\begin{split}
&=  1 + \sum^{\ell-1}_{k=0} x^{k+1} \prod^k_{r=0}\eta^{\ell}_r\\
&= 1 + \sum^{s_1}_{j=0} x^{2j+1} \prod^{2j}_{r=0}\eta^{\ell}_r + \sum^{s_2}_{j=1} x^{2j} \prod^{2j-1}_{r=0}\eta^{\ell}_r\\
% &= 1 + \sum^{s_1}_{j=0} (-c)^jx \prod^{2j}_{r=0}\eta^{\ell}r + \sum^{s_2}_{j=1} (-c)^j \prod^{2j-1}_{r=0}\eta^{\ell}_r\\
&=1  + \sum^{s_2}_{j=1} (-c)^j \prod^{2j-1}_{r=0}\eta^{\ell}_r
+ x\sum^{s_1}_{j=0} (-c)^j \prod^{2j}_{r=0}\eta^{\ell}_r.
\end{split}
\end{equation}

\begin{equation}\label{eg-ab-lj}
\begin{split}
a^{\ell}_j &= \prod^{2j}_{r=0}\eta^{\ell}_r, \quad 0 \le j \le s_1   \\
b^{\ell}_j &= \prod^{2j-1}_{r=0}\eta^{\ell}_r,  \quad 1 \le j \le s_2, \\
b^{\ell}_0 & = 1,
\end{split}
\end{equation}

\begin{equation}
n(s_1,x) = \sum^{s_1}_{j=0} a^{\ell}_j x^j, \quad
d(s_2,x) = \sum^{s_2}_{j=0} b^{\ell}_j x^j,
\end{equation}

\begin{equation}
\begin{split}
%&= \sum^{s_2}_{j=0} b^{\ell}_j \cdot (-c)^j +x \sum^{s_1}_{j=1} a^{\ell}_j \cdot (-c)^j\\
&=  d( s_2,-c) + x \cdot n(s_1,-c).
\end{split}
\end{equation}

\begin{equation}
\beta(\ell,c) = \frac{n(s_1,-c)}{d(s_2,-c)}
=\frac{\sum\limits^{s_1}_{j=0} a^{\ell}_j (-c)^j }{ \sum\limits^{s_2}_{j=0} b^{\ell}_j (-c)^j},
\end{equation}

\begin{equation*}
g_\ell(x) = \frac{P_\ell(x)}{P_\ell(-x)} = \frac{1 + \beta(\ell,c)x}{1-\beta(\ell,c)x} = \phi\left(\beta(\ell,c) x\right).
\end{equation*}

§.§ Procedure for Computing $\beta(\ell,c)$

The flow chart calculating the auxilary parameters $c$ and $\beta$ is illustrated in Figure <ref>.

Flow chart of computing the key auxilary parameters $c$ and $\beta$.

There are two key tasks for computing $\beta(\ell,c)$:

* determine the coefficients $a^{\ell}_j$ and $b^{\ell}_j$;

* calculate the values of polynomials $n(s_1,-c)$ and $d( s_2, -c)$ when $c$ is determined.

For the coefficients of interest, we can use iterative algorithms to generate them. For the computations of polynomials, we can use the famous rule by  Jiushao Qin and Horne independently [32, 33, 34], see Algorithm <ref> in Appendix <ref>.

There are two typical ways to compute the coefficients $a^{\ell}_j$ and $b^{\ell}_j$ for the polynomials $n(s_1,-c)$ and $d(s_2,-c)$.

1) Parallel Iterative method:
According to the definition of $a^{\ell}_j$, we can find the following iterative formulas for the coefficients $a^{\ell}_j$:

\begin{equation} \label{eq-i-a}
\left\{
\begin{split}
&= a^{\ell}_{j}\cdot\eta^{\ell}_{2j+1}\cdot
\eta^{\ell}_{2j+2}, \quad 0 \le j \le s_1 -1; \\
a^{\ell}_0  &= \frac{1}{2}.
\end{split}
\right.
\end{equation}

Similarly, we can obtain

\begin{equation} \label{eq-i-b}
\left\{
\begin{split}
&= b^{\ell}_{j} \cdot \eta^{\ell}_{2j}\cdot \eta^{\ell}_{2j+1}, \quad 0\le  j \le s_2 -1 \\
b^{\ell}_0  &= 1.
\end{split}
\right.
\end{equation}

2) Alternative Iterative Method:
We can also obtain the following alternative iterative formula for the coefficients $a^{\ell}_j$ and $b^{\ell}_j$

\begin{equation} \label{eq-ai-ab}
\left\{
\begin{split}
b^{\ell}_0 &= 1, a^{\ell}_0 = \frac{1}{2},\\
b^{\ell}_{j+1} &=a^{\ell}_j \cdot \eta^{\ell}_{2j+1}, \quad 0 \le j \le s_1 -1\\
a^{\ell}_{j+1} &=b^{\ell}_{j+1}\cdot \eta^{\ell}_{2j+2},\quad 0 \le j \le s_1 -1\\
b^{\ell}_{s_2} &= a^{\ell}_{s_1}\cdot \eta^{\ell}_{2s_1+1}, \quad \hbox{for even $\ell$ }.
\end{split}
\right.
\end{equation}

It is easy to find that the computations in alternative iterative formula (<ref>) are more efficient than those in (<ref>) and (<ref>) since only $\eta^{\ell}_{2j+1}$ and $\eta^{\ell}_{2j+2}$ will be computed.

The iterative algorithms for generating $a^{\ell}_j$ and $b^{\ell}_j$ are shown in Figure <ref>. We can find some interesting facts from this figure:

* The structures of the graph for  odd $\ell$ and even $\ell$ are different. If $\ell$ is even, then $s_2 = s_1 + 1 = \ell/2$ and we need an extra iterative step from $a^{\ell}_{s_1}$ to $b^{\ell}_{s_2}$ if compared with the case of odd $\ell$.

* For the parallel iterative method, the paths for computing $\set{a^{\ell}_j}^{s_1}_{j=0}$ and $\set{b^{\ell}_j}^{s_2}_{j=0}$ are  $a^{\ell}_0 \rightarrow a^{\ell}_1 \rightarrow  \cdots \rightarrow a^{\ell}_j \rightarrow \cdots \rightarrow a^{\ell}_{s_1}$ and  $b^{\ell}_0 \rightarrow b^{\ell}_1 \rightarrow  \cdots \rightarrow b^{\ell}_j \rightarrow \cdots \rightarrow b^{\ell}_{s_2}$ respectively. It is obvious that these two paths are independent and the computations involves lots of the products $\eta^{\ell}_{2j+1}\eta^{\ell}_{2j+2}$ and $\eta^{\ell}_{2j}\eta^{\ell}_{2j+1}$.

* For the alternative iterative method, the path for computing the coefficients of interest is $b^{\ell}_0 \rightarrow a^{\ell}_0 \rightarrow b^{\ell}_1 \rightarrow a^{\ell}_1 \rightarrow \cdots \rightarrow b^{\ell}_j \rightarrow a^{\ell}_j \rightarrow \cdots \rightarrow b^{\ell}_{s_1} \rightarrow a^{\ell}_{s_1} \rightarrow b^{\ell}_{s_2}$ when $\ell$ is even or $b^{\ell}_0 \rightarrow a^{\ell}_0 \rightarrow b^{\ell}_1 \rightarrow a^{\ell}_1 \rightarrow \cdots \rightarrow b^{\ell}_j \rightarrow a^{\ell}_j \rightarrow \cdots \rightarrow b^{\ell}_{s_1} \rightarrow a^{\ell}_{s_1}$ when $\ell$ is odd.

[$\ell$ is odd, $s_2 = s_1=(\ell-1)/2$.]

\xymatrix{
a^{\ell}_0\ar[r]^{\eta^{\ell}_1\eta^{\ell}_2}\ar[dr]_{\eta^{\ell}_1} &
a^{\ell}_1\ar[r]^{\eta^{\ell}_3\eta^{\ell}_4}\ar[dr]_{\eta^{\ell}_3} &
a^{\ell}_2\ar[r]^{\eta^{\ell}_5\eta^{\ell}_6}\ar[dr]_{\eta^{\ell}_5} &
<EMAIL_ADDRESS><EMAIL_ADDRESS>&
a^{\ell}_j\ar[rr]^{\eta^{\ell}_{2j+1}\eta^{\ell}_{2(j+1)}}\ar[drr]^{\eta^{\ell}_{2j+1}} &
b^{\ell}_0\ar[r]_{\eta^{\ell}_0\eta^{\ell}_1}\ar[u]^{\eta^{\ell}_0} &
b^{\ell}_1\ar[r]_{\eta^{\ell}_2\eta^{\ell}_3}\ar[u]^{\eta^{\ell}_2} &
b^{\ell}_2\ar[r]_{\eta^{\ell}_4\eta^{\ell}_5}\ar[u]^{\eta^{\ell}_4} &
<EMAIL_ADDRESS>}\ar[u]^{\eta^{\ell}_6} &
b^{\ell}_j\ar[rr]_{\eta^{\ell}_{2j}\eta^{\ell}_{2j+1}}\ar[u]^{\eta^{\ell}_{2j}} &

[$\ell$ is even, $s_2 = s_1 + 1=\ell/2$.]

\xymatrix{
a^{\ell}_0\ar[r]^{\eta^{\ell}_1\eta^{\ell}_2}\ar[dr]_{\eta^{\ell}_1} &
a^{\ell}_1\ar[r]^{\eta^{\ell}_3\eta^{\ell}_4}\ar[dr]_{\eta^{\ell}_3} &
a^{\ell}_2\ar[r]^{\eta^{\ell}_5\eta^{\ell}_6}\ar[dr]_{\eta^{\ell}_5} &
<EMAIL_ADDRESS><EMAIL_ADDRESS>&
a^{\ell}_j\ar[rr]^{\eta^{\ell}_{2j+1}\eta^{\ell}_{2(j+1)}}\ar[drr]^{\eta^{\ell}_{2j+1}} &
b^{\ell}_0\ar[r]_{\eta^{\ell}_0\eta^{\ell}_1}\ar[u]^{\eta^{\ell}_0} &
b^{\ell}_1\ar[r]_{\eta^{\ell}_2\eta^{\ell}_3}\ar[u]^{\eta^{\ell}_2} &
b^{\ell}_2\ar[r]_{\eta^{\ell}_4\eta^{\ell}_5}\ar[u]^{\eta^{\ell}_4} &
<EMAIL_ADDRESS>}\ar[u]^{\eta^{\ell}_6} &
b^{\ell}_j\ar[rr]_{\eta^{\ell}_{2j}\eta^{\ell}_{2j+1}}\ar[u]^{\eta^{\ell}_{2j}} &

Iterative processes for calculating the coefficients $a^{\ell}_j$ and $b^{\ell}_j$ correspond to paths on graph.

With the help of $a^{\ell}_j$ and $b^{\ell}_j$, we can find the expression of $\beta(\ell,c)$ with fixed $\ell$. For $1\le \ell \le 6$, the $\beta(\ell,c)$ are illustrated in Table <ref>.
Obviously, $\beta(1,c)\equiv 1/2$ does not depend on the constant $c$. The series expansion  of $\beta(\ell,c)$ at $c=0$ are also demonstrated in Table <ref>.

Expressions of $\beta(\ell,c)$ for $1\le \ell\le 6$

$\ell$     $\beta(\ell,c)=\frac{n(s_1,-c)}{d(s_2,-c)}$
Taylor series of $\beta(\ell,c)$ at $c = 0$     Domain of convergence

$1$     $
\frac{1}{2}$        $\frac{1}{2}$      $[0, +\infty)$

$2$     $ \frac{\frac{1}{2}}{1-\frac{1}{12}c}$        $ \frac{1}{2} +  \frac{c}{24} + \frac{c^2}{288} + \frac{c^3}{3456} + \frac{c^4}{41472} + \frac{c^5}{497664} + \frac{c^6}{5971968} + \mathcal{O}(c^7)$      $[0, 12)$

$3$     $ \frac{\frac{1}{2}-\frac{1}{120}c}{1-\frac{1}{10}c}$        $\frac{1}{2}+\frac{c}{24}+\frac{c^2}{240}+\frac{c^3}{2400}+  \frac{c^4}{24000}
+   \frac{c^5}{240000} + \frac{c^6}{2400000}  +  \mathcal{O}\left(c^7\right)$     $[0,10)$

$4$     $ \frac{\frac{1}{2}-\frac{1}{84}c}{1-\frac{3}{28}c+\frac{1}{1680}c^2}$        $\frac{1}{2}+\frac{c}{24}+\frac{c^2}{240}+\frac{17 c^3}{40320}+\frac{241 c^4}{5644800}+ \frac{41c^5}{9483264}
+ \frac{29063c^6}{66382848000} +\mathcal{O}\left(c^7\right)$     $[0, \frac{12650}{1281})=[0,9.8751)$

$5$     $ \frac{\frac{1}{2}-\frac{1}{72}c+\frac
{1}{30240}c^2}{1 -\frac{1}{9}c + \frac{1}{1008} c^2}$        $\frac{1}{2}+\frac{c}{24}+\frac{c^2}{240}+\frac{17 c^3}{40320}+\frac{31 c^4}{725760}+\frac{1583 c^5}{365783040}+ \frac{2887c^6}{6584094720} +
\mathcal{O}\left(c^7\right)$     $[0, \frac{2349}{238})=[0,9.8697)$

$6$    $ \frac{\frac{1}{2}-\frac{1}{66}c + \frac{1}{15840}c^2}{1-\frac{5}{44}c+\frac{1}{792}c^2-\frac{1}{665280}c^3}$        $\frac{1}{2}+\frac{c}{24}+\frac{c^2}{240}+\frac{17 c^3}{40320}+\frac{31 c^4}{725760}+\frac{691 c^5}{159667200}+\frac{388151 c^6}{885194956800}+\mathcal{O}\left(c^7\right)$     $[0, \frac{10294}{1043})=[0,9.8696)$

§.§ Impacts of parameters $\ell$ and $c$ on $\beta(\ell,c)$

It is obvious that the transition matrix $\mat{G}^{k+1}_k(\ell,\tau)$ depends on the order parameter $\ell$ by the value of $\beta(\ell,c)$ essentially. Some expressions of $\beta(\ell,c)$ for $1\le \ell \le 6$ are demonstrated in Table <ref>. Note that for different $\ell$, the function $\beta(\ell,c)$ will have different domain of convergence (DOC) for parameter $c$.

Figure <ref> illustrates the $\beta(\ell,c)$ for different $\ell$ and $c$. Figure <ref> shows  that for $c\le 1.0\times 10^{-2}$ we have $\beta(\ell,c)\approx \beta(1,c) = \frac{1}{2}$. However, for $c\ge 0.1$, $\beta(\ell,c)$ varies with $c$ significantly. Moreover, for $\ell\ge 4$, $\beta(\ell,c)$ has  stable values for any $c$, which implies that $\ell \ge 4$ will be a satisfactory order parameter. Note that $\ell$ is a discrete variable and we use continuous lines for visualization in Figure <ref>.
On the other hand, Figure <ref> shows that for $\beta(\ell,c)$ differs a lot when $c > 2.0$. Moreover, for $\ell \ge 4$, the $\beta(\ell,c)$ coincides well for any $c$. In general, it is a good choice for us to set $\ell = 4$ if we prefer high order algorithm.

[$\beta(\ell,c)$ varies with $\ell$ for fixed $c$]

[$\beta(\ell,c)$ varies with $c$ for fixed $\ell$]

Impact of $\ell$ and $c$ on $\beta(\ell,c)$

§.§ $2\ell$-th order SDS  for LTI-QKDE

The diagonal Padè approximation of order $\ell$ for the transition matrix  for the LTI-QKDE  is

\begin{equation} \label{eq-gMat}
=  \phi\left( \beta(\ell,c)\frac{\tau}{2}  \mat{\Omega}\right) = \frac{\mat{I} + \beta(\ell,c)\frac{\tau}{2}  \mat{\Omega}}{\mat{I} - \beta(\ell,c)\frac{\tau}{2}  \mat{\Omega}}
\end{equation}

\begin{equation}
c = \norm{\vec{\omega}}^2\tau^2/4.
\end{equation}

Proof: This lemma holds obviously since

\begin{equation*}
\left(\tau\mat{\Omega}/2\right)^2 =  -c \mat{I}, \quad c = \norm{\vec{\omega}}^2\tau^2/4
\end{equation*}

and Lemma <ref> can be applied.

For any time step $\tau$ and time-invariant vector $\vec{\omega}$, the transition matrix  for the LTI-QKDE is

\begin{equation}  \label{eq-GtqExp}
\begin{split}
\mat{G}(\ell,\tau)
&= g_\ell\left(\frac{\tau}{2}\mat{\Omega}  \right)\\
&= \frac{1}{1+\alpha(\ell,c)}  \left[ \left( 1 - \alpha(\ell,c)\right)\mat{I}
+\tau \beta(\ell,c)\mat{\Omega}\right]\\
&= \cos\delta(\ell,\vec{\omega}, \tau)   \mat{I} + \sin \delta(\ell,\vec{\omega},\tau)  \hat{\mat{\Omega}},
\end{split}
\end{equation}

in which

\begin{align}
&\alpha(\ell, c) = c [\beta(\ell,c)]^2, \\
&\delta(\ell, \vec{\omega},\tau) = 2 \arctan (\beta(\ell,c)\norm{\vec{\omega}}\tau/2),     \\
&\hat{\mat{\Omega}} =\mat{\Omega}/\norm{\vec{\omega}}, \\
&\hat{\mat{\Omega}}^2 = -\mat{I}.
\end{align}

Proof: Let

\begin{align}
y = \tau\beta(\ell,c)/2, \quad
\gamma = \norm{\vec{\omega}}, \quad
\delta  = 2 \arctan (y\gamma),
\end{align}

then  $\mat{\Omega}^2 = -\gamma^2 \mat{I}$ and

\begin{equation}
\alpha(\ell,c) = y^2 \gamma^2 =  c [\beta(\ell,c)]^2.
\end{equation}

Let $\mat{M} = \mat{\Omega}$, then this theorem follows from the
Lemma <ref> directly. Q.E.D.

For any $\ell\in \mathbb{N}$, constant $\vec{\omega}\in \ES{R}{3}{1}$ and any $\tau\in \mathbb{R}$, the transition matrix $\mat{G}(\ell,\tau)$
is an orthogonal, invertible  and symplectic transformation with $2\ell$-th order accuracy, i.e.,

\begin{equation}\label{eq-Gq-1st-sp}
\inv{\mat{G}(\ell,\tau)} = \mat{G}(\ell,-\tau) = \trsp{\mat{G}(\ell,\tau)}
%&\trsp{\mat{G}(\ell,\tau)} \mat{J} \mat{G}(\ell,\tau)  = \mat{J} + \mathscr{O}(\sin \delta(\ell, \mat{\omega},\tau))
\end{equation}

for the LTI-QKDE.

Proof: For a time-invariant  vector $\vec{\omega}$, the function $ \delta(\ell,\vec{\omega},\tau)$ is an odd function of $\tau$. With the help of $ \hat{\mat{\Omega}}^2 = -\mat{I}$ and $\trsp{\hat{\mat{\Omega}}} = -\hat{\mat{\Omega}}$, we immediately obtain

\begin{align*}
\trsp{\mat{G}(\ell,\tau)}
&= \trsp{[\cos \delta(\ell,\vec{\omega},\tau)\mat{I} +\sin \delta(\ell,\vec{\omega},\tau)\hat{\mat{\Omega}}]}\\
&= \cos \delta(\ell,\vec{\omega},\tau)\mat{I} -\sin \delta(\ell,\vec{\omega},\tau)\hat{\mat{\Omega}} \\
&= \cos (-\delta(\ell,\vec{\omega},\tau))\mat{I} +\sin (-\delta(\ell,\vec{\omega},\tau))\hat{\mat{\Omega}}\\
&= \cos \delta(\ell,\vec{\omega},-\tau))\mat{I} +\sin \delta(\ell,\vec{\omega},-\tau)\hat{\mat{\Omega}}\\
&= \mat{G}(\ell,-\tau),
\end{align*}

which implies that the transition matrix is invertible. Moreover, $\mat{G}(\ell,\tau)$ is orthogonal by
Lemma <ref>. In consequence, $\trsp{\mat{G}(\ell,\tau)} = \mat{G}(\ell,-\tau) = \inv{\mat{G}(\ell,\tau)}$.

Since the $\mat{G}(\ell,\tau)$ is constructed from the Padè approximation $g_\ell(\cdot)$ directly according to Theorem <ref>, it must be a SDS with $2\ell$-th  order accuracy. Q.E.D.

§.§ $2\ell$-th order SDS for LTV-QKDE

For the time-varying vector $\vec{\omega}(t)$,  the LTV-QKDE is a non-autonomous Hamiltonian system essentially and the corresponding SGA  can be obtained via the concept of extended phase space and time-centered difference scheme  (See [30], Chapter 5).
Let $\bar{t}_k = t[k] +\tau/2,
\vec{\omega}_k =\vec{\omega}(\bar{t}_k),
\mat{\Omega}_k = \mat{\Omega}(\vec{\omega}_k)$ with the similar procedure as we do
in finding the transition matrix for SDS of the LTI-QKDE.
However, we have to treat the discrete time sequence $t_0, t_0 + \frac{\tau}{2}, t_0 + \frac{3\tau}{2}, \cdots, t_0 + \frac{(2k+1)\tau}{2}, t_0 + \frac{(2k+3)\tau}{2}, \cdots$ and the corresponding symplectic transition matrices $\mat{G}^{k+1}_k(\cdot)$ at these discrete time. Thus we have to treat the initial condition $\quat{q}(t_0)$ carefully since the length of interval $[t_0, t_0 + \frac{\tau}{2}]$ is not equal to that of $[t_0 + \frac{(2k+1)\tau}{2}, t_0 + \frac{(2k+3)\tau}{2}]$. In order to avoid this phenomena and the fractional sampling problem, we use the following modified version of the discrete  vectors $\vec{\omega}_k$ and matrices $\mat{\Omega}_k$

\begin{equation} \label{eq-omegaA-LTV}
\vec{\omega}_k =\vec{\omega}(t[k]), \quad \mat{\Omega}_k = \mat{\Omega}(\vec{\omega}_k),
\end{equation}

thus we can obtain the $k$-th transition matrix at time $t[k]$

\begin{equation}
\begin{split}
\mat{G}^{k+1}_k(\ell, \tau)
&= g_{\ell}\left(\tau\mat{\Omega}_k/2\right)
= \phi\left(\beta(\ell,c_k)\tau\mat{\Omega}_k/2 \right)
\end{split}
\end{equation}

for the LTV-QKDE. By Lemma <ref>, it is easy to prove the following theorem.

For the LTV-QKDE $\cfrac{\dif \quat{q}}{\dif t} = \frac{1}{2}\mat{\Omega}[\vec{\omega}(t)]\quat{q}$, let

\begin{align}
c_k &= \tau^2\norm{\vec{\omega}_k}^2/4, \\
\beta_k &= \beta(\ell,c_k), \\
\alpha_k &= \alpha(\ell,c_k) = c_k \beta_k^2,
\end{align}

then the transition matrices $\mat{G}^{k+1}_k(\ell,\tau): \quat{q}[k]\mapsto \quat{q}[k+1]$ can be expressed by

\begin{equation}
\begin{split}
&\mat{G}^{k+1}_k(\ell, \tau)\\
&= \frac{1}{1+\alpha(\ell,c_k)}\left[(1-\alpha(\ell,c_k)\mat{I} + \tau \beta(\ell,c_k)\mat{\Omega}_k\right]\\
&= \cos [\delta\left(\ell, \vec{\omega}_k,\tau\right) ] \mat{I}
+\sin [\delta\left(\ell, \vec{\omega}_k,\tau\right)] \hat{\mat{\Omega}}(\vec{\omega}_k),
\end{split}
\end{equation}

which generates a  SDS with $2\ell$-th order accuracy for $k= 0, 1, 2, \cdots$.

Note that  $\mat{G}^{k+1}_k(\ell, \tau)$ is also  symplectic and orthogonal. However,
$\mat{G}^{k+1}_k(\ell, -\tau) \neq \inv{[\mat{G}^{k+1}_k(\ell, \tau)]}$ since $\vec{\omega}(t)$ is time-dependent and $\vec{\omega}(t[k]-\tau) \neq \vec{\omega}(t[k]+\tau)$ in general.

§ SYMPLECTIC GEOMETRIC ALGORITHMS FOR QKDE

§.§ Algorithm for Computing the $\beta(\ell,c)$

For the practical computaiton of $\beta(\ell,c)$, it is necessary to find $\eta^{\ell}_k$ frequently. For the purpose of clearity, we design a simple procedure $\ProcName{Eta}$ to compute $\eta^{\ell}_k$, see <ref>.

Compute the Parameter $\eta^{\ell}_k$
[1]
Parameter $\ell\in \mathbb{N}$ for order and integer $k$.
The value of $\eta^{\ell}_k$.
$\eta^{\ell}_k \gets \cfrac{\ell-k}{(2\ell - k)(k+1)}$;

With the help of <ref> in Appendix <ref>, we can desgin two Fast and Stable Iterative Algorithms for generating $\beta(\ell,c)$ (FSIA), see <ref> and <ref>.

Parallel Fast and Stable Iterative Algorithm for Generating $\beta(\ell,c)$

[1]
Positive integer $\ell$ and auxiliary parameter $c$.
Key auxiliary parameter $\beta$.
PFsiaGenBeta$\ell, c$
$s_1 \gets \mfloor{(\ell-1)/2}, s_2 \gets \mfloor{\ell/2}$;
Allocate memory for  $s_1+1$ coefficients $a^{\ell}_j$.
Allocate memory for  $s_2+1$ coefficients $b^{\ell}_j$.
$a^{\ell}_0 \gets \frac{1}{2}$;
($ j\gets 0; j \le s_1-1; j \gets j + 1$)
$a^{\ell}_{j+1} \gets a^{\ell}_j \cdot\ProcName{Eta}(\ell, 2j+1)\cdot\ProcName{Eta}(\ell, 2j+2)$;
$b^{\ell}_0 \gets 1$;
($ j\gets 0; j \le s_2-1; j\gets j+1$)
$b^{\ell}_{j+1} \gets b^{\ell}_j \cdot \ProcName{Eta}(\ell,2j)\cdot \ProcName{Eta}(\ell, 2j+1)$;
$n \gets \ProcName{Polynomial}(a^{\ell}_0, \cdots, a^{\ell}_{s_1},-c)$;
$d \gets \ProcName{Polynomial}(b^{\ell}_0, \cdots, b^{\ell}_{s_2},-c)$;
$\beta \gets n/d$;

Alternative Fast and Stable Iterative Algorithm for Generating $\beta(\ell,c)$

[1]
Positive integer $\ell$ and auxiliary parameter $c$.
Key auxiliary parameter $\beta$.
AFsiaGenBeta$\ell, c$
$s_1 \gets \mfloor{(\ell-1)/2}, s_2 \gets \mfloor{\ell/2}$;
Allocate memory for  $s_1+1$ coefficients $a^{\ell}_j$.
Allocate memory for  $s_2+1$ coefficients $b^{\ell}_j$.
Set $b^{\ell}_0 \gets 1, a^{\ell}_0 \gets \frac{1}{2}$;
($ j\gets 0; j \le s_1-1; j \gets j + 1$)
$b^{\ell}_{j+1} \gets a^{\ell}_j \cdot\ProcName{Eta}(\ell, 2j+1)$;
$a^{\ell}_{j+1} \gets b^{\ell}_{j+1} \cdot\ProcName{Eta}(\ell, 2j+2)$;
($\ell \mod 2 \equiv 0$)
$b^{\ell}_{s_2} \gets a^{\ell}_{s_1}\cdot \ProcName{Eta}(\ell,2s_1+1)$;
$n \gets \ProcName{Polynomial}(a^{\ell}_0, \cdots, a^{\ell}_{s_1},-c)$;
$d \gets \ProcName{Polynomial}(b^{\ell}_0, \cdots, b^{\ell}_{s_2},-c)$;
$\beta \gets n/d$;

§.§ Algorithm for Computing the Symplectic Transition Matrix

The algorithm for computing the symplectic transition matrix $\mat{G}(\ell,\tau)$ is summarized in <ref>.

Compute the Symplectic Transition Matrix for LTV-QKDE
[1]
Angular velocity $\vec{\omega}\in \ES{R}{3}{1}$, order parameter $\ell$ and time step $\tau$.
Symplectic transition matrix $\mat{G}$ for the LTV-QKDE $\fracode{\quat{q}}{t} = \frac{1}{2}\mat{\Omega}(\vec{\omega}(t))\quat{q}$ at time $t$ with $2\ell$-th order accuracy.
SpTranMatQkdeLTV$\ell, \tau, \vec{\omega}$
\mat{\Omega}
\gets\begin{bmatrix}
0 & -\omega_1 & -\omega_2  & -\omega_3 \\
\omega_1 & 0 & \omega_3 & -\omega_2 \\
\omega_2 & -\omega_3 & 0 & \omega_1 \\
\omega_3 & \omega_2 & -\omega_1 & 0
\end{bmatrix};
$c \gets \tau^2\norm{\vec{\omega}}^2/4$;
$\beta \gets \ProcName{AFsiaGenBeta}(\ell,c)$;
$\alpha \gets c \beta^2$;
$\mat{G}(\ell, \tau) \gets \cfrac{(1-\alpha)\mat{I} + \tau\beta \mat{\Omega}}{1+\alpha}$;
$\mat{G}(\ell, \tau)$;

§.§ Algorithm for Solving the LTI-QKDE

The $2\ell$-th order ESGA for LTI-QKDE based on Padé approximation, $\ProcName{EoEsgaQkdeLTI}$ for
brevity[“Eo" means even order since $2\ell$ must be an even integer. We use EoEsgaQkdeLTI instead of EoEsgaLtiQkde since the EoEsgaQkdeLTI and its counterpart  EoEsgaQkdeLTV can be merged into EoEsgaQkde with function overloading in object-oriented programming.],  is given in Algorithm <ref> by the explicit expression of $\mat{G}(\ell, \tau)$  obtained.
Particularly, if $\ell = 1$, then the $\mat{G}(1, \tau)$ is the same with its counterpart in the ESGA-I algorithm proposed in [29] since  $\beta(1,c) \equiv \frac{1}{2}$ for any positive number $c$.

$2\ell$-th Order Explicit Symplectic Geometric Algorithm for the LTI-QKDE
[1]
The constant angular velocity vector $\vec{\omega}\in \ES{R}{3}{1}$, time step $\tau$, order parameter $\ell$, time interval $[t_0, t_f]$
and the initial state quaternion $\quat{q}_0 $
at time $t_0$.
Numerical solution to the LTI-QKDE
$\fracode{\quat{q}}{t} = \frac{1}{2}\mat{\Omega}(\vec{\omega})\quat{q}$ for $ t_0 \le t \le t_f$ with $2\ell$-th order SDS.
EoEsgaQkdeLTI$\ell, \tau, \vec{\omega}, \quat{q}_0, t_0, t_f$
$t[0] \gets t_0$,
$\quat{q}[0] \gets \quat{q}_0, n \gets \mfloor{(t_f-t_0)/\tau}$;
$\mat{G}(\ell,\tau) = \ProcName{SpTranMatQkdeLTV}(\ell, \tau, \vec{\omega})$;
($k\gets 0; k < n; k \gets k + 1$)
$\quat{q}[k+1] \gets \mat{G}(\ell, \tau) \quat{q}[k]$;
$t[k+1] \gets t[k] + \tau$;
The sequence $\set{\quat{q}[k]}$;

Note that for computing $\beta(\ell,c)$ in Step 6 of <ref>,  we can use $\ProcName{PFsiaGenBeta}(\ell,c)$ or $\ProcName{AFsiaGenBeta}$ for $\ProcName{FsiaGenBeta}$ according to the practical engineering requirements. Moreover, if a fixed $\ell$ is preferred, say $\ell = 4$, we can choose a specific polynomial from Table I. This strategy also holds for <ref> in the following subsection.

We remark that for real-time applications, it is not necessary to store the sequences $\set{t[k]}$ and $\set{\quat{q}[k]}$. It is enough for us to store the current state $\quat{q}[k]$ at time $t[k]$, which implies that we just need to store two real numbers.  Obviously, this remark holds for the <ref> EoEsgaQkdeLTV in the following subsection.

§.§ Algorithm for Solving the LTV-QKDE

The even order ESGA for LTV-QKDE based on Padè approximation, EoEsgaQkdeLTV for brevity, is presented in <ref>. It should be noted that we can share some memories for the
time-dependent parameters $\vec{\omega}_k, \mat{\Omega}_k, c_k, \beta_k, \alpha_k$ and $\mat{G}^{k+1}_k(\ell,\tau)$ since they are local variables for omputing the sequence $\set{\quat{q}[k]}$.

By comparison, EoEsgaQkdeLTV differs the EoEsgaQkdeLTI from two points: the angular velocity at time $t[k]$ should be computed and the computation of $\mat{G}$ is moved into the loop.

$2\ell$-th Order Explicit Symplectic Geometric Algorithm for the LTV-QKDE
[1]
The time-varying vector $\vec{\omega}(t)\in \ES{R}{3}{1}$, positive integer $\ell$, time interval $[t_0,t_f]$, initial quaternion
$\quat{q}_0$  and time step $\tau$.
Numeric solution $\set{\quat{q}[k]}$ to the LTV-QKDE
$\fracode{\quat{q}}{t} = \frac{1}{2}\mat{\Omega}(\vec{\omega}(t))\quat{q}$ for $ t_0\le t\le t_f$
with $2\ell$-th order accuracy.
EoEsgaQkdeLTV$\ell, \tau, \vec{\omega}(t), \quat{q}_0,t_0, t_f$
$t[0] \gets t_0, \quat{q}[0] \gets \quat{q}_0, n\gets \mfloor{(t_f-t_0)/\tau}$;
($k\gets 0; k < n; k \gets k + 1$)
$\vec{\omega}_k \gets \trsp{[\omega_1(t[k]), \omega_2(t[k]), \omega_3(t[k])]}$;
$\mat{G}^{k+1}_k(\ell,\tau) = \ProcName{SpTranMatQkdeLTV}(\ell, \tau, \vec{\omega}_k)$;
$\quat{q}[k+1] \gets \mat{G}^{k+1}_k(\ell, \tau) \quat{q}[k]$;
$t[k+1] \gets t[k] + \tau$;
The sequence $\set{\quat{q}[k]}$;

Note that for the steps in the loop, we can use the notations
$\vec{\omega}_k$ and $\mat{G}^{k+1}_k$ for understanding (see Theorem <ref>). For the implementation of $\ProcName{EoEsgaQkdeLTV}$, we can use
the notations $\vec{\omega}$ and $\mat{G}$ in order to share memories for the tempary (local) variables in the sense of computer programming.

§.§ Analysis of Computational Complexity

$\tic{*}{(\cpvar{expr})}$ and
$\tic{+}{(\cpvar{expr})} $
be the times of multiplication and addition in some operation expression $\cpvar{expr}$. The time  complexity vector of computation (TCVC)  for $\cpvar{expr}$ is defined by

\begin{equation}
\timer(\cpvar{expr}) = [\tic{*}{(\cpvar{expr})}, \tic{+}{(\cpvar{expr})}].
\end{equation}

Note that we just list two components of $\timer(\cpvar{expr})$ here since for the problems and algorithms which contain massive matrix-vector operations, the multiplication and addition operations for real numbers are fundamental and essential.

Similarly, we use $T(\cpvar{expr})$ to denote the computation time for $\cpvar{expr}$. We also use $T_*(\cpvar{expr})$ and $T_+(\cpvar{expr})$ to represent the computation time for the multiplications and additions involved in $\cpvar{expr}$.
Suppose that the time units for multiplication and addition are  $u_1$ and $u_2$ respectively. Let

\begin{equation}
\vec{u} = [u_1, u_2]
\end{equation}

be the vector of time units. Then the time for computing $\cpvar{expr}$ will be

\begin{equation}
\begin{split}
& = \braket{\timer(\cpvar{expr})}{\vec{u}} \\
%&= [T^*(\cpvar{expr}), T^+(\cpvar{expr})]\begin{bmatrix}
%u_1 \\
%u_2\end{bmatrix} \\
&= \tic{*}{(\cpvar{expr})} u_1 + \tic{+}{(\cpvar{expr})}u_2\\
&= T_*(\cpvar{expr}) + T_+(\cpvar{expr})
\end{split}
\end{equation}

if only multiplication and addition are essential for the total time consumed.

For an algorithm named with $\ProcName{Alg}$, its TCVC is defined by

\begin{equation}
\begin{split}
\tic{}{(\ProcName{Alg})}
&= \sum_{\cpvar{expr}\in \ProcName{Alg}} \tic{}{(\cpvar{expr})}\\
&=[\tic{*}{(\ProcName{Alg})},\tic{+}{(\ProcName{Alg})} ].
\end{split}
\end{equation}

The time needed for  algorithm Alg can be measured by

\begin{equation}
T(\ProcName{Alg}) = \braket{\tic{}{(\ProcName{Alg})}}{\vec{u}}= \tic{*}{(\ProcName{Alg})} u_1 + \tic{+}{(\ProcName{Alg})}u_2
\end{equation}

theoretically with acceptable accuracy.

It is easy to check that the TCVC for Algorithm <ref> for computing $\eta^{\ell}_k= \ProcName{Eta}(\ell,k)$ is

\begin{equation}
\tic{}{(\ProcName{Eta})}
= \sum_{ \cpvar{expr}\in \ProcName{Eta} }
\tic{}{(\cpvar{expr})} =  [3, 3].
\end{equation}

The TCVC for computing the polynomial $p_s(x) = \sum^s_{i=0}c_i x^i$ via Horne's rule is

\begin{equation}
\tic{}{(\ProcName{Polynomial})}
= \sum_{ \cpvar{expr}\in \ProcName{Polynomial} }
\tic{}{(\cpvar{expr})} =  [s, s]
\end{equation}

Simple algebraic operations show that The TCVC for <ref> and <ref> are

\begin{equation}
\begin{split}
&\tic{}{(\ProcName{PFsiaGenBeta})}  \\
&= \sum_{ \cpvar{expr}\in \ProcName{PFsiaGenBeta} }
\tic{}{(\cpvar{expr})}  \\
&= s_1[11,9] + s_2[11,8] + [5,3]  \\
&= [11(s_1 + s_2)+5, 9s_1 + 8s_2 + 3]  \\
\begin{array}{ll}
\left[11\ell-6, (17\ell-11)/2 \right],   &  \text{for}~\ell~\text{is~odd};\\
\left[11\ell-6, (17\ell-12)/2 \right],       &  \text{for}~\ell~\text{is~even};
\end{array}
\right.\\
&= [11\ell -6, (17\ell -12 + \mod(\ell,2))/2], \quad \forall \ell \in \mathbb{N}
\end{split}
\end{equation}

\begin{equation}
\begin{split}
&\tic{}{(\ProcName{AFsiaGenBeta})}  \\
&= \sum_{ \cpvar{expr}\in \ProcName{AFsiaGenBeta} }
\tic{}{(\cpvar{expr})}  \\
&=[6\ell-1, 5\ell-3], \quad \forall \ell\in \mathbb{N}
\end{split}
\end{equation}

respectively. Obviously, the time complexity of $\ProcName{AFsiaGenBeta}$ is lower than that of $\ProcName{PFsiaGenBeta}$ since there is only one loop in the former. However, if we take parallel programming, the two loops in PFsiaGenBeta can be executed at the same time. In this way, the difference of the time complexity  will be ignorable.

The TCVC for $\ProcName{SpTranMatQkde}$ can be expressed by

\begin{equation}
\begin{split}
&\tic{}{(\ProcName{SpTranMatQkde})}  \\
&= \sum_{ \cpvar{expr}\in \ProcName{SpTranMatQkde} }
\tic{}{(\cpvar{expr})}  \\
&=[6\ell+29, 5\ell+6], \quad \forall \ell\in \mathbb{N}.
\end{split}
\end{equation}

This implies that the computational cost is constant for fixed $\ell$ and varies with $\ell$ linearly when $\ell$ increases.

Time and Space Consumptions of EoEsgaQkdeLTV

Step     TCVC      Numbers/matrices to be stored

2         —           Setting initial values

3         —           $k\in \mathbb{Z}, t[k]\in \mathbb{R}$

4         —          $\vec{\omega}\in \ES{R}{3}{1}$

5        $ [6\ell+29, 5\ell+6]$       $c, \beta, \alpha\in \mathbb{R}, \mat{G}\in \ES{R}{4}{4}$

6        $ [16,12] $       $\quat{q}[k+1]\in \ES{R}{4}{1}$

7       $[0,1]$     —

8       —     —

9       —     —

\begin{equation}
n = \mfloor{(t_f-t_0)/\tau},
\end{equation}

then there are $n$ iterations for the loops in <ref> and  <ref>. Table <ref> shows that the
TCVC for EoEsgaQkdeLTV is

\begin{equation}
\begin{split}
&=\sum_{\cpvar{expr}\in \ProcName{EoEsgaQkdeLTV}} \tic{}{(\cpvar{expr})}\\
&= [6\ell n + 45n, 5\ell n + 19n]
\end{split}
\end{equation}

and the total time consumed (without optimization) can be measured by

\begin{equation}
\begin{split}
&= (6\ell n +45n )u_1 + n(5\ell n+19n)u_2\\
%&= ((6u_1 + 5u_2)\ell + (45u_1 + 19u_2)) n  \\
&= \mathcal{O}(n)
\end{split}
\end{equation}

for fixed order parameter $\ell$, time units $u_1$ and $u_2$.
By comparison, the TCVC for EoEsgaQkdeLTI is

\begin{equation}
\begin{split}
&=\sum_{\cpvar{expr}\in \ProcName{EoEsgaQkdeLTI}} \tic{}{(\cpvar{expr})}\\
&= [16n+6\ell+29, 13n+ 5\ell+6]
\end{split}
\end{equation}

and the time consumed will be

\begin{equation}
\begin{split}
&= (16n+6\ell+29)u_1 + (13n+ 5\ell+6)u_2\\
&= \mathcal{O}(n).
\end{split}
\end{equation}

Therefore, the time complexity for both $\ProcName{EoEsgaQkdeLTI}$ and $\ProcName{EoEsgaQkdeLTV}$ are linear.
Moreover, it is necessary to store the global variables $t_0, t_f, \tau$ and local variables $\omega_1, \omega_2, \omega_3,  \mat{\Omega},  \alpha, \beta, \mat{G}, k$ and $\quat{q}[k]$. The storage sharing mechanics implies that  only $45$ real numbers as well as a pointer to $\vec{\omega}(t)$ should be stored in this algorithm for realtime applications. Obviously, the space complexity is  constant, i.e., $\mathcal{O}(1)$.

We remark that we can set $\beta(\ell,c)$ according to the Table <ref> without invoking the <ref> or <ref>. By this way,  the time  complexity of <ref> and <ref> will still be linear.

§ PERFORMANCE EVALUATION

The fundamental aspects of verification and validation of ESGA for QKDE are discussed in detail in [29]. Our emphasis here lies in the following issues:

* Analysis of accuracy of $\ProcName{EoEsgaQkdeLTI}$;

* Verification and evaluatation of  $\ProcName{EoEsgaQkdeLTV}$;

* Comparison of the NS and AS for some special time-varying $\vec{\omega}(t)$ and $\quat{q}[0]$;

* Verification of $\ProcName{EoEsgaQkdeLTV}$ with MATLAB Simulink for general case.

§.§ Measure of NS and AS

Let $\mat{G}_{\mathrm{AS}}(k,\tau)$ and $\mat{G}_{\mathrm{NS}}(k,\tau)$ be the AS and NS to the single step transition matrix at time $t[k]$ respectively.
Obviously, we have

\begin{equation}
\mat{G}_{\mathrm{AS}}(k,\tau) = \mat{G}^{k+1}_k(\tau), \quad\mat{G}_{\mathrm{NS}}(k,\tau) = \mat{G}^{k+1}_k(\ell,\tau).
\end{equation}

The error operator of step $k$ is defined by

\begin{equation}
\begin{split}
\mat{E}_\mathcal{G}[k]
&= \mat{G}_{\mathrm{AS}}(k,\tau) - \mat{G}_{\mathrm{NS}}(k,\tau) \\
&= \mat{G}^{k+1}_k(\tau) - \mat{G}^{k+1}_k(\ell,\tau).
\end{split}
\end{equation}

For the fixed $\tau
$ and  $\ell$, the absolute error of the NS  for the $k$-th solution   $\quat{q}[k]$ is

\begin{equation}
\begin{split}
&=  \norm{\quat{q}^{\mathrm{NS}}[k] - \quat{q}^{\mathrm{AS}}[k]}\\
\end{split}
\end{equation}

The maximum absolute error for solving $\quat{q}(t)$ on interval $[t_0, t_f]$ can be expressed by

\begin{equation}
\begin{split}
&= \max_{t[k]\in [t_0, t_f]} E_k(\ell,\tau) =  \\
&=  \max_{0\le k \le \mfloor{(t_f-t_0)/\tau}} \norm{\quat{q}^{\mathrm{NS}}[k] - \quat{q}^{\mathrm{AS}}[k]}.
\end{split}
\end{equation}

§.§ Analysis of Accuracy of $\ProcName{EoEsgaQkdeLTI}$ for the LTI-QKDE

For the LTI-QKDE, we can obtain the AS to
the equation (<ref>) for constant $\vec{\omega}$. Since $\mat{\Omega}^2 = -\norm{\vec{\omega}}^2 \mat{I}$ and $\hat{\mat{\Omega}} = \mat{\Omega}/\norm{\vec{\omega}}$, Lemma <ref> shows that

\begin{align*}
\quat{q}(t[k+1])
&= \exp\left((t[k+1]-t[k])\mat{\Omega}/2\right)\quat{q}(t[k]) \\
&= \left(\mat{I} \cos (\norm{\vec{\omega}}\tau/2)
+ \hat{\mat{\Omega}} \sin (\norm{\vec{\omega}}\tau/2)\right)\quat{q}(t[k]).
\end{align*}

\begin{equation}
\quat{q}[k+1] = \mat{G}_{\mathrm{AS}}(k,\tau)\quat{q}[k]
\end{equation}

\begin{equation}
\begin{split}
\mat{G}_{\mathrm{AS}}(k,\tau)
&=\mat{G}^{k+1}_k(\tau) =\mat{G}(\tau)
= \exp\left(\tau\mat{\Omega}/2\right)\\
&= \mat{I} \cos (\norm{\vec{\omega}}\tau/2)
+ \hat{\mat{\Omega}} \sin (\norm{\vec{\omega}}\tau/2)
\end{split}
\end{equation}

is the AS to the transition matrix for any $k$. On the other hand,

\begin{equation}
\mat{G}^{k+1}_k(\ell,\tau) = \mat{G}(\ell,\tau)
\end{equation}

for the LTI-QKDE. Consequently,

\begin{equation}
\begin{split}
\mat{E}_\mathcal{G}
%&= \mat{G}_{\mathrm{AS}}(\tau) - \mat{G}_{\mathrm{NS}}(\tau) \\
%&= \mat{G}^{k+1}_k(\tau) - \mat{G}^{k+1}_k(\ell,\tau)\\
&= \mat{G}(\tau) - \mat{G}(\ell,\tau)\\
&= \mat{I}\left[ \cos \frac{\norm{\vec{\omega}}\tau}{2} -
\frac{1-\alpha(\ell,c)}{1+\alpha(\ell,c)}\right]\\
+ \hat{\mat{\Omega}} \left[
\sin  \frac{\norm{\vec{\omega}}\tau}{2} - \frac{\norm{\vec{\omega}}\tau\beta(\ell,c)}{1+\alpha(\ell,c)}
\right]
\end{split}
\end{equation}

for any $k$.

\begin{equation}
x = \norm{\vec{\omega}}\tau
\end{equation}

then  $x$ must be a pure scalar parameter without any unit and we immediately have

\begin{equation}
\begin{split}
c &= x^2/4,\\
\beta(\ell,c)&= \beta(\ell,x^2/4),\\
\alpha(\ell,c)&= c\beta(\ell,x)^2 = x^2[\beta(\ell,x^2/4)]^2/4.
\end{split}
\end{equation}

\begin{equation}
\left\{
\begin{split}
f_1(\ell,x) &= \cos \frac{x}{2} -
\frac{1-\frac{x^2}{4}\left[\beta\left(\ell,\frac{x^2}{4}\right)\right]^2}{1+\frac{x^2}{4}\left[\beta\left(\ell,\frac{x^2}{4}\right)\right]^2},\\
f_2(\ell,x) &= \sin \frac{x}{2} -
\frac{x\beta\left(\ell, \frac{x^2}{4}\right)}{1+\frac{x^2}{4}\left[\beta\left(\ell,\frac{x^2}{4}\right)\right]^2},
\end{split}
\right.
\end{equation}

\begin{equation}
\mat{E}_\mathcal{G} = f_1(\ell,x)\mat{I}
+ f_2(\ell,x)\hat{\mat{\Omega}}.
\end{equation}

\begin{equation}\label{eq-ErrMatNormEst}
\begin{split}
\norm{\mat{E}_\mathcal{G}}
\norm{f_1(\ell,x)\mat{I}
+ f_2(\ell,x)\hat{\mat{\Omega}}} \\
&\le \norm{f_1(\ell,x)\mat{I}}
+ \norm{f_2(\ell,x)\hat{\mat{\Omega}}}\\
%&= \abs{f_1(\ell,x)}\norm{\mat{I}} +
%\abs{f_2(\ell,x)} \norm{\hat{\mat{\Omega}}}\\
&=\abs{f_1(\ell,x)}\norm{\mat{I}} +
\abs{f_2(\ell,x)} \norm{\mat{\Omega}}/\norm{\vec{\omega}}\\
&= 2(\abs{f_1(\ell,x)} + \abs{f_2(\ell,x)})\\
\end{split}
\end{equation}

for sufficiently small $x = \norm{\vec{\omega}}\tau$ and any $k$.
According to Table <ref> and (<ref>), we can obtain Table <ref>. It is clear that the upper bound $2(\abs{f_1(\ell,x)} + \abs{f_2(\ell,x)})$ for
$\norm{\mat{E}_\mathcal{G}} $ decreases rapidly when $\ell$ increases.

Furthermore, for the LTI-QKDE we have the single step computational error

\begin{align*}
&\quat{q}^{\mathrm{AS}}[k+1] - \quat{q}^{\mathrm{NS}}[k+1]\\
&=\mat{G}_{\mathrm{AS}}(\tau)\quat{q}^{\mathrm{AS}}[k] - \mat{G}_{\mathrm{NS}}(\tau)\quat{q}^{\mathrm{NS}}[k]\\
&=\left(\mat{G}_{\mathrm{AS}}(\tau) - \mat{G}_{\mathrm{NS}}(\tau)\right)\quat{q}^{\mathrm{AS}}[k]+\mat{G}_{\mathrm{NS}}(\tau)\left(\quat{q}^{\mathrm{AS}}[k]-\quat{q}^{\mathrm{NS}}[k]\right)
\end{align*}

and its estimation

\begin{align*}
&\norm{\quat{q}^{\mathrm{AS}}[k+1] - \quat{q}^{\mathrm{NS}}[k+1]}\\
&\le \norm{(\mat{G}_{\mathrm{AS}}(\tau) - \mat{G}_{\mathrm{NS}}(\tau))\quat{q}^{\mathrm{AS}}[k]} + \norm{\mat{G}_{\mathrm{NS}}(\tau)(\quat{q}^{\mathrm{AS}}[k]-\quat{q}^{\mathrm{NS}}[k])}\\
&\le \norm{\mat{G}_{\mathrm{AS}}(\tau) - \mat{G}_{\mathrm{NS}}(\tau)}\norm{\quat{q}^{\mathrm{AS}}[k]}
&=\norm{\mat{G}_{\mathrm{AS}}(\tau) - \mat{G}_{\mathrm{NS}}(\tau)} + \norm{\quat{q}^{\mathrm{AS}}[k]-\quat{q}^{\mathrm{NS}}[k]}\\
&= \norm{\mat{E}_\mathcal{G}} + \norm{\quat{q}^{\mathrm{AS}}[k]-\quat{q}^{\mathrm{NS}}[k]}
\end{align*}

since we have $\norm{\quat{q}^{\mathrm{AS}}[k]} \equiv 1$ for any  integer $k$ and $\mat{G}_{\mathrm{NS}}(\tau)$ is an orthogonal matrix.  For $k=0$, we have

E_0(\ell,\tau) = \norm{\quat{q}^{\mathrm{AS}}[0]-\quat{q}^{\mathrm{NS}}[0]} = 0

because of $\quat{q}^{\mathrm{AS}}[0] =\quat{q}^{\mathrm{NS}}[0] = \quat{q}[0]$.
In consequence,

\begin{equation}
\begin{split}
&= \norm{\quat{q}^{\mathrm{AS}}[k]-\quat{q}^{\mathrm{NS}}[k]} \\
&\le \norm{\mat{E}_\mathcal{G}} + E_{k-1}(\ell,\tau)\\
&\le k \norm{\mat{E}_\mathcal{G}}
\end{split}
\end{equation}

for any $k$ by induction. Therefore,

\begin{equation}
\begin{split}
&\le 2 k_{max} (\abs{f_1(\ell,x)} + \abs{f_2(\ell,x)})\\
& = 2 \mfloor{\frac{t_f - t_0}{\tau}} (\abs{f_1(\ell,x)} + \abs{f_2(\ell,x)}) \\
\end{split}
\end{equation}

This implies that the maximum absolute error $E_{max}(\ell,\tau)$ has an upper bound with order $\mathcal{O}(\tau^{2\ell})$ for fixed length of time.

Accuracy of $\mat{G}(\ell,\tau)$ for LTI-QKDE: $\norm{\mat{E}_\mathcal{G}} \le 2(\abs{f_1(x)} + \abs{f_2(x)})$

$\ell$      $f_1(\ell,x)$      $f_2(\ell,x)$      $2(\abs{f_1(x)} + \abs{f_2(x)})=\mathcal{O}(x^{2\ell+1})$

$ 1$        $-\frac{x^4}{192}+\mathcal{O}\left(x^6\right)$        $\frac{x^3}{96}+\mathcal{O}\left(x^5\right)$       $\frac{x^3}{48} + \mathcal{O}\left(x^4\right)=\mathcal{O}(x^{2\times 1 +1})$

$ 2$        $-\frac{x^6}{46080}+\mathcal{O}\left(x^7\right)$        $\frac{x^5}{23040}+\mathcal{O}\left(x^7\right)$       $\frac{x^5}{11520}+  \mathcal{O}\left(x^6\right)=\mathcal{O}(x^{2\times 2 +1})$

$ 3$         $-\frac{x^8}{25804800}+\mathcal{O}\left(x^{10}\right)$       $\frac{x^7}{12902400}+\mathcal{O}\left(x^9\right)$       $\frac{x^7}{6451200}+\mathcal{O}\left(x^8\right)=\mathcal{O}(x^{2\times 3 +1})$

$ 4$         $-\frac{x^{10}}{26011238400}+\mathcal{O}\left(x^{12}\right)$       $\frac{x^9}{13005619200}+\mathcal{O}\left(x^{11}\right)$       $\frac{x^9}{6502809600} + \mathcal{O}\left(x^{10}\right)=\mathcal{O}(x^{2\times 4 +1})$

$ 5$         $-\frac{x^{12}}{41201801625600}+\mathcal{O}\left(x^{13}\right)$       $\frac{x^{11}}{20600900812800}+\mathcal{O}\left(x^{12}\right)$       $\frac{x^{11}}{10300450406400} + \mathcal{O}\left(x^{12}\right)=\mathcal{O}(x^{2\times 5 +1})$

For the order parameter $\ell \in \set{1, 2, \cdots, 10}$ and  time step
$\tau \in \set{10^{-6}, 10^{-5}, 10^{-4}, 10^{-3}, 10^{-2},10^{-1},  0.2, 0.4, 0.6, 0.8}$ (seconds), Figure <ref> demonstrates the impacts of $\ell$ and $\tau$ on the maximum absolute error $E_{max}(\ell,\tau)$ for the LTI-QKDE.

[Effect of $\ell$ on $E_{max}(\ell,\tau)$ for fixed time step $\tau$.]

[Effect of $\tau$ on $E_{max}(\ell,\tau)$ for fixed order parameter $\ell$.]

$E_{max}(\ell,\tau)$ of LTI-QKDE for different order and time step with a long time span $[t_0, t_f] = [0, 2000]$, constant angular velocity $\vec{\omega} = \trsp{[\pi\sin \frac{\pi}{8}, -\frac{\pi}{3}\cos\frac{\pi}{8}, -2\sin\frac{\pi}{3}]}$   and initial state $\quat{q}_0=\trsp{[1,0,0,0]}$.

Figure <ref> shows that if $\tau\le 10^{-3}$ second, we can set $\ell \ge  2$ since we have $E_{max}(\ell,\tau)\sim 10^{-13}$ for $\ell \ge 2$ instead of $E_{max}(\ell,\tau)\sim 10^{-7}$ for $\ell = 1$. If we set $E_{max}(\ell,\tau)\le 10^{-5}$, we should set $ \tau \le 0.007$ second for $\ell =1$ and $\tau \le 0.11$ for $\ell\ge 2$. Moreover, for $\ell \ge 4$, $E_{max}(\ell,\tau)$ will be small enough which implies the time step can be set large enough ($0< \tau < 1$) and the requirement for the performance of hardware will be easily satisfied.

§.§ Simulation of $\ProcName{EoEsgaQkdeLTV}$ for the LTV-QKDE

§.§.§ A special LTV-QKDE with AS

For the constant parameters
$\omega^{\ell}_0, \xi$, initial state

$$\quat{q}[0] = \trsp{[\cos (\xi/2), 0, \sin(\xi/2), 0]}$$

and  time-varying angular velocity

\begin{equation}  \label{eq-omega-special}
\vec{\omega}(t) = \begin{bmatrix}
-\omega^{\ell}_0(1-\cos \xi)\\
-\omega^{\ell}_0 \sin \xi \sin (\omega^{\ell}_0 t)\\
\omega^{\ell}_0 \sin \xi \cos (\omega^{\ell}_0 t)
\end{bmatrix},
\end{equation}

it is easy to check that the corresponding AS to the QKDE is

\begin{equation*}
\quat{q}^{\mathrm{AS}}(t) = \begin{bmatrix}
e_0(t)\\ e_1(t) \\ e_2(t) \\ e_3(t) \end{bmatrix}=
\begin{bmatrix}
\cos(\xi/2) \\  0 \\ \sin (\xi/2) \cos(\omega^{\ell}_0 t) \\
\sin(\xi/2)\sin(\omega^{\ell}_0 t)
\end{bmatrix}.
\end{equation*}

For the discrete version, we have

$$\quat{q}^{\mathrm{AS}}[k] =\quat{q}^{\mathrm{AS}}(t_0 + k\tau) = \mat{G}^k_{k-1}(\tau)\quat{q}^{\mathrm{AS}}[k-1].$$

On the other hand, the NS by the EoEsgaQkdeLTV is

\begin{equation}
\quat{q}^{\mathrm{NS}}_{\ell,\tau}[k]
= \begin{bmatrix}
e^{\mathrm{NS}}_0[k]\\ e^{\mathrm{NS}}_1[k] \\ e^{\mathrm{NS}}_2[k] \\ e^{\mathrm{NS}}_3[k]
\end{bmatrix}
= \mat{G}^{k}_{k-1}(\ell,\tau)\quat{q}^{\mathrm{NS}}_{\ell,\tau}[k-1].
%e_0(t_0+ k\tau)\\ e_1(t_0+ k\tau) \\ e_2(t_0+ k\tau) \\ e_3(t_0+ k\tau) \end{bmatrix}
\end{equation}

Effect of $\ell$ and $\tau$ on $E_{max}(\ell,\tau)$ for LTV-QKDE

Figure <ref> demonstrates a numeric example of LTV-QKDE with AS. Here we set    the time span $t_f - t_0 = 2000$ seconds, $\omega^{\ell}_0 = 2\pi$,$\xi = \pi/80$, the initial state $\quat{q}_0=\trsp{[\cos\frac{\xi}{2},  0 ,\sin\frac{\xi}{2} \cos(\omega^{\ell}_0 t), \sin\frac{\xi}{2}\sin(\omega^{\ell}_0 t)]}$ and $\vec{\omega}(t)$ according to (<ref>). Obviously, larger order parameter $\ell$ implies smaller $E_{max}(\ell,\tau)$ for fixed $\tau$.
If we set $\tau = 10^{-2}$ and need  $E_{max}(\ell,\tau)\sim 10^{-5}$, then any $\ell \ge 1$ is acceptable.
However, if we set $\tau = 10^{-1}$ and $E_{max}(\ell,\tau)\le 10^{-5}$, then we shoud set $\ell \ge 2$.
For $\ell \ge 4$ and  $E_{max}(\ell,\tau) \le 10^{-4}$,
a big step $\tau = 0.8$ is satisfactory. The figure shows that for fixed $E_{max}(\ell,\tau)$, we can set large step for high order algorithms. In other words, with configurable order parameter $\ell$ we can choose proper hardware for practical flight control systems by setting the step $\tau$ for sampling.  Moreover, for fixed $\tau$, a larger $\ell$ means better manoeuvrability  in flight since it traces the variation of $\vec{\omega}(t)$ better. Obviously, this observation  may be valuable for designing aircrafts, spacecrafts, torpedoes and so on.

§.§.§ Simulation of LTV-QKDE with MATLAB Simulink

It is well known that there is no AS for a LTV system such as the LTV-QKDE. In addition to the NS, simulation is another method to explore the solution to  LTV-QKDE.

Figure <ref> illustrates the diagram for simulating LTV-QKDE with MATLAB Simulink.

[MATLAB Simulink diagram.]

[Error between Pseudo AS to LTV-QKDE and NS with MATLAB Simulink]

$\|\quat{q}^{sim}-\quat{q}^{NS}\|$ of LTV-QKDE for $\ell \in \{1,2,3,4,5\}$ and time step $\tau = 10^{-5}$
with a time span $[t_0, t_f] = [0, 100]$,  time-varying angular velocity $\vec{\omega}(t) = \trsp{[-\omega^{\ell}_0\cos(\xi t)\exp(-\omega^{\ell}_0t),
\omega^{\ell}_0\cos(\xi t) \cos(\omega^{\ell}_0 t)]}$, and initial state $\quat{q}_0=\trsp{[\cos(\xi/2), 0 ,\sin(\xi/2),0]}$ where $\omega^{\ell}_0 = 2\pi$ and $\xi = \pi/80$.

§ CONCLUSIONS

The essential observations for solving the QKDE with high order explicit symplectic geometric algorithms include four aspects:

* the general LTV-QKDE can be modeled by the linear Hamiltonian system with infinitesimal symplec structure;

* the symplectic Padè approximation can be used to construct the SDS for QKDE with $2\ell$-th order accuracy;

* the symplectic transform (single-step transition operator) can be simplified by the Padè-Cayley lemma;

* the parameters for the SDS can be obtained with parallel and alternative iterative methods;

The performance of algorithms — EoEsgaQkdeLTV and EoEsgaQkdeLTI —  are attractive due to the following reasons:

* there are no accumulative computational errors in the sense of long term time since they are symplectic;

* the order parameter $\ell$ is configurable, which benifits  the designer of practical enginnering systems such as aircrafts, spacecrafts and so on according to the acceptable angular velocity;

* the time  complexity of computation for the algorithms proposed   is linear while the space complexity of computation is constant, thus our algorithms are appropriate for realtime applications;

* the maximum absolute error when solving the LTI-QKDE is  controlled by the time span $t_f-t_0$, time step $\tau$ and the order parameter $\ell$ through the relation $E_{max}(\ell,\tau) = \mathcal{O}((t_f-t_0)\tau^{2\ell})$, and the simulation result shows that this relation also holds for some LTV-QKDE.

Since the LTI-QKDE is a special case of LTV-QKDE and the EoEsgaQkdeLTV surpasses the algorithms ESGA-I and ESGA-II proposed in [29], it is enough for us to use the EoEsgaQkdeLTV algorithm for solving QKDE.

§ ACKNOWLEDGEMENT

This work was supported by the
Hainan Provincial Natural Science Foundation of China under grant number 2019RC199.

§ HAMILTONIAN SYSTEM

W. R. Hamilton introduced the canonical differential equations [35, 36]

\begin{equation*}
\fracode{p_i}{t} = -\fracpde{H}{q_i}, \quad \fracode{q_i}{t} = \fracpde{H}{p_i}, \quad i =1, 2, \cdots, N
\end{equation*}

for problems of geometrical optics, where $p_i$ are the generalized momentums,
$q_i$ are the generalized displacements and $H = H(p_1, \cdots, p_N, q_1, \cdots, q_N)$ is the Hamiltonian, viz., the total energy of the system. Let
$\vec{p} =\trsp{[p_1, \cdots, p_N]}\in \ES{R}{N}{1}$,
$\vec{q} = \trsp{[q_1, \cdots, q_N]}\in \ES{R}{N}{1}$, and
$\vec{z} =\trsp{[p_1, \cdots, p_N, q_1, \cdots, q_N]}= \trsp{[\trsp{\vec{p}}, \trsp{\vec{q}}]}\in \ES{R}{2N}{1}$,
then $H = H(\vec{p}, \vec{q}) = H(\vec{z})$ can be specified by  $\vec{z}$ in the $2N$-dimensional phase space. Since the gradient of $H$ is

=\trsp{\left[ \fracpde{H}{p_1}, \cdots, \fracpde{H}{p_N}, \fracpde{H}{q_1}, \cdots, \fracpde{H}{q_N}\right]} \in \ES{R}{2N}{1}.

Then the canonical equation is equivalent to

\begin{align} \label{eq-Heq}
\fracode{\vec{z}}{t} =  \inv{\mat{J}_{2N}} \cdot H_{\vec{z}}(\vec{z}), \quad
\mat{J}_{2N} = \begin{bmatrix} \mat{O}_N & \mat{I}_N \\ -\mat{I}_N & \mat{O}_N\end{bmatrix}
\end{align}

where $\mat{I}_N$ is the $N$-by-$N$ identical matrix, $\mat{O}_N$ is the $N$-by-$N$ zero matrix and $\mat{J}_{2N}$
is the $2N$-th order standard symplectic matrix [37, 30]. For simplicity, the subscripts in $\mat{I}_{N}, \mat{I}_{2N}$ and $\mat{J}_{2N}$ may be omitted. Any system which can be described by  (<ref>) is called an Hamiltonian system. There are some fundamental properties for the canonical equation of Hamiltonian system [35, 38, 39, 40]:

(i) it is invariant under the symplectic transform (phase flow);

(ii) the evolution of the system is the evolution of symplectic transform;

(iii) the symplectic symmetry and the total energy of the system can be preserved simultaneously and automatically.

§ COMPUTATION OF POLYNOMIALS

For the polynomial
$p_s(x) = \sum\limits^s_{i=0} c_i x^i = c_0 + c_1 x + \cdots + c_s x^s
we have the <ref> for computing this polynomial with linear computational complexity of time.

Computing Polynomials by Horner's Rule

[1]
The degree $s$, coefficients $c_0, c_1, \cdots, c_s$ and variable $x$ of the polynomial to be calculated.
$p_s(x) = \sum\limits^s_{i=0} c_i x^i$
Polynomial$c_0, \cdots, c_s, x$
$T \gets c_s$;
($i\gets1; i \le  s; i\gets i + 1$)
$T \gets T \cdot x + c_{s-i}$;
$p_s(x) \gets T$;

[1]
Bong Wie and Peter M Barba.
Quaternion feedback for spacecraft large angle maneuvers.
Journal of Guidance, Control, and Dynamics, 8(3):360–365, May

[2]
Bong Wie, H Weiss, and A Arapostathis.
Quarternion feedback regulator for spacecraft eigenaxis rotations.
Journal of Guidance, Control, and Dynamics, 12(3):375–380, May

[3]
Bong Wie and Jianbo Lu.
Feedback control logic for spacecraft eigenaxis rotations under slew
rate and control constraints.
Journal of Guidance, Control, and Dynamics, 18(6):1372–1379,
Nov. 1995.

[4]
Bong Wie, David Bailey, and Christopher Heiberg.
Rapid multitarget acquisition and pointing control of agile
Journal of Guidance, Control, and Dynamics, 25(1):96–104, Jan.

[5]
J Kuipers.
Quaternions and Rotation Sequences: A Primer with Applications
to Orbits, Aerospace and Virtual Reality.
Princeton University Press, 2002.

[6]
Bernard Friedland.
Analysis strapdown navigation using quaternions.
IEEE Transactions on Aerospace and Electronic Systems,
AES-14(5):764–768, Sept. 1978.

[7]
Anthony Kim and MF Golnaraghi.
A quaternion-based orientation estimation algorithm using an inertial
measurement unit.
In Position Location and Navigation Symposium, 2004. PLANS
2004, pages 268–272. IEEE, Apr. 2004.

[8]
R.M. Rogers.
Applied Mathematics in Integrated Navigation Systems.
AIAA, 2007.

[9]
Yong Min Zhong, She Sheng Gao, and Wei Li.
A quaternion-based method for SINS/SAR integrated navigation
IEEE Transactions on Aerospace and Electronic Systems,
48(1):514–524, Jan. 2012.

[10]
Joseph M Cooke, Michael J Zyda, David R Pratt, and Robert B McGhee.
NPSNET: Flight simulation dynamic modeling using quaternions.
Presence: Teleoperators & Virtual Environments, 1(4):404–420,

[11]
David Allerton.
Principles of flight simulation.
John Wiley & Sons, 2009.

[12]
Dominic J Diston.
Computational Modelling and Simulation of Aircraft and the
Environment, Volume 1: Platform Kinematics and Synthetic Environment.
John Wiley & Sons, 2009.

[13]
Jon S. Berndt, Tony Peden, and David Megginson.
JSBSim, open source flight dynamics software library.

[14]
J. S. Yuan.
Closed-loop manipulator control using quaternion feedback.
IEEE Journal on Robotics and Automation, 4(4):434–440, Aug.

[15]
Janez Funda, Russell H Taylor, and Richard P Paul.
On homogeneous transforms, quaternions, and computational efficiency.
IEEE Transactions on Robotics and Automation, 6(3):382–388,
Jun. 1990.

[16]
Jack CK Chou.
Quaternion kinematic and dynamic differential equations.
IEEE Transactions on Robotics and Automation, 8(1):53–64, Feb.

[17]
O-E Fjellstad and Thor I Fossen.
Position and attitude tracking of AUV's: A quaternion feedback
IEEE Journal of Oceanic Engineering, 19(4):512–518, Oct. 1994.

[18]
Mohammad Zamani, Jochen Trumpf, and Robert Mahony.
Minimum-energy filtering for attitude estimation.
IEEE Transactions on Automatic Control, 58(11):2917–2921, Nov.

[19]
Raul Murartal, J. M. M. Montiel, and Juan D. Tardos.
ORB-SLAM: A Versatile and Accurate Monocular SLAM System.
IEEE Transactions on Robotics, 31(5):1147–1163, 2015.
arXiv:1502.00956v2 [cs.RO].

[20]
CH Robotics content library: Understandiing quaternions.

[21]
Richard Szeliski.
Computer Vision: Algorithms and Applications.
Texts in Computer Science. Springer, London, 2010.
Available on line: <http://szeliski.org/Book/>.

[22]
John Vince.
Quaternions for Computer Graphics.
Springer, 2011.

[23]
Kit Ian Kou and Yong-Hui Xia.
Linear quaternion differential equations: Basic theory and
fundamental results.
arXiv:1510.02224v5 [math.CA]. available at:
<https://arxiv.org/abs/1510.02224>, September 7, 2017.

[24]
Alfred C Robinson.
On the use of quaternions in simulation of rigid-body motion.
Technical report, Wright Air Development Center, Dec. 1958.

[25]
ShunJin Wang and Hua Zhang.
Algebraic dynamics algorithm: Numerical comparison with Runge-Kutta
algorithm and symplectic geometric algorithm.
Science in China Series G: Physics Mechanics and Astronomy,
50(1):53–69, Feb. 2007.

[26]
John C Butcher.
Implicit Runge-Kutta processes.
Mathematics of Computation, 18(85):50–64, Jan, 1964.

[27]
Arieh Iserles.
A First Course in the Numerical Analysis of Differential
Cambridge University Press, 1996.

[28]
James M Varah.
On the efficient implementation of implicit Runge-Kutta methods.
Mathematics of Computation, 33(146):557–561, 1979.

[29]
Hong-Yan Zhang, Zi-Hao Wang, Lu-Sha Zhou, Qian-Nan Xue, Long Ma, and Yi-Fan
Explicit symplectic geometric algorithm for quaternion kinematical
differential equation.
IEEE/CAA Journal of Automatica Sinica, 5(2):479–488, Mar.

[30]
Kang Feng and Mengzhao Qin.
Symplectic Geometric Algorithms for Hamiltonian systems.
Springer, 2010.

[31]
Walter Rudin.
Real Analysis and Complex Analysis.
McGraw-Hill, 3rd edition, 1987.

[32]
Jiu-Shao Qin.
Available from
<http://www-history.mcs.st-andrews.ac.uk/Printonly/Qin_Jiushao.html> and

[33]
Allan Borodin.
Horner's Rule is Uniquely Optimal, pages 45–58.
Academic Press, 12 1971.
In book: Theory of Machines and Computations.

[34]
Donald E. Knuth.
The Art of Computer Programming: Seminumerical Algorithms, pp.
467-469, 1998., volume 2.
Addison-Wesley, 3rd edition, 1998.

[35]
Vladimir I Arnold.
Mathematical Methods of Classical Mechanics, volume 60 of Graduate Texts in Mathematics.
Springer, Apr. 1989.

[36]
Vladimir I Arnold, Valery V Kozlov, and Anatoly I Neishtadt.
Mathematical Aspects of Classical and Celestial Mechanics,
volume 3.
Springer, 2007.

[37]
Kang Feng.
On difference schemes and symplectic geometry.
In Proceeding of the 1984 Beijing Symposium on differential
geometry and differential equations., pages 42–58. Science Press, 1984.

[38]
Ling Hua Kong, Jia Lin Hong, Lan Wang, and Fang Fang Fu.
Symplectic integrator for nonlinear high order Schrödinger
equation with a trapped term.
Journal of computational and applied mathematics,
231(2):664–679, Sept. 2009.

[39]
JM Franco and I Gómez.
Construction of explicit symmetric and symplectic methods of
Runge–Kutta–Nyström type for solving perturbed oscillators.
Applied Mathematics and Computation, 219(9):4637–4649, Jan.

[40]
David M Hernandez.
Fast and reliable symplectic integration for planetary system
N-body problems.
Monthly Notices of the Royal Astronomical Society,
458(4):4285–4296, Jun. 2016.