Grassmannian spectral shooting

We present a new numerical method for computing the pure-point spectrum associated with the linear stability of coherent structures. In the context of the Evans function shooting and matching approach, all the relevant information is carried by the flow projected onto the underlying Grassmann manifold. We show how to numerically construct this projected flow in a stable and robust manner. In particular, the method avoids representation singularities by, in practice, choosing the best coordinate patch representation for the flow as it evolves. The method is analytic in the spectral parameter and of complexity bounded by the order of the spectral problem cubed. For large systems it represents a competitive method to those recently developed that are based on continuous orthogonalization. We demonstrate this by comparing the two methods in three applications: Boussinesq solitary waves, autocatalytic travelling waves and Ekman boundary layer.

For such problems the matching condition is typically a discriminant known as the Evans function (Evans [34], Alexander, Gardner and Jones [2]) or miss-distance function (Pryce [85], Greenberg and Marletta [39]).It measures the degree of (possibly transversal) intersection of the stable and unstable solution subspaces satisfying the longitudinally separated far field boundary data.The stable subspace decays in the direction of wave propagation whilst the unstable subspace decays in the opposite direction.Equivalently, the Evans function is the determinant of the set of solution vectors spanning both subspaces.With this end-goal matching condition in mind, the problem boils down to how to numerically construct the solution subspaces in a robust fashion, as well as where to match longitudinally.This is especially difficult for large scale problems.These might either emerge from high order systems, or more specifically, we envisage the stability of nonlinear travelling waves with multi-dimensional structure; for example wrinkled fronts travelling in a fixed longitudinal direction.The transverse structural information can be projected onto finite dimensional transverse basis, generating a large linear spectral problem posed on the one-dimensional longitudinal coefficient function set (see Ledoux, Malham, Niesen and Thümmler [66]).Hitherto such large problems could not be solved by shooting and matching and were resolved by projecting the whole problem onto a finite basis and solving the resulting large algebraic eigenvalue problem.However recently, in the context of the Evans function, Humpherys and Zumbrun [49] proposed continuous orthogonalization as a viable approach to help make large scale problems amenable to shooting and matching.Here we provide our own answer.
We propose the new Grassmann Gaussian elimination method (GGEM) which resolves several numerical problems all in one, in particular it: (1) Evolves the solution along the underlying Grassmann manifold avoiding representation singularities.(2) Retains analyticity in the spectral parameter.
(3) Allows for matching anywhere in the longitudinal computational domain.(4) Has polynomial complexity, operations are of the order of the size of the system cubed.(5) Naturally evolves the solution in what in practice is the optimal coordinate representation (generated by optimal partial pivoting).
That we have a numerical method that avoids representation singularities for nonautonomous systems is new.The idea is as follows.Given data that lies in a suitable coordinate patch of the Grassmann manifold, pullback to the Stiefel manifold.Evolve the solution one steplength along the Stiefel manifold either directly using a Runge-Kutta method or a Lie group method.Then project onto a suitable and possibly different coordinate patch of the Grassmann manifold using optimal Gaussian elimination.In this last step, the practical method we propose picks a quasi-optimal coordinate patch in which to best represent the solution in the Grassmann manifold (see below and Section 5).
Note that on first inspection the Stiefel manifold is the direct natural setting for the stable and unstable solution subspaces.Afterall in each case we have to construct the full set of solutions to a large system of differential equations satisfying the correct respective asymptotic boundary conditions in the far field.Each solution set represents a curve in the Stiefel manifold of dimension commensurate with the size of the solution set.That the spectral problems are linear means that all the relevant spectral information can be reconstructed from the flow in the corresponding Grassmann manifold.That the matching condition is determinental, means that we only need the Grassmann flow information-see Martin and Hermann [71] and Brockett and Byrnes [18] where this reduction was first considered for autonomous linear control problems (in practice we will also need to retain a complex scalar field to ensure analytic dependence on parameters).This reduction is crucial because long-range integration along the Stiefel manifold has been problematic (due to multiple distinct exponential growth and decay rates) and one of the simplest resolutions in the Evans function context was to use Plücker coordinates-whilst ignoring the Plücker relations (more on these below).
Second, retaining analyticity away from the essential spectrum is standard for any shooting method; we prove analyticity in Section 6.This allows for a global search for eigenvalues in that region by numerically computing the change in argument of the Evans function round any closed contour.Invoking the argument principle, this integer value represents the number of zeros of the Evans function, and hence the number of eigenvalues counting multiplicity, inside the contour; see Brockett and Byrnes [18], Alexander, Gardner and Jones [2] and Ying and Katz [99].
Third, that our method allows matching anywhere in the longitudinal domain is new.We provide substantive numerical evidence.Previously, other than in trivial cases, most numerical practioners used the common-sense rule of thumb of integrating the spectral problem from both ends of the longitudinal domain and matching at a point roughly centered on the front (which is assumed to be localized).When solving the linear problem with Plücker coordinates, it was first important to rescale for the far field spatial behaviour to neutralize its total exponential growth.Integrating from the far field initial conditions (a subset of the spatial eigenvectors), the solution remains roughly constant until the coefficient matrix starts to reveal its nonautonmous character due to the integration step impinging on the front.Accuracy is retained whilst integrating through the front, but thereafter the problem becomes stiff.The issue is that the numerical methods cannot resolve the simultaneous exponential growth and decay character associated with the other far-end stable and unstable subspaces.
Fourth, having polynomial complexity is essential and we provide here a new alternate.After Humpherys and Zumbrun [49] introduced their continuous orthogonalization method in this context, which also has polynomial complexity, any new numerical spectral shooting method should have this property and also be competitive.Previous successful methods used Plücker coordinates, ignoring the quadratic Plücker relations (see Section 2).They solved the flow for the corresponding linear vector field in the higher dimensional Plücker embedding space.Details of this Plücker coordinate or compound matrix approach can be found in, for example, Alexander and Sachs [3], Brin [16,17] and Allen and Bridges [5].Unfortunately the number of Plücker coordinates typically grows exponentially with the order of the original system, and so this approach cannot be used for medium to large order systems.However the continuous orthogonalization method of Humpherys and Zumbrun, and our method, are especially suited to large scale problems.
Fifth, our method for dynamic practical optimal coordinate representation is new.Given data on the Stiefel manifold, for example generated by advancing the solution one steplength along the Stiefel manifold, how can we project down onto the Grassmann manifold using the best representation patch possible?The idea is as follows.The natural map projection from the Stiefel to the Grassmann manifolds is a linear fractional map (Milnor and Stasheff [74]; Martin and Hermann [71]).This map represents the action of equivalencing by transformations whose rank matches that of the Stiefel manifold (this takes us from the space of frames to the space of spaces spanning those frames).The Stiefel manifold has a non-square matrix representation.Projection onto the Grassmann manifold corresponds to equivalencing by a full rank submatrix of the non-square Stiefel matrix-this renders the corresponding submatrix as the identity matrix.We are free to choose which submatrix to equivalence by, each distinct choice corresponds to the matrix representation of a coordinate patch on the Grassmann manifold.We can use Gaussian elimination, via elementary column operations, to render any given full rank submatrix of the Stiefel matrix as the identity matrix.The key is to try to pick the full rank submatrix which has the largest determinant-this corresponds to choosing the Grassmannian patch that gives the best representation for the projection from the Stiefel to Grassmann manifold.Ideally we would check the size of every full rank submatrix of the Stiefel matrix and equivalence by the one with the largest determinant.However this is an NP problem (equivalent to using the Plücker coordinates described above).We provide a practical solution of polynomial complexity.The method maximises the pivot used at each step of the Gaussian elimination process.In the current context, we call it quasi-optimal Gaussian elimination.
Our paper is organised as follows.In Section 2, we provide a tailored review of Grassmann manifolds and their representation.We then show, in Section 3 how flows generated by linear vector fields on the Stiefel manifold, produce a natural flow on the underlying Grassmann manifold that is decoupled from the flow through the remaining fibres.We subsequently show how this leads to using Riccati systems to resolve spectra, but that singularities that develop in the Riccati flows present spectral matching problems.In Section 4 we introduce two new practical approaches to avoiding these representation singularities.One is the idea behind our main method, the Grassmann Gaussian elimination method.The other is a modification of the Riccati approach that changes the coordinate patch when deemed necessary.Also in this section we show the connection between the Riccati and continuous orthogonalization approaches.We present our proposed Grassmann Gaussian elimination method fully in Section 5, including details of how in practice to choose the quasi-optimal Grassmannian coordinate representation patch.We review the Evans function in Section 6 and discuss further simple practical numerical refinements that retain analyticity and prevent potential numerical overflow.We then implement and compare all the competing numerical methods in Section 7 in three distinct applications.Finally in Section 8 we conclude and present future directions.
2. Review: Grassmann manifolds 2.1.Grassmann and Stiefel manifolds.A k-frame is a k-tuple of k ≤ n linearly independent vectors in C n .The Stiefel manifold V(n, k) of k-frames is the open subset of C n×k of all k-frames centred at the origin.The set of k dimensional subspaces of C n form a complex manifold Gr(n, k) called the Grassmann manifold of k-planes in C n (see Steenrod [88,p. 35] or Griffiths and Harris [40,p. 193]).
The fibre bundle π : V(n, k)→Gr(n, k) is a principle fibre bundle.For each y in the base space Gr(n, k), the inverse image π −1 (y) is homeomorphic to the fibre space GL(k) which is a Lie group-see Montgomery [76, p. 151].The projection map π is the natural quotient map sending each k-frame centered at the origin to the k-plane it spans-see Milnor In other words, U i is the set of k-planes Y ∈ Gr(n, k) such that the k × k submatrix of one, and hence any, matrix representation of Y is nonsingular (representing a coordinate patch labelled by i).
Any element of U i has a unique matrix representation y i • whose ith k × k submatrix is the identity matrix.For example, if i = {1, . . ., k} then any element of U {1,...,k} can be uniquely represented by a matrix of the form , where ŷi,j ∈ C for i = k + 1, . . ., n and j = 1, . . ., k. Conversely, a n × k matrix of this form represents a k-plane in U i .Each coordinate patch U i is an open, dense subset of Gr(n, k) and the union of all such patches covers Gr(n, k).For each i, there is a bijective map ϕ i : U i → C (n−k)k given by Each ϕ i is thus a local coordinate chart for the coordinate patch U i of Gr(n, k).
Since u i,i ′ represents the transformation between representative patchs and depends holomorphically on y i • , we deduce ϕ i •ϕ −1 i ′ is holomorphic.Note that Gr(n, k) has a structure of a complex manifold (see Griffiths and Harris [40,p. 194]).Further the unitary group U(n) acts continuously and surjectively on Gr(n, k).Hence Gr(n, k) is compact and connected.Lastly, the general linear group GL(n) acts transitively on Gr(n, k) and it is a homogeneous manifold isomorphic to GL(n)/GL(n − k) × GL(k) (see Chern [22, p. 65] or Warner [97, p. 130]).

