How to vectorize this code?

Can anyone please help me to write a vectorized version of following code assignment. I tried by using meshgrid instead of two outer loops of x and y values:

=============================

clc
clear all
close all
r1=[0;0];
r2=[1;1];
r3=[-1;-1];
for x=-2:0.03:2
     for y=-2:0.03:2
         X=[x;y];
for j=1:30
     f=[X(1)^3 - X(2), X(2)^3 - X(1)]';
    Jf=[3*X(1)^2, -1; -1, 3*X(2)^2];
    X=X-Jf\f;
end
format long
    X;
    if norm(X-r1)<1e-8 % if the distance between X and r1 is less than 10^-8
     c='b';        % use the color red
elseif norm(X-r2)<1e-8
     c='m';        % use the color blue
elseif norm(X-r3)<1e-8
     c='r';        % use the color green
else               % if not close to any of the roots
     c='y';        % use the color black
    end
plot(x,y,'.','Color',c,'Markersize',40);
hold on
    end
     end

=============================================

Thanx

Answer 1

[Roger Stafford]

Your outer two for-loops can be vectorized, but the innermost one which involves the iteration process cannot - at least as far as I am aware.

 [X,Y] = meshgrid(-2:0.03:2);
 for j = 1:30
   D = 9*X.^2.*Y.^2-1; % Avoid matrix division
   N = 6*X.^2.*Y.^3+2*X.^3; % Use direct expression for newton-raphson formula
   X = (6*X.^3.*Y.^2+2*Y.^3)./D;
   Y = N./D;
   tr = (X-0).^2+(Y-0).^2<1e-16; % Color the dots red for true tr
   tb = (X-1).^2+(Y-1).^2<1e-16; % Color the dots blue for true tb
   tg = (X+1).^2+(Y+1).^2<1e-16; % Color the dots green for true tg
   ty = ~(tr|tb|tg);  % Color the other dots black
 end
 plot(X(tr),Y(tr),'r.',etc....

By the way, there are six more roots to your equations. They are all complex numbers on the unit circle in the complex plane. For example, (x,y) = (1i,-1i) is a root.

[Athar]

thank you Roger it is revealing answer to me. But I am unable to complete the plot(X(tr),Y(tr),'r.',etc.... part. What I get is a single red dot in the centre. Please complete this. What about using scatter(), would it speed up the fig making?

[Roger Stafford]

Don't try to plot the black points. They are far, far out because the newton-raphson algorithm is not converging for them but rather diverging. That is true for both the vectorized method and your nested for-loop method. It only works when the initial points are within some reasonable neighborhood of the solution. To attempt a plot of all the diverging points will cause the scaling of the plot to be so huge as to make all valid solution points look as if they are packed right at the origin. Look at the values in your plot axes to see this. The last few lines should be:

   tr = (X-0).^2+(Y-0).^2<1e-16; % Color the dots red for true tr
   tb = (X-1).^2+(Y-1).^2<1e-16; % Color the dots blue for true tb
   tg = (X+1).^2+(Y+1).^2<1e-16; % Color the dots green for true tg
 end
 plot(X(tr),Y(tr),'ro',X(tb),Y(tb),'bo',X(tg),Y(tg),'go')

[Roger Stafford]

I would think it would be far more informative to plot the original positions, rather than final positions, of the points which converged to one of your three solutions. You already know where the final points are to be found but you don't know where they came from.

You can do this by saving a copy of the original X and Y and using the 't' logical arrays with them. This will make plots of the "neighborhoods" I mentioned previously.

 [X,Y] = meshgrid(-2:0.03:2);
 X0 = X; Y0 = Y;
 .....
 plot(X0(tr),Y0(tr),'r*',X0(tb),Y0(tb),'b*',X0(tg),Y0(tg),'g*')

Answer 2

[Walter Roberson]

Look at scatter()