terminate loop diverging to inf

I am trying to iterate a basins of attraction map, where the values diverge to either 0 or inf. When they diverge to inf the loop is continuous. I would like to terminate it at a particularly large value and store the index value as p(k,kk)=2 and then continue iterations, and similarly when it diverges to 0 to store as p(k,kk)=1

Any help is massively appreciated !

NT = 1024;
N = 1024;
r = 0.9;
omega = (sqrt(5)+1)/2;
phi = 2*pi*omega;
eps = 1e-6;
x0_vec = linspace(-1.0,1.8,N); y0_vec = linspace(-1.0,1.8,N);
z0_vec = complex(x0_vec,  y0_vec);
P = zeros(N,N);

for k = 1:length(y0_vec)
    y0 = y0_vec(k);

    for kk = 1:length(x0_vec)
        x0 = x0_vec(kk);
z = z0_vec;
  
for j = 1:NT

    z(j+1) = z(j)^2 + r*exp(1i*phi)*z(j);
    %stop loop continuing to inf

if abs(z(j+1)) <= eps
    disp('Winning attractor: at 0');
    P(k,kk) = 1;    
else
  
    disp('Winning attractor: at infinity');
    P(k,kk) = 2;
end
end
end
 end

[dpb]

MAGICNUMBEER=99999999999999999999;
....
%loop iterator here...

if newvalue>=MAGICNUMBER
  P(k,kk)=2
  break                             % breaks innermost loop only...
end

Answer 1

[Torsten]

I'm not sure, but I guess you meant something like this:

NT = 1024;
N = 2048;
r = 0.9;
omega = (sqrt(5)+1)/2;
phi = 2*pi*omega;
eps = 1e-6;
x0_vec = linspace(-1.0,1.8,N); y0_vec = linspace(-1.0,1.8,N);
P = zeros(N,N);
for i = 1:numel(x0_vec)
    for j = 1:numel(y0_vec)
      z0 = complex(x0_vec(i),y0_vec(j));
      for k = 1:NT
        z0=z0^2+r*exp(1i*phi)*z0;
        if abs(z0) <= eps
           %disp('Winning attractor: at 0');
           P(i,j)=1; 
           break
        end
        if abs(z0)  > 1e8
           %disp('Winning attractor: at infinity');
           P(i,j)=2; 
           break
        end
      end
    end
end
contourf(x0_vec,y0_vec,P)
colorbar

sum(P(:)==1)

ans = 1242217

sum(P(:)==2)

ans = 2952087

sum(P(:)==1)+sum(P(:)==2)

ans = 4194304

numel(x0_vec)*numel(y0_vec)

ans = 4194304

[Reece]

This works perfectly, thank you!

Answer 2

[Sulaymon Eshkabilov]

There is one critical point is the eps (MATLAB built in var ==> not to use) = 1e-6 is not met. The smallest abs value of z is 0.012; therefore, EPS_lim (different name instead of eps) = 0.0012; Why not ot use vectorization instead of the loop. The [for end] loop is very slow and just displays the phrases in disp() for over million times, e.g.:

NT = 1024;
N = 1024;
r = 0.9;
omega = (sqrt(5)+1)/2;
phi = 2*pi*omega;
EPS_lim = 0.012;
x0_vec = linspace(-1.0,1.8,N); 
y0_vec = linspace(-1.0,1.8,N);
z0_vec = complex(x0_vec,  y0_vec);
P = zeros(N,N);
z = z0_vec;
z(2:end+1) = z(1:end).^2 + r*exp(1i*phi)*z(1:end);
IDX1 = abs(z(2:end)) <= EPS_lim;
P(:,IDX1)=1;
IDX2 = P==0;
P(IDX2)=2;
disp('Number of Winning attractor: at 0');

Number of Winning attractor: at 0

sum((P(:)==1))

ans = 7168

disp('Winning attractor: at infinity');

Winning attractor: at infinity

sum((P(:)==2))

ans = 1041408

[Walter Roberson]

By the way, matlab eps is a function rather than a variable. You can pass a parameter to it, and the result is the epsilon relative to the input value.