A minimal-variable symplectic method for isospectral flows

TL;DR

This paper introduces a simple, second-order numerical integrator for isospectral flows that is intrinsically defined on quadratic Lie algebras and symmetric matrices, preserving spectral properties and Lie--Poisson structure.

Contribution

It presents a minimal-variable, structure-preserving integrator that is computationally efficient and applicable to a broad class of isospectral flows, improving upon existing methods.

Findings

The integrator is isospectral for general flows.

It preserves Lie--Poisson structure when Hamiltonian.

It is simple, second-order, and structure-preserving.

Abstract

Isospectral flows are abundant in mathematical physics; the rigid body, the the Toda lattice, the Brockett flow, the Heisenberg spin chain, and point vortex dynamics, to mention but a few. Their connection on the one hand with integrable systems and, on the other, with Lie--Poisson systems motivates the research for optimal numerical schemes to solve them. Several works about numerical methods to integrate isospectral flows have produced a large varieties of solutions to this problem. However, many of these algorithms are not intrinsically defined in the space where the equations take place and/or rely on computationally heavy transformations. In the literature, only few examples of numerical methods avoiding these issues are known, for instance, the \textit{spherical midpoint method} on $\SO(3)$. In this paper we introduce a new minimal-variable, second order, numerical integrator for isospectral flows intrinsically defined on quadratic Lie algebras and symmetric matrices. The algorithm is isospectral for general isospectral flows and Lie--Poisson preserving when the isospectral flow is Hamiltonian. The simplicity of the scheme, together with its structure-preserving properties, makes it a competitive alternative to those already present in literature.

Equations

$$\begin{array}[]{ll}&\dot{W}=[B(W),W],\quad W\in V\subset\mathfrak{gl}(n,\mathbb{C})\\ &W(0)=W_{0}.\end{array}$$

$$\begin{array}[]{ll}&\dot{W}=[\nabla H(W)^{\dagger},W],\quad W\in\mathfrak{g}\subset\mathfrak{gl}(n,\mathbb{C})\\ &W(0)=W_{0},\end{array}$$

$$N(V)=\{g\in G\mid g^{-1}Vg\subseteq V\}$$

$$\mathfrak{n}(V)=\{\xi \in \mathfrak{g}\mid [\xi ,V]\subseteq V\}$$

$$\mathfrak{c}(V)=\{\xi \in \mathfrak{g}\mid [\xi ,V]=0\}$$

$$W\in \mathfrak{g}⟺W^{†}J+JW=0.$$

$$W^{†}J+JW=0,$$

$$WJ+J{W}^{†}=0.$$

$$[\mathfrak{g},\mathfrak{g}]\subset \mathfrak{g},[\mathfrak{g},{\mathfrak{g}}^{\perp }]\subset {\mathfrak{g}}^{\perp },[{\mathfrak{g}}^{\perp },{\mathfrak{g}}^{\perp }]\subset \mathfrak{g}.$$

$$N\in {\mathfrak{g}}_0⟺N^{†}P+PN=0.$$

$$\boxed{\begin{array}[]{ll}&W_{k}=(Id-\frac{h}{2}B(\widetilde{W}))\widetilde{W}(Id+\frac{h}{2}B(\widetilde{W}))\\ &W_{k+1}=(Id+\frac{h}{2}B(\widetilde{W}))\widetilde{W}(Id-\frac{h}{2}B(\widetilde{W})),\end{array}}$$

$$Y=(Id-\frac{h}{2}B(X))X(Id+\frac{h}{2}B(X))$$

$$\frac{d}{dh}X$$

$$X(0)$$

$$W_{k+1}=\text{mbox}Cay{\left(\frac{h}{2}B(\tilde{W})\right)}^{-1}{W}_k\text{mbox}Cay\left(\frac{h}{2}B(\tilde{W})\right).$$

$$\begin{array}[]{llll}&X=-h(W_{k}+\frac{1}{2}X)B(\widetilde{W})\\ &Y=hB(\widetilde{W})(W_{k}+\frac{1}{2}Y)\\ &K=\frac{h}{2}B(\widetilde{W})(X+K)\\ &\widetilde{W}=W_{k}+\frac{1}{2}(X+Y+K)\\ &W_{k+1}=W_{k}+h[B(\widetilde{W}),\widetilde{W}],\end{array}$$

$$Y+K=hB(\tilde{W})\tilde{W}$$

$$X=-h\tilde{W}B(\tilde{W})+\frac{h^2}{2}B(\tilde{W})\tilde{W}B(\tilde{W}).$$

$$W_{k+1}=W_k+h\sum _{i=1}^s{b}_i[B({\tilde{W}}_i),{\tilde{W}}_i],$$

$$H(W)=\frac{1}{2}\mathrm{Tr}(({\mathcal{I}}^{-1}W{)}^{†}W),W\in \mathfrak{so}(n),$$

$$\left.\begin{array}[]{ll}\dot{W}=-[\mathcal{I}^{-1}W,W]\\ W(0)=W_{0}.\end{array}\right.$$

