Introduction
In Trust-Region-Dogleg and Newton-Raphson - A Quick Comparison on a 5th order polynomial, we saw how the Newton-Raphson method produces a fractal pattern in the complex plane when finding roots of a polynomial. The fractal patterns show how Newton is unstable near the boundary of the basins of attraction for the roots. In this post, we will look at taking a step that is stable at every step according to the fixed point theorem. I am anticipating that this will produce a much smoother pattern in the complex plane.
polynomial setup
Given the polynomial
\[P(z) = z^5 + z^2 - z + 1 = 0\]
find a root starting from anywhere in the interval of -2 to 2 and -2i to 2i.
Common code:
clc;clear; close all;
% setup polynomial function
coefficients = [1 0 0 1 -1 1];
polynomialFunction = @(z) polyval(coefficients, z);
% setup the derivative of the polynomial function
coefficients_dot = polyder(coefficients);
polynomialDerivativeFunction = @(z) polyval(coefficients_dot, z);
% roots of the polynomial. roots uses a eigenvalues of a companion matrix approach.
% rootsOfPolynomial = roots(coefficients);
% exact roots
rootsOfPolynomial(1) = -1i;
rootsOfPolynomial(2) = 1i;
rootsOfPolynomial(3) = (2^(2/3)*(3^(1/3) - 3^(5/6)*1i)*(69^(1/2) + 9)^(1/3))/12 + (2^(2/3)*3^(1/3)*(9 - 69^(1/2))^(1/3))/12 + (2^(2/3)*3^(5/6)*(9 - 69^(1/2))^(1/3)*1i)/12;
rootsOfPolynomial(4) = (2^(2/3)*(3^(1/3) + 3^(5/6)*1i)*(69^(1/2) + 9)^(1/3))/12 + (2^(2/3)*3^(1/3)*(9 - 69^(1/2))^(1/3))/12 - (2^(2/3)*3^(5/6)*(9 - 69^(1/2))^(1/3)*1i)/12;
rootsOfPolynomial(5) = -(2^(2/3)*3^(1/3)*((69^(1/2) + 9)^(1/3) + (9 - 69^(1/2))^(1/3)))/6;
% Calculate a grid of starting points using the ndgrid format. x (real part)
% follow rows (1st dimension) and y follows the columns (2nd dimension).
% ndgrid is format is transpose of meshgrid format. meshgrid format is used
% for plotting surfaces.
nStartingPointReal = 400;
nStartingPointImaginary = 401;
startingPointReal = linspace(-2,2,nStartingPointReal)';
startingPointImaginary = linspace(-2,2+eps(2),nStartingPointImaginary)*1i;
startingPoint = startingPointReal + startingPointImaginary; % using singleton expansion
maxIterations = 100;
functionTolerance = 1e-14;
divideByZeroTolerance = 1000; % 1000*eps(p)Newton-Raphson Iteration
Using Newton’s method starting anywhere in the interval -2 to 2 and -2i to 2i produces the following pattern.
% setup for iteration
isConverged = false(nStartingPointReal, nStartingPointImaginary);
isDivideByZero = false(nStartingPointReal, nStartingPointImaginary);
whileFlag = true;
iteration = 0;
z = startingPoint;
% check if the starting value was already converged.
% polynomial is less then functionTolerance.
p = polynomialFunction(z(~isConverged));
isConverged(~isConverged) = abs(p) < functionTolerance;
% begin iteration
p = polynomialFunction(z(~isConverged));
while whileFlag
% increment iteration number
iteration = iteration + 1;
p = polynomialFunction(z(~isConverged));
% Calculate the derivative of the polynomial.
% Check for divide by 0. If p_dot is close to machine precision
% of p, then its a divide by 0 error.
% Stop all further iteration.
p_dot = polynomialDerivativeFunction(z(~isConverged));
isDivideByZeroIteration = abs(p_dot) < divideByZeroTolerance*eps(abs(p));
isDivideByZero(~isConverged) = isDivideByZeroIteration;
isConverged(~isConverged) = isDivideByZeroIteration;
% calculate Newton iteration on unconverged points.
z(~isConverged) = z(~isConverged) - p(~isDivideByZeroIteration)./p_dot(~isDivideByZeroIteration);
% check for convergence
p = polynomialFunction(z(~isConverged));
isConverged(~isConverged) = abs(p) < functionTolerance;
% set while flag
% all roots have converged or reached max iterations
whileFlag = iteration <= maxIterations || all(all(~isConverged));
end
% successful roots
isRoot = isConverged & ~isDivideByZero;
% which root did newton's method converge too?
% calculate the difference between the converged root and the roots of the polynomial
% the smallest distance corresponds to the root index
zDiff = abs(z - reshape(rootsOfPolynomial,1,1,[])); % using singelton expansion
[rootErrors,rootIndex] = min(zDiff,[],3); % minimum along 3rd dimension
rootIndex(~isRoot) = 0;
% create the fractal image
% imagesc uses the meshgrid format where
% x is the 2nd dimension
% y is the first dimension
% therefore, transpose rootIndex because its in the ndgrid format
figure;
imagesc(startingPointReal', imag(startingPointImaginary), rootIndex')
c = prism(numel(rootsOfPolynomial)+1);
colormap(c);
xlabel("Real Axis"); ylabel("Imaginary Axis")
ylabel(colorbar('Ticks',1:numel(rootsOfPolynomial)),"Root Index")
hold all;
plot(rootsOfPolynomial,'ko','MarkerSize',10,'LineWidth',3)
title("Newton-Raphson Convergence for z^5+z^2-z+1")
Newton-Raphson Iteration with step sized derived from the fixed point theorem
A fixed point is
\[x_{k+1} = g(x_k)\]
In the context of newton’s method, applying the fixed point theorem yields the following:
\[x_{k+1} = x_k - \alpha \frac{P(x_k)}{P'(x_k)}\]
where \(\alpha\) is a step size that is normally 1. The convergence of the fixed point iteration is guaranteed if the following condition is satisfied:
\[|g'(x_k)| < 1\]
The derivative of \(g\) is given by:
\[g'(x_k) = 1 - \alpha \left( \frac{P'(x_k)P'(x_k) - P(x_k)P''(x_k)}{P'(x_k)^2} \right)\]
To ensure convergence, we can choose \(\alpha\) such that \(|g'(x_k)| < 1\). This can be achieved by setting \(\alpha\) to a value that satisfies the inequality. We need the following to be true:
\[-1 < 1 - \alpha \left( \frac{P'(x_k)P'(x_k) - P(x_k)P''(x_k)}{P'(x_k)^2} \right) < 1 \]
Algorithmically, this is what will lead to our choice of \(\alpha\):
- Compute \(g'(x_k)\) using the formula above.
- If \(|g'(x_k)| < 1\), then set \(\alpha = 1\) (standard Newton step).
- If \(|g'(x_k)| \ge 1\), then choose \(\alpha\) to satisfy the inequality by the following:
- If \(g'(x_k) > 1\), then set \(\alpha = \frac{2}{\frac{P'(x_k)P'(x_k) - P(x_k)P''(x_k)}{P'(x_k)^2}}\) to ensure \(g'(x_k) < 1\).
- If \(g'(x_k) < -1\), then set \(\alpha = \frac{2}{\frac{P'(x_k)P'(x_k) - P(x_k)P''(x_k)}{P'(x_k)^2}}\) to ensure \(g'(x_k) > -1\).