Lie-Poisson methods for isospectral flows

The theory of isospectral flows comprises a large class of continuous dynamical systems, particularly integrable systems and Lie--Poisson systems. Their discretization is a classical problem in numerical analysis. Preserving the spectra in the discrete flow requires the conservation of high order polynomials, which is hard to come by. Existing methods achieving this are complicated and usually fail to preserve the underlying Lie--Poisson structure. Here we present a class of numerical methods of arbitrary order for Hamiltonian and non-Hamiltonian isospectral flows, which preserve both the spectra and the Lie--Poisson structure. The methods are surprisingly simple, and avoid the use of constraints or exponential maps. Furthermore, due to preservation of the Lie--Poisson structure, they exhibit near conservation of the Hamiltonian function. As an illustration, we apply the methods to several classical isospectral flows.


Introduction
Lie-Poisson systems and isospectral flows are two well-studied classes of dynamical systems. The former appear as Poisson reductions of Hamiltonian systems for which the configuration and symmetry space is a Lie group (see the monograph [17] and references therein). The classical example is the free rigid body as viewed by Poincaré [26]. The latter, isospectral flows, appear as Lax formulations of integrable systems (see the survey papers [31,5,30] and references therein). The classical example is the Toda lattice as viewed by Flaschka [29,8].
The study of numerical methods for the two classes of systems are by now classical subjects in numerical analysis. The motivation for such schemes came through the strong connection between matrix factorizations in numerical linear algebra and isospectral flows (see the survey papers [6,23]). This was initiated by the remarkable discovery that the iterative QR-algorithm for computing eigenvalues is a discretization of the (non-periodic) Toda flow [28,7].
The general form of an isospectral flow iṡ Here, [·, ·] denotes the matrix commutator, S is a linear subspace of the Lie algebra gl(n, C), and the function B : S → n(S) maps into the normalizer algebra n(S) (see Section 3 below for details). The most studied setting is when S = Sym(n, R) is the space of symmetric real matrices, for which the normalizer is the Lie algebra of skew-symmetric real matrices n(S) = so(n). Another setting is when S = g is a Lie subgroup of gl(n, C), for which the normalizer is the subalgebra itself n(S) = g.
Let us now discuss the connection between isospectral flows and Lie-Poisson systems. The predominant example connecting the two is Manakov's n-dimensional rigid body [16].
Recall that a Lie-Poisson system evolves on the dual g * of a Lie algebra g. Given a Hamiltonian function H on g * , the flow W (t) ∈ g * is given bẏ where the operator ad * is defined by Without loss of generality we may assume that g is a subalgebra of gl(n, C). To identify gl(n, C) * with gl(n, C) we use the Frobenius inner product where W † denotes the conjugate transpose. In this way we also identify g * with the subspace g ⊂ gl(n, C). Next we extend the Hamiltonian to all of gl(n, C) by taking it to be constant on the affine spaces given by translations of the orthogonal complement of g. Then dH corresponds to ∇H. From the definition (3) and the identification of g * with g we get ad * W (M ) = Π [W † , M ], where Π is the orthogonal projection gl(n, C) → g. We thus arrive at an explicit formulation of the Lie-Poisson system (2), namelẏ Now, the key observation is that if the representation of g as a subalgebra of gl(n, C) is closed under conjugate transpose, then equation (4) becomes the isospectral flowẆ = [∇H(W ) † , W ]. (5) Such a representation is possible if and only if g is a reductive Lie algebra (see Section 2-3 below for details). Thus we arrive at the statement that Lie-Poisson systems for any reductive Lie algebra can be viewed as isospectral flows. Recall that most classical Lie algebras are reductive, for example gl(n, C), gl(n, R), sl(n, C), sl(n, R), u(n), su(n), so(n), and sp(n).
An interesting consequence of equation (5) is that whenever the function B(W ) in the isospectral flow (1) can be written as B(W ) = ∇H(W ) † , then it can be extended to a Lie-Poisson system on gl(n, C) (or possibly a smaller reductive Lie algebra containing S). Indeed, just extend the Hamiltonian function H to be constant of the affine fibres orthogonal to S. In this way we obtain an extended system foliated into invariant affine subspaces generated by S. The Toda flow is an example where this construction is possible (see Section 5.2 below).
The key feature of isospectral flows is, of course, that the eigenvalues of W are preserved. Equivalently, given any analytic function f extended to matrices, the function F (W ) = Tr(f (W )) is a first integral regardless of the choice of B(W ) in (1). From the perspective of Lie-Poisson systems (5), this means that F (W ) is a Casimir function associated with the Lie-Poisson structure (3). Although there are infinitely many Casimir functions, only a finite number of them can be functionally independent.
In this paper, we develop spectral preserving numerical methods for flows of the form (1) which, in the case of Hamiltonian isospectral flows (5), also preserves the Lie-Poisson structure. There already exist at least four ways to achieve this: • If the Hamiltonian can be written as a sum of explicitly integrable Hamiltonians one can use splitting method (see [22] and references therein). • The Lie-Poisson system on g * g can be extended to a constrained canonical Hamiltonian system on T * G T G ⊂ T GL(n, C). One can then use the symplectic RATTLE method (or higher order versions of it) for the constrained system (see [13,21]).
• One can use symplectic Lie group methods on T * G as developed in [3].
These methods rely on an invertible mapping between the Lie algebra and (an identity neighbourhood of) the Lie group, such as the exponential map (works in general) or the Cayley map (works for quadratic Lie groups). • One can, in some cases, use collective symplectic integrators, which rely on Clebsch variables originating from a Hamiltonian action of G on a symplectic vector space (see [18,19] for details). Compared to these methods our approach is: (i) simpler since the algorithms are formulated directly on the algebra g ⊂ gl(n, C); (ii) free of constraints; (iii) free of algebra-to-group maps, such as the exponential or Cayley map; (iv) generic as they apply to any isospectral Hamiltonian flow. Furthermore, through the framework of Poisson reduction (cf. [17]) our methods are directly related to classical symplectic Runge-Kutta methods (or partitioned symplectic Runge-Kutta methods). Therefore they merit the designation Isospectral Symplectic Runge-Kutta (IsoSyRK) methods.
The paper is organized as follows. In Section 2 we give the definitions of the new methods and we state our main results. In Section 3 we develop a discrete reduction theory for isospectral Lie-Poisson flows. These results are instrumental in Section 4, where we specialize our construction to symplectic Runge-Kutta methods. All numerical examples are given in Section 5.

