All pastes #2133600 Raw Edit

Unnamed

public text v1 · immutable
#2133600 ·published 2012-03-29 18:23 UTC
rendered paste body
 
%TO VISUALISE THE FRACTAL

%The polynomial and it's derivative are decribed.
f=@(x) x^3 - 1;                 
g=@(x) 3*x^2;                 

%The density of the data is decided by altering the size of n (Wait roughly
%10 seconds for n=1000).
n=1000;
p=zeros(n,n); %Define an initial nxn free matrix for a value of p for each point on the 
z=p;          %graph that represent convergence.
l=p;
a=p; 
solmap=zeros(n,n,3); %Create three nxn matricies.
x=linspace(-2.2,2.2,n); %Create n evenly distributed points along the x axis between -2.2 and 2.2
y=x;  %Create n evenly distributed points along the y axis between -2.2 and 2.2
[X,Y] = meshgrid(x,y); %Create a xy cartesian coordinate system to plot on.

for k=1:n %Iterate between 1 and n
    for j=1:n %Iterate between 1 and n
        z(k,j)=x(j)+1i*y(n+1-k); %Create imaginary values to insert into original polynomials.
        m=(z(k,j)); %Convert these values into a new variable.
        while abs(f(m)) > 1e-6 %Contraint that causes the iteration to stop at a specific accuracy.
            m=m-((f(m)/g(m))); %The Newton Raphson formula.
            a(k,j)=imag(m); %Variable 'a' represents the imaginary part of 'm' which is later transfered to 'y' values.
           
            if a(k,j)<1e-6 && a(k,j)>-1e-6 %Contraint that causes the iteration to stop at a specific accuracy.
                a(k,j)=0; %As above
            end
            p(k,j)=p(k,j)+1; %Used to count the number of iterations.
        end
    end
end
maxp=max(max(p)); %Identifies the point with the largest number of iterations.
for k=1:n %Iterate between 1 and n
    for j=1:n %Iterate between 1 and n
        if a(k,j) == 0 %If all the elements are zero used the first solmap matrix generated.
            solmap(k,j,1)=1-(p(k,j)/maxp)^0.4;
        elseif a(k,j)> 0 %If the elements are greater than zero use the second matrix.
            solmap(k,j,2)=1-(p(k,j)/maxp)^0.4;
        else
            solmap(k,j,3)=1-(p(k,j)/maxp)^0.4; %Likewise with the third.
        end
    end
end

pcolor(X,Y,p); %Create a graph with a color bar, where the darkest spots are points of convergence
title('Newton Fractals for a thrid order polynomial')%i.e. where the roots are.
colorbar;
colormap hot;
shading interp; axis equal tight;

figure
imagesc(solmap);shading interp;axis equal tight off; %Create another graph that shows the same idea.