0001 function img = NewtonFractal(P)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 NITER = 100;
0013 threshold = .001;
0014 pixelnum = 1000;
0015 colorArr = [7,3];
0016
0017 colorArr(1,1) = 1;colorArr(1,2) = 0;colorArr(1,3) = 0;
0018
0019 colorArr(2,1) = 0;colorArr(2,2) = 1;colorArr(2,3) = 0;
0020
0021 colorArr(3,1) = 0;colorArr(3,2) = 0;colorArr(3,3) = 1;
0022
0023 colorArr(4,1) = 1;colorArr(4,2) = 1;colorArr(4,3) = 0;
0024
0025 colorArr(5,1) = 1;colorArr(5,2) = 0;colorArr(5,3) = 1;
0026
0027 colorArr(6,1) = 0;colorArr(5,2) = 1;colorArr(5,3) = 1;
0028
0029 colorArr(7,1) = 1;colorArr(6,2) = 1;colorArr(6,3) = 1;
0030
0031
0032 [xs,ys] = meshgrid(linspace(-1,1,pixelnum), linspace(-1,1,pixelnum));
0033
0034
0035 solutions = xs(:) + 1i*ys(:);
0036
0037 select = 1:numel(xs);
0038
0039
0040 niters = NITER*ones(numel(xs), 1);
0041
0042
0043 Proots = roots(P);
0044 if isempty(Proots)
0045 disp('Polynomial has no roots');
0046 return;
0047 end
0048
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
0056
0057 oldi = solutions(select);
0058
0059
0060 solutions(select) = oldi - polyval(P,oldi) ...
0061 ./ polyval(Pderivative,oldi);
0062
0063
0064 differ = (oldi - solutions(select));
0065
0066 converged = abs(differ) < threshold;
0067
0068 problematic = isnan(differ);
0069
0070
0071 niters(select(converged)) = iteration;
0072
0073
0074 niters(select(problematic)) = NITER+1;
0075
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
0091 tmp = abs(repmat(solutions(i1,i2), size(Proots))-Proots);
0092 rootIndex = find(tmp<threshold);
0093 if ~isempty(rootIndex)
0094
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