Matlab: Newton-Fractal

[Spalding]

Hi FractalForum,

I am currently writing a code in Matlab to create Newton fractals. I am in school so, the code needn't to be perfect.

Code:

clc;
clear all;
close all;

funk = inline('x^3-1');                                                   
abl = inline('3*x^2');                                                      

n = 200;                                                                  
k = zeros(n,n);                                                             
l = zeros(n,n);                                                             
m = zeros(n,n);                                                             
x = linspace(-1,1,n);                                                       
y = linspace(-1,1,n);                                                       
[X,Y] = meshgrid(x,y);                                                      
tol = 0.0005;                                                               
maxIter = 40;                                                               

for u=1:n
   for v=1:n
        l(u,v)=x(v)+1i*y(n+1-u);                                            
        t=(l(u,v));                                                         
        while abs(funk(t)) > tol && k(u,v) < maxIter                        
            t=t-((funk(t))/(abl(t)));                                       
            m(u,v)=imag(t);                                                 
            if m(u,v)<tol && m(u,v)>-tol                                    
               m(u,v)=0;                                                    
            end
            k(u,v)=k(u,v)+1;                                               
        end
    end
end



pcolor(X,Y,m);
colormap jet;                                                               

                                                                          
shading interp;                                                             

I want to:

a) Choose as many colors as roots + black for divergence.
b) For each number in the visible plane, determine whether the method converges and if yes to which root. Color the point according to the result.
c) The brightness should display the number of iterations needed. (e.g. red for the root and light red => high number of iterations, dark red => low number of itertions

The code above:
pcolor(X,Y,m) solves question number b)
pcolor(X,Y,k) solves question number c)

But how can I combine both and how is it possible to include a) ?

Thank you

[superheal]

If you know the roots beforehand, lets say 3 roots then you will know where newton's converged simply be measuring the distance from each root. You have to pick the root with the smallest distance.
Knowing the roots will also help you pick the colors. Lets say root 1 is red, root 2 is green, root 3 is blue.

If you are using an arbitrary function with unknown roots, the best thing you can do is to create an array of different colors (many if possible) and depending of the final value of the newton's iteration pick a color from that array for the root. I usually combine the argument and the norm of the complex number into a single index pointing to that array. Obviously with this way you gonna have some times different roots having the same color.

Now for the shading part, after you pick the root color you can use the number of iterations with linear interpolation to get the final color.

Black         -------------- X Gray Shade--------------------->White
0 Iterations -------------- X Iteration    ---------------------> Max Iteration reached on the current image

after you get the shade, combine it with root color 50%-50%

final_red = shade_red *0.5 + root_red * 0.5
final_green = shade_green *0.5 + root_green * 0.5
final_blue = shade_blue *0.5 + root_blue * 0.5

[Spalding]

Thank you, I tried.
The following code defines 'red', 'blue' and 'green'
and draw each root in another color.
But how can i measure the distance to the root if I don't know the root?
e.g. f(x) = x^3-1  
I know that 3 roots exists, but I don't want to solve the problem beforehand.
So I wrote the following code:

Code:

myMap = [1, 0, 0
    0, 1, 0
    0, 0, 1];

pcolor(X,Y,m);
colormap(myMap);   

The matrix k contains the number of iterations depending of the start value.
How can I use these values to create a shading color ?
In the code above maxIter = 40 so, for 40 it should be black
while the color should be white when the number of iterations is 0.

The color white is given by the triple:
white = [1, 1, 1];

Is it like 100/maxIter*k(u,v)*[1, 1, 1] ?
So I can implement it in my code like this:

Code:

s = zeros(n,n);
s(u,v) = (1/maxIter)*k(u,v)

so my shading 'percentage' is [1,1,1]*s(u,v)

[superheal]

Well since you cant know the roots, the idea is to make a hashtable along with a hashfunction.
myMap is that table, so you need a way to convert the final (real, imaginary) number into an index number between 1 and 3 (or 0 and 2 w/e).
So my proposal is to keep a list of converged values.

