Home > K25104 > NonlinearSystems > HalleyFractal.m

HalleyFractal

PURPOSE ^

Calculates the basins of attractions when Halley's method

SYNOPSIS ^

function img = HalleyFractal(P)

DESCRIPTION ^

 Calculates the basins of attractions when Halley's 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 = HalleyFractal(P)
0002 % Calculates the basins of attractions when Halley's method
0003 %   for root finding is used on a polynomial of degree at most 7.
0004 % Input arguments:
0005 %   P, vector specifying the coefficents of the polynomial starting
0006 %     with the coefficient of the highest power
0007 % Output arguments:
0008 %   img, each pixel is coloured according to the root it converges to
0009 %     and shaded by the number of iterations necessary, the
0010 %     more iterations, the darker shade
0011 
0012 NITER = 100;        % maximum number of iterations
0013 threshold = .001;   % convergence criterion
0014 pixelnum = 800;     % 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 %  lengthpixelnum * 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 
0049 % calculate the coefficients of the derivative
0050 Pderivative = zeros(length(P) - 1,1);
0051 for it = 1:length(P)-1
0052     Pderivative(it)=(length(P)-it)*P(it);
0053 end
0054 % calculate the coefficients of the second derivative
0055 P2derivative = zeros(length(Pderivative) - 1,1);
0056 for it = 1:length(Pderivative)-1
0057     P2derivative(it)=(length(Pderivative)-it)*Pderivative(it);
0058 end
0059 
0060 for iteration = 1:NITER
0061     % each iteration considers the entire grid minus the grid pints
0062     % where convergence has occured
0063     oldi = solutions(select);
0064 
0065     solutions(select) = oldi - 2 .* polyval(P,oldi) .* ...
0066         polyval(Pderivative,oldi) ./ (2.*polyval(Pderivative,oldi).^2 ...
0067         - polyval(P,oldi).*polyval(P2derivative,oldi));
0068 
0069     %check for convergence or NaN
0070     differ = (oldi - solutions(select));
0071     %logical array marking converged grid points
0072     converged = abs(differ) < threshold;
0073     %logical array marking problematic grid points
0074     problematic = isnan(differ);
0075 
0076     % if converence occured update the necessary number of iterations
0077     niters(select(converged)) = iteration;
0078     % for problematic grid points set the number of iterations to the
0079     % maximum + 1
0080     niters(select(problematic)) = NITER+1;
0081     %remove indices of converged or problematic points
0082     select(converged | problematic) = [];
0083 end
0084 
0085 Max = max(niters);
0086 niters = reshape(niters,size(xs));
0087 solutions = reshape(solutions,size(xs));
0088 
0089 A = zeros(pixelnum, pixelnum);
0090 B = uint8(round(A * 255));
0091 RowCol = size(solutions);
0092 rows = RowCol(1);
0093 cols = RowCol(2);
0094 for i1 = 1:rows
0095     for i2 =1:cols
0096         % to which root did the method converge
0097         tmp = abs(repmat(solutions(i1,i2), size(Proots))-Proots);
0098         rootIndex = find(tmp<threshold);
0099         if ~isempty(rootIndex)
0100         %color associated with roots and rate of convergence
0101         B(i1,i2,1)=colorArr(rootIndex,1) * ...
0102             (1-(niters(i1,i2)/(Max))) * 255;
0103         B(i1,i2,2)=colorArr(rootIndex,2) * ...
0104             (1-(niters(i1,i2)/(Max))) * 255;
0105         B(i1,i2,3)=colorArr(rootIndex,3) * ...
0106             (1-(niters(i1,i2)/(Max))) * 255;
0107         end
0108     end
0109 end
0110 img = image(B);
0111 axis off;
0112 end

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