实验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,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时分别对应函数321,,f f f
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)) % | % |
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)]; %保存迭代点
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)];
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] 】
实验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:m
A(:,i)=A(:,i)+x(j)*A(:,j);
end
end
hold 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的列向量指标。
%有两种不同的绘制方式，由参数l和p控制。
%l=1,向量用红色粗线条绘制；l=2,向量用蓝色细线条绘制。p是控制箭头的大小参数。
a=-10;b=10;k=50;
x=linspace(a,b,k);y=x;
[X,Y]=meshgrid(x,y);axis([a b a b a b]);
xlabel('x1','fontsize',14),ylabel('x2','fontsize',14),zlabel('x3','fontsize',1 4)
if l==1 %第一种绘图方式
for i=r
h1=plot3([0 A(1,i)],[0 A(2,i)],[0 A(3,i)],'rs','linewidth',3); %画线段
u(1)=A(1,i);v(1)=A(2,i);w(1)=A(3,i);u(2)=eps; v(2)=eps; w(2)=eps;
h=quiver3([0 A(1,i)],[0 A(2,i)],[0,A(3,i)],u,v,w,p); %画箭头
set(h,'linewidth',1,'color','red'), axis([a b a b a b]),hold on,
text(A(1,i),A(2,i),A(3,i),['\leftarrowA',int2str(i)],'fontsize',14); %文字标注
end
Z=zxyplane(X,Y,A(:,1)',A(:,2)',0,0,0); %计算平面的z坐标
mesh(X,Y,Z), axis([a b a b a b]),hidden off, %画平面
elseif l==2 %第二种绘图方式，结构同上
for i=r
h1=plot3([0 A(1,i)],[0 A(2,i)],[0 A(3,i)],'o--');
set(h1,'linewidth',2,'color','blue')
u(1)=A(1,i);v(1)=A(2,i);w(1)=A(3,i);u(2)=eps;v(2)=eps; w(2)=eps;
h=quiver3([0 A(1,i)],[0 A(2,i)],[0,A(3,i)],u,v,w,p);
set(h,'linewidth',1,'color','blue'), axis([a b a b a b]),hold on,
text(A(1,i),A(2,i),A(3,i),['\leftarrowA',int2str(i)],'fontsize',14);
end
end 】
zxyplane.m （计算平面的z坐标）
【 function z=zxyplane(x,y,a,b,x0,y0,z0)
%向量a,b叉积构成平面的法矢t,平面过（x0,y0,z0）点
t(1)=det([a(2) a(3);b(2) b(3)]);
t(2)=-det([a(1) a(3);b(1) b(3)]);
t(3)=det([a(1) a(2);b(1) b(2)]);
if t(3)~=0
z=z0-(t(1)/t(3))*(x-x0)-(t(2)/t(3))*(y-y0);
k=find(z<=-8&z>=8);z(k)=NaN;
elseif t(3)==0
disp('t(3)=0,this programme can not finish the work you want to do ') end 】
8.1.2应用、思考与练习
1. 观察极大线性无关组的意义
【 B=rref(A), C=B(1:2,:), A(:,1:2)*C, A 】
2. 平面四连杆机构的设计
3. 用Matlab 做线性代数题（矩阵的符号演算）
◆ 试题 设T T T T a a a a a )4,,10,3(,)1 ,11,0(,)3,1,7,2(,)2,0,4,1(4321=-===，求 a 的值使得 1. a 4不能由a 1、a 2、a 3线性表示。 2. a 4可以由a 1、a 2、a 3线性表示，并写出表达式。 【 syms a
a1=[1;4;0;2]; a2=[2;7;1;3];a3=[0;1;-1;1];a4=[3;10;a;4]; A=[a1,a2,a3,a4]
for i=2:4 %行初等变换 A(i,:)=A(i,:)-A(1,:)*A(i,1); end
A(2,:)=A(2,:)/A(2,2); for i=3:4
A(i,:)=A(i,:)-A(2,:)*A(i,2); end A 】
8.2 实验(Ⅱ)：矩阵的相似化简
8.2.1 实验与观察：矩阵的特征―相似标准形的作用
1. 逼近直线的迭代点列
2. 估计直线－特征值、特征向量
【 [P,D] = eig(A) 】
【 x=linspace(a,b,30); [pc,lamda]=eig(A),pc=-pc; z1=pc(2,1)/pc(1,1)*x;
plot(x,z1,'linewidth',3) 】
3.特征值和特征向量决定迭代性质?
【 x0=[1 1]'; A=[1/5,99/100;1,0];
[P,D]=eig(A);
y0=inv(P)*x0;y=y0;
for i=1:50
y=[D^i*y0,y];
end
x=P*y;
plot(y(1,:),y(2,:),'o',x(1,:),x(2,:),'*'),legend('Y','X=PY') 】
4. 观测程序说明
zxy8_2.m的源代码如下：
【 clear，clf
a=-20*100;b=-a;c=a;d=b;p=0.1; %设定画图范围
n=100;
A=[1/5,99/100;1,0]; %设定矩阵
axis([a b c d]),grid on,hold on
button=1
while button==1
[xi,yi,button]=ginput(1); %用鼠标选初始点
plot(xi,yi,'o'),hold on,
X0=[xi;yi];X=X0;
for i=1:n
X=[A*X,X0]; %用这种方式作迭代，并画图
h=plot(X(1,1),X(2,1),'R.',X(1,1:2),X(2,1:2),'r:');hold on
quiver([X(1,2),1]',[X(2,2),1]',[X(1,1)-X(1,2),0]',[X(2,1)-
X(2,2),0]',p)
set(h,'MarkerSize',6), grid,
end
end
pause
x=linspace(a,b,30); %画最大特征值所对应的特征向量所决定的直线[pc,lamda]=eig(A),pc=-pc;
z1=pc(2,1)/pc(1,1)*x;
z2=pc(2,2)/pc(1,2)*x;
h=plot(x,z1),set(h,'linewidth',2) 】
8.2.2 应用、思考与练习
1. 植物基因的分布、杂交育种问题
3.高维线性离散动力系统
参考程序zxy8_3.m . (计算图8.6和8.7)
【 clear,clf
a=-2;b=-a;c=a;d=b;p=0.0001*abs(a); %设定画图范围
n=30;
A=[2,0;0,1/2] %设定矩阵
axis([a b c d]),hold on
button=1
while button==1
[xi,yi,button]=ginput(1); %用鼠标选初始点
plot(xi,yi,'o'),hold on,
X0=[xi;yi];X=X0;
for i=1:n
X=[A*X,X0]; %用这种方式作迭代，并画图
h=plot(X(1,1),X(2,1),'R.',X(1,1:2),X(2,1:2),'r:');hold on
quiver([X(1,2),1]',[X(2,2),1]',[X(1,1)-X(1,2),0]',[X(2,1)-
X(2,2),0]',p)
set(h,'MarkerSize',6), grid,
end
end 】
3. 主成分分析和线性变换
下面是相应的程序：
【x1=[-0.931 -3.931 20.069 -6.931 -4.931 20.931 -3.931 -12.931 28.069 0.069 -16.931 16.069 -11.931 -8.931 -7.931 0.069 0.069 -16.931 -8.931 -3.931 14.069 11.069 -15.931 21.069 -10.934 7.069 15.069 8.069 16.069];
x2=[-0.379 -0.379 12.621 -4.379 -3.379 -11.379 -1.379 -3.379 6.621 2.621 -9.379 3.621 -7.379 -2.379 -3.379 3.621 0.379 -6.379 -3.379 -0.379 2.621 7.621 -7.379 4.621 -3.379 4.621 9.621 4.621 5.621];
X=[x1',x2'];C=cov(X) ％计算协方差矩阵，注意数据需按列向量存入X 中。
[P,latent,explained] = pcacov(C) ％主成分分析，P即为所需的变换矩阵。】
【 clf, a=-20;b=30;c=-20;d=30;
t=linspace(a,b,20); z1= P(2,1)/P(1,1)*t; z2=P(2,2)/P(1,2)*t;hold on,
plot(x1,x2,'o',[a b],[0 0],':',[0 0],[c d],'b:',t,z1,t,z2),
xlabel('x1'),ylabel('x2') , axis([a b c d]),axis('equal') 】
zxy8_5.m
【 clear, clf,a=-20;b=30,n=8，x=linspace(a,b,n);y=x;
[x1,x2]=meshgrid(x);[x4,x3]=meshgrid(x); %生成网格线
A=[0 ,2.2;1 0.2]; a=pi/6;
for k=1:n %作网格线的线性变换
for l=1:n
z=A*[x1(k,l),x2(k,l)]'; %垂直线的的变换
y1(k,l)=z(1);y2(k,l)=z(2); %将变换后的坐标适当排列以便作图
z=A*[x3(k,l),x4(k,l)]'; %水平线的的变换
y3(k,l)=z(1); y4(k,l)=z(2);
end
end
t=linspace(a,b,30);z1=5*cos(t);z2=15+5*sin(t); %产生一个椭圆
z3=A*[z1;z2]; %对椭圆的变换
figure(1) %在第一个图形窗口画x1-x2平面的内容
plot(x1,x2,x3,x4,z1,z2),axis('equal')
figure(2) %在第二个图形窗口画y1-y2平面的内容
实验9 概率统计实验
9.1 实验（I）：Galton钉板试验
9.1.1 实验与观察: Galton钉板模型和二项分布
1. 动画模拟Calton钉板试验
【 rand('seed',1), u=rand(1,6) 】
【 rand('seed',2),u=rand(1,6) 】
观察程序zxy9_1.m
【 clear，clf, m=100;n=5;y0=2; %设置参数。
ballnum=zeros(1,n+1);
p=0.5;q=1-p;
for i=n+1:-1:1 %创建钉子的坐标x,y 。
x(i,1)=0.5*(n-i+1);y(i,1)=(n-i+1)+y0;
for j=2:i
x(i,j)=x(i,1)+(j-1)*1;y(i,j)=y(i,1);
end
end
mm=moviein(m); % 动画开始，模拟小球下落路径。
for i=1:m
s=rand(1,n); %产生n个随机数。
xi=x(1,1);yi=y(1,1);k=1;l=1; % 小球遇到第一个钉子。
for j=1:n
plot(x(1:n,:),y(1:n,:),'o',x(n+1,:),y(n+1,:),'.-'), % 画钉子的位置。
axis([-2 n+2 0 y0+n+1]),hold on
k=k+1; % 小球下落一格。
if s(j)>p
l=l+0; %小球左移。
else
l=l+1; %小球右移。
end
xt=x(k,l);yt=y(k,l); %小球下落点的坐标。
h=plot([xi,xt],[yi,yt]);axis([-2 n+2 0 y0+n+1]) %画小球运动轨
迹。
xi=xt;yi=yt;
end
ballnum(l)=ballnum(l)+1; %计数。
ballnum1=3*ballnum./m;
bar([0:n],ballnum1),axis([-2 n+2 0 y0+n+1]) %画各格子的频率。
mm(i)=getframe; %存储动画数据。
hold off
end
movie(mm,1) %播放动画1次。】
2. 用二项分布描述Galton钉板模型
【 X = binornd(5,0.5,1,10) 】
3. 数学期望和平均收益
【 n=5;p=0.5; m=5000;
rand('seed',5); R = binornd(n,p,1,m); %模拟投球。
f=[5, 1, 0.2, 0.2, 1, 5]; %奖品的价
值向量。
s=[];
for k=1:n+1 %计算
第0～n格奖品价值。
u=[];
u=find(R==(k-1)); %计算落入第k-1格的小球
下标，并存于向量u中。
s=[f(k)*length(u),s]; %计算相应的奖品价值，并
存于向量s中。
end
mean_return=sum(s)/m %计算一次抽奖的平均回
报。】
{ mean_return = 0.7506 }
【f=[5, 1, 0.2, 0.2, 1, 5]; x=[0:1:n];
pu=binopdf(x,n,p);EF=sum(f.*pu) 】
9.1.2 应用、思考和练习
1. 二项分布的应用模型
◆电力供应问题。
提示：
【x=linspace(90,160,10); r = binornd(200,0.6,1,960); hist(r,x) ;
[n,x]=hist(r,x); 】
◆废品问题。
下面的程序可作为参考。
【clear,clf ,n=50;r=linspace(0,1,n);s=linspace(0,10,n);
[R,S]=meshgrid(r,s); x=[0 1 2 3 4 5];p = binopdf(x,5,0.5);
for k=1:n
f=[S(k,l) S(k,l)*R(k,l) S(k,l)*R(k,l)^2 S(k,l)*R(k,l)^2 S(k,l)*R(k,l) S(k,l)];
for l=1:n
z(k,l)=sum(f.*p);
if z(k,l)>=1 z(k,l)=NaN; end
end
end
contour(R,S,z)】
3. 单服务台定长服务时间排队系统的计算机模拟
4. 随机变量的模拟：逆概率法
mu=0;sigma=1;a=mu-
4*sigma;b=mu+4*sigma;
t=linspace(a,b,30);
f1=normcdf(t,mu,sigma);
for i=1:3
hold on
u=rand(1);
f=norminv(u,mu,sigma);
plot(t,f1,[0 0],[0,b]),hold on,
plot(0,u,'rO'),plot([0,f],[u,u]);plot([
f,f],[u,0]);
plot(f,0,'rO'), axis([-4 4 0 1])
end 】
9.2 实验(Ⅱ) ：热轧机的调整
9.2.1实验与观察: 控制粗轧的浪费
1. 用正态分布描述热轧机模型
2. 调整额定长度使浪费最小
.观察程序
zxy9_2.m
【clf,clear, m=20;l=3;mu=3.5;sigma=0.3; mu1=3.5;sigma1=0.9;
a=0; b=max([mu+4*sigma,mu1+4*sigma1]); %设定坐标的范围。
subplot(2,2,1),
for k=1:m
rand('seed',1), x=normrnd(mu,sigma);
plot([0,x],[k,k],'linewidth',5),hold on, axis([0
mu+4*0.6 -2 m+5])
end
plot([l,l],[-2,m+5],'r-');hold on,axis([a b -2 m+5])
subplot(2,2,3),
t=linspace(a,b,50);f=normpdf(t,m u,sigma);
y0=max(f)+0.1;plot(t,f),hold on,axis([a b 0 y0]),
plot([l,l],[0,y0],'r-',[mu,mu],[0,y0],'r:');hold on
t=linspace(a,b,30);f=normpdf(t,mu,sigma);plot(t,f)
subplot(2,2,2)
for k=1:m
rand('seed',1),x=normrnd(mu1,sigma1);
plot([0,x],[k,k],'linewidth',5),hold on,axis([a b -2 m+5])
end
plot([l,l],[-2,m+5],'r-');hold on
subplot(2,2,4),
t=linspace(a,b,50);f=normpdf(t,mu1,sigma1);y0=max(f)+0.
1;
plot(t,f),hold on,axis([a b 0 y0])
plot([l,l],[0,y0],'r-',[mu1,mu1],[0,y0],'r:');hold
on 】
下面是观察与实验程序zxy9_3.m，对照上一节给出的算法，这一程序也是不难理解的。
zxy9_3.m
【 clf,clear
l=2;sigma=0.2;
n=10000;m=50;a=2.2;b=3;
mu=linspace(a,b,m);
for k=1:m
x=normrnd(mu(k),sigma,1,n); %模拟轧钢。
k_chenpin=find(x>=l); %求可轧制的成品材的下标。
k_feipin=find(x w_chenpin=x(k_chenpin)-l; %计算浪费量。 w_feipin=x(k_feipin); %计算浪费量。 if length(k_chenpin)==0 waste(k)=NaN; else waste(k)=(sum(w_chenpin)+sum(w_feipin))/length(k_chenpi n); end end [wmin,i]=min(waste); %求最小浪费量及其下标。 [mu(i) wmin] plot(mu,waste,'.- ',mu(i),wmin,'ro'),set(gca,'fontsize',15) text(mu(i),wmin,['\bullet\leftarrow The Min Value is ',num2str(wmin),' at \it{mu}= ',num2str(mu(i))]); 】 9.2.2应用、思考与练习 1. 随机优化：确定热轧机的额定长度 2. 二维正态分布: 如何制定胖和瘦的标准？ 题的思路，这种思路在实践中经常被采用，真正要确定正常体重标准可能要复杂得多。 ◆。 【 clear,clf n=1000;x=normrnd(170,4.5,1,n); y=0.36*x+normrnd(0,7,1,n); a=min(x);b=max(x);c=min(y);d=max(y);xx=linspace(a,b,20); yy=0.36*xx;plot(x,y,'ko'),hold on, plot(xx,yy,'k-','linewidth',5),grid, axis([a b c d]),axis('equal'), xlabel('身高X');ylabel('体重Y'); 】 ◆（2）观察 程序zxy9_4.m （画图9.8） 【 clear,clf,n=1000;m=15;x=normrnd(170,4.5,1,n); y=0.36*x+normrnd(0,7,1,n);a=min(x);b=max(x);c=min(y);d= max(y); h1=(b-a)/m;h2=(d-c)/m; xx=linspace(a-h1,b+h1,m);yy=linspace(c-h2,d+h2,m); [X,Y]=meshgrid(xx,yy); for k=1:m %计算频数。 for l=1:m j=find(x>=(X(k,l)-h1)&x<=(X(k,l)+h1)); h=find(y(j)>=(Y(k,l)-h2)&y(j)<=(Y(k,l)+h2)); z(k,l)=length(h)/n; end end mu=[170 0.36*170]; %计算二维正态分布密度。 C=[4.5^2 0.36*4.5^2; 0.36*4.5^2, (0.36*4.5)^2+7^2]; detC=det(C); for k=1:m for l=1:m u=[X(k,l),Y(k,l)]'-mu';s=2*pi*sqrt(detC); z1(k,l)=s^-1*exp((-0.5)*u'*C^-1*u); end end subplot(2,2,1),surf(X,Y,z),axis([a b c d 0 0.15]), set(gca,'fontsize',15); subplot(2,2,2),contour(X,Y,z),axis([a b c d ]),axis('equal'), subplot(2,2,3),surf(X,Y,z1),axis([a b c d 0 0.005]), subplot(2,2,4),contour(X,Y,z1),axis([a b c d ]),axis('equal')】 3.用线性回归方法确定正常体重标准 观察： ◆（1）参考程序 zxy9_4.m中的模拟部分，。 【alpha=0.01; polytool(x',y',1,alpha) 】 ◆（3) 用多元线性回归指令regress做体重预测。 ◆对模拟的100对身高体重数据，运行下面的程序，了解指令regress 的用法。 【[b,bint,r,rint,stats] = regress(y',[ones(100,1) x']); b,bint,stats rcoplot(r,rint),set(gca,'fontsize',15) 】 9.3实验（Ⅲ）参数估计和假设检验 9.3.1实验与观察: 极大似然估计 1.极大似然估计原理: 如何决定废品率？ 观察 ◆（1） 【 clear,p=0.04; n=10; %设定产品的档次，抽样次数。 for k=1:n %抽样n次。 r(k)=rand(1,1); %产生随机数。 if r(k)<=p x(k)=1; %抽样得到废品。 elseif r(k)>p x(k)=0; %抽样得到正品。 end end I=[1:n];[I;x] 】 ◆（2） 2.实验观察的参考程序 实验观察的主要程序为zxy9_5.m，其源代码如下，该程序是不难读懂的。 zxy9_5.m 【 clear,clf,n=50; p=0.04; rand('seed',1),r=rand(1,n); k=+find(r<=p); x=zeros(1,n); x(k)=ones(1,length(k)); p_estimate=sum(x)/n, t=[0.01 0.02 0.03 0.04 0.05 0.06]; %t=linspace(0.01,0.5,40) Lt=t.^sum(x).*(1-t).^(n-sum(x)); %lnLt=sum(x)*log(t)+(n-sum(x))*log(1-t); set(gca,'fontsize',16), plot(t,Lt,'o:'), [Lmax,I]=max(Lt); tmax=t(I); text(t(I),Lmax,['\fontsize{16}\leftarrow\it Lmax=' ,num2str(Lmax)]) text(t(I),0.95*Lmax,['\fontsize{16}\it{ at p}= ',num2str(tmax)]) xlabel('p','fontsize',16);ylabel('L(p)','font size',16); 】 9.3.2 应用、思考与练习 1. 用Matlab符号演算求解极大似然估计 ◆练习：用Matlab符号演算指令求解9.3.1节中废品问题的似然方程获得极大似然估计。 【syms p %未知参数为p,所以作为符号变量处理，用 syms指令说明。 clear,clf,n=50; %产生50个样本。 p=0.04; %设定真实参数。 x=zeros(1,n); %令x全为0 。 rand('seed',1),r=rand(1,n); k=+find(r<=p); %找出废品的下标。 x(k)=ones(1,length(k)); %在废品下标处改x为1；x为50个样本观察 值。 %观察似然函数和似然方程的一般表达式。 L=sym('p^sx*(1-p)^(n-sx)') %正确写出似然函数，L是符号p的函 数。 likely_equ=diff(L,'p') %对p进行符号求导，得到似 然方程。 %观察含具体数值的似然函数和似然方程 sx=sum(x), p='p'; %代入s x的具体数值。 Lp=subs(L) %将具体的数值代入似然函数中。 likely_equ=diff(Lp,'p') %求似然方程。】 【 %求似然方程的符号解，得到极大似然估计 s=solve('p^sx*sx/p*(1-p)^(n-sx)-p^sx*(1-p)^(n-sx)*(n- sx)/(1-p)=0','p') sx,n %看看具体的数值。 sp=subs(s) %对已获得的样本，观察极大似然估计的具体数 值。】 2. 水库入库径流的分布估计 7月上旬径流数据 356 258 222 208 163 342 501 501 782 225 630 2305 931 485 503 501 422 101 280 1807 922 390 466 211 922 444 233 370 788 802 219 524 470 1097 1160 702 566 222 630 7月中旬径流数据 98 262 117 1687 291 1318 292 716 254 519 270 273 275 274 374 147 345 70 940 440 2839 141 699 324 900 311 870 596 187 2231 111 949 303 888 328 459 370 1360 1320 7月下旬径流数据 69 133 392 596 4518 1051 336 867 541 1733 149 266 324 1365 891 918 1751 219 513 438 1033 1217 1290 247 2360 1023 453 1622 1272 1383 1217 1530 1724 703 299 638 548 1200 1220 zxy9_6.m 【 clear,clf, Q=[ 98 262 117 1687 291 1318 292 716 254 519 270 273 275 274 374 147 345 70 940 440 2839 141 699 324 900 311 870 596 187 2231 111 949 303 888 328 459 370 1360 1320]; m=50; as=linspace(0.6561, 2.2477,m); bs=linspace(215.4687, 637.4421,m); [A B]=meshgrid(as,bs); for i=1:m for j=1:m ax=A(i,j);bx=B(i,j); y(i,j) = sum(log(gampdf(Q,ax,bx))); end end [y0,s]=max(y); [y00,s0]=max(y0); a_est=A(s0,s(s0));b_est=B(s0,s(s0)); meshc(A,B,y),hold on plot3(a_est,b_est,y00,'.','markersize',25); text(a_est,b_est,y00+0.06*abs(y00),['(a_e_s_t,b_e_s_t,L _m_a_x)=']); text(a_est,b_est,y00+0.03*abs(y00),['(',num2str(a_est), ',',num2str(b_est),',',num2str(y00),')']) figure(2) [h,n0]=hist(Q,5); h0=max(h); a=min(Q)- 100;b=max(Q); x=linspace(a,b,30); hist(Q,5); hold on y = gampdf(x,a_est,b_est); plot(x,h0*y/max(y),'r-','linewidth',5); 】 3. 数学建模竞赛: 零件的参数设计 zxy9_7.m 【 clear,clf, n=1000; A=0.01/3;B=0.05/3;C=0.1/3; %x=[0.075,0.375,0.125,0.1198,1.252,12.5,0.77];%A=[B B B C C B B]; x=[0.1 0.3 0.1 0.1 1.5 16 0.75]; %设置标定值。 A=[B C C C C C B]; %设置容差。 X1=normrnd(x(1),A(1)*x(1),n,1); %模拟零件参数。 X2=normrnd(x(2),A(2)*x(2),n,1);X3=normrnd(x(3),A(3)*x(3 ),n,1); X4=normrnd(x(4),A(4)*x(4),n,1);X5=normrnd(x(5),A(5)*x(5 ),n,1); X6=normrnd(x(6),A(6)*x(6),n,1);X7=normrnd(x(7),A(7)*x(7 ),n,1); Y=z9_5fun(X1,X2,X3,X4,X5,X6,X7); %模拟y的样本。 k=find(imag(Y)==0);Y1=Y(k); %如果产生了复数，剔除复 数。 histfit(Y1) %直方 图和正态密度拟合。 dh=0.001; %求数值 导数。 for i=1:7 for j=1:7 xx(:,j)=x(j)*ones(2,1); end xx(:,i)=linspace(x(i),x(i)+dh,2)'; F(:,i)=z9_5fun(xx(:,1),xx(:,2),xx(:,3),xx(:,4),xx(:,5), xx(:,6),xx(:,7)); g(:,i)=diff(F(:,i),1)/dh; %求数值导数。 end mu=F(1,1),EY=mean(Y1), %均值对比。 R=(A.*x).^2; sigma=sqrt(sum((g.^2).*R)),DY=std(Y1), %均方差值对比. [h,p,ci]=ttest(Y1,mu,0.01,0) %均值的t－检验。】 zxy9_6fun.m（计算函数的子程序） 【 function y=zxy9_5fun(x1,x2,x3,x4,x5,x6,x7) y1=174.42*(x1./x5).*(x3./(x2-x1)).^0.85; y2=(1-0.36*(x4./x2).^(-0.56)).^(1.5) ; y3=(1-2.62*y2.*(x4./x2).^1.16)./(x6.*x7); y=y1.*sqrt(y3); 】 实验10.数值仿真 10.1知识要点与背景：单自由度阻尼系统 2.观察程序 zxy10_1.m （图10.1(a)） 【 clear;clf; global c w x0(1)=1;x0(2)=0;w=10;n=3;tspan=linspace(0,4,100);cc=[0.1 0.4 0.7 1]; xx=[];hold on，，xlabel('t'),ylabel('x'), for i=1:length(cc); c=cc(i); [t,x]=ode45('zxy10_1f ',tspan,x0); %求微分方程的数值解。 xx=[x(:,1),xx]; text(t(10),x(10,1),['\leftarrow c=',num2str(c)],'fontsize',15) plot(t,x(:,1)), end 】 zxy10_1f.m od45指令中的微分方程组函数子程序 【 function dx=zxy10_2f(t,x) global c w u=0;dx=zeros(2,1); dx(1)=x(2);dx(2)=u-c*w*x(2)-w^2*x(1); 】 zxy10_2.m （图10.1(b)） 【 global c w x0(1)=1;x0(2)=0;w=10;n=3;tspan=linspace(0,4,1500);cc=1:-1/10:0;
相关主题
文本预览