list =[]

converged_value = iterate()

if ( list is empty) {
   add converged_value to list
   color_index = 0
}
else {
   if (converged_value+-error is in the list) {
       color_index = get index of found item
   }
   else {
      add converged_value to list
      color_index = list size - 1
   }

}

now for the shading part
white_red = 1
white_green = 1
white_blue = 1

black_red = 0
black_green = 0
black_blue = 0

shade_red = white_red + (current_iteration / max_iteration)*(black_red - white_red)
shade_green = white_green + (current_iteration / max_iteration)*(black_green - white_green)
shade_blue = white_blue + (current_iteration / max_iteration)*(black_blue - white_blue)


And finally

final_red = myMap[color_index][0] * 0.5 + shade_red * 0.5
final_green = myMap[color_index][1] * 0.5 + shade_green * 0.5
final_blue = myMap[color_index][2] * 0.5 + shade_blue * 0.5

I think matlab uses 1 offset indexes so change it accordingly

Last thoughts, if you want to handle more roots add more colors to myMap table and just to be sure use color_index % length(myMap) when you acess it so you wont run out of it. It will cyrcle back the same colors at that point though.

[Spalding]

Sorry, but i am confused.
myMap is the matrix which defines the colors. So the first line is 'red', 2nd is 'green', 3rd is 'blue'.
In my code, I built a matrix m with the final (real, imaginary) number.
and after finishing the whole for-loop, I draw the whole matrix.

Now I need to built a new 3x1 vector for each start value?
So I need to write it into the newton-algorithm and draw each point after it is approximated?

Really, I have no idea how to change it.

[superheal]

Confused in what way?

If the newton iteration converges to a number that must be a root right?
So we need to keep every root into a list. If we converge into the same root again and that root was already found, we will know its index on the list.
If the root was not found then we add it to the list.
Since we can never converge exactly to the same number due to floating point errors, when you check the list for the existence
of a root you need to check if the distance norm(converged_value - current_root) < epsilon.
So every number that will converge to that root will have the same color.
That was the idea.

[Spalding]

Well, I did it like this:

Code:

r = roots(p);
white_red = 1;
white_green = 1;
white_blue = 1;
black_red = 0;
black_green = 0;
black_blue = 0;

iter=1; 


for u=1:n
   for v=1:n
        l(u,v)=x(v)+1i*y(n+1-u);                                            
        t=(l(u,v));                                                         
        while abs(funk(t)) > tol && iter < maxIter                         
            t=t-((funk(t))/(abl(t)));                                       
            iter=iter+1;                                                        
            if abs(r(1)-funk(t)) > tol
                shade_red = white_red + (iter/maxIter)*(black_red - white_red);
                final_red = 0.5*[1, 0, 0] + 0.5*shade_red*[1, 0, 0];
            end
            if abs(r(2)-funk(t)) > tol
                shade_green = white_green + (iter/maxIter)*(black_green - white_green);
                final_green = 0.5*[0, 1, 0] + 0.5*shade_green*[0, 1, 0];
            end
            if abs(r(3)-funk(t)) > tol
                shade_blue = white_blue + (iter/maxIter)*(black_blue - white_blue);
                final_blue = 0.5*[0, 0, 1] + 0.5*shade_blue*[0, 0, 1];
            end
        end
    end
end

now I just need something to plot the pixel? right?

[superheal]

I think you confused the R,G,B color format.

Color red in RGB is: 1, 0, 0
Color blue in RGB is: 0, 0, 1
Color green in RGB is: 0, 1, 0
Color white in RGB is: 1, 1, 1
Color black in RGB is: 0, 0, 0

I am doing the calculations for every color channel R,G,B seperatelly.

Also:
this abs(r(2)-funk(t)) > tol

might need to change to abs(r(2) - t) < tol for every root checking.
and Instead of just checking those 3 roots maybe do it in a loop for the entire length of the vector r.
If one of them matches the criterion then you can break the loop.