%
% Example 2 by 2 matrix polynomial pseudospectrum
% in the Lagrange basis, using Mandelbrot polynomials
%
% RMC 13.12.2010 after an idea by Piers Lawrence
% and after discussion with Nic Fillion
%
%
% This first section just sets up the matrix polynomial of degree
% 2^level - 1, formed by interpolating at n+2 = 2^(level-1)+1 points
% giving 2 by 2 matrices at each interpolation point.
%
level = 6;
p = 1;
rts = [-1];
for i=2:level,
    % Evaluate the matrix polynomials at x=0, x=-2, and x=all roots
    % previous.
    n = 2^(i-1) - 1; % n roots available
    tau = [0, rts', -2 ]; % n+2 evaluation points
    gamma = genbarywts( tau, 1 );
    gamma = gamma/norm(gamma,inf);
    A = zeros( 2*(n+2)+2, 2*(n+2)+2 ); % extra block at end
    B0 = eye( 2*(n+2)+2, 2*(n+2)+2 );
    B0(end-1:end, end-1:end ) = zeros(2,2);
    A(1:2, 1:2) = [ 0, 0; 0, 0 ];
    A(1:2, end-1:end) = [ 0, -1; 1, 1];
    A(end-1:end,1:2) = [gamma(1), 0; 0, gamma(1)];
    for j=1:n,
        A( 2 + 2*j-1:2+2*j, 2+2*j-1:2+2*j ) = [ tau(j+1), 0 ; 0, tau(j+1) ];
        A( 2+2*j-1:2+2*j, end-1:end ) = [ 0, -1; 1, 0 ];
        A( end-1:end, 2+2*j-1:2+2*j ) = [ gamma(j+1), 0; 0, gamma(j+1) ];
    end;
    A(end-3:end-2,end-3:end-2) = [-2, 0; 0, -2];
    A(end-3:end-2,end-1:end) = [2, -1; 1, -1];
    A(end-1:end,end-3:end-2) = [gamma(end),0; 0, gamma(end)];
    rts = eig( A, B0 );
    k = find( abs(rts) < 3 );
    rts = rts(k);
end;
%
% At this point, we have our matrix, interpolated at the nonlinear
% eigenvalues of the previous level matrix polynomial.

% Compute the weights for B as the 2-norms of the polynomial coefficients
alpha = zeros( 1, n+2 );
for i=1:n+2,
    alpha(i) = norm( A( 2*i-1:2*i, end-1:end ), 2 );
end;

% Now evaluate the matrix polynomial at a bunch of points and compute 
% the 2-norm of its inverse, by simple svd computation.
nsamp = 400;
x = linspace(-2.5,0.7,nsamp);
y = linspace(-2,2,nsamp);
% z(j,k) = x(j) + i*y(k)
V = zeros(nsamp,nsamp);
for j=1:nsamp,
    for k=1:nsamp,
        z = x(j)+sqrt(-1)*y(k);
        w = prod( z - tau );
        P = gamma(1)*[ 0, -1; 1, 1]/(z-0) + gamma(end)*[2, -1; 1, -1]/(z+2);
        B = abs(alpha(1)*gamma(1)/z) + abs(alpha(end)*gamma(end)/(z+2));
        for ell=1:n,
            P = P + gamma(ell+1)*[ 0, -1; 1, 0 ]/(z-tau(ell+1));
            B = B + abs( alpha(ell+1)*gamma(ell+1)/(z-tau(ell+1)) );
        end;
        S = svd( w*P );
        V(k,j) = S(2)/(abs(w)*B);
    end;
end;
% Choose among the following options by commenting out as appropriate.
% The first two lines, with level=2, were used to generate the figure
% in the text.
%[C,h] = contour( x, y, log(V)/log(10), [-2,-1.5,-1,-0.75, -0.5, -0.25], 'k-' );
[C,h] = contour( x, y, log(V)/log(10), [-22,-18,-14,-10,-6,-2], 'k-' );
clabel( C, h );
%contour( x, y, log(V)/log(10),[0,-2,-4], 'k-' );
%[C,h] = contour( x, y, log(V)/log(10), 'k-' );
%clabel( C, h );
axis('square')
