Let $A\in\mathbb{C}^{n\times n}$ and $\varepsilon>0$. The $\varepsilon$-pseudospectrum of $A$ is the set $$\sigma_{\varepsilon}(A)=\left\{z\in\mathbb{C} \middle| \|(A-zI)^{-1}\|\geq\frac{1}{\varepsilon}\right\}.$$ If, for some $z\in\mathbb{C}$, $A-zI$ is not invertible (that is $z$ is an eigenvalue of $A$), we put $\|(A-zI)^{-1}\|:=\infty$. So the spectrum is always contained in the pseudospectrum. The norm here is the operator norm or induced norm of a matrix (see wiki).

I would like to plot the level sets $$\left\{z\in\mathbb{C} \middle| \|(A-zI)^{-1}\|=\frac{1}{\varepsilon}\right\}$$ for various values of $\varepsilon>0$ and given $A\in\mathbb{C}^{n\times n}$ of large order $n$.

There is a package called EigTool for MATLAB which works nicely. However, I did not find anything similar for Mathematica. Is there any (hopefully straightforward) way to effectively plot the level sets in Mathematica?


I still haven't gotten around to implementing the Arnoldi-based algorithm of Toh and Trefethen, so what follows is straightforward, but not very efficient code for visualizing $\varepsilon$-pseudospectra:

mat = ExampleData[{"Matrix", "WEST0067"}];

With[{c = Norm[mat], n = Length[mat]}, 
     DensityPlot[-c Log10[First[
                 SingularValueList[mat - (x + I y) IdentityMatrix[n, SparseArray],
                                   -1, Method -> "Arnoldi", Tolerance -> 0]]],
                 {x, -2, 2}, {y, -2, 2}, ColorFunction -> "SolarColors"]]

pseudospectrum of  sparse matrix

Compare with the "spectral portrait" here.

It is quite easy to modify the code above to use ContourPlot[] instead of DensityPlot[], since all that matters is the plotting of the (logarithm of the) greatest singular value of the resolvent matrix $(\mathbf A-\lambda\mathbf I)^{-1}$ (equivalently, the least singular value of its inverse). Here's the example used in the EigTool download page:

(* Grcar matrix, https://arxiv.org/abs/1203.2390 *)
grcar[r : _Integer?Positive : 3, n_Integer?Positive] := 
     SparseArray[{{j_, k_} /; j == k + 1 :> -1, {j_, k_} /; 0 <= k - j <= r :> 1}, {n, n}]

mat = grcar[32]; eigs = Eigenvalues[N[mat]];

With[{c = Norm[mat], n = Length[mat]}, 
     ContourPlot[-c Log10[First[
                 SingularValueList[mat - SparseArray[Band[{1, 1}] -> x + I y, {n, n}],
                                   -1, Method -> "Arnoldi", Tolerance -> 0]]],
                 {x, -2, 4}, {y, -4, 4}, AspectRatio -> Automatic,
                 Contours -> Range[1, 6, 1/2], ContourShading -> None,
                 ContourStyle -> Black,
                 Epilog -> {Directive[Black, AbsolutePointSize[5]], Point[ReIm /@ eigs]}]]

pseudospectral contours for Grcar matrix

Coloring and other sundry styling is left to the interested reader.

  • $\begingroup$ Notes: 1. The matrix norm (c in the code above) is merely a scaling factor, and can be left out if desired. 2. Jeff Bryant did some previous work here. $\endgroup$ – J. M. will be back soon Jan 28 '17 at 18:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.