%
%   PROOTS.M   Demo for Working with Polynomials and Roots of Polynomials with Matlab
%
%   -->  Consider the following 5th order polynomial:
%            f(x) = x^5 - 3x^4 - 10x^3 + 10x^2 +44x + 48
%
%   This file works with this function in a number of ways.  f(x) is evaluated and
%   plotted for a range of x.  We find the roots/factors of the polynomial. This 
%   demo also shows how to multiply, divide, and find derivatives of polynomials.
%
%   File prepared by J. R. White, UMass-Lowell  (last update:  Nov. 2017)
%

      clear all,   close all
%
%   let's first define the desired polynomial by storing its coefficients in a vector
%   (note that the coefficients are in the order of decreasing powers of x)
      a = [1 -3 -10 10 44 48]
%
%   we can evaluate f(x) at any desired x value (say x = 3.5)
      disp('f(3.5) = '),  polyval(a,3.5)
%
%   we can also evaluate and plot it for a whole range of x values
      x = linspace(-3,5,101);   f = polyval(a,x);
      subplot(2,1,1),plot(x,f,'LineWidth',2),grid,axis([-3 5 -200 200]);
      title('PRoots:  Plot of f(x) = x^5 - 3x^4 - 10x^3 + 10x^2 + 44x + 48')
      xlabel('x value'),ylabel('f(x)')
%
%   Based on this figure, it looks like there are only 3 real roots.
%   However, we know that a 5th order polynomial has 5 roots.  Thus the other
%   two must be complex conjugates.  Let's use Matlab's roots command to see.
      r = roots(a)
%   Notice that the roots are stored in a column vector, and that there are 
%   two complex roots along with the three seen in the plot.
%
%   Given the factors of a polynomial in a vector, we can construct the 
%   actual polynomial coefficients with Matlab's poly command.
      p = poly(r)
%   Notice that the elements of p are identical to the original a vector.
%
%   Let's evaluate and plot the derivative of f(x)  (on the same plot as above)
      ap = polyder(a);  fp = polyval(ap,x);
      subplot(2,1,2),plot(x,fp,'LineWidth',2),grid,axis([-3 5 -200 200]);
      title('PRoots:  Plot of f''(x)')
      xlabel('x value'),ylabel('f''(x)')
%
%   We can also multiply and divide polynomials easily in Matlab.
%   To show this, let's 
%     1. extract the last two roots from the r vector -- call this rr
%     2. form a 2nd order polynomial from these roots  -- call this pp
%     3. divide the original p by pp  -- with Q and R as results (R should be zero)
%     4. multiply Q by pp to give back the original polynomial --  call this P
%     5. modify pp so that it is not an exact factor (here RR will not be zero)
      rr = r(4:5)                                      % step 1
      pp = poly(rr)                                    % step 2
      [Q,R] = deconv(p,pp)                             % step 3
      P = conv(Q,pp)                                   % step 4
      ppp = pp;  ppp(3) = 1,  [QQ,RR] = deconv(p,ppp)  % step 5
%
%   Well, there are many other things, but this is enough for now!!!
%