Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Matlab script to test NAClab's root finder 'uvFactor'
% 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)');
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.