$$(\mathcal{I}^{-1}W)_{ij}=\left\{\begin{array}[]{ll}\frac{W_{ij}}{i},\quad i=1,\ldots,5,j=1,\ldots,10\\ \frac{W_{ij}}{11-i},\quad i=6,\ldots,10,j=1,\ldots,10\end{array}\right.$$

$$(W_0{)}_{ij}=1/10\text{for}i<j\text{and}{W}_0^{†}=-W_0$$

$$\dot{W}=[[N,W],W],$$

$$\begin{array}[]{ll}&w_{k}=\widetilde{w}+\frac{h}{2}\widetilde{w}\times B(\widetilde{w})-\frac{h^{2}}{4}B(\widetilde{w})(B(\widetilde{w})\cdot\widetilde{w})\\ &w_{k+1}=\widetilde{w}-\frac{h}{2}\widetilde{w}\times B(\widetilde{w})-\frac{h^{2}}{4}B(\widetilde{w})(B(\widetilde{w})\cdot\widetilde{w}),\end{array}$$

$$\begin{array}[]{ll}&w_{k+1}=w_{k}+h\dfrac{\sqrt{w_{k+1}}\sqrt{w_{k}}(w_{k+1}+w_{k})}{|w_{k+1}+w_{k}|}\times B\left(\dfrac{\sqrt{w_{k+1}}\sqrt{w_{k}}(w_{k+1}+w_{k})}{|w_{k+1}+w_{k}|}\right).\end{array}$$

$$\begin{array}[]{ll}&w_{k+1}=w_{k}+h\dfrac{w_{k+1}+w_{k}}{2}\times B\left(\dfrac{w_{k+1}+w_{k}}{2}\right).\end{array}$$

$${\dot{w}}_i=w_i\times (w_{i-1}+w_{i+1}),$$

$$\dot{w}=w\times {\partial }_{xx}w,$$

$$H(w_1,\dots ,w_N)=\sum _{i=1}^N{w}_i\cdot w_{i+1}.$$

$$w(x)=(\cos (2\pi x^2)\sin (2\pi x^3),\sin (2\pi x^2)\sin (2\pi x^3),\cos (2\pi x^3)),$$

Full text

1. Introduction

The numerical integration of isospectral flows is a classical subject of study in numerical analysis [6, 9]. The interest in this problem is motivated by the numerical simulation of integrable systems, which are deeply related to isospectral flows via the Lax pair formulation. The quasi-periodic dynamics of integrable systems depends on the presence of a large number of first integrals. In the Lax pair formulation, some of these first integrals can be presented as a linear combination of the eigenvalues of the dynamical variable. Therefore, the preservation of the spectrum of the dynamical variable is a key feature of a numerical scheme applied to isospectral flows, in order to expect the right qualitative behaviour of the discrete approximate solutions [6]. Furthermore, as a special case, Lie–Poisson systems on the dual of reductive Lie algebras can be seen as isospectral flows [18]. A reductive Lie algebra is defined as the direct sum of a semisimple Lie algebra and an abelian Lie algebra. In this paper, the crucial property of real semisimple Lie algebras is that they can be represented as matrix Lie algbras which are closed under conjugate transpose [11, Prop. 6.28]. Moreover, any real matrix Lie algebra which is closed under conjugate transpose is reductive [11, Prop. 1.56].111In view of [11, Sec. I.8], there is no restriction in looking at real matrix Lie algebras, since the complex ones can be seen as real matrix Lie algebras of double dimension. Because of this, Lie–Poisson systems on the dual of a reductive Lie algebra can be equivalently seen as isospectral flows of the form (2) below. Lie–Poisson systems originate from the Poisson reduction of canonical Hamiltonian systems on the cotangent bundle of a Lie group [14]. Classical examples of Lie–Poisson systems are the rigid body [21], the heavy top and the incompressible Euler equations [1]. It is known from the backward error analysis, that Lie–Poisson preserving numerical schemes are superior to standard methods when applied to Lie–Poisson systems, especially for long-time simulations. As state-of-the-art, well established theories on numerical methods for both isospectral and Lie–Poisson systems exist in the literature (see for example [4, 9, 20]). For Lie–Poisson systems various symplectic algorithms have been developed (see [15] for a recent survey). However, few examples of numerical schemes that are intrinsically defined in the space where the dynamics takes place are known (e.g. [16]). This issue often causes a lack of efficiency for these schemes, which rely on group to algebra maps (e.g. the matrix exponential or the Cayley map) or a large number of unknowns. Before presenting our result, let us introduce the mathematical setup used throughout the paper.

Isospectral flows are first order ODEs of the form:

[TABLE]

Here, $[\cdot,\cdot]$ denotes the matrix commutator, $V$ is a linear subspace of the Lie algebra $\mathfrak{gl}(n,\mathbb{C})$, and the function $B\colon V\to\mathfrak{n}(V)$ maps $V$ into its $\mathfrak{gl}(n,\mathbb{C})$-normalizer algebra $\mathfrak{n}(V)$222$\mathfrak{n}(V)$ is the largest Lie subalgebra of $\mathfrak{gl}(n,\mathbb{C})$ such that $[\mathfrak{n}(V),V]\subset V$ (see Definition 1 in Section 2).. The most studied case in literature is when $V=\operatorname{Sym}(n,\mathbb{R})$ is the space of symmetric real matrices, for which the normalizer is the Lie algebra of skew-symmetric real matrices $\mathfrak{n}(V)=\mathfrak{so}(n)$. For Lie–Poisson systems on the dual of a reductive Lie algebra, we have that $V=\mathfrak{g}^{*}$ is the dual of a reductive Lie subalgebra of $\mathfrak{gl}(n,\mathbb{C})$, for which the $\mathfrak{gl}(n,\mathbb{C})-$normalizer is $\mathfrak{n}(V)=\mathfrak{g}_{0}\oplus\mathfrak{c}(\mathfrak{g})$ (see Definition 2 and Lemma 2, in Section 2). Throughout the paper, we identify $\mathfrak{g}^{*}$ with $\mathfrak{g}$, via the Frobenius inner product $\langle A,B\rangle=\operatorname{Tr}(A^{\dagger}B)$, where † is the conjugate transpose. Via this identifications, Lie–Poisson systems on the dual of a reductive Lie algebra $\mathfrak{g}$ take the form:

[TABLE]

for $H:\mathfrak{g}\rightarrow\mathbb{C}$, a smooth function called Hamiltonian.

A class of numerical methods to solve (1)-(2), called Isospectral Symplectic Runge-Kutta (IsoSyRK), has been introduced in [18]. In the case of Lie–Poisson systems these schemes are symplectic. In this paper, we focus on the IsoSyRK associated to the implicit midpoint method, which turns out to have a specially nice structure. On the one hand, we provide a simpler proof (avoiding the use of the B-series theory) that for the implicit midpoint method, the respective IsoSyRK defined in [18] is isospectral for any $B=B(W)$. On the other hand, we derive a simpler scheme reducing the number of unknowns up to minimality, revealing an intrinsic relation between the implicit midpoint method and the Cayley transform. The resulting integrator, although implicit, is second order, isospectral and symplectic when the isospectral flow is Lie–Poisson. The scheme is also intrinsically defined on $V$, for a large class of isospectral flows (see section 2 for details) and, when the isospectral flow is Lie–Poisson, it preserves the coadjoint orbits. Furthermore, only one evaluation of $B(\cdot)$ and two matrix multiplications per iteration are required, making the scheme very efficient. In the last section of this paper, we show some numerical examples of our scheme and we compare it with the spherical midpoint method, which is another minimal-variable Lie–Poisson integrator on $\mathbb{R}^{3}$. Finally, we show how our scheme looks on $\mathfrak{sl}(2,\mathbb{R})$, defining what we call the hyperbolic midpoint method.

Acknowledgements. The author was supported by EU Horizon 2020 grant No 691070, by the Swedish Foundation for International Cooperation in Research and Higher Eduction (STINT) grant No PT2014-5823, by the Swedish Foundation for Strategic Research grant ICA12-0052, and by the Swedish Research Council (VR) grant No 2017-05040. The author would like to thank Klas Modin for the support and the enlightening discussions during the work on this paper.

2. Main result

Let us consider an isospectral flow of the form (1). In order to present our result, we need a short detour on some concepts and basic results on Lie algebras. As already mentioned in section 1, for (1) to be well defined, we have to require $B(\cdot)$ to take values in the $\mathfrak{gl}(n,\mathbb{C})$-normalizer algebra of $V$. We recall here the definition of the normalizer Lie group and normalizer Lie algebra:

Definition 1**.**

Let $G$ be a Lie group and $\mathfrak{g}$ its Lie algebra. Furthermore, let $V\subseteq\mathfrak{g}$ be a linear subspace. Then the two sets

[TABLE]

are respectively called the $G$-normalizer and the $\mathfrak{g}$-normalizer of $V$. Notice that $N(V)$ is a subgroup of $G$ and $\mathfrak{n}(V)$ is a Lie subalgebra of $\mathfrak{g}$.

A related concept to normalizer is the centralizer Lie algebra.

Definition 2**.**

Let $\mathfrak{g}$ be a Lie algebra and let $V\subseteq\mathfrak{g}$ be a linear subspace. Then the set

[TABLE]

is called the $\mathfrak{g}$-centralizer of $V$. Notice that $\mathfrak{c}(V)$ is a Lie subalgebra of $\mathfrak{g}$.

We now recall the definition of a $J-$quadratic Lie algebra.

Definition 3**.**

A Lie subalgebra $\mathfrak{g}$ of $\mathfrak{gl}(n,\mathbb{C})$ is called $J-$quadratic Lie algebra if there exists an invertible matrix $J$, such that

[TABLE]

Lemma 1**.**

Let $\mathfrak{g}$ be a Lie subalgebra of $\mathfrak{gl}(n,\mathbb{C})$ such that there exists a matrix $J$ for which

[TABLE]

Then $\mathfrak{g}=\mathfrak{g}^{\dagger}$ implies $J^{2}\in\mathfrak{c}(\mathfrak{g})$. Moreover, if $J$ is invertible and $J^{2}\in\mathfrak{c}(\mathfrak{g})$, then $\mathfrak{g}=\mathfrak{g}^{\dagger}$.

Proof.

Suppose $\mathfrak{g}=\mathfrak{g}^{\dagger}$. Then, for all $W\in\mathfrak{g}$, both the following identities hold:

[TABLE]

The second of these implies that $WJ^{2}+JW^{\dagger}J=0$ and the first one $JW^{\dagger}J+J^{2}W=0$. Subtracting these identities, we get $[W,J^{2}]=0$. Hence $J^{2}\in\mathfrak{c}(\mathfrak{g})$.

Now assume that $\mathfrak{g}$ is $J-$quadratic and $J^{2}\in\mathfrak{c}(\mathfrak{g})$. Then, for all $W\in\mathfrak{g}$ we have $0=JW^{\dagger}J+J^{2}W=JW^{\dagger}J+WJ^{2}$ and hence $WJ+JW^{\dagger}=0$, being $J$ invertible. Therefore $\mathfrak{g}$ and $\mathfrak{g}^{\dagger}$ are defined by the same identity and they coincide. ∎

Lemma 2**.**

Let $\mathfrak{g}$ be a Lie subalgebra of $\mathfrak{gl}(n,\mathbb{C})$ such that $\mathfrak{g}=\mathfrak{g}^{\dagger}$. Then the $\mathfrak{gl}(n,\mathbb{C})-$normalizer of $\mathfrak{g}$ is $\mathfrak{n}(\mathfrak{g})=\mathfrak{g}_{0}\oplus\mathfrak{c}(\mathfrak{g})$, where $\mathfrak{g}_{0}$ is the semisimple ideal of $\mathfrak{g}$ such that $\mathfrak{g}=\mathfrak{g}_{0}\oplus\mathfrak{z}(\mathfrak{g})$, for $\mathfrak{z}(\mathfrak{g})$ the center of $\mathfrak{g}$.333Such a decomposition always exists being $\mathfrak{g}$ reductive by [11, Prop. 1.56].

Proof.

Let $\mathfrak{g}^{\perp}$ be the orthogonal complement of $\mathfrak{g}$ in $\mathfrak{gl}(n,\mathbb{C})$ with respect to the Frobenius inner product. It is not hard to check that the following properties hold:

[TABLE]

Hence, if $A\in\mathfrak{n}(\mathfrak{g})\cap\mathfrak{g}^{\perp}$ it must be $[\mathfrak{g},A]=0$. Therefore, $A$ has to be in $\mathfrak{c}(\mathfrak{g})$. Moreover, we have that the following inclusions always hold $\mathfrak{z}(\mathfrak{g})\subset\mathfrak{c}(\mathfrak{g})\subset\mathfrak{n}(\mathfrak{g})$. Therefore, $\mathfrak{n}(\mathfrak{g})=\mathfrak{g}_{0}\oplus\mathfrak{c}(\mathfrak{g})$, being $\mathfrak{g}_{0}$ centerless. ∎

Notice that, since always $\mathfrak{n}(\mathfrak{g})^{\dagger}=\mathfrak{n}(\mathfrak{g}^{\perp})$, we have that $\mathfrak{n}(\mathfrak{g}^{\perp})=\mathfrak{g}_{0}\oplus\mathfrak{c}(\mathfrak{g})^{\dagger}$, whenever $\mathfrak{g}=\mathfrak{g}^{\dagger}$. In conclusion, we have the following corollary:

Corollary 1**.**

Let $\mathfrak{g}$ be a $J-$quadratic Lie subalgebra of $\mathfrak{gl}(n,\mathbb{C})$ such that $J^{2}\in\mathfrak{c}(\mathfrak{g})$. Then the $\mathfrak{gl}(n,\mathbb{C})-$normalizer of $\mathfrak{g}$ and $\mathfrak{g}^{\perp}$ are respectively $\mathfrak{n}(\mathfrak{g})=\mathfrak{g}_{0}\oplus\mathfrak{c}(\mathfrak{g})$ and $\mathfrak{n}(\mathfrak{g}^{\perp})=\mathfrak{g}_{0}\oplus\mathfrak{c}(\mathfrak{g})^{\dagger}$. In particular, under the identification of $\mathfrak{g}^{*}$ with $\mathfrak{g}$ via the Frobenius inner product, any Lie–Poisson system on $\mathfrak{g}^{*}$ can be written in the form (2).

Notice that by [11, Prop. 1.56] any $\mathfrak{g}$ such as in Corollary 1 is a reductive Lie algebra. As mention in the introduction, Lie–Poisson systems on the dual of reductive Lie algebras can be written in the form (2). In particular, this is true for Lie–Poisson systems on the dual of $\mathfrak{g}\oplus\mathfrak{Z}$, where $\mathfrak{Z}$ is an Abelian Lie algebra. We can now state the main result of this paper.

Theorem 1**.**

Let $W_{k}\in D\subset V$, for a domain $D$ in the linear subspace $V\subset\mathfrak{gl}(n,\mathbb{C})$. Assume that the normalizer splits as $\mathfrak{n}(V)=\mathfrak{g}_{0}\oplus\mathfrak{c}(V)$, for some Lie algebra $\mathfrak{g}_{0}$, which satisfies

[TABLE]

for some constant matrix $P$. Furthermore, let $B:D\subset\mathfrak{gl}(n,\mathbb{C})\rightarrow\mathfrak{n}(V)$ be continuously differentiable. Then, for some $h>0$, there exists $\widetilde{W}\in V$ such that the numerical scheme $W_{k}\rightarrow W_{k+1}$, implicitly defined by:

[TABLE]

is a second order isospectral integrator for (1), for any $k\geq 0$.444Here $Id$ denotes the $n\times n$ identity matrix.. Moreover, when (1) is a Lie-Poisson system on $\mathfrak{gl}(n,\mathbb{C})^{*}$ or on the dual of some $J-$quadratic Lie algebra $\mathfrak{g}$ such that $J^{2}\in\mathfrak{c}(\mathfrak{g})$ (or even on $\mathfrak{g}\oplus\mathfrak{Z}$, where $\mathfrak{Z}$ is an Abelian Lie algebra), then (5) is a Lie-Poisson integrator for (1) which preserves the coadjoint orbits in $\mathfrak{g}^{*}$.

Remark 1**.**

The main contribution of Theorem 1, with respect to the results presented in [18], is that the scheme (5) is a minimal-variable isospectral (Lie–Poisson) integrator. Minimal-variable means here that the only unknown is $\widetilde{W}$, which lives in a vector space of dimension $\dim(V)$. Hereafter, we will refer to the scheme (5) as the isospectral minimal midpoint. Moreover, the proof of the properties of (5), unlikely to [18, Cor. 1], does not require any application of the B-series theory and reveals a deep connection with the Cayley transform (see the proof of Lemma 4). The latter is a quite interesting fact because the Cayley transform arises as a necessary consequence of the use of the implicit midpoint scheme and not, as it has always appeared in literature, as a prescribed choice to construct a certain numerical scheme. We also emphasize that the condition (4) and the ones on $J$ in Theorem 1 to get a Lie–Poisson integrator is slightly more general than the one considered in [18, Thm. 1-2].

We will give the proof of Theorem 1 in some lemmas.

Lemma 3**.**

Let $B:D\subset\mathfrak{gl}(n,\mathbb{C})\rightarrow\mathfrak{gl}(n,\mathbb{C})$ continuously differentiable in the domain $D$. Then, for every $Y\in D$, there exist $\overline{h}>0$ such that the equation

[TABLE]

has a solution $X\in\mathfrak{gl}(n,\mathbb{C})$ for any $0\leq h<\overline{h}$.

Proof.

In order to get a solution to (6), we consider the function $F_{h}(X):=Y+\frac{h}{2}[B(X),X]+\frac{h^{2}}{4}B(X)XB(X)$, such that (6) is equal to $X=F_{h}(X)$. In order to determine $\overline{h}$, we consider the initial value problem:

[TABLE]

Since $\dfrac{\partial F_{h}(X)}{\partial h}$ is continuous, if we prove that the operator $(Id-DF_{h}(X))$ is continuous and invertible, then the Peano existence theorem will ensure a solution $X(h)$, for any $h$ in some interval $[0,\overline{h})$, for $\overline{h}>0$. Indeed, $DF_{h}(X)=hG(DB(X),B(X),X,h)$, where $G$ is polynomial in its variables and $DB(X)$ is continuous by hypothesis. Hence, $DF_{h}(X)\rightarrow 0$, for $h\rightarrow 0$, therefore there exist some $\overline{h}>0$ such that $(Id-DF_{h}(X))$ is invertible, for any $0\leq h<\overline{h}$. ∎

Lemma 4**.**

Let $B:D\subset\mathfrak{gl}(n,\mathbb{C})\rightarrow\mathfrak{gl}(n,\mathbb{C})$ continuously differentiable in the domain $D$ and let $0\leq h<\overline{h}$ as in lemma 3. Then, for every $W_{k}\in D$, the numerical scheme $W_{k}\rightarrow W_{k+1}$ is isospectral. Moreover, if $D\subset V$, for $D$ domain in the linear subspace $V\subset\mathfrak{gl}(n,\mathbb{C})$, and $B:D\subset\mathfrak{gl}(n,\mathbb{C})\rightarrow\mathfrak{n}(V)$, where $\mathfrak{n}(V)=\mathfrak{g}_{0}\oplus\mathfrak{c}(V)$, for some Lie algebra $\mathfrak{g}_{0}$ which satisfies (4), then $W_{k+1}\in V$. Furthermore, when $V$ is a $J-$quadratic Lie algebra such that $J^{2}\in\mathfrak{c}(\mathfrak{g})$, then $W_{k+1}\in\mathcal{O}_{W_{k}}\subset V$, where $\mathcal{O}_{W_{k}}$ is the coadjoint orbit of which $W_{k}$ belongs.

Proof.

Clearly $W_{k+1}=(Id+\frac{h}{2}B(\widetilde{W}))(Id-\frac{h}{2}B(\widetilde{W}))^{-1}W_{k}(Id-\frac{h}{2}B(\widetilde{W}))(Id+\frac{h}{2}B(\widetilde{W}))^{-1}$ and hence $W_{k+1}$ and $W_{k}$ are similar. Furthermore, we notice that $(Id-\frac{h}{2}B(\widetilde{W}))(Id+\frac{h}{2}B(\widetilde{W}))^{-1}=\mbox{Cay}(\frac{h}{2}B(\widetilde{W}))$, where Cay is the Cayley transform. Therefore we have that:

[TABLE]

Assuming $\mathfrak{n}(V)=\mathfrak{g}_{0}\oplus\mathfrak{c}(V)$, for some Lie algebra $\mathfrak{g}_{0}$ which satisfies (4), by [6, Lemma IV.8.7], $\mbox{Cay}\left(\frac{h}{2}B(\widetilde{W})\right)$ is in the normalizer group $N(V)$ of $V$ and therefore $W_{k+1}$ is in $V$ as well. When $V=\mathfrak{g}^{*}$ for $\mathfrak{g}$ a $J-$quadratic Lie algebra such that $J^{2}\in\mathfrak{c}(\mathfrak{g})$, the transformation (7) coincides with the coadjoint action of $G$ on $\mathfrak{g}^{*}$, where $G$ is the respective connected component to the identity of a Lie group with Lie algebra $\mathfrak{g}$. Therefore, (7) fixes the coadjoint orbits. ∎

Remark 2**.**

We point out that the equation (7) reveals an interesting relation between the Cayley transform and the implicit midpoint method. Indeed, the fact that the Cayley transform appears as a consequence of the reduction of the implicit midpoint from $T^{*}GL(n,\mathbb{C})$ to $\mathfrak{gl}(n,\mathbb{C})^{*}$ may indicate a deeper, perhaps canonical, relation between symplectic Runge–Kutta methods and the Cayley transform. The former are associated to conservation of quadratic first integrals of ODEs and the latter to transforming quadratic Lie algebras in quadratic Lie groups.

Corollary 2**.**

Let the hypotesis of Lemma 3 hold. Then, if $\mathfrak{n}(V)=\mathfrak{g}_{0}\oplus\mathfrak{c}(V)$ and $\mathfrak{g}_{0}$ is a compact Lie algebra, there exists $\overline{h}>0$ independent from $k$ such that the scheme (5) has solution.

Proof.

Since $\mathfrak{g}_{0}$ is a compact Lie algebra, the associate connected Lie group $G$ is compact and therefore the orbits of the action (7) are compact. Hence, we can find a minimum $\overline{h}>0$ in Lemma 3 independent from the iteration $k\geq 0$. ∎

Lemma 5**.**

Let $B:D\subset\mathfrak{gl}(n,\mathbb{C})\rightarrow\mathfrak{gl}(n,\mathbb{C})$ continuously differentiable in the domain $D$ and let $0\leq h<\overline{h}$ as in lemma 3. Then, for every $W_{k}\in D$, the numerical scheme $W_{k}\rightarrow W_{k+1}$ in (5) descends from the method defined in [18, Def. 1] associated with the implicit midpoint method. In particular, if $B=\nabla H^{\dagger}$ for some function $H:D\subset\mathfrak{gl}(n,\mathbb{C})^{*}\rightarrow\mathbb{R}$, then the method is Lie–Poisson in $\mathfrak{gl}(n,\mathbb{C})^{*}$.

Proof.

Consider the second order method as defined in [18, Def. 1] associated with the implicit midpoint method:

[TABLE]

for $k\geq 0$ with unknowns $X,Y,K\in\mathfrak{gl}(n,\mathbb{C})$. It is not hard to check that the following identities hold:

[TABLE]

Applying these to (8) we get, after some computations, the scheme (5). In [18, Thm. 3], it has been proven that when $B(W)=\nabla H(W)^{\dagger}$, for some functions $H:\mathfrak{gl}(n,\mathbb{C})^{*}\rightarrow\mathbb{R}$, the method is a Lie-Poisson integrator in $\mathfrak{gl}(n,\mathbb{C})^{*}$. The scheme $W_{k}\rightarrow W_{k+1}$ defined in (5) coincides with (8), but with the elimination of the intermediate variables $X,Y,K$. Therefore, (5) is a Lie–Poisson integrator. ∎

Proof.

$[Theorem~{}\ref{thm1}]$. The proof simply follows from the lemmas. Lemma 3 says that the method (5) has solution for $h$ sufficiently small. Under the assumptions of Lemma 4 proves the isospectrality of the scheme and its intrinsic restriction $V$, when $\mathfrak{n}(V)=\mathfrak{g}_{0}\oplus\mathfrak{c}(V)$, for some Lie algebra $\mathfrak{g}_{0}$ which satisfies (4). Finally, in Lemma 5 is shown that the scheme descends from the isospectral midpoint method defined in [18] and therefore it is a second order Lie–Poisson integrator in $\mathfrak{gl}(n,\mathbb{C})^{*}$, when (1) is. Putting together Lemma 4, Lemma 5 and Corollary 1, we have that when a $V$ is the dual of a $J-$quadratic Lie algebra such that $J^{2}\in\mathfrak{c}(\mathfrak{g})$ (possibly plus a commutative Lie algebra) and (1) is Lie–Poisson, the scheme (5) is a Lie–Poisson integrator on $V$, which preserves the coadjoint orbits. ∎

Remark 3**.**

We notice that the isospectral minimal midpoint (5) is somehow similar to the modified implicit midpoint rule introduced in [4]. However, in their scheme, $\widetilde{W}$ was set to be $\frac{W_{k+1}+W_{k}}{2}$ which does not hold in general while solving the isospectral minimal midpoint (5). In fact, even though the scheme in [4] is isospectral, it is not symplectic.

Remark 4**.**

The isospectral minimal midpoint (5) can be derived in a different way, as proposed in [15]. The construction there is more general and (5) can be recovered choosing as a retraction map the Cayley transform instead of the exponential map. This surprising connection opens up a question about a geometrical description of the methods proposed in [18, Def. 1]. Let us consider for any $s=1,2,\dots$ and $a$ $s\times s$ real matrix, a retraction map $\tau_{a}:\mathfrak{g}^{\oplus s}\rightarrow G^{\times s}$. Then, similarly to [15], for each $i=1,\dots,s$, it is implicitly defined by the differential of the retraction map $d\tau_{a}$ a discrete map $W_{k}\mapsto\widetilde{W}_{i}\ni\mathfrak{g}$. Finally, we define our integrator $\Psi_{h}:W_{k}\mapsto W_{k+1}$ as:

[TABLE]

for some real numbers $b_{i}$. The question is whether, for any $s$-stage symplectic Runge–Kutta method, there exists a retraction map $\tau_{a}:\mathfrak{g}^{\oplus s}\rightarrow G^{\times s}$, such that any Lie–Poisson integrator defined in [18] can be obtained in this way.

3. Numerical examples

In this section we present some applications of the isospectral minimal midpoint (5) on isospectral flows and Lie–Poisson systems found in literature. We also compare our method in the case of $\mathfrak{so}(3)\cong(\mathbb{R}^{3},\times)$ with the spherical midpoint method, showing that the isospectral minimal midpoint (5) has the same computational cost. Finally, we show explicitly how the isospectral minimal midpoint (5) looks on $\mathfrak{sl}(2,\mathbb{R})$, applying it to the the point vortex equations on the hyperbolic plane. For the example considered in this section, we plot the variation of the first integrals of (1). As expected, we get exact conservation (up to round-off errors) of the Casimir functions and, when the flow is Hamiltonian, near conservation of the Hamiltonian.

3.1. The generalized rigid body

A classical example among Hamiltonian isospectral systems is the generalized rigid body. It represents a class of completely integrable systems on $\mathfrak{so}(n)$, for every $n\geq 1$, [13]. The Hamiltonian is given by

[TABLE]

where $\mathcal{I}\colon\mathfrak{so}(n)\rightarrow\mathfrak{so}(n)$ is a symmetric positive definite inertia tensor. The equations of motion are then

[TABLE]

We discretize this system for $n=10$. Our implementation uses Newton iterations for the non-linear system.555Fixed-point iteration could be used as well, no numerical issues have arisen in our experiments with it. The inertia tensor is given by

[TABLE]

and we use the stepsize $h=0.1$. The initial conditions are given by

[TABLE]

As shown in Figure 1, the Hamiltonian is nearly conserved and the Casimir functions are conserved up to the accuracy of the Newton iterations.

3.2. The Brockett flow

In this section we specify the isospectral minimal midpoint (5) for the Brockett flow, or double bracket flow:

[TABLE]

where $N,W$ are $n\times n$ self-adjoint complex matrices. In [3], Brockett shows that for $N$ diagonal matrix with distinct entries and $W_{0}$ self-adjoint matrix with distinct eigenvalues, for $t\rightarrow\infty$, $W(t)$ converges exponentially fast to a diagonal matrix with the eigenvalues sorted accordingly to the order of the entries of $N$. In figure 2, we plot the eigenvalue variation for a randomly generated666Throughout the whole paper, a randomly generated matrix is understood as a matrix whose components are defined as pseudorandom values drawn from the standard uniform distribution on the open interval $(0,1)$, generated via the MATLAB function rand. self-adjoint initial matrix $W_{0}$ of dimension $10\times 10$ and $N=diag(1,2,\dots,10)$. The asymptotic stationarity of the eigenvalues variation reflects the fact that $W$ is close to a diagonal matrix.

3.3. Lie–Poisson systems on $(\mathbb{R}^{3},\times)$

On $(\mathbb{R}^{3},\times)$ the isospectral minimal midpoint (5) can be written as:

[TABLE]

for $\widetilde{w},w_{k},w_{k+1}\in\mathbb{R}^{3}$ and $B:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$. We want to compare the isospectral minimal midpoint (14) with another minimal-variable symplectic integrator on $\mathbb{R}^{3}$ introduced in [17], i.e. the spherical midpoint method:

[TABLE]

Remark 5**.**

Let $B(\cdot)$ be orthogonal with respect to the rays, i.e. $B(w)\cdot w=0$, for every $w\in\mathbb{R}^{3}$. It is immediate to check that then (14) coincides with the classical midpoint scheme:

[TABLE]

In [17] it is shown that also (15) coincides with (16) when $B(\cdot)=\nabla H(\cdot)$, for some Hamiltonian function $H:\mathbb{R}^{3}\rightarrow\mathbb{R}$ constant on the rays (which implies $B(\cdot)$ to be orthogonal to the rays). In this case, (16) is known to be symplectic, whereas this fails for general Hamiltonian $H$. Therefore, (14) can be seen as the second order correction of (16) to be symplectic for any Hamiltonian $H$.

Let us now consider the two schemes (14) and (15). Both methods are implicit and therefore an implicit solver has to be used. Here we show that they exhibit the same computational cost. The example we consider is the Heisenberg spin chain on $\mathbb{R}^{3N}$. For this one has to extend both the isospectral minimal midpoint (14) and the spherical midpoint (15) to direct products of $\mathbb{R}^{3}$ (see [17] and [18]).

The Heisenberg spin chain of micromagnetics is defined as:

[TABLE]

where $w_{i}\in\mathbb{S}^{2}$, for $i=1,\ldots,N$ and $w_{N+1}=w_{1}$. It corresponds, up to scaling, to spatial discretization of the Landau-Lifschitz PDE:

[TABLE]

for $w:\mathbb{S}^{1}\rightarrow\mathbb{S}^{2}$ smooth. We notice that (17) is a Lie–Poisson system on $\mathbb{R}^{3N}$, with Hamiltonian:

[TABLE]

Clearly, to get a good approximation of (18), $N$ has to be large. In figure 3 we show the average time cost for time-step with respect to the number of spin particles for both the isospectral minimal midpoint (14) and the spherical midpoint (15). We conclude that the complexity grows similarly.

In terms of conservation properties, the two schemes exactly preserve the linear invariants and nearly conserve the the quadratic first integrals of the form $\sum_{i,j}w_{i}^{\dagger}Aw_{j}$. Moreover, the spherical midpoint has the advantage of exactly conserving all the quadratic first integrals of the form $w_{i}^{\dagger}Aw_{i}$, for some square matrix $A$, whereas the isospectral minimal midpoint (14) exactly conserves only the quadratic invariants $w_{i}^{\dagger}w_{i}$. In figure 4 we compare the isospectral minimal midpoint (14) and the spherical midpoint (15) for the initial data given as am equispaced discretization in $N=100$ points of the closed spherical curve:

[TABLE]

for $x\in[0,1]$.

We can conclude from figure 4 that the spherical midpoint performs slightly better than the isospectral minimal midpoint (14). However, both the schemes show the desired conservation properties due to their symplecticity.

Finally, in figure 5, we present the error diagram for the isospectral minimal midpoint (14) and the spherical midpoint (15). The curves in figure 5 show the expected second order of the schemes, with no significant difference.

3.4. Lie–Poisson systems on $\mathfrak{sl}(2,\mathbb{R})^{*}$

In this section we specify the isospectral minimal midpoint (5) in the case of the Lie algebra $\mathfrak{sl}(2,\mathbb{R})$, i.e. the $2\times 2$ matrices with zero trace. The first observation is that $\mathfrak{sl}(2,\mathbb{R})\cong\mathfrak{sp}(2,\mathbb{R})$, which is a non-compact $J-$quadratic Lie algebra with respect to $J=\begin{bmatrix}0&-1\\ 1&0\end{bmatrix}.$ From this, it is straightforward to see that, when $B:\mathfrak{gl}(2,\mathbb{R})\rightarrow\mathfrak{sl}(2,\mathbb{R})$ and $W_{k}\in\mathfrak{sl}(2,\mathbb{R})$, $\widetilde{W}$ in (5) is also in $\mathfrak{sl}(2,\mathbb{R})$. We notice that this is no longer true for $\mathfrak{sl}(n,\mathbb{R})$, for $n>2$. On the other hand, any element in $\mathfrak{sl}(2,\mathbb{R})$ can be written as a vector in $\mathbb{R}^{3}$, via the vector spaces isomorphism:

[TABLE]

for any $x,y,z\in\mathbb{R}$. In this coordinates, we can express the isospectral minimal midpoint (5) for $\mathfrak{sl}(2,\mathbb{R})$ as:

[TABLE]

for $\widetilde{w},w_{k},w_{k+1}\in\mathbb{R}^{3}$ and $B:\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$, where :

[TABLE]

and

[TABLE]

for any $a,b,c\in\mathbb{R}^{3}$.

We notice that the map (19) is a Lie algebra isomorphism from $\mathfrak{sl}(2,\mathbb{R})$ to $(\mathbb{R}^{3},\times_{L})$, where $a\times_{L}b=2L(a\times b)$, for any $a,b\in\mathbb{R}^{3}$. Furthermore, the tensor $L$ defines also the hyperbolic inner product $a\cdot_{L}b=a\cdot(Lb)$, for any $a,b\in\mathbb{R}^{3}$. We recall that the coadjoint orbits in $\mathfrak{sl}(2,\mathbb{R})^{*}$ are hyperboloids of the form $\{x^{2}+y^{2}-z^{2}=const\}$. Hence, in analogy with the spherical midpoint method in [16], we will call (20) the hyperbolic midpoint method.

We illustrate an application of the hyperbolic midpoint method (20) on the point vortex equations on the hyperbolic plane (see for example [7, 8, 19]). The interest in this equations is motivated also by the studies on ideal hydrodynamics on hyperbolic spaces [5, 10], and in particular on the Euler equations, for which the point vortices can be seen as a finite dimensional approximation [2]. These equations are a Lie–Poisson system on $(\mathfrak{sl}(2,\mathbb{R})^{*})^{N}\cong(\mathbb{R}^{3},\times_{L})^{N}$, with initial values on the coadjoint orbit determined by the equations $w_{i}\cdot_{L}w_{i}=-1$, for $i=1,2,\dots,N$. The Hamiltonian is given by:

[TABLE]

The equations of motion are then:

[TABLE]

Equations (22) constrain the vortices to move on the hyperboloid $x^{2}+y^{2}-z^{2}=-1$. Furthermore, the $SL(2,\mathbb{R})$ symmetry of (22) gives the conservation of the momentum vector:

[TABLE]

Equations (22) and their $SL(2,\mathbb{R})$-relative equilibria have been studied in [7, 8, 19]. In particular, for two and three vortices most of the stability issues have been worked out in [19]. However, unlike to point vortex equations on a sphere [12], it is still unknown a general result on the stability of a relative equilibrium of point vortices on the hyperbolic plane.

Here we present the results of some numerical simulations of (22) with the hyperbolic midpoint method (20). In particular, we consider as initial values some of the relative equilibria found in [8, 19]. Let us take two different initial values in $(\mathbb{R}^{3})^{3}$, $w_{0}^{1}$ and $w_{0}^{2}$, whose columns represent the initial position of three vortices with strengths respectively equal to $\Gamma^{1}$ and $\Gamma^{2}$, as defined here below:

[TABLE]

[TABLE]

The initial condition $w_{0}^{1}$ is an equilateral relative equilibrium, whereas $w_{0}^{2}$ is a geodesic relative equilibrium, as defined in [19]. The fist initial condition is known to be stable, whereas it is not known for the second one. However, from figure 6 we can see that both the initial conditions evolve in close trajectories, which proves numerically the stability for both of them.

We conclude showing in figure 7 the conservation properties for the hyperbolic midpoint method (20), concerning the first integrals (23) and (21).

Acknowledgements The author would like to thank Klas Modin for the support and the enlightening discussions during the work on this paper.