Collocation Methods for Exploring Perturbations in Linear Stability Analysis

Eigenvalue analysis is a well-established tool for stability analysis of dynamical systems. However, there are situations where eigenvalues miss some important features of physical models. For example, in models of incompressible fluid dynamics, there are examples where linear stability analysis predicts stability but transient simulations exhibit significant growth of infinitesimal perturbations. This behavior can be predicted by pseudo-spectral analysis. In this study, we show that an approach similar to pseudo-spectral analysis can be performed inexpensively using stochastic collocation methods and the results can be used to provide quantitative information about instability. In addition, we demonstrate that the results of the perturbation analysis provide insight into the behavior of unsteady flow simulations.

1. Introduction. This study is concerned with a refined understanding of the classic problem of stability of dynamical systems. Let Let γ = γ(x, 0) represent a small perturbation of u (s) . Suppose the perturbed quantityû(x, 0) := u (s) (x) + γ(x, 0) is taken as an initial condition for (1.1), for which integration leads to a solutionû(x, t) = u (s) (x) + γ(x, t). Ifû(x, t) reverts to u (s) (x) (γ(x, t) → 0) as t increases, then the steady solution is said to be stable; otherwise it is unstable. In typical applications, f (u, t) = f α (u, t) depends on a parameter α, as does the resulting steady solution u (s) α , and we are interested in the set of values of α for which u (s) α is stable. Spatial discretization of (1.1) leads to a discrete version of it, which has the form where u and f (u, t) are finite-dimensional vectors of size n u , the size of the spatial discretization. For finite-element discretization, M is a mass matrix. As above, we wish to know if a steady solution u (s) to (1.2) is stable.
Linear stability analysis addresses this question by examining the eigenvalues of the algebraic system Jv = λM v, (1.3) where J = ∂f ∂u (u (s) ) is the Jacobian matrix of f with respect to u, evaluated at u (s) ; see, for example, [13,Ch. 1]. A necessary condition for stability of u (s) is that all eigenvalues λ of (1.3) have negative real part. If any eigenvalue has positive real part, then there exists an arbitrary small perturbation γ such that if u (s) + γ is used as an initial condition for (1.2), the integrated solution will not revert to u (s) .
A problematic aspect of linear stability analysis is that it fails to account for transient effects that may take a long time to resolve. In particular, it may happen that the solution of the system (1.1) with initial condition u(x, 0) = u (s) (x) + γ(x), consisting of a small perturbation of a steady solution, exhibits large growth over a significant period of time even if u (s) is linearly stable. This is discussed for models of flow in [21,Sections 2.3,4.1], [23,Sections 20,22]. It can be explained using pseudospectra: the -pseudospectrum of the Jacobian matrix, defined for M = I in (1.3), is the set of eigenvalues of J + E for E ≤ . (A generalization to forms of M considered in the present study is discussed in [11].) Transient growth is exhibited when some elements of this set protrude into the right-half of the complex plane [23].
Our aim in this study is to develop and explore a simple procedure to study the sensitivity of the eigenvalues of (1.3) when the dynamical system comes from models of incompressible flow. As observed in [11], it is not practical to compute pseudospectra for the large-scale systems that arise in this setting. This difficulty is addressed in [11] by projecting such systems into invariant subspaces of (shifted versions of) J −1 M , which have smaller dimension and for which computation of pseudospectra is feasible. It is shown in [11] that these pseudospectra estimates provide interior bounds on pseudospectra of (1.3) as well as insight into transient growth of solutions.
In this work, we develop a complementary approach to study the sensitivity of the eigenvalues of (1.3) for models of incompressible flow. The methodology derives from a two-fold procedure: 1. Introduce a simple way to construct perturbed versions of the eigenvalue problem (1.3) using spatial perturbations that depend on a finite number of parameters. 2. Approximate the critical eigenvalues of the perturbed problem using a surrogate function defined by interpolation. This requires the solution of a relatively small number of perturbed eigenvalue problems determined from a special set of parameter values, using sparse-grid methods [1,22]. The surrogate function interpolates the critical eigenvalues obtained from these eigenvalue problems and provides a means of approximating the critical eigenvalues for an additional set of perturbed problems. The surrogate function is very inexpensive to evaluate. As a result, it is possible to generate many samples of (approximate) eigenvalues in order to gain an understanding of the effects of perturbation. We apply this technique to the eigenvalue problems arising from stability analysis of the incompressible Navier-Stokes equations.
An outline of the remainder of the paper is as follows. In Section 2, we describe the collocation strategy and show in detail how it is developed for the Navier-Stokes equations. In Section 3, we describe two benchmark problems we use to test the methodology and show how the perturbed eigenvalues behave with respect to Reynolds numbers and sizes of perturbation, and in Section 4, we demonstrate that the behavior of perturbed eigenvalues predicts the behavior of transient solutions obtained from perturbed flow conditions. Finally, in Section 5, we make some concluding remarks.
2. Approach. In this section, we describe the methodology we will use to explore the sensitivity to perturbation of the eigenvalue problem (1.3), which is based on sampling. We first outline the approach in general terms in Section 2.1, and then we continue in Section 2.2 with a more detailed statement of how the ideas are applied to a specific benchmark problem, the incompressible Navier-Stokes equations.
2.1. General approach. Let u (s) be a steady solution to (1.2) and let δ be a small perturbation of u (s) . We will specify δ = δ(ξ) to depend on a vector of parameters ξ := (ξ 1 , ξ 2 , . . . , ξ m ) T with δ(0) = 0, and we will explore a perturbed eigenvalue problemĴ with the aim of understanding the impact of the perturbation δ on the eigenvalues {λ}. One way to defineĴ is to evaluate the Jacobian at the perturbed velocity, J(u (s) , δ) := J(u (s) + δ). In this study, which concerns the incompressible Navier-Stokes equations, we will insist that the perturbation is not dissipative. Details on the structure of the perturbation and its parameter dependence are given in Section 2.2.
Remark 2.1. We call attention here to an important aspect of the issue under study. Classic linear stability analysis concerns the sensitivity of the steady solution u (s) to perturbation. Our (different) concern here, like that of [23], is the sensitivity of the eigenvalues λ to perturbation, and in particular whether the conclusions reached from stability analysis predict behavior. To highlight this distinction, we use different symbols for perturbation depending on context: γ is used for perturbations arising in linear stability analysis, and δ for perturbations of eigenvalue problems as in (2.1).
Given the eigenvalue problem (2.1), let g(ξ) := rightmost eigenvalue of (2.1), (2.2) where, if there is a complex conjugate pair of rightmost eigenvalues, g(ξ) can be taken to be the eigenvalue with positive imaginary part. One way to explore the sensitivity of (1.3) is by sampling ξ, that is, to evaluate g(ξ) for a large set of sample values of ξ. If this function is very sensitive, that is, if small changes in δ(ξ) lead to large changes in g(ξ), then linear stability analysis may not provide an accurate assessment of stability; conversely, if g is not sensitive to perturbation, then linear stability analysis is likely to yield insight. The point of view here is that the study of perturbation is done by sampling a large number of nearby problems. A potential downside is that this approach requires the solution of many eigenvalue problems (2.1), one for each choice of ξ and resulting δ(ξ), which tends to incur a high computational cost. To reduce this expense, instead of evaluating the function of (2.2) (by solving an eigenvalue problem), we will replace g(ξ) with an approximation, a surrogate function g (I) (ξ), which is inexpensive to compute and therefore can be evaluated cheaply for many samples of ξ. For this, we will use the method of collocation designed to construct approximations to functions on high-dimensional spaces [1,22]. This entails evaluation of g(ξ) at a relatively small number of special points, {ξ (1) , ξ (2) , . . . , ξ (n ξ ) }. The surrogate function is then taken to be the polynomial interpolant of g, where { k (ξ)} are multidimensional Lagrange interpolation polynomials, For the interpolation points, we use sparse grids derived from the extrema of onedimensional Chebyshev polynomials [1].
Remark 2.2. It might happen that there are multiple eigenvalues of (1.3) with the same rightmost real part and different imaginary parts. In this case, g(ξ) of (2.2) would also be multi-valued or nearly so, and the ideas presented here would need to be applied to each of the rightmost eigenvalues. As long as there are not too many such values, this would have minimal impact on costs.

2.2.
Application to the Navier-Stokes equations. We will explore these ideas when the dynamical system (1.1) comes from the incompressible Navier-Stokes equations, and we now describe a way to specify a perturbation δ(ξ) for this benchmark problem for use in (2.1). To this end, consider the Navier-Stokes equations For fixed time t ∈ (0, ∞), the weak formulation of (2.4) is to find u(·, t) ∈ H 1 E , p(·, t) ∈ L 2 (D) such that (2.5) Linear stability analysis uses a linearized form of the first (momentum) equation of (2.4)-(2.5). Given a steady velocity field u (i.e., u t = 0), consider a perturbation u + γ. Substitution of this perturbed velocity into the quadratic term from (2.4) gives where the approximation on the right is made under the assumption that γ is small. Addition of the diffusion operator and specification of a perturbed weak formulation then leads to a trilinear form associated with the linearized convection-diffusion operator, (2.6) Mixed finite-element discretization of (2.5) uses finite-dimensional subspaces X h 0 ⊂ H 1 E0 and Y h ⊂ L 2 (D) together with X h E ⊂ H 1 E containing functions that interpolate the Dirichlet boundary data at element nodes lying in ∂D D . We will assume that this discretization is div-stable [15,Sect. 2.2]. The discrete weak formulation is to find be a discrete steady solution to (2.7), i.e., [ u where the aim is to find eigenvalues λ h and associated eigenfunctions satisfying (2.8) Here, we have linearized around a steady flow velocity field u (s) h satisfying (2.7).

Remark 2.3.
A complete discussion of the development of the trilinear form a(·, ·; ·) of (2.6) and the derivation of (2.8) is given in [10,]. This form also arises from use of Newton's method for solving the nonlinear system of equations arising from implicit time discretization of (2.5).
Let the dimensions of X h 0 and Y h be n u and n p , respectively. Let u (s) be the vector of coefficients of the steady finite-element solution u Here, F = F (u (s) ) is the matrix of order n u derived from the bilinear form a(·, ·; u (s) h ), B and B T are matrix representations of negative-divergence and gradient operators, respectively (B is of size n p × n u ), and Q is a velocity mass matrix, also of order n u .
Remark 2.4. The matrix on the right side of (2.9) is singular, and the resulting infinite eigenvalue can lead to instability in eigenvalue computations [19]. This can be avoided by replacing the matrix by −Q αB T αB 0 , which leaves the finite eigenvalues intact and maps the infinite eigenvalue to 1/α, see [4]. 1 We want to explore the sensitivity of our modified eigenvalue problem to perturbation. For this, we add a small perturbation The perturbation is defined in two steps. First, we specify a discretely divergence-free vector field δ h , that is, one satisfying This ensures that the perturbed velocity u h + δ h would be appropriate as an initial condition for testing the stability of u h . Second, the field δ h is used to generate a nondissipative perturbation . This means that the perturbation does not introduce any damping effects associated with numerical diffusion. To illustrate the construction, we will suppose that T h denotes a subdivision of D ⊂ R 2 into triangular or rectangular elements. The extension to three-dimensional problems is perfectly straightforward.
Note that the local construction (2.11) generates a discontinuous velocity field so that in general δ h ∈ H 1 (D) d .
Claim 2.2. Let the perturbation operator on element k ∈ T h be given by h is skew-adjoint on T h , and we need to apply Green's theorem on each element: where the last equality follows from the fact that δ is skew-adjoint. Summation over all the elements establishes the same property for c h .
The perturbed variant of (2.8) is This leads to the perturbed matrix eigenvalue problem (2.1) so that in particular N is a skew-symmetric matrix, N T = −N , for all parameter values ξ independent of the boundary conditions of the flow problem. It remains to specify the finite element function φ h used in Claim 2.1 to define the vector field δ h . Following [20], we take φ h (x, ξ) ∈ H 1 (D) to be a parameter-dependent scalar potential specified using a covariance function C(x (1) , x (2) ) for x (i) ∈ D. In particular, given C, let C := C(x, x) be the covariance matrix of order n consisting of the vertices in the subdivision associated with X h 0 , so that C ij = C (x i , x j ). Now let φ be an n-dimensional zero-mean stationary random process with covariance matrix C, i.e., C = E(φφ T ), where "E" refers to expected value. If C = σ 2 V ΘV T is an eigenvalue-eigenvector decomposition (scaled by the variance), then φ can be defined using a discrete Karhunen-Loève (KL) expansion where the eigenvector v j is the jth column of V and {ξ j } n j=1 are uncorrelated random variables with zero mean and unit variance [18,Section 5.4]. It is often the case that many of the eigenvalues are small and some of the terms in (2.15) can be removed without significant loss of accuracy. We will choose m < n such that n j=m+1 θ j n j=1 θ j ≤ 5/100, and, in the sequel, ξ := (ξ 1 , . . . , ξ m ) T will represent an m-dimensional vector of parameters and φ(ξ) is defined using the truncated KL-expansion This ξ-dependent coefficient vector (of length n) then characterizes a piecewise-defined linear or bilinear function φ h (ξ). For the computational results described in Section 3, we take the smooth covariance function where c 1 and c 2 are correlation lengths. We will also assume that {ξ j } in (2.16) are mutually independent, with each satisfying a truncated Gaussian distribution with range [−3, 3], so that ξ j has the density function A Matlab implementation of this distribution is given in [3] and described in [2]. We note two differences between this formulation of perturbations and traditional approaches based on pseudospectra. First, the perturbed eigenvalue problem (2.1), as specified in (2.13), is restricted to have a structure determined from that of the original problem, so that the resulting perturbed eigenvalues are closer in form to structured pseudoeigenvalues [23,Ch. 50]. Moreover, the perturbation itself derives explicitly from the nonlinear term u · ∇ u in (2.4). Indeed, if the original problem (1.1) were linear so that the Jacobian J did not depend on u, then the eigenvalues of J might still be sensitive to perturbation but the ideas discussed here would not give insight into this. Despite these limitations, as will be shown in Sections 3-4, the behavior of the perturbed eigenvalues gives insight into transient growth of solutions and other features of transient solvers.
We conclude this section with an analytic result bounding the size of the eigenvalue perturbation in proportion to the perturbation size.
Theorem 2.1. If δ h is the perturbation defined using (2.11) and λ is the rightmost eigenvalue of (1.3), then the eigenvalueλ of the perturbed problem (2.1) closest to λ Proof. Let M denote the matrix on the right side of (2.13), and let F and F = F + N denote the unperturbed and perturbed matrices on the left sides of (2.9) and (2.13), respectively; here N = N (ξ) = N (ξ) 0 0 0 . We are interested in the eigenvalue problems Fv = λMv and Fv =λMv. M can be factored as The velocity mass matrix Q and Schur complement Thus, we can consider the unperturbed and perturbed standard eigenvalue problems then we can use the Bauer-Fike Theorem [12,Ch. 7] to explore the eigenvalue perturbation: given an eigenvalue λ of A, We seek a bound on E 2 = Ẽ 2 , whereẼ = LN L −T ; the equality here follows from the fact that D is unitary. The structure of L leads tõ This holds for all α = 0, and so it also holds in the limit as α → 0, giving Ẽ 2 ≤ L −1 N L −T 2 . Since both N and L −1 N L −T are skew-symmetric, it follows that E 2 ≤ ρ(L −1 N L −T ), the spectral radius.
Thus, we require a bound on the Rayleigh quotient |(y,N y)| (y,Qy) . For this, we use (2.14) together with a standard bound on |c h | (see [10, p. 243 The assertion then follows from the inverse estimate 3. Benchmark problems and structure of eigenvalues. We will illustrate these ideas for two benchmark problems. In this section, we specify the problems and their features of interest, the eigenvalues associated with linear stability analysis and the effect of perturbation of these eigenvalues. Each of these is a model of flow through a channel for which there are inflow and outflow boundaries. The position of the outflow boundary is far enough downstream that the flow is fully developed. The spatial approximation is done using Q 2 -P −1 (biquadratic velocity; discontinuous linear pressure) mixed approximation [10, Section 3.3.1], implemented in the ifiss software package [8,9]. Unless otherwise specified, the discretization is done on a uniform grid with element width h = 1/32, which gives 64 elements along the vertical interval [−1, 1]. This corresponds to "grid level" = 6 in ifiss, with element width h = 2/2 . For both benchmark problems, we will explore the stability of solutions obtained for choices of the viscosity parameter near the critical values for linear stability. We determined these critical values experimentally, as described for example in [7]. This flow problem exhibits a pitchfork bifurcation [6]: as the viscosity decreases through a critical value (approximately ν = 1/220.5), the rightmost eigenvalue of (2.9), which is real, changes from negative (indicating linear stability) to positive (instability). Figure 3.2 shows the ten rightmost eigenvalues for three values of ν in this range and the rightmost eigenvalues (detail in the inset) for each choice, whose values are also identified on the right.
We explored the sensitivity of the rightmost negative eigenvalues using the perturbed eigenvalue problem (2.13). This was derived using correlation lengths c 1 = c 2 = 2 in (2.17), which resulted in a finite expansion (2.16) with m = 19 terms. The surrogate function g (I) of (2.3) used to estimate eigenvalues was constructed from a two-level sparse grid on the m-dimensional parameter space, which in turn resulted in n ξ = 761 sparse grid nodes. Thus, it is necessary to solve 760 eigenvalue problems, that is, find the rightmost eigenvalues of 760 perturbed systems (2.13), one for each  sparse-grid node. (One of the sparse-grid nodes is ξ = 0, which corresponds to an unperturbed system.) Once these are available, the estimates of eigenvalues for other choices of ξ can be obtained by evaluating g (I) . We implemented the sparse-grid interpolation using the matlab toolbox spinterp [16,17].
The dependence of the estimated perturbed eigenvalues on the perturbation size is illustrated in Figure 3.3. We found that insight can be provided using small values of the standard deviation in (2.16), and in these tests we used σ = βh 2 for β = . . Also note that the critical eigenvalues move to the left with mesh refinement, indicating that stability of the discrete systems is enhanced with mesh refinement, although it is also clear that a limiting value is approached with refinement.

Flow around a square obstacle.
In this case, the domain is a rectangular duct containing a square obstacle, with boundary conditions For this example, we used a level-6 stretched grid with local refinement near the obstacle. A representative steady solution that retains the reflectional symmetry is shown in Figure 3.5. In this case, there is a symmetry-breaking Hopf bifurcation for ν ≈ 1/186; that is, for ν in this range there is a complex conjugate pair of rightmost eigenvalues whose real parts change from negative to positive as ν is reduced. values of ν, the rightmost unperturbed eigenvalue (center of the set of perturbations) has negative real part, showing that the associated steady solution is stable, but for the close-to-critical value ν = 1/185.6, some of the perturbed eigenvalues have a positive real part. As for the step problem, it is readily seen that the magnitude of the perturbations does not depend on ν and varies linearly with σ.

Unsteady flow simulations.
In this section, we explore the connection between time integration of the Navier-Stokes equations and the eigenvalue perturbation results in the previous section. We will do this by computing time-accurate solutions of the Navier-Stokes equations using the adaptive (stabilized) Trapezoidal Rule (sTR) time stepping methodology built into ifiss. The suitability of sTR for long-time integration is discussed in [14]. Full details of the ifiss implementation of sTR can be found in section 10.2.3 of [10]. (Stabilization is based on time step averaging, which prevents the "ringing" to which TR is susceptible for stiff systems.) In what follows, we present results obtained from a nonlinear version of the integrator, denoted We model the laboratory scenario computationally via a two-stage process. 1. We start from a quiescent state and a tiny time step (1e-9). The inflow profile is smoothly ramped up to a fully developed flow using an exponential startup. The sTR2 integration is then carried out for 330 time steps with a relatively tight accuracy tolerance (i.e., a bound on an estimate of local truncation error), 3e-5. The number of steps taken is arbitrary but needs to be chosen large enough so that the reference flow is visually steady. More precisely, when this phase is complete, the instantaneous acceleration a(t), time steps, stopping prematurely only if the time reaches T * =1e14 -which we interpret as reaching a "computational steady-state" -at which point the adaptive time-stepping routine is taking very large time steps, see Figure  4.5 below, and the acceleration a(T * ) will almost certainly be smaller than unit roundoff. The mean vorticity ω(t), or the average vertical velocity at the outflow, provides a convenient way of assessing the degree of departure from the reflectionally symmetric flow (for which ω=0). At the conclusion of the first phase of the time integration, a pseudo-steady flow is computed for each of the three values of ν. The evolution of the mean vorticity and the acceleration visualized in Figure 4.1 shows that a symmetric flow is established for ν = 1/250 before the interruption is made after 330 time steps; this corresponds to time T ≈ 32.
Moving on to the second stage, we show results for the subcritical cases ν = 1/210 and 1/220, for three representative flow perturbations, each of which derives from a particular collocation point ξ (k) used in (2.3). For the first of these, no perturbation is made (this corresponds to the point ξ ≡ 0) and the integration simply continues from the first-stage stopping point T . The other two are representatives of a "benign" perturbation and a "lively" perturbation and have the spatial structure shown in The evolution of the flow after the restart for ν = 1/210 is depicted in Figure 4.3, where mean vorticity is shown in both actual and logarithmic scales. The unperturbed flow is perfectly stable; the sTR1 integrator reaches the end time (T * =1e14) at time step 396, 66 steps after the restart. The distinctive jumps in the acceleration are associated with the stabilization of the integrator, which has the effect of periodically injecting a small amount of dissipation into the flow. The benign perturbation, which respects the reflectional symmetry, has no effect on the long-term flow evolution. In contrast, the lively perturbation excites visible instability at about time step 390, 60 steps after the restart. But (as seen in particular from the acceleration), the size of the perturbation is not big enough to stop the long-term evolution to the symmetric flow at the designated end time T * . The growth in vorticity for the stable examples toward  Figure 4.4. The unperturbed case is just about stable: the sTR1 integrator reaches the end time 75 steps after the restart. A virtually identical evolution is evident when the perturbation to the flow is benign. (As in the previous example, roundoff effects lead to some growth in the vorticity after a steady solution is obtained.) The time evolution for the lively perturbation is noticeably different, however. In this case, the sTR1 integrator rejects time step 415 (85 steps after the restart) and the computational flow evolves to a numerically noisy solution where the magnitude of the oscillation is of the order of the time-stepping accuracy.
These observations are substantiated in Figure 4.5, which shows the history of the time step sizes chosen by the adaptive integrator. For each of the plots in this figure, the switch from the first to the second stage is identified by a vertical dotted line. When either no perturbation or a benign perturbation is made, the time step sizes rapidly increase because the integration goes to a steady state for the subcritical  bation is added at the interrupt point. The different outcomes for the two benchmark problems are representative of the difference between a pitchfork bifurcation (for the step) and a Hopf bifurcation (for the obstacle) [5], [10, p. 343].
To study the flow breakdown mechanism in detail, the second phase of the time integration is computed with a very small accuracy tolerance (1e-9) using the unstabilised TR1 integrator. 4 In all cases discussed below the time integrator is run for 2500 steps after the interrupt. Figure 4.9 shows the evolution of the mean vorticity and the acceleration using this refined strategy for each value of the viscosity parameter, when no perturbation is done. In the super-critical case of ν = 1/200 (bottom), there is a fast breakdown to the vortex-shedding solution. (Note that the evolution is plotted against physical time in this figure.) For both the subcritical (ν = 1/175) and near-critical (ν = 1/185.6) cases, there are long delays (until t ≈ 6e7 and t ≈ 1.1e4, respectively) after the interrupt, after which numerical instability kicks in and (as in the preceding section) generates a numerically noisy solution. The onset of instability is dramatically later for the subcritical case.
We explore the breakdowns in more depth in Figure 4.10, which shows magnified images of the noisy solution measures at the time they become unsteady. These images show that the magnitudes of the numerical oscillations (of order 1e-8 in the subcritical case and 1e-7 in the near-critical case) are comparable to the time-stepping accuracy. Even when no explicit perturbation is done, time accuracy plays a role in long-term simulation to compute steady solutions in near-critical regimes. Continuing this exploration of subcritical cases, we now consider the effects of perturbations at the interrupt. As in the previous section, we look at one perturbation that respects the reflectional symmetry of the flow solution in Figure 3.5 and is expected to be "benign" and one that breaks the reflectional symmetry and so is expected to be "lively". The results for ν = 1/175 are shown in Figure 4.11 and those for ν = 1/185.6 are in Figure 4.12. In these figures, the vertical scaling for the mean vorticities are now set to be equal in order to discern differences for the two viscosity values. These images should be compared with those corresponding to analogous experiments with no perturbation in Figures 4.9-4.10. In particular, for the sub-critical viscosity ν = 1/175 with either type of perturbation, after a long delay, the solution moves away from a steady state. This is not surprising, since the same phenomenon occurs when no perturbation is done. The onset of periodic behavior for the perturbed data is slightly earlier than for no perturbation (and earlier still for  the lively perturbation), but the magnitude of the oscillations is small. The results for ν = 1/185.6 ( Figure 4.12) bear some similarity to these -most notably, the behavior for the benign perturbation is virtually identical to that for no perturbation (middle of Figure 4.9). The structure of the oscillations for the near-critical viscosity is more like that for the super-critical viscosity (for both perturbations as well as without perturbation, compare the images for lively perturbation in Figure 4.12 with the images in Figures 4.9-4.10). In contrast, for the sub-critical viscosity ν = 1/175, the structure of the oscillations is is more like that arising when no perturbation is done. But for the lively perturbation and ν = 185.6, the onset of unstable behavior is significantly earlier (bottom of Figure 4.12). Finally, when we check to see what happens when the perturbation is significantly larger (of the order of the perturbation made in computing the pseudo-eigenvalues in Figure 3.6) we observe that there is a big difference in the time-stepping behavior in any case where the perturbation is not benign. This is illustrated by the results shown in the bottom plot in Figure 4.13. In this case the size of the lively perturbation is

Flow acceleration
Lively perturbation big enough to destabilize the integrator and a noisy periodic solution is computed. This mirrors the vortex-shedding solution that is computed in the unstable case but has an amplitude that is too small to be seen when plotted.
5. Concluding remarks. Our aims in this study were twofold. First, we have developed a new approach to assess the stability of dynamical systems by constructing perturbed systems based on collocation methods. This is reminiscent of methods for computing pseudospectra, but it has the advantage that the process of sampling (approximate) spectra is significantly less costly. Second, we compared the results of such assessments with the performance of time-stepping computations for a nontrivial application, the incompressible Navier-Stokes equations. In particular, for two benchmark problems, we examined the behavior of a stable integration scheme for simulating transient behavior for values of the viscosity in the system that are "sub-critical", nearly critical (very slightly smaller than the critical value), and super-critical.
In general, we found that the predictions of instability made by the collocation

Flow acceleration
Lively perturbation method were consistent with the behavior of integrators: in the nearly critical regime (of parameter values, viscosity in this case), there is more sensitivity to perturbation than in the sub-critical regime, and outcomes are qualitatively like those for super-critical parameters. We also note that making such assessments is complicated somewhat by the delicate nature of computations in regimes at or near stability limits. Eigenvalues and pseudoeigenvalues are not the sole determining factor affecting stability; the form of the perturbation also plays a significant role.