function [Poly_Rts, no_Real_Rts, no_Cmplx_Rts, MN, ...
Real_Rts, Real_Rts_Mult, ...
Poly_Rts_Cmplx_Abs_Distnt, Poly_Rts_Cmplx_Conj_Mult, ...
Poly_Rts_Cmplx_Angle_Distnt, Poly_Rts_Cmplx_Angle_Mult] = ...
cmplx_roots_sort_Det ( Poly_Rts, str_Rev )
%
% Sundar Krishnan (sundark100@yahoo.com)
% Release : R14, R13 (not R12 since R12 does not have strfind)
% Date : June-July 2005
%
% Functional Description :
% ------------------------
% Imagine a function that can pair the Complex Conjugate Roots
% of a Polynomial, and also tell us how many Roots are Real,
% and how many are Complex ; imagine further, that this function
% should also be able to give us the multiplicity of the roots ;
% and imagine further, that this function should be able to do all this
% even if not all the Complex Roots of the Polynomial are Conjugate Pairs.
% Information about multiplicity of roots is important when we need
% to raise a Polynomial to Fractional Powers.
%
% One (feeble) immediate answer could be that Matlab's standard cplxpair
% could do (part of) the functions.
% But, I have observed that as the degree of a Polynomial increases,
% the results obtained by Matlab's roots has so much variation
% that the standard cplxpair will result in an error !
% The standard Matlab roots () function does not return
% n exactly identical roots if the input Polynomial has
% n Repeated Roots, specially when n is above 10.
% So, firstly, we need to modify the tolerances of cplxpair.
% My cplxpair_Mod.m is a modified version of the Standard R14 Matlab's
% cplxpair.m.
%
% So, in addition to computing the number of Real and Complex Roots,
% and their multiplicities, we also need this new function
% to be able to perform the foll additional fns :
% a) optionally reverse the listing order of the roots
% b) List the Pure Real Roots in the beginning, and then the
% Complex (Conjugate-Paired) Roots for both the NORmal and the REVerse
% sorting cases.
%
% Such a function is this Programme cmplx_roots_sort_Det.m.
% This function makes calls to the lighter version cmplx_roots_sort.m,
% to Real_Roots_Multicty.m and to cplxpair_Mod.m.
% When cplxpair_Mod.m throws up an error as in the case of the
% presence of Non-Conjugate Complex Roots, the case is dealt with
% in the catch block of code by considering individual roots.
%
% An other Programme of mine, Poly_POWER.m, uses cmplx_roots_sort_Det.m
% to get the multiplicity of roots.
%
% We also know that multiple roots form the GCD of a Polynomial and it's
% Derivative. My Programme Poly_GCD.m calculates the GCD of 2 Polynomials.
% So, the product of the multiple roots of a Polynomial, as obtained from
% this Programme cmplx_roots_sort_Det.m, can be cross-checked
% with the GCD of the Polynomial and it's Derivative.
%
% If this suite of Programmes (in this zip file) don't work for you
% directly when downloaded into an independent directory, you could
% try loading these files in the Matlab dir : ...\toolbox\matlab\specfun
%
% &&&&&&&&&&&&&&&&&
%
% Author's Legal Disclaimer :
% ---------------------------
% I share my code, developed with lots of efforts, and often based
% on lots of studies, in good faith.
% Everyone is free to use this code at their own risk provided
% they include my name for this part in their works.
% However, I am still learning many facets of this subject ;
% there may be still some open unanswered technical questions for me
% on this subject and possibly in this programme.
% I am NOT responsible for any technical or legal complications
% directly or indirectly, arising out of any direct or indirect use
% of my codes or my interpretations.
% Any apparent TECHNICAL fault in interpreting any external references,
% or others' replies in any Public User Group Forums, is my own,
% not theirs. However, I am NOT LEGALLY responsible for any of my
% technical or other interpretation of others' comments in ANY Forum.
%
% &&&&&&&&&&&&&&&&&
%
% % Usage Egs :
% % For each case below, the "Test Code" is as follows :
% % Test Code :
% len_Poly_Rts_Cmplx = length ( Poly_Rts_Cmplx ) ;
% Poly_T_1 = [1] ;
% for i = 1 : len_Poly_Rts_Cmplx
% Poly_T_1 = conv ( Poly_T_1, [1, -Poly_Rts_Cmplx(i)] ) ;
% end
% roots_Poly_T_1 = roots ( Poly_T_1 )
% clc
% [Poly_Rts, no_Real_Rts, no_Cmplx_Rts, MN, ...
% Real_Rts, Real_Rts_Mult, ...
% Poly_Rts_Cmplx_Abs_Distnt, Poly_Rts_Cmplx_Conj_Mult, ...
% Poly_Rts_Cmplx_Angle_Distnt, Poly_Rts_Cmplx_Angle_Mult] = ...
% cmplx_roots_sort_Det ( roots_Poly_T_1 )
%
%
% Poly_T_2 = Poly_POWER (Poly_T_1, 2)
% roots_Poly_T_2 = roots ( Poly_T_2 )
% clc
% [Poly_Rts, no_Real_Rts, no_Cmplx_Rts, MN, ...
% Real_Rts, Real_Rts_Mult, ...
% Poly_Rts_Cmplx_Abs_Distnt, Poly_Rts_Cmplx_Conj_Mult, ...
% Poly_Rts_Cmplx_Angle_Distnt, Poly_Rts_Cmplx_Angle_Mult] = ...
% cmplx_roots_sort_Det ( roots_Poly_T_2 )
%
% % &&&&&&&&&&
%
% % Case 1 : (1 x 2, 1, 1) Real ; (3, 3, 2) Complex Pairs
% % OK for square - gives (1x4, 1x2, 1x2) Real, (6, 6, 4) Complex Pairs
% clc, KK = 0.1 ;
% Poly_Rts_Cmplx = KK * [ -2, -2, +7, -10 ...
% -13 + 4.02j, -13 - 4.02j, -13 + 4.02j, -13 - 4.02j, ...
% 4 + 13j, 4 - 13j, 4 + 13j, 4 - 13j, ...
% 4 + 13j, 4 - 13j, ...
% 2 - 1j, 2 + 1j, 2 - 1j, 2 + 1j, ...
% 2 - 1j, 2 + 1j ] ;
%
% % &&&&&&&&&&
%
% % Case 2a : With an additional root at 0 :
% % (1 x 2, 1, 1, 1) Real ; (3, 3, 2) Complex Pairs
% % OK for square - gives (1x4, 1x2, 1x2, 1x2) Real, (6, 6, 4) Complex Pairs
% clc, KK = 0.1 ;
% Poly_Rts_Cmplx = KK * [ -2, -2, +7, -10, 0, ...
% -13 + 4.02j, -13 - 4.02j, -13 + 4.02j, -13 - 4.02j, ...
% 4 + 13j, 4 - 13j, 4 + 13j, 4 - 13j, ...
% 4 + 13j, 4 - 13j, ...
% 2 - 1j, 2 + 1j, 2 - 1j, 2 + 1j, ...
% 2 - 1j, 2 + 1j ] ;
%
% % Case 2b : OK with KK = 10 also
%
% % &&&&&&&&&&
%
% % Case 3a : With multiple roots at 0 :
% % (1x2, 1, 1, 1x2) Real ; (2, 2, 2) Complex Pairs
% % OK for square - gives (1x4, 1x2, 1x2, 1x2) Real, (6, 6, 4) Complex Pairs
% clc, KK = 0.1 ;
% Poly_Rts_Cmplx = KK * [ -2, -2, +7, -10, 0, 0, ...
% -9 + 7.02j, -9 - 7.02j, -9 + 7.02j, -9 - 7.02j, ...
% 7 + 9j, 7 - 9j, 7 + 9j, 7 - 9j, ...
% 2 - 1j, 2 + 1j, 2 - 1j, 2 + 1j ] ;
%
% % Case 3b : OK with KK = 10 also
%
% % &&&&&&&&&&
%
% % Case 4 : (1 x 4) Real ; (3, 3, 2) Complex Pairs
% % OK for square - gives (1 x 8) Real, (6, 6, 4) Complex Pairs
% clc, KK = 1 ;
% Poly_Rts_Cmplx = KK * [ -2, -2, -2, -2, ...
% -13 + 4.02j, -13 - 4.02j, -13 + 4.02j, -13 - 4.02j, ...
% 4 + 13j, 4 - 13j, 4 + 13j, 4 - 13j, ...
% 4 + 13j, 4 - 13j, ...
% 2 - 1j, 2 + 1j, 2 - 1j, 2 + 1j, ...
% 2 - 1j, 2 + 1j ] ;
%
% % Case 4b : However, with KK = 0.1, the foll fails at the square operation.
% % The multiplicity of Real Roots obtained is 5 and 3 instead of 8.
%
% % &&&&&&&&&&
%
% % Case 5a : Multiplicity of 4 Real and 5 Complex Pairs
% clc, KK = 1 ;
% Poly_Rts_Cmplx = KK * [ -4, -4, -4, -4, ...
% 4 + 13j, 4 - 13j, 4 + 13j, 4 - 13j, ...
% 4 + 13j, 4 - 13j, ...
% 4 + 13j, 4 - 13j, 4 + 13j, 4 - 13j ] ;
%
% % Case 5b : With KK = 0.1, for the Square operation,
% % 2 Real Roots shift to complex form !
%
% % &&&&&&&&&&
%
% % Case 6a : Also here, for the Square operation,
% % 2 Real Roots shift to complex form !
% clc, KK = 0.1 ;
% Poly_Rts_Cmplx = KK * [ -5, -5, -5, -5, ...
% 4 + 13j, 4 - 13j, 4 + 13j, 4 - 13j, ...
% 4 + 13j, 4 - 13j, ...
% 4 + 13j, 4 - 13j, 4 + 13j, 4 - 13j ] ;
%
% % Case 6b, 6c : With Real Roots = -6 or -7 instead of -5 above,
% % for the Square operation, 1 Real Root shifts to complex form !
%
% % &&&&&&&&&&
%
% % Case 7 : cplxpair_Mod fails rather badly here :
% % Imag parts not conjugates (Test 2)
% clc, KK = 0.1 ;
% Poly_Rts_Cmplx = KK * [ -8, -8, -8, -8, ...
% 4 + 13j, 4 - 13j, 4 + 13j, 4 - 13j, ...
% 4 + 13j, 4 - 13j, ...
% 4 + 13j, 4 - 13j, 4 + 13j, 4 - 13j ] ;
%
% *******************
% fprintf ( '\nSTART of cmplx_roots_sort_Det \n' ) ; % for TESTing
% 1) Inits :
if nargin < 2
str_Rev = 'NOR' ; % Default is 'NORmal' ordering of roots.
end
rad_deg = 180 / pi ;
deg_rad = pi / 180 ;
% &&&&&&&&&&&&&&&&&
% 2) Convert the roots to a Col Vector :
Poly_Rts = Poly_Rts (:) ;
% &&&&&&&&&&&&&&&&&
% 3) Initial Common Procedure for NORMAL and Reciprocal Roots :
Poly_Rts_Real = [ ] ;
Poly_Rts_Cmplx = [ ] ;
len_roots = length ( Poly_Rts ) ;
if len_roots > 40
fprintf ( '\nlen_roots = %d \n', len_roots ) ;
wrnmsg_1 = strcat ( 'VERY HIGH DEG > 40 ;', ...
' RESULTS MAY NOT BE ACCURATE.' ) ;
warning (wrnmsg_1) ;
fprintf ( '\nDue to errors at high degrees in Matlab''s' ) ;
fprintf ( '\nconv and roots, this Programme may not work well' ) ;
fprintf ( '\nfor very high Polynomial Degrees like for eg, > 40.' ) ;
fprintf ( '\nPausing ... Press any key. \n' ) ;
pause
end
% 3b) Compute Empirical Abs Value of Roots, MN :
MN_mean = mean ( mean ( abs( Poly_Rts ) ) ) ;
MN_max = max ( max ( abs( Poly_Rts ) ) ) ;
MN = MN_mean * ( 1 + 0.5 * ( MN_max - MN_mean ) / MN_mean ) ;
if len_roots > 30 & MN < 0.8
fprintf ( '\nlen_roots = %d ; MN = %d \n', len_roots, MN ) ;
wrnmsg_2 = strcat ( 'HIGH DEG > 30 & VERY LOW ABS VALUE MN < 0.8 ;',...
' RESULTS MAY NOT BE ACCURATE.' ) ;
warning (wrnmsg_2) ;
fprintf ('\nlen_roots = %d ; Mean Abs Value = %d \n', len_roots, MN ) ;
fprintf ('\nPausing ... Press any key. \n' ) ;
pause
end
% &&&&&&&&&&&&&&&&&
% 4) Eliminate the "superfluous" negligible imag parts created by roots.
% For this, we use an empirical pure_real_fctr to test imag / real ratio
% and empirical imag_fctr to test the imag value.
% Empirically, it seems that the value of pure_real_fctr is based on :
% a) both len_roots and MN if MN <= 3.
% b) only MN if MN > 3.
% Calculate pure_real_fctr and imag_fctr :
% 4a) 1st level of calcs for pure_real_fctr based on len_roots :
pure_real_fctr = 0.005 ;
imag_fctr = 1 ;
if len_roots > 6
if len_roots <= 10
pure_real_fctr = 0.006 ;
elseif len_roots <= 20
pure_real_fctr = 0.007 ;
elseif len_roots <= 25
pure_real_fctr = 0.008 ;
elseif len_roots <= 30
pure_real_fctr = 0.01 ;
elseif len_roots <= 35
pure_real_fctr = 0.025 ;
elseif len_roots <= 40
pure_real_fctr = 0.0475 ;
else
pure_real_fctr = 0.06 ;
% Even though the Pgm will not work well for Deg > 40.
end
end % if len_roots > 6
% 4b) 2nd level of calcs for pure_real_fctr and imag_fctr based on MN :
if MN <= 0.25
pure_real_fctr = pure_real_fctr * 4 ;
imag_fctr = 0.25 ;
elseif MN <= 0.5
pure_real_fctr = pure_real_fctr * 2 ;
imag_fctr = 0.75 ;
elseif MN <= 0.75
pure_real_fctr = pure_real_fctr * 1.5 ;
imag_fctr = 1.5 ;
% For 0.75 < MN < 3, pure_real_fctr * 1 & imag_fctr = 1.
elseif MN >= 3
pure_real_fctr = 4 * 0.0006 * MN ;
imag_fctr = 4 ;
end
% &&&&&&&&&&&&&&&&&
% 5) The foll empirical constants are required to re-combine
% genuine multiple cases that may appear as different.
Diff_Ang_Alwd = 2.5 * deg_rad ; % 0.5 degrees = 0.0087 radians
% 2.5 degrees = 0.0436 radians
% 3 degrees = 0.0523 radians
if MN <= 0.25
Diff_Abs_Alwd = 0.1 ; % Based on 90 % match in Abs Values
elseif MN <= 0.5
Diff_Abs_Alwd = 0.07 ; % Based on 93 % match in Abs Values
elseif MN <= 1
Diff_Abs_Alwd = 0.06 ; % Based on 94 % match in Abs Values
elseif MN <= 5
Diff_Abs_Alwd = 0.055 ; % Based on 94.5 % match in Abs Values
elseif MN <= 10
Diff_Abs_Alwd = 0.0525 ; % Based on 94.75 % match in Abs Values
else % ie, for > 10
Diff_Abs_Alwd = 0.05 ; % Based on 95 % match in Abs Values
end
% &&&&&&&&&&&&&&&&&
% 6) Separate the Pure Real roots :
% In cplxpair_Mod.m too, the Pure Real roots are separated out.
% But, the value of tol used for that in cplxpair_Mod.m
% is significantly lower than what I have used here,
% or rather, what we actually need due to the inaccurate nature of
% roots and conv.
% Yet, I have not modified tol very drastically in cplxpair_Mod.m
% The main reason is that I did not want to trample the original code
% of cplxpair_Mod.m too drastically ; secondly, tol does not affect us
% much in the current context of working with cmplx_roots_sort_Det.m
% or cmplx_roots_sort.m since we are making calls to cplxpair_Mod.m
% only with inputs as Pure Complex Conjugates.
% So, the real challenge in cplxpair_Mod.m in the current context is
% tol_Real of cplxpair_Mod.m, not tol.
temp_Real_Rts = [ ] ;
temp_Cmplx_Rts = [ ] ;
for k = 1 : len_roots
% Note that here the code is slightly different from that
% followed in Real_Roots_Multicty for multiplicity of Real Roots.
% In this code section, we are only interested in
% seggeregating the real roots, not (yet) in counting the multiplicity.
if ( abs ( Poly_Rts (k) ) < 1000 * eps ) | ...
( abs ( imag ( Poly_Rts (k) ) / real ( Poly_Rts (k) ) ) ...
< pure_real_fctr ...
& abs ( imag ( Poly_Rts (k) ) ) < imag_fctr * pure_real_fctr )
% using "|" and "&" for compatibility with R12
Poly_Rts (k) = real ( Poly_Rts (k) ) ;
% The foll arrays temp_Real_Rts and temp_Cmplx_Rts are reqd
% in the CATCH block below for processing when we don't have
% Complex Roots in Conjugate Pairs, and hence, the call to
% cmplx_roots_sort, and thence to cplxpair_Mod will fail.
temp_Real_Rts = [ temp_Real_Rts, real( Poly_Rts (k) ) ] ;
else
temp_Cmplx_Rts = [ temp_Cmplx_Rts, Poly_Rts(k) ] ;
end
end
% &&&&&&&&&&&&&&&&&
% 7) 7 to 14 :
% We put the most common scenario - that of Complex Roots being in
% Conjugate Pairs - in the main try block.
% The reason is that the call to cplxpair_Mod, either via cmplx_roots_sort,
% or directly below, may fail if the Complex Roots are not in
% Conjugate Pairs
try % try - catch 1
% 8) After the true Real Roots have been seggreagted above,
% we need to sort to get the Real Roots pushed up in the list.
% Call the simple version of my sort programme :
[ Poly_Rts, no_Real_Rts, no_Cmplx_Rts ] = ...
cmplx_roots_sort ( Poly_Rts ) ;
% Note that cmplx_roots_sort will in turn call cplxpair_Mod,
% but with the input degree being only for the Complex Roots part
% (Poly_Rts_Cmplx).
% &&&&&&&&&&&&&&&&&
% 9) Separate the Real and Complex Roots :
for index = 1 : len_roots
if isreal (Poly_Rts(index))
Poly_Rts_Real = [ Poly_Rts_Real, Poly_Rts(index) ] ;
else
Poly_Rts_Cmplx = [ Poly_Rts_Cmplx, Poly_Rts(index) ] ;
end
end
% &&&&&&&&&&&&&&&&&
% 10) We are relying on a successful call of cplxpair (cplxpair_Mod)
% ie, this try block is for the case of Complex Roots being in
% Conjugate Pairs.
% Note that this direct call to cplxpair_Mod is also only for the
% Complex Roots part (Poly_Rts_Cmplx).
Poly_Rts_Cmplx = cplxpair_Mod(Poly_Rts_Cmplx) ;% if not OK, error here.
Poly_Rts_Cmplx_Abs = abs ( Poly_Rts_Cmplx ) ;
[Poly_Rts_Cmplx_Abs, ind_Cmplx] = sort ( Poly_Rts_Cmplx_Abs ) ;
[Poly_Rts_Real_Abs, ind_Real ] = sort ( Poly_Rts_Real ) ;
% Note : NOT by abs for Real Roots above.
% &&&&&&&&&&&&&&&&&
% 11) For 'NOR(MAL)'_Roots :
% Normal Sorted Order with Real first, then Complex Roots.
Poly_Rts_Real = Poly_Rts_Real(ind_Real) ;
Poly_Rts_Cmplx = Poly_Rts_Cmplx(ind_Cmplx) ;
% &&&&&&&&&&&&&&&&&
% Now, the "EXTENDED PART" with count (multiplicity)
% of Roots in Poly_Rts:
% 12) Count of Multiplicity of Real Roots :
% The implicit calculation of mplcity_fctr within Real_Roots_Multicty
% is good enough for normal cases ; however, when the variance
% of Real Roots is so large as in the case of the results of roots,
% we need to supply a larger mplcity_fctr ;
% mplcity_fctr = pure_real_fctr seems to work for many cases.
% (I tried mplcity_fctr = 0.5 * pure_real_fctr, but was not good enough
% for a case multiplicity of 4 Real Roots.)
[Real_Rts, Real_Rts_Mult] = ...
Real_Roots_Multicty ( Poly_Rts_Real, pure_real_fctr ) ;
% Note that we have already counted the no of Real and Complex Roots
% in the call to cmplx_roots_sort.
if no_Real_Rts ~= sum (Real_Rts_Mult)
errmsg_1 = 'no_Real_Rts and sum (Real_Rts_Mult) don''t match.' ;
no_Real_Rts
sum_Real_Rts_Mult = sum (Real_Rts_Mult)
error (errmsg_1) ;
end
% &&&&&&&&&&&&&&&&&
% 13) Now, the Multiplicity of Complex Roots :
% 13a) Get the multiplicity based on the "ABSolute" value
% of Complex Roots:
% Note that this multiplicity will show a factor of x 2 because of
% Conjugate Pairs.
% Even though multiplicity of Real Roots requires a higher
% mplcity_fctr = pure_real_fctr, the multiplicity of Absolute value of
% genuine Conjugate-Pair Complex Roots seems to work with the
% intrinsic value of mplcity_fctr as provided in Real_Roots_Multicty.
% So, it seems that we don't need to "boost" mplcity_fctr
% while checking for multiplicity of Conjugate-paired Complex Roots ?
% However, to be sure, I have boosted here too by 0.25 * pure_real_fctr.
[Poly_Rts_Cmplx_Abs_Distnt, Poly_Rts_Cmplx_Abs_Mult] = ...
Real_Roots_Multicty ( Poly_Rts_Cmplx_Abs, 0.25 * pure_real_fctr ) ;
if no_Cmplx_Rts ~= sum (Poly_Rts_Cmplx_Abs_Mult)
errmsg_2 = strcat ( 'no_Cmplx_Rts and', ...
'sum (Poly_Rts_Cmplx_Abs_Mult) don''t match.');
no_Cmplx_Rts
sum_Rts_Cmplx_Abs_Mult = sum (Poly_Rts_Cmplx_Abs_Mult)
error (errmsg_2) ;
end
if any ( mod ( Poly_Rts_Cmplx_Abs_Mult, 2 ) ) == 1
errmsg_3='The Multiplicity of Complex Roots are NOT ALL Even Nos.';
error (errmsg_3) ;
end
% ^^^^^^^^^
% 13b) The Angle Part
Poly_Rts_Cmplx_Angle_Distnt = [ ] ;
Poly_Rts_Cmplx_Angle_Mult = [ ] ;
Poly_Rts_Cmplx_Conj_Mult = [ ] ;
for i = 1 : length ( Poly_Rts_Cmplx_Abs_Mult )
% Get the individual angles within a set of Poly_Rts_Cmplx_Abs_Mult
ang_ind_roots = angle ( Poly_Rts_Cmplx ( ...
sum ( Poly_Rts_Cmplx_Abs_Mult ( 1 : i-1 ) ) + 1 : ...
sum ( Poly_Rts_Cmplx_Abs_Mult ( 1 : i ) ) ) ) ;
% Get the multiplicity of angles in the above set of angles.
% Note that we need to step by 2 in the try block
% because of Conjugate Pairs for Complex Roots.
mplcity_fctr = min ( 0.005, ...
0.06 * mean ( mean ( abs( ang_ind_roots ) ) ) ) ;
[ang_ind_roots, ang_ind_roots_Mult] = ...
Real_Roots_Multicty ...
( ang_ind_roots (2 : 2 : end), mplcity_fctr ) ;
len_ang_ind_roots_Mult = length(ang_ind_roots_Mult) ;
Poly_Rts_Cmplx_Angle_Distnt = [ Poly_Rts_Cmplx_Angle_Distnt, ...
ang_ind_roots ] ;
Poly_Rts_Cmplx_Angle_Mult = [ Poly_Rts_Cmplx_Angle_Mult, ...
ang_ind_roots_Mult ] ;
end
% ^^^^^^^^^
% 13c) The Kernel Part of the Counting of Complex Roots :
% Note that for the same Poly_Rts_Cmplx_Abs_Distnt, we could have
% different angle pairs. So, there is a need to further split the above
% Poly_Rts_Cmplx_Abs_Mult value.
i = 1 ;
while i <= length ( Poly_Rts_Cmplx_Abs_Mult )
% Initial No of Roots within Individual Poly_Rts_Cmplx_Abs_Mult :
no_ind_roots = Poly_Rts_Cmplx_Abs_Mult(i) ;
% No need to split further if the multiplicity is only 2.
if no_ind_roots >= 4
% Get the individual angles
% within a set of Poly_Rts_Cmplx_Abs_Mult
ang_ind_roots = angle ( Poly_Rts_Cmplx ( ...
sum ( Poly_Rts_Cmplx_Abs_Mult ( 1 : i-1 ) ) + 1 : ...
sum ( Poly_Rts_Cmplx_Abs_Mult ( 1 : i ) ) ) ) ;
% Get the multiplicity of angles in the above set of angles.
% Note that we need to step by 2 in the try block
% because of Conjugate Pairs for Complex Roots.
mplcity_fctr = min ( 0.005, ...
0.06 * mean ( mean ( abs( ang_ind_roots ) ) ) ) ;
[ang_ind_roots, ang_ind_roots_Mult] = ...
Real_Roots_Multicty ( ang_ind_roots (2 : 2 : end), ...
mplcity_fctr ) ;
len_ang_ind_roots_Mult = length(ang_ind_roots_Mult) ;
% If ang_ind_roots_Mult has more than 1 column, it means that
% we need to split.
if len_ang_ind_roots_Mult > 1
% Current Starting Length BEFORE Split
% (under the current i) :
len_Abs_Distnt = length ( Poly_Rts_Cmplx_Abs_Distnt ) ;
% len_Abs_Distnt and len_Abs_Mult should always match.
% If they don't, then, it means that there may be
% some problem with the Programme
len_Abs_Mult = length ( Poly_Rts_Cmplx_Abs_Mult ) ;
if len_Abs_Distnt ~= len_Abs_Mult
len_Abs_Distnt, len_Abs_Mult
errmsg_4A = ...
strcat ( 'len_Abs_Distnt and len_Abs_Mult', ...
' do not match - Check Pt 4A.' ) ;
error ( errmsg_4A ) ;
end
% The index change of the terms BEYOND
for k = len_Abs_Distnt : -1 : i + 1
Poly_Rts_Cmplx_Abs_Distnt ...
( k + len_ang_ind_roots_Mult - 1) = ...
Poly_Rts_Cmplx_Abs_Distnt ( k ) ;
end
% To insert the multiplicity of the Current term
for kk = 1 : len_ang_ind_roots_Mult - 1
Poly_Rts_Cmplx_Abs_Distnt (i + kk) = ...
Poly_Rts_Cmplx_Abs_Distnt (i) ;
end
% The index change of the terms BEYOND
for r = len_Abs_Mult : -1 : i + 1
Poly_Rts_Cmplx_Abs_Mult ...
( r + len_ang_ind_roots_Mult - 1) = ...
Poly_Rts_Cmplx_Abs_Mult ( r ) ;
end
% To insert the multiplicity of the Current term
Poly_Rts_Cmplx_Abs_Mult (i) = 2 * ang_ind_roots_Mult (1) ;
for rr = 2 : len_ang_ind_roots_Mult
Poly_Rts_Cmplx_Abs_Mult (i + rr - 1) = ...
2 * ang_ind_roots_Mult (rr) ; % 2 for Conj Pairs
end
len_Abs_Distnt = length ( Poly_Rts_Cmplx_Abs_Distnt ) ;
len_Abs_Mult = length ( Poly_Rts_Cmplx_Abs_Mult ) ;
if len_Abs_Distnt ~= len_Abs_Mult
len_Abs_Distnt, len_Abs_Mult
errmsg_5A = ...
strcat ( 'len_Abs_Distnt and len_Abs_Mult', ...
' do not match - Check Pt 5A.' ) ;
error ( errmsg_5A ) ;
end
end % if len_ang_ind_roots_Mult > 1
end % if no_ind_roots >= 4
i = i + 1 ;
end % while i <= length ( Poly_Rts_Cmplx_Abs_Mult )
% ^^^^^^^^^
% 13d) The multiplicity of Conjugate Roots is half :
Poly_Rts_Cmplx_Conj_Mult = Poly_Rts_Cmplx_Abs_Mult / 2 ;
% ^^^^^^^^^
% 13e-1) Poly_Rts_Cmplx_Conj_Mult and Poly_Rts_Cmplx_Angle_Mult
% should match.
if any ( Poly_Rts_Cmplx_Conj_Mult - Poly_Rts_Cmplx_Angle_Mult )
Poly_Rts_Cmplx_Conj_Mult, Poly_Rts_Cmplx_Angle_Mult
errmsg_6A = strcat ( 'Poly_Rts_Cmplx_Conj_Mult and', ...
' Poly_Rts_Cmplx_Angle_Mult', ...
' do not match - Check Pt 6A.' ) ;
error ( errmsg_6A ) ;
end
% 13e-2)A Redundant Check :
% Lengths of Poly_Rts_Cmplx_Abs_Distnt and Poly_Rts_Cmplx_Angle_Distnt
% should match.
if length (Poly_Rts_Cmplx_Abs_Distnt) ~= ...
length (Poly_Rts_Cmplx_Angle_Distnt)
Poly_Rts_Cmplx_Abs_Distnt, Poly_Rts_Cmplx_Angle_Distnt
errmsg_7A = strcat ( 'Lengths of Poly_Rts_Cmplx_Abs_Distnt', ...
' and Poly_Rts_Cmplx_Angle_Distnt', ...
' do not match - Check Pt 7A.' ) ;
error ( errmsg_7A ) ;
end
% ^^^^^^^^^
% 13f) Due to the errors creeping in from cplxpair, we may still have
% genuine multiple cases shown separately.
% So, we need to iron them out, and show true multiplicity.
% For this, we use Diff_Abs_Alwd (within 95 % t0 90 %)
% and Diff_Ang_Alwd (within 2.5 degrees) in the lists above
% to identify genuine multiples.
i = 1 ;
while i < length (Poly_Rts_Cmplx_Abs_Distnt)
% for will NOT work : ie, instead of :
% for i = 1 : length (Poly_Rts_Cmplx_Abs_Distnt) - 1
j = i + 1 ;
while j <= length (Poly_Rts_Cmplx_Abs_Distnt)
if abs ( Poly_Rts_Cmplx_Abs_Distnt(i) - ...
Poly_Rts_Cmplx_Abs_Distnt(j) ) ...
/ abs ( Poly_Rts_Cmplx_Abs_Distnt(i) ) ...
<= Diff_Abs_Alwd ...
& abs ( Poly_Rts_Cmplx_Angle_Distnt(i) - ...
Poly_Rts_Cmplx_Angle_Distnt(j) ) ...
<= Diff_Ang_Alwd
Poly_Rts_Cmplx_Conj_Mult(i) = ...
Poly_Rts_Cmplx_Conj_Mult(i) + ...
Poly_Rts_Cmplx_Conj_Mult(j) ;
Poly_Rts_Cmplx_Conj_Mult(j) = [ ] ;
Poly_Rts_Cmplx_Abs_Distnt(j) = [ ] ;
Poly_Rts_Cmplx_Angle_Distnt(j) = [ ] ;
else
% Note Tricky here !
% We need to move "forward" in j
% ONLY if we have NOT found a match.
% If a match is found we have effectively moved forward
% because of "shrinking".
j = j + 1 ;
end % if abs ...
end % while j <= length (Poly_Rts_Cmplx_Abs_Distnt)
i = i + 1 ;
end % while i < length (Poly_Rts_Cmplx_Abs_Distnt)
% Since we have already checked above that Poly_Rts_Cmplx_Conj_Mult
% and Poly_Rts_Cmplx_Angle_Mult match, we simply assign as :
Poly_Rts_Cmplx_Angle_Mult = Poly_Rts_Cmplx_Conj_Mult ;
% &&&&&&&&&&&&&&&&&
% 14) Now, (re-)construct Poly_Rts_Real, Poly_Rts_Cmplx,
% and thus Poly_Rts.
% 14 a) Real Roots first (The Real Roots portion could be common
% for try and catch blocks) :
Poly_Rts_Real = [ ] ;
for i = 1 : length (Real_Rts)
Poly_Rts_Real = [ Poly_Rts_Real, ...
Real_Rts(i) * ones( 1, Real_Rts_Mult(i) ) ] ;
end
% ^^^^^^^^^
% 14 b) Complex Roots next (Note that the Complex Roots portion
% cannot be common for try and catch blocks) :
Poly_Rts_Cmplx = [ ] ;
for i = 1 : length (Poly_Rts_Cmplx_Abs_Distnt)
% Why am I not using exp (j * ...) ?
% Ans : Because j is already used above !
temp = [ ] ;
cs = cos ( Poly_Rts_Cmplx_Angle_Distnt(i) ) ;
sn = sin ( Poly_Rts_Cmplx_Angle_Distnt(i) ) ;
temp = Poly_Rts_Cmplx_Abs_Distnt(i) * complex( cs, sn ) ;
temp = [ temp, conj(temp) ] ; % this way only for Conj Roots
Poly_Rts_Cmplx = [ Poly_Rts_Cmplx, ...
repmat( temp, 1, Poly_Rts_Cmplx_Angle_Mult(i) ) ] ;
end
% &&&&&&&&&&&&&&&&&
catch % try - catch 1
fprintf ( '\n\n **** CATCH Block in cmplx_roots_sort_Det **** \n' ) ;
fprintf ( '\n ie, All Complex Roots are NOT as Conjuagte Pairs \n\n' );
% 15) 15 to 21 :
% If the error was in the call in cplxpair_Mod
% (in the call through cmplx_roots_sort), we need to reconcile to
% computing the multiplicity of roots by taking each individual root.
% One reason could be that we could have a mixed case :
% some Complex Conjugates Pairs, some Complex Non-Conjugates
% and some Real Roots etc.
if ~isempty ( strfind ( lasterr, 'cplxpair_Mod' ) ) ...
|| ~isempty ( strfind ( lasterr, 'cplxpair_mod' ) ) % Reqd for R13 !!!
% strfind not in R12 !
% Here, we make use of temp_Real_Rts and temp_Cmplx_Rts
% calculated above - precisely for this scenario.
% 16) Get Poly_Rts_Real_Abs and Poly_Rts_Cmplx_Abs and
[Poly_Rts_Real_Abs, ind_Real] = sort( temp_Real_Rts ) ;
% Note : NOT by abs for Real Roots above.
Poly_Rts_Cmplx_Abs = abs ( temp_Cmplx_Rts ) ;
[Poly_Rts_Cmplx_Abs, ind_Cmplx] = sort ( Poly_Rts_Cmplx_Abs ) ;
% &&&&&&&&&&&&&&&&&
% 17) For 'NOR(MAL)'_Roots :
% Normal Sorted Order with Real first, then Complex Roots.
Poly_Rts_Real = temp_Real_Rts (ind_Real) ;
Poly_Rts_Cmplx = temp_Cmplx_Rts (ind_Cmplx) ;
% &&&&&&&&&&&&&&&&&
% 18) Count of Multiplicity of Real Roots :
[Real_Rts, Real_Rts_Mult] = ...
Real_Roots_Multicty ( Poly_Rts_Real, pure_real_fctr ) ;
% In the Catch Block, we do not have no_Real_Rts and no_Cmplx_Rts
% from cmplx_roots_sort as in the Try Block.
no_Real_Rts = sum (Real_Rts_Mult) ;
% &&&&&&&&&&&&&&&&&
% 19) Now, the Multiplicity of Complex Roots :
% 19a) Get the multiplicity based on the "ABSolute" value
% of Complex Roots:
[Poly_Rts_Cmplx_Abs_Distnt, Poly_Rts_Cmplx_Abs_Mult] = ...
Real_Roots_Multicty ( Poly_Rts_Cmplx_Abs, ...
0.25 * pure_real_fctr ) ;
no_Cmplx_Rts = sum (Poly_Rts_Cmplx_Abs_Mult) ;
% ^^^^^^^^^
% 19b) The Angle Part
Poly_Rts_Cmplx_Angle_Distnt = [ ] ;
Poly_Rts_Cmplx_Angle_Mult = [ ] ;
Poly_Rts_Cmplx_Conj_Mult = [ ] ;
for i = 1 : length ( Poly_Rts_Cmplx_Abs_Mult )
% Get the individual angles within a set
% of Poly_Rts_Cmplx_Abs_Mult
ang_ind_roots = angle ( Poly_Rts_Cmplx ( ...
sum ( Poly_Rts_Cmplx_Abs_Mult ( 1 : i-1 ) ) + 1 : ...
sum ( Poly_Rts_Cmplx_Abs_Mult ( 1 : i ) ) ) ) ;
% Get the multiplicity of angles in the above set of angles.
% Note that here, we need to step by 1 only because
% we are not dealing with Conjugate Pairs for Complex Roots.
mplcity_fctr = min ( 0.005, ...
0.06 * mean ( mean ( abs( ang_ind_roots ) ) ) ) ;
[ang_ind_roots, ang_ind_roots_Mult] = ...
Real_Roots_Multicty ...
( ang_ind_roots (1 : 1 : end), mplcity_fctr ) ;
len_ang_ind_roots_Mult = length (ang_ind_roots_Mult) ;
Poly_Rts_Cmplx_Angle_Distnt = [ Poly_Rts_Cmplx_Angle_Distnt,...
ang_ind_roots ] ;
Poly_Rts_Cmplx_Angle_Mult = [ Poly_Rts_Cmplx_Angle_Mult, ...
ang_ind_roots_Mult ] ;
end % for i = 1 : length ( Poly_Rts_Cmplx_Abs_Mult )
% ^^^^^^^^^
% 19c) The Kernel Part of the Counting of Complex Roots :
% Note that for the same Poly_Rts_Cmplx_Abs_Distnt, we could have
% different angle pairs.
% So, there is a need to further split the above
% Poly_Rts_Cmplx_Abs_Mult value.
i = 1 ;
while i <= length ( Poly_Rts_Cmplx_Abs_Mult )
% Initial No of Roots within Individual values of
% Poly_Rts_Cmplx_Abs_Mult :
no_ind_roots = Poly_Rts_Cmplx_Abs_Mult(i) ;
% No need to split further if the multiplicity is only 1.
if no_ind_roots >= 2
% Get the individual angles within a set of
% Poly_Rts_Cmplx_Abs_Mult :
ang_ind_roots = angle ( Poly_Rts_Cmplx ( ...
sum ( Poly_Rts_Cmplx_Abs_Mult ( 1 : i-1 ) ) + 1 : ...
sum ( Poly_Rts_Cmplx_Abs_Mult ( 1 : i ) ) ) ) ;
% Get the multiplicity of angles in the above set
% of angles.
% Note that here, we need to step by 1 only because we
% are not dealing with Conjugate Pairs for Complex Roots.
mplcity_fctr = min ( 0.005, ...
0.06 * mean ( mean ( abs( ang_ind_roots ) ) ) ) ;
[ang_ind_roots, ang_ind_roots_Mult] = ...
Real_Roots_Multicty ( ang_ind_roots (1 : 1 : end), ...
mplcity_fctr ) ;
len_ang_ind_roots_Mult = length (ang_ind_roots_Mult) ;
% If ang_ind_roots_Mult has more than 1 column, it means
% that we need to split.
if len_ang_ind_roots_Mult > 1
% Current Starting Length BEFORE Split
% (under the current i) :
len_Abs_Distnt = length ( Poly_Rts_Cmplx_Abs_Distnt ) ;
% len_Abs_Distnt and len_Abs_Mult should always match.
% If they don't, then, it means that there may be
% some problem with the Programme
len_Abs_Mult = length ( Poly_Rts_Cmplx_Abs_Mult ) ;
if len_Abs_Distnt ~= len_Abs_Mult
len_Abs_Distnt, len_Abs_Mult
errmsg_1B = ...
strcat ( 'len_Abs_Distnt and len_Abs_Mult', ...
' do not match - Check Pt 1B.' ) ;
error ( errmsg_1B ) ;
end
% The index change of the terms BEYOND
for k = len_Abs_Distnt : -1 : i + 1
Poly_Rts_Cmplx_Abs_Distnt ...
( k + len_ang_ind_roots_Mult - 1) = ...
Poly_Rts_Cmplx_Abs_Distnt ( k ) ;
end
% To insert the multiplicity of the Current term
for kk = 1 : len_ang_ind_roots_Mult - 1
Poly_Rts_Cmplx_Abs_Distnt (i + kk) = ...
Poly_Rts_Cmplx_Abs_Distnt (i) ;
end
% The index change of the terms BEYOND
for r = len_Abs_Mult : -1 : i + 1
Poly_Rts_Cmplx_Abs_Mult ...
( r + len_ang_ind_roots_Mult - 1) = ...
Poly_Rts_Cmplx_Abs_Mult ( r ) ;
end
Poly_Rts_Cmplx_Abs_Mult(i) = ang_ind_roots_Mult(1) ;
% To insert the multiplicity of the Current term
% This section took quite some time.
for rr = 2 : len_ang_ind_roots_Mult
Poly_Rts_Cmplx_Abs_Mult (i + rr - 1) = ...
ang_ind_roots_Mult (rr) ;
end
len_Abs_Distnt = length ( Poly_Rts_Cmplx_Abs_Distnt ) ;
len_Abs_Mult = length ( Poly_Rts_Cmplx_Abs_Mult ) ;
if len_Abs_Distnt ~= len_Abs_Mult
len_Abs_Distnt, len_Abs_Mult
errmsg_2B = ...
strcat ( 'len_Abs_Distnt and len_Abs_Mult', ...
' do not match - Check Pt 2B.' ) ;
error ( errmsg_2B ) ;
end
end % if len_ang_ind_roots_Mult > 1
end % if no_ind_roots >= 2
i = i + 1 ;
end % while i <= length ( Poly_Rts_Cmplx_Abs_Mult )
% ^^^^^^^^^
% 19d-1) Poly_Rts_Cmplx_Abs_Mult and Poly_Rts_Cmplx_Angle_Mult
% should match.
if any ( Poly_Rts_Cmplx_Abs_Mult - Poly_Rts_Cmplx_Angle_Mult )
Poly_Rts_Cmplx_Abs_Mult, Poly_Rts_Cmplx_Angle_Mult
errmsg_3B = strcat ( 'Poly_Rts_Cmplx_Abs_Mult and', ...
' Poly_Rts_Cmplx_Angle_Mult', ...
' do not match - Check Pt 3B.' ) ;
error ( errmsg_3B ) ;
end
% 19d-2) A Redundant Check :
% Lengths of Poly_Rts_Cmplx_Abs_Distnt and
% Poly_Rts_Cmplx_Angle_Distnt should match.
if length (Poly_Rts_Cmplx_Abs_Distnt) ~= ...
length (Poly_Rts_Cmplx_Angle_Distnt)
Poly_Rts_Cmplx_Abs_Distnt, Poly_Rts_Cmplx_Angle_Distnt
errmsg_4B = strcat ( 'Lengths of Poly_Rts_Cmplx_Abs_Distnt',...
' and Poly_Rts_Cmplx_Angle_Distnt', ...
' do not match - Check Pt 4B.' ) ;
error ( errmsg_4B ) ;
end
% ^^^^^^^^^
% 19e) Let us use Diff_Abs_Alwd (within 95 % t0 90 %)
% and Diff_Ang_Alwd (within 2.5 degrees) in the lists above
% to find genuine multiples.
i = 1 ;
while i < length (Poly_Rts_Cmplx_Abs_Distnt)
j = i + 1 ;
while j <= length (Poly_Rts_Cmplx_Abs_Distnt)
if abs ( Poly_Rts_Cmplx_Abs_Distnt(i) - ...
Poly_Rts_Cmplx_Abs_Distnt(j) ) ...
/ abs ( Poly_Rts_Cmplx_Abs_Distnt(i) ) ...
<= Diff_Abs_Alwd ...
& abs ( Poly_Rts_Cmplx_Angle_Distnt(i) - ...
Poly_Rts_Cmplx_Angle_Distnt(j) ) ...
<= Diff_Ang_Alwd
Poly_Rts_Cmplx_Angle_Mult(i) = ...
Poly_Rts_Cmplx_Angle_Mult(i) + ...
Poly_Rts_Cmplx_Angle_Mult(j) ;
Poly_Rts_Cmplx_Angle_Mult(j) = [ ] ;
Poly_Rts_Cmplx_Abs_Distnt(j) = [ ] ;
Poly_Rts_Cmplx_Angle_Distnt(j) = [ ] ;
else
% Note Tricky here !
% We need to move "forward" in j
% ONLY if we have NOT found a match.
% If a match is found we have effectively moved
% forward because of "shrinking".
j = j + 1 ; %
end % if abs ...
end % while j <= length (Poly_Rts_Cmplx_Abs_Distnt)
i = i + 1 ;
end % while i < length (Poly_Rts_Cmplx_Abs_Distnt)
% &&&&&&&&&&&&&&&&&
% 20) Now, (re-)construct Poly_Rts_Real, Poly_Rts_Cmplx,
% and thus Poly_Rts.
% 20 a) Real Roots first (The Real Roots portion could be common
% for try and catch blocks) :
Poly_Rts_Real = [ ] ;
for i = 1 : length (Real_Rts)
Poly_Rts_Real = [ Poly_Rts_Real, ...
Real_Rts(i) * ones( 1, Real_Rts_Mult(i) ) ] ;
end
% ^^^^^^^^^
% 20 b) Complex Roots next (Note that the Complex Roots portion
% cannot be common for try and catch blocks) :
Poly_Rts_Cmplx = [ ] ;
for i = 1 : length (Poly_Rts_Cmplx_Abs_Distnt)
% Why am I not using exp (j * ...) ?
% Ans : Because j is already used above !
temp = [ ] ;
cs = cos ( Poly_Rts_Cmplx_Angle_Distnt(i) ) ;
sn = sin ( Poly_Rts_Cmplx_Angle_Distnt(i) ) ;
temp = Poly_Rts_Cmplx_Abs_Distnt(i) * complex( cs, sn ) ;
% This is NOT for the catch block !
% temp = [ temp, conj(temp) ] ;
Poly_Rts_Cmplx = ...
[Poly_Rts_Cmplx, repmat( temp, 1, Poly_Rts_Cmplx_Angle_Mult(i) )] ;
end
% &&&&&&&&&&&&&&&&&
else % If the error in the try block was some error other than
% that which originated from cplxpair_Mod :
% 21) Exit with that other different error in try block.
error (lasterr) ;
end % if ~isempty ( strfind ( lasterr, 'MATLAB:cplxpair_Mod' ) )
end % catch of try - catch 1
% &&&&&&&&&&&&&&&&&
% 22) Reverse Sorting :
if ~any ( str_Rev ~= 'REV' )
str_Rev
% For Reciprocal ie, 'REV(ERSE)' Roots :
% Here is where we reverse the sorted order
% within Real and within Complex, respectively :
Poly_Rts_Real = fliplr( Poly_Rts_Real ) ;
Poly_Rts_Cmplx = fliplr( Poly_Rts_Cmplx ) ;
Real_Rts = fliplr ( Real_Rts ) ;
Real_Rts_Mult = fliplr ( Real_Rts_Mult ) ;
Poly_Rts_Cmplx_Abs_Distnt = fliplr ( Poly_Rts_Cmplx_Abs_Distnt ) ;
% Poly_Rts_Cmplx_Abs_Mult = fliplr ( Poly_Rts_Cmplx_Abs_Mult ) ;
Poly_Rts_Cmplx_Angle_Distnt = fliplr ( Poly_Rts_Cmplx_Angle_Distnt ) ;
Poly_Rts_Cmplx_Angle_Mult = fliplr ( Poly_Rts_Cmplx_Angle_Mult ) ;
% Note that for Non-Conjugate Complex Roots, ie, for code in the
% Catch Block, Poly_Rts_Cmplx_Conj_Mult is an empty matrix.
Poly_Rts_Cmplx_Conj_Mult = fliplr ( Poly_Rts_Cmplx_Conj_Mult ) ;
end
% &&&&&&&&&&&&&&&&&
% 23) Finally, the overall Poly_Rts :
% This is common to both NOR and REV orders (respectively) :
Poly_Rts = [ Poly_Rts_Real, Poly_Rts_Cmplx ] ;
% &&&&&&&&&&&&&&&&&
% fprintf ( '\nEND of cmplx_roots_sort_Det \n' ) ; % for TESTing
% *******************
