eigentools: A Python package for studying differential eigenvalue problems with an emphasis on robustness

1 Department of Physics and Astronomy, Bates College 2 Department of Mathematics, MIT 3 School of Natural Sciences, Institute for Advanced Study 4 CIERA, Northwestern University 5 Department of Astrophysical and Planetary Sciences, University of Colorado, Boulder 6 School of Mathematics and Statistics, University of Sydney 7 Department of Engineering Sciences and Applied Mathematics, Northwestern University DOI: 10.21105/joss.03079


Summary and Statement of need
Linear stability analysis of partial differential equations (PDEs) is a fundamental tool in chaotic dynamics, fluid dynamics, biophysics, and many other scientific disciplines. These analyses require solving eigenvalue problems. However, studying these eigenvalues is not without significant peril: discretized PDEs are poorly conditioned and many of the numerical eigenvalues reported by standard solvers are unreliable. Additionally, it is far from trivial to start with a set of PDEs and construct a matrix discretization in the first place. In order to solve these problems, we present eigentools, a Python package that extends the eigenvalue problem (EVP) capabilities of the Dedalus Project  to provide a complete analysis toolkit for EVPs. eigentools provides a convenient, parallelized interface for both modal and non-modal stability analyses for nearly arbitrary systems of PDEs. It provides a toolkit that makes it very easy to take a model, find robust eigenvalues and eigenmodes, and find critical parameter values for stability. The only thing a user needs to do is find a state to linearize about. One constructs a Dedalus EigenvalueProblem object for the PDE linearized about the chosen background, and passes that to eigentools. eigentools provides robust spurious eigenvalue rejection (Boyd, 2001), spectrum and eigenmode visualization, ε−pseudospectra, and the ability to project a given eigenmode onto an 2-or 3-dimensional domain and save it as a Dedalus-formatted HDF5 file to use as an initial condition for a non-linear simulation of the same system.

Robust Eigenvalues and Finding Critical Parameters
The core eigentools class, Eigenproblem, provides a simple interface to accurately solve a Dedalus eigenvalue problem using sparse or dense methods. A common use of eigenvalue problems is stability analysis, in which one characterizes the exponential growth rate γ of a mode via (typically) the temporal part of the eigenvalue, e γt with γ ∈ R. When γ > 0, the mode is unstable. Eigenproblem allows the user to choose how this "growth" is defined via custom functions of the eigenvalues returned by the solver. This permits instability to correspond to positive real parts of the eigenvalue (e.g. e σt , σ = γ + iω and γ, ω ∈ R), negative imaginary parts of the eigenvalue (e.g. e i (kx−ωt) , ω = ω r + iω i , ω r , ω i ∈ R), or any other choice. Because every eigenvalue has a corresponding growth rate, Eigenproblem returns an overall growth rate defined as the robust eigenvalue with the highest growth rate. Eigenproblem also has simple tools to plot all components of eigenmodes corresponding to a selected eigenvalue. In order to find robust eigenvalues, eigentools performs mode rejection by solving the same problem twice, once at a user specifiable multiple (1.5 by default) of the original resolution. In order to ascertain which modes are good, we calculate a figure of merit called the inverse drift ratio one of two ways (for details, see Chapter 7 of Boyd, 2001). For simple problems with only one mode family, one can use the ordinal method in which the eigenvalues are compared in sorted order. However, many important problems have multiple wave families. By increasing the resolution, the number of resolved modes for each family increases; because of this, one must compare the drift ratios of the nearest eigenvalues between the two resolutions. The right panel of Figure 1 shows both ordinal and nearest drift ratios; by the ordinal criterion, all eigenvalues would be rejected, despite the fact that many eigenvalues are robust. This shows that for this problem, one must use the nearest method. The MRI spectrum at the critical parameters. Inverse drift ratios for modes shown in the spectrum. Those below 10 6 are rejected according to nearest (blue) and ordinal (orange) criteria. For this problem, the nearest criterion must be used.
One of the original motivations for eigentools was to quickly and easily find critical parameters for eigenvalue stability problems. Critical parameters are those at which the fastest growing mode has zero growth rate. In order to do so, eigentools provides an object CriticalFinder, which allows users to specify an Eigenproblem and a 2-dimensional tuple of parameters. The user then provides a grid of points for those two parameters, and CriticalFinder finds the maximum growth rate for the EVP at each point in the grid, exploiting MPI parallelism on multiprocessor systems. It then interpolates to find the zero crossings of one parameter, and finally minimizes over the remaining parameter to find approximate critical values. CriticalFinder also provides simple visualization tools, root polishing algorithms to further refine the critical values, and the ability to save and load the grid of eigenvalues. Figure 1 demonstrates three core features of eigentools: the ability to find critical parameters, the ability to use sparse and dense eigenvalue solvers, and the ability to reject spurious eigenvalues. In the left panel, the growth rate of the magnetorotational instability (defined as the positive real part of the eigenvalue σ) is plotted on a 20 × 20 grid of magnetic Reynolds number Rm and wavenumber Q, finding the critical values Rm c = 4.88, Q = 0.747 (Clark & Oishi, 2017b); in Figure 1, we used 4 cores each performing 100 sparse eigenvalue solves finding the 15 modes with σ closest to zero. The middle panel shows the spectrum at the critical parameters; this was solved using a dense eigenvalue solver to find all modes. The unstable mode is a rotationally modified Alfvén wave highlighted in red. Finally, the rightmost panel shows a plot of the inverse drift ratio 1/δ for both ordinal and nearest comparisons.

