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.