Matlab script to test NAClab's root finder 'uvFactor'

polytest.m

% Matlab script by Stefan Rickli
% 
% This Matlab script calls uvFactor of the NAClab toolbox by Zhonggang Zeng
% Get the Toolbox at http://homepages.neiu.edu/~zzeng/naclab.html

%% user variables
polydeg_min = 16;
polydeg_max = 23;

%% init
f1 = figure; s1 = subplot(1,1,1);
hold on;
f2 = figure; s2 = subplot(1,1,1);
hold on;
f3 = figure; s3 = subplot(1,1,1);

labels1 = {};
labels2 = {};
vertscat1 = [];
vertscat2 = [];
ratios = [];

%% root finding
for ii = polydeg_min:polydeg_max
    coeffs = fliplr(poly(1:ii));
    ratios(end+1) = coeffs(end)/coeffs(1);
    fprintf('\n\n\nRatio between constant coeff and highest order coeff: %d:%d = %d\n', coeffs(end), coeffs(1), coeffs(end)/coeffs(1));
    coeffs_perturbed = [coeffs(1:end-1),0];
    
    try
        f1 = uvFactor(coeffs,1e-9);
        f2 = uvFactor(coeffs_perturbed,1e-9);
        
        roots =             -f1(:,2)./f1(:,1);
        roots_perturbed =   -f2(:,2)./f2(:,1);
        
        scatter(s1, real(roots),imag(roots));
        scatter(s2, real(roots_perturbed),imag(roots_perturbed));
        
        labels1{end+1} = ['deg(p_2) = ', num2str(ii)];
        labels2{end+1} = ['deg(p_3) = ', num2str(ii)];
        
        vertscat1(end+1) = max(abs(imag(roots)));
        vertscat2(end+1) = max(abs(imag(roots_perturbed)));
    catch err
        disp(getReport(err,'extended'));
    end
end

%% plots
ylim(s1, [-0.15,0.15]);
title(s1, 'p_2(x) for different max poly degree');
legend(s1, labels1);

ylim(s2, [-0.15,0.15]);
title(s2, 'p_3(x) (perturbed) for different max poly degree');
legend(s2, labels2);

n = polydeg_min:polydeg_max;
plot(s3, n, -0.01*log10(abs(ratios)), n, vertscat1, n, vertscat2);
legend(s3, '-0.01*log10(abs(coefficient ratios))', 'vertical scattering p_2(x)', 'vertical scattering p_3(x)');