Main results
A Runge-Kutta method is defined by its Butcher tableau (cf. [9]) for i, j = 1, . . . , s, then the corresponding Runge-Kutta method is symplectic when applied to canonical Hamiltonian systems on R 2n [27]. However, directly applying a symplectic Runge-Kutta method to the Hamiltonian isospectral flow (5) does not yield a Poisson integrator. Nor does it, in general, preserve the isospectral property, as is well known.
Theorem 1. The method in Definition 1 fulfills the following properties: (1) It has the same order as the underlying Runge-Kutta method.
(3) It is equivariant with respect to Lie algebra morphisms; if A : gl(n, C) → gl(n, C) is a linear invertible mapping fulfilling for all X, Y ∈ gl(n, C) then the following diagram commutes It is a Lie-Poisson integrator if the isospectral flow is Hamiltonian, i.e., of the form (5). Furthermore, if the isospectral flow is Hamiltonian, it is equivariant with respect to linear Lie-Poisson isomorphisms B : gl(n, C) * → gl(n, C) * . (5) It restricts to a Lie-Poisson integrator for any Lie subalgebra g ⊂ gl(n, C) defined by where J 2 = cI for some c ∈ R\{0}. (6) It restricts to a Lie-Poisson integrator for any Lie subalgebra given by arbitrary intersections of gl(n, R), sl(n, C), and Lie algebras of the form (8). (7) It extends to a Lie-Poisson integrator for direct products of Lie algebras of the form in item (6). (8) It restricts to an isospectral integrator on the orthogonal complement g ⊥ ⊂ gl(n, C) of any Lie algebra g of the form in item (6), provided that B restricts to a mapping B : g ⊥ → g.
Proof. The theorem is a combination of results proved in Theorem 3, Corollary 1, Theorem 4, and Theorem 6 below.
We also have an analogous, albeit slightly weaker, result for partitioned symplectic Runge-Kutta methods, such as defined by two Butcher tableaux (cf. [9]) If, for i, j = 1, . . . , s, the coefficients in the tableaux fulfill then the corresponding partitioned Runge-Kutta method is symplectic when applied to canonical Hamiltonian systems on R 2n .
Theorem 2. The method in Definition 2 fulfills the following properties: (1) It has the same order as the underlying partioned Runge-Kutta method.
(3) It is equivariant with respect to Lie algebra morphisms; if A : gl(n, C) → gl(n, C) is a linear invertible mapping fulfilling for all X, Y ∈ gl(n, C) A[X, Y ] = [AX, AY ], then the following diagram commutes It is a Lie-Poisson integrator if the isospectral flow is Hamiltonian, i.e., of the form (5). Furthermore, if the isospectral flow is Hamiltonian, it is equivariant with respect to linear Lie-Poisson isomorphisms B : gl(n, C) * → gl(n, C) * . (5) If it restricts to a Lie-Poisson integrator for a Lie subalgebra g ⊂ gl(n, C) defined by W ∈ g ⇐⇒ W † J + JW = 0, where J 2 = cI for some c ∈ R\{0} and b i = 0 for i = 1, . . . , s, then the two Butcher tableaux coincide (it is a standard Runge-Kutta method).
Proof. The theorem is a combination of results proved in Theorem 3, Corollary 1, Theorem 5, and Theorem 6 below.

