隐函数、方程求根、不动点和迭代

隐函数、方程求根、不动点和迭代
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,clf
a=-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时分别对应函数
function y=zxy7_3f(x,s)
if s==1
   y=(4-x.*x);
elseif s==2
   y=4./(1+x);
elseif s==3
   y=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))    
%              |
%              |
%    (x(1),y(1))                         
  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,clear
 s=2;            %  选择函数：s=1、2、3时分别对应函数f1、f2、f3
if 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=x
   y1=zxy7_3f(x,s);    %计算曲线y=f(x)
   c=0;d=max((1+0.3)*max(abs(y1)),5.2);
end
clear 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<=5
      pause                       % 按任意键逐次观察前5次迭代的蛛网图
   end
   x00=xt(2); y00=yt(2);   % 将本次迭代的终点作为下一次迭代的起点。
   z=[z,xt(1)];                 %保存迭代点
end                               】
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:4
   x(1)=x0;                                %初值
   for j=2:n
      x(j)=r*x(j-1)*(1-x(j-1));                    %迭代
   end  
   s=[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)])    ％标注
end
plot(s,xx,'k.','markersize',15);        %调整点的大小获得较好的视觉效果
%plot(s,xx,'k.','markersize',1);   xlabel('参数r'),   ylabel('迭代序列x')        】
zxy7_5.m   (用这个程序可以画出蛛网迭代图(图7.10))
【   clf
r=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:n         
   xt(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)];
end                               】
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中
end           
movie(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:100
   xt(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]              】