求迭代分形的MATLAB程序

function

cycle6(px,py,r,n)

%px,py和r分别是大圆圆心坐标和半径,n是迭代次数

%例子:

%cycle6(0,0,1,4)axis

equal

t=linspace(-pi,pi)

plot(r*cos(t)+px,r*sin(t)+py)

theta=pi/3*(0:5)

newr=r/3

newpx=px+2*newr*cos(theta)

newpy=py+2*newr*sin(theta)

if

n>1

for

i=1:6

hold

oncycle6(newpx(i),newpy(i),newr,n-1)

end

end

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集"}]

你看看行不行

---