下面的代碼用來畫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
end

function 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

請注意，我只是搬運，並沒有版權。如果有任何版權問題我不負法律責任。===================================================尋找過程中我發現了幾段特別牛的代碼，c++直接畫點陣圖，貼兩張上來看看

===========================================知乎，最近你是怎麼了，發圖各種發不上來？

---

這是我以前寫的代碼

%% 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)

---

推薦閱讀：

相关文章