terminate loop diverging to inf

조회 수: 1 (최근 30일)
Reece
Reece 2023년 11월 24일
편집: Torsten 2023년 11월 24일
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
  댓글 수: 1
dpb
dpb 2023년 11월 24일
MAGICNUMBEER=99999999999999999999;
....
%loop iterator here...
if newvalue>=MAGICNUMBER
P(k,kk)=2
break % breaks innermost loop only...
end

댓글을 달려면 로그인하십시오.

채택된 답변

Torsten
Torsten 2023년 11월 24일
편집: Torsten 2023년 11월 24일
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
  댓글 수: 1
Reece
Reece 2023년 11월 24일
This works perfectly, thank you!

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

Sulaymon Eshkabilov
Sulaymon Eshkabilov 2023년 11월 24일
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
  댓글 수: 1
Walter Roberson
Walter Roberson 2023년 11월 24일
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.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Environment and Settings에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by