Public interface
Simple spectral portrait
Pseudospectra.spectralportrait
— Functionspectralportrait(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.
Importing a matrix or operator
Pseudospectra.new_matrix
— Functionnew_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 ifA
is not large?:real_matrix::Bool
: treatA
as unitarily equivalent to a real matrix?:verbosity::Int
: obvious:eigA
: eigenvalues ofA
, if already known:proj_lev
: projection level (seepsa_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
: distributeZ
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
).
Setting up graphics subsystem
Pseudospectra.setpsplotter
— Functionsetpsplotter(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.
Pseudospectra.setgs
— Functionsetgs(; headless=false, savefigs=true) => gs
Construct a GUIState
for subsequent use by Pseudospectra functions.
Assumes plotting package has been chosen via setpsplotter()
.
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!
— Functiondriver!(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 bynew_matrix
gs::GUIState
: object handling graphical outputopts::Dict{Symbol,Any}
::ax
, axis limits (overrides value stored inps_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
.
Pseudospectra computation
Pseudospectra.psa_compute
— Functionpsa_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. fromschur()
npts
: grid will havenpts × npts
nodesax
: axis on which to plot[min_real, max_real, min_imag, max_imag]
eigA
: eigenvalues of the matrix, usu. also produced byschur()
. Pass an empty vector if unknown.S
: 2nd matrix, if this is a generalized problem arising from an original rectangular matrix.opts
: aDict{Symbol,Any}
holding options. Keys used here are as follows:
Key | Type | Default | Description |
---|---|---|---|
:levels | Vector{Real} | auto | log10(ϵ) for the desired ϵ levels |
:recompute_levels | Bool | true | automatically recompute ϵ levels? |
:real_matrix | Bool | eltype(A)<:Real | is the original matrix real? (Portrait is symmetric if so.) This is needed because T could be complex even if A was real. |
:proj_lev | Real | ∞ | The 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_equal | Bool | false | force the grid to be isotropic? |
:threaded | Bool | false | distribute 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 gridx
: the x coordinates of the grid linesy
: the y coordinates of the grid lineslevels
: the levels used for the contour plot (if automatically calculated)Tproj
: the projected matrix (an alias toT
if no projection was done)eigAproj
: eigenvalues projected ontoalgo
: a Symbol indicating which algorithm was usedinfo
: flag indicating where automatic level creation fails:
info | Meaning |
---|---|
0 | No error |
-1 | No levels in range specified (either manually, or if matrix is too normal to show levels) |
-2 | Matrix is so non-normal that only zero singular values were found |
-3 | Computation cancelled |
Eigen/Pseudo-mode computation and plotting
Pseudospectra.modeplot
— Functionmodeplot(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
.
Other computations
Pseudospectra.psa_radius
— Functionpsa_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
Pseudospectra.psa_abscissa
— Functionpsa_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
Pseudospectra.numerical_range
— Functionnumerical_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.
Pseudospectra.numerical_abscissa
— Functionnumerical_abscissa(A)
Compute the numerical abscissa of a matrix A
, ω(A)
.
Uses eigvals()
. ω(A)
provides bounds and limiting behavior for opnorm(exp(t*A))
.
Other plots
Pseudospectra.surfplot
— Functionmake a surface plot of the current spectral portrait
Pseudospectra.mtxexpsplot
— Functionmtxexpsplot(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
.
Pseudospectra.mtxpowersplot
— Functionmtxpowersplot(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.
- Trefethen1999L.N.Trefethen, "Computation of pseudospectra," Acta Numerica 8, 247-295 (1999).