Reduction theory for isospectral Lie-Poisson integrators
Let us consider Lie-Poisson systems of the form (5). Conditions under which the flow remains in a linear subspace S of gl(n, C) are intrinsically connected to the gl(n, C)-normalizer of S. We recall here its definition: Definition 3. Let G be a Lie group and g its Lie algebra. Furthermore, let S ⊆ g be a linear subspace. Then the two sets are respectively called the G-normalizer and the g-normalizer of S. Notice that N (S) is a subgroup of G and n(S) is a Lie subalgebra of g.
We now give some examples of Definition 3. We first give the following definition.
Definition 4. Let g ⊆ sl(n, C) be a Lie algebra and J ∈ GL(n, C). Then g is said to be a J-quadratic Lie algebra if A † J + JA = 0, for any A ∈ g.
Remark 2. If S = g is a Lie subalgebra of gl(n, C) one may ask under which conditions the isospectral Hamiltonian system (5) coincide with the Lie-Poisson system on g. Recall that ∇H is the gradient of H with respect to the Frobenius inner product. It is not a restriction to assume that ∇H(W ) ∈ g for all W ∈ gl(n, C) (since we can extend H to be constant on the affine complements of g). Due to the conjugate transpose on ∇H(W ) it is not, however, enough that ∇H(W ) ∈ g; instead we need ∇H(W ) † ∈ g. A sufficient condition for this to be true is that g is closed under conjugate transpose: g † ⊂ g. Such g are, up to representation, the semisimple Lie algebras ( [14], Prop. 6.28). This means that, after the identification of the dual of the Lie algebra g with itself (using the Frobenius inner product), a Lie-Poisson system on a semisimple Lie algebra g coincides with a Lie-Poisson system on gl(n, C) restricted to g. In fact, slightly more is true: due to the bracket in the right hand side of (5) it is enough that This holds when g is a reductive Lie algebra, i.e., the direct sum of a semisimple Lie algebra and an abelian Lie algebra.
A Lie-Poisson systems on a Lie algebra g can be viewed as the Lie-Poisson reduction of a canonical Hamiltonian system on T * G with a G-symmetric Hamiltonian. Going backwards, one may 'unreduce' any Lie-Poisson system to a canonical Hamiltonian system on the cotangent bundle of the corresponding Lie group. Our 1 Orthogonal complements are taken with respect to the Frobenius inner product. We always have that n(S ⊥ ) = n(S) † . objective is to show that the Isospectral Symplectic Runge-Kutta methods (cf. Section 2) originates as a "discrete Lie-Poisson reduction" of symplectic Runge-Kutta methods.
We thus proceed by extending the equations (5) to a canonical system on T * GL(n, C). To do so, one needs the momentum map (cf. [17]) associated with the right action of GL(n, C) on T * GL(n, C): for G, Q ∈ GL(n, C) and P ∈ T * Q GL(n, C). The momentum map for this (Hamiltonian) action is given by This momentum map provides a left-invariant Hamiltonian function H(Q, P ) = H(Q † P ) on T * GL(n, C), i.e., an Hamiltonian function invariant with respect to the left action The fact that the momentum map is a Poisson map between T * GL(n, C) and gl(n, C) * means that a symplectic map in Φ : T * GL(n, C) → T * GL(n, C) which is equivariant with respect to the action (11) descends to a corresponding map φ : gl(n, C) * → gl(n, C) * . In terms of numerical integrators, this means that a GL(n, C)-equivariant symplectic integrator on T * GL(n, C) induces a Poisson integrator on gl(n, C) * . As we shall see, this is precisely how the isospectral symplectic Runge-Kutta methods come about.
Using the momentum map (11), the canonical Hamiltonian system on T * GL(n, C) is given byQ where H is the same Hamiltonian as in (5). We now translate the condition of staying on S from (5) to (12).
Consider a solution (Q(t), P (t)) of Hamilton's equations (12) for a given initial point (Q(0), P (0)) and let S ⊆ gl * (n, C) be a linear subspace as before. Then there exists a time T > 0 such that the following three statements are equivalent: , since S is a linear space and (Q † P )(t) ∈ S, for any 0 ≤ t ≤ T . But this means that ∇H(Q † P ) † has to be in n(S), for any 0 ≤ t ≤ T .
2) ⇒ 1) For 0 ≤ t ≤ T , we have that: which proves the statement, since N (S) ⊇ exp(n(S)). 2) ⇒ 3) Let G ∈ GL(n, C) such that GQ(t) † ∈ N (S). Then we have: 3) ⇒ 2) By the formula above we have Although the statements in Proposition 1 are equivalent for the exact flow, they are different after discretization. Indeed, in order to understand the conditions for our isospectral symplectic Runge-Kutta methods to preserve the flow on S we need the definition of weak and strong first integrals.
Definition 5. Let M be a smooth manifold and N ⊂ M a smooth submanifold. Consider the following dynamical system on N : with X a smooth vector field on N and z 0 ∈ N . Assume further that X can be extended on a ε−neighbourhood N ε of N in M . Then a differentiable function I : N ε → C is said to be a weak, respectively, strong first integral of (13) if dI(z), X(z) = 0 for all z ∈ N dI(z), X(z) = 0 for all z ∈ N ε .
In numerical analysis it is often the case that integration schemes on a submanifold N actually depends on how N is embedded in a larger (vector) space M . That is, the integration scheme is not intrinsic to N (for example evaluations of the vector field outside of N may occur). In this situation, the difference between strong and weak first integrals is essential. Indeed, for non-intrinsic methods one can at best expected to conserve strong first integrals. Motivated by this we make the following: Since S is a linear space the natural way to extend ∇H † is to take it to be constant on the affine complements of S. With this extension the gradient of the Hamiltonian requires only an orthogonal projection of W to S.
Under Assumption 1 our Proposition 1 says that (Q † P ) ∈ S is determined by weak first integrals of the Hamiltonian system (12) provided that the gradient of the Hamiltonian is in n(S). In fact, having (Q † P ) ∈ S is equivalent to [∇H(Q † P ) † , Q † P ] ∈ S which in general is not true for Q † P in an ε−neighbourhood of S. Instead an equivalent formulation corresponding to strong first integrals is given by the third statement, which says that there exists a fixed matrix G such that GQ † ∈ N (S). Therefore only the numerical methods that have GQ † ∈ N (S) as a discrete invariant correspond to integrators that preserve S. In particular, if N (S) is a quadratic Lie group one can expect symplectic Runge-Kutta methods to yield a discrete flow that preserves S since they preserve general quadratic first integrals. On the other hand, the same cannot be expected from symplectic partitioned Runge-Kutta methods, since they only preserve special (bilinear) quadratic first integrals.
We summarize our findings in the following theorem.
Theorem 3. Consider a Lie-Poisson system of the form (5) evolving on a linear subspace S ⊂ gl(n, C). Let Φ h : T * GL(n, C) → T * GL(n, C) be a symplectic numerical method for the corresponding canonical Hamiltonian system (12) obtained by extension from S in accordance with Assumption 1.
(1) If Φ h is equivariant with respect to the action (11), i.e., then it descends to a Lie-Poisson integrator φ h on gl(n, C).
then φ h restricts to an integrator on S.
Based on the results in Theorem 3, we can now generalize the results to a general B(·), i.e., to isospectral flows that are not necessarily Hamiltonian. This extension requires that the underlying method can be expanded in a B-series or P-series (cf. [9] for definitions and notation). 2 Consider first the generalization of Assumption 1: Assumption 2. Let S ε be a ε−neighbourhood of S in gl(n, C). We assume that B(·) can be extended to S ε such that B(W ) ∈ n(S) for all W ∈ S ε .
Then, based on Theorem 3, we have the following result. Corollary 1. Consider an isospectral flow of the form (1) evolving on a linear subspace S ⊂ gl(n, C). Let Φ h : T * GL(n, C) → T * GL(n, C) be a symplectic Bseries (or P-series) method for the corresponding system: obtained by extension from S in accordance with Assumption 2.
(1) If Φ h is equivariant with respect to the action (11), i.e., then it descends to an isospectral integrator φ h on gl(n, C).
then φ h restricts to an integrator on S.
Proof. From Theorem 3 we know that for B(W ) = ∇H(W ) † when we solve (14) with a symplectic integrator the discrete flow is isospectral for W := Q † P . Therefore, the (truncated) modified equation is of the forṁ for some modified Hamiltonian H. On the other hand, since Φ h is a symplectic B-series method, the right hand side in (15) is a B-series whose coefficients satisfy  (15) is of the form for A k n , B k n homogeneous polynomials of degree n for each k. We claim that to get the modified equation for a general B we just replace in (16) (H (k) ) † with B (k−1) (which is possible since k ≥ 1). Indeed, this follows since the coefficients of a symplectic B-series are uniquely determined by Hamiltonian vector fields [9, Thm IX.9.10]. Therefore, we conclude that a symplectic B-series integrator applied to (14), for a general B, is isospectral for W = Q † P with a modified equation of the formQ for some B(·) obtained by replacing in (16) the (H (k) ) † with B (k−1) . In the case when Φ h is a symplectic P-series, the proof is repeated similarly, using instead [9, Thm IX.10.3, Lem IX.10.6, Thm IX.10.8].

