# Plotting resolvent of the matrix

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

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}]

• You can also define A like this: A=RotateLeft[IdentityMatrix[5]]. Commented Jul 8, 2022 at 4:57

Here is a solution using the default Norm:

resolvent=Inverse[z*ii-A];
aux[zz_?NumericQ]:=Block[{z=zz},Norm[resolvent]];
Plot3D[aux[x+I*y],{x,-2,2},{y,-2,2},PlotPoints->100]


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}] *)


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] &)]


(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.

• 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. Commented Jul 9, 2022 at 10:27
• 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. Commented Jul 9, 2022 at 14:09
• Thanks, this seems to be the paper you are referring to; will have a look. Commented Jul 9, 2022 at 15:03