Stochastic Galerkin methods for linear stability analysis of systems with parametric uncertainty^††thanks: This work was supported by the U.  S.  National Science Foundation under grant DMS1913201.

The identification of instability in large-scale dynamical system is important in a number of applications such as fluid dynamics, epidemic models, pharmacokinetics, analysis of power systems and power grid, or quantum mechanics and plasma physics. A steady solution $u$ is stable, if when in a transient simulation it is introduced with a small perturbation as initial data and the simulation reverts to $u$, and it is unstable otherwise. This is of fundamental importance since unstable solutions may lead to inexplicable dynamic behavior. Linear stability analysis entails computing the rightmost eigenvalue of the Jacobian evaluated at $u$, and thus it leads to solution of, in general, large sparse generalized eigenvalues problems see, e.g. [4, 6, 8, 13, 16, 27] and the references therein. Typically, a complex pair of rightmost eigenvalues leads to a Hopf bifurcation, and a real rightmost eigenvalue may lead to a pitchfork bifurcation. The analysis is further complicated if the parameters in the systems are functions of one or more random variables. This is quite common in many real-world applications, since the precise values of model coefficients or boundary conditions are often not known. A popular method for this type of problems is Monte Carlo, which is known for its robustness but also slow convergence. In this study, we use spectral stochastic finite element methods [12, 17, 33, 34] with the main focus on the so-called stochastic Galerkin method, for the linear stability analysis of Navier–Stokes equation with stochastic viscosity. Specifically, we consider the parameterized viscosity given in the form of generalized polynomial chaos (gPC) expansion. In the first step, we apply the algorithms developed in [18, 30], see also [24], to find a gPC expansion of the solution of the Navier–Stokes equation. In the second step, we use the expansions of the solution and viscosity to set up a generalized eigenvalue problem with a nonsymmetric matrix operator, and in the assessment of linear stability of this problem we identify the gPC expansions of the rightmost eigenvalue. The main contribution in this study is development of stochastic Galerkin method for nonsymmetric eigenvalue problems. Our approach is based on inexact Newton iteration: the linear systems with Jacobian matrices are solved using GMRES, for which we also develop several preconditioners. The preconditioners are motivated by our prior work on (truncated) hierarchical preconditioning [32, 19], see also [2]. For an overview of literature on solving eigenvalue problems in the context of spectral stochastic finite element methods we refer to [1, 19, 29] and the references therein. Recently, Hakula and Laaksonen [15] studied crossing of eigenmodes in the stochastic parameter space, and Elman and Su [9] developed a low-rank inverse subspace iteration. However, to the best of our knowledge, there are only a few references addressing nonsymmetric stochastic eigenvalue problems: by Sarrouy et al. [25, 26], but there is no discussion of efficient solution strategies, and also by Sonday et al. [28], who studied distribution of eigenvalues for the Jacobian in the context of stochastic Galerkin method. Most recently, the authors with collaborators also compared surrogate learning strategies based on a sampling method, Gaussian process regression and a neural network in [31].

A study of linear stability of Navier–Stokes equation under uncertainty was conducted by Elman and Silvester [5]. The study was based on a judiciously chosen perturbation of the state variable and a stochastic collocation method was used to characterize the rightmost eigenvalue. Our approach here is different. We consider parametric uncertainty (of the viscosity), and the solution strategy is based on the stochastic Galerkin method. In fact, also the variant of the collocation method used here is based on the stochastic Galerkin projection (sometimes called a nonintrusive stochastic Galerkin method in the literature, see [33, Chapter 7] for a discussion). From this perspective, our study can be viewed as an extension of the setup from [6] to problems with viscosity given in the form of stochastic expansion and their efficient solution using stochastic Galerkin method. However, more importantly, we illustrate that the inexact methods for stochastic eigenvalue problems proposed recently by Lee and Sousedík [19] can be also applied to problems with nonsymmetric matrix operator¹¹1Specifically, the methods based on inexact Newton iteration, since in our experience the stochastic Rayleigh quotient and inverse iteration methods are less effective for nonsymmetric problems.. This in general allows to perform a linear stability analysis for other types of problems as well. We do not address eigenvalue crossing here, which is a somewhat delicate task for gPC-based techniques. We assume that the eigenvalue of interest is sufficiently separated from the rest of the spectrum and no crossing occurs. This is often the case for outliers and other eigenvalues that may be of interest in applications. A suitability of the algorithm we propose in this study can be assessed, e.g., by running first a low-fidelity (quasi-)Monte Carlo simulation. We also note that from our experience with the problem at hand indicates that the rightmost eigenvalue remains relatively well separated from the rest of the spectrum, and it is less prone to switching unlike the other eigenvalues, even with moderate values of coefficient of variation, which in turn allows its efficient characterization by a gPC expansion.