Output and creation of initial conditions
eigentools can output eigenmodes in the Dedalus HDF5 data format so they can easily be used as initial conditions for non-linear simulations. Figure 2 highlights this capability. We solve an EVP at Ra = 10 6 for Rayleigh-Benard convection between two no-slip plates using eigentools at a resolution of n z = 32, and select the most unstable mode. We then project that mode on a 2-D domain, write it to disk, and load the data into a Dedalus initial value problem (IVP) solver using the full, non-linear equations for Rayleigh-Benard convection. Using Dedalus's ability to change parameters and resolutions on the fly, we then run IVP with a resolution of (512, 64) until it reaches a non-linear steady state. The IVP in Figure 2 was run in parallel on 32 cores. In the right panel of Figure 2, we see excellent agreement between the growth rate from the non-linear IVP and the initial eigenvalue until non-linearity begins to become important around t ≈ 0.01.

Pseudospectra
The eigentools package can also solve for the ε-pseudospectra (Trefethen & Embree, 2005) of generalized eigenvalue problems of the form using the recent algorithm given by Embree & Keeler (2017). To our knowledge, this is the first open-source system for computing ε−pseudospectra for arbitrary generalized EVPs, including those with singular M matrices. Such M arise quite commonly in the solution of temporally differential-algebraic equations, which occur any time there are algebraic equations in the system such as linear constraints (e.g. ∇·u = 0) or boundary conditions. The pseudospectrum shows the sensitivity of the eigenvalue λ to a parameter ε for bounded perturbations in the underlying operators The notion of pseudospectra relies on a particular choice of operator norm; in fluid dynamics, this is often the energy inner product. eigentools allow the user to implement any norm useful for their particular problem, provided it is induced by an inner product on the underlying vector space V : Oishi et al., (2021). eigentools: A Python package for studying differential eigenvalue problems with an emphasis on robustness. Journal of Open Source Software, 6(62), 3079. https://doi.org/10.21105/joss.03079 The pseudospectrum identifies the robust parts of the spectrum. Pseudospectra have numerous applications, including non-modal stability analysis in fluid dynamics (Schmid & Henningson, 2012), locating topological states (Loring, 2015), and helping to understand the unusual properties of PT -symmetric quantum mechanics (Krejčiřík et al., 2015).  Figure 3 shows an example pseudospectrum, its corresponding spectrum, and four representative eigenvectors for the classic Orr-Sommerfeld problem in hydrodynamic stability theory Trefethen et al., 1993). As a twist on the standard problem, we demonstrate Dedalus and eigentools ability to solve the problem using the standard Navier-Stokes equations linearized about a background velocity, rather than in the traditional, single fourthorder equation for wall-normal velocity. This is not possible without using the generalized eigenvalue pseudospectra algorithm implemented above. Note that for the four eigenvectors, we plot u and w, the stream wise and wall-normal directions, respectively, rather than w and η, the vorticity as would be the case in the reduced Orr-Sommerfeld form. The solid and dashed lines represent the real and imaginary parts of the eigenvectors, respectively.

Example
Here we present a script that computes the spectra and pseudospectra for the Orr-Sommerfeld problem, producing a simplified version of the center plot in Figure 3. The first block of code sets up the Navier-Stokes equations in Dedalus, making use of its expressive substitution mechanism.

Related Work
There are a few other packages dedicated to the automatic construction of eigenvalue problems, including Chebfun, which can also produce pseudospectra. Chebfun, while itself released under the standard 3-Clause BSD license, is written in the proprietary MATLAB language. For computing spectra and pseudospectra for existing matrices, the venerable EigTool package is another open-source option, also written in MATLAB. It does not feature parallelism nor the ability to construct matrices from PDEs. EigTool has also been ported to the open-source Julia language in the Pseudospectra.jl package.