Isospectral symplectic Runge-Kutta methods
In this section we specialize Theorem 3 to the symplectic Runge-Kutta and partitioned Runge-Kutta methods. As a result, we obtain the novel numerical schemes for isospectral (Lie-Poisson) systems presented in Section 2 above.  (1) is for i, j = 1, . . . , s. Recall that the method is symplectic, i.e., the discrete flow is a symplectic map, if b i a ij + b j a ji = b i b j for any i, j = 1, . . . , s. (1) The symplectic integrator Φ h descends to a Lie-Poisson integrator φ h on gl(n, C) * gl(n, C) for the isospectral Hamiltonian system (5). Furthermore, the map φ h is completely constructive as an implicit integration scheme (see below for specific formulas).
(2) If S is an invariant subspace of (5) (as described above), then φ h preserves S in the cases S = sl(N, C), S = g, and S = g ⊥ , for g a J-quadratic Lie subalgebra.
The schemes obtained in Theorem 4 are the following: 1. S = sl(n, C) or S = gl(n, C).
, for i, j = 1, . . . , s, where the unknowns are X i , Y i , K ij for i, j = 1, . . . , s and the last two lines are explicit.
. . , s, where the unknowns are X i , K ij for i, j = 1, . . . , s and the last two lines are explicit. The last line is also equivalent to 3. S = g ⊥ , g ⊆ sl(n, C) J-quadratic.
, for i, j = 1, . . . , s, where the unknowns are X i , K ij for i, j = 1, . . . , s and the last two lines are explicit. The last line is also equivalent to Proof of Theorem 4. hej (1) For S := sl(n, C) we have that n(S) = gl(n, C) and N (S) = GL(n, C). Therefore the hypotheses of Theorem 3 are trivially satisfied. To get the explicit construction, we look at the argument of the gradient of the Hamiltonian which suggests to define for i, j = 1, . . . , s. The equations for X i , Y i are straightforward (consider the equations of the Runge-Kutta method (18) and take the transpose of the first equation and multiply by P n , and multiply the second equation by Q † n , respectively).
To get the equations for K ij , we first transpose the first equation of (18), then we multiply it (indexed now by j ) by h 2 a ij K P j and sum over j . We thereby get Multiplying the second equation of (18) (indexed now by j ) by h 2 a ij (K Q j ) † and then summing over j , we obtain for i, j = 1, . . . , s.

