Public interface

spectralportrait(A::AbstractMatrix; npts=100) => Plots object

compute pseudospectra of matrix A and display as a spectral portrait.

Pseudospectra are computed on a grid of npts by npts points in the complex plane, including a neighborhood of the spectrum. Contour levels are log10(ϵ) where ϵ is the inverse resolvent norm. This is a convenience wrapper for simple cases; see the Pseudospectra package documentation for more elaborate interfaces.

new_matrix(A::AbstractMatrix, opts::Dict{Symbol,Any}=()) -> ps_data

process a matrix into the auxiliary data structure used by Pseudospectra.

Options

• :direct::Bool: force use of a direct algorithm?
• :keep_sparse::Bool: use sparse matrix code even if A is not large?
• :real_matrix::Bool: treat A as unitarily equivalent to a real matrix?
• :verbosity::Int: obvious
• :eigA: eigenvalues of A, if already known
• :proj_lev: projection level (see psa_compute)
• :npts: edge length of grid for computing and plotting pseudospectra
• :arpack_opts::ArpackOptions: (see type description)
• :levels::Vector{Real}: contour levels (if automatic choice is not wanted)
• :ax::Vector{Real}(4): bounding box for computation [xmin,xmax,ymin,ymax]
• :scale_equal::Bool: force isotropic axes for spectral portraits?
• :threaded::Bool: distribute Z values over Julia threads?

new_matrix(A, opts::Dict{Symbol,Any}=()) -> ps_data

process a linear operator object into the auxiliary data structure used by Pseudospectra.

There must be methods with A for eltype, size, and mul!. It is up to the user to make sure that mul! is consistent with any options passed to the iterative solver (see documentation for xeigs).

setpsplotter(plotter::Symbol=:default)

Select a plotting package for use with Pseudospectra.

Currently :Plots and :PyPlot are implemented. Defaults to :Plots unless PyPlot is already imported without Plots.

setgs(; headless=false, savefigs=true) => gs

Construct a GUIState for subsequent use by Pseudospectra functions.

Assumes plotting package has been chosen via setpsplotter().

driver!(ps_data, opts, gs=defaultgs(); revise_method=false)

Compute pseudospectra and plot a spectral portrait.

If using an iterative method to get some eigenvalues and a projection, invokes that first.

Arguments

• ps_data::PSAStruct: ingested matrix, as processed by new_matrix
• gs::GUIState: object handling graphical output
• opts::Dict{Symbol,Any}:
  • :ax, axis limits (overrides value stored in ps_data).
  • other options passed to redrawcontour, arnoldiplotter!

Note that many options stored in ps_data by new_matrix() influence the processing.

When revising a spectral portrait (revise_method==true), the following entries in opts also apply:

• :arpack_opts::ArpackOptions,
• :direct::Bool.

psa_compute(T,npts,ax,eigA,opts,S=I) -> (Z,x,y,levels,info,Tproj,eigAproj,algo)

Compute pseudospectra of a (decomposed) matrix.

Uses a modified version of the code in ^{[Trefethen1999]}. If the matrix T is upper triangular (e.g. from a Schur decomposition) the solver is much more efficient than otherwise.

Arguments

• T: input matrix, usu. from schur()
• npts: grid will have npts × npts nodes
• ax: axis on which to plot [min_real, max_real, min_imag, max_imag]
• eigA: eigenvalues of the matrix, usu. also produced by schur(). Pass an empty vector if unknown.
• S: 2nd matrix, if this is a generalized problem arising from an original rectangular matrix.
• opts: a Dict{Symbol,Any} holding options. Keys used here are as follows:

Notes:

• Projection is only done for square, dense matrices. Projection for sparse matrices may be handled (outside this function) by a Krylov method which reduces the matrix to a projected Hessenberg form before invoking psa_compute.
• This function does not compute generalized pseudospectra per se. They may be handled by pre- and post-processing.

Outputs:

• Z: the singular values over the grid
• x: the x coordinates of the grid lines
• y: the y coordinates of the grid lines
• levels: the levels used for the contour plot (if automatically calculated)
• Tproj: the projected matrix (an alias to T if no projection was done)
• eigAproj: eigenvalues projected onto
• algo: a Symbol indicating which algorithm was used
• info: flag indicating where automatic level creation fails:

modeplot(ps_data, pkey [,z])

Extract and plot an eigenmode or pseudo-eigenmode for the matrix encapsulated in the Pseudospectra object ps_data. Use the value z if provided or prompt for one. If pkey is 1, find the pseudoeigenmode for z; otherwise find the eigenmode for the eigenvalue closest to z.

psa_radius(A,ϵ [,d]) -> r,z

Compute ϵ-pseudospectral radius for a dense matrix.

Quadratically convergent two-way method to compute the ϵ-pseudospectral radius r of a dense matrix A. Also returns a vector z of points where the pseudospectrum intersects the circle of radius r. Uses the "criss-cross" algorithm of Overton and Mengi.

The ϵ-pseudospectral radius is

   maximum(abs(z)) for z s.t. minimum(σ(A-zI)) == ϵ

Optional arg:

• d: eigenvalues of A, if known in advance

Keyword args:

• `verbosity: 0 for quiet, 1 for noise

psa_abscissa(A,ϵ [,d]) -> α,z

Compute ϵ-pseudospectral abscissa for a dense matrix.

Quadratically convergent two-way method to compute the ϵ-pseudospectral abscissa α of a dense matrix A. Also returns a vector z of points where the pseudospectrum reaches the abscissa. Uses the criss-cross algorithm of Burke et al.

The ϵ-pseudospectral abscissa is

   maximum(real(z)) for z s.t. minimum(σ(A-zI)) == ϵ

Optional arg:

• d: eigenvalues of A, if known in advance

Keyword args:

• `verbosity: 0 for quiet, 1 for noise

numerical_range(A, nstep=20) -> Vector{Complex}

Compute points along the numerical range of a matrix.

Note: this solves an eigensystem for each point, so may be expensive.

numerical_abscissa(A)

Compute the numerical abscissa of a matrix A, ω(A).

Uses eigvals(). ω(A) provides bounds and limiting behavior for norm(expm(t*A)).

make a surface plot of the current spectral portrait

mtxexpsplot(ps_data,dt=0.1,nmax=50; gs::GUIState = defaultgs(), gradual=false)

plot the evolution of ∥e^(tA)∥.

This is useful for analyzing linear initial value problems ∂x/∂t = Ax.

mtxpowersplot(ps_data, nmax=50; gs::GUIState = defaultgs(), gradual=false)

plot norms of powers of a matrix ∥A^k∥

This is useful for analyzing iterative linear algebra methods.