%% FINDING THE SMALLEST POSITIVE ROOT OF A POLYNOMIAL EQUATION FOR A
% In a previous blog at autarkaw.wordpress.com (June 10), we set up a
% polynomial equation that would allow us to find the longest mast 
% that can be setup before it buckles under its own weight.
% In this blog, we will find the root of this equation.
% The problem is given at
% http://numericalmethods.eng.usf.edu/blog/mast_of_length.pdf
% and we are solving Exercise 1 of the pdf file.

%% TOPIC
% Smallest positive root of a polynomial equations with infinite terms

%% SUMMARY

% Language : Matlab 2008a; 
% Authors : Autar Kaw; 
% Mfile available at
% http://numericalmethods.eng.usf.edu/blog/rootinfinite.m; 
% Last Revised : July 4, 2010
% Abstract: This program shows you how to find the smallest positive 
% real root of a polynomial equations with infinite terms
clc
clear all

%% INTRODUCTION

disp('ABSTRACT')
disp('   This program shows you how to')
disp('   find the smallest positive real root')
disp('   of a polynomial equations with infinite terms')
disp(' ')
disp('AUTHOR')
disp('   Autar K Kaw of http://autarkaw.wordpress.com')
disp(' ')
disp('MFILE SOURCE')
disp('   http://numericalmethods.eng.usf.edu/blog/rootsinfinite.m')
disp('  ')
disp('PROBLEM STATEMENT')
disp('   http://numericalmethods.eng.usf.edu/blog/mast_of_length.pdf')
disp(' ')
disp('LAST REVISED')
disp('   July 4, 2010')
disp(' ')

%% INPUTS
% prespecified tolerance, eps
eps=0.000001;
% maximum number of terms of polynomial
nmax=100;

%% DISPLAYING INPUTS

disp('INPUTS')
fprintf(' The max number of terms of the polynomial chosen, nmax= %g',nmax)
fprintf('\n The prespecified tolerance, eps= %g',eps)
disp('  ')
disp('  ')

%% CODE
for N=2:1:nmax
    % The above looop is to see how many terms we should take of the 
    % infinite polynomial
  aa(1)=-3.0/8.0;
    % Setting up the polynomial via recursive relations
for i=2:1:N
  aa(i)=-3*aa(i-1)/(4*i*(3*i-1));
end
% Since it is a polynomial of order N, 
% there are N+1 coefficients
% To set up the polynomial for MATLAB
% N+1 th coefficient is the constant term
% N th coefficient is the term of order 1
% and so on till 1st coefficient is of order N
bb(N+1)=1;
for i=1:1:N
    bb(N-i+1)=aa(i);
end

% Finding all the roots of the Nth order polynomial
abc=roots(bb);

% Finding the first real positive root so that it 
% would be used as the starting minimum value available 
for i=1:1:N
    if isreal(abc(i))==true & abc(i)>0
        minval=abc(i);
        break;
    end
end

% Finding the smallest positive real root

for i=1:1:N
    if isreal(abc(i))==true & abc(i)>0
        if (abc(i) < minval)
            minval=abc(i);
        end
    end
end
% Checking if prespecified tolerance is met
if N>2
    absea=abs((minval-previous)/minval)*100;
    if absea<=eps 
        terms_needed=N;
        break;
    end
end
previous=minval;
end

%% DISPLAYING OUTPUTS

disp('OUTPUTS')
fprintf(' The number of terms used in the polynomial is %g',terms_needed)
fprintf('\n The smallest positive real root is %g',minval)
disp('  ')