I'm trying to reproduce resolvent plot of circulant matrix, page 59 of Mark Embree's slides

plot of matrix resolvent

Straightforward solution below, relying on ComplexPlot3D is below. It works for some matrices, but for this one it fails. Is there a more robust way to plot the resolvent?

n = 5;
A = IdentityMatrix[n];
A = A[[2 ;;]]~Join~{First[A]};
ii = IdentityMatrix[n];
A // MatrixForm
ComplexPlot3D[Norm@Inverse[z*ii - A], {z, -2 - 2 I, 2 + 2 I}]
  • 3
    $\begingroup$ You can also define A like this: A=RotateLeft[IdentityMatrix[5]]. $\endgroup$
    – user293787
    Commented Jul 8, 2022 at 4:57

3 Answers 3


Here is a solution using the default Norm:


For matrix,maybe use "Frobenius" norm.

norm = Norm[Inverse[z*ii - A], "Frobenius"];
Block[{z = x + I*y}, 
 Plot3D[norm, {x, -2, 2}, {y, -2, 2}, AxesLabel -> {"Re z", "Im z"}]]
(* ComplexPlot3D[norm, {z, -2 - 2 I, 2 + 2 I}] *)

enter image description here


It's a bit wasteful to invert just to take the 2-norm afterwards. Instead, recall that the 2-norm of $\mathbf A^{-1}$ is the reciprocal of the smallest singular value of $\mathbf A$, which is about the same amount of effort as computing the 2-norm (largest singular value) of $\mathbf A$. To accentuate the positions of the singularities (i.e. the eigenvalues), you can then take a logarithm:

mat = SparseArray[ToeplitzMatrix[UnitVector[5, 5], UnitVector[5, 2]]];
Plot3D[-Log[First[SingularValueList[SparseArray[Band[{1, 1}] -> x + I y, {5, 5}] - mat,
                                    -1, Tolerance -> 0]]], {x, -2, 2}, {y, -2, 2}, 
       BoxRatios -> Automatic, ClippingStyle -> None, 
       ColorFunction -> (ColorData["SolarColors", #3] &)]

resolvent of circulant matrix

(I also use this identity in this answer involving the $\varepsilon$-pseudospectrum.)

It should perhaps also be noted that Embree's example also happens to be the companion matrix of the polynomial $z^5-1$, which should explain why the eigenvalues are arranged that way.

  • $\begingroup$ Since you bring up efficiency: For every normal matrix mat (and OP's matrix is normal) there is an orthonormal basis relative to which it is diagonal and therefore evals=Eigenvalues[mat]; maxsval[z_]:=Max[Map[1/Abs[z-#]&,evals]] computes the max singular value of the resolvent, which is hard to beat. What is an efficient algorithm for a non-normal matrix? Perhaps a Schur decomposition (cf. this question of mine) would help but I have not thought it through. $\endgroup$
    – user293787
    Commented Jul 9, 2022 at 10:27
  • $\begingroup$ As I recall, Toh and Trefethen had an Arnoldi-based method for matrices that you could use after a preliminary Schur decomposition. I haven't experimented with that approach yet, tho. $\endgroup$ Commented Jul 9, 2022 at 14:09
  • $\begingroup$ Thanks, this seems to be the paper you are referring to; will have a look. $\endgroup$
    – user293787
    Commented Jul 9, 2022 at 15:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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