Plücker embedding.
There is a natural map, the Plücker map, ; here P n denotes the complex projective space of dimension n.See Griffiths and Harris [40] or Coskun [26]  A natural coordinatization of P k C n is through the determinants of all the k × k submatrices of Y , normalized by a chosen minor characterized by an index i, hence These minor determinants-the Plücker coordinates-are not all independent, indeed, they satisfy quadratic relations known as the Plücker relations (which may themselves not all be independent).The image of the Plücker map p is thus the subspace of P ( n k )−1 cut out by the quadratic Plücker relations.

Grassmannian flows
3.1.Tangent space decomposition.Recall that we can consider the Stiefel manifold as a principle fibre bundle π : V(n, k) → Gr(n, k).Our goal here is to characterize the induced decomposition of the tangent space T Y V(n, k) for Y ∈ V(n, k).
We can decompose the tangent space T Y V(n, k) into horizontal and vertical subspaces (see, for example, Montgomery [76, p. 149]) The horizontal subspace H Y is associated with the tangent space of the Grassmannian base space, while the vertical subspace V Y is associated with the fibres homeomorphic to GL(k).Let us choose the coordinate patch representation U i for Gr(n, k) for some i = {i 1 , . . ., i k } ⊂ {1, . . ., n}.Let P i denote the projection matrix of size n × n that contains zeros everywhere except at positions (i l , i l ) for l = 1, . . ., k where it contains ones.Note that one can additively decompose any given tangent vector V = P i V + P i • V .Hence we have 3.2.Fibre bundle flow.Suppose we are given a vector field on the Stiefel manifold where we now think of a, b, c and d as functions of x, u and ŷ.
Proof.Using that Y = y i • u the ordinary differential system Y ′ = V (x, Y ) becomes where A i represents the submatrix obtained by restricting the matrix A(x, y i • u) to its ith columns.Applying the projections P i • and P i to both sides of this equation, respectively generates the equations for ŷ and u shown.Note that ŷ is the projection of y i • onto its i • th rows, as well as its image under the coordinate chart ϕ i .

