terminate loop diverging to inf

2 ビュー (過去 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 Center および File ExchangeEnvironment and Settings についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by