Using then
a jj K ij for i, j = 1, . . . , s.
Therefore the K ij depend completely on the K ij and so we can neglect them, obtaining the desired equations for the K ij . Finally, to get the equation for W n+1 , we multiply the third one of (18) transposed with the fourth one of (18) and we get Using the symplecticity of the method, the last term becomes Therefore, Now substituting the equations found for X i , Y i , K ii , K ii we get the desired equation for W n+1 . (2) Symplectic Runge-Kutta methods preserve exactly the strong quadratic first integrals of a dynamical system. In particular, when S is one of the spaces stated in the theorem, they preserve N (S) = {Q ∈ GL(n, C)|Q † JQ = J}. Therefore, by Theorem 3, they descend to an integrator on S. It is also easy to check that, if we assume B † to be in n(S), we get Y i = −J −1 X † i J. Moreover, from the definition of K ij and K ij and the equations: (a ij X j + a jj K ij )) for i, j = 1, . . . , s, we get also that Remark 3. We stress that our methods are not intrinsically formulated on S. That is, they depend on how S is embedded as a subspace in gl(n, C). Therefore, there is no hope to present the schemes above only in terms of the matrix commutator.   of a symplectic partitioned s-stage Runge-Kutta method, let Φ h : T * GL(n, C) → T * GL(n, C) denote the corresponding integrator map for the system (12). Then: (1) The symplectic integrator Φ h descends to a Lie-Poisson integrator φ h on gl(n, C) * gl(n, C) for the isospectral Hamiltonian system (5). Furthermore, the map φ h is completely constructive as an implicit integration scheme (see below for a specific formula).
(2) If S = sl(N, C) is an invariant subspace of (5) (as described above), then φ h preserves S. (3) If S = g, and S = g ⊥ , for g a J-quadratic Lie subalgebra and b i = 0, then φ h preserves S (for general Hamiltonians on S extended to gl(n, C)) if and only if a ij = a ij , for i, j = 1, . . . , s.
The scheme obtained in Theorem 5 is the following: S = sl(n, C) or S = gl(n, C).
for i, j = 1, . . . , s, where the unknowns are X i , Y i , K ij for i, j = 1, . . . , s and the last two lines are explicit.
Proof of Theorem 5. hej (1) The proof is, mutatis mutandis, identical to the one of the previous theorem.
We have just to change accordingly the following definitions: and pointing out the following identity: Finally we just use the condition of symplecticity for partitioned Runge-Kutta methods. (2) Follows directly from the formula for W k+1 .  (1) and (5) are not affine equivariant; linear equivariance is the best we can expect. The linear isomorphisms that leave equations (1) and (5) invariant in form are, respectively, Lie algebra isomorphisms and Lie-Poisson isomorphisms. Indeed, consider a Lie algebra isomorphism A : gl(n, C) → gl(n, C). Applying A to equation (1) gives which shows the invariance in form of equation (1) to Lie algebra isomorphism. In particular, we have the identity Via the identification gl(n, C) * gl(n, C) as previously explained, it is easy to check that the adjoint operator A * : gl(n, C) * → gl(n, C) * acts on the coadjoint represen- , for X ∈ gl(n, C) and Y ∈ gl(n, C) * . In particular, A * is a Lie-Poisson map for equation (5): which leaves equation (5) invariant in form. We thus obtain the following identity: For any map B and Hamiltonian H, let φ h (B) and φ h (H), respectively, denote integrators as in Theorem 4 or Theorem 5. Now, the numerical scheme φ h (B) is Lie equivariant, and, correspondingly, The identities (19) and (20) show that the right hand sides of equations (21), (22) have the same form. Therefore, it is enough to prove equation (21).
Proof. Let us consider equation (21) for the partitioned symplectic Runge-Kutta schemes. The same conclusion for the symplectic Runge-Kutta method will follow straightforwardly from this. We want to check equation (21) for any W n ∈ gl(n, C) and A as above. The right hand side is , for i, j = 1, . . . , s which is exactly the left hand side of (21).

