MCS 471 Lecture Seven

Roots of Polynomials

We defined "synthetic division" which is useful in Newton's method when applied to polynomials.

MATLAB calls the synthetic division "deconvolution". To divide the polynomial x^3 + x^2 - 3*x - 3 by x - 2 we execute the following commands:

>> p = [1 1 -3 -3]; d = [1 -2];
>> [q,r] = deconv(p,d)

q =

     1     3     3


r =

     0     0     0     3

So we find that

  x^3 + x^2 - 3*x - 3 = (x^2 + 3*x + 3)*(x - 2) + 3,

where r contains p(2) = 3. To compute the derivative of p at 2 we go one step further:

>> [q1,r1] = deconv(q,d)

q1 =

     1     5


r1 =

     0     0    13

The result of the synthetic division shows that

  x^2 + 3*x + 3 = (x + 5)(x-2) + 13

and if we relate this to the original polynomial p, we can write

  x^3 + x^2 - 3*x - 3 = (x + 5)*(x - 2)^2 + 13*(x-2) + 3

If we compute the derivative of the last expresssion and evaluate at 2, we find that p'(2) = 13, which confirms the 13 we see in r1.

Below is the code for a simple MATLAB program to perform Newton's method on a polynomial:

function [q,x] = newton(p,x,eps,N)
%
%  Does one step with Newton's method using synthetic division
%  on the polynomial p, starting at x.
%  On return is the quotient and a new approximation for the root.
%  The iteration stops until the accuracy eps is reached.
%  No more than N steps are performed.
%
%  Examples :
%     p = [1 0 -1]; x = 0.9;
%     [q,x] = newton(p,x,eps,4)
%     p = [1 0 1]; x = 0.1 + 0.9*i;
%     [q,x] = newton(p,x,eps,4)
%
n = size(p,2);
d = [1 -x];
[q,r] = deconv(p,d);
fprintf('  real(x)               imag(x)                 dx        f(x) \n');
for i = 1:N
   [q1,r1] = deconv(q,d);
   dx = r(n)/r1(n-1);
   x = x - dx;
   d(2) = -x;
   [q,r] = deconv(p,d);
   fprintf('%.15e %.15e  %.2e  %.2e\n', real(x), imag(x), dx, r(n));
   if (abs(dx) < eps) | (abs(r(n)) < eps)
      fprintf('succeeded after %d steps\n', i);
      return;
   end
end
fprintf('failed requirements after %d steps\n', N);

The output for the two examples:

>> p = [1 0 1]; x = 0.9 + 0.9*i;
>> [q,x] = newton(p,x,eps,4)
  real(x)               imag(x)                 dx        f(x) 
1.005555555555556e+00 0.000000000000000e+00  -1.06e-01  1.11e-02
1.000015346838551e+00 0.000000000000000e+00  5.54e-03  3.07e-05
1.000000000117761e+00 0.000000000000000e+00  1.53e-05  2.36e-10
1.000000000000000e+00 0.000000000000000e+00  1.18e-10  0.00e+00
succeeded after 4 steps

q =

     1     1


x =

     1
>> p = [1 0 1]; x = 0.1 + 0.9*i; 
>> [q,x] = newton(p,x,eps,4)    
  real(x)               imag(x)                 dx        f(x) 
-1.097560975609757e-02 9.987804878048780e-01  1.11e-01  2.56e-03
1.274517695672162e-05 9.999402989079396e-01  -1.10e-02  1.19e-04
-7.609680927669747e-10 1.000000001700982e+00  1.27e-05  -3.40e-09
-1.294393212231021e-18 1.000000000000000e+00  -7.61e-10  0.00e+00
succeeded after 4 steps

q =

   1.0000            -0.0000 + 1.0000i


x =

  -0.0000 + 1.0000i

Observe the quadratic convergence in both cases.