Riccati flow.
A natural decoupling of the flow on the base space Gr(n, k) from the flow on the fibres GL(k) occurs when the vector field V is linear, i.e. when The following corollary is immediate from Theorem 1.
Corollary 1.If the vector field V is linear so that A = A(x) only, then: (1) The flow on the base space Gr(n, k) decouples from the flow evolving through the fibres GL(k)-the flow in the fibres is slaved to that in the base space; (2) For a fixed coordinate patch U i index by i, in the coordinate chart variables ŷ = ϕ i • y i • , the flow is governed by the Riccati differential system: Suppose we are required to determine the flow generated by a linear nonautonomous vector field defined on the Stiefel manifold V(n, k).The first conclusion in the corollary implies that all the relevant information is carried in the flow in the Grassmann manifold Gr(n, k), and the flow through the fibres GL(k) can be completely determined a-posteriori from the Grassmannian flow.The second conclusion suggests that if we fix a coordinate patch, then the flow in the Grassmannian can be determined from the solution to the Riccati system for ŷ.If required, we can solve the differential system for u in Theorem 1 to determine Y , thereby completely resolving the flow generated by the linear vector field V on the Stiefel manifold.Provided ŷ remains finite, this approach in fact works.
The problem is that, though Y = y i • u must be globally finite as it is generated by a linear vector field (with globally smooth coefficients), the Riccati solution ŷ itself can become singular in a finite integration interval.Of course simultaneously the determinant of u itself becomes zero.The solution on the Grassmannian does not become singular.The issue is representation.(Note that the flow on GL(k) is linear but rank is not preserved because its coefficients depend on the Riccati flow.) The Riccati flow is a flow in a given fixed coordinate chart indexed by i, which is chosen at the start of integration.Given an initial element in V(n, k), we pick a (good) coordinate patch U i for Gr(n, k), this fixes the Grassmannian representation y i • .Projecting onto the i • th rows of y i • , or equivalently looking at the image under the coordinate chart ϕ i , generates ŷ ∈ C (n−k)k .The Riccati flow is the flow in the Euclidean chart image space C (n−k)k .Each coordinate patch U i is dense in Gr(n, k).Therefore in numerical computations, the Riccati solution ŷ is likely to leave and return to the patch U i across any discrete integration step that staddles a representation singularity (generating a large but finite solution ŷ the other side).
With this in mind, Schiff and Shnider [90] suggested the following method that integrates through singularities in the Riccati flow.Fix a coordinate patch U i with index i.The general linear group GL(n) acts transitively on V(n, k) (and also Gr(n, k)): the left Lie group action Λ : GL The Möbius Lie group action of S, we have: Thus, given data ŷ0 in the Euclidean chart image space C (n−k)k , pullback to the Lie group GL(n), to the identity element I n , using the Möbius Lie group action map µ ŷ0 .Advance the solution across one integration step in the Lie group generating the element S ∈ GL(n).Schiff and Shnider used a Neumann/Runge-Kutta method to do this, but a Lie group method could also be used.Now push forward to ŷ = µ ŷ0 • S ∈ C (n−k)k using the Möbius Lie group action map.This takes you back to the chart corresponding to the original coordinate patch U i .
This method integrates through singularities in the Riccati flow.However, we are still left with another associated practical problem.If the linear vector field V depends on a parameter we wish to vary, the singularities in the Riccati flow can drift in the domain of integration and impinge on the matching position.

Practical Grassmann integration
Our goal in this section is construct numerical methods that integrate the flow associated with the push forward of the linear vector field V onto the Grassmannian manifold, whilst avoiding the representation singularities that occur in the Riccati approach.The solution is to change patch continuously in some optimal fashion or whenever the coordinate patch becomes a poor representation.
The following diagram helps map out the two new strategies we suggest.
Continuous optimal patch evolution.The idea behind the Grassmann Gaussian elimination method (GGEM) is this.Given data in Gr(n, k) that lies in a given patch Y 0 ∈ U i identified by i, pullback to the Stiefel manifold using the identity map-note Advance the solution one integration step along the Stiefel manifold V(n, k), using say a classical Runge-Kutta method, thus generating the element Y ∈ V(n, k).
Now with the next step solution on the Stiefel manifold, we use quasi-optimal Gaussian elimination with partial pivoting (QOGE) to decompose Y = y (i ′ ) • u and project onto the coordinate patch U i ′ producing the element y (i ′ ) • ∈ U i ′ .The quasioptimal Gaussian elimination process (described in Section 5) picks out a suitable coordinate patch U i ′ to represent the solution, which may be different than the original patch U i .
To ensure we numerically remain with the Stiefel manifold, rather than use a Runge-Kutta method, we might use a Lie group method as follows (see Munthe-Kaas [77]).Pullback the data Y 0 ∈ V(n, k) from the Stiefel manifold to the general linear group GL(n), via the action map Λ Y0 .The corresponding element in GL(n) is naturally the identity element I n .Subsequently pullback, via the exponential map, to the zero element in the corresponding Lie algebra, i.e. o ∈ gl(n).Evolve the solution on the Lie algebra to σ ∈ gl(n) using the Magnus expansion (Magnus [68]; also see Iserles, Munthe-Kaas, Nørsett and Zanna [54]).Pushforward from gl(n) to GL(n), via the exponential map, producing the Lie group element S = exp σ ∈ GL(n).Now pushfoward to the Stiefel manifold V(n, k) via the Lie group action map Λ Y0 .For more details on Lie group methods for Stiefel manifolds, see Krogstad [61] and Celledoni and Owren [20].

4.2.
Riccati flow with patch swapping.Since we are interested in constructing the flow generated by the push forward of the vector field V onto the Grassmann manifold, we can avoid singularities in any given Riccati system chart flow associated with a given Grassmannian coordinate patch, by simply changing patch when the solution representation appears to become poor.In particular, we could either change the coordinate patch when the: • Norm ŷ ∞ becomes too large (the easier and preferred approach we take); • Determinant det u becomes too small (this involves constructing u as integration proceeds and we want to avoid carrying information unnecessarily).Hence if say ŷ ∞ exceeds a prescribed tolerance at the end of one step, to change patch, we apply the quasi-optimal Gaussian elimination method (QOGE) to the matrix y i • = ϕ −1 i • ŷ.This identifies a new patch and index i ′ to use for the next set of successive steps until ŷ ∞ becomes too large again, and so forth.

Drury-Oja flow.
There is a close connection between the Riccati flow in Corollary 1 and the continuous orthogonalization method of Humpherys and Zumbrun [49], proved by direct comparison.
Humpherys and Zumbrun [49] derive this flow by a QR-decomposition of the solution Y = QR ∈ V(n, k) to the linear, spectral, globally bounded flow Y ′ = A(x) Y .We can think of this continuous orthogonalization method as generating an approximate flow on the Grassmann manifold whilst evolving the coordinatization (see Edelman, Arias and Smith [33]; Bindel, Demmel and Friedman [11]).It is also known as a Drury-Oja flow (see Drury [32], Oja [80], Yan, Helmke and Moore [98], Bridges and Reich [15], Dieci and Van Vleck [30] and Hairer, Lubich and Wanner [42, p. 136]).The determinant of the upper triangular matrix, det R, grows exponentially.Thus, since it is nonzero in the far field, we know it remains nonzero in the whole integration interval R. Consequently, Q = Y R −1 is globally finite on R, i.e. there are no singularities in the Drury-Oja flow.
There is a natural su(n) Lie algebra action on the Stiefel manifold of orthonormal k-frames, V 0 (n, k), generated by the map (σ, We could use this to construct a numerical method that preserves orthonormality of Q.

Grassmann Gaussian elimination method
5.1.Algorithm.The Grassmann Gaussian elimination method using a Lie group method on the Stiefel manifold (GGEM-LG), proceeds as follows: (1) Suppose initially we are given data y i • m (x m ) in the coordinate patch U im .(2) Across the integration interval [x m , x m+1 ], compute σ m using the Magnus expansion-we recommend the fourth order Magnus expansion.4) Apply quasi-optimal Gaussian elimination (QOGE) (outlined below) to Y m+1 ; this generates the solution y i • m+1 (x m+1 ) in the coordinate chart U im+1 , and the rank k matrix U m+1 we have effectively equivalenced by.
Across the integration interval [x m , x m+1 ], we could generate Y m+1 from y i • m (x m ) by solving the flow on the Stiefel manifold V(n, k) using a classical Runge-Kutta step (GGEM-RK).Both algorithms are summarized in the following diagram.
How much of U m+1 do we retain at each step?We address this in Section 6.4.