Numerical examples
In this section we demonstrate the Isospectral Symplectic Runge-Kutta methods on some Hamiltonian isospectral flows often seen in the literature. As expected, we obtain near conservation of the Hamiltonian (owing to the symplectic quality) and exact conservation (up to round-off errors) of the Casimir functions (owing to the isospectral quality). 3 5.1. The generalized rigid body. The core example among Hamiltonian isospectral systems is the generalized rigid body. It is known that in any dimension n it forms a complete integrable system in so(n), as proved by Manakov [16]. The Hamiltonian is given by where I : so(n) → so(n) is a symmetric positive definite inertia tensor. The equations of motion are thenẆ = −[I −1 W, W ] W (0) = W 0 . We discretize this system for n = 10 with the method in Theorem 4 and with the Butcher tableau corresponding to the implicit midpoint method. Our implementation uses Newton iterations for the non-linear system. The inertia tensor is given by and we use the stepsize h = 0.1. The initial conditions are given by (W 0 ) ij = 1/10 for i < j and W † 0 = −W 0 As shown in Figure 2, the Hamiltonian is nearly conserved and the Casimir functions are conserved up to the accuracy of the Newton iterations.
The Casimirs of the generalized rigid body only constitutes n first integrals, and they are therefore not enough to obtain the integrability. The additional, non-Casimir first integrals are not exactly preserved by our methods. However, from backward error analysis combined with KAM theory (see e.g. [9]), one obtains that

