A minimal-variable symplectic method for isospectral flows

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 spherical midpoint method on so(3)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathfrak {s}}}{{\mathfrak {o}}}(3)$$\end{document}. 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.


Introduction
The numerical integration of isospectral flows is a classical subject of study in numerical analysis [7,10]. 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 [7]. 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 [12,Prop. 6.28]. Moreover, any real matrix Lie algebra which is closed under conjugate transpose is reductive [12,Prop. 1.56]. 1 Because of this, Lie-Poisson systems on the dual of a reductive Lie algebra can be equivalently seen as isospectral flows of the form (1.2) below. Lie-Poisson systems originate from the Poisson reduction of canonical Hamiltonian systems on the cotangent bundle of a Lie group [15]. 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,10,20]). For Lie-Poisson systems various symplectic algorithms have been developed (see [6] 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: Here, [·, ·] denotes the matrix commutator, V is a linear subspace of the Lie algebra gl(n, C), and the function B : V → n(V ) maps V into its gl(n, C)-normalizer algebra n(V ). 2 The most studied case in literature is when V = Sym(n, R) is the space of symmetric real matrices, for which the normalizer is the Lie algebra of skew-symmetric real matrices n(V ) = so(n). For Lie-Poisson systems on the dual of a reductive Lie algebra, we have that V = g * is the dual of a reductive Lie subalgebra of gl(n, C), for which the gl(n, C)−normalizer is n(V ) = g 0 ⊕ c(g) (see Definition 2 and Lemma 2, in Section 2). Throughout the paper, we identify g * with g, via the Frobenius inner product A, B = Tr(A † B), where † is the conjugate transpose. Via this identifications, Lie-Poisson systems on the dual of a reductive Lie algebra g take the form: for H : g → C, a smooth function called Hamiltonian.
A class of numerical methods to solve (1.1)-(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(·) 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 R 3 . Finally, we show how our scheme looks on sl(2, R), defining what we call the hyperbolic midpoint method.

Main result
Let us consider an isospectral flow of the form (1.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.1) to be well defined, we have to require B(·) to take values in the gl(n, 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 g its Lie algebra. Furthermore, let V ⊆ g be a linear subspace. Then the two sets are respectively called the G-normalizer and the g-normalizer of V . Notice that N (V ) is a subgroup of G and n(V ) is a Lie subalgebra of g.
A related concept to normalizer is the centralizer Lie algebra. Definition 2 Let g be a Lie algebra and let V ⊆ g be a linear subspace. Then the set is called the g-centralizer of V . Notice that c(V ) is a Lie subalgebra of g.
We now recall the definition of a J -quadratic Lie algebra. Definition 3 A Lie subalgebra g of gl(n, C) is called J -quadratic Lie algebra if there exists an invertible matrix J , such that (2.1)

Lemma 1
Let g be a Lie subalgebra of gl(n, C) such that there exists a matrix J for which Then g = g † implies J 2 ∈ c(g). Moreover, if J is invertible and J 2 ∈ c(g), then g = g † .
Proof Suppose g = g † . Then, for all W ∈ g, both the following identities hold: The second of these implies that W J 2 + J W † J = 0 and the first one J W † J + J 2 W = 0. Subtracting these identities, we get [W , J 2 ] = 0. Hence J 2 ∈ c(g). Now assume that g is J −quadratic and J 2 ∈ c(g). Then, for all W ∈ g we have 0 = J W † J + J 2 W = J W † J + W J 2 and hence W J + J W † = 0, being J invertible. Therefore g and g † are defined by the same identity and they coincide.

Lemma 2
Let g be a Lie subalgebra of gl(n, C) such that g = g † . Then the gl(n, C)−normalizer of g is n(g) = g 0 ⊕ c(g), where g 0 is the semisimple ideal of g such that g = g 0 ⊕ z(g), for z(g) the center of g. 3 Proof Let g ⊥ be the orthogonal complement of g in gl(n, C) with respect to the Frobenius inner product. It is not hard to check that the following properties hold: Hence, if A ∈ n(g) ∩ g ⊥ it must be [g, A] = 0. Therefore, A has to be in c(g). Moreover, we have that the following inclusions always hold z(g) ⊂ c(g) ⊂ n(g). Therefore, n(g) = g 0 ⊕ c(g), being g 0 centerless.
Notice that, since always n(g) † = n(g ⊥ ), we have that n(g ⊥ ) = g 0 ⊕ c(g) † , whenever g = g † . In conclusion, we have the following corollary: Then the gl(n, C)−normalizer of g and g ⊥ are respectively n(g) = g 0 ⊕ c(g) and n(g ⊥ ) = g 0 ⊕c(g) † . In particular, under the identification of g * with g via the Frobenius inner product, any Lie-Poisson system on g * can be written in the form (1.2).
Notice that by [12,Prop. 1.56] any 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 (1.2). In particular, this is true for Lie-Poisson systems on the dual of g ⊕ Z, where Z is an Abelian Lie algebra. We can now state the main result of this paper.
Assume that the normalizer splits as n( is a second order isospectral integrator for (1.1), for any k ≥ 0. 4 . Moreover, when (1.1) is a Lie-Poisson system on gl(n, C) * or on the dual of some J -quadratic Lie algebra g such that J 2 ∈ c(g) (or even on g ⊕ Z, where Z is an Abelian Lie algebra), then (2.3) is a Lie-Poisson integrator for (1.1) which preserves the coadjoint orbits in g * .

Remark 1
The main contribution of Theorem 1, with respect to the results presented in [18], is that the scheme (2.3) is a minimal-variable isospectral (Lie-Poisson) integrator. Minimal-variable means here that the only unknown is W , which lives in a vector space of dimension dim(V ). Hereafter, we will refer to the scheme (2.3) as the isospectral minimal midpoint. Moreover, the proof of the properties of (2.3), 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 (2.2) 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.
has a solution X ∈ gl(n, C) for any 0 ≤ h < h.
Proof In order to get a solution to (2.4), we consider the function In order to determine h, we consider the initial value problem: ) is continuous and invertible, then the Peano existence theorem will ensure a solu- Lemma 4 Let B : D ⊂ gl(n, C) → gl(n, C) continuously differentiable in the domain D and let 0 ≤ h < h as in Lemma 3. Then, for every W k ∈ D, the numerical scheme where Cay is the Cayley transform. Therefore we have that: Assuming n(V ) = g 0 ⊕ c(V ), for some Lie algebra g 0 which satisfies (2.2), by [7, Lemma IV.8.7], Cay h 2 B( W ) is in the normalizer group N (V ) of V and therefore W k+1 is in V as well. When V = g * for g a J -quadratic Lie algebra such that J 2 ∈ c(g), the transformation (2.5) coincides with the coadjoint action of G on g * , where G is the respective connected component to the identity of a Lie group with Lie algebra g. Therefore, (2.5) fixes the coadjoint orbits.

Remark 2
We point out that the equation (2.5) 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 * G L(n, C) to gl(n, 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. Proof Since g 0 is a compact Lie algebra, the associate connected Lie group G is compact and therefore the orbits of the action (2.5) are compact. Hence, we can find a minimum h > 0 in Lemma 3 independent from the iteration k ≥ 0.
for k ≥ 0 with unknowns X , Y , K ∈ gl(n, C). It is not hard to check that the following identities hold:

Proof. (Theorem 1)
The proof simply follows from the lemmas. Lemma 3 says that the method (2.3) has solution for h sufficiently small. Under the assumptions of Lemma 4 proves the isospectrality of the scheme and its intrinsic restriction V , when n(V ) = g 0 ⊕ c(V ), for some Lie algebra g 0 which satisfies (2.2). 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 gl(n, C) * , when (1.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 ∈ c(g) (possibly plus a commutative Lie algebra) and (1.1) is Lie-Poisson, the scheme (2.3) is a Lie-Poisson integrator on V , which preserves the coadjoint orbits.

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

Remark 4
The isospectral minimal midpoint (2.3) can be derived in a different way, as proposed in [6]. The construction there is more general and (2.3) 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, . . . and a s × s real matrix, a retraction map τ a : g ⊕s → G ×s . Then, similarly to [6], for each i = 1, . . . , s, it is implicitly defined by the differential of the retraction map dτ a a discrete map W k → W i g. Finally, we define our integrator Ψ h : W k → W k+1 as: for some real numbers b i . The question is whether, for any s-stage symplectic Runge-Kutta method, there exists a retraction map τ a : g ⊕s → G ×s , such that any Lie-Poisson integrator defined in [18] can be obtained in this way.

Numerical examples
In this section we present some applications of the isospectral minimal midpoint (2. For the example considered in this section, we plot the variation of the first integrals of (1.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.

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 so(n), for every n ≥ 1, [14]. The Hamiltonian is given by where I : so(n) → so(n) is a symmetric positive definite inertia tensor. The equations of motion are thenẆ We discretize this system for n = 10. Our implementation uses Newton iterations for the non-linear system. 5 The inertia tensor is given by As shown in Fig. 1, the Hamiltonian is nearly conserved and the Casimir functions are conserved up to the accuracy of the Newton iterations.

The Brockett flow
In this section we specify the isospectral minimal midpoint where N , W are n × 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 → ∞, W (t) converges exponentially fast to a diagonal matrix with the eigenvalues sorted accordingly to the order of the entries of N . In Fig. 2, we plot the eigenvalue variation for a randomly generated 6 self-adjoint initial matrix W 0 of dimension 10 × 10 and N = diag (1, 2, . . . , 10). The asymptotic stationarity of the eigenvalues variation reflects the fact that W is close to a diagonal matrix.

Lie-Poisson systems on (R 3 , ×)
On (R 3 , ×) the isospectral minimal midpoint (2.3) can be written as: for w, w k , w k+1 ∈ R 3 and B : R 3 → R 3 . We want to compare the isospectral minimal midpoint (3.6) with another minimal-variable symplectic integrator on R 3 introduced in [17], i.e. the spherical midpoint method: (3.7) Remark 5 Let B(·) be orthogonal with respect to the rays, i.e. B(w) · w = 0, for every w ∈ R 3 . It is immediate to check that then (3.6) coincides with the classical midpoint scheme: In [17] it is shown that also (3.7) coincides with (3.8) when B(·) = ∇ H (·), for some Hamiltonian function H : R 3 → R constant on the rays (which implies B(·) to be orthogonal to the rays). In this case, (3.8) is known to be symplectic, whereas this fails for general Hamiltonian H . Therefore, (3.6) can be seen as the second order correction of (3.8) to be symplectic for any Hamiltonian H .
Let us now consider the two schemes (3.6) and (3.7). 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 R 3N . For this one has to extend both the isospectral minimal midpoint (3.6) and the spherical midpoint (3.7) to direct products of R 3 (see [17,18]). The Heisenberg spin chain of micromagnetics is defined as: where w i ∈ S 2 , for i = 1, . . . , N and w N +1 = w 1 . It corresponds, up to scaling, to spatial discretization of the Landau-Lifschitz PDE: for w : S 1 → S 2 smooth. We notice that (3.9) is a Lie-Poisson system on R 3N , with Hamiltonian: Clearly, to get a good approximation of (3.10), N has to be large. In Fig. 3 we show the average time cost for time-step with respect to the number of spin particles for both the isospectral minimal midpoint (3.6) and the spherical midpoint (3.7). 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 i, j w † i Aw j . Moreover, the spherical midpoint has the advantage of exactly conserving all the quadratic first integrals of the form w † i Aw i , for some square matrix A, whereas the isospectral minimal midpoint (3.6) exactly conserves only the quadratic invariants w † i w i . In Fig. 4 we compare the isospectral minimal midpoint (3.6) and the spherical midpoint (3.7) for the initial data given as am equispaced discretization in N = 100 points of the closed spherical curve: for x ∈ [0, 1]. We can conclude from Fig. 4 that the spherical midpoint performs slightly better than the isospectral minimal midpoint (3.6). However, both the schemes show the desired conservation properties due to their symplecticity. Finally, in Fig. 5, we present the error diagram for the isospectral minimal midpoint (3.6) and the spherical midpoint (3.7). The curves in Fig. 5 show the expected second order of the schemes, with no significant difference.

Lie-Poisson systems on sl(2, R) *
In this section we specify the isospectral minimal midpoint (2.3) in the case of the Lie algebra sl(2, R), i.e. the 2 × 2 matrices with zero trace. The first observation is that sl(2, R) ∼ = sp(2, R), which is a non-compact J -quadratic Lie algebra with respect to J = 0 −1 1 0 . From this, it is straightforward to see that, when B : gl(2, R) → sl(2, R) and W k ∈ sl(2, R), W in (2.3) is also in sl(2, R). We notice that this is no longer true for sl(n, R), for n > 2. On the other hand, any element in sl(2, R) can be written as a vector in R 3 , via the vector spaces isomorphism: for any x, y, z ∈ R. In this coordinates, we can express the isospectral minimal midpoint (2.3) for sl(2, R) as: for any a, b, c ∈ R 3 . We notice that the map (3.11) is a Lie algebra isomorphism from sl(2, R) to 2L(a × b), for any a, b ∈ R 3 . Furthermore, the tensor L defines also the hyperbolic inner product a · L b = a · (Lb), for any a, b ∈ R 3 . We recall that the coadjoint orbits in sl(2, 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 (3.12) the hyperbolic midpoint method.
We illustrate an application of the hyperbolic midpoint method (3.12) on the point vortex equations on the hyperbolic plane (see for example [8,9,19]). The interest in this equations is motivated also by the studies on ideal hydrodynamics on hyperbolic spaces [5,11], 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 (sl(2, R) * ) N ∼ = (R 3 , × L ) N , with initial values on the coadjoint orbit determined by the equations w i · L w i = −1, for i = 1, 2, . . . , N . The Hamiltonian is given by: The equations of motion are then: (3.14) Equations (3.14) constrain the vortices to move on the hyperboloid x 2 + y 2 − z 2 = −1. Furthermore, the SL(2, R) symmetry of (3.14) gives the conservation of the momentum vector: Equations (3.14) and their SL(2, R)-relative equilibria have been studied in [8,9,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 [13], 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 (3.14) with the hyperbolic midpoint method (3.12). In particular, we consider as initial values some of the relative equilibria found in [9,19]. Let us take two different initial values in (R 3 ) 3 , w 1 0 and w 2 0 , whose columns represent the initial position of three vortices with strengths respectively equal to Γ 1 and Γ 2 , as defined here below: The initial condition w 1 0 is an equilateral relative equilibrium, whereas w 2 0 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 Fig. 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 Fig. 7 the conservation properties for the hyperbolic midpoint method (3.12), concerning the first integrals (3.15) and (3.13).