The program consists of 3 .m files:

- The main program: qrange.m
- The subprogram for plotting ellipses: ellipse.m
- The subprogram for plotting circles: circle.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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); % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%