用Mathematica7绘制分形图Mandelbrot集合的程序问题，急！！

cx,cy的意思就是参数Z的实部和虚部的值。这个程序中令cx=1,cy=0,则迭代方程还是Zn+1 =Zn^k+C，n趋向正无穷时就有z=z^n+c了

这程序设c的模小于2，迭代次数小于50次

用了个Block函数将z作为局部分量和全局变量区分清楚

通过画密度函数画出c的图像

实际上，这个程序可简化成

M1[x_,y_]:=

Block[{z,k=0},

z=x+y*I

While[(Abs[z]<2.0)&&(k<50),++k

z=z^2+(x+y*I)]

Return[k]

]

M2[px_,py_,p_]:=

Block[{t},

t=DensityPlot[M1[xx,yy],{xx,-1.5,0.5},{yy,-1.2,1.2},p,

Mesh->False,ColorFunction->Hue]

Return[t]

]

Mandelbrot=M2[x,y,{PlotPoints->120,PlotLabel->"Mandelbrot集"}]

你看看行不行

曼德勃罗特集的计算与显示

曼德勃罗特集是易并行计算的一个典型例子.采用分治技术,并行算法设计时分为静态

任务分配和动态任务分配(可用work-pool or processor farm).

1.测试用例

其中底部的数据(real_min, imag_min) to (real_max, imag_max)表示复平面窗口,real_min表示

实部最小值, imag_min表示虚部最小值,real_max表示实部最大值, imag_max表示虚部最大

值.

(-0.84950, 0.21000) to (-0.84860, 0.21090) (0.32000, -0.45000) to (0.50000, 0.05000) (0.26304, 0.00233) to (0.26329, 0.00267)

(-0.63500, 0.68000) to (-0.62500, 0.69000) (-0.46510, -0.56500) to (-0.46470, -0.56460) (-1.50000, -1.00000) to (0.50000, 1.00000)

2.曼德勃罗特集的计算

显示曼德勃罗特集是处理位映射图像的一个例子.首先要对图像进行计算,且计算量很

大.曼德勃罗特集是复数平面中的点集,当对一个函数迭代计算时,这些点将处于准稳定状

态(quasi-stable),即将会增加或减少,但不会超过某一限度.通常该函数为:

zk+1 = zk

2 + c

式中zk+1是复数z = a + bi的第k+1次迭代,c是确定该点在复平面中位置的复数值.z的初始

值为0.迭代将一直进行下去,直到z的幅值(向量长度,这里为22ba+)大于2或者迭代次

数已经达到某种任意的规定的限度.

化简计算.

z2 = a2 - b2 + 2abi

用zreal表示z的实部,zimag表示z的虚部.则

zreal = a2 - b2 + creal

zimag = 2ab + cimag

3.顺序代码

structure complex {

float real

float image

}

对一点值的计算并返回迭代次数的例程:

int cal_pixel(complex c)

{ int count, max

complex z

float temp, lengthsq

max = 256

z.real = 0

z.imag = 0

count = 0

do {

temp = z.real * z.real - z.imag * z.imag + c.real

z.imag = 2 * z.real * z.imag + c.imag

z.real = temp

lengthsqc = z.real * z.real + z.imag * z.imag

count++

} while ((lengthsq <4.0) &&(count <max))//直到z的幅值大于2或者迭代次数达到max

return count

}

因此,所有给定的曼德勃罗特点必将处在以原点为中心,半径为2的圆中.

计算和显示这些点的代码需要对坐标系统进行一定的缩放来与显示区域的坐标系统相匹配.

假设显示高度为disp_height,宽度为disp_width.而点在显示区域中的位置为(x, y).如果显

示复数平面的这个窗口具有最小值(real_min, imag_min)和最大值(real_max, imag_max),则每

个点需用以下系数加以缩放.

c.real = real_min + x * (real_max - real_min) / disp_width

c.imag = imag_min + y * (imag_max - imag_min) / disp_height

设置scale_real = (real_max - real_min) / disp_width

scale_imag = (imag_max - imag_min) / disp_height

缩放代码:

for (x = 0x <disp_widthx++)

for (y = 0y <disp_heighty++) {

c.real = real_min + ((float) x * scale_real)

c.imag = imag_min + ((float) y * scale_imag)

color = cal_pixel(c)

display(x, y, color)

}

4.并行代码

1)静态任务分配

假定显示区域为640 × 480,在一个进程中要计算10行.即将10 × 640像素变为一组,共48

个进程.为代码如下:

主进程Master

for ( i = 0, row = 0i <48i++, row = row + 10)

send (&row, pi)

for (i = 0i <(480 * 640)i++)

recv(&c, &color, PANY)

display(c, color)

}

从进程Slave (process i)

recv(&row, Pmaster)

for (x = 0x <disp_widthx++)

for (y = rowy <(row + 10)y++) {

c.real = real_min + ((float) x * scale_real)

c.imag = imag_min + ((float) y * scale_imag)

color = cal_pixel(c)

send(&c, &color, Pmaster)

}

改进:成组发送数据.减少通信启动次数.先将结果保存在数组中,然后以一个消息将整个

数组发送给主进程.主进程可用一个通配符以任意顺序接收来自从进程的消息.

2)动态任务分配—工作池/处理器农庄(work pool / processor farm)

动态负载平衡可以用工作池方法实现.在我们的问题中,像素集(更确切应该是坐标集)构成

了任务.任务数是固定的,要计算的像素数在计算前是已知的.各个处理器从工作池中请求

下一个像素子集的坐标.

主进程

count = 0

row = 0

for (k = 0k<procnok++) {

send(&row, pk, data_tag)

count++

row++

}

do {

recv(&slave, &r, color, PANY, result_tag)

count--

if (row 0)

从进程

recv(y, Pmaster, ANYTAG, source_tag)//接收Pmaster发送的第y行的点

while(source_tag == data_tag) { //判断是否还有消息需要处理

c.imag = imag_min + ((float) y * scale_imag)

for (x = 0x <disp_widthx++) {

c.real = real_min + ((float) x * scale_real)

color = cal_pixel(c)

}

send(&i, &y, color, Pmaster, result_tag)//将所计算的第y行点的color发给Pmaster

recv(y, &r, Pmaster, source_tag)//接收下一任务

} //如果退出while循环, 则一定是source_tag = termiate_tag, 表明没有任务,程序终止

---