The program consists of 3 .m files:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function z = qrange(q,a,b,c)
%
% Function qrange(q,a,b,c)
%
% This is a Matlab program developed by Chi-Kwong Li to plot the ellipses
% and cirle that form the boundary of the q-numerical range of a 3 x 3
% diagonal matrix % diag(a,b,c) when 0 < q < 1, based on the result in
% [H. Nakazato, The boundary of therange of a contrained sesquilinear form,
% LAMA 40 (1996), 37-43.]
% Please feel free to use or modify it for your research or study.
% It would be nice if you would aknowledge the origin of the program
% if you use it.
%
% If a, b, c are not collinear, the following will be plotted:
% (i) three ellipses and the triangle whose vertices are their foci,
% (ii) the smallest circle containing the classical numerical range
% of diag(a,b,c) with its center translated to contain the 3 ellipses.
%
% The covex boundary of the q-numerical range will come from the ellipses
% and the circle. (See the paper of Nakazato for details.)
%
% To run this program, one also needs the programs:
% ellipse.m and circle.m developed by Chi-Kwong Li
%
%
if abs(a-b) + abs(b-c) == 0
z = 'You have a scalar matrix, whose q-numerical range is a singleton!';
else
%
% We identify the case when we have a scalar matrix above.
% In the next few cases, we have to generate circles and ellipses.
% Here n is the number of grid points when we plot the ellipses and circles.
%
n = 50;
A = [real(a) imag(a) 1/2 ; real(b) imag(b) 1/2; real(c) imag(c) 1/2];
if det(A) == 0
%
% Here we treat the case when a,b,c are collinear.
%
[X1, Y1] = ellipse(q*a, q*b, q, n);
[X2, Y2] = ellipse(q*b, q*c, q, n);
plot(X1,Y1,X2,Y2);
else
%
% Finally, we treat the case when a,b,c are in general positions.
% The next few commands compute the center cw and the radius r of the
% smallest circle that contiains a,b,c, and shift the center to qw.
%
w = (inv(A)*[(abs(a))^2 (abs(b))^2 (abs(c))^2]')/2;
r = sqrt(w(3) + w(1)^2 + w(2)^2);
s = q*r;
cw = w(1)+i*w(2);
X0 = q*w(1);
Y0 = q*w(2);
qw = X0+i*Y0;
%
[X1, Y1] = circle(qw, r, n);
[X2, Y2] = circle(qw, s, n);
%
% Here we determine the 3 ellipses needed to plot the boundary.
%
a1 = q*a;
b1 = q*b;
c1 = q*c;
[X3, Y3] = ellipse(a1,b1,q,n);
[X4, Y4] = ellipse(c1,b1,q,n);
[X5, Y5] = ellipse(a1,c1,q,n);
%
% Now determine W(A) and qW(A) shifted to the new center.
%
Z6 = [a b c a] + (q-1)*cw*ones(1,4);
X6 = real(Z6);
Y6 = imag(Z6);
Z7 = [a1 b1 c1 a1];
X7 = real(Z7);
Y7 = imag(Z7);
plot(X0,Y0,'x',X1,Y1,X3,Y3,X4,Y4,X5,Y5,X7,Y7);
end
% We choose to display big circle and its center, the 3 ellipses and the
% triangle whose vertice are their foci.
% One can choose to plot any of the figures determined above. For example:
%
% plot(X0,Y0,'x',X1,Y1,X2,Y2,X3,Y3,'g',X4,Y4,'g',X5,Y5,'g',X6,Y6,X7,Y7);
% pause
% plot(X1,Y1,X2,Y2,X3,Y3,X4,Y4,X5,Y5,X6,Y6);
%
axis('square');
title('The q-numerical range of A.');
xlabel('real axis');
ylabel('imaginary axis');
grid;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [X, Y] = ellipse(u,v,e,n)
%
% Function ellipse(u,v,e,n)
%
% This is a Matlab program developed by Chi-Kwong Li
% to plot an ellipse on the complex point with the complex numbers
% u and v as foci, and eccentricity 0 < e < 1, and 2n grid points.
%
z = (u-v)/2;
if z == 0
X = real(u);
Y = imag(u);
else
%
c = abs(z);
r = angle(z);
a = c/e;
b = sqrt(a^2 - c^2);
z1 = (u+v)/2;
p = real(z1);
q = imag(z1);
for k = 1:n
t = k*pi/n;
X(k) = a*cos(r)*cos(t)-b*sin(r)*sin(t) + p;
Y(k) = a*cos(t)*sin(r)+b*sin(t)*cos(r) + q;
X(n+k) = -a*cos(r)*cos(t)+b*sin(r)*sin(t) + p;
Y(n+k) = -a*cos(t)*sin(r)-b*sin(t)*cos(r) + q;
X(2*n+1) = X(1);
Y(2*n+1) = Y(1);
end
end
% plot(X,Y);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [X, Y] = circle(w,r,n)
%
% Function circle(w,r,n)
%
% This is a Matlab porgram developed by Chi-Kwong Li
% to plot a circle on the complex point with center w and radius r
% and 2n grid points.
%
%
w1 = real(w);
w2 = imag(w);
for k = 1:n
t = k*pi/n;
X(k) = w1 + r*cos(t);
Y(k) = w2 + r*sin(t);
X(n+k) = w1 - r*cos(t);
Y(n+k) = w2 - r*sin(t);
X(2*n+1) = X(1);
Y(2*n+1) = Y(1);
end
%
% plot(X,Y);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%