Home > K25104 > NonlinearSystems > NewtonFractal.m

NewtonFractal

PURPOSE ^

Calculates the basins of attractions when the Newton method for root

SYNOPSIS ^

function img = NewtonFractal(P)

DESCRIPTION ^

 Calculates the basins of attractions when the Newton method for root
   finding is used on a polynomial of degree at most 7.
 Input arguments:
   P, vector specifying the coefficents of the polynomial starting with 
     the coefficient of the highest power
 Output arguments:
   img, each pixel is coloured according to the root it converges to and 
     shaded by the number of iterations necessary, the more iterations, 
     the darker shade

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function img = NewtonFractal(P)
0002 % Calculates the basins of attractions when the Newton method for root
0003 %   finding is used on a polynomial of degree at most 7.
0004 % Input arguments:
0005 %   P, vector specifying the coefficents of the polynomial starting with
0006 %     the coefficient of the highest power
0007 % Output arguments:
0008 %   img, each pixel is coloured according to the root it converges to and
0009 %     shaded by the number of iterations necessary, the more iterations,
0010 %     the darker shade
0011 
0012 NITER = 100;        % maximum number of iterations
0013 threshold = .001;   % convergence criterion
0014 pixelnum = 1000;    % resolution of image
0015 colorArr = [7,3];   % seven colours specified by their RGB values
0016 %RED
0017 colorArr(1,1) = 1;colorArr(1,2) = 0;colorArr(1,3) = 0;
0018 %GREEN
0019 colorArr(2,1) = 0;colorArr(2,2) = 1;colorArr(2,3) = 0;
0020 %BLUE
0021 colorArr(3,1) = 0;colorArr(3,2) = 0;colorArr(3,3) = 1;
0022 %YELLOW
0023 colorArr(4,1) = 1;colorArr(4,2) = 1;colorArr(4,3) = 0;
0024 %WHITE
0025 colorArr(5,1) = 1;colorArr(5,2) = 0;colorArr(5,3) = 1;
0026 %CYAN
0027 colorArr(6,1) = 0;colorArr(5,2) = 1;colorArr(5,3) = 1;
0028 %RED
0029 colorArr(7,1) = 1;colorArr(6,2) = 1;colorArr(6,3) = 1;
0030 
0031 % generate grid over the square [-1, 1] x [-1, 1]
0032 [xs,ys] = meshgrid(linspace(-1,1,pixelnum), linspace(-1,1,pixelnum));
0033 % grid points interpreted as complex numbers, which is an array of
0034 % length pixelnum * pixelnum
0035 solutions = xs(:) + 1i*ys(:);
0036 % array of indices of grid points under consideration
0037 select = 1:numel(xs);
0038 % for each grid point initialize the necessary number of iterations
0039 % to the maximum
0040 niters = NITER*ones(numel(xs), 1);
0041 
0042 % calculate the roots of the polynomial
0043 Proots = roots(P);
0044 if isempty(Proots)
0045     disp('Polynomial has no roots');
0046     return;
0047 end
0048 % calculate the coefficients of the derivative
0049 Pderivative = zeros(length(P) - 1,1);
0050 for it = 1:length(P)-1
0051     Pderivative(it)=(length(P)-it)*P(it);
0052 end
0053 
0054 for iteration = 1:NITER
0055     % each iteration considers the entire grid minus the grid points
0056     % where convergence has occured
0057     oldi = solutions(select);
0058 
0059     % in newton's method we have z_{i+1} = z_i - p(z_i) / p'(z_i)
0060     solutions(select) = oldi - polyval(P,oldi) ...
0061         ./ polyval(Pderivative,oldi);
0062 
0063     % check for convergence or NaN (in case of a division by zero)
0064     differ = (oldi - solutions(select));
0065     % logical array marking converged grid points
0066     converged = abs(differ) < threshold;
0067     % logical array marking problematic grid points
0068     problematic = isnan(differ);
0069 
0070     % if converence occured update the necessary number of iterations
0071     niters(select(converged)) = iteration;
0072     % for problematic grid points set the number of iterations to the
0073     % maximum + 1
0074     niters(select(problematic)) = NITER+1;
0075     %remove indices of converged or problematic points
0076     select(converged | problematic) = [];
0077 end
0078 
0079 Max = max(niters);
0080 niters = reshape(niters,size(xs));
0081 solutions = reshape(solutions,size(xs));
0082 
0083 A = zeros(pixelnum, pixelnum);
0084 B = uint8(round(A * 255));
0085 RowCol = size(solutions);
0086 rows = RowCol(1);
0087 cols = RowCol(2);
0088 for i1 = 1:rows
0089     for i2 =1:cols
0090         % to which root did the method converge
0091         tmp = abs(repmat(solutions(i1,i2), size(Proots))-Proots);
0092         rootIndex = find(tmp<threshold);
0093         if ~isempty(rootIndex)
0094         % color associated with roots and rate of convergence
0095         B(i1,i2,1)=colorArr(rootIndex,1) * ...
0096             (1-(niters(i1,i2)/(Max))) * 255;
0097         B(i1,i2,2)=colorArr(rootIndex,2) * ...
0098             (1-(niters(i1,i2)/(Max))) * 255;
0099         B(i1,i2,3)=colorArr(rootIndex,3) * ...
0100             (1-(niters(i1,i2)/(Max))) * 255;
0101         end
0102     end
0103 end
0104 img = image(B);
0105 axis off;
0106 end

Generated on Mon 18-Jan-2016 10:25:49 by m2html © 2005