Analytic and numerical aspects of isospectral flows Amandeep Kaur Wolfson College DAMTP, Center for Mathematical Sciences University of Cambridge April 2016 Abstract In this thesis we address the analytic and numerical aspects of isospectral flows. Such flows occur in mathematical physics and numerical linear algebra. Their main structural feature is to retain the eigenvalues in the solution space. We explore the solution of Isospectral flows and their stochastic counterpart using explicit generalisation of Magnus expansion. In the first part of the thesis we expand the solution of Bloch–Iserles equa- tions, the matrix ordinary differential system of the form X ′ = [N,X2], t ≥ 0, X(0) = X0 ∈ Sym(n), N ∈ so(n), where Sym(n) denotes the space of real n× n symmetric matrices and so(n) denotes the Lie algebra of real n× n skew-symmetric matrices. This system is endowed with Poisson structure and is integrable. Various important properties of the flow are discussed. The flow is solved using explicit Magnus expansion and the terms of expansion are represented as binary rooted trees deducing an explicit formalism to construct the trees recursively. Unlike classical numerical methods, e.g. Runge–Kutta and multistep methods, Magnus expansion respects the isospectrality of the system, and the shorthand of binary rooted trees reduces the computational cost of the exponentially growing terms. The desired structure of the solution (also with large time steps) has been displayed. Having seen the promising results in the first part of the thesis, the technique has been extended to the generalised double bracket flow X ′ = [[N,X] + M,X], t ≥ 0, X(0) = X0 ∈ Sym(n), where N ∈ diag(n) and M ∈ so(n), which is also a form of an Isospectral flow. In the second part of the thesis we define the generalised double bracket flow and discuss its dynamics. It is noted that N = 0 reduces it to an integrable flow, while for M = 0 it results in a gradient flow. We analyse the flow for various non-zero values of N and M by assigning different weights and observe Hopf bifurcation in the system. The discretisation is done using Magnus series and the expansion terms have been portrayed using binary rooted trees. Although this matrix system appears more complex and leads to the tri-colour leaves; it has been possible to formulate the explicit recursive rule. The desired structure of the solution is obtained that leaves the eigenvalues invariant in the solution space. Keywords: Isospectral flow, ordinary differential equations, eigenvalues, Mag- nus expansion, Lie group, Lie algebra, binary trees, Cayley transform, double bracket flow, Toda lattice, QR algorithm, Runge–Kutta. Dedicated to my late grandmother, Veeran Wali. Acknowledgements I never thought that writing the acknowledgements section in my thesis would be more difficult than any other section. I express the deepest appreciation to each and everyone who helped me in these four years at Cambridge. First of all, I would like to express my gratitude to my supervisor Arieh Iserles for giving me the opportunity to work in DAMTP and to pursue research in Mathematics in such a beautiful University. He truly is a great mathematician and his suggestions have always been valuable to me. Also, my adviser Carola- Bibiane Scho¨nlieb has always been a great support. I thank her for motivating me from time to time. My officemates and friends Evangelos Papoutsellis and Pranav Singh created a beautiful working environment. I have learnt a lot having both of them as my friends. Evangelos and I had an exchange program for teaching me Greek in exchange of every Hindi(or Punjabi) word he learnt! I thank all the academic and supporting staff at the department for being helpful in various ways. I acknowledge all the visitors during my stay at Cam- bridge; professors, postdocs, and students, the discussions with them were highly knowledgeable and fruitful. I sincerely thank Wolfson college for accepting me as a graduate student and for providing me with accommodation for four years. The porters and all other staff have been very helpful to make this a home away from home. I appreciate my friends Zab, Varun, Rav, Angad, Deepyanti, Liza, Adarsh for all the fun and support, and the residents of Toda house for being so understanding. A special word of gratitude is due to Nuri Purswani for all her support. She has been a great friend and stood with me in all the tragic moments I had to go through during my time at Cambridge. I also thank my friends Parminder Kaur and Ankur Dutt for their support, specially in the writing period. I am indebted to Mugdha who allowed me to stay at her place. I acknowledge Shiromani Gurdwara Prabandhak Committee–Cambridge Com- monwealth Trust (SGPC–CCT) for the financial support. Also, I am grateful to my friend Nuri my brother Davinder who contributed to the funding in the final year of my PhD. I dedicate this thesis to my beloved grandmother who died just before I started my journey to Cambridge. I am indebted to my cute nephews Jaiveen, Ganeev and Fateh for making me smile all the time. I thank Shavinder and Nishu (my sisters-in-law) for their love and support. This thesis would have been impossible without the support of my loving par- ents Jaspal Singh and Kulwant Kaur and my brothers Davinder and Kuljinder. Their unconditional love, patience and support in every aspect help me to move forward. I thank my brothers for being there at every single step helping me cross all the hurdles. A special thanks to my parents for nurturing me. They have always been a source of encouragement and inspiration. And I am greatful to someone special, my friend, guide and protector, for always being with me... So Satgur pyara mere naal hai !!! 8 Statement of Originality This dissertation is the result of my own work and includes nothing which is the outcome of work done in collaboration except where specifically indi- cated in the text. I further state that no substantial part of my dissertation has already been submitted, or, is being concurrently submitted for any such degree, diploma or other qualification at the University of Cambridge or any other University or similar institution except as declared in the Preface and specified in the text Contents 1 Introduction 13 2 Mathematical Preliminaries 17 2.1 Isospectral Flows and their applications . . . . . . . . . . . . . . . . . . . . 17 2.2 Lie groups and Lie algebra . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.3 Some examples of Lie groups and algebras . . . . . . . . . . . . . . . . . . . 26 2.4 Magnus expansion and binary rooted trees . . . . . . . . . . . . . . . . . . . 27 3 Bloch–Iserles equations 33 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.2 An expansion of the solution . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.3 Representation by binary trees . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.4 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4 The dynamic generalised double bracket flow 51 4.1 Double bracket flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2 Analysing the generalised double bracket flow . . . . . . . . . . . . . . . . . 53 4.2.1 Analysing taking 2× 2 matrices . . . . . . . . . . . . . . . . . . . . . 53 4.2.2 An explicit solution taking 2× 2 matrices . . . . . . . . . . . . . . . 56 4.3 Some interesting phase portraits . . . . . . . . . . . . . . . . . . . . . . . . 59 4.4 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 5 Discretisation of generalised double bracket flow 67 5.1 Expansion of the solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.2 Representation by binary trees . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.2.1 Constructing the elements (5.7) and (5.10) . . . . . . . . . . . . . . 74 5.2.2 Growing trees using the algorithm . . . . . . . . . . . . . . . . . . . 79 6 Conclusion 91 Bibliography 93 11 CONTENTS 12 Chapter 1 Introduction Isospectral flows are matrix systems of ordinary differential equations of the form X ′ = [B(X), X], t ≥ 0, X(0) = X0 ∈ Sym(n), (1.1) where B(X) : Sym(n) → so(n). Their main structural feature is that they preserve the eigenvalues of the solution matrix. Isospectral flows occur in many important applications. First and the best known example is the Toda lattice, a one-dimensional lattice of particles whose motion is described by a nearest-neighbour interaction of an exponential type. It can be used to model a wide range of particle systems, ranging from the hard-sphere limit to the atomic case [Tod81, Tod89]. Another important example is the QR flow. The QR method for finding the eigenvalues of a matrix can be executed as an isospectral flow at unit intervals. QR flow is the generalization of non-periodic Toda flow. Such flows were first investigated by Symes [Sym82] and subsequently in [Nan82, DNT83, DRTW91, Chu84, CD89, Lag91] and elsewhere. Other well known examples include eigenvalue problems and inverse eigenvalue problems for symmetric Toeplitz matrices [FNO87]. Note that if we let B(X) = [N,X], in (1.1) where N ∈ Sym(n) then it leads to the double-bracket flows. Double bracket flows are isospectral flows given by the equations X ′ = [[N,X], X], t ≥ 0, X(0) = X0 ∈ Sym(n) (1.2) where N ∈ Sym(n). They were introduced by Brockett [Bro91] and Chu and Driessel [CD90]. Double bracket flows were discretised and then solved by Iserles by the method of Magnus series [Ise02] and were generalized for more parameters [BI05]. Also, methods based on Magnus expansion are proposed for the numerical integration of the double- bracket flow and a bound on the convergence domain is provided by Casas [Cas04]. It is obvious that we can discretise isospectral flows by traditional numerical methods (e.g. Runge–Kutta and multistep), but, once n ≥ 3, these methods cannot respect the 13 Introduction isospectrality of the system, i.e. the numerical solution changes the eigenvalues [CIZ99]. Isospectrality is essential for applications ranging from classical mechanics, like Toda flows and N -body systems, to linear algebra, like QR flows and inverse eigenvalue problems, so we need to solve systems of the form (1.1) by a method that respects it [CIZ97, CIZ99, Zan98]. In this thesis we solve the given isospectral flows (namely, Bloch–Iserles equations (3.1) and generalised double bracket flow (4.1)) using the method of Magnus series. For the solution, first we convert the isospectral flow to a Lie-group flow and then translate it into a Lie-algebraic equation. This method preserves the isospectrality and gives the desired structure of the solution with large time steps. We show that the solution of an isospectral flow can be represented in the form X(t) = eΩ(t)X0e −Ω(t), where instead of computing X at the first place, we obtain the Taylor expansion of Ω. Note that this ensures automatically that the numerical solution, being similar to X0, is isospectral. We will see that the Taylor expansion of Ω can be formed algorithmically from the elements of the flow and linear combinations of their commutators and(or) anti-commutators. Our goal is to determine the rules for finding the terms of Ω to an arbitrary accuracy. The terms are represented by binary rooted trees and an algorithm is formed to construct the next tree by recursion and to calculate the coefficient of each tree. This lays the foundations to a more general setting, namely the explicit representation of the solution when B(X) can be represented in a finite “alphabet”. The representation as binary trees is very important because otherwise, as the number of terms in each iteration grows exponentially, the complexity of manual computation becomes prohibitive. By indexing the terms in the expansion with a subset of binary trees, it is convenient to derive explicit recurrence relations. Also, it is remarkable that the skew-symmetry and Jacobi identity obeyed by the commutator help us to reduce the number of the terms by cancelling or writing certain terms as the linear combination of other terms. We organise the thesis as follows. In chapter 2, we discuss the theory of isospectral flows and their applications ranging from classical mechanics to linear algebra. We also discuss the classical numerical methods, isospectral methods, and give introduction to Lie group and Lie algebra along with illustrating some examples. Towards the end of the chapter, we explain Magnus expansion which is used in the later chapters for discretisation, and define some basic concepts of graph theory for constructing the binary rooted trees for the expansion terms. Chapter 3 consists of expansion of the solution of Bloch–Iserles equations, a matrix ordinary differential system of the form X ′ = [N,X2], t ≥ 0, X(0) = X0 ∈ Sym(n), N ∈ so(n), where Sym(n) denotes the space of real n×n symmetric matrices and so(n) denotes the Lie algebra of real n × n skew-symmetric matrices. We discuss the motivational properties of the system. This system is an isospectral flow, and it has been seen that the system is endowed with Poisson structure. We represent the solution of isospectral flow 14 in the form X(t) = eΩ(t)X0e −Ω(t), where instead of computing X at the first place, we obtain the Taylor expansion of Ω. The numerical solution being similar to X0, ensures the isospectrality. The solution is computed using Magnus series and the Taylor expansion is formed from the elements X0 and N , and the linear combination of their commutators and anti-commutators. We also introduce curly bracket, {X0} = NX0 +X0N , that means the matrix N features implicitly in the bracket. Thus, each term of the expansion is made of commutators and curly brackets with some expression X0. This representation further contributes to the shorthand of binary rooted trees by preventing bicolour leaves. It has been shown graphically that the solution with the Magnus series preserves the spectrum of the solution space. The terms of the expansion are represented using binary rooted trees that helps to reduce the manual computation which has been growing exponentially with each iteration. Moreover, we deduce an explicit formalism to construct the trees and their coefficients recursively. Results are analysed by plotting the error graphs of the solution and the graphs for error in eigenvalues, these plots are generated in MATLAB. In chapter 4, we introduce the generalised double bracket flows, the matrix system of ordinary differential equations of the form X ′ = [[N,X] + M,X], t ≥ 0, X(0) = X0 ∈ Sym(n) where N ∈ diag(n) and M ∈ so(n). Analysis is done for 2 × 2 and 3 × 3 matrices and the dynamics are represented as phase portraits. We see that the phase plots, computed with MATLAB result in spirals. It is observed that the system undergoes a Hopf bifurcation that gives birth to the limit cycles. We further discuss the special cases where we see that the different weightage of N and M influence the result. We also discuss the case with identical diagonal elements of the matrix N . We perform these experiments taking 3× 3 matrices. In chapter 5, we extend the method of Magnus series, to generalised double bracket flow. This system seems to be more challenging to be discretised and for the solution to be deduced as an explicit representation. Due to the fact that it is going to have three matrices appearing implicitly in the Taylor expansion together with commutators. We obtain the Taylor expansion and represent the terms as binary rooted trees. Although this matrix system appears more complex and leads to the tri-colour leaves; it has been possible to formulate the explicit recursive rule. In the case of this matrix system, we have two types of trees and an interplay between them, which makes it more complicated. A step by step algorithm is developed and recurrence formula is defined that enables us to compute the trees and coefficients explicitly. We compare the solution and error in eigenvalues against ode45 method. In the concluding chapter 6, we briefly summarise the work presented in this dissertation. 15 Introduction 16 Chapter 2 Mathematical Preliminaries In this chapter we assemble and briefly explain, the mathematical objects and the keywords that form a basis for the work presented in this dissertation. We start with explaining the mathematical term, isospectral flow in section 2.1 along with its applications and discuss about the classical methods and isospectral methods. This is followed by the introduction of Lie groups, Lie algebras, and their examples in section 2.2 and section 2.3 respectively. Later in section 2.4 we explain the Magnus expansion and the basic definitions of graph theory for the formation of binary rooted trees. Main references to this chapter are [IMKNZ99] [Ise99] [Zan98]. 2.1 Isospectral Flows and their applications Isospectral flows are the matrix system of ordinary differential equations of the form X ′ = [B(t,X), X], t ≥ 0 X(0) = X0, (2.1) where X,B(t,X) ∈ Vn×n, set of n× n matrices with the entries in V = R or C. The functions B and X which satisfy the matrix differential equation (2.1) are called a Lax pair [Lax68] and [A,B] = AB−BA is the commutator or Lie bracket. It is easy to verify that X(t) = Q(t)X0Q(t) −1, t ≥ 0, (2.2) is the solution of (2.1), where Q(t) is the solution of Q ′ = B(t, QX0Q −1)Q, t ≥ 0, Q(0) = I. (2.3) To verify this isospectral deformation: We begin with associating the differential equation Q ′ = B(t,X)Q, t ≥ 0, Q(0) = I, (2.4) 17 Mathematical Preliminaries where I is the identity matrix, to (2.1). Let us construct the matrix functionQ(t)−1X(t)Q(t), we observe that, from (2.1) and (2.4), together with d dt Q−1 = −Q−1Q′Q−1, (2.5) we get d dt Q(t)−1X(t)Q(t) = O (2.6) where O is the zero matrix. Thus, we can say that Q(t)−1X(t)Q(t) is time independent. Therefore, it must equal its initial condition substituting t = 0. Hence, X(t) = Q(t)X0Q(t) −1 (2.7) is a similarity transformation that leaves the eigenvalues same. In the following chapters, while discussing about Bloch–Iserles equations and double bracket flows, we shall see that there are many important examples where B is skew-symmetric. Then the equation (2.4) becomes an orthogonal flow. In the case of orthogonal flow, Q−1 = QT and if X0 is symmetric, X(t) also remains symmetric. Isospectral flows occur in various applications. First and a well known example is the Toda lattice. Toda lattices, introduced by Morikazu Toda (1967) are one-dimensional lattices of par- ticles whose motion is described by a nearest-neighbour interaction of an exponential type and can be used to model a continuum of flows, ranging from the hard-sphere limit to the atomic case [Tod81, Tod89]. Consider the equations of motion m d2xi dt2 = d dr ψ(xi+1 − xi)− d dr ψ(xi − xi−1), i = 1, 2, ..., (2.8) where m is the mass of all the particles, xi is the displacement of the ith particle and ψ(r) is the interaction potential. Introduce the momenta pi = mx ′ i and generalised coordinates qi = xi, i = 1, 2, ..., and let ψ be the exponential potential ψ(r) = e−r + r, (2.9) such that ψ ′ (r) = −e−r + 1, the equations of motion become the Hamiltonian system q ′ i = 1 m pi, p ′ i = e −(qi−qi−1) − e−(qi+1−qi) (2.10) 18 2.1. Isospectral Flows and their applications assume m = 1 the corresponding Hamiltonian function is H(p,q) = 1 2 ∑ i p2i + ∑ (e−(qi+1−qi) − 1), (2.11) where p = {pi} and q = {qi}. To change this Hamiltonian system into isospectral form we introduce the change of variables as done by [Tod81] [Fla74] αi = 1 2 e−(qi+1−qi)/2, (2.12) βi = pi, (2.13) so that α ′ i = αi(βi − βi+1), (2.14) β ′ i = 2(α 2 i−1 − α2i ), (2.15) we get the Lax form, X ′ = [B,X], where X =  β1 α1 0 ... 0 α1 β1 α2 . . . ... 0 α2 . . . . . . 0 ... . . . . . . . . . αd−1 0 . . . 0 αd−1 βd  , and B =  0 −α1 0 ... 0 α1 0 −α2 . . . ... 0 α2 . . . . . . 0 ... . . . . . . . . . −αd−1 0 . . . 0 αd−1 0  . (2.16) Another important example is the QR flow. QR flow is the generalization of non- periodic Toda flow. The connection between Toda flow and QR algorithm was observed by Symes in [Sym82]. Let f is an analytic function in an open domain containing the 19 Mathematical Preliminaries spectrum of X0. Then we can define B(X) = fl(X)− fu(X), (2.17) where the subscripts l and u denote the lower and upper triangular parts of f(X), re- spectively. Let us assume X is symmetric and QΛQT is its factorization, here Q is an orthogonal matrix and Λ is a diagonal matrix. So we can say that f(X) = Qf(Λ)QT is also symmetric, and hence, B(X) is skew symmetric. QR flow is originated from the fact that for a symmetric positive definite initial condition X0 there exists a connection between the isospectral flow (2.1) and the iterates of the familiar QR method for finding the eigenvalues of a matrix. QR method consists in finding the decomposition Xk = QkUk (2.18) where Qk is an orthogonal matrix and Uk is an upper triangular matrix, with positive diagonal entries (Golub and van Loan 1989), k = 1, 2, · · ·n. The next matrix in the sequence is then obtained as Xk+1 = UkQk (2.19) and, when X0 is symmetric and positive definite, the method converges at an exponential speed to a diagonal matrix whose entries are the eigenvalues of X0. It is easy to verify that Xk+1 = Q T kXkQk, (2.20) rendering the connection with the isospectral flow is less obscure. If X0 is symmetric and positive definite, choosing f(x) = log x, the matrices X(1), X(2), ..., interpolate the iterates X1, X2, ..., obtained this time with QR method for eigenvalues. As we said earlier this connection was discovered by [Sym82] and then later investegated by [Nan82, DNT83, DRTW91, Chu84, CD89, Lag91] and many others. Other well known examples are Eigenvalue problems and inverse eigenvalue problems for symmetric Toeplitz matrices [FNO87]. There are existing numerical methods that are most commonly used to solve ODEs. Clearly, it is possible that we discretise the isospectral flows by traditional numerical meth- ods (e.g. Runge–Kutta and multistep), but these methods cannot respect the isospectrality of the system; once n ≥ 3, the numerical solution changes the eigenvalues [CIZ99]. We need to solve these flows by a method that respects it, since isospectrality is essential for applications ranging from classical mechanics, like Toda flows and N -body systems, to linear algebra, like QR flows and inverse eigenvalue problems [CIZ97, CIZ99, Zan98]. Also, the detailed analysis done in [Zan98] shows the negative results for classical ODE methods. Now, while we talk that the behaviour of the classical ODE methods in the preservation 20 2.2. Lie groups and Lie algebra of conservation laws is generally poor, the question that comes to our mind is: Is it possible to devise numerical methods that are isospectral? The answer is positive. Consider again the lax pair (2.1) X ′ = [B(t,X), X], t ≥ 0, with X0, a symmetric initial condition and B(X), a skew-symmetric matrix function. If we solve the system (2.1) with standard ODE method, we are bound to lose isospectrality. Recalling the proof of isospectrality of the system (2.1) from above, we see that Q(t) is the solution of Q ′ = B(t,X)Q, t ≥ 0, Q(0) = I, which is an orthogonal flow because B(X) is skew-symmetric and X(t) can be obtained by similarity transformation, X(t) = Q(t)X0Q(t) T . Hence, the main idea of the isospectral methods is to solve the orthogonal flow instead and then approximate X(t). As we know that the majority of numerical methods for ODEs do not preserve quadratic conservation laws, an easy way to construct orthogonal methods is by means of projection methods [LD94], evaluation of QR factorization can be seen in [GVL89] and [Dem90]. There are other examples of Runge–Kutta schemes [CIZ97][CIZ97] and [DRVV94], Cayley transforms [GVL89] and [DVV99], Modified Gauss– Legendre schemes, Semi-explicit methods [CIZ97] and the Lie group methods; Magnus expansion and Fer expansion. We discuss Lie group and Lie algebra in more detail in the following section. 2.2 Lie groups and Lie algebra In this section we assemble the elements of Lie groups and Lie algebras. We briefly state background theory and introduce differentiable manifolds, Lie groups and their properties referring mainly to [IMKNZ99]. Lie groups and Lie algebras have originated by Sophus Lie (1842-1899) while solving dif- ferential equations by quadrature, using symmetry methods. In the twentieth century an abstract view of Lie group theory emerged. This formulation simplifies mathematical anal- ysis and has become popular in Mathematics. However, the abstract theory concentrates on understanding mathematical structures rather than exposing applications in solving differential equations. Therefore, it is not clearly known to most applied mathematicians that Lie groups are really very useful in applied and computational mathematics. Tra- ditionally, numerical integration of ordinary differential equations(ODEs) is basically the concept of solving the initial value problems of the form x ′ = f(t, x), t ≥ 0, x(0) = x0, x(t) ∈ Rn (2.21) 21 Mathematical Preliminaries where f is vector field on Rn×Rn. Traditional methods for numerical integration, like Euler forward method, Runge-Kutta and multistep methods iterate step by step by adding a vector on each iteration, xn+1 = xn + han, (2.22) where an is calculated by the numerical method which we use and h is the stepsize. We can say that classical integrators are formulated using a set of basic notations given by translations on Rn to advance the numerical solution. A major motivation for Lie group methods is the possibility of replacing the domain Rn with more general configuration spaces and replacing translation on Rn by more general families of basic motions on the domain. To use Lie group methods and Lie algebra we should have an abstract view of differential manifolds, which is the domain on which differential equations evolve. Definition 2.2.1. A d-dimensional manifold M is a d-dimensional smooth surface M⊂ Rn for some n ≥ d. Definition 2.2.2. Let M be a d-dimensional manifold and suppose that ρ(t) ∈ M is a smooth curve such that ρ(0) = p. A tangent vector at p is defined as a = dρ(t) dt |t=0 . The set of all tangents at p is called the tangent space at p and denoted by TM|p. It has the structure of a d-dimensional linear space: if a, b ∈ TM|p then a + b ∈ TM|p and αa ∈ TM|p for any real α. The collection of all tangent spaces at all points p ∈ M is called the tangent bundle of M and denoted by TM = ⋃p∈M TM|p. Definition 2.2.3. A (tangent) vector field on M is a smooth function F : M → TM such that F (p) ∈ TM|p for all p ∈M. The collection of all vector fields onM is denoted by X(M). Addition and scalar multiplication of vector fields are defined pointwise in a natural way as (F + G)(p) = F (p) + G(p) and (αF )(p) = α(F (p)). If F,G ∈ X(M) then also F +G ∈ X(M) and αF ∈ X(M) for all real α. Definition 2.2.4. Let F be a tangent vector field on M. By a differential equation (evolving) on M we mean a differential equation of the form y′ = F (y), t ≥ 0, y(0) ∈M, (2.23) where F ∈ X(M). Whenever convenient, we allow F in (2.23) to be a function of time, F = F (t,y). The flow of F is the solution operator Ψt,F : M→M such that y(t) = Ψt,F (y(0)) 22 2.2. Lie groups and Lie algebra solves (2.23). Lemma 2.2.5. Given two vector fields F,G on Rn, the commutator H = [F,G] can be computed component wise at a given point y ∈ Rn as Hi(y) = n∑ j=1 { Gj(y) ∂Fi(y) ∂yj − Fj(y)∂Gi(y) ∂yj } . (2.24) Lemma 2.2.6. If F,G ∈ X(M) then H = [F,G] ∈ X(M). Lemma 2.2.7. Two flows Ψs,F and Ψt,G commute if and only if [F,G] = 0. From (2.24) one may prove the following important features of the commutator which should be familiar in the special case (which we will encounter soon again) of a commutator of two matrices. Definition 2.2.8. A Lie algebra is a linear space V equipped with a Lie bracket, a bilinear, skew-symmetric mapping [·, ·] : V × V → g that obeys identities ((2.25)− (2.28)) from Lemma 2.2.9. Lemma 2.2.9. The commutator of vector fields satisfies the identities [F,G] = −[G,F ], (skew symmetry), (2.25) [αF,G] = α[F,G], for α ∈ R, (2.26) [F +G,H] = [F,H] + [G,H], (bilinearity), (2.27) 0 = [F, [G,H]] + [G, [H,F ]] + [H, [F,G]], (Jacobi’s identity). (2.28) Definition 2.2.10. A Lie algebra of vector fields is a collection of vector fields which is closed under linear combination and commutation. In other words, letting g denote the Lie algebra, B ∈ g ⇒ αB ∈ g for all α ∈ R. B1, B2 ∈ g ⇒ B1 +B2, [B1, B2] ∈ g. Given a collection of vector fields B = {B1, B2, . . .}, the least Lie algebra of vector fields containing B is called the Lie algebra generated by B. 23 Mathematical Preliminaries Definition 2.2.11. A Lie algebra homomorphism is a linear map between two Lie alge- bras, ϕ : g× h, satisfying the identity ϕ([v, w]g) = [ϕ(v), ϕ(w)]h, v, w ∈ g. An invertible homomorphism is called an isomorphism. Definition 2.2.12. A Lie group is a differential manifold G equipped with a product · : G × G → G satisfying p · (q · r) = (p · q) · r ∀p, q, r ∈ G (associativity), ∃I ∈ G such that I · p = p · I = p ∀p ∈ G (identity element), ∀p ∈ G ∃p−1 ∈ G such that p−1 · p = I (inverse), The maps (p, r) 7→ p · r and p 7→ p−1 are smooth functions (smoothness). Definition 2.2.13. An action of a Lie group G on a manifold M is a smooth map Λ : G ×M→M satisfying Λ(I,y) = y ∀y ∈M. Λ(p,Λ(r,y)) = Λ(p · r,y) ∀p, r ∈ G, ∀y ∈M. (2.29) If this relation does hold only in a local sense, for all elements p and r sufficiently close to the identity I ∈ G, we say that Λ is local action. Definition 2.2.14. Let G be a Lie group and g its Lie algebra. The exponential mapping exp : g→ G is defined as exp(a) = σ(1) where σ(t) ∈ G satisfies the differential equation σ′(t) = aσ(t), σ(0) = I. Definition 2.2.15. Let p ∈ G and let σ(t) be a smooth curve on G such that σ(0) = I and σ′(0) = b ∈ g. The adjoint representation is defined as Adp(b) = d dt pσ(t)p−1|t=0. (2.30) The derivative of Ad with respect to the first argument is denoted ad. Let ρ(s) be a smooth curve on G such that ρ(0) = I and ρ′(0) = a. Now we know that: Ada(b) = d ds Adρ(s)(b)|s=0 = [a, b]. (2.31) The following formulae show that Ad is both a linear group action (of G on g) and also 24 2.2. Lie groups and Lie algebra that for a fixed argument p it is a Lie-algebra isomorphism of g onto itself: Adp(a) ∈ g for all p ∈ G, a ∈ g (2.32) Adp ◦Adq = Adpq (2.33) Adp(a+ b) = Adp(a) + Adp(b) (2.34) Adp([a, b]) = [Adp(a),Adp(b)]. (2.35) Note that according to (2.34) both Adp and ada are linear in their second argument, hence they may be regarded as matrices acting on the linear space g. This gives meaning to the following important formula relating Ad, ad and the exponential mapping: Adexp(a) = expm(ada). (2.36) Definition 2.2.16. A real matrix Lie group is a smooth subset G ⊆ Rn×n, closed under matrix products and matrix inversion. We let I ∈ G denote the identity matrix. Definition 2.2.17. The Lie algebra g of a matrix Lie group G is the linear subspace g ⊆ Rn×n consisting of all matrices of the form g = { A ∈ Rn×n : A = dρ(s) ds |s=0 } , where ρ(s) ∈ G is a smooth curve such that ρ(0) = I. The space g is closed under matrix additions, scalar multiplication and the matrix commutator [A,B] = AB −BA. (2.37) Complex matrix Lie groups and algebras are defined similarly. Definition 2.2.18. A differential equation on a matrix Lie group is an equation of the form Y ′ = A(t, Y )Y, t ≥ 0, Y (0) ∈ G, (2.38) where A : R× G → g and AY is the usual matrix product between A ∈ g and Y ∈ G. This can be verified that this is the special case of the general form of a differential equation on a manifold, where M = G is a matrix Lie group and the action Λ is taken to be the left (matrix) multiplication in G, Λ(R, Y ) = RY. We find λ∗(A)Y = AY. Since g is defined as the collection of all tangent directions at I ∈ G and matrix multipli- 25 Mathematical Preliminaries cation by Y is an invertible mapping, we see that any tangent at Y can be written in the form AY and all differential equations on G can be written in the form (2.38). Definition 2.2.19. The exponential mapping expm : g→ G is defined as expm(A) = ∞∑ j=0 Aj j! . (2.39) The adjoint representation, Ad, and its derivative, ad, are defined as AdP (A) = PAP −1 (2.40) adA(B) = AB −BA = [A,B]. (2.41) 2.3 Some examples of Lie groups and algebras We introduce here some concrete examples of Lie groups and algebras. In each case it is easy to verify that all the axioms of a group or an algebra, as the case might be, are fulfilled. • The set of all real n × n nonsingular matrices is a (multiplicative) Lie group, the general linear group GL(n). The corresponding Lie algebra is the set Rn×n of all n× n real matrices which, we denote by gl(N). The general linear group and algebra can be defined over other fields than R in which case we communicate this in the second argument. For example, GL(n,C) consists of all nonsingular n× n complex matrices. • All members of GL(n) with unit determinant form the special linear group SL(n). Its Lie algebra, sl(N), consists of all matrices in gl(N) with zero trace. • All the matrices X ∈ SL(4) such that XJXT = J, where J =  1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 −1  form the Lorenz group SO(3, 1). Its Lie algebra so(3, 1) is made out of all F ∈ gl(4) such that FJ + JF T = O. • N ×N real orthogonal matrices form the orthogonal group O(n), whose Lie algebra so(N) consists of N ×N skew-symmetric matrices. 26 2.4. Magnus expansion and binary rooted trees The set SO(n) = SL(N)∩O(N), consisting of N ×N real orthogonal matrices with a unit determinant, is the special orthogonal group. Its Lie algebra is so(N). This is not contradictory; this is to be noted, we never claimed that two different Lie groups must have different Lie algebras. If G is a Lie group and GId its connected component such that I ∈ GId (precisely the situation with O(n) and SO(n), respectively) then they produce the same Lie algebra. • The set of all (2N)× (2N) real matrices X such that XJXT = J , where J = ON IN −IN ON  , is the symplectic group and is denoted by Sp(N). (The Jacobian of the flow of a Hamiltonian ODE system evolves in Sp(N).) The corresponding Lie algebra, sp(N), consists of F ∈ gl(2N) such that FJ + JF T = O. • As an example of complex Lie groups, we mention the unitary group U(N ;C) of all the N ×N complex unitary matrices: X ∈ U(N ;C) if and only if XXH = I. The Lie algebra corresponding to the U(N ;C) is the set u(N ;C) of all the skew-Hermitian matrices in gl(N ;C). • Similarly, for the case of O(n) and SO(n), we obtain the special unitary group in- tersecting U(N ;C) with SL(N ;C). Its Lie algebra, su(N ;C), is composed of N ×N complex skew-Hermitian and traceless matrices In this dissertation we use Magnus expansion, which is a Lie group expansion, to dis- cretise Bloch–Iserles equations and generalised double bracket flow. Let us have a brief look at the Magnus expansion in the following section. 2.4 Magnus expansion and binary rooted trees In this section we explain the Magnus expansion [Mag54] which we apply in the next chapters to discretise Bloch–Iserles equations and generalised double bracket flow. The main references for this section are [Mag54], [Ise99] and [IN99] other references if any, are mentioned explicitly. Magnus expansion provides an exponential representation of the solution of linear ordi- nary differential equation. Consider a matrix differential equation Y ′(t) = A(t)Y (t), t ≥ 0, Y (0) = Y0, (2.42) 27 Mathematical Preliminaries where A is a n× n matrix. For n = 1, solution of (2.42) is Y (t) = e ´ t 0 A(ξ)dξY0, t ≥ 0. (2.43) For n > 1 and if A is not a constant matrix then the above expression (2.43) is no longer the solution of the problem. Magnus [Mag54] proposed the solution to the matrix initial value problem and he ex- pressed the solution in the form Y (t) = eΩY0. (2.44) The expression for Ω is (by Felix Hausdorff [Hau06]) given by Ω′ = ∞∑ m=0 Bm m! admΩA, t ≥ 0, Ω(0) = 0, (2.45) where Bm,m ∈ Z are Bernoulli’s numbers and admΩ is an iterated commutator and we have ad0ΩA = A, ad 1 ΩA = [Ω, A], ad 2 ΩA = [Ω, [Ω, A]], ..., ad m ΩA = [Ω, ad m−1 Ω A], where [Ω, A] = ΩA−AΩ. Magnus observed in [Mag54] that Ω can be written as a linear combination of multiple integrals. Employing Picard’s iteration, Ω0(t) = 0, Ωs+1(t) = ˆ t 0 ∞∑ m=0 Bm m! admΩs(ξ)A(ξ). it is possible to show that Ω(t) = ˆ t 0 A(ξ)dξ − 1 2 ˆ t 0 ˆ ξ1 0 [A(ξ2), A(ξ1)]dξ2dξ1 + 1 12 ˆ t 0 ˆ ξ1 0 ˆ ξ1 0 [A(ξ3), [A(ξ2), A(ξ1)]]dξ3dξ2dξ1 + 1 4 ˆ t 0 ˆ ξ1 0 ˆ ξ1 0 [[A(ξ3), A(ξ2)], A(ξ1)]dξ3dξ2dξ1 + · · · . Calculating further terms get out of hand because it is growing exponentially with fourfold integrals and three nested commutators. Therefore, Iserles and Nørsett [IN99] proposed an alternative way to use binary rooted trees as a shorthand for expansion terms. Before going to the tree representation of the Magnus expansion, let us have a 28 2.4. Magnus expansion and binary rooted trees brief view of basic concepts of graph theory [Har69]. • Let V = {v1, v2, ..., vr} be a finite set of distinct vertices and E = V × V a set of edges. Then G = 〈V,E〉 is a graph. • The graph is said to be connected if there is a path between any two vertices. • It is a tree if exactly one path links every two vertices. • The ordered set {(vsl , vtl) : l = 1, 2, ..., r} of edges is a path from vi ∈ V to vj ∈ V, i 6= j, if s1 = i, tl = sl+1, l = 1, 2, ..., r − 1 and tr = j. • The pair T = (G,w), where G is a tree and w ∈ V is its root, is called a rooted tree. There exist a natural partial order on T : we say that vi ≺ vj if vi precedes vj in the unique path extending from the root w to vj . In that case vi is the ancestor of vj , while vj is the successor of vi. • If vi ≺ vj and there is no vk ∈ V such that vi ≺ vk ≺ vj , we say that vi is the parent of vj and vj is the child of vi. Childless vertices are called leaves. • If each vertex in a rooted tree has at most two children, T is called a binary tree. If each vertex has either exactly two children or is a leaf, T is said to be a strictly binary tree. Constructing the binary rooted trees for the expansion [Ise99], we commence by assign- ing to A a trivial tree t A, t@τ1 τ2 [Hτ1 , Hτ2 ], and tτ1 ˆ t 0 Hτ1(ξ)dξ. A is an integrable n× n matrix function, we define a map τ → Hτ from T, a subset of binary rooted trees into n× n matrix functions. Thus, Ω can be written as Ω = tt − 1 2 tt@ t tt + 1 12 tt@ t tt @t t t + 1 4 tt@ t tt@ t tt + · · · . 29 Mathematical Preliminaries Following the tree algorithm proposed in [Ise99], writing the Magnus expansion in the form Ω(t) = ∞∑ r=1 ∑ τ∈Tr α(τ)Hτ (t), where the trees τs are of the form τ = tt@ τ1 t@τ2 t@ τn t , (2.46) and the constant α(τ) is given by α(τ) = Bs s! s∏ i=1 α(τi), s ∈ N, where Bs is the sth Bernoulli number and the trees τ1, τ2, · · · are featured earlier in the expansion. We construct the next level trees, Ω = tt − 1 2 tt@ t tt + 1 12 tt@ t tt @t t t (2.47) + 1 4 tt@ t tt@ t tt − 1 8 tt@ t tt@ t tt@ t tt − 1 24 tt@ t tt @t t t@t t t − 1 24 tt@ t tt@ t tt @t t t · · · . The Magnus expansion absolutely converges for every t ≥ 0 such that ˆ t 0 A(ξ)dξ ≤ ˆ 2pi 0 dξ 4 + ξ[1− cot( ξ2)] ≈ 1.086868702. 30 2.4. Magnus expansion and binary rooted trees In the next chapter we discretise the Bloch–Iserles equations using Magnus expansion that helps preserving the isospectrality of the solution space and gives the desired structure underlying our approach. Moreover, the shorthand of binary rooted trees helps to reduce the computational cost. We see that the solution terms in the Taylor expansion using the Mangus series, grow exponentially and manual computation becomes very expensive. Our approach of using binary rooted trees makes it possible to formulate the explicit representation of the solution. 31 Mathematical Preliminaries 32 Chapter 3 Bloch–Iserles equations This chapter is based on my paper “On solving an isospectral flow” [Kau16]. In this chapter we expand the solution of the matrix ordinary differential system, namely Bloch– Iserles equations, of the form X ′ = [N,X2], t ≥ 0, X(0) = X0 ∈ Sym(n), N ∈ so(n), where Sym(n) denotes the space of real n×n symmetric matrices and so(n) denotes the Lie algebra of real n × n skew-symmetric matrices. The flow is solved using explicit Magnus expansion, which respects the isospectrality of the system. We represent the terms of expansion as binary rooted trees and deduce an explicit formalism to construct the trees recursively. 3.1 Introduction We are concerned with the discretization of the matrix differential equation X ′ = [N,X2], t ≥ 0, X(0) = X0 ∈ Sym(n), N ∈ so(n). (3.1) The system (3.1) is known as the Bloch–Iserles (BI) equations. It is isospectral (preserves the eigenvalues of X(t)), is endowed with a Poisson structure and is integrable as proved by Bloch and Iserles [BI06]. We discretise this system using a similar approach, for instance in [Ise02] and [Cas04]. However, solving BI is much more complicated since it contains X2 in the expression. The above system is of interest for a number of reasons. Firstly, we can easily verify that it can be written in the form X ′ = [N,X]X +X[N,X], t ≥ 0, X(0) = X0 ∈ Sym(n). (3.2) For N ∈ so(n) and X ∈ Sym(n), we have [N,X] ∈ Sym(n), therefore it is a special case 33 Bloch–Iserles equations of a congruent flow X ′ = A(X)X +XAT (X), t ≥ 0, X(0) = X0 ∈ Sym(n), (3.3) where A : Sym(n) → M(n), where M(n) is the set of real n × n matrices, is sufficiently smooth. It is easy to verify that X(t) = V (t)X0V T (t), where V ′ = A(V X0V T )V, V (0) = I. That means the solution is an outcome of the general linear group GL(n) acting on Sym(n) by congruence. That proves that the signature of X(t) is constant [HJ91]. Another interesting aspect of the given set of equations is that they are dual to the generalized rigid body equations M ′ = [Ω,M ], t ≥ 0, M(0) ∈ so(n), where M = ΩJ + JΩ, J ∈ Sym(n) therefore Ω ∈ so(n) [Man76]. Also, it is clear that (3.1) can be rewritten in the form X ′ = [XN +NX,X], t ≥ 0, X(0) = X0 ∈ Sym(n), N ∈ so(n). Since XN +NX ∈ so(n) for X ∈ Sym(n), N ∈ so(n), it follows that the system (3.1) is indeed isospectral. There are traditional numerical methods (e.g. Runge–Kutta and multistep) to discre- tise (3.1), but, as we discussed earlier, this is unfortunate that once n ≥ 3, these methods cannot respect the isospectrality of the system and hence we are bound to lose the struc- tural feature of the system, i.e. the numerical solution changes the eigenvalues [CIZ99]. Isospectrality is essential for applications ranging from classical mechanics, like Toda flows and N -body systems, to linear algebra, like QR flows and inverse eigenvalue problems, so we need to solve (3.1) by a method that respects it [CIZ97, CIZ99, Zan98]. In this chapter we solve the given isospectral flow using the method of Magnus series. We represent the solution of (3.1) in the form X(t) = eΩ(t)X0e −Ω(t), where instead of comput- ing X at the first place, we obtain the Taylor expansion of Ω and in each step an orthogonal matrix Q(t) = eΩ(t) is evaluated. Approximating X(t) by X(t) = eΩ(t)X0e −Ω(t), ensures automatically that the numerical solution, being similar to X0, is isospectral. We will see that the Taylor expansion of Ω can be formed algorithmically from X0 and N and linear combinations of their commutators and anti-commutators. Our purpose is to determine the rules for finding the terms of Ω to an arbitrary accuracy. For the solution, first we convert the isospectral flow to a Lie-group flow and then translate it into a Lie-algebraic equation. It is observed that this method preserves the isospectrality (Figure 3.2) and gives the desired structure of the solution with large time steps. In section 3.2 we solve the given system of differential equations using the Magnus expansion to obtain the Taylor expansion of Ω. We will see that the number of terms in each iteration grow exponentially therefore the complexity of manual computation becomes prohibitive. Finally, in section 34 3.2. An expansion of the solution 3.3 we represent the terms by binary rooted trees and form an algorithm to obtain the recursive formula that is used to construct the next generation trees and to calculate the coefficient of each tree. This technique is the inauguration to the general setting, i.e. the explicit representation of the solution of (3.1) when B(X) can be represented in a finite “alphabet”, namely X0 and N . Once we construct the trees, it is comparatively easier to translate them into terms containing X0 and N . Also, it is remarkable that the skew- symmetry and Jacobi identity obeyed by the commutator help us to reduce the number of the terms by cancelling or writing certain terms as the linear combination of other terms. 3.2 An expansion of the solution As stated above, the Bloch–Iserles system can be rewritten in the form X ′ = [B(X), X], t ≥ 0, X(0) = X0 ∈ Sym(n). with B(X) = NX + XN, where B(X) : Sym(n) → so(n). This system is seen to be isospectral and it is standard to verify that X(t) = Q(t)X0Q T (t), t ≥ 0, (3.4) where Q(t) ∈ SO(n) is the solution of Q′(t) = (Q(t)X0QT (t)N +NQ(t)X0QT (t))Q(t), Q(0) = I. (3.5) In a similar way as Magnus [Mag54] did for linear equations, our idea is to represent the solution of (3.5) in the form Q(t) = eΩ(t), where Ω′ = ∞∑ 0 Br r! adrΩ(e ΩX0e −ΩN +NeΩX0e−Ω), Ω(0) = 0. (3.6) Here Bm, m ∈ Z are Bernoulli numbers and adrΩ is an iterated commutator defined by ad0ΩA = A, ad 1 ΩA = [Ω, A], ad 2 ΩA = [Ω, [Ω, A]], ..., ad m ΩA = [Ω, ad m−1 Ω A], where [Ω, A] = ΩA−AΩ. 35 Bloch–Iserles equations Now, taking Ω(t) = ∑∞ m=0 Ωmt m gives Ω′(t) = ∞∑ m=0 (m+ 1)Ωm+1t m (3.7) and this implies ∞∑ m=0 (m+ 1)Ωm+1t m = ∞∑ r=0 Br r! adrΩ(t)(e Ω(t)X0e −Ω(t)N +NeΩ(t)X0e−Ω(t)). (3.8) Comparing coefficients of t0, t1, t2... we get the values of Ω1,Ω2,Ω3... as follows Ω1 = NX0 +X0N, Ω2 = 1 2 (N [Ω1, X0] + [Ω1, X0]N), Ω3 = 1 3 (N [Ω2, X0] + [Ω2, X0]N) + 1 6 (N [Ω1, [Ω1, X0]] + [Ω1, [Ω1, X0]]N) + 1 6 [Ω2,Ω1]. We denote NX + XN = {X}, thereby rewriting the above values of Ω1,Ω2,Ω3 . . . . in a more succinct manner as Ω1 = {X0}, (3.9) Ω2 = 1 2 {[Ω1, X0]}, Ω3 = 1 3 {[Ω2, X0]}+ 1 6 {[Ω1, [Ω1, X0]]}+ 1 6 [Ω2,Ω1], etc. Note that X ∈ Sym(n), N ∈ so(n) implies that {X} ∈ so(n). Introducing the curly bracket, i.e. featuring N implicitly in NX + XN = {X} is important mainly for two reasons. Firstly, this helps to understand the recurrence of the terms in the expansion, else XN +NX gets convoluted with other terms (either by getting cancelled or by adding up). Later in section 3.3, while defining the structure of binary rooted trees; it prevents from getting bicolored leaves, for instance, bicolored leaves in [Ise02]. Thus, Ω(t) = t{X0}+ 1 2 t2{[{X0}, X0]} 36 3.2. An expansion of the solution +t3( 1 6 {[{[{X0}, X0]}, X0]}+ 1 6 {[{X0}, [{X0}, X0]]} + 1 12 [{[{X0}, X0]}, {X0}]) + . . . . (3.10) In next term there are three nested commutators. These are getting increasingly com- plex with each iteration and it is clear that the number of such terms is growing expo- nentially. In this chapter we will represent these terms in terms of rooted trees, similar to the case of double-bracket flows in [Ise02]. This builds upon an idea of Iserles and Nørsett [IN99] to use binary rooted trees as a shorthand for expansion terms. 10−4 10−3 10−2 10−6 10−4 10−2 100 102 104 ∆t l 2 e r r o r Error for Ω[3] Expected error O((∆t)3) Figure 3.1: Global error on logarithmic scale across an interval [0,1] with different time steps, after truncating the Taylor expansion up to third order terms. Before we follow the procedure to simplify the above expansion, preliminary error graph of this method and error graph for eigenvalues of this method, as compared to the MAT- LAB ode45 solver with built-in parameters, are presented in Figure 3.1 and Figure 3.2 respectively. In Figure 3.1 we display error in the solution of (3.1) in the interval [0,1] for a range of different step sizes ∆t. The error plot is generated by truncating the expansion up to order three and is compared against the theoretically expected error of O((∆t)3). The experiments were performed on random 25× 25 matrices. Note that in Figure 3.1 we are interested in the case when ∆t→ 0 (asymptotic limit), which corresponds to the left part 37 Bloch–Iserles equations 10−3 10−2 10−1 100 101 10−15 10−10 10−5 100 t l ∞ er ro r in ei g en va lu es Solution with ode45 Magnus with ∆t ≈ 1/10 Magnus with ∆t ≈ 1/160 Figure 3.2: Error plot showing absolute error of eigenvalues of the two methods, Magnus expansion and ode45, on a logarithmic scale. of the figure. In Figure 3.2 we calculate absolute error of eigenvalues of the two methods on a logarithmic scale. It is clearly seen that our method preserves the correct eigenvalues to machine accuracy (we note here that the machine epsilon 10−16 corresponds to relative error; however, we are calculating the absolute error in the graph). Despite a large time- step in the Magnus method, the error in eigenvalues stays very close to machine precision while the solution obtained using ode45 quickly strays away in terms of eigenvalues as time increases. This favourable behavior of the Magnus method is to be expected from the principles underlying our approach. Also, we have computed numerically the solution of the system for random 3×3 matrices using Lie group method using Magnus expansion. In the Figure 3.3, the phase portraits (X1,2, Xk,l) are displayed for (k, l) = (1, 1), (1, 3), (2, 2), (2, 3), (3, 3), with random initial condition, here by Xk,l we mean the kl th element of the matrix X. It has been proved in [BI06] that the system is Lie-Poisson and in [BBI+09] and [LT06] that it is integrable. So, it is clear that the behaviour of the solution displays the regularity and the solution curve lies on invariant tori. This is an indication of integrable Lie-Poisson structure. Similar behaviour is obtained for other randomly calculated matrices as well. 38 3.2. An expansion of the solution −1.5 −1 −0.5 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 −1.5 −1 −0.5 0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 −1.5 −1 −0.5 0 0.5 1 1.5 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 −1.5 −1 −0.5 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 Figure 3.3: The phase portraits (X1,2, Xk,l) for (k, l) = (1, 1), (1, 3), (2, 2), (2, 3), (3, 3), respectively for (3.1), with a random initial condition. Here, by Xk,l we mean the kl th element of the matrix X. 39 Bloch–Iserles equations 3.3 Representation by binary trees Once we look at the expansion Ω, we note that each term is written in just two ‘letters’ X0 and N , hence belongs to a free structure generated by them: The matrix N features implicitly in the bracket {Z} = NZ + ZN . Thus we can say that each term is made of commutators and curly brackets with some expression Z, which itself has been formed from X0 and N . We attempt to find the expansion of the solution in Taylor series of the form Ω(t) = ∞∑ r=1 ∑ τ∈Tr α(τ)Hτ (3.11) where Tr is the set of all binary trees of power r : A tree τ is of power r ≥ 1 if r is the greatest integer such that Hτ = O(t r), Hτ is an expression constructed from X0 and N according to rules implicit in the structure of the tree τ which will be explained next, and α is a scalar coefficient. Using rooted trees as a shorthand for expansion terms is an approach introduced by [IN99], that leads to a framework that elucidates the structure of individual terms and their relationship. While constructing trees, we commence by assigning X0 to a single node, i.e. a trivial tree, t X0. We define a function τ → Hτ from T = ⋃∞ r=1 Tr, a subset of binary rooted trees into n×n matrix functions by letting H s= X0 and, by induction, t@τ1 τ2 [Hτ1 , Hτ2 ], and tτ1 {Hτ1}. where Hτ1 and Hτ2 are already constructed expansion terms. For example, tt@ t tt ⇒ tt@ {X0} X0 ⇒ t[{X0}, X0]⇒ {[{X0}, X0]} ⇒ {[Ω1, X0]}. In particular (3.10) can be written as 40 3.3. Representation by binary trees Ω = tt t + 1 2 tt@ t tt t2 + 1 6 tt@ t tt@ t tt t3 (3.12) + 1 6 tt@ t tt @t t t t3 + 1 12 @t t ttt@ t tt t3 + · · · . To explore the general rules underlying this correspondence between trees and expansion terms, let us have again a look at the equation ∞∑ m=0 (m+ 1)Ωm+1t m = ∞∑ r=0 Br r! adrΩ(t) ( eΩ(t)X0e −Ω(t)N +NeΩ(t)X0e−Ω(t) ) . We know that eΩX0e −Ω = AdΩX0 = eadΩX0 = ∞∑ n=0 1 n! adnΩX0 Thus ∞∑ m=0 (m+ 1)Ωm+1t m = ∞∑ r=0 Br r! adrΩ(t) ( N ∞∑ n=0 1 n! adnΩ(t)X0 + ∞∑ n=0 1 n! adnΩ(t)X0N ) = ∞∑ r=0 Br r! adrΩ(t) ∞∑ n=0 1 n! tad n Ω(t)X0 = B0 0! ∞∑ n=0 1 n! tad n Ω(t)X0 + B1 1! ∞∑ n=0 1 n! t@Ω(t) t adnΩ(t)X0 + B2 2! ∞∑ n=0 1 n! t@Ω(t) t@ Ω(t) tad n Ω(t)X0 + · · · . (3.13) where 41 Bloch–Iserles equations adnΩ(t)X0 = t@@ Ω(t) ad n−1 Ω(t)X0 . Clearly, in (3.13) we see that each tree here can be represented in the form Ts,n 3 τ = t@τ1 t@ τ2 t@τs t t@κ1 t@ κ2 t@κn t , (3.14) where s ∈ {0, 1, 2, 3, . . .}, n ∈ {0, 1, 2, 3, . . .}. Here, the trees τ1, τ2, . . . τs, κ1, κ2, . . . κn have been featured earlier in the expansion, where τi ∈ Tpi , and κj ∈ Tqj , i = 1, 2, ..., s, j = 1, 2..., n and p1 + p2 + ...+ ps + q1 + q2 + ...+ qn + 1 = r. For s = 0 (for example, in the first four trees in (3.12)), the structure of the trees becomes T0,n 3 τ = tt@ κ1 t@κ2 t@ κn t . (3.15) For n = 0 (for example, the fifth tree), this becomes Ts,0 3 τ = t@τ1 t@ τ2 t@τs t t . (3.16) It is also possible to deduce the explicit form of the constant α by substituting the value of Ω from (3.11) in the form of α(τ) and Hτ in (3.13) and simplifying it for α(τ). Set α ( tt) = 1. (3.17) Let τ ∈ Tr, r ∈ N and suppose that α(τi) and α(κj) are known for i = 1, 2, . . . , s, j = 1, 2, . . . , n. Then α(τ) = 1 r Bs s! 1 n! s∏ i=1 α(τi) n∏ j=i α(κj), s, n ∈ N, (3.18) where Bs is the sth Bernoulli number. 42 3.3. Representation by binary trees Now we have the general pattern for recursion. Suppose that T1,T2, . . . ,Tr−1 are known and also the coefficient α in these sets is known. To construct Tr we note that every τ is of the form (3.14) for some s ∈ {0, 1, 2, 3, . . . , r− 1} and n ∈ {0, 1, 2, 3, . . . , r− 1}. For every such s and n we consider all the partitions p1 + p2 + ...+ ps + q1 + q2 + ...+ qn + 1 = r. For every partition we construct the tree τ in (3.14) and use (3.18) to determine the coefficient α. The trees which correspond to zero terms would be eliminated. Moreover, some trees can be replaced by linear combinations of other trees. Let us start from T1 (i) T1: For s = 0, n = 0, we have τ11 = tt , α(τ11 ) = 11 · 1 · 10! = 1. (ii) T2: (1) For s = 0, n = 1, we can have only one possibility, κ1 = tt τ21 = tt@ t tt , α(τ21 ) = 1 2 · 1 · 1 1! = 1 2 ; (2) For s = 1, n = 0 : τ1 = tt τ22 = @t t ttt , vanishing tree, discard. (iii) T3: (1) s = 0, n = 1 : κ1 = tt@ t tt τ˜31 = tt@ t tt@ t tt , α(τ˜31 ) = 1 3 · 1 · 1 1! · 1 2 = 1 6 ; (2) s = 0, n = 2 : κ1 = κ2 = tt τ˜32 = tt@ t tt @t t t , α(τ˜32 ) = 1 3 · 1 · 1 2! · 1 · 1 = 1 6 ; 43 Bloch–Iserles equations (3) s = 1, n = 1 : κ1 = τ1 = tt τ˜33 = @t t ttt @ t tt , α(τ˜33 ) = 1 3 · (−1 2 ) · 1 1! · 1 · 1 = −1 6 ; (4) s = 1, n = 0 : τ1 = tt@ t tt τ˜34 = @t t ttt@ t tt , α(τ˜34 ) = 1 3 · (−1 2 ) · 1 0! · 1 · 1 2 = − 1 12 ; (5) s = 2, n = 0 : τ1 = τ2 = tt τ˜35 = @t tt t @ t tt t , vanishing tree, discard. Before we proceed further, let us clean up the set T3. We clearly see that τ˜33 is nothing but τ˜34 with opposite sign. The two trees can be aggregated into τ˜ 3 4 say, with the coefficient replaced by α(τ˜34 )−α(τ˜33 ). After trivial rotations (corresponding to commutation) we obtain three trees in the set T3, τ31 = tt@ t tt@ t tt , α(τ31 ) = 1 6 ; τ32 = tt@ t tt @t t t , α(τ32 ) = 1 6 ; 44 3.3. Representation by binary trees τ33 = @t t ttt@ t tt , α(τ33 ) = 1 12 . Using the tree formalism, we are now constructing the next “generation” of terms. (i) T4: (1) s = 0, n = 1 : i. κ1 = tt@ t tt@ t tt : τ˜41 = tt@ t tt@ t tt@ t tt , α(τ˜41 ) = 1 24 ; ii. κ1 = tt@ t tt @t t t : τ˜42 = tt@ t tt@ t tt @t t t , α(τ˜42 ) = 1 24 ; iii. κ1 = @t t ttt@ t tt : τ˜43 = tt@ t t@t t t t@t t t , α(τ˜43 ) = 1 48 ; (2) s = 0, n = 2 : i. κ1 = tt, κ2 = tt@ t tt : τ˜44 = tt@ t tt @t t t@t t t , α(τ˜44 ) = 1 16 ; 45 Bloch–Iserles equations ii. κ1 = tt@ t tt , κ2 = tt : τ˜45 = tt@ @ t tt @t t t @ t tt , α(τ˜45 ) = 1 16 ; (3) s = 0, n = 3 : κ1 = κ2 = κ3 = tt τ˜46 = tt@ t tt @t t t @ t tt , α(τ˜46 ) = 1 24 ; (4) s = 1, n = 0 : i. τ1 = tt@ t tt@ t tt : τ˜47 = @t t ttt@ t tt@ t tt , α(τ˜47 ) = − 148 ; ii. τ1 = tt@ t tt @t t t : τ˜48 = @t t ttt@ t tt @t t t , α(τ˜48 ) = − 148 ; iii. τ1 = @t t ttt@ t tt : τ˜49 = @t t tt@ tt tt@t t t , α(τ˜49 ) = − 196 ; (5) s = 1, n = 1 : i. τ1 = tt@ t tt , κ1 = tt : 46 3.3. Representation by binary trees τ˜410 = @@ t tt t@t t t t@t t t , vanishing tree, discard. ii. τ1 = tt , κ1 = tt@ t tt : τ˜411 = @t t ttt @ t tt@ t tt , α(τ˜411) = − 1 16 ; (6) s = 1, n = 2 : κ1 = κ2 = τ1 = tt : τ˜412 = @t t ttt @ t tt @t t t , α(τ˜412) = − 1 16 ; (7) s = 2, n = 0 : i. τ1 = tt , τ2 = tt@ t tt : τ˜413 = @t tt t @ t tt t@ t tt , α(τ˜413) = 1 96 ; ii. τ1 = tt@ t tt , τ2 = tt : τ˜414 = @@ t tt @ t tt tt@t t t , vanishing tree, discard. (8) s = 2, n = 1 : κ1 = τ1 = τ2 = tt : 47 Bloch–Iserles equations τ˜415 = @t tt t @ t tt t@ t tt , α(τ˜415) = 1 48 ; (9) s = 3, n = 0 : τ1 = τ2 = τ3 = tt : τ˜416 = @t tt t @ t tt @t t t t , α(τ˜416) = 0; vanishing tree, discard. Three trees are vanishing here. Also, with trivial rotations, few trees get cancelled or merge into other trees obeying the anti symmetry, for example, τ˜49 is −τ˜413 which is also equal to −(−τ˜415), similarly τ˜47 is nothing but τ˜411 with opposite sign and also τ˜48 = −τ˜412 and the trees τ˜43 , τ˜ 4 4 and τ˜ 4 5 satisfy Jacobi identity. Tiding up T4 and translating every tree in terms of commutators and curly brackets, the Taylor expansion of Ω(t) becomes Ω(t) = t{X0} + 1 2 t2{[{X0}, X0]} +t3( 1 6 {[{[{X0}, X0]}, X0]}+ 1 6 {[{X0}, [{X0}, X0]]} + 1 12 [{[{X0}, X0]}, {X0}]) +t4( 1 24 {[{[{[{X0}, X0]}, X0]}, X0]} + 1 24 {[{[{X0}, [{X0}, X0]]}, X0]}+ 1 24 {[{X0}, [{[{X0}, X0]}, X0]]} + 1 12 {[{[{X0}, X0]}, [{X0}, X0]]}+ 1 24 {[{X0}, [{X0}, [{X0}, X0]]]} + 1 24 [{[{[{X0}, X0]}, X0]}, {X0}] + 1 24 [{[{X0}, [{X0}, X0]]}, {X0}]) + . . . . In explanation of the simplification exercise of the above terms: As we stated that the terms get simplified obeying the anti symmetry and Jacobi identity; By anti symmetry we mean [A,B] = −[B,A], which is an easy exercise to verify for the reader, whereas A,B,C are said to satisfy the Jacobi identity if [A, [B,C]] + [B, [C,A]] + [C, [A,B]] = 0. As we 48 3.3. Representation by binary trees translate from the trees, τ˜43 → {[[{[{X0}, X0]}, {X0}], X0]} → {[[Ω2,Ω1], X0]}, τ˜44 → {[{X0}, [{[{X0}, X0]}, X0]]} → {[Ω1, [Ω2, X0]]}, τ˜45 → {[{[{X0}, X0]}, [{X0}, X0]]} → {[Ω2, [Ω1, X0]]} and hence use the Jacobi identity for simplification. Also, one may observe that by co- incidence the tree formation seems identical with [Ise02] by replacing @ t d → tt . 10−4 10−3 10−2 10−6 10−4 10−2 100 102 104 ∆t l 2 e r r o r Error for Ω[3] Error for Ω[4] Expected error O((∆t)4) Expected error O((∆t)3) Figure 3.4: Global error on logarithmic scale across an interval [0,1] with different time steps, after truncating the Taylor expansion up to fourth order terms. It seems quite tempting to introduce such a correspondence between the brackets {X} and [., N ] for the whole tree exercise. With this replacement we get correct terms up to third order but it gives different values when we reach next generation terms. The reason for this fact is that the difference between {X} and [., N ] is not only the brackets but the different signs, + and −, which appear implicitly in the expression and surely make difference as the complexity increases. Moreover, the number of independent fourth order 49 Bloch–Iserles equations terms in [Ise02] is 8, whereas the number of independent terms in Ω4 in BI system is 7. This explains that the fourth order terms obtained from the two approaches are different. Therefore, no such correspondence for the tree structure has been found between the two. We see that once a tree formulation is obtained, translating it into commutators and curly brackets costs less labour as compared to the complexity of manual computation. After translating the trees into mathematical expressions, we plot the error graph. Trun- cating the expansion up to fourth order terms, the error graph of the solution as compared to the MATLAB ode45 solver with built-in parameters, is shown in Figure 3.4. The error plot is generated by comparing against the theoretically expected error of O((∆t)3) and O((∆t)4). The experiments were performed on random 25 × 25 matrices. Clearly, this plot shows that the Magnus method is a fourth order method. 3.4 Concluding remarks One thus observes that the BI equations have a number of remarkable properties which motivate us to expand the equations. It has been seen in Figure 3.3 that these equations have Lie-Poisson structure. Figure 3.2 shows that the discretisation of the BI equations using Magnus expansion preserves the eigenvalues of the solution matrix. Also, it is observed in Figure 3.4 that the Lie group method using Magnus expansion is a fourth order method, here by fourth order we mean the truncation of Ω(t) upto the fourth power of t. By employing the shorthand of binary rooted trees for expansion terms, the computation is made affordable. This also lays a foundation to the explicit representation of the solution of the BI system. 50 Chapter 4 The dynamic generalised double bracket flow The subject matter of this chapter is generalised double bracket flow, that is, the matrix system of ordinary differential equations of the form X ′ = [[N,X] +M,X], t ≥ 0, X(0) = X0 ∈ Sym(n) (4.1) where N ∈ diag(n) and M ∈ so(n). We are concerned with analysing and discretising (4.1). Clearly, we see that it is of the form X ′ = [B(X), X], t ≥ 0, X(0) = X0 ∈ Sym(n), where B(X) : Sym(n) → so(n) is given by B(X) = [N,X] + M . This is a form of an isospectral flow. Further we observe that when M = 0, it reduces to the double bracket flow, which is explained in detail in section 4.1. Moreover, when N = 0, this reduces to an integrable system. Now the question arises as to how it behaves when both N and M are non-zero. In this chapter we analyse the dynamics of (4.1) for different values of N and M . We present the phase portraits for various non-zero values of N and M . Before proceeding further let us have a brief look at the double bracket flow in the following section. 4.1 Double bracket flow As discussed above, when M = 0, (4.1) reduces to the double bracket flow (dbf) X ′ = [[N,X], X], t ≥ 0, X(0) = X0 ∈ Sym(n). 51 The dynamic generalised double bracket flow Dbf was introduced by Brockett in [Bro91], and Chu and Driessel in [CD90]. The double bracket equations possess a number of features. First of all, it is a particular example of an isospectral flow [Ise02]. X ′ = [A(t,X), X], X(0) = X0 ∈ Sym(n), where A ∈ so(n), the Lie-algebra of n×n real skew-symmetric matrices. Then there exists a matrix function Q(t) ∈ SO(n) such that X(t) = Q(t)X0Q T (t) and X(t) has the same eigenvalues as X0 for all t ≥ 0. The second feature of dbf is that it is a gradient system, with a global Lyapunov function, therefore it is assured of convergence to a fixed point of the flow as t→ 0 [Bro91]. Given the potential function ψ(X) = 1 2 ‖X −N‖2F where ‖.‖F is the Frobenius norm and X ranges across all symmetric matrices orthogonally similar to X0, then it is shown in [CD90] that ∇ψ(X) = −[[N,X], X] and thus dbf is precisely gradient system X ′ = −∇ψ(X). Hence, the double-bracket flow minimizes ψ. X∞ = limt→∞X(t) is a local minimizer of ψ in [Bro91]. We can use the double-bracket flow to diagonalize real symmetric matrices, and to find their eigenvalues. It has been shown by Brockett that when N is a real diagonal matrix and both X0 and N have distinct eigenvalues, then X(t) tends exponentially to a diagonal matrix as t → +∞ and the eigenvalues are sorted accordingly to the diagonal entries of N . There are other applications including sorting lists and solving certain linear program- ming problems [Blo90]. In the double bracket flow, different choices of N correspond to special realization processes. Therefore, the appealing part of the dbf is also the flexibility caused by its dependence on the matrix N . If N = diag(1, 2, . . . , n) and Y is tridiagonal, then dbf gives the Toda flow on tridiagonal matrices [Blo90]. Double-bracket equation has been well understood from a theoretical point of view. When we talk about efficiently computing the solutions of dbf, we know there are several numerical methods, including the family of Lie-group methods [IMKNZ99] and [Zan97]. The idea is to write the solution in the form X(t) = Q(t)X0Q T (t) and, instead of computing X directly, in each step an orthogonal matrix Qk+1 is evaluated, so that Xk+1 = Qk+1XkQ T k+1. The matrix Qk+1 is chosen as the solution of the initial value problem Q ′ k+1 = A(t, Qk+1XkQ T k+1)Qk+1, Qk+1(kh) = I 52 4.2. Analysing the generalised double bracket flow which is a particular example of a Lie-group flow. In [Ise02], the dbf is represented in the form X(t) = exp(Ω(t))X0exp(−Ω(t)) and the Taylor expansion of Ω is constructed explicitly identifying individual expansion terms with certain rooted trees with bicolour leaves. In the following section, we proceed to the analysis of generalised double bracket flow (4.1), and we compute their solutions in the next chapter. 4.2 Analysing the generalised double bracket flow In this section we analyse the given system (4.1). As mentioned earlier, when N = 0, (4.1) reduces to an integrable system and when M = 0 it reduces to the double bracket flow. We will see how the dynamics change with different values of M and N . We also do the stability analysis taking examples of 2× 2 and 3× 3 matrices. 4.2.1 Analysing taking 2× 2 matrices Taking X0 =  a b b c  , N =  0 0 0 0  and M =  0 −1 1 0  we get X =  12(a+ c)− b sin 2t+ 12(a− c) cos 2t b sin 2t+ 12(a− c) cos 2t b sin 2t+ 12(a− c) cos 2t 12(a+ c) + b sin 2t− 12(a− c) cos 2t  which is an oscillatory matrix. Taking X =  a b b c  , N =  d1 0 0 d2  and M =  0 α −α 0  we get X ′ =  2b((d1 − d2)b+ α) (c− a)((d1 − d2)b+ α) (c− a)((d1 − d2)b+ α) −2b((d1 − d2)b+ α)  that is a ′ = 2b((d1 − d2)b+ α) (4.2) b ′ = (c− a)((d1 − d2)b+ α) c ′ = −2b((d1 − d2)b+ α) 53 The dynamic generalised double bracket flow Solving the equations (4.2) a ′ + c ′ = 0 or a+ c = constant In particular, we assume a+ c = 0, a = −c then a ′ = 2b((d1 − d2)b+ α) (4.3) b ′ = −2a((d1 − d2)b+ α) ⇒ a ′ b′ = b −a or −aa′ = bb′ or a2 + b2 = γ2 where γ2 = √ a20 + b 2 0. Therefore X =  a √ γ2 − a2√ γ2 − a2 −a  Taking a = γ cos θ, we have X = √ a20 + b 2 0  cos θ sin θ sin θ − cos θ  We can see that this is a circle. To find the stability of the equations, we will use the Jacobian matrix Now, fixed points of these equations are (i) b = 0, a = c, (ii) b = −αd1−d2 . We know that equilibrium is stable if all eigenvalues have negative real part, it is unstable if at least one eigenvalue has positive real part. Jacobian is given as J =  0 2((d1 − d2)b+ α) + 2b(d1 − d2) 0 −((d1 − d2)b+ α) (c− a)(d1 − d2) ((d1 − d2)b+ α) 0 −2((d1 − d2)b+ α)− 2b(d1 − d2) 0  54 4.2. Analysing the generalised double bracket flow Here, the determinant of Jacobian matrix is always zero. So, at least one of the eigenvalues is zero. Then the origin is not an isolated fixed point. Finding the eigenvalues of the Jacobian matrix at: (i) b = 0, a = c J =  0 2α 0 −α 0 α 0 −2α 0  Characteristic equation is given by JI − λI = 0∣∣∣∣∣∣∣∣∣∣∣ −λ 2α 0 −α −λ α 0 −2α −λ ∣∣∣∣∣∣∣∣∣∣∣ = 0 λ = 0, ±2αι, where ι = √−1. We see that system is unstable at the origin. It is stable for λ < 0. The nature of this equilibrium point says non hyperbolic fixed point is center. (ii) b = − αd1−d2 J =  0 −2α 0 0 (c− a)(d1 − d2) 0 0 2α 0  Characteristic equation is given by JI − λI = 0∣∣∣∣∣∣∣∣∣∣∣ −λ −2α 0 0 (c− a)(d1 − d2)− λ 0 0 2α −λ ∣∣∣∣∣∣∣∣∣∣∣ = 0 λ = 0, λ = (c− a)(d1 − d2) For stability we should have λ = (c− a)(d1 − d2) < 0, That means one of c− a and d1 − d2 should be negative for stability. This is a non hyperbolic fixed point. The dynamics of the model system (4.2) is showing stable limit cycle. The solution is presented as time series in Figure 4.1 and 2D and 3D phase plots showing the por- 55 The dynamic generalised double bracket flow traits with respect to (a, b), (a, c), (b, c) and (a, b, c) respectively in the Figures 4.2(i), 4.2(ii), 4.2(iii) and 4.2(iv). These plots are generated taking the initial condition X0 =  1 1 1 2  . Time 0 5 10 15 20 25 30 a; b; c -3 -2 -1 0 1 2 3 4 5 Figure 4.1: Showing time series of the model system (4.2) with d1 = 4, d2 = 5, α = 4. 4.2.2 An explicit solution taking 2× 2 matrices We also attempt to find an explicit solution for (4.1) taking 2× 2 matrices. Having again a look at the above system (4.2), we denote, δ = d1 − d2 6= 0 therefore rewriting (4.2), a′ = 2b(δb+ α) b′ = (c− a)(δb+ α) 56 4.2. Analysing the generalised double bracket flow a 0 0.5 1 1.5 2 2.5 3 b -1.5 -1 -0.5 0 0.5 1 1.5 a 0 0.5 1 1.5 2 2.5 3 c 0 0.5 1 1.5 2 2.5 3 (i) (ii) b -1.5 -1 -0.5 0 0.5 1 1.5 c 0 0.5 1 1.5 2 2.5 3 3 2 a 1 0-1.5 -1 -0.5 0 b 0.5 1 3 2 1 0 1.5 c (iii) (iv) Figure 4.2: Phase plots (a, b), (a, c), (b, c) and (a, b, c), respectively, of the model system (4.2) with d1 = 4, d2 = 5, α = 4. 57 The dynamic generalised double bracket flow c′ = −2b(δb+ α) ∴ a′ + c′ = 0 and a+ c = const. Let p(t) = a(t)− c(t) , q(t) = δb(t) + α; [or b(t) = q(t)−αδ ] ∴ p′ = a′ − c′ = 4b(δb+ α) = 4( q(t)− α δ ).q = 4 δ q(q − α) Similarly, q′ = δb′ = δ(c− a)(δb+ α) = δpq The general solution with MAPLE and trivial simplifications is p(t) = − c1[(64α 2 − 16c21)e2c2+2+c1 − 1] δ[(8αetc1+c2 + 1)2 − 16c21e2(tc1+c2)] q(t) = 4c21e tc1+c2 (8αetc1+c2 + 1)2 − 16c21e2(tc1+c2) (i) When c1 > 0 lim t→∞ q(t) = 0, limt→∞ p(t) = − c1 δ ( c1 6= 4α) (ii) When c1 < 0 lim t→∞ q(t) = 0, limt→∞ p(t) = c1 δ For ec2 → c2, we get the solution as p(t) = − c1[(64α 2 − 16c21)c2e2tc1 − 1] δ[(8αc2etc1 + 1)2 − 16c21c22e2tc1 ] q(t) = 4c21c2e tc1 (8αc2etc1 + 1)2 − 16c21c22e2tc1 58 4.3. Some interesting phase portraits At t = 0, the solution is p(0) = −c1[64α 2c22 − 16c21c22 − 1] δ[(8αc2 + 1)2 − 16c21c22] q(0) = 4c21c2 (8αc2 + 1)2 − 16c21c22 ∴ c1 = zeroes of Z2 = λ2p20 + 8αq0 − 4q20 ∴ c1 = ± √ λ2p20 + 8αq0 − 4q20 c2 = q0 2 ±λp0 √ λ2p20 + 8αq0 − 4q20 + λ2p20 + 4αq0 − 4q20 Hence, we have three cases, firstly, when we have λ2p20+8αq0−4q20 ≥ 0 this implies p(t), q(t)→ 0. In the second case, if λ2p20 + 8αq0 − 4q20 < 0, the system reduces to a limit cycle. And thirdly, in case λ2p20 + 8αq0 = 4q 2 0, this results in a Hopf bifurcation and vice versa. In Hopf bifurcation the real parts of a pair of complex conjugate eigenvalues of the least- stable equilibrium point increase through zero as a control parameter is varied through a critical value. A time-periodic solution may arise after bifurcation. 4.3 Some interesting phase portraits In the last section of this chapter, we present some phase portraits showing the remark- able behaviour of generalised double bracket flows. We have computed numerically the solution of the system for random 3 × 3 matrices using Lie group method using Mag- nus expansion. In the Figures 4.3–4.6, the phase portraits (X1,2, Xk,l) are displayed for (k, l) = (1, 1), (1, 3), (2, 2), (2, 3), (3, 3), with random initial condition. Consider the ma- trices N =  d1 0 0 0 d2 0 0 0 d3  and M =  0 h g −h 0 l −g −l 0  We discuss the following cases- We fix the matrices choosing d1 = 1, d2 = 2, d3 = 3 and h = 2, g = 3, l = 1. Setting a parameter p where p ∈ (0, 1), we generate the phase portraits for p ·N and (1− p) ·M for different values of p. Plots are generated for p = 0.9, 0.7, 0.1, and 0.0001, respectively in the Figures 4.3, 4.4, 4.5 and 4.6. In figure 4.6, each string is again a bunch of strings, it is shown in the last plot in Figure 4.6 by focusing on the enlarged image in one of the 59 The dynamic generalised double bracket flow plots. Similar behaviour (as in 4.6) is seen for other values of the form p = 10−n where n = {1, 2, 3, · · · }. We see Hopf bifurcation occurs in a periodic orbit of an autonomous flow. In this case the invariant curve corresponds to an invariant torus for the flow and attracting periodic orbits on the circle correspond to mode-locked periodic motion on the torus, whilst dense orbits correspond to quasi-periodic motion. 60 4.3. Some interesting phase portraits -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.45 0.5 0.55 0.6 0.65 0.7 0.75 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0 0.2 0.4 0.6 0.8 1 1.2 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 Figure 4.3: The phase portraits (X1,2, Xk,l) for (k, l) = (1, 1), (1, 3), (2, 2), (2, 3), (3, 3), respectively, with a random initial condition. Here by Xk,l we mean the kl th element of the matrixX. These plots are generated choosing d1 = 1, d2 = 2, d3 = 3 and h = 2, g = 3, l = 1 and p = 0.9. 61 The dynamic generalised double bracket flow -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.5 0.55 0.6 0.65 0.7 0.75 0.8 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 0.2 0.4 0.6 0.8 1 1.2 1.4 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 Figure 4.4: The phase portraits (X1,2, Xk,l) for (k, l) = (1, 1), (1, 3), (2, 2), (2, 3), (3, 3), respectively, with a random initial condition. Here by Xk,l we mean the kl th element of the matrixX. These plots are generated choosing d1 = 1, d2 = 2, d3 = 3 and h = 2, g = 3, l = 1 and p = 0.7. 62 4.3. Some interesting phase portraits -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 0 0.5 1 1.5 2 2.5 3 3.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 Figure 4.5: The phase portraits (X1,2, Xk,l) for (k, l) = (1, 1), (1, 3), (2, 2), (2, 3), (3, 3), respectively, with a random initial condition. Here by Xk,l we mean the kl th element of the matrixX. These plots are generated choosing d1 = 1, d2 = 2, d3 = 3 and h = 2, g = 3, l = 1 and p = 0.1. 63 The dynamic generalised double bracket flow -1 -0.5 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 3.5 -1 -0.5 0 0.5 1 1.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 -1 -0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 -1 -0.5 0 0.5 1 1.5 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 -1 -0.5 0 0.5 1 1.5 0 0.5 1 1.5 2 2.5 3 -0.7 -0.69 -0.68 -0.67 -0.66 -0.65 -0.64 -0.63 0.45 0.46 0.47 0.48 0.49 0.5 Enlarged image of the string Figure 4.6: The phase portraits (X1,2, Xk,l) for (k, l) = (1, 1), (1, 3), (2, 2), (2, 3), (3, 3), respectively, with a random initial condition. Here by Xk,l we mean the kl th element of the matrixX. These plots are generated choosing d1 = 1, d2 = 2, d3 = 3 and h = 2, g = 3, l = 1 and p = 0.0001. 64 4.4. Concluding remarks 4.4 Concluding remarks Thus, we observe the phase portraits of the generalised double bracket flow for various values of M and N . We observe that, for Mij = 0, X(t) → Xˆ (this is the case of dbf), and for Mij  1, X(t) → pi as the system has Hopf bifurcation and gives birth to the limit cycles. Limit cycles have been used to model the behaviour of many oscillatory sys- tems. The Hopf bifurcation underlies many spontaneous oscillations such as airfoil flutter and other wind-induced oscillations (e.g., those that caused the Tacoma-Narrows bridge collapse) in structural engineering systems, vortex shedding in fluid flow around a solid body at sufficiently high stream velocity, LCR oscillations in electrical circuits, relaxation oscillations (e.g., as described by the Van der Pol oscillator), the periodic firing of neurons in nervous systems (e.g., in the FitzHugh-Nagumo equation modelling these phenomena), oscillations in autocatalytic chemical reactions (e.g., the Belousov-Zhabotinsky reaction) as described by the Brusselator and similar models, oscillations in fish populations (as described by predator-prey models), periodic fluctuations in the number of individuals suffering from an infectious disease (as described by epidemic models), etc. For references see, [MM12] [VDM85] [Sac64]. Having seen the interesting applications of Hopf bifurcation, we could study the analysis of this system more in detail in the sequel. In the next chapter we move towards discretising the flow using Magnus series. 65 The dynamic generalised double bracket flow 66 Chapter 5 Discretisation of generalised double bracket flow This chapter extends the method of Magnus expansion to Lie-algebraic equations in the generalised double bracket flow (gdbf) (4.1). The idea is to write the solution of (4.1) in the form X(t) = eΩ(t)X0e −Ω(t), where instead of computing X at the first place, we obtain the Taylor expansion of Ω. Our goal is to determine the rules for finding the terms of Ω to an arbitrary accuracy. Following the same solution structure as we did for BI equations, first we convert the isospectral flow to a Lie-group flow and then translate it into a Lie-algebraic equation. This method preserves the isospectrality and gives the desired structure of the solution with large time steps. We solve the given system of differential equations using the Magnus expansion to obtain the Taylor expansion of Ω. The terms are represented by binary rooted trees with tri-colour leaves and an algorithm is formed to construct the next tree by recursion and to calculate the coefficient of each tree. The explicit representation of the solution of the given equation is derived and B(X) is represented in a finite “alphabet”. The representation as binary trees is very important because, as the number of terms in each iteration grows exponentially, the complexity of manual computation becomes prohibitive. By indexing the terms in the expansion with a subset of binary trees, it is convenient to derive explicit recurrence relations. We organise this chapter as follows. In section 5.1 we expand the solution of (4.1) using Magnus series and represent it in the form X(t) = eΩ(t)X0e −Ω(t), which ensures that the solution structure remains isospectral. Taylor expansion of Ω has been formed algorithmically from X0, N , M and linear combinations of their commutators. The graph for solution error and the error in eigenvalues are displayed. We determine the precise rules underpinning this process and obtain explicit numerical procedure for the approximation of Ω to arbitrary accuracy. The procedure to construct Taylor expansion of Ω is similar using the Magnus expansion and it employes the terminology of binary rooted trees. This 67 Discretisation of generalised double bracket flow is much more complicated and delicate because of the three “alphabet”, X0, N , M , and their commutators. In the next section 5.2 the terms are represented by binary rooted trees with tri-colour leaves and we acquire a rule to construct the next tree by recursion and to calculate the coefficient of each tree. Although this matrix system appears more complex and leads to the tri-colour leaves; it has been possible to formulate the explicit recursive rule. In the case of this matrix system, we have two types of trees and an interplay between them, which makes it more convoluted and difficult. A step by step algorithm is developed and recurrence formula is defined that enables us to compute the trees and coefficients explicitly. In section 5.2.1 we assemble all the information into a well defined numerical algorithm. Using this algorithm next generation trees are constructed and their coefficient are calculated in section 5.2.2. Translating the trees into X0, N , M and their commutators gives us the required expansion of Ω. 5.1 Expansion of the solution We know (see [IMKNZ99]) that no classical numerical methods can respect the Lie–group structure. Exceptionally, in case of orthogonal flows, symplectic Runge–Kutta methods preserve orthogonality [LD94]. However, because of being implicit, it is expensive to use those methods. We use the method which respects the structure. Rewriting (4.1) to discretise using Magnus expansion X ′ = [[N,X] +M,X], t ≥ 0, X(0) = X0 ∈ Sym(n) X ′ = [B(X), X], B(X) = [N,X] + M where B(X) : Sym(n) → so(n). It is standard to verify that X(t) = Q(t)X0Q T (t), t ≥ 0, where Q(t) ∈ SO(n) is the solution of Q′(t) = B(Q(t)X0QT (t))Q(t), Q(0) = I. (5.1) In a similar way as we did for Bloch–Iserles equations in chapter 2, our idea is to represent the solution of the linear equation (5.1) in the form Q(t) = eΩ(t), where Ω′ = ∞∑ 0 Br r! adrΩ([N, e ΩX0e −Ω] +M), Ω(0) = 0. (5.2) 68 5.1. Expansion of the solution Here Bm, m ∈ Z are Bernoulli numbers and adrΩ is an iterated commutator defined by ad0ΩA = A, ad 1 ΩA = [Ω, A], ad 2 ΩA = [Ω, [Ω, A]], ..., ad m ΩA = [Ω, ad m−1 Ω A], where [Ω, A] = ΩA−AΩ. Now, taking Ω(t) = ∑∞ m=0 Ωmt m gives Ω′(t) = ∞∑ m=0 (m+ 1)Ωm+1t m (5.3) and this implies ∞∑ m=0 (m+ 1)Ωm+1t m = ∞∑ r=0 Br r! adrΩ(t)([N, e ΩX0e −Ω] +M). (5.4) Comparing coefficients of t0, t1, t2... we get the values of Ω1,Ω2,Ω3... as follows Ω1 = [N,X0] +M, Ω2 = 1 2 [N, [Ω1, X0]], Ω3 = 1 3 [N, [Ω2, X0]] + 1 6 [N, [Ω1, [Ω1, X0]]] − 1 6 [Ω1,Ω2]. It is time to analyse the error graphs of the solution that is obtained using Lie group method. Before simplifying the above expansion, we present the preliminary error graph and error graph for eigenvalues of this method, as compared to the MATLAB ode45 solver with built-in parameters, in Figure 5.1 and Figure 5.2, respectively. In Figure 5.1 we display error in the solution of gdbf in the interval [0,1] for a range of different step sizes ∆t. The error plot is generated by truncating the expansion up to order three and is compared against the theoretically expected error of O((∆t)3). The experiments were performed on random matrices. In Figure 5.2 we calculate absolute error of eigenvalues of the two methods on a logarithmic scale. It is clearly seen that our method preserves the correct eigenvalues to machine accuracy. We calculated the error in eigenvalues with small and large time steps. Despite a large time-step in the Magnus method, the error in eigenvalues stays very close to machine precision while the solution obtained using ode45 quickly strays away in terms of eigenvalues as time increases. Which is the expected behaviour of Magnus method from the principles underlying our approach. 69 Discretisation of generalised double bracket flow "t 10-3 10-2 10-1 l 2 er ro r 10-6 10-4 10-2 100 102 Error for +[3] Expected error O ! ("t)3 " Figure 5.1: Global error on logarithmic scale across an interval [0,1] with different time steps, after truncating the Taylor expansion up to third order terms. t 10-3 10-2 10-1 100 101 l 1 er ro r in ei ge nv al ue s 10-15 10-10 10-5 100 Solution with ode45 Magnus with "t : 1/10 Magnus with "t : 1/160 Figure 5.2: Error plot showing absolute error of eigenvalues of the two methods on a logarithmic scale. 70 5.2. Representation by binary trees 5.2 Representation by binary trees Let us again look at the expansion Ω, we note that each term is written in just three ‘letters’ X0, N and M , hence belongs to a free structure generated by them. This seems relatively complex than the expansion of Bloch–Iserles equations in chapter 2, because there is an extra ‘letter’ M and it appears implicitly in every expression and making it more complex. To understand the expansion terms and the tree structure in a better way we tend to separate the commutators in every term when it comes with +M e.g. [., X +M ] = [., X] + [.,M ], that for sure is going to increase the number of trees for every term. We attempt to find the expansion of the solution in Taylor series of the form Ω(t) = ∞∑ r=1 tr ∑ τ∈Tr α(τ)Hτ (5.5) where Tr is the set of all binary trees of power r, Hτ is an expression constructed from X0, N and M according to rules implicit in the structure of the tree τ which will be explained next, and α is a scalar constant. We use binary rooted trees as a shorthand for expansion terms, an approach introduced by [IN99], also used for BI equations in chapter 3, that leads to a groundwork illuminating the structure of individual terms and their relationship. We choose to assign to the three ‘letters’, X0, N and M , the leaves in three different colours namely, black, white, and black-n-white respectively, t X0. and d N, and M. We define a function τ → Hτ from T = ⋃∞ r=1 Tr, a subset of binary rooted trees into n× n matrix functions by letting H s= X0 and, by induction, @ τ1 τ2 [Hτ1 , Hτ2 ], where Hτ1 and Hτ2 are already constructed expansion terms. We delve into the general convention for the correspondence between trees and expan- sion terms. Let us have again a look at the equation ∞∑ m=0 (m+ 1)Ωm+1t m = ∞∑ r=0 Br r! adrΩ(t) ( [N, eΩ(t)X0e −Ω(t)] +M ) . 71 Discretisation of generalised double bracket flow We know that eΩX0e −Ω = AdΩX0 = eadΩX0 = ∞∑ n=0 1 n! adnΩX0 Thus ∞∑ m=0 (m+ 1)Ωm+1t m = ∞∑ r=0 Br r! adrΩ(t) ( [N, ∞∑ n=0 1 n! adnΩ(t)X0] +M ) = ∞∑ r=0 Br r! adrΩ(t)[N, ∞∑ n=0 1 n! adnΩ(t)X0] + ∞∑ r=0 Br r! adrΩ(t)M = ∞∑ r=0 Br r! adrΩ(t) ∞∑ n=0 1 n! @ d adnΩ(t)X0 + ∞∑ r=0 Br r! adrΩ(t)M = B0 0! ( ∞∑ n=0 1 n! @ d adnΩ(t)X0 + M ) + B1 1! ( ∞∑ n=0 1 n! @ Ω(t)@ d adnΩ(t)X0 + @ Ω(t) M ) + B2 2! ( ∞∑ n=0 1 n! @ Ω(t)@ Ω(t)@ d adnΩ(t)X0 + @ Ω(t)@ Ω(t) M) + · · · . (5.6) where adnΩ(t)X0 = @ @ Ω(t) ad n−1 Ω(t)X0 . If we compare the coefficients of t0, t1, t2... we get the following representation 72 5.2. Representation by binary trees Ω1 = B0 0! 1 0! @ d X0 +M = @ d t + and 2Ω2 = B0 0! @ N @ Ω1 X0 + B1 1! @ Ω1 @ N X0 + B1 1! @ Ω1 M ︸ ︷︷ ︸ = @ N @ Ω1 X0 + B1 1! @ Ω1 [N,X0] +M︸ ︷︷ ︸ = @ N @ Ω1 X0 + B1 1! ( @ Ω1 Ω1 = 0 ) Therefore, Ω2 = 1 2 @ N @@ [N,X0]X0 + 1 2 @ N @ M X0 = 1 2 @ d @@ d t t + 1 2 @ d @ t Similarly, Ω3 = 1 6 @ d @@ d t@ t@d t + 1 6 @ d @@ d t@ t + 1 6 @ d @@ @ d t @ t@d t + 1 6 @ d @ @ t@d t + 1 6 @ d @@ @ d t @ t + 1 6 @ d @ @ t − 1 12 @@ @ d t @ d @ t@ d t − 1 12 @@ @ d t @ d @ t 73 Discretisation of generalised double bracket flow − 1 12 @ @ d @ t@ d t − 1 12 @ @ d @ t Clearly, we see that each tree here can be represented in the form Ts,n 3 τ = @ τ1 @ τ2 @ τs @ d @κ1 @ κ2 @ κn t , or Ts 3 τ = @ τ1 @ τ2 @ τs . (5.7) where s ∈ {0, 1, 2, 3, . . .}, n ∈ {0, 1, 2, 3, . . .}. Here, the trees τ1, τ2, . . . τs, κ1, κ2, . . . κn have been featured earlier in the expansion, where τi ∈ Tpi , and κj ∈ Tqj , i = 1, 2, ..., s, j = 1, 2..., n and p1 + p2 + ...+ ps + q1 + q2 + ...+ qn + 1 = r. It is also possible to deduce the explicit form of the constant α. Set α ( @ d t ) = 1 (5.8) and α ( ) = 1. (5.9) Let τ ∈ Tr, r ∈ N and suppose that α(τi) and α(κj) are known for i = 1, 2, . . . , s, j = 1, 2, . . . , n. Then α(τ) = 1 r Bs s! 1 n! s∏ i=1 α(τi) n∏ j=i α(κj), s, n ∈ N, (5.10) where Bs is the sth Bernoulli number. We have the elements set to generate the trees for the Taylor expansion. In the next subsection, we construct the recursive rule (5.7) and (5.10) step by step in an alternative way. 5.2.1 Constructing the elements (5.7) and (5.10) We may write the dexp equation as Ω ′ = dexp−1Ω ([N, e adΩX0] +M), t ≥ 0,Ω(t0) = O where eadΩ = Σ∞m=0( 1 m!ad m Ω ). We attempt to find the expansion of the solution in Taylor series of the form (5.5). As defined earlier in section 5.2, while constructing trees, we 74 5.2. Representation by binary trees commence by assigning X0 to a single node, i.e. a trivial tree, t X0. and d N, and M. We define a function τ → Hτ from T = ⋃∞ r=1 Tr, a subset of binary rooted trees into n× n matrix functions by letting H s= X0 and, by induction, @ τ1 τ2 [Hτ1 , Hτ2 ], where Hτ1 and Hτ2 are already constructed expansion terms. Let us define Pn(t) = ad n ΩX0, n ∈ Z+, Qs(t) = ad s Ω([N, e adΩX0] +M), s ∈ Z+, the Taylor expansions of Pn and Qs can also be written as nested commutators and their linear combinations, Pn(t) = ∞∑ k=n tk ∑ τ∈Pn,k β(τ)Hτ , Qs(t) = ∞∑ l=s tl ∑ τ∈Qs,l γ(τ)Hτ , Here, Hτ is an expression constructed from X0, N and M in the structure of the tree τ , and β(τ) and γ(τ) are scalars. The index sets Pn,k and Qs,l are formal collections of binary rooted trees. Our goal is to establish the recursive rule from Pn. We commence by establishing the rules to form the sets Pn,k, and the coefficient map β. Since P0 = X0, we obtain P0,0 = s ;β( s) = 1 Thus substituting (5.5), P1 = [Ω(t), P0(t)] = ∞∑ m=1 tm ∑ τ∈T α(τ)[Hτ , X0], therefore, for every k ≥ 1, τ ∈ P1,k ⇔ τ = @ κ1 t where κ1 ∈ Tk. 75 Discretisation of generalised double bracket flow In this case β(τ) = α(κ1). Here, we suppose that we have already determined the sets Pm,k and the underlying coefficient map β, for m ≤ n − 1 and k ≥ m, we will prove for m = n Pn = adΩPn−1 = [Ω, Pn−1] Pn = [ ∞∑ m=1 tm ∑ α(κ1)Hκ1 , ∞∑ k=n−1 tm ∑ β(κ2)Hκ2 ] = ∞∑ m=1 ∞∑ k=n−1 tm+k ∑∑ α(κ1)β(κ2)[Hκ1 , Hκ2 ], therefore, τ ∈ Pn,k ⇔ τ ∈ @ κ1 κ2 , where κ1 ∈ Tk1 , κ2 ∈ Pn−1,k2 , k1 + k2 = k β(τ) = α(κ1)β(κ2) Using induction and noting the fact that P0 = X0 can now be used to express in terms of trees from T1,T2, . . . ,Tn−1 Pn,k 3 τ = @ κ1 @ κ2 @ κn t ⇔ κl ∈ Tk,l, l = 1, 2, . . . , n, ∑ kl = k, β(τ) = n∏ i=1 α(κi). After expressing the Pn in T, we proceed to do the same for the functions Qs. We commence by noting eadΩX0 = ∞∑ n=0 1 n! adnΩX0 = ∞∑ n=0 1 n! Pn Therefore, Q0 = [N, e adΩ(t)X0] +M = ∞∑ r=0 1 r [N,Pn] +M = ∞∑ l=0 tl l∑ r=0 1 n! ∑ β(τ)[N,Hτ ] +M. hence, ∑ γ(τ)Hτ = l∑ r=0 1 n! ∑ β(τ)[N,Hτ ] +M. 76 5.2. Representation by binary trees We thus deduce that, for every l ∈ Z+, Q0,0 3 τ = @ d κ1 + , where κ1 ∈ ⋃l n=0 Pn,l and γ(τ) = 1 n!β(κ1). We can reformulate above as, Q0,l ∈ τ = @ d @κ1 @ κ2 @ κn t + , where κi ∈ Tk, ∑r i=l ki = l β(τ) = 1 n! n∏ i=1 α(κi). Continuing in the similar manner we get, Q1(t) = [Ω(t), Q0(t)] = ∞∑ m=0 ∞∑ l=0 tm+l ∑∑ α(τ1)γ(τ2)[Hτ1 , Hτ2 ], therefore, Q1,l 3 τ = @ τ1 @ d @κ1 @ κ2 @ κn t + @ τ1 , where κi ∈ Tk, qi + ∑s i=1 ki = l and γ(τ) = 1 n! s∏ i=1 α(τi) n∏ j=1 α(κj). (5.11) proceeding in inductive way using the fact, Qs = [Ω(t), Qs−1(t)], s ∈ N, we deduce that Qs,l 3 τ = @ τ1 @ τ2 @ τs @ d @κ1 @ κ2 @ κn t or @ τ1 @ τ2 @ τs , 77 Discretisation of generalised double bracket flow hence we obtained the tree structure (5.7), where τi ∈ Tqi , κj ∈ Tkj , ∑s i=1 qi+ ∑n j=1 kj = l and γ(τ) = 1 n! s∏ i=1 α(τi) n∏ j=1 α(κj). Now we work towards the derivation of (5.10). Having derived Qs, s ∈ Z+, we note that dexp−1Ω ([N, e ad Ω X0] +M) = ∞∑ s=0 Bs s! Qs. Integrating, we deduce from the dexpinv Ω(t) = ∞∑ s=0 Bs s! ˆ t t0 Qsdx = ∞∑ s=0 Bs s! ∞∑ m=s tm+1 m+ 1 ∑ γ(τ)Hτ ∞∑ m=1 tm m m−1∑ s=0 Bs s! ∑ γ(τ)Hτ Thus, ∑ α(τ)Hτ = tm m m−1∑ s=0 Bs s! ∑ γ(τ)Hτ , m ∈ N Thus we deduce that Tm = ⋃m−1 s=0 Qs,m−1, m ∈ N Let τ ∈ Tm,m ∈ N then ∃s ∈ 0, 1, · · · ,m− 1 such that τ ∈ Qs,m−1, we deduce α(τ) = 1 m Bs s! γ(τ). Also, from (5.11), γ(τ) = 1n! ∏s i=1 α(τi) ∏n j=1 α(κj). Together which result in the equation (5.10) α(τ) = 1 r Bs s! 1 n! s∏ i=1 α(τi) n∏ j=i α(κj), s, n ∈ N. Let us take an example, taking m = 1, T1 = Q0,0 Q0,0 = @ d τ1 + , τ1 ∈ P0,0, γ(τ) = β(τ1) P0,0 = { t}, β( t) = 1 therefore, T1 = Q0,0 = {@ d t , } α{ @ d t } = 1, α{ } = 1 This is the set of first order trees. We proceed in the similar way to obtain the next generation trees. Let us jump to the next section to grow the binary rooted trees using 78 5.2. Representation by binary trees this algorithm. 5.2.2 Growing trees using the algorithm Now we have the general pattern for recursion. Suppose that T1,T2, . . . ,Tr−1 are known and also the coefficient α in these sets is known. To construct Tr we note that every τ is of the form (5.7) for some s ∈ {0, 1, 2, 3, . . . , r−1} and n ∈ {0, 1, 2, 3, . . . , r−1}. For every such s and n we consider all the partitions p1 + p2 + ...+ ps + q1 + q2 + ...+ qn + 1 = r. For every partition we construct the tree τ in (5.7) and use (5.10) to determine the coefficient α. The trees which correspond to zero terms would be eliminated. Moreover, some trees can be replaced by linear combinations of other trees. Let us start from T1 (i) T1: For s = 0, n = 0, (1) τ11 = @ d t , α(τ11 ) = 1 1 · 1 · 10! = 1. (2) τ12 = , α(τ12 ) = 11 · 1 · 10! = 1. (ii) T2: (1) For s = 0, n = 1, i. When κ1 = @ d t , τ21 = @ d @@ d t t , α(τ21 ) = 1 2 · 1 · 1 1! = 1 2 ; ii. When κ1 = , τ22 = @ d @ t , α(τ22 ) = 1 2 · 1 · 1 1! = 1 2 ; (2) For s = 1, n = 0 : i. When τ1 = @ d t , τ23 = @ @ @ d t @ d t , vanishing tree, discard. τ24 = @ @ d t , α(τ24 ) = 1 2 −1 2 · 1 · 1 0! = −1 4 79 Discretisation of generalised double bracket flow ii. When τ1 = τ25 = @ @ d t , α(τ25 ) = 1 2 · −1 2 · 1 · 1 0! = −1 4 τ26 = @ , vanishing tree, discard. Before we proceed further, let us clean up the set T2. We clearly see that τ24 is nothing but τ25 with opposite sign. Therefore, both the trees get cancelled as the coefficient is also the same. Therefore, T2 contains, τ21 = @ d @@ d t t , α(τ21 ) = 1 2 ; τ22 = @ d @ t , α(τ22 ) = 1 2 . (iii) T3 : (1) when s = 0, n = 1 i. When κ1 = @ d @@ d t t τ31 = @ d @@ d t@ t@d t , α(τ31 ) = 1 6 ; ii. When κ1 = @ d @ t τ32 = @ d @@ d t@ t , α(τ32 ) = 1 6 ; (2) When s = 0, n = 2 i. When κ1 = κ2 = @ d t 80 5.2. Representation by binary trees τ33 = @ d @@ @ d t @ t@d t , α(τ33 ) = 1 6 ; ii. κ1 = , κ2 = @ d t τ34 = @ d @ @ t@d t , α(τ34 ) = 1 6 ; iii. κ1 = @ d t , κ2 = τ35 = @ d @@ @ d t @ t , α(τ35 ) = 1 6 ; iv. κ1 = κ2 = τ36 = @ d @ @ t , α(τ36 ) = 1 6 ; (3) For s = 1, n = 0 i. τ1 = @ d @@ d t t τ37 = @ @ @ d @ d t@t@t d , α(τ37 ) = − 1 12 ; τ38 = @ @ @ d @t@t d , α(τ38 ) = − 1 12 ; 81 Discretisation of generalised double bracket flow ii. τ1 = @ d @ t τ39 = @ @ @ d @ d t@t , α(τ39 ) = − 1 12 ; τ310 = @ @ @ d @t , α(τ310) = − 1 12 ; (4) For s = 1, n = 1 i. τ1 = κ1 = @ d t τ311 = @ @ @ d t @ d @ t@ d t , α(τ311) = − 1 6 ; ii. τ1 = @ d t , κ1 = τ312 = @ @ @ d t @ d @ t , α(τ312) = − 1 6 ; iii. τ1 = , κ1 = @ d t τ313 = @ @ d @ t@ d t , α(τ313) = − 1 6 ; iv. τ1 = κ1 = τ314 = @ @ d @ t , α(τ314) = − 1 6 ; (5) For s = 2, n = 0 i. When τ1 = τ2 = @ d t 82 5.2. Representation by binary trees τ315 = @ @ @ d t @@ @ d t @ d t , vanishing tree, discard; τ316 = @ @ @ d t @@ @ d t , α(τ316) = 1 36 ; ii. τ1 = , τ2 = @ d t τ317 = @ @@ @ d t @ d t , vanishing tree, discard; τ318 = @ @@ @ d t , α(τ318) = 1 36 ; iii. τ1 = @ d t , τ2 = τ319 = @ @ @ d t @ @ d t , α(τ319) = 1 36 ; τ320 = @ @ @ d t @ , vanishing tree, discard; iv. τ1 = τ2 = τ321 = @ @ @ td , α(τ321) = 1 36 ; τ322 = @ @ , vanishing tree, discard. Cleaning up the set T3. We clearly see that τ311 is nothing but τ37 with opposite sign. The two trees can be aggregated into τ37 say, with the coefficient replaced by α(τ 3 7 ) − α(τ311). Similarly, τ312 = −τ39 , τ313 = −τ38 and τ314 = −τ310. Proceeding further we see that τ315, τ317, τ320, τ 3 22 are vanishing whereas, τ 3 19 is nothing but τ 3 16 with opposite sign also τ 3 21 is negative of τ318. Therefore, simplifying the above generation, the trees having same coefficient get cancelled and the others get convoluted with adding or subtracting the respective 83 Discretisation of generalised double bracket flow coefficients. After trivial rotations (corresponding to commutation) we obtain the set T3, 1 6 @ d @@ d t@ t@d t , 1 6 @ d @@ d t@ t , 1 6 @ d @@ @ d t @ t@d t , 1 6 @ d @ @ t@d t 1 6 @ d @@ @ d t @ t , 1 6 @ d @ @ t , 1 12 @@ @ d @ d t@t@t d , 1 12 @@ @ d @t@t d , 1 12 @@ @ d @ d t@t , 1 12 @@ @ d @t . Next we generate the set T4 of fourth order trees. We get 102 trees with all the possible combinations but these get reduced with the simplifications. Some of trees are zero trees and few get cancelled or added up to the others because of anti symmetry and Jacobi identity. Also, when we have s = 3, we get zero coefficient with respect to those trees because the Bernoulli’s number B3 equals zero. Hence we lose the trees with coefficient containing B3. After all the possible simplifications we are left with the following thirty eight trees with respective coefficients. Therefore, T4 contains, 1 48 @ d @@ t d @ t@ d @ t@ d t , 1 48 @ d @@ t d @ t@ d @ t , 1 48 @ d @@ t d @ @ d t @ t@ td , 1 48 @ d @@ t d @ @ @ t td , 1 48 @ d @ t@ d @@ @ t @ td , 1 48 @ d @ t@ d @ @ t , 1 48 @ d @ t@@ @ d t @ d t@@ d t , 1 48 @ d @ t@@ @ d t @ d t@ , 84 5.2. Representation by binary trees 1 48 @ d @ t@@ @ d t@@ d t , 1 48 @ d @ t@@ @ d t@ , 1 16 @ d @@@ d t @ t@d @ t@ td , 1 16 @ d @@@ @ d @ t@d t@ t@ td , 1 16 @ d @@@ d t @ t@d @ t , 1 16 @ d @ @ d @ t @ t@ td , 1 16 @ d @ @ t@d @ t@ td , 1 16 @ d @@@ d @ t @ t@ td , 1 16 @ d @ @ t@d @ t , 1 16 @ d @@@ d @ t @ t , 1 24 @ d @@@ d t @@ @ td @ t@ td , 1 24 @ d @ @ @ t , 1 24 @ d @@@ d t @@ @ td @ t , 1 24 @ d @@@ d t @ @ t@d t , 1 24 @ d @@@ d t @ @ t , 85 Discretisation of generalised double bracket flow 1 24 @ d @ @@ @ td @ t@ td , 1 24 @ d @ @ @ t@d t , 1 24 @ d @ @@ @ td @ t , 1 24 @ @ @ d t @ d t @ @ d t @ @ d t , 1 24 @ @ d t @ @ d t @ @ d t , 1 24 @ @ @ d t @ d t @ @ d t @ , 1 24 @ @ d t @ @ d t @ , 1 24 @ @ @ d t @ d t @ @ d t @ @ d t , 1 24 @ @ d t @ @ d t @ @ d t , 1 24 @ @ @ d t @ d t @ @ @ d t , 1 24 @ @ @ d t @ @ @ d t , 1 24 @ @ @ d t @ d t @ @ d t @ , 1 24 @ @ @ d t @ @ d t @ , 1 24 @ @ @ d t @ d t @ @ , 1 24 @ @ @ d t @ @ . Translating the trees in terms of commutators, the Taylor expansion of Ω(t) becomes, 86 5.2. Representation by binary trees Ω(t) = t([N,X] +M) +t2( 1 2 [N, [[N,X], X]] + 1 2 [N, [M,X]] +t3( 1 6 [N, [[N, [[N,X], X]], X]] + 1 6 [N, [[N, [M,X]], X]] + 1 6 [N, [[N,X], [[N,X], X]]] + 1 6 [N, [M, [[N,X], X]]] + 1 6 [N, [[N,X], [M,X]]] + 1 6 [N, [M, [M,X]]] + 1 12 [[N, [[N,X], X]], [N,X]] + 1 12 [[N, [M,X]], [N,X]] + 1 12 [[N, [[N,X], X]],M ] + 1 12 [[N, [M,X]],M ]) +t4( 1 24 [N, [[N, [[N, [[N,X], X]], X]], X]] + 1 24 [N, [[N, [[N, [M,X]], X]], X]] + 1 24 [N, [[N, [[N,X], [[N,X], X]]], X]] + 1 24 [N, [[N, [M, [[N,X], X]]], X]] + 1 24 [N, [[N, [[N,X], [M,X]]], X]] + 1 24 [N, [[N, [M, [M,X]]], X]] + 1 48 [N, [[[N, [[N,X], X]], [N,X]], X]] + 1 48 [N, [[[N, [M,X]], [N,X]], X]] + 1 48 [N, [[[N, [[N,X], X]],M ], X]] + 1 48 [N, [[[N, [M,X]],M ], X]] + 1 16 [N, [[N,X], [[N, [[N,X], X]], X]]] + 1 16 [N, [[N, [[N,X], X]], [[N,X], X]]] + 1 16 [N, [[N,X], [[N, [M,X]], X]]] + 1 16 [N, [[N, [M,X]], [[N,X], X]]] + 1 16 [N, [M, [[N, [[N,X], X]], X]]] + 1 16 [N, [[N, [[N,X], X]], [M,X]]] + 1 16 [N, [M, [[N, [M,X]], X]]] + 1 16 [N, [[N, [M,X]], [M,X]]] + 1 24 [N, [[N,X], [[N,X], [[N,X], X]]]] + 1 24 [N, [M, [M, [M,X]]]] + 1 24 [N, [[N,X], [[N,X], [M,X]]]] + 1 24 [N, [[N,X], [M, [[N,X], X]]]] + 1 24 [N, [[N,X], [M, [M,X]]]] + 1 24 [N, [M, [[N,X], [[N,X], X]]]] + 1 24 [N, [M, [M, [[N,X], X]]]] + 1 24 [N, [M, [[N,X], [M,X]]]] + 1 24 [[N, [[N, [[N,X], X]], X]], [N,X]] + 1 24 [[N, [[N, [[N,X], X]], X]],M ] + 1 24 [[N, [[N, [M,X]], X]], [N,X]] + 1 24 [[N, [[N, [M,X]], X]],M ] + 1 24 [[N, [[N,X], [[N,X], X]]], [N,X]] + 1 24 [[N, [[N,X], [[N,X], X]]],M ] + 1 24 [[N, [M, [[N,X], X]]], [N,X]] + 1 24 [[N, [M, [[N,X], X]]],M ] 87 Discretisation of generalised double bracket flow + 1 24 [[N, [[N,X], [M,X]]], [N,X]] + 1 24 [[N, [[N,X], [M,X]]],M ] + 1 24 [[N, [M, [M,X]]], [N,X]] + 1 24 [[N, [M, [M,X]]],M ]) + · · · . Once we obtain a tree formulation, translating it into commutators costs less labour as compared to the complexity of manual computation of the terms. Subsequent to converting the trees into mathematical expressions, we plot the error graph. Truncating the expansion up to fourth order terms, the error graph of the solution as compared to the MATLAB ode45 solver with built-in parameters, is shown in Figure 5.3. The error plot is generated by comparing against the theoretically expected error of O((∆t)4). The experiments were performed on random initial matrices. "t 10-4 10-3 10-2 l 2 er ro r 10-12 10-10 10-8 10-6 10-4 10-2 100 Error for +[4] Expected error O ! (t)4 " Figure 5.3: Global error on logarithmic scale across an interval [0,1] with different time steps, after truncating the Taylor expansion up to fourth order terms. We thus observe that Lie group methods have the generality to be applied to various isospectral flows. Although the matrix system (4.1) appears more complex and leads to the tri-colour leaves; it has been possible to formulate the explicit recursive rule. This matrix system has two types of trees and an interplay between them, which makes it more difficult, yet we successfully develop a step by step algorithm and define a recurrence formula that enables us to compute the trees and coefficients explicitly. Figure 5.2 shows that the discretisation of the gdb equations using Magnus expansion preserves the eigenvalues of 88 5.2. Representation by binary trees the solution matrix. Also, it is observed in Figure 5.3 that the Lie group method using Magnus expansion is a fourth order method; here by fourth order we mean the truncation of Ω(t) upto the fourth power of t. By employing the shorthand of binary rooted trees for expansion terms, the computation is made affordable. This also lays a foundation to the explicit representation of the solution of the generalised double bracket flow. 89 Discretisation of generalised double bracket flow 90 Chapter 6 Conclusion In this thesis we considered solving isospectral flows, the ordinary differential equations of the form X ′ = [B(X), X], t ≥ 0, X(0) = X0, possessing the characteristic to leave the spectrum of the initial matrix X0 invariant in the solution space X. Evidently, these flows emerge in a number of remarkable applications ranging from classical mechanics to linear algebra. Traditional numerical methods do not preserve the isospectrality of the system. Therefore, we emphasise on solving the isospectral flows using Lie group methods which preserves the isospectrality of the solution space and gives the desired structure of the solution with large time steps. The emphasis is also given on presenting the terms using shorthand of binary rooted trees that enables us to derive a recursive relation to generate text trees and calculate their coefficients. That also helps us to reduce the manual computational cost. Also, using anti symmetry and Jacobi identity to reduce the trees to the independent terms, further reduces the cost. We observe in chapter 3 that BI equations consist of phenomenal properties and hence, motivate us to expand the equations. It has been seen that these equations have Lie- Poisson structure. The comparison of the error of the solution against ode45 after trun- cating up to third order and fourth order terms gives the desired result. Discretisation of the BI equations using Magnus expansion preserves the eigenvalues of the solution matrix. This favourable behavior of the Magnus method is to be expected from the principles un- derlying our approach. The shorthand of binary rooted trees for expansion terms, makes the computation affordable. This also lays a foundation to the explicit representation of the solution of the BI system. Motivated by the promising results of Lie group methods for BI equations, we applied the procedure on generalised double bracket flow and observed the favourable results. Though the flow in chapter 5 is much more complicated and has to result in tri-colour leaves while constructing trees, yet one observes that it has been possible to derive an explicit representation of the solution using Magnus expansion and binary rooted trees and preserve the desired structure leaving the eigenvalues invariant. 91 Conclusion Whereas in chapter 4, we observe the dynamical structure of the flow. It is observed that the system has Hopf bifurcation that gives birth to the limit cycles. There are a number of issues that can be addressed in future. An attempt can be made to formulate the rule to count the trees. It has been known that the number of trees having 2n nodes is calculated using Witt–Birkhoff theorem 12n ∑ d|n µ ( n d )( 2d d ) , where µ is Mo¨bius function, and the number of ordered trees of degree n + 1 is given by the Catalan numbers C(n) = (2n)!n!(n+1)! (see [MKK03] [Ise02] [Bou08]). Upper limit of independent trees for BI equations can be known using 12n ∑ d|n µ ( n d )( 2d d ) . Whereas in the case of gdbf, we have two types of trees, and an interplay between them, so it does not seem to be easy to determine the number of independent trees in each iteration but might be possible to formulate the rule. Lie group methods are geometric integrators and arise in a number of applications. Numerical implementation of Lie group methods has been seen in [Ise02] [Zan98] [IMKNZ99] [Ise99]. In this thesis the technique has been extended to nonlinear cases and we believe that this would have generality to be applied to a variety of problems in the same class. Also, the generalised form of BI equations X ′ = [N,Xm], t ≥ 0, X(0) = X0 ∈ Sym(n), N ∈ so(n), for different values of m can be studied, which for sure would have complexity in handling large number of brackets in the expansion depending upon how large is the power of X. 92 Bibliography [BBI+09] A. M. Bloch, V. Brnzanescu, A. Iserles, J. E. Marsden, and T. S. Ratiu, A Class of Integrable Flows on the Space of Symmetric Matrices, Communica- tions in Mathematical Physics 290 (2009), 399–435. (Cited on page 38). [BI05] A. M. Bloch and A. Iserles, Aspects of generalized double-bracket flows, Proc. Centre de Recherche Montreal, American Mathematical Society (2005), 65– 76. (Cited on page 13). [BI06] A.M. Bloch and A. Iserles, On an isospectral Lie–Poisson system and its Lie- algebra, Found. Comput. Math. (2006), 121–144. (Cited on pages 33 and 38). [Blo90] A. M. Bloch, Steepest descent, linear programming and Hamiltonian flows, Contemp. Math. AMS 114 (1990), 77–88. (Cited on page 52). [Bou08] N. Bourbaki, Lie groups and Lie algebras, vol. 7, Springer Science & Business Media, 2008. (Cited on page 92). [Bro91] R.W. Brockett, Dynamical systems that sort lists, diagonalize matrices, and solve linear programming problems, Linear Algebra and its applications 146 (1991), 79–91. (Cited on pages 13 and 52). [Cas04] F. Casas, Numerical integration methods for the double-bracket flow, Journal of Computational and Applied Mathematics 166 (2004), no. 2, 477 – 495. (Cited on pages 13 and 33). [CD89] M. Chu and K. R. Driessel, Can real symmetric Toeplitz matrices have ar- bitrary real spectra?, Tech. report, Idaho State University, 1989. (Cited on pages 13 and 20). [CD90] M.T. Chu and K.R. Driessel, The projected gradient method for least squares matrix approximations with spectral constraints, SIAM Journal on Numerical Analysis 27 (1990), no. 4, 1050–1060. (Cited on pages 13 and 52). [Chu84] M. Chu, The generalized Toda flow, the QR algorithm and the center manifold theory, SIAM Journal on Algebraic Discrete Methods 5 (1984), no. 2, 187–201. (Cited on pages 13 and 20). 93 BIBLIOGRAPHY [CIZ97] M.P. Calvo, A. Iserles, and A. Zanna, Numerical solution of isospectral flows, Mathematics of computation 66 (1997), no. 220, 1461–1486. (Cited on pages 14, 20, 21, and 34). [CIZ99] M P. Calvo, A. Iserles, and A. Zanna, Conservative methods for the Toda lattice equations, IMA journal of numerical analysis 19 (1999), no. 4, 509– 523. (Cited on pages 13, 14, 20, and 34). [Dem90] J. W. Demmel, Matrix Computations; Second Edition (Gene Golub and Charles F. Van Loan), SIAM Review 32 (1990), no. 4, 690–691. (Cited on page 21). [DNT83] P. Deift, T. Nanda, and C. Tomei, Ordinary differential equations and the symmetric eigenvalue problem, SIAM Journal on Numerical Analysis 20 (1983), no. 1, 1–22. (Cited on pages 13 and 20). [DRTW91] P. Deift, S. Rivera, C. Tomei, and D. Watkins, A monotonicity property for Toda-type flows, SIAM Journal on Matrix Analysis and Applications 12 (1991), no. 3, 463–468. (Cited on pages 13 and 20). [DRVV94] Luca Dieci, Robert D. Russell, and Erik S. Van V., Unitary integrators and applications to continuous orthonormalization techniques, SIAM Journal on Numerical Analysis 31 (1994), no. 1, 261–281. (Cited on page 21). [DVV99] Luca Dieci and Erik S Van V, Computation of orthonormal factors for funda- mental solution matrices, Numerische Mathematik 83 (1999), no. 4, 599–620. (Cited on page 21). [Fla74] H. Flaschka, The Toda lattice I, Phys. Rev. B9 (1974). (Cited on page 19). [FNO87] S. Friedland, J. Nocedal, and M. Overton, The formulation and analysis of numerical methods for inverse eigenvalue problems, SIAM Journal on Numer- ical Analysis 24 (1987), no. 3, 634–667. (Cited on pages 13 and 20). [GVL89] G.H. Golub and C.F. Van Loan, Matrix Computations 2nd Edn Johns Hopkins University Press, 1989. (Cited on page 21). [Har69] F. Harary, Graph Theory, Addison-Wesley, Reading MA, 1969. (Cited on page 28). [Hau06] F. Hausdorff, Die symbolische Exponentialformel in der Gruppentheorie, Berichte der Sa¨chsischen Akademie der Wißenschaften, Math. Phys. Klasse 58 (1906), 19–48. (Cited on page 28). 94 BIBLIOGRAPHY [HJ91] R.A. Horn and C.R. Johnson, Topics in Matrix Analysis, Cambridge Univer- sity Press, Cambridge, 1991. (Cited on page 34). [IMKNZ99] A. Iserles, H.Z. Munthe-Kaas, S.P. Nørsett, and A. Zanna, Lie-group methods, Acta Numerica (1999). (Cited on pages 17, 21, 52, 68, and 92). [IN99] A. Iserles and S.P. Nørsett, On the Solution of linear differential equations in Lie groups, Phil. Trans Royal Soc. A 357 (1999), 983–1019. (Cited on pages 27, 28, 37, 40, and 71). [Ise99] A. Iserles, Expansions that grow on trees, Notices Amer. Math. Soc 49 (1999), 430–440. (Cited on pages 17, 27, 29, and 92). [Ise02] A. Iserles, On the discretization of double-bracket flows, Foundations of Com- putational Mathematics 2 (2002), no. 3, 305–329. (Cited on pages 13, 33, 36, 37, 49, 50, 52, 53, and 92). [Kau16] A. Kaur, On solving an isospectral flow, Journal of Computational and Ap- plied Mathematics 308 (2016), 263 – 275. (Cited on page 33). [Lag91] J. Lagarias, Monotonicity properties of the Toda flow, the QR-flow, and sub- space iteration, SIAM Journal on Matrix Analysis and Applications 12 (1991), no. 3, 449–462. (Cited on pages 13 and 20). [Lax68] P. D. Lax, Integrals of Nonlinear Equations of Evolution and Solitary Waves, Commun. Pure Appl. Math. 21 (1968), 467–490. (Cited on page 17). [LD94] E. S. V. Vleck L. Dieci, R. D. Russell, Unitary Integrators and Applications to Continuous Orthonormalization Techniques, SIAM Journal on Numerical Analysis 31 (1994), no. 1, 261–281. (Cited on pages 21 and 68). [LT06] L. Li and C. Tomei, C.: The complete integrability of a Lie-Poisson system proposed by Bloch and Iserles, Intern. Math. Res. Notes 64949 (2006). (Cited on page 38). [Mag54] W. Magnus, On the exponential solution of differential equations for a linear operator, Comm. Pure Appl. Math. 7 (1954), 649–673. (Cited on pages 27, 28, and 35). [Man76] S.V. Manakov, Note on the integration of Euler’s equations of the dynamics of an n-dimensional rigid body, Funct. Anal. Appl. 10 (1976), 328–329. (Cited on page 34). [MKK03] H. Munthe-Kaas and S. Krogstad, On enumeration problems in liebutcher theory, Future Generation Computer Systems 19 (2003), no. 7, 1197 – 1205, 95 BIBLIOGRAPHY Selected papers on Theoretical and Computational Aspects of Structural Dy- namical Systems in Linear Algebra and Control. (Cited on page 92). [MM12] J. E. Marsden and M. McCracken, The Hopf bifurcation and its applications, vol. 19, Springer Science & Business Media, 2012. (Cited on page 65). [Nan82] T. Nanda, Isospectral flows on band matrices, Ph.D. thesis, New York Uni- versity, 1982. (Cited on pages 13 and 20). [Sac64] R. J. Sacker, On invariant surfaces and bifurcation of periodic solutions of ordinary differential equations, Tech. report, DTIC Document, 1964. (Cited on page 65). [Sym82] W. W. Symes, The QR algorithm and scattering for the finite nonperiodic Toda lattice, Physica D Nonlinear Phenomena 4 (1982), 275–280. (Cited on pages 13, 19, and 20). [Tod81] M. Toda, Theory of Nonlinear Lattices, Springer-Verlag, Berlin, 1981. (Cited on pages 13, 18, and 19). [Tod89] , Nonlinear Waves and Solitons, Mathematics and Its Applications (Japanese Series) (1989). (Cited on pages 13 and 18). [VDM85] J. Van Der Meer, The hamiltonian Hopf bifurcation, Springer, 1985. (Cited on page 65). [Zan97] A. Zanna, Lie-group methods for isospectral flows, University of Cam- bridge Department of Applied Mathematics and Theoretical Physics-Report- DAMTP NA (1997). (Cited on page 52). [Zan98] A. Zanna, On the numerical solution of isospectral flows, Ph.D. thesis, Uni- versity of Cambridge, 1998. (Cited on pages 14, 17, 20, 34, and 92). 96