Public interface

Simple spectral portrait

Pseudospectra.spectralportraitFunction
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.

source

Importing a matrix or operator

Pseudospectra.new_matrixFunction
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?
source
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).

source

Setting up graphics subsystem

Pseudospectra.setpsplotterFunction
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.

source
Pseudospectra.setgsFunction
setgs(; headless=false, savefigs=true) => gs

Construct a GUIState for subsequent use by Pseudospectra functions.

Assumes plotting package has been chosen via setpsplotter().

source

Driver function

After a matrix has been ingested into a PSAStruct and the graphics subsystem has been established, the following function will compute pseudospectra and plot a spectral portrait:

Pseudospectra.driver!Function
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.
source

Pseudospectra computation

Pseudospectra.psa_computeFunction
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:
KeyTypeDefaultDescription
:levelsVector{Real}autolog10(ϵ) for the desired ϵ levels
:recompute_levelsBooltrueautomatically recompute ϵ levels?
:real_matrixBooleltype(A)<:Realis the original matrix real? (Portrait is symmetric if so.) This is needed because T could be complex even if A was real.
:proj_levRealThe proportion by which to extend the axes in all directions before projection. If negative, exclude subspace of eigenvalues smaller than inverse fraction. ∞ means no projection.
:scale_equalBoolfalseforce the grid to be isotropic?
:threadedBoolfalsedistribute computation over Julia threads?

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:
infoMeaning
0No error
-1No levels in range specified (either manually, or if matrix is too normal to show levels)
-2Matrix is so non-normal that only zero singular values were found
-3Computation cancelled
source

Eigen/Pseudo-mode computation and plotting

Pseudospectra.modeplotFunction
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.

source

Other computations

Pseudospectra.psa_radiusFunction
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
source
Pseudospectra.psa_abscissaFunction
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
source
Pseudospectra.numerical_rangeFunction
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.

source
Pseudospectra.numerical_abscissaFunction
numerical_abscissa(A)

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

Uses eigvals(). ω(A) provides bounds and limiting behavior for opnorm(exp(t*A)).

source

Other plots

Pseudospectra.mtxexpsplotFunction
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.

source
Pseudospectra.mtxpowersplotFunction
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.

source
  • Trefethen1999L.N.Trefethen, "Computation of pseudospectra," Acta Numerica 8, 247-295 (1999).