5.2.
Quasi-optimal Gaussian elimination.To project the endpoint solution Y m+1 in the Stiefel manifold onto a quasi-optimal coordinate patch, say U im+1 , we use a Gauss-Jordan approach.This entails using elementary column operations and optimal pivoting, until a specified k ×k submatrix of Y m+1 becomes the identity matrix (we naturally assume k ≥ 2).To be more explicit about what we mean by optimal pivoting, we outline the procedure: (1) Look for the largest term in magnitude in Y (1) := Y m+1 .Nominate this term as the pivot term and suppose this occurs in row i 1 , column j 1 .Use elementary column operations with this term, i1,j1 , as the pivot to render all the other (k − 1) terms in that row equal to zero.Finally, use the elementary column operation of scalar multiplication to render the pivot term itself equal to one.Let us call the resulting n × k matrix Y (2) ; we can write (see for example Meyer [73]) (1) , where the elementary column operations performed are encoded in the elementary matrix U (1) .To be precise, we can write where for ℓ ∈ {1, . . ., k}/{j 1 } the elementary matrix E (1) ℓ encodes the following elementary column operation on column c ℓ : For ℓ = j 1 the elementary matrix E (1) j1 encodes the elementary column operation c j1 → c j1 /p i1j1 .In summary, the matrix Y (2) has value one at position (i 1 , j 1 ) and otherwise zeros in row i 1 ; and we have det (2) In the (n − 1) × (k − 1) submatrix of Y (2) identified by excluding row i 1 and column j 1 , look for the largest term in magnitude.Again nominate this term as the pivot term and suppose this occurs in row i 2 (which will be distinct from i 1 ), and column j 2 (which is distinct from j 1 ).Here i 2 and j 2 refer to the row and column relative to the original n × k matrix Y (2) .Use elementary column operations with this term, i2,j2 , as the pivot to render the terms in row i 2 and columns {1, . . ., k}/{j 1 , j 2 } of Y (2) equal to zero.Again use the elementary column operation of scalar multiplication to render the pivot term itself equal to one.Let us call the resulting n × k matrix Y (3) ; we can write (2) , where the elementary column operations performed on Y (2) are encoded in the elementary matrix U (2) .Again, to be precise, we can write where for ℓ ∈ {1, . . ., k}/{j 1 , j 2 } the elementary matrix E (2) ℓ encodes the following elementary column operation on column c ℓ : Note that in our expression for U (2) above we could have either j 1 < j 2 or j2 encodes the elementary column operation c j2 → c j2 /p i2j2 .In summary, the resulting matrix Y (3) has value one at positions (i 1 , j 1 ) and (i 2 , j 2 ) and otherwise zeros in row i 1 , and zeros in row i 2 in columns {1, . . ., k}/{j 1 , j 2 }; and we have det U (2) = (p i2j2 ) −1 .
(3) Continue this process.Focus on the (n − 2) × (k − 2) submatrix of Y (3)  identified by excluding the rows i 1 , i 2 and columns j 1 , j 2 ; look for the largest term in magnitude.Nominate this term as the pivot term-suppose it occurs in row i 3 and column j 3 , relative to the original n × k matrix Y (3) , and so forth.On completing the final kth step in this process, we will have The final n × k matrix Y (k) will have ones in positions (i ℓ , j ℓ ), for ℓ = 1, . . ., k.In row i ℓ it will have zeros in columns {1, . . ., k}/{j 1 , . . ., j ℓ }.We set (4) Perform column swaps in Y (k) so that column j 1 becomes column 1, column j 2 becomes column 2, and so forth so that finally column j k is forced to be column k.These column swaps can be encoded in the elementary matrix Σ with det Σ = (−1) #{swaps} .The resulting matrix Ỹ (k) is given by The i m+1 × {1, . . ., k} submatrix of Ỹ (k) given by is lower triangular with ones on the diagonal.Finally we set In practice of course, we do not compute L −1 , but continue performing elementary column operations on Ỹ (k) so that the submatrix L becomes the identity matrix, thus generating y i • m+1 .Hence we have effectively performed the decomposition Note that if Y m+1 ∈ V(n, k) then the entries in y i • m+1 and U m+1 , produced as a result of this process, will all be finite.Further if Y m+1 depends analytically on a parameter, then the product y i • m+1 • U m+1 , naturally does as well.Lastly we remark that we could have performed alternative elementary column operations of the form i,ℓ c j (we do not include the scalar multipication operations here) with the result that det

5.3.
Complexity.The complexity of the quasi-optimal Gaussian elimination algorithm, dominated by the search for the largest elements in the successively decreasing submatrices of Y m+1 , is of order nk 2 .The method is a practical approach to maximize the determinant of the k × k submatrix removed from Y m+1 .It will not in general choose the submatrix with the largest determinant-hence the label quasi-optimal.This could be achieved by searching through all the k × k submatrices of Y m+1 , i.e. all the Plücker coordinates, but this has complexity of order n choose k.An interesting question here is whether there is an efficient way to use the Plücker relations to reduce this complexity?6. Spectral problems 6.1.Linear Stiefel flow.Consider the linear spectral problem on R: We assume there exists a subdomain Ω ⊆ C containing the right-half complex plane, such that for λ ∈ Ω there exists exponential dichotomies on R − and R + with the same Morse index k in each case (see Henry [43] and Sandstede [88]).Let Y − (x; λ) ∈ V(n, k) denote the matrix whose columns are solutions to the spectral problem and which span the unstable manifold section at x ∈ [−∞, +∞).Let Y + (x; λ) ∈ V(n, n − k) denote the matrix whose columns are the solutions which span the stable manifold section at x ∈ (−∞, +∞].

