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 3So 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 13The result of the synthetic division shows that
x^2 + 3*x + 3 = (x + 5)(x-2) + 13and 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) + 3If 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.0000iObserve the quadratic convergence in both cases.