Evolution of error in Hamiltonian
Evolution of error in Casimirs (eigenvalues) Figure 2. Evolution of errors for the generalized rigid body in so (10). The Casimir functions correspond to the 10 eigenvalues (which occur in pairs). The Hamiltonian is given by (23). The data for the simulation are given by: stepsize h=0.1; inertia tensor I = diag(1, . . . , 10); initial conditions (W 0 ) ij = 1/10 if i < j, the additional integrals are nearly conserved (just as the Hamiltonian is nearly conserved).

5.2.
The (periodic) Toda lattice. Among Hamiltonian integrable systems the Toda lattice is perhaps the best known and most studied example. It represents a system of particles interacting pairwise with exponential forces. The equations of motion are determined by the Hamiltonian where (q i , p i ) are canonical coordinates of the n particles. Independently, Hénon [10], Flaschka [8] and Manakov [16] proved that the Toda system is integrable when q n = q n+1 (periodic boundary conditions). This is most easily seen by providing a Lax pair formulation. Indeed, by the following change of variables one obtains an equivalent isospectral floẇ where In these coordinates the canonical Hamiltonian is simply H(L) = 2 Tr(L 2 ). So far, the mapping B(·) is defined only for matrices of the form L above. A natural extension to any matrix W ∈ gl(n, C) is Next, in order to extend (24) to a Hamiltonian isospectral flow on gl(n, C) of the form in (5), we notice that we can take as a new Hamiltonian the function The flow (5)  Since the H(L) is one of the Casimir functions of the flow, it is preserved up to the iteration tolerance, as shown in Figure 3. We notice that in general our methods does not exactly preserve the zero entries of L (although they are nearly preserved). This is because the normalizer of the subspace of the symmetric matrices with the form of L is not J-quadratic for n > 3.

