ppt02迭代分形图形

u=[uuu]; subplot(3,3,n+1)，plot(u(:,1),u(:,2))，axis([-0.5,0.5,0,1]) end toc 计算出的各点实数坐标如下，为了对比，后面的复数结果也列出。 u = x y u = x + iy
0 0 %o点 0 0 0.33333 %a点 0 + 0.33333i 0 0.33333 %a点 0 + 0.33333i -0.16667 0.62201 %d点 -0.16667 + 0.62201i 0 0.33333 %a点 0 + 0.33333i 0 0.66667 %b点 0 + 0.66667i 0 0.66667 %b点 0 + 0.66667i 0.16667 0.95534 %e点 0.16667 + 0.95534i 0 0.66667 %b点 0 + 0.66667i 0 1 %c点 0 + 1i 3.1.2 树程序2 以复数表示点的坐标,作图方法相同,计算过程略快，程序可读性更强。 u=[0,i]; subplot(3,3,1); plot(u)， axis([-0.5,0.5,0,1]) for n=1:7 uuu=[]; for I=0:length(u)/2-1 diff=(u(2*I+2)-u(2*I+1))/3; p1=u(2*I+1)+diff; p2=p1+diff; lp=p1+diff*(sqrt(3)/2+1/2*i); rp=p2+diff*(sqrt(3)/2-1/2*i); uu=[u(2*I+1);p1;p1;lp;p1; p2;p2;rp;p2;u(2*I+2)]; uuu=[uuu;uu]; end
end 3.2.2 雪花程序2 所有的边上新点坐标同时计算,减少一重for循环 figure new=[0,i,sqrt(3)/2+i/2,0]; subplot(3,3,1); plot(new), axis equal; for k=1:8; old=new; n=length(old)-1; diff = (old(2:end) - old(1:end-1)) / 3; new(1:4:4*n-3) = old(1:end-1); new(2:4:4*n-2) = old(1:end-1) + diff; new(3:4:4*n-1) = new(2:4:4*n-2) + diff*(0.5+0.5*sqrt(3)*i); new(4:4:4*n) = old(1:end-1)+ diff + diff; new(4*n+1)=old(n+1); subplot(3,3,k+1); plot(new), axis equal; end 3.2.3 雪花程序3 相似移动，将图形缩小，移动，旋转组合成一条边，再将这条边旋转， 移动，形成另外两边。 a=[0,i]; b=a*(-0.5 - 0.86603i)+i; c=(a-i)*(-0.5 + 0.86603i); subplot(3,3,1), plot([a,b,c]); axis equal for k=1:8 subplot(3,3,k+1) a1=a/3; a=[a1, a1*(0.5 + 0.86603i)+i/3,... (a1-i/3)*(0.5 - 0.86603i)+i*2/3, a1+2*i/3]; b=a*(-0.5 - 0.86603i)+i;
。
1.2 例2 Logistic模型 xn+1 = µ(xn − x2 n ) (0 < µ < 4, 0 < x < 1) 给定x初值,对不同的µ值循环计算x，共循 环250次,循环计算150次后开始画图。 注解： 计算所得的x值是矩阵，行标对应循环次数，列标对 应µ值。使用矢量化编程以后，对所有的µ值同时计算一次新 的x值，得到矩阵x中的一列元素。下面两个程序不同之处在于 是否保留所有的x值。预先设置一个矩阵存放x的元素可以节约 扩充矩阵元素的时间。在Workspace中可以查看这两个程序中x的 区别，设置间断点以后运行程序，可以看到程序中各个语句的 运行结果。
3.分形图形画法
编程不仅仅是语句的运用......
1 1 1
3.1 分形树 将一条线段三等分，在等分点 上各画一条长度为原线段长度 的三分之一的线段，该线段与 原线段成固定的夹角。依次做 下去，得到如下图形。 3.1.1 树程序1(运行时间120秒) 用二列的实数矩阵(x,y)表示点 的坐标,按类似一笔画的方式 作图, 。参看右边电影演示。
uu=[m, i/3+m*(sqrt(3)*0.5+0.5i),m+i/3, ... 2i/3+m*(sqrt(3)*0.5-0.5i), m+2i/3]; subplot(3,3,k); plot(uu) u=uu; end
3.2 科契(Koch)雪花曲线 将等边三角形的一条边三等 分，将中间的一段去掉，代 之以一个更小的等边三角形的 两条边，依次做下去，将得到 图1.5 所示的雪花图。 3.2.1 雪花程序1 用复数表示点的坐标，对每边计算4个新点坐标。书《理论...模拟》p41 new=[0,i,sqrt(3)/2+i/2,0]; subplot(2,4,1); plot(new) axis equal for kkk=1:7; old=new; n=length(old)-1; for jj=0:n-1; diff=(old(jj+2)-old(jj+1))/3; new(4*jj+1)=old(jj+1); new(4*jj+2)=old(jj+1)+diff; new(4*jj+3)=new(4*jj+2)+diff*((1+sqrt(3)*i)/2); new(4*jj+4)=old(jj+1)+2*diff; end new(4*n+1)=old(n+1); subplot(2,4,kkk+1), plot(new), axis equal;
方法：用定性分析预测，借助算法、模型、参数，特例验算。 算得快：计算速度快。 方法：避免不必要或重复的计算，采用节省时间的算法。 写得好：简洁易懂，可读性强。 方法：程序集中在一个文件。加注解，模块化编程。排版规范， 标点正确。 不信！！！ 请看对比,画出同样图形的两个程序
1. 矢量化编程
“Life is too short to spend writing DO loops...” “生命诚可贵，少做DO循环......”
%程序1 u=2.6:0.001:4; x=0.6; for j=1:150, x=u.*(x-x.^2); end for i=1:100 x=u.*(x-x.^2); plot(u,x,’r.’) hold on; end 程序运行时间 0.28秒。不保留旧 的X 值，而是直接用它画图。
%程序2 u=2.6:0.001:4; X=ones(250,1401); X(1,:)=0.6*X(1,:); for j=1:250 X(j+1,:)=u.*(X(j,:)-X(j,:).^2); end plot(u,X(150:end,:),’r.’) 运行时间0.15秒，比程序1快。 保 留 所 有 X 值 ， 每 次 计 算 的X 值 生成矩阵的一行元素，最后用矩 阵X 的后150行作图，程序可读性 强。
u=uuu; subplot(3,3,n+1); end
plot(u)， axis([-0.5,0.5,0,1])
3.1.3 树程序3 以矩阵元素赋值的方法计算新点坐标，减少一重循环。 new=[0,i]; subplot(3,3,1); plot(new); axis([-0.5,0.5,0,1 ]); for k=1:8 old=new; n=length(old)/2-1; diff=(old(2:2:end)-old(1:2:end-1))/3; p1=old(1:2:end-1)+diff; p2=p1+diff; lp=p1+diff*(sqrt(3)/2+1/2*i); rp=p2+diff*(sqrt(3)/2-1/2*i); new(2:10:10*n+2)=p1; %第一次循环计算第2,12,22,...点 new(3:10:10*n+3)=p1; %第一次循环计算第3,13...点 new(4:10:10*n+4)=lp; new(5:10:10*n+5)=p1; new(6:10:10*n+6)=p2; new(7:10:10*n+7)=p2; new(8:10:10*n+8)=rp; new(9:10:10*n+9)=p2; new(10:10:10*(n+1))=old(2:2:end); %10,20,... new(1:10:10*n+1)=old(1:2:end-1); %1,11,21.... subplot(3,3,k+1)， plot(new) axis([-0.5,0.5,0,1 ]); end 3.1.4 树程序4(运行时间0.21秒) 每个图形缩小为原来的1/3,再移动位置, 添加到原图形。参看右边的电影。 u=[0,i]; subplot(3,3,1); plot(u) for k=2:8 m=u/3;
0.5
0.5
0.5
0 −0.5 1
0
0.5
0 −0.5 1
0
0.5
0 −0.5 1
0
0.5
0.5
0.5
0.5
0 −0.5 1
0
0.5
0 −0.5 1
0
0.5
0 −0.5 1
0
0.5
0.5
0.5
0.5
0 −0.5
0
0.5
0 −0.5
0
0.5
0 −0.5
0
0.5
tic theta=pi/6; u =[0,0;0,1]; subplot(3,3,1)， plot(u(:,1),u(:,2)), axis([-0.5,0.5,0,1]) for n=1:7 uuu=[]; for I=0:(length(u)/2-1) p1=(u(2*I+1,:)*2+u(2*I+2,:))/3 ; p2=(u(2*I+1,:)+u(2*I+2,:)*2)/3; lp=[cos(theta),-sin(theta);sin(theta),cos(theta)]*... [(u(2*I+2,1)-u(2*I+1,1));u(2*I+2,2)-u(2*I+1,2)]/3 ; lp=p1+lp’; rp=[cos(theta),sin(theta);-sin(theta),cos(theta)]*... [(u(2*I+2,1)-u(2*I+1,1));u(2*I+2,2)-u(2*I+1,2)]/3; rp=p2+rp’; uu=[u(2*I+1,:);p1;p1;lp;p1;p2;p2;rp;p2;u(2*I+2,:)]; uuu=[uuu;uu]; end
迭代是编程中常用的技巧，许多分形图形可以通过迭代来画，本节通过介 绍分形图形的画法来学习迭代的技巧，理解编程的要求，学习程序分析工 具Profiler的用法。
Βιβλιοθήκη Baidu
第二章
迭 代 –分 形 图形
编程能力很重要:
是基本功，要勤学苦练，熟能生巧，要耐心细致，一丝不苟......
好程序应该这样: 算得对：语法正确，计算结果正确。
2. 递归调用
例1. 归调用就是程序调用自己，以此实现循环计算的功能。尽 管计算n阶阶乘可以运用指令prod(1:n), 为了举例,下面编制一个示意 性的递归程序来计算阶乘。在后面将这种方法画更复杂的图形。 %program fact.m function output=fact(n) if n==1 output=1; else output=n*fact(n-1); end
%中止递归的条件 %递归调用的语句
调用时,输入fact(n),例如计算5的阶乘就是fact(5)。 例2. 用递归法画Sierpinski triangle。 调用时在指令窗口输入三角形 三 个 顶 点 和 迭 代 次 数 。N1=7; A1=-1; A2=sqrt(3)*i; A3=1; sjx(A1,A2,A3,N1)。 function sjx(A1,A2,A3,N1) triangle1=[A1,A2,A3,A1]; B1=(A1+A2)/2; B2=(A1+A3)/2; B3=(A2+A3)/2; triangle2=[B1,B2,B3,B1]; N2=N1-1; if N2>0; sjx(A1,B1,B2,N2); sjx(B1,A2,B3,N2); sjx(B2,B3,A3,N2); subplot(2,3,N2) fill([real(triangle1)],[imag(triangle1)],’g’); hold on fill([real(triangle2)],[imag(triangle2)],’w’) end
1.3 运用程序分析工具profil 运用profile可以查看各部分程序的运行时间，便于有针对性地改进程序。 从两个程序的分析报告可以看到，程序1重复作图100次，所用时间与程 序2全部运行时间相当。所以程序1 比程序2更费时间。 >>profile on >>bug >>profile report >>profile off %启动Profile %运行被检测的程序 %生成检测报告 %关闭profile
迭代并不是一定要通过for或while循环来完成。矢量化编程可以替代for循 环，发挥MATLAB进行矩阵运算的特点。 1.1 例1 画图形sin nπx (n = 1, 2, 3, 4; x = 0 : 0.001 : 1) x=[0:0.001:1]; y=sin(pi*[1:4]’*x); plot(x,y) 说明：x是矢量，函数值y是矩阵, plot用y的每行对X作图(新用法)