The rest of the paper is organized as follows. In Section 2 we recall the basic elements of linear stability analysis and link it to an eigenvalue problem for a specific model given by Navier–Stokes equation, in Section 3 we introduce the stochastic eigenvalue problem, in Section 4 we formulate the sampling methods and in Section 5 the stochastic Galerkin method and inexact Newton iteration for its solution, in Section 6 we apply the algorithms to linear stability analysis of Navier–Stokes equation with stochastic viscosity, in Section 7 we report the results of numerical experiments, and finally in Section 8 we summarize and conclude our work.

Following [6], let us consider the dynamical system

where $f:\mathbb{R}^{n}\times\mathbb{R}{\color[rgb]{0,0,0}\rightarrow}\mathbb{R}^{n}$ is a nonlinear mapping, $u\in\mathbb{R}^{n}$ is the state variable (velocity, pressure, temperature, deformation, etc.), in the finite element setting $M\in\mathbb{R}^{n\times n}$ is the mass matrix, and $\nu$ is a parameter. Let $u$ denote the steady-state solution to (1), i.e., $u_{t}=0$. We are interested in the stability of$~{}u$: if a small perturbation $\delta(0)$ is introduced to $u$ at time $t=0$, does $\delta(t)$ grow with time, or does it decay? For a fixed value of $\nu$, linear stability of the steady-state solution is determined by the spectrum of the eigenvalue problem

where $J=\frac{\partial f}{\partial u}(u(\nu),\nu)$ is the Jacobian matrix of $f$ evaluated at $\nu$. The eigenvalues have a general form $\lambda=\alpha+i\beta$, where $\alpha=\operatorname{Re}\lambda$ and $\beta=\operatorname{Im}\lambda$. Then $e^{\lambda t}=e^{\alpha t}e^{i\beta t}$, and since $\left|e^{\lambda t}\right|=\left|e^{\alpha t}\right|$, there are in general two cases: if $\alpha<0$ the perturbation decays with time, and if $\alpha>0$ the perturbation grows. We refer, e.g., to [4, 13] and the references therein for a detailed discussion. That is, if all the eigenvalues have strictly negative real part, then $u$ is a stable steady solution, and if some eigenvalues have nonnegative real part, then$~{}u$ is unstable. Therefore, a change of stability can be detected by monitoring the rightmost eigenvalues of (2). A steady-state solution may lose its stability in one of two ways: either the rightmost eigenvalue of (2) is real and passes through zero from negative to positive as $\nu$ varies, or (2) has a complex pair of rightmost eigenvalues and they cross the imaginary axis as $\nu$ varies, which leads to a Hopf bifurcation with the consequent birth of periodic solutions of (1).

Consider a special case of (1), the time-dependent Navier–Stokes equations governing viscous incompressible flow,

subject to appropriate boundary and initial conditions in a bounded physical domain$~{}D$, where $\nu$ is the kinematic viscosity, $\vec{u}$ is the velocity and $p$ is the pressure. Properties of the flow are usually characterized by the Reynolds number

where $U$ is the characteristic velocity and $L$ is the characteristic length. For convenience, we will also sometimes refer to the Reynolds number instead of the viscosity. Mixed finite element discretization of (3) gives the following Jacobian and the mass matrix, see [6] and [8, Chapter $8$] for more details,

where $n_{x}=n_{u}+n_{p}$ is the number of velocity and pressure degrees of freedom, respectively, $n_{u}>n_{p}$, $F\in\mathbb{R}^{n_{u}\times n_{u}}$ is nonsymmetric, $B\in\mathbb{R}^{n_{p}\times n_{u}}$ is the divergence matrix, and the velocity mass matrix $G^{n_{u}\times n_{u}}$ is symmetric positive definite. The matrices are sparse and $n_{x}$ is typically large. We note that the mass matrix$~{}M$ is singular, and (2) has an infinite eigenvalue. As suggested in [3], we replace the singular mass matrix$~{}M$ with the nonsingular, shifted mass matrix

which is symmetric but in general indefinite, and it maps the infinite eigenvalues of (2) to $\sigma^{-1}$ and leaves the finite ones unchanged. Then, the generalized eigenvalue problem (2) can be transformed into an eigenvalue problem

The flow is considered stable if $\operatorname{Re}\lambda<0$, and we wish to detect the onset of instability, that is to find when the rightmost eigenvalue crosses the imaginary axis. Efficient methods for finding the rightmost pair of complex eigenvalues of (2) (or (6)) were studied in [6]. Here, our goal is different. We consider parametric uncertainty in the sense that the parameter $\nu\equiv\nu(\xi)$, where $\xi$ is a set of random variables and which is given by the so-called generalized polynomial chaos expansion. To this end, we first recast the eigenvalue problem (6) in the spectral stochastic finite element framework, then we show how to efficiently solve it, and finally we apply the stability analysis to Navier–Stokes equation with stochastic viscosity.

Let $\left(\Omega,\mathcal{F},\mathcal{P}\right)$ be a complete probability space, that is $\Omega$ is the sample space with $\sigma$-algebra$~{}\mathcal{F}$ and probability measure$~{}\mathcal{P}$. We assume that the randomness in the mathematical model is induced by a vector $\xi:\Omega\mapsto\Gamma\subset\mathbb{R}^{m_{\xi}}$ of independent random variables $\xi_{1}(\omega),\dots,\xi_{m_{\xi}}(\omega)$, where $\omega\in\Omega$. Let $\mathcal{B}(\Gamma)$ denote the Borel $\sigma$-algebra on$~{}\Gamma$ induced by $\xi$ and $\mu$ the induced measure. The expected value of the product of measurable functions on$~{}\Gamma$ determines a Hilbert space $T_{\Gamma}\equiv L^{2}\left(\Gamma,\mathcal{B}(\Gamma),\mu\right)$ with inner product

where the symbol $\mathbb{E}$ denotes the mathematical expectation.

In computations we will use a finite-dimensional subspace $T_{p}\subset T_{\Gamma}$ spanned by a set of multivariate polynomials $\left\{\psi_{\ell}(\xi)\right\}$ that are orthonormal with respect to the density function $\mu$, that is $\mathbb{E}\left[\psi_{k}\psi_{\ell}\right]=\delta_{k\ell}$, and $\psi_{1}$ is constant. This will be referred to as the gPC basis [34]. The dimension of the space$~{}T_{p}$, depends on the polynomial degree. For polynomials of total degree$~{}p$, the dimension is $n_{\xi}=\binom{m_{\xi}+p}{p}$.

Suppose that we are given an affine expansion of a matrix operator, which may correspond to the Jacobian matrix in (2), as

where each$~{}K_{\ell}\in\mathbb{R}^{n_{x}\times n_{x}}$ is a deterministic matrix, and $K_{1}$ is the mean-value matrix $K_{1}=\mathbb{E}\left[K(\xi)\right]$. The representation (8) is obtained from either Karhunen-Loève or, more generally, a stochastic expansion of an underlying random process; a specific example is provided in Section 7.

We are interested in a solution of the following stochastic eigenvalue problem: find a specific eigenvalue $\lambda(\xi)$ and possibly also the corresponding eigenvector$~{}v(\xi)$, which satisfy in$~{}D$ almost surely (a.s.) the equation

where $K(\xi)\in\mathbb{R}^{n_{x}\times n_{x}}$ is a nonsymmetric matrix operator, $M_{\sigma}\color[rgb]{0,0,0}{\in\mathbb{R}^{n_{x}\times n_{x}}}$ is a deterministic mass matrix, $\lambda(\xi)\in\mathbb{\mathbb{C}}$ and $v(\xi)\in\mathbb{\mathbb{C}}^{n_{x}}$ along with a normalization condition

which is further specified in Section 5.

We will search for expansions of the eigenpair ${(\lambda(\xi),v(\xi))}$ in the form

where $\lambda_{k}\in\mathbb{\mathbb{C}}$ and $v_{k}\in\mathbb{\mathbb{C}}^{n_{x}}$ are coefficients corresponding to the basis$~{}\left\{\psi_{k}\right\}$. We note that the series for $\lambda(\xi)$ in (11) converges for $n_{\xi}\rightarrow\infty$ in the space $T_{\Gamma}$ under the assumption that the gPC polynomials provide its orthonormal basis and provided that $\lambda(\xi)$ has finite second moments see, e.g., [10] or [33, Chapter 5]. Convergence analysis of this approximation for self-adjoint problems can be found in [1].

Both Monte Carlo and stochastic collocation are based on sampling. The coefficients are defined by a discrete projection

