Computing the Nth Roots of a Number

In Matlab when we want to compute the square root of a number we use the 'sqrt()' command. For example, if z = 9 we find the square root of z by typing the

           >> sqrt(z)

command giving us the result

           ans =
               3

That result, however, provides only the 'principle' square root of z. As we know the number z = 9 has two square roots, i.e., +3 and -3. Likewise, if we want the cube root of z (also called the "3rd root of z") we type

           >> z^(1/3)

command giving us the principle result

           ans =
               2.0801

Again, the 3rd root of z has three solutions. They are:

           2.0801
           -1.0400 + 1.8014i
           -1.0400 - 1.8014i.
           (1)

The following Matlab code computes all of the n nth roots of the number(s) Num. If input argument Num has more than one element, then it must be a row vector. The elements in Num can be real-valued of complex-valued. Argument n is the order of the roots. If you want a complex-plane plot of the roots on Num that option, using the string 'Y', is also available. In such a plot the principle root is plotted as a blue square and the remaining roots are plotted with red circles.

For example, using the commands

           Num = 9;
           n=3;
           [nth_Roots] = roots_nth(Num, n, 'Y')

produces the results in Eq. (1) and the plot shown below.

Using the command

           [nth_Roots] = roots_nth(Num, n)

computes the n roots of Num, but no complex-plane plot is generated.

function [nth_Roots] = roots_nth(Num, n, Plot)
% ROOTS_NTH computes the nth roots of Num.
%
% Inputs:
%   Num is the number(s) for which the nth roots will be computed.
%     Num can be real-valued or complex-valued.
%     If Num has more than one element, then Num must be a row vector.
%   n must be a positive integer.
%   Plot is a string equal to 'Y' if a complex plane plot of 
%        the n complex roots of each element in Num is desired. 
%        Principle root is plotted as a blue square.  Remaining roots 
%        are red dots.  If input Plot string does not exist, then no 
%        plot is generated.
%
% Output:
%   nth_Roots is a column vector (with n rows) giving the nth roots of Num.
%
% Example:
%     clear, clc               % clear all variables
%     Num = [2+j*11, 3*exp(j*pi)]; % Two numbers
%     n = 3; % find the 3rd roots
%     [nth_Roots] = roots_nth(Num, n, 'Y')
%
%   Results are:
%       nth_Roots =
%
%          2.0000 + 1.0000i   0.7211 + 1.2490i
%         -1.8660 + 1.2321i  -1.4422 + 0.0000i
%         -0.1340 - 2.2321i   0.7211 - 1.2490i
%
%   [R. Lyons, Jan. 2013]
 
% Input argument check
if (nargin<2),
       disp('Not enough input arguments!!')
       return
    end;
    if (nargin==2)
        Plot == 'N'
    end
   
Num_of_Columns = max(size(Num)); % How many elements are in Num?
 
for Column_Loop = 1:Num_of_Columns;
    Mag(Column_Loop) = abs(Num(Column_Loop));
    Angle_Degrees(Column_Loop) = angle(Num(Column_Loop))*180/pi;
    for k = 1:n
        Root_Mag = Mag(Column_Loop)^(1/n);
        Root_Angle_Deg(k) = Angle_Degrees(Column_Loop)/n + (k-1)*360/n;
        nth_Roots(k,Column_Loop) = Root_Mag*exp(j*Root_Angle_Deg(k)*pi/180);
    end
    
    % Plot roots, on complex plane(s), if necessary
    if Plot == 'Y'
        figure
            % Plot unit circle
            Num_Pts = 200;     % Number of circle points
            Index = 0 : Num_Pts-1;
            Alpha = 1/Num_Pts;
            Points = Root_Mag*exp(j*2*pi*Alpha*Index);
            plot(Points, ':k','markersize',5)
            grid on
        hold on
 
        % Plot principle root (square blue dot) and list results in plot
            axis([-1.1*Root_Mag, 1.1*Root_Mag, -1.1*Root_Mag, 1.1*Root_Mag]);
            plot(Root_Mag*cos(Root_Angle_Deg(1)*pi/180),  ...
                Root_Mag*sin(Root_Angle_Deg(1)*pi/180), 'bs');
            text(-Root_Mag/2, 3*Root_Mag/4, ['Magnitudes = ',num2str(Root_Mag)]);
            text(-Root_Mag/2, 3*Root_Mag/4-Root_Mag/10, ...
                ['Principle Angle (deg.) = ',num2str(Root_Angle_Deg(1))]);
    
        % Plot remaining roots as red dots
            for k = 2:n
                plot(Root_Mag*cos(Root_Angle_Deg(k)*pi/180),  ...
                    Root_Mag*sin(Root_Angle_Deg(k)*pi/180), 'ro');
                text(-Root_Mag/2, 3*Root_Mag/4-k*Root_Mag/10, ...
                    ['             Angle ',num2str(k),' (deg.) = ', ...
                    num2str(Root_Angle_Deg(k))]);
            end % End k loop
            xlabel('Real'); ylabel('Imag.');
        hold off
    end % End 'if Plot == 'Y''
 
end % End Column_Loop