《数学实验与Matlab》程序2

实验7.隐函数、方程求根、不动点和迭代7.1知识要点与背景7.1.1 隐函数存在定理与四连杆机构的运动7.1.2 不动点和函数迭代7.2实验与观察7.2.1 隐函数的存在定理的可视化1. 隐函数为什么存在？【 clf,plot(Y(:,40),z1(:,40),'.-');hold on,xlabel('x'),ylabel('y'),plot([-5,5],[0,0],'r-'),legend('z=F(x0,y)','z=0'); 】观察程序zxy7_1.m【 clear,clfa=-5;b=5;c=-5;d=5;n=60;u=[15 25 40];x=linspace(a,b,n);y=linspace(c,d,n);[X,Y]=meshgrid(x,y);z1=Y.^3./(2+0.2*sin(X.*Y))+X.^2-4*X; z2=zeros(size(X));surf(X,Y,z1),shading flat,hold on,mesh(X,Y,z2),hidden off,xlabel('x','fontsize',16);ylabel('y','fontsize',16);zlabel('z','fontsize',16);r0=(abs(z1-z2)<=0.7); %处理交线zz=r0.*z1; yy=r0.*Y; xx=r0.*X;plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'r.','markersize',12);plot3(X(:,u),Y(:,u),z1(:,u),'b.-','markersize',10); 】2. 如何决定隐函数－非线性方程的求根zxy7-2.m【 global p %说明全局变量p，P用于本程序中和函数子程序zxy7_2f.m间传递参数m=100;x=linspace(-5,5,m); %确定隐函数自变量的范围y0=-4.6; %第一个方程的初值y=[];f=[];for k=1:m%开始循环 p=x(k);%第k 个方程的自变量x(k)通过全局变量p 传递到zxy7_2f.m 中 [y1,fval,exitflag,output] = fzero('zxy7_2f',y0); %求第k 个方程y0=y1; %将第k 方程的解作为下一个方程的初值 y=[y,y1];f=[f,fval]; %保存计算结果end %循环结束plot(x(1:m),y(1:m),'r.-'), %绘制隐函数图形axis([-5 5 -5 3]),grid on 】zxy7_2f.m【 function z=zxy7_2f(y)global p %在这里也要对应说明全局变量p ，使得可以获得外界传递来的P 值 x=p;z=y.^3/(2+0.2*sin(x *y))+x^2-4*x; %计算函数 】7.2.2.用蛛网图观察不动点迭代观察程序： 下面的程序可以绘制这三个函数迭代的蛛网图。

zxy7_3f.m【%计算问题3中的三个函数，s=1、2、3时分别对应函数321,,f f ffunction y=zxy7_3f(x,s)if s==1y=(4-x.*x);elseif s==2y=4./(1+x);elseif s==3y=x-(x.^2+x-4)./(2*x+1);end 】zxyplot7_3.m【 %zxyplot7_3 画一次迭代的蛛网图, 改变p 可调节箭头的大小function out = zxyplot7_3(x,y,p)% 已知有向折线段的始点为(x(1),y(1)),终点为(x(2),y(2)),画出有向折线段% (x(1),y(2))――>(x(2),y(2))% |% |u(1)=0; v(1)=(y(2)-y(1)); u(2)=eps; v(2)=eps;h=quiver([x(1) x(1)],[y(1) y(2)],u,v,p);set(h,'color','red');hold on,u(1)=(x(2)-x(1)); v(1)=0; u(2)=eps; v(2)=eps;h=quiver([x(1) x(2)],[y(2) y(2)],u,v,p); set(h,'color','red');plot([x(1) x(1) x(2)],[y(1) y(2) y(2)],'r.-') 】主程序zxy7_3.m【% 绘制迭代的蛛网图(对问题3的三个函数)clf,clears=2; % 选择函数：s=1、2、3时分别对应函数f1、f2、f3if s==1 %对于函数f1,决定坐标轴的范围。

(以便得到较好的图形显示效果) a=-5;b=5;x00=2.1;y00=0; %初始点x=linspace(a,b,80); y0=x; y1=zxy7_3f(x,s);c=-(1+0.3)*max(abs(y1));d=(1+0.3)*max(abs(y1));elseif s==2|s==3 %对于函数f1,决定坐标轴的范围，将函数限制在同一范围内a=0;b=5; %显示，以便进行观察和比较x00=4;y00=0; %初始点x=linspace(a,b,80);y0=x; %计算直线y=xy1=zxy7_3f(x,s); %计算曲线y=f(x)c=0;d=max((1+0.3)*max(abs(y1)),5.2);endclear y;y=[y0;y1];plot(x,y,'linewidth',2),legend('y=x','y=f2'),hold on, % 画图plot([a b],[0,0],'k-',[0 0],[c d],'k-'),axis([a,b,c,d]), % 画坐标轴z=[];for i=1:10 % 画蛛网图，迭代过程为n=10次xt(1)=x00; yt(1)=y00; %决定始点坐标xt(2)=zxy7_3f(xt(1),s); yt(2)=zxy7_3f(xt(1),s); %决定终点坐标 zxyplot7_3(xt,yt,0.6) % 画蛛网图if i<=5pause % 按任意键逐次观察前5次迭代的蛛网图endx00=xt(2); y00=yt(2); % 将本次迭代的终点作为下一次迭代的起点。