The evaluations of coefficients in (12) entail solving a set of independent deterministic eigenvalue problems at a set of sample points$~{}\left\{\xi^{(q)}\right\}$, $q=1,\dots,n_{MC}$ or $n_{q}$,

In the Monte Carlo method, the sample points$~{}\xi^{(q)}$, $q=1,\dots,n_{MC},$ are generated randomly following the distribution of the random variables$~{}\xi_{i}$, $i=1,\dots,m_{\xi}$, and moments of solution are computed by ensemble averaging. In addition, the coefficients in (11) could be computed using Monte Carlo integration as²²2In numerical experiments, we avoid this approximation of the gPC coefficients and directly work with the sampled quantities.

For stochastic collocation, which is used here as the so-called nonintrusive stochastic Galerkin method, the sample points$~{}\xi^{(q)}$, $q=1,\dots,n_{q},$ consist of a predetermined set of collocation points, and the coefficients $\lambda_{k}$ and $v_{k}$ in the expansions (11) are determined by evaluating (12) in the sense of (7) using numerical quadrature as

where $\mathbb{\xi}^{(q)}$ are the quadrature (collocation) points and $w^{(q)}$ are quadrature weights. Details of the rule we use in our numerical experiments are discussed in Section 7, and we refer to monograph [17] for a detailed discussion of quadrature rules.

The stochastic Galerkin method is based on the projection

Let us introduce the notation

In order to formulate an efficient algorithm for eigenvalue problem (9) with nonsymmetric matrix operator using the stochastic Galerkin formulation, we introduce a separated representation of the eigenpair using real and imaginary parts,

where $v_{\operatorname{Re}}(\xi),v_{\operatorname{Im}}(\xi)\in\mathbb{R}^{n_{x}}$ and $\lambda_{\operatorname{Re}}(\xi),\lambda_{\operatorname{Im}}(\xi)\in\mathbb{R}$. Then, replacing $v(\xi)$ and $\lambda(\xi)$ in eigenvalue problem (9) results in

Expanding the terms in (17) and collecting the real and imaginary parts yields a system of equations that can be written in a separated representation as

and the normalization condition corresponding to the separated representation in (18) is taken as

Now, we seek expansions of type (11) for $v_{\operatorname{Re}}(\xi)$, $v_{\operatorname{Im}}(\xi)$, $\lambda_{\operatorname{Re}}(\xi)$, and $\lambda_{\operatorname{Im}}(\xi)$, that is

where the symbol $\otimes$ denotes the Kronecker product, $\Psi(\xi)=[\psi_{1}(\xi),\ldots,\psi_{n_{\xi}}(\xi)]^{T}$, $\bar{\lambda}_{s}=[\lambda_{s,1},\ldots,\lambda_{s,n_{\xi}}]^{T}\in\mathbb{R}^{n_{\xi}}$, and $\bar{v}_{s}=[v_{s,1}^{T},\ldots,v_{s,n_{\xi}}^{T}]^{T}\in\mathbb{R}^{n_{x}n_{\xi}}$ for $s=\operatorname{Re},\operatorname{Im}$.

Let us consider expansions (20) as approximations to the solution of (18)–(19). Then, we can write the residual of (18) as

and the residual of (19) as

Imposing the stochastic Galerkin orthogonality conditions (14) and (15) to$~{}\widetilde{F}$ and$~{}\widetilde{G}$, respectively, we get a system of nonlinear equations

where

and

We will use Newton iteration to solve system of nonlinear equations (21). The Jacobian matrix$~{}DJ(\bar{v}_{\operatorname{Re}},\bar{v}_{\operatorname{Im}},\bar{\lambda}_{\operatorname{Re}},\bar{\lambda}_{\operatorname{Im}})$ of (21) can be written as

where

and

However, for convenience in the formulation of the preconditioners presented later, we formulate Newton iteration with rescaled Jacobian matrix$~{}\widehat{DJ}(\bar{v}_{\operatorname{Re}}^{(n)},\bar{v}_{\operatorname{Im}}^{(n)},\bar{\lambda}_{\operatorname{Re}}^{(n)},\bar{\lambda}_{\operatorname{Im}}^{(n)})$, which at step$~{}n$ entails solving a linear system

followed by an update

Linear systems (26) are solved inexactly using a preconditioned GMRES method. We refer to Appendix A for the details of evaluation of the right-hand side and of the matrix-vector product, and to [18] for a discussion of GMRES in a related context.