6.2.
Matching.The values of spectral parameter λ ∈ Ω for which the columns of Y − and columns of Y + are linearly dependent on R are pure-point eigenvalues.The Evans function D(λ) is the measure of the degree linear dependence between the two basis sets Y − and Y + , i.e. of the degree of transversal intersection between the unstable and stable manifolds (see Alexander, Gardner and Jones [2]; Nii [79]): It is analytic in Ω.In practice we drop the non-zero, scalar exponential prefactor and evaluate the Evans function at a matching point x * .
There are other matching criteria measuring the degree of intersection between subspaces that do not use the determinant: for example computing the angle between subspaces as suggested by Björck and Golub [13] or computing the smallest eigenvalue as suggested by Hutson [51] and Ixaru [55].Both these latter techniques might be important for large systems when computing the determinant could be an unstable process, indeed, we investigate them in this context in Ledoux, Malham, Niesen and Thümmler [66].However in both cases the magnitude of a function of the spectral parameter is computed, whose zeros correspond to eigenvalues.Hence we must search for touchdowns to zero in the complex parameter spectral plane which can be problematic.For the examples we consider here, which are not too large, using the determinant suffices.6.3.Initialization.We construct the n×k matrix Y − 0 (λ) whose columns are the k eigenvectors of A(−∞; λ) corresponding to eigenvalues with a positive real part (see Humpherys and Zumbrun [49] and also Humpherys, Sandstede and Zumbrun [50] for how to preserve analyticity with respect to the spectral parameter λ).In practice, integration starts at x = ℓ − for some suitable, usually negative, value of ℓ − .Analogously we construct the n × (n − k) matrix Y + 0 (λ) whose columns are the n − k eigenvectors of A(+∞; λ) corresponding to eigenvalues with a negative real part.Again, in practice, we integrate backwards from x = ℓ + for some suitable, usually large and positive, value of ℓ + .6.4.GGEM matching and analyticity.To compute Y ± (x * ; λ) we start with Y ± 0 (λ) at x = ℓ ± and integrate centrally towards x = x * .Our goal in this section is to show that we only need the determinant of the rank k transformations in the GGEM method described at the beginning of Section 5, to retain analyticity for the Evans function in Ω.Since the procedure is the same in both intervals [ℓ ± , x * ] we will describe it for the generic interval [ℓ, x * ] for Y (x; λ) ∈ V(n, k) starting with value Y 0 (ℓ; λ) at x = ℓ.Suppose we use M successive computation subintervals [x m , x m+1 ] in [ℓ, x * ] where x m := ℓ + m (x * − ℓ)/M .We label the nodal solution values at x = x m as Y m (λ).
Let S m,m+1 (λ) denote an approximation to the flow-map across [x m , x m+1 ] to the linear system Y ′ = A(x; λ) Y .We assume that S m,m+1 (λ) preserves analytic dependency on λ, so that the next step solution value Y m+1 (λ) = S m,m+1 (λ) Y m (λ) analytically depends on λ if Y m (λ) does.For example, in the case of GGEM-LG then S m,m+1 (λ) = exp σ m and most straightforward Magnus based integrators will naturally preserve analyticity with respect to λ. Similarly most simple Runge-Kutta methods used to generate S m,m+1 (λ), or directly generate the next step solution value Y m+1 (λ), will preserve analyticity.
Our numerical procedure would proceed as follows.Across [x 0 , x 1 ] we have where in the last step we applied QOGE to the n × k matrix S 0,1 (λ) where we applied QOGE to S 1,2 (λ) y i • 1 (x 0 ; λ).Repeating this argument across the subsequent intervals [x m , x m+1 ] for m = 2, . . ., M − 1, we get the following approximation Ŷ (x * ; λ) to Y (x * ; λ): Returning to the separate interval calculations on [ℓ ± , x * ], we see that by the procedure just outlined we can generate the solution approximations The number of integration steps M can of course be different in each interval.Hence, after dropping the exponential prefactor and fixing the matching point to be x * , the Evans function D(λ; x * ) can be approximated by D(λ; x * ) where This is an analytic function of λ.Indeed as we hinted previously, at each computation step x = x m we need only store y i • m (x m ; λ) and the value which gets updated at each step by simply multiplying the previous step value by the complex scalar determinental factor for current step.Hence to preserve analyticity for the Evans function using GGEM we must generate an approximate flow on Gr(n, k) ⊗ C.
6.5.Scaled GGEM.The scalar determinental factor just described, that we update at each step, grows exponentially.This would be tempered by the scalar exponential prefactor in the definition of the Evans function.An accurate practical procedure here is as follows (to be applied with due care).When integrating in the interval [ℓ − , x * ], at each step, divide the scalar determinental factor in GGEM by exp (µ , where h is the stepsize, and the µ − i (λ) are the (spatial) eigenvalues, with positive real part, of A(−∞; λ).When integrating in the interval [ℓ + , x * ], at each step, divide the scalar factor by exp h , where the µ + i (λ) are the eigenvalues, with negative real part, of A(+∞; λ).To see that this normalization is appropriate, we recall the Plücker coordinates of Section 2. After applying the optimal Gaussian elimination algorithm to Y m+1 the i • m+1 th row elements of y i • m+1 are themselves (n − k)k Plücker coordinates; normalized by det U m+1 .The i • m+1 th row elements of y i • m+1 and det U m+1 , can be used to reconstruct the remaining Plücker coordinates through the homogeneous, quadratic Plücker relations.Hence the Plücker coordinates, or complete set of k × k minors, of Y m+1 and y i • m+1 differ by a factor det U m+1 .It is well known that if the original vector field on V(n, k) is linear, then the Plücker coordinates corresponding to Y m+1 satisfy a (larger) linear system of equations.In the left far field the coordinates thus grow exponentially, in fact with growth rate µ − 1 (λ) + • • • + µ − k (λ); hence our recommendation to divide by the exponential factor suggested (with an analogous argument for the right far field).See Alexander, Gardner and Jones [2], Alexander and Sachs [3], Brin [16,17] or Allen and Bridges [5] for more details.

Applications
We present some numerical results for three different applications.The three applications reduce to the solution of a system showing multiple distinct exponential growth and decay rates in the stable and unstable subspaces, respectively.We show that our approach resolves this numerical obstacle successfully and can compete with the continuous orthogonalization method of Humpherys and Zumbrun [49].7.1.Algorithms.We implement six different algorithms as follows.
(1) Riccati-RK: Riccati method with fixed coordinatization with the flow of the Riccati vector field approximated by the classical fourth order Runge-Kutta method.We generically chose the coordinatization labelled by i − = {1, . . ., k} and i + = {k + 1, . . ., n} for the left-hand and right-hand intervals, respectively.Hence if Y − 0 (λ) and Y + 0 (λ) denote the unstable and stable subspaces of A(−∞; λ) and A(+∞; λ), respectively, then we set where (Y ) i,i ′ denotes the i × i ′ submatrix of Y .We integrate the Riccati equation outlined in Corollary 1 in the two intervals and evaluate the modified Evans function Provided neither Riccati flow becomes singular, this Evans function is analytic in the spectral parameter λ.
(2) Möbius-Magnus: Uses the Schiff and Shnider approach to integrate through singularities, combined with a Lie group method to advance the solution on the general linear group, as described at the end of Section 3. The same generic fixed coordinate charts are used as for the Riccati-RK method above.Over the integration interval [x m , x m+1 ] with an equidistant mesh stepsize h, we advance the solution on the Lie algebra using the fourth order Magnus method with the two Gauss-Legendre points (see Iserles, Marthinsen and Nørsett [52]) x [1] We then compute the Möbius map ŷm+1 = µ ŷm • exp σ m to advance the solution in the fixed Grassmannian chart-for the left-hand interval i − = {1, . . ., k} while for the right-hand interval i + = {k + 1, . . ., n}.We evaluate the same Evans function as for the Riccati-RK method above.
(3) GGEM-RK: Scaled Grassmann Gaussian elimination method, with the classical fourth order Runge-Kutta method used to advance the solution on the Stiefel manifold, as described in Sections 5 and 6.5.We evaluate the Evans function D(λ; x * ) in Section 6.
(4) GGEM-LG: Same as GGEM-RK but with a fourth order Magnus method used to advance the solution on the Stiefel manifold instead, i.e.Y m+1 = exp(σ m ) y im where σ m is generated as for the Möbius-Magnus method above.
(5) Riccati-QOGE: Riccati method with coordinate swapping as described in Section 4.2.We have chosen to implement the method in the following form.At each integration step we advance the solution on the Stiefel manifold using the Magnus method (we could also use a Runge-Kutta method here).We apply elementary column operations to the resulting solution matrix to convert the pre-determined rows indexed by i from the previous step to the identity matrix.Then if ŷ ∞ is less than or equal to a tolerance size, we keep this index i for the next step.If it is greater, we apply QOGE at the end of the next step after advancing the solution on the Stiefel manifold, thus generating a new index.As for GGEM-RK and GGEM-LG, we update the scalar determinental factor at each step (produced by the elementary column operations with the pre-determined index or QOGE).We divide the scalar determinental factor by the scalar exponential factors, as described for the scaled GGEM method.We evaluate the same Evans function also.
(6) CO-RK: Continuous orthogonalization method of Humpherys and Zumbrun with the classical fourth order Runge-Kutta method used to advance the solution on the Stiefel manifold of orthonormal frames.The initial conditions Q ± 0 (λ) for the Q-matrices are obtained by QR-factorization of Y ± 0 (λ).From Humpherys and Zumbrun [49], to ensure analyticity we must also solve the scalar problems (det x det R ± .The Evans function is then given by

Boussinesq system.
As the first test system, we consider the Boussinesq system studied by Humpherys and Zumbrun [49].The (good) Boussinesq equation, expressed in a co-moving frame moving to the right with wave speed c, is given by It has solitary wave solutions of the form ū(x) ≡ 3  2 x , where |c| < 1.These waves are stable when 1/2 < |c| < 1 and unstable when |c| < 1/2.
If we consider small perturbations about the travelling wave ū we generate a linear spectral problem of the form Y ′ = A(x; λ)Y , where When the spectral parameter λ lies in the right-half complex plane the eigenvalues of A(±∞; λ) spectrally separate into two growth and two decay modes, i.e. k = 2.We used ℓ ± = ±8 in our experiments.
In Figure 1 we show the Evans function computed along the real axis from λ = 0 to λ = 0.2 for the unstable pulse with c = 0.4.The Riccati-RK (left plot) and CO-RK (right plot) methods detect a zero of the Evans function near λ = 0.155, indicating an unstable eigenvalue there.An accurate value of the eigenvalue can be found by using a standard root-finding method.This yields the values in Table 1.The Riccati-RK and the CO-RK methods both converge to the same eigenvalue when the number of steps N increases.As a check, the Matlab ode45 solver was used with a relative tolerance 10 −8 and absolute tolerance 10 −10 to integrate the systems, leading to the same resulting eigenvalue: λ = 0.15543141 for both methods.
Table 1.Zero of the Evans function for the Boussinesq problem computed with the Riccati-RK method (with fixed coordinate patches identified by i = − {1, 2} over [−8, 0] and i + = {3, 4} over [0, 8]) and the CO-RK method.Here N is the number of (equidistant) steps used in the mesh.

N
Riccati  Function evaluation for the Riccati vector field requires three matrix-matrix multiplications.This is the same number of matrix-matrix multiplications needed to evaluate the Drury-Oja vector field.However, the matrices in the Drury-Oja vector field are n × k and n × (n − k), respectively, while the matrices in the Riccati vector fields have smaller dimension: (n − k) × k.Because of the smaller dimension of the systems to be integrated, our Riccati approach is faster than the continuous orthogonalization problem.For example, to construct Figure 1 the Evans function was evaluated at 200 distinct λ values between λ = 0 and λ = 0.2.Using the fourth-order Runge-Kutta method with N = 512 steps, this required 33 seconds for the CO-RK method, while the Riccati-RK method needed 24 seconds (Matlabimplementation, CPU 2.4GHz).
In Figure 2 we compare the error in the eigenvalue and efficiency of computation for all six methods, when we match at x * = 0. We see that the methods that use the Magnus expansion to advance the solution on the Stiefel manifold are the most accurate for a given stepsize.They are also the most efficient, delivering the best accuracy for given computational effort.The Riccati-RK method does not suffer from singularities for the chosen fixed patches when matching at x * = 0, at least for the range of values of the spectral parameter in the vicinity of the eigenvalue (as well as the origin and anywhere inbetween).However, if we change the matching point to x * = ℓ + = +8 there are singularities in the Riccati-RK solution (as a result of a vanishing determinant of u − ).In particular, a singularity appears around x = 2 for λ equal to the eigenvalue (see Figure 4).Hence we compare the remaining five methods in Figure 3 in this case.We see that using the Möbius-Magnus method to integrate through a singularity does not introduce loss in accuracy.The resulting Evans function can have poles, as seen in Figure 4, which appear at λ-values where the Riccati equation has a singularity at the matching point.This means that in  some cases the matching point should be chosen rather carefully in order not to have the poles interfering with the eigenvalue(s).When applying GGEM-RK the Evans function is analytic and the choice of the matching point is less important.Figure 5 shows the Evans function obtained when the GGEM-RK evolves y i from ℓ − = −8 to ℓ + = +8.To construct the plot in Figure 5, the quasi-optimal Gaussian elimination process was applied at each step in the integration.However it is clear from the right plot in Figure 5, that multiple successive steps can be integrated in the same coordinate patch.For example, between x = −8 and x = −3 the coordinate patch does not change.Performing the whole quasi-optimal Gaussian elimination process only when a certain criterion is satisfied, reduces the computing time.Using the Riccati-QOGE method, we change the coordinatization when ŷ− ∞ > 2. This generates an Evans function very similar to that in Figure 5.As seen in Figure 6 the quasi-optimal Gaussian elimination process is then performed only two times for λ equal to the eigenvalue.We compare the error in computing the eigenvalue for different choices of matching point-in fact for x * anywhere in the interval [−8, 8]-for all six methods in Figure 7.We see that the most accurate and robust methods are the GGEM-LG and Riccati-QOGE methods.Some methods, such as the Riccati-RK method as discussed already, break down when singularities impinge on the matching pointthe singularities in the Evans function are observed in Figure 8. Generally we also see in Figure 7 that the GGEM-LG and Riccati-QOGE methods outperform the CO-RK method in terms of accuracy.
Overall, we observe in this example that when computing the eigenvalue, those methods based on the Magnus expansion are superior in accuracy and efficiency.Note that for GGEM-RK, the quasi-optimal Gaussian elimination process is an additional nk 2 operation.However, to ensure analyticity for the CO-RK method, there are two additional matrix-matrix multiplications in the equations for det R ± (operational cost kn 2 ) required at each step.7.3.Autocatalytic fronts.As a second example, we study travelling waves in a model of autocatalysis in an infinitely extended medium Here u(x, t) is the concentration of the reactant and v(x, t) is the concentration of the autocatalyst.We suppose (u, v) approaches the stable homogeneous steady state (0, 1) as x → −∞, and the unstable homogeneous steady state (1, 0) as x → +∞.The diffusion parameter δ is the ratio of the diffusivity of the reactant to that of the autocatalyst and m is the order of the autocatalytic reaction.The speed of the co-moving reference frame is c.The system is globally well-posed for smooth initial data and any finite δ > 0 and m ≥ 1.
From Billingham and Needham [9] we know that a unique heteroclinic connection between the unstable and stable homogeneous steady states exists for wavespeeds c ≥ c min .The unique travelling wave for c = c min converges exponentially to the homogeneous steady states and is computed by a simple shooting algorithm (see Balmforth, Craster and Malham [8]).The resulting travelling wave for δ = 0.1 and m = 9 is shown in Figure 9.The stability of the travelling wave of velocity c can be deduced from the location of the spectrum of the eigenvalue problem Y ′ = A(x; λ)Y , where where ū and v represent the travelling wave solution.
The pulsating instability occurs when δ < 1 is sufficiently small and m is sufficiently large (see Metcalf, Merkin and Scott [72] and Balmforth, Craster and Malham [8]).For δ fixed and m increasing, a complex conjugate pair of eigenvalues crosses into the right-half λ-plane signifying the onset of instability via a Hopf bifurcation.Figure 10 shows the onset of this instability for δ = 0.1 as m is increased from 8 to 9 (see Aparicio, Malham and Oliver [7]).The figure shows the zero contour lines of the real and imaginary parts of the Evans function.Solid lines correspond to zero contours of the real part of D(λ), dashed lines to the imaginary part of D(λ).We see that a complex-conjugate pair of eigenvalues crosses into the right-half plane, indicating the onset of instability.Figure 10 was constructed using the Riccati-RK method with the fixed coordinate patches identified by i − = {1, 2} from x = −10 to x * = −7, and i + = {3, 4} from x = 10 to x * = −7.The matching point x * = −7 is chosen roughly centred on the wavefront.
We compare the Riccati-RK, Möbius-Magnus, CO-RK and GGEM-LG methods in Figure 11 where we plot the absolute error in the eigenvalue vs the stepsize (upper panel) and also vs cputime (lower panel).The eigenvalue in question is that in the first quadrant in Figure 10 for δ = 0.1 and m = 9. Figure 11 was generated as follows.Starting with an initial guess lying within a small square around the eigenvalue, we iterated a standard root finding algorithm until we arrived in a square (containing the eigenvalue) which was smaller than a preset tolerance.We see in Figure 11 that the Riccati-RK method produces a slightly better error for a given stepsize, and is marginally more efficient than the GGEM-LG method.The CO-RK method produces a larger error for a given computational effort.This is not surprising, as again, the matrices in the Drury-Oja vector field are twice as big (4 × 2) as the ones in the Riccati vector fields (2 × 2).
When we match at x * = −7, there is little to distinguish the Riccati-RK, Möbius-Magnus, CO-RK and GGEM-LG methods.We compare all four methods for different matching positions x * in Figure 12, which was generated using the same root finding criteria as for Figure 11, except all the methods used N = 256 steps.Note that the errors in the Möbius-Magnus and GGEM-LG methods are uniform for any matching values in the range [ −10, 10].The CO-RK method error doesn't vary that much either and is slightly larger.Note that no values are plotted for the CO-RK method at the matching points x * = 8, 10.In these cases the classical Runge-Kutta method applied to the Drury-Oja vector field is unstable for N = 256 steps for some λ-values close to the eigenvalue.This problem is resolved by increasing the number of steps to N = 512.For a range of matching positions roughly in [−10, 0], there are no singularities of the Riccati-RK method in the left and right-hand integration intervals for values of the spectral parameter λ close to the eigenvalue.Indeed for this range of matching positions the Riccati-RK method delivers the best accuracy.However for matching positions outside this range, for values of the spectral parameter λ not far from the eigenvalue, the Riccati-RK solution does have a singularity for some matching points (which we can see in the contour plots in Figure 13).This makes the eigenvalue-searching algorithm failindicated by no error points for those matching position values.We also do not show the error for the Riccati-RK method for the matching points x * = −8, −10, as  We show three different closed contours in the first quadrant of the complex λ-plane in each of the left panels.In the top two left panels, the contour encloses the eigenvalue in that quadrant (found in Figure 13 for δ = 0.1, m = 9).In the corresponding three panels on the right we show how the argument (in multiples of 2π) of the Evans function D(λ) changes as we perform a complete circuit of the contour.The Evans function was computed using the GGEM-LG method matching at x * = +14.The number of patch changes performed for each fixed λ value are indicated.
there are singularities in the Evans function close to the real axis for these matching points (again see Figure 13).This means that we cannot for example, apply the argument principle in the first quadrant, though starting sufficiently close to the eigenvalue we can still use the Riccati-RK method as part of a root-finding algorithm to determine the eigenvalue.
Figure 13 shows the contour lines of |D(λ)| for δ = 0.1 and m = 9, close to the eigenvalue in the first quadrant, when using the Riccati-RK, Möbius-Magnus, CO-RK and GGEM-LG methods, respectively, and matching at three positions x * = −8, 0, +8.We see that the CO-RK and GGEM-LG methods show the least sensitivity to the choice of matching position and produce smooth contour plots for all three matching points.The contour plots across shape and scale look very similar for both these methods.By contrast the Riccati-RK and Möbius-Magnus methods appear to develop singularities close to the eigenvalue when the matching position x * is 0 or +8.
Lastly in Figure 14 we demonstrate the argument principle for counting zeros of the Evans function inside closed contours.We computed the Evans function using the GGEM-LG method and matched at x * = +14.As expected, if the closed contour in the complex λ-plane encloses the eigenvalue, then the change in the argument of the Evans function around the complete contour is one (once we have accounted for the 2π factor in the argument principle).We also show, for each fixed λ value, the number of patch changes that occured as we integrated from x = −14 through to x = +14.7.4.Ekman boundary layer.The third test system is a boundary layer flow over a flat plate which is infinitely extended in the x and y direction and rotates around the half infinite z-axis with a given rotational speed.Linear stability of the Ekman boundary layer has been investigated in Allen and Bridges [6] and Allen [4] using the compound matrix method.The flow is governed by the continuity equation u x + v y + w z = 0, and the Navier-Stokes equations in a co-rotating frame Here R e , R o and E k denote the Reynolds, Rossby and Ekman numbers, respectively.
After non-dimensionalization and setting R e = R o , E k = 1, the linear stability of the boundary layer is determined by the eigenvalues λ of the linear problem Y ′ = A(z; λ)Y , where (see Allen [4, p. 176 Here the parameters γ and ǫ represent the radial and angle components, respectively, of a polar coordinate parameterization of horizontal wavenumbers associated with the x and y directions-see Allen and Bridges [6] for more details.We choose the fixed coordinate patch identified by i + = {1, 2, 3} to integrate the corresponding Riccati equation from z = ℓ + to z * = 0.The boundary condition for  the rigid wall at z * = 0 as given in Allen and Bridges is We thus compute the Evans function: We computed neutral curves, i.e. curves in the ǫ-γ plane where Re(λ) = 0, using the Riccati-RK method with i + = {1, 2, 3} to compute ŷ(0; λ) ∈ C 3×3 and consequently the Evans function D(λ; z * ).For continuation of the curves we used the Matlab package MatCont which uses pseudo-arclength continuation (Dhooge, Govaerts and Kuznetsov [29]).Figures 15 and 16 show the neutral curves which match those in Allen and Bridges [6] and Allen [4].The integration of the Riccati system has been done with the Matlab ODE-solver ode23s from z = 10 to z * = 0 (as in Allen and Bridges) with absolute and relative tolerances 10 −6 and 10 −4 .The stable subspace of A(+∞; λ) was constructed using the Matlab eigenvalue-solver eig (we also used direct formulae for the eigenvectors to construct analytic bases for the stable subspace but this did not significantly change the overall performance).
For comparison we also implemented the CO-RK and GGEM-LG methods.We tested the performance of all three methods, in each case evaluating the Evans function on a 20 × 20 grid for λ in the complex plane.In Figures 17 we present contour plots of |D(λ; 0)|, and see that all methods find a root at λ ≈ 0.002−0.117i.The computation times for a 2.4Ghz machine were: 134 seconds for Riccati-RK, 150 seconds for GGEM-LG and 189 seconds for CO-RK.As a comprehensive check, we also implemented the compound matrix method (i.e.Plücker coordinates, of which there are 20), described in Allen and Bridges, for this performance test.As expected, since this method involves integrating a linear system of order 20, it was an order of magnitude slower (while giving the same results).

Concluding discussion
We have shown that the new scaled Grassmann Gaussian elimination method as well as the Riccati method with quasi-optimal patch swapping, compete with the continuous orthogonalization method for computing the Evans function.Both new methods deliver superior accuracy for the same computational cost when combined with Lie group Magnus integration to advance the solution.Moreover, as hoped, numerically these new methods appear to be robust in the sense that they are insensitive to the choice of the matching position in the computational domain.We now outline several directions in which we plan to use and extend these methods.
One of the main goals we have had in mind in this paper is that of large scale spectral problems, in particular the stability of travelling waves with a multi-dimensional structure.There is recent research extending the Evans function approach in this direction-see Deng and Nii [27], Gesztesy, Latushkin and Makarov [36] and Gesztesy, Latushkin and Zumbrun [37].From a numerical perspective we have, together with Niesen, implemented some of the methods we propose in this paper in a multi-dimensional context.In particular, it is well known in autocatalysis and combustion that planar travelling fronts can be unstable to transverse perturbations and develop into steadily propagating travelling fronts with wrinkles.In Ledoux, Malham, Niesen and Thümmler [66] we show that the wrinkled fronts themselves develop an instability as a diffusion parameter is further increased.
For large scale problems the Lie group methods we propose using the Magnus expansion may become prohibitive.This is because of the effort required to compute the matrix exponential-see Moler and Van Loan [75], Celledoni and Iserles [19], Munthe-Kaas and Zanna [78] and Iserles and Zanna [53].For the examples we considered this was not an issue.However it remains to be seen if such Lie group methods will be cost effective for larger problems-the methods we proposed based on Runge-Kutta integration such as GGEM-RK can be used as they scale favourably with system size.
The constructs and Grassmannian reductions we have considered in this paper, it turns out, have their origins in the control theory literature dating back to the early seventies, in particular in the pioneering papers of Hermann and Martin [44,45,46,47,48], Martin and Hermann [71] and Brockett and Byrnes [18].We also found Bittanti, Laub and Willems [12], Lafortune and Winternitz [63], Rosenthal [87], Shayman [92] and Zelikin [100] particularly useful resources.A future direction we would like to explore is whether there are any applications of the numerical methods we have outlined here to practical non-autonomous control problems?
Riccati methods in particular also have their origins in the quantum chemistry literature also dating back to the early seventies-a recent survey of these numerical methods can be found in Chou and Wyatt [24].However also see Light and Walker [67], Johnson [56], Hutson [51] and Gray and Manopoulous [39].In particular the log-derivative and R-propagation methods correspond to special choices of Grassmannian patch in the Riccati methods we mention above.Prüfer methods, for which we can think of the patch evolving, originate even further back; see Prüfer [84] and Pryce [85].
Of course, our quasi-optimal Gaussian elimination process for choosing a suitable representative patch was inspired by the Schubert cell decomposition of the Grassmann manifold; see for example Billey [10], Griffiths and Harris [40], Kleiman and Laksov [58], Kresch [60], Postnikov [83], Sottile [93] and, in a somewhat different vein, Kodama [59].Since the Grassmann manifold is the disjoint union of Schubert cells, the question is, can we express the flow on the Grassmann manifold as a flow on Schubert cells (see Griffiths and Harris and also Ravi, Rosenthal and Wang [86])?Can we construct the corresponding flow on the cohomological ring of Schubert cycles (Chern [22]; Fulton [35])?Ernst Hairer for inviting us and providing so much support and enthusiasm.We are also indebted to the facilities at the Isaac Newton Institute which were invaluable.We also thank Tom Bridges, Jitse Niesen, Jacques Vanneste and Antonella Zanna for stimulating discussions on this work.Veerle Ledoux is a postdoctoral fellow of the Fund of Scientific Research-Flanders (F.W.O.-Vlaanderen).Vera Thümmler was supported by CRC 701: Spectral Structures and Topological Methods in Mathematics.

Figure 1 .
Figure 1.The Evans function of the Boussinesq system for the unstable pulse having wave speed c = 0.4.The left plot is generated by the Riccati-RK method with fixed coordinate patches identified by i − = {1, 2} over [−8, 0] and i + = {3, 4} over [0, 8].The right plot shows the result of the CO-RK method.

Figure 4 .Figure 5 .
Figure 4.The Evans function when the Möbius-Magnus method is applied (left panel) with the matching point as x * = +8.The entries of the solution ŷ− , passing through the singularity when integrating from −8 to 8, for λ = 0.15543141 (right panel).

Figure 6 .Figure 7 .Figure 8 .
Figure 6.The Evans function when the Riccati-QOGE method is applied (left panel).The entries of y i • are shown for λ = 0.15543141 (right panel).The criterion used for swapping to a new coordinate patch was ŷ− ∞ > 2.

Figure 10 .
Figure 10.Zero contour lines of the real (solid) and imaginary (dashed) parts of the Evans function for the autocatalysis problem with δ = 0.1 and m = 8 (left panel) and m = 9 (right panel).

Figure 11 .
Figure 11.Error in the eigenvalue in the first quadrant when δ = 0.1 and m = 9 for different methods, vs stepsize (upper panel) and vs cputime (lower panel), matching at x * = 0.

Figure 12 .
Figure 12.Error in the eigenvalue in the first quadrant when δ = 0.1 and m = 9, for different methods, for different matching points.The number of steps in the equidistant mesh is N = 256.

Figure 14 .
Figure 14.We show three different closed contours in the first quadrant of the complex λ-plane in each of the left panels.In the top two left panels, the contour encloses the eigenvalue in that quadrant (found in Figure13for δ = 0.1, m = 9).In the corresponding three panels on the right we show how the argument (in multiples of 2π) of the Evans function D(λ) changes as we perform a complete circuit of the contour.The Evans function was computed using the GGEM-LG method matching at x * = +14.The number of patch changes performed for each fixed λ value are indicated.

Figure 15 .
Figure 15.Neutral curves for rigid wall, R e fixed (values indicated).
and Stasheff [74, p. 56].2.2.Representation.Following the exposition in Griffiths and Harris [40], any k-plane in C n can be represented by an n × k matrix of rank k, say Y ∈ C n×k .Any two such matrices Y and Y ′ represent the same k-plane element of Gr(n, k) if and only if Y ′ = Y u for some u ∈ GL(k) (the k-dimensional subspace elements are invariant to rank k closed transformations mapping k-planes to k-planes).Let i = {i 1 , . . ., i k } ⊂ {1, . . ., n} denote a multi-index of cardinality k.Let Y i • ⊂ C n denote the (n − k)-plane in C n spanned by the vectors {e j : j ∈ i} and for more details.If we change the basis, the basis for the image changes by the determinant of the transformation matrix.Hence the map is a point in P k C n .We can recover Y from its image Y 1 ∧ . . .∧ Y k as the set of all vectors v such that v∧Y 1 ∧. ..∧Y k = 0. Further, a point of P k C n is in the image of p if and only if its representation as a linear combination of the basis elements of k C n , consisting of all possible distinct wedge products of a k-dimensional basis in C n , is completely decomposable.Hence the image of p is a subvariety of P k C n of completely decomposable elements.It can also be realized as follows.
GL(k).Let a, b, c and d denote the i × i, i × i • , i • × i and i • × i • submatrices of A, respectively.
k) where A ∈ gl(n).Fixing a coordinate patch U i for Gr(n, k) for some i = {i 1 , ..., i k } we can always decomposeY ∈ V(n, k) into Y = y i • u,where u ∈