Matlab中是否有更簡便的方法生成這樣的矩陣M?
下面的代碼用來畫Mandelbrot set
其中生成矩陣M這一步感覺有些麻煩,還得需要id1, id2。而Mathematica中可以用這樣比較簡潔的寫法:
M = Table[mandel[r + i*I], {i, -1, 1, 0.005}, {r, -2, 0.5, 0.005}]
(* mandel是自定義的函數 *)Matlab中有沒有更方便的做法呢?
function M=mandelTest()
d=0.005;
M = zeros(length(-1:d:1), length(-2:d:0.5));
id1 = 1;
for i = -1:d:1
id2 = 1;
for r = -2:d:0.5
M(id1, id2) = mandel(complex(r,i));
id2 = id2 + 1;
end
id1 = id1 + 1;
end
endfunction n = mandel(z)
n = 0;
c = z;
for n=0:99
if abs(z)&>2
break
end
z = z^2+c;
end
end% 畫圖語句 :imshow(mandelTest(), [])
% mandel 是自定義的函數
M = arrayfun(@mandel, (-1:.005:1)+1i*(-2:.005:.5));
謝邀。
Mathematica裡面mandel是你自定義的函數吧。利用Table能在一行內實現兩重循環,這是Mathematica比較兇殘的地方。
Matlab裡面,我覺得如果要簡化代碼,則不應該用循環,而是用向量化運算來解決問題。大體給個思路:
首先,for n = 0:100 這個循環肯定是需要的。於是,為了實現向量化運算,要把mandelTest裡面的內容搬到mandel裡面。也就是隻用一個函數來實現你的功能。定義一個矩陣Z0,包含所有你要考察的點。那麼每次迭代為,Z = Z.^2+Z0。另外定義一個矩陣N,存放迭代次數n,可以這樣更新,N((~M).*(abs(Z)&>2)) = n。這裡M是上一步裡面abs(Z)&>2的結果。這樣循環完,整個矩陣就出來了。這個思路有兩個問題,一個是效率比較低,每次循環要計算所有的點,另一個是如果點密度比較大,則內存可能不夠。對於前者,可以考慮在每次循環只計算還沒發散的點,這可通過M來篩選。對於後者,可以把待求點陣分解為小區域,分別求解。手機不方便寫代碼,先給個我的建議吧,有空再來詳細補充一下。
題主可以用空間來換代碼簡潔度,試著看看meshgrid這個函數。
另外,mathematica作為函數式編程的一員,其代碼簡潔度本身就不是matlab能比擬的雖然收到了邀請,這個我也不大會... 不敢瞎寫,不是一個方向....
發揮 matlab 矩陣計算的優勢
通過 meshgrid 和 arrayfun 兩個函數的應用來「消滅」 for 循環,從而實現簡化。d=0.005
xspace = -1:d:1;yspace = -1:d:1;[x, y] = meshgrid ( xspace, yspace);M = arrayfun ( @(z) Interaction(z), x+y*sqrt(-1) );
Interaction 需要自己編寫,就是對複數 x+y*sqrt(-1) 反覆迭代,判斷是否落在圓形區域內,類似樓主程序中的 mandel 函數。
Mandelbrot集,上帝的戒指-function Mandelbrot(iter,pixel)
% MANDELBROT(ITER,PIXEL) draws the Mandelbrot set;
% ITER gives the number of iterations and
% PIXEL indicates the number of points on the x-axis.
% MANDELBROT uses iter = 23, pixel = 400.
% This file was generated by students as a partial fulfillment
% for the requirements of the course "Fractals", Winter term
% 2004/2005, Stuttgart University.
% Author : Jochen Schurr
% Date : Nov 2004
% Version: 1.0
% default setting
switch nargin
case 0
iter=23;
pixel=400;
end
% decompose x- and y-axis according to the chosen
% proportion of r = 3:4.
r = 3/4;
x = linspace(-2.5,1.5,pixel);
y = linspace(-1.5,1.5,round(pixel*r));
% create the matrix C containing all pixel
% to be analysed.
[Re,Im] = meshgrid(x,y);
C = Re + i * Im;
% The matrix B assigns a number to every pixel which is
% characteristic with respect to the velocity of divergence.
% Therefore B and C have the same dimension.
B = zeros(round(pixel*r),pixel);
Cn = B; % C0 = 0+0i
for l = 1:iter
Cn = Cn.*Cn + C; % performing the map
B = B + (abs(Cn)&<2); % if |cn| &> 2 the sequence diverges
end;
% plot settings
imagesc(B);
colormap(jet);
axis equal
axis off
這是我以前寫的代碼
%% define the region
x = -2:0.005:0.5;
y = -1:0.005:1;
%% broadcast
[x,y] = meshgrid(x, y);
sz = size(x);
%% initialize the iterates
z0 = x + 1i*y;
z = zeros(sz);
c = zeros(sz);
j = true(sz);
%% Mandelbrot iteration
depth = 100;
for k = 1:depth
z(j) = z(j).^2 + z0(j);
j = abs(z) &< 2;
c(j) = k;
end
%% show
image(c)
Elapsed time is 1.394212 seconds.
@lx其實可以藉助bsxfun的廣播功能實現meshgrid,避免產生中間變數
參照上述幾個答案,改了一下,能夠動態顯示整個生成過程:
function MandelbrotPlot(maxDepth, N)
% MANDELBROTPLOT show the generation process of Madelbrot set
x = linspace(-2.5, 1.5, N);
y = linspace(-1.5, 1.5, N);
z = bsxfun(@plus, x, y*i); % use the Binary Singleton Expansion func
% to generate "x+y*i" instead of using MESHGRID
c0 = z;
n_z = zeros(size(z));
index = zeros([size(z), maxDepth]);
for k = 1:maxDepth
n_z = n_z + (abs(z)&<2);
z = z.^2 + c0;
index(:,:,k) = n_z; % index records every iteration state of n_z
end
for n = 1:maxDepth
imagesc(index(:, :, n))
colormap hsv
%axis equal
axis off
pause(0.15)
end
end
&>&> MandelbrotPlot(100, 1000)
推薦閱讀: