## An Effective Algorithm for Computing the Numerical Range

### Carl C. Cowen and Elad Harel

In this paper, we collect a few fairly well known facts about the
numerical range and assemble them into an effective algorithm for
computing the numerical range of an n x n matrix. Particular
attention is paid to the case in which a line segment is embedded
in the boundary of the numerical range, a case in which multiplicity
is present. The result is a parametrization of the curve that forms
the boundary of the numerical range.

A Matlab implementation of the algorithm is appended below.

Download the
paper in PDF (Acrobat) format.

If you cannot use a PDF document, send email to ccowen@iupui.edu

To get a copy of the program that finds the numerical range
use your browser to save a (text file) copy of the source for this page.

% An Effective Algorithm for Calculating the Numerical Range
% Carl C. Cowen and Elad Harel
% Purdue University
%
% This script finds the numerical range of an n x n matrix by
% finding the real and imaginary parts of rotates of the matrix
% and finding the associated boundary point of that rotate by
% finding the largest eigenvalue of the real part and using the
% corresponding eigenvector's contribution to the numerical range.
% Multiplicity of the largest eigenvalue, as occurs in a normal
% matrix, is handled by plotting the end points of the corresponding
% segment in the boundary of the numerical range.
%
A=input('For what matrix do you want the numerical range? ')
nm=ceil(norm(A));
th=[0:.01:6.29];
k=1;
w=zeros(1,630);
for j=1:630
Ath=(exp(i*(-th(j))))*A;
Hth=(Ath+Ath')/2;
[r e]=eig(Hth);
e=real(diag(e));
m=max(e);
s=find(e==m);
if size(s,1)==1
w(k)=r(:,s)'*A*r(:,s);
%
% This is the point of the numerical range contributed by
% v_t=r(:,s) when the eigenspace of Hth (H_t) is one dimensional.
%
else
Kth=i*(Hth-Ath);
pKp=r(:,s)'*Kth*r(:,s);
%
% The matrix Q described above is r(:,s)
%
[rr ee]=eig(pKp);
ee=real(diag(ee));
mm=min(ee);
sm=find(ee==mm);
temp=rr(:,sm(:,1))'*r(:,s)'*A*r(:,s)*rr(:,sm(:,1));
w(k)=temp(1,1);
%
% This is the point of the numerical range contributed by
% v_t^- = r(:,s)*rr(:,sm(:,1))
%
k=k+1;
mM=max(ee);
sM=find(ee==mM);
temp=rr(:,sM(:,1))'*r(:,s)'*A*r(:,s)*rr(:,sM(:,1));
w(k)=temp(1,1);
%
% This is the point of the numerical range contributed by v_t^+
%
end
k=k+1;
end
fill(real(w),imag(w),'y')
axis([-nm,nm,-nm,nm])
axis('equal')

Back to Carl Cowen's Home Page