z=[z,xt(1)]; %保存迭代点7.2.3 简单和复杂：二次函数的迭代和混沌观察程序zxy7_4.m【 clear,clf,n0=100;n=150;x0=0.1;hold on,s=[];xx=[];for r=2.6:0.001:4%for r=0:0.3:4x(1)=x0; %初值for j=2:nx(j)=r*x(j-1)*(1-x(j-1)); %迭代ends=[r*ones(1,n+1-n0);s]; %在固定的r处画出n以后的迭代值xn,xn+1，…xx=[x(n0:n);xx]; xs=max(x(n0:n));text(r-0.1,xs+0.05,['\it{r}=',num2str(r)]) ％标注endplot(s,xx,'k.','markersize',15); %调整点的大小获得较好的视觉效果%plot(s,xx,'k.','markersize',1); xlabel('参数r'), ylabel('迭代序列x') 】zxy7_5.m (用这个程序可以画出蛛网迭代图(图7.10))【 clfr=2.7; %r=3.2;r=3.9; n=15;x00=0.2;y00=0;a=0;b=1;x=linspace(a,b,50);plot(x,x),axis([a b a b]),hold on,theaxes=axis;y=r*x.*(1-x);plot(x,y),z=[];for i=1:nxt(1)=x00; yt(1)=y00;xt(2)=r*xt(1).*(1-xt(1)); yt(2)=xt(2);zxyplot7_4(xt,yt,0.6)x00=xt(2); y00=yt(2);z=[z,xt(2)];7.3.应用、思考与练习7.3.1四连杆输出角的运动规律和动画模拟1. 确定四杆机构的转角关系2. 动画模拟四杆机构的运动zxy7_6.m【 x=-8:0.5:8; %定义曲面[XX,YY]=meshgrid(x);r=sqrt(XX.^2+YY.^2)+eps;Z=sin(r)./r;surf(Z); %画出祯theAxes=axis; %保存坐标值，使得所有帧都在同一坐标系中fmat=moviein(20); %创建一个动画的矩阵，保存20祯for j=1:20; %循环创建动画数据surf(sin(2*pi*j/20)*Z,Z) %画出每一步的曲面axis(theAxes) %使用相同的坐标系fmat(:,j)=getframe; %拷贝祯到矩阵fmat中endmovie(fmat,10) %演示动画10次%这很有趣】7.3.2轨道飞行器的拦截7.3.3怎样保证或加速迭代序列的收敛1. 函数越平坦，迭代越快吗？2. 如何构造迭代函数使之具有较快的收敛速度？3. 关于迭代的收敛性和收敛速度的定理zxy7_7.m【clf,x=linspace(a,b,50);y=x; plot(x,y),axis([a b a b]),hold on,theaxes=axis;theaxes=axis;button=1;x1=[];y1=[];while button==1[xi,yi,button]=ginput(1); h=plot(xi,yi,'+'),axis(theaxes);x1=[xi,x1];y1=[yi,y1];end[p,c]=polyfit(x1,y1,4);ypolyfit=polyval(p,x);plot(x,ypolyfit),axis(theaxes);x00=(b-a)/2;y00=0;for i=1:100xt(1)=x00; yt(1)=y00; xt(2)=polyval(p,xt(1)); yt(2)=xt(2);zxyplot7_4(xt,yt,0.6); x00=xt(2); y00=yt(2); z=[z,xt(2)];end 】7.3.4.混沌有哪些特点？1.Feigenbaum普适常数2. 周期窗口3. 混沌对初值的敏感性4. 其它迭代函数7.4 非线性方程组求解zxy7_8f.m【 function f=zxy7_7f(x)f=[sin(x(1))+x(2)+x(3)^2*exp(x(1))-4;x(1)+x(2)*x(3); x(1)*x(2)*x(3)]; 】【 options=optimset('Display','off'); %若取'Display'为iter则显示中间迭代结果 [x,fval]=fsolve('zxy7_8f',[1 1 1],options) 】【 syms x1 x2 x3 a c %用syms 对每个符号变量进行说明 f1='a*sin(x1)+c*x2+x3^2*exp(x1)-a'; %象这样定义各个方程f2='x1+x2*x3'; f3='x1*x2*x3';[x1,x2,x3]=solve(f1,f2,f3); %用solve指令求解solution=[x1,x2,x3] 】实验8.线性代数实验8.1 实验(Ⅰ)：用Matlab学线性代数8.1.1实验与观察：向量组的线性关系和解线性方程组1. 用线性组合的方式产生向量组【 clear n=3; m=2; a=-10; b=10;rand('seed',32), A = unifrnd(a,b,[n,m]) 】【 x = unifrnd(-1,1,[1,m]), A(:,3)=x(1)*A(:,1)+x(2)*A(:,2) 】2.Gauss消元法和向量组的线性关系的判定【 B=rref(A) 】【 r=2;m=2;s=5;X=[-B(1:2,r+1:m+s) ;eye(m+s-r)] %基础解系】【 r1=rank(A(:,1:2)), r2=rank(A(:,1:3)) 】【 B=rref(A(:,1:3)) 】【 B=rref(A);C=B(1:2,:) 】3 . 观察程序zxy8_1.m （主程序，产生向量并画图）【 clf,n=3;m=2;s=5; % 确定相关的参数a=-10;b=10;rand('seed',32),A = unifrnd(a,b,[n,m]), %产生m个n维向量（生成向量）r=[1:m]; l=1;p=1.2;zxy8_1plot(A,r,l,p) %将向量画出for i=m+1:m+s %产生组合向量x = unifrnd(-1,1,[1,m]);A(:,i)=zeros(n,1);for j=1:mA(:,i)=A(:,i)+x(j)*A(:,j);endendhold on,r=[m:m+s];l=2;p=0.5;zxy8_1plot(A,r,l,p) %用另一种方式画组合向量】其中调用了画图子程序zxy8_1plot【function out=zxy8_1plot(A,r,l,p)%A为若干3维向量拼成的矩阵, 绘制这些向量；r为向量，指明要绘制A的列向量指标。