隱函數(shù)、方程求根、不動點和迭代
7.1知識要點與背景
? 7.1.1?? 隱函數(shù)存在定理與四連桿機構的運動
? 7.1.2? 不動點和函數(shù)迭代
7.2實驗 與觀察
7.2.1?? 隱函數(shù)的存在定理的可視化
1. 隱函數(shù)為什么存在？
【??? 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. 如何決定隱函數(shù)－非線性方程的求根
? zxy7-2.m
【?? global p???????????? %說明全局變量p，P用于本程序中和函數(shù)子程序zxy7_2f.m間傳遞參數(shù)
m=100;x=linspace(-5,5,m);????? %確定隱函數(shù)自變量的范圍
y0=-4.6;??????????????????????????????? %第一個方程的初值
y=[];f=[];
for k=1:m???? ???????%開始循環(huán)
?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???????????? ?????? %循環(huán)結束
plot(x(1:m),y(1:m),'r.-'),??? %繪制隱函數(shù)圖形
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;?? %計算函數(shù)???????????????? 】
7.2.2.用蛛網(wǎng)圖觀察不動點迭代
??????? 觀察程序： 下面的程序可以繪制這三個函數(shù)迭代的蛛網(wǎng)圖。
zxy7_3f.m
【%計算問題3中的三個函數(shù)，s=1、2、3時分別對應函數(shù)
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?? 畫一次迭代的蛛網(wǎng)圖, 改變p可調節(jié)箭頭的大小
?????? 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
【? %? 繪制迭代的蛛網(wǎng)圖(對問題3的三個函數(shù))
clf,clear
?s=2;??????????? %? 選擇函數(shù)：s=1、2、3時分別對應函數(shù)f1、f2、f3
if s==1???????? %對于函數(shù)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?????? %對于函數(shù)f1,決定坐標軸的范圍，將函數(shù)限制在同一范圍內(nèi)?? 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???????????????????????????????? % 畫蛛網(wǎng)圖，迭代過程為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)????????????? % 畫蛛網(wǎng)圖????
?? if i<=5
????? pause?????????????????????? % 按任意鍵逐次觀察前5次迭代的蛛網(wǎng)圖
?? end
?? x00=xt(2); y00=yt(2);?? % 將本次迭代的終點作為下一次迭代的起點。
?? z=[z,xt(1)];???????????????? %保存迭代點
end?????????????????????????????? 】
7.2.3 簡單和復雜：二次函數(shù)的迭代和混沌
???
? 觀察程序
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('參數(shù)r'),?? ylabel('迭代序列x')??????? 】
zxy7_5.m?? (用這個程序可以畫出蛛網(wǎng)迭代圖(圖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四連桿輸出角的運動規(guī)律和動畫模擬
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);??? %創(chuàng)建一個動畫的矩陣，保存20禎
for j=1:20;?????????????? %循環(huán)創(chuàng)建動畫數(shù)據(jù)
?? 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. 函數(shù)越平坦，迭代越快嗎？
2. 如何構造迭代函數(shù)使之具有較快的收斂速度？
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普適常數(shù)
2. 周期窗口
3. 混沌對初值的敏感性
4. 其它迭代函數(shù)
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]????????????? 】