5.3.
The Euler equations on a sphere. Let us briefly mention a beautiful approach for spatial discretization of the incompressible Euler equations on a sphere, which leads to a finite dimensional Hamiltonian isospectral flow. For a full account we refer to the publication [24]. 4 On the 2-sphere S 2 the hydrodynamical Euler equations for an incompressible, inviscid, and homogeneous fluid can be formulated in terms of vorticity of the velocity vector. The formulation isω where the vorticity function ω is a smooth function on S 2 with zero mean, ∆ −1 is the inverse of the Laplace-Beltrami operator on the sphere (which is invertible because the kernel consists only of the constant functions), and {·, ·} is the Poisson bracket between functions. This system is an infinite dimensional Lie-Poisson system on the dual of the algebra of divergence free vector field on S 2 . The Casimir functions are given by where f : R → R is any smooth function. Thus, there are infinitely many independent first integrals. (Although not enough integrals for the system to be integrable.) In geometric quantization theory (cf. [1]) the space of smooth functions is replaced by a Hilbert space of linear operators, and the Poisson bracket {·, ·} is replaced by the commutator [·, ·]. The aim is to obtain a construction such that, in some sense, [·, ·] approximates {·, ·}. In his PhD thesis, Hoppe [11] gave an explicit quantization of (C ∞ (S 2 ), {·, ·}) in terms of the finite dimensional Lie algebras su(n, C), such that [·, ·] → {·, ·} as n → ∞. This naturally leads to a spatial discretization of the vorticity equation (25) by simply replacing {·, ·} by [·, ·] and then working out what the corresponding discrete Laplacian ∆ n should be. An explicit formula for ∆ n was given by Hoppe and Yau [12]. The resulting spatially discretized equations thus becomeẆ where W ∈ su(n). This is an isospectral Hamiltonian system with respect to the Hamiltonian H(W ) = 1 2 Tr((∆ −1 n W ) † W ). In our paper [24] we develop and further explore a fully discrete version of (26) based on the isospectral methods in this paper. We thus obtain a discrete flow that preserves all the underlying structure of the Euler equations: conservation of Casimirs and the Lie-Poisson structure. In particular, conservation of Casimirs is essential for numerical studies of the long-time behavior of (25) and of the mechanisms behind the inverse energy cascade exclusive to 2D turbulence. 5.4. Point vortices on a sphere and the Heisenberg spin chain. As stated above, the symplectic isospectral Runge-Kutta methods are readily extended to product spaces. For example we can deal with (su(2) * ) n , where n is the number of vortices or spin particles in the point-vortices equation [25] or, respectively, the Heisenberg spin chain [20].
The Hamiltonian for point-vortex dynamics is where W 1 , . . . , W n thought of a vectors in R 3 are the positions of the point vortices and Γ 1 , . . . , Γ n the respective strengths. For the Heisenberg spin chain the Hamiltonian is where W 1 , . . . , W n are the spins of the particles and W 1 = W n+1 , see for example [20]. For these systems a new first integral arises, due to the SU (2) symmetry of the Hamiltonians for any G ∈ SU (2). The corresponding first integrals are given by the (weighted) sum of the vortices/spins As before, the Casimirs are conserved and the Hamiltonian is nearly conserved. In addition, the extra integral M is conserved up to machine precision, as can be seen in Figure 4. It can be cast as an isospectral flow (1) on Sym(n, R) with B(W ) = N W + W N . Its interest lies in its integrable structure, which is fundamentally different from that of the Toda lattice and the generalized rigid body. The Bloch-Iserles flow can be extended to a Hamiltonian isospectral flow on gl(n, R) * gl(n, R) such that Sym(n, R) is an invariant subspace, just as the Toda flow in Section 5.  The evolution of energy and the Casimirs (eigenvalues) are given in Figure 5. The integrable structure of the flow is revealed as quasi-periodicity in projections of the phase diagram, as seen in Figure 6.
5.6. The Toeplitz inverse eigenvalue problem. In this section we demonstrate that the methods introduced can be applied also to non-Hamiltonian systems. To this end, consider Chu's flow on symmetric real matrices, which is of the form (1) with Notice that if W ∈ Sym(n, R) then B(W ) ∈ so(n). The Toeplitz inverse eigenvalue problem reads as follows. Given a certain set of eigenvalues, find a symmetric Toeplitz matrix with that prescribed spectra (recall that a Toeplitz matrix is a matrix with constant elements on the diagonals). In [15], H.J. Landau established that, for any given spectra, there exists a symmetric Toeplitz matrix with those eigenvalues. Towards a practical algorithm, Chu [5] instead proved that fixed points of the isospectral flow with B(W ) as above are symmetric Toeplitz matrices, provided the eigenvalues are distinct. Chu's flow is particularly interesting from a numerical point of view because there exist periodic orbits. Thus, the flow does not always converge to a fixed point. However, the periodic orbits are unstable and because of the floating point drift in numerical methods Chu's flow in practice always converge to a symmetric Toeplitz matrix when the starting point has distinct eigenvalues [32].
A qualitatively better simulation of Chu's flow, which preserves the periodic orbits, can be obtained by restriction to centrosymmetric matrices. A matrix is said to be centrosymmetric if it is invariant with respect to a rotation of the components of π grade. In other words a matrix A is centrosymmetric if The set of centrosymmetric matrices of dimension n is a Lie algebra which we denote by Centro(n). In particular we have that the symmetric Toeplitz matrices are centrosymmetric and B(W ) is centrosymmetric when W is symmetric and centrosymmetric [32]. Therefore, for the Toeplitz problem, Chu's flow can be restricted to the symmetric-centrosymmetric matrices. With this restriction the periodic orbits are numerically preserved, and therefore the simulation of the flow is more realistic. In fact, in order to avoid the drifting out from the periodic orbits, it is necessary to respect the centrosymmetric symmetry of the original flow in the discrete approximation.
If S denotes the linear subspace of symmetric-centrosymmetric matrices, then by Theorem 3 the isospectral symplectic Runge-Kutta methods descend to an isospectral integrator on the symmetric-centrosymmetric matrices, provided that B(W ) is in the normalizer of S, which is so(n) ∩ Centro(n).
We use the midpoint based numerical scheme of Theorem 4 for Chu's flow with n = 4 and stepsize h = 0.1. The initial conditions, proposed in [32], are where N and W are n × n self-adjoint complex matrices. In [4], Brockett shows that for a diagonal N with distinct entries and W 0 a self-adjoint matrix with distinct eigenvalues, W (t) converges exponentially fast to a diagonal matrix with the eigenvalues sorted accordingly to the order of the entries of N . There are interesting connections between the Brockett flow and information theory. Indeed, the Brockett flow can be viewed as a gradient flow, with respect to the Fisher-Rao information metric, of a relative entropy functional on the statistical manifold of multivariate Gaussian distributions [23,Sec. 3.4.3]. We apply the isospectral midpoint method with h = 0.1. In Figure 9 we plot the eigenvalues and the components variation for a randomly generated self-adjoint initial matrix W 0 of dimension 3 × 3 and N = diag(1, 2, 3). Figure 9 displays the exponential convergence to a similar diagonal matrix.

Evolution of components
Evolution of error in Casimirs (eigenvalues)