If $A$ is a square matrix, the resolvent of $A$ for a particular complex number $z$ is the matrix $(zI-A)^{-1}$. The $2$-norm of the resolvent is a quantity of interest in many applications. For example, if $|(zI-A)^{-1}| = 1/\varepsilon$ for some quantity $\varepsilon$, then there is a matrix $E$ with norm $|E| = \varepsilon$ such that $z$ is an eigenvalue of $A+E$. This is the starting point of the theory of pseudospectra [1].

In particular, suppose all the eigenvalues of $A$ are in the left half of the complex plane, so that $A$ is stable in the sense that all solutions of the differential equation $\frac{du}{dt} = Au$ decay to zero as $t \to \infty$. How small a perturbation matrix $E$ might make $A$ unstable? The answer is $|E| = \varepsilon$, where $1/\varepsilon$ is the maximum of $|(zI-A)^{-1}|$ as $z$ ranges over the imaginary axis. Therefore in a number of fields such as control theory, there is special interest in the values taken by the norm of the resolvent on the imaginary axis.

Let's compute this function with Chebfun. As an example we take the matrix

A = [-1 3 5 2; -3 -2 4 6; -5 -4 -2 1; -2 -6 -1 3]
A =
    -1     3     5     2
    -3    -2     4     6
    -5    -4    -2     1
    -2    -6    -1     3

A has two pairs of eigenvalues near the imaginary axis:

format short, format compact
eig(A)
ans =
  -0.7688 + 8.9660i
  -0.7688 - 8.9660i
  -0.2312 + 1.5019i
  -0.2312 - 1.5019i

Suppose $z=x+iy$. It takes Chebfun a fraction of a second to compute a chebfun representing $|(zI-A)^{-1}|$ as a function of $y$, with $x=0$. Here is that calculation and a plot of the result:

I = eye(size(A));
nr = @(y) 1/min(svd(1i*y*I-A));
f = chebfun(nr,[-25,25],'vectorize');
LW = 'linewidth';
plot(f,LW,1.6), grid on

The maximum of $f$ is this,

format long
maxf = max(f)
maxf =
   6.227545522966220

and the distance to instability is the reciprocal of this quantity,

dist_sing = 1/maxf
dist_sing =
   0.160576907276254

Let us consider another example matrix, and this time, let's make an anonymous function to construct the chebfun.

normfun = @(A) chebfun(@(y) 1/min(svd(1i*y*eye(size(A))-A)),...
   1.5*norm(A)*[-1,1],'vectorize');

Here is a $5\times5$ matrix which we take to be complex, to break the symmetry:

B =  [ -3-2i   1+1i    -1i      0   -1+1i
           0  -2-3i    -1i     1i   -2-1i
          1i      0  -2-4i  -2-1i    2-1i
           0      1     1i  -2-4i      1i
        1-2i      0      1      1   -2-3i ];
format short, eig(B)
ans =
  -5.3054 - 3.2003i
  -0.6662 - 0.8209i
  -0.3296 - 4.5158i
  -2.9797 - 3.2972i
  -1.7191 - 4.1659i

And here is its resolvent norm plot:

fB = normfun(B);
plot(fB,LW,1.6), grid on
title(['maximum = ' num2str(max(fB))]);

Here are 12 random $6\times6$ complex matrices, all with rightmost eigenvalue having real part $-0.25$:

rng(1)
for j = 1:12
    N = 6;
    A = randn(N) + 1i*randn(N) + 2i*diag(randn(N,1));
    abscissa = max(real(eig(A)));
    A = A - (abscissa+0.25)*eye(N);
    subplot(4,3,j)
    plot(normfun(A),LW,1)
    axis([-10 10 0 8]), drawnow
end

References

  1. L. N. Trefethen and M. Embree, Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators, Princeton U. Press, 2005.