数学实验报告
姓名： 学号： 班级：
姓名： 学号： 班级：
姓名： 学号： 班级：
1、整数规划类型（选自课本P154，上机练习题第7题）
【实验问题】：某轰炸机群轰炸敌人军事目标。已知该目标有四个要害部位，只要摧毁其中之一即可达到目标。为完成此项任务的汽油消耗量限制为48000升、重型炸弹48枚、轻型炸弹32枚。飞机携带重型炸弹时每升汽油可飞行2千米，带轻型炸弹时每升汽3千米。又已知每架飞机每次只能装载一枚炸弹，每出发轰炸一次除来回路程汽油消耗（空载时每升汽油可飞行4千米）外，起飞和降落每次各消耗100升。有关数据如表所示：
为了使摧毁敌方军事目标的可能性最大，试确定飞机轰炸的方案，要求建立这个问题的线性规划模型。
【问题分析】
【总结和体会】
汽油限制为48000升，故有：
537.5X11+560 X12+605 X13+650 X14+462.5 X21+480 X22+515 X23+550 X24
48000---------1
再考虑炸弹数量限制，重型炸弹为48枚，轻型为32枚，即有：
X11 +X12 +X13 +X1448 -----------------------------------------------------------------------------------2
X21 +X22 +X23 +X2432------------------------------------------------------------------------------------3
另外，炸弹数量定大于零，故有：
Xij
0 （i=1,2；j=1，2,3）-------------------------------------------------------------------------------4
537.5X11+560 X12+605 X13+650 X14+462.5 X21+480 X22+
515 X23+550 X2448000
s.t. X11 +X12 +X13 +X1448
X21 +X22 +X23 +X2432
Xij
0 （i=1,2；j=1，2,3）
2、多目标类型（选自）
【实验问题】某工厂生产A1、A2、A3三种产品以满足市场的需要，该厂每周生产的时间为40h,其数据如下表所示。
（1） 请确定最佳时间分配方案，使得该工厂一周内获得的利润最大。
（2） 由于能源紧缺，现要求每周能耗不得超过20t标准煤，三种产品单位能耗如下表：
【问题分析】
假设三种产品每周分配的小时数分别为x1、x2、x3，
则该厂一周内得到的总利润f1(x)为：f1(x)= 500x1+400 x2+600x3 .
该厂一周内的总能耗f2(x)为:f2(x)= 0.48 x1+0.65 x2+0.42 x3
由约束条件，一周生产时间不超过40小时，故：x1+x2+x340 1
再考虑销售过程，由于产品数量应不超过本周最大销量才能使成本最低，故有：
40 x1500 2
25 x2 400 3
15 x3600 4
由约束条件 1 2 3 4 可建立如下模型：
Max f1(x)= 500x1+400 x2+600x3
s.t. x1+x2+x340
40 x1500
25 x2 400
15 x3600
x1、x2、x3
0
Matlab 程序代码为：
a=[-500,-400,-600];
b=[1,1,1;
40,0,0;
0,25,0;
0,0,15;
-1,0,0;
0,-1,0;
0,0,-1];
c=[40;700;800;500;0;0;0];
[x,fal]=linprog(a,b,c)
程序运行结果为：
Optimization terminated.
x = 6.6667
0.0000
33.3333
fal = -2.3333e+004
即 ：目标函数的最优解为：x1 =6.6667 ,x2 =0.0000 ,x3 =33.3333 ;
最大利润为23333元
（2）
Max f1(x)= 500x1+400 x2+600x3
Min f2(x)= 0.48 x1+0.65 x2+0.42 x3
s.t. x1+x2+x340
40 x1500
25 x2 400
15 x3600
x1、x2、x30
此问题中，目标函数一个为求极大，一个为求极小，将其均转换成求极小，编写M-函数文件fun.m如下：
function f=fun(x)
f(1)=-500*x(1)-400*x(2)-600*x(3);
f(2)=0.48*x(1)+0.65*x(2)+0.42*x(3);
然后设定期望目标goal,我们假设f2(x)乘以1500后的值显然比f1(x)的值大，而f2(x)20,故f1(x)30000,我们设定f1(x)的期望目标为30000，此值不能达到，希望求解过程尽量接近该值，于是设定初始点x0=[0;0;0],期望目标goal=[-30000,20].权重为期望目标的绝对值，然后再根据相应的约束确定A和b,调用fgoalattain求解该多目标规划问题，程序代码为：
x0=[0;0;0];
lb=[0;0;0];
goal=[-30000;20];
weight=abs(goal);
a=[1,1,1;
40,0,0;
0,25,0;
0,0,15];
b=[40 700 800 500];
[x,favl,attainfactor]=fgoalattain(@objfun,x0,goal,weight,a,b,[],[],lb,[])
运行结果为：
Local minimum possible. Constraints satisfied.
fgoalattain stopped because the predicted change in the objective function
is less than the default value of the function tolerance and constraints
were satisfied to within the default value of the constraint tolerance.
<stopping criteria details>
Active inequalities (to within options.TolCon = 1e-006):
lower upper ineqlin ineqnonlin
1 1
4
x =
1.4479
5.2188
33.3333
favl =
1.0e+004 *
-2.2811 0.0018
attainfactor =
0.2396
【总结和体会】
分析第二道实验题目时，第一问属于单目标规划问题，相对容易解决 ; 第二问属于多目标规划，由于需要同时考虑利润和能耗，利润越大，能耗也就越大，能耗越小，利润也随之减小，使问题难度增大。但经过仔细分析知道，我们不妨令能耗乘以一个较大数使之大于利润的最高期望值，于是设定初始点x0=[0;0;0],期望目标goal=[-30000,20].权重为期望目标的绝对值，然后再根据相应的约束确定A和b，调用Matlab 函数 fgoalattain求解该多目标规划问题。
第二篇：matlab数学实验实例
《数学实验与Matlab》程序
周晓阳
华中科技大学数学系
我将程序按书中的顺序排列，这样便于读者利用。
本书程序均通过了调式。可直接拷贝到命令窗口运行或编制M文件运行。
如出现问题，可能是中英文标点的缘故（排版有可能使用了中文标点），请将中文标点换为英文标点试试。
实验１．矩阵运算与Matlab命令
1.1 知识要点与背景：知识要点和背景：日常矩阵及其运算
【 A=[4 2 3;1 3 2;1 3 3;3 2 2], % 表1-1、表1-2的数据分别写成矩阵形式
B=[35 20 60 45;10 15 50 40;20 12 45 20] 】
【 C=A*B %矩阵乘法，求各订单所对应的原材料和劳动力 。 】
【 whos % 查看Matlab工作空间中变量及其规模 】
1.2实验与观察：矩阵和Matlab语言
1.2.1 向量的生成和运算
【 x=linspace(0,4*pi,100); %将[0,4π]区间100等分，产生了一个100维向量
y=sin(x); %计算函数值，产生了一个与x同维的100维函数向量y
y1=sin(x).^2; %计算函数向量，注意元素群运算
y2=exp(-x).*sin(x);
%以x为横坐标，y为纵坐标画函数的图用不同的线型将函数曲线绘制在一个图上
plot(x,y,'-',x,y1,'-',x,y2,'.-') 】
1. 向量的创建
◆直接输入向量。
【x1=[1 2 4],x2=[1,2,1],x3=x1' 】
◆冒号创建向量 。
【 x1=3.4:6.7
x2=3.4:2:6.7
x3=2.6:-0.8:0 】
◆生成线性等分向量。
【 x=linspace(0,1,5) 】
2. 向量的运算
【 y=sin(x) 】
【 y1=sin(x).^2; y2=exp(-x).*sin(x); 】
1.2.2.矩阵创建和运算
1.创建矩阵
（1）数值矩阵的创建
◆直接输入法创建简单矩阵。
【 A=[1 2 3 4; 5 6 7 8; 9 10 11 12] 】
【 B=[-1.3,sqrt(3);(1+2)*4/5,sin(5);exp(2),6] 】
（2）符号矩阵的创建
◆
【 syms a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 …
b11 b12 b13 b14 b21 b22 b23 b24 b31 b32 b33 b34
A1=[a11 a12 a13 a14 ;a21 a22 a23 a24; a31 a32 a33 a34],
B1=[b11 b12 b13 b14 ;b21 b22 b23 b24; b31 b32 b33 b34] 】
2.矩阵的运算
【 C=A1+B1,D=A1-B1 】
【 syms c
cA=c*A1 】
【 C=A1*B1 】
{ ??? Error using ==> sym/mtimes, Inner matrix dimensions must agree. }
【 A2=A1(:,1:3), B1 】
【 G=A2*B1 】
【 g11=A2(1,:)*B1(:,1) 】
【 A, A_trans=A' 】
【 H=[1 2 3 ; 2 1 0 ; 1 2 3 ], K=[1 2 3 ; 2 1 0 ; 2 3 1]
h_det=det(H), k_det=det(K),H_inv=inv(H),K_inv=K^-1 】
【 A=[3 0 1; 1 1 0;0 1 4];
B=inv(A-2*eye(3))*A, B=(A-2*eye(3))\A 】
3.分块矩阵：矩阵的裁剪、分割、修改与抽取
（1）
【 A=[1 0 1 1 2;0 1 -1 2 3;3 0 5 1 0;2 3 1 2 1], vr=[1,3];vc=[1,3];
A1=A(vr,vc) %取出A的1、3行和1、3列的交叉处元素构成新矩阵A1 】
◆将上面的矩阵A分为四块，并把它们赋值到矩阵B中，观察运行后的结果。
【 A11=A(1:2,1:2),A12=A(1:2,3:5),A21=A(3:4,1:2),A22=A(3:4,3:5)
B=[A11 A12;A21 A22] 】
A =
2 0 5
4 2 1
0 -1 2
B =
1 2 4 -1
5 3 1 0
-1 0 2 3
C =
-3 4 18 13
13 14 20 -1
-7 -3 3 6
（2）矩阵的修改和提取
◆ 【 A=[1 0 1 1 2;0 1 -1 2 3;3 0 5 1 0;2 3 1 2 1]
A(1,:)=[0 0 0 0 0]; A 】
◆ 观察：
【 B(:,[2,4])=[ ] %删除矩阵B的第2、4列 】
（3）矩阵元素的抽取
4.生成特殊矩阵
。
◆
【 y1=rand(1,5), y2=rand(1,5),
rand('seed',3), x1=rand(1,5), rand('seed',3), x2=rand(1,5) 】
5. 常用矩阵函数
6. 数据的简单分析
◆
【 rand('seed',1);A=rand(3,6),
Asort=sort(A), Amax=max(A), Asum=sum(A) 】
1.2.3 Matlab工作环境和编程
2.Matlab的基本设计
1.3应用、思考与练习
1.3.1 关系矩阵
1.3.2 投入产出
1.3.3 循环比赛的名次
【 A=[0 1 1 0; 0 0 1 1; 0 0 0 1; 1 0 0 0],
e=ones(4,1); c=A*e;
s=c' 】
★ 画矩阵结构图的gplot指令。
◆（3）
【 clf, A=[0 1 1 0;0 0 1 1;0 0 0 1;1 0 0 0]; xy=[0 1;0 0;-1 –0.5;1 –0.5];
graphy_plot(A,xy,1,0.5), % gplot(A,xy) 】
1.3.4 参考程序
graphy_plot.m
【 function y=graphy_plot(A,xy,l,p)
%画矩阵的有向结构图。 A为邻接矩阵，xy为顶点坐标，l控制参数，l=0，画无向图；
%l~=0,画有向图。p为控制箭头大小的参数。
a=-max(abs(xy(:,1)))*1.1;b=max(abs(xy(:,1)))*1.1;
c=-max(abs(xy(:,2)))*1.1;d=max(abs(xy(:,2)))*1.1;
if l=0
gplot(A,xy),axis([a b c d]),hold on,
elseif l~=0
U=[];V=[];X=[];Y=[];
n=length(A(:,1)) ;
for i=1:n
k=find(A(i,:)~=0);
m=length(k);
if(m~=0)
for j=1:m
u(1)=(xy(k(j),1)-xy(i,1)); v(1)=(xy(k(j),2)-xy(i,2));
u(2)=eps; v(2)=eps; U=[u;U];V=[v;V];
X=[[xy(i,1) xy(k(j),1)];X]; Y=[[xy(i,2) xy(k(j),2)];Y];
end
text(xy(i,1),xy(i,2),['\bullet\leftarrow\fontsize{16}\it{V}',…
um2str(i)]); hold on,
end
end
gplot(A,xy),axis([a b c d]),hold on,
h=quiver(X,Y,U,V,p);set(h,'color','red');hold on,
plot(xy(:,1),xy(:,2),'k.','markersize',12),hold on,
end , hold off 】
实验２．函数的可视化与Matlab作
2.1 实验与观察：函数的可视化
2.1.1 Matlab二维绘图命令
1.周期函数与线性p－周期函数
◆ 观察 ：
【 clf, x=linspace(0,8*pi,100);
F=inline('sin(x+cos(x+sin(x)))');
y1=sin(x+cos(x+sin(x))); y2=0.2*x+sin(x+cos(x+sin(x)));
plot(x,y1,'k:',x,y2,'k-') legend('sin(x+cos(x+sin(x))','0.2x+sin(x+cos(x+sin(x)))',2) 】
2. plot指令：绘制直角坐标的二维曲线
3. 图形的属性设置和屏幕控制
【 h=plot([0:0.1:2*pi],sin([0:0.1:2*pi])); grid on
set(h,'LineWidth',5,'color','red'); set(gca,'GridLineStyle','-','fontsize',16) 】
◆设置y坐标的刻度并加以说明，并改变字体的大小。
【 h=plot([0:0.1:2*pi],sin([0:0.1:2*pi]));grid on，
set(gca,'ytick',[-1 -0.5 0 0.5 1]), set(gca,'yticklabel','a|b|c|d|e'),
set(gca,'fontsize',20) 】
4. 文字标注指令
【 plot(x,y1,'b',x,y2,'k-') ,
set(gca,'fontsize',15,'fontname','times New Roman'), %设置轴对象的字体为times
% New Roman，字体的大小为15
title(' \it{Peroid and linear peroid function}'); %加标题，注意文字用单引号' ' 加上
%斜杠'\'后可输入不同的设置，例如it{…}表示花括号里的文字为斜体；如果有多项设置，
%则可用\…\…\…连续输入。
xlabel('x from 0 to 8*pi it{t}\'); ylabel('\it{y}'); %说明坐标轴
text(x(49),y1(50)-0.4,'\fontsize{15}\bullet\leftarrowThe period function {\itf(x)}');
%在坐标(x(49),y1(50)-0.4）处作文字说明， 各项设置用"\"隔开。
%\fontsize{15}\bullet\leftarrow的意义依次是：\字体大小=15 \ 画圆点 \左箭头
text(x(14),y2(50)+1,'\fontsize{15}The linear period function {\itg(x)}\rightarrow
\bullet') %与上一语句类似，用右箭头 】
◆观察指令legend和num2str的用法：在同一张图上画出
， 这里
, 并进行适当的标注。
zxy2_2.m
【 clf, t=0:0.1:3*pi;alpha=0:0.1:3*pi;
plot(t,sin(t),'r-');hold on; plot(alpha,3*exp(-0.5*alpha),'k:');
set(gca,'fontsize',15,'fontname','times New Roman'),
xlabel('\it{t(deg)}');ylabel('\it{magnitude}');
title(' \it{sine wave and {\it{Ae}}^{-\alpha{\itt}}wave}'); %注意\alpha的意义
text(6,sin(6),'\fontsize{15}The Value \it{sin(t)} at {\itt}=6\rightarrow\bullet', 'HorizontalAlignment','right'),
%上面的语句是一整行，如果要写成两行，必须使用续行号 … ，例如要在“ bullet',”
%后换行，需写“bullet', …”后才能换行。
% 'HorizontalAlignment','right' 表示箭头所指的曲线对象在 文字的右边。
text(2,3*exp(-0.5*2),['\fontsize{15}\bullet\leftarrow The Value of \it{3e}^{-0.5 \it{t}}=',num2str(3*exp(-0.5*2)),' at \it{t} =2 ']);
%num2str的用法：['string1' ,num2str,'string2']，注意方括号的使用。
%legend('\itsin(t)','{\itAe}^{-\alphat}') % 请结合图形观察此命令的使用 】
运行结果如图2.6所示。
5. 图形窗口的创建和分割
◆观察：
【 clf，b=2*pi;x=linspace(0,b,50);
for k =1:9
y=sin(k*x);
subplot(3,3,k),plot(x,y),axis([0,2*pi,-1,1])
end 】
x=linspace(0,4*pi,150);
y1=cos(x);
y=cos(cos(cos(cos(cos(cos(cos(x)+sin(x)))+sin(x)))+sin(x)));
plot(x,y,x,y1,'r-'),grid
2.1.2多元函数的可视化与空间解析几何（三维图形）
本节通过高等数学的几个例子观察Matlab的三维绘图功能和技巧。
1. 绘制二元函数
◆观察：绘制 
的图象，作定义域的裁剪。
◆（1）观察meshgrid指令的效果。
【 a=-0.98;b=0.98;c=-1;d=1;n=10;
x=linspace(a,b,n); y=linspace(c,d,n);
[X,Y]=meshgrid(x,y);
plot(X,Y,'+') 】
★三维绘图指令mesh、meshc、surf。
◆（2）做函数的定义域裁剪，观察上述三维绘图指令的效果。
程序zxy2_4.m
【 clear,clf,
a=-1;b=1;c=-15;d=15;n=20;eps1=0.01;
x=linspace(a,b,n);y=linspace(c,d,n);
[X,Y]=meshgrid(x,y);
for i=1:n %计算函数值z ，并作定义域裁剪
for j=1:n
if (1-X(i,j))<eps1|X(i,j)-Y(i,j)<eps1 %if语句这样用
z(i,j)=NaN; %作定义域裁剪，定义域以外的函数值为NaN
else
z(i,j)=1000*sqrt(1-X(i,j))^-1.*log(X(i,j)-Y(i,j));
end
end
end
zz=-20*ones(1,n);plot3(x,x,zz),grid off,hold on %画定义域的边界线
mesh(X,Y,z) %绘图，读者可用meshz, surf, meshc在此替换之
xlabel('x'),ylabel('y'),zlabel('z'), box on %把三维图形封闭在箱体里 】
◆运行了zxy2_4.m 以后，有关向量存储在工作空间中，在此基础上，观察上述等值线绘制指令的运行效果。
【 [cs,h]=contour(X,Y,z,15); clabel(cs,h,'labelspacing',244) 】
2. 三元函数可视化: slice指令
◆ 观察： 绘制三元函数
的可视化图形。
【 clf,x=linspace(-2,2,40); y=x; z=x;
[X,Y,Z]=meshgrid(x,y,z); w=X.^2+Y.^2+Z.^2;
slice(X,Y,Z,w,[-1,0,1],[-1,0,1],[-1,0,1]),colorbar 】
3. 空间曲线及其运动方向的表现：plot3和quiver指令
【 clf, t=0:0.1:1.5;
Vx=2*t;Vy=2*t.^2;Vz=6*t.^3-t.^2;
x=t.^2;y=(2/3)*t.^3;z=(6/4)*t.^4-(1/3)*t.^3; %由速度得到曲线
plot3(x,y,z,'r.-'),hold on %画飞行轨迹
%算数值梯度，也就是重新计算数值速度矢量，这只是为了编程的方便，不是必须的
Vx=gradient(x);Vy=gradient(y);Vz=gradient(z);
quiver3(x,y,z,Vx,Vy,Vz),grid on %画速度矢量图
xlabel('x'),ylabel('y'),zlabel('z') 】
2.2应用、思考和练习
2.2.1 线性p周期函数
zxy2_3_f.m
【 function f=zxy2_3_f(x)
f=sin(x+cos(x)); 】
zxy2_3.m
【 clear，clf
a=-8;b=12;n=300;xx=linspace(a,b,n);
h=zxy2_3_f(xx);
S(1)=0;
for i=2:n
S(i)=S(i-1)+quad('zxy2_3_f',xx(i-1),xx(i));
end
subplot(1,2,1),plot(xx,S,'k-'),axis([a,b,-1.5,9])
subplot(1,2,2),plot(xx,[h;zeros(1,length(xx))],'k-'),axis([a,b,-1.5,1.5]) 】
2.2.2 平面截割法和曲面交线的绘制
◆用平行截面法讨论由曲面 
构成的马鞍面形状。
zxy2_6.m （ 平行截割法示例 ， 本程序的绘制两曲面交线方法可以套用）
【 clf, a=-20;eps0=1;
[x,y]=meshgrid(-10:0.2:10); %生成平面网格
v=[-10 10 -10 10 -100 100]; %设定空间坐标系的范围
colormap(gray) %将当前的颜色设置为灰色
z1=(x.^2-2*y.^2)+eps; %计算马鞍面函数z1=z1(x,y)
z2=a*ones(size(x)); %计算平面 z2=z2(x,y)
r0=abs(z1-z2)<=eps0;
%计算一个和z1同维的函数r0,当abs(z1-z2)<=eps时r0 =1；当abs(z1-z2)>eps0时，r0 =0。
%可用mesh(x,y,r0)语句观察它的图形，体会它的作用，该方法可以套用。
zz=r0.*z2;xx=r0.*x;yy=r0.*y; %计算截割的双曲线及其对应的坐标
subplot(2,2,2), %在第2图形窗口绘制双曲线
h1=plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'+');
set(h1,'markersize',2),hold on,axis(v),grid on
subplot(2,2,1), %在第一图形窗口绘制马鞍面和平面
mesh(x,y,z1);grid,hold on;mesh(x,y,z2);
h2=plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'.'); %画出二者的交线
set(h2,'markersize',6),hold on,axis(v),
for i=1:5 %以下程序和上面是类似的，通过循环绘制一系列的平面去截割马鞍面
a=70-i*30; %在这里改变截割平面
z2=a*ones(size(x)); r0=abs(z1-z2)<=1; zz=r0.*z2;yy=r0.*y;xx=r0.*x;
subplot(2,2,3),
mesh(x,y,z1);grid,hold on;mesh(x,y,z2);hidden off
h2=plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'.'); axis(v),grid
subplot(2,2,4),
h4=plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'o');
set(h4,'markersize',2),hold on,axis(v),grid on
end 】
2.2.3 微分方程的斜率场
◆ 绘制微分方程 
的斜率场，并将解曲线画在图中，观察斜率场和解曲线的关系。
zxy2_5.m ( 绘制一阶微分方程的斜率场和解曲线)
【 clf,clear %清除当前所有图形窗口的图像，清除当前工作空间的内存变量。
a=0;b=4;c=0;d=4;n=15;
[X,Y]=meshgrid(linspace(a,b,n),linspace(c,d,n)); %生成区域中的网格。
z=X.*Y; %计算斜率函数。
Fx=cos(atan(X.*Y));Fy=sqrt(1-Fx.^2); %计算切线斜率矢量。
quiver(X,Y,Fx,Fy,0.5),hold on,axis([a,b,c,d])
%在每一网格点画出相应的斜率矢量，0.5是控制矢量大小的控制参数，可以调整。
[x,y]=ode45('zxy2_5f',[0,4],0.4); %求解微分方程。
%zxy2_5f.m是方程相应函数f(x,y)的程序，单独编制；[x0,xs]=[0,4]为求解区间；
%y0=0.4为初始值；输出变量x,y分别为解轨线自变量和因变量数组。
plot(x,y,'r.-') %画解轨线 】
zxy2_5f.m （微分方程的函数子程序）
【 function dy=zxy2_5f(x,y)
dy=x.*y; 】
2.2.4颜色控制和渲染及特殊绘图指令
1.地球表面的气温分布（sphere指令）
◆
【 [a,b,c]=sphere(40);t=max(max(abs(c)))-abs(c);surf(a,b,c,t);
axis('equal'),colormap('hot'), shading flat,colorbar 】
2.旋转曲面的生成：柱面指令cylinder和光照控制指令surfl
◆
【 x=0:0.1:10;z=x;y=1./(x.^3-2.*x+4);
[u,v,w]=cylinder(y);surfl(u,v,w,[45,45]);
shading interp 】
3.若干特殊图形
◆ 运行下面程序，了解各指令的用法和效果。
【 x=[1:10]; y=[5 6 3 4 8 1 10 3 5 6];
subplot(2,3,1),bar(x,y),axis([1 10 1 11])
subplot(2,3,2),hist(y,x),axis([1 10 1 4])
subplot(2,3,3),stem(x,y,'k'),axis([1 10 1 11])
subplot(2,3,4),stairs(x,y,'k'), axis([1 10 1 11])
subplot(2,3,5), x = [1 3 0.5 5];explode = [0 0 0 1];pie(x,explode)
subplot(2,3,6),z=0:0.1:100; x=sin(z);y=cos(z).*10;
comet3(x,y,z) 】
实验３．函数式－直接确定型模型
3.1知识要点和背景：函数 — 直接确定性模型
3.2实验与观察：插值与拟合
3.2.1 插值方法与多项式拟合的概念
3.2.2 用Matlab作插值和拟合
3.2.3 用鼠标选节点 观察插值、拟合的效果
3.2.4 观察程序说明
zxy3_1.m
【 clf,a=-1;b=1;n=100;
%用内联函数inline命令定义函数，
%在后面可直接用于函数g的计算,要改变函数做实验，可按此格式重新定义g
g=inline('x^2-x^4'); xx=linspace(a,b,n);
for i=1:n
gx(i)=gxx(i)); % 前面已经用inline命令定义了g，可以这样用g计算函数值
end
ymin=min(gx)*0.8;ymax=max(gx)*1.2;%分四个界面画图g的图形，以便于结果比较
subplot(2,2,1),
plot(xx,gx,'--'),grid,hold on,axis([a b ymin ymax]),title('近邻插值')
subplot(2,2,2),
plot(xx,gx,'--'),grid,hold on,axis([a b ymin ymax]),title('线性插值')
subplot(2,2,3),
plot(xx,gx,'--'),grid,hold on,axis([a b ymin ymax]),title('样条插值')
subplot(2,2,4),plot(xx,gx,'--'),
grid,hold on,axis([a b ymin ymax]),title('多项式拟合')
%用鼠标在屏幕上选点[x,y,button] = ginput(n)，可套用下面程序的格式
button=1;
x1=[a];y1=[gx(1)];
while button==1
[xi,yi,button]=ginput(1);
subplot(2,2,1),h=plot(xi,yi,'ro') %在4个图形窗口画点
subplot(2,2,2),h=plot(xi,yi,'ro')
subplot(2,2,3),h=plot(xi,yi,'ro')
subplot(2,2,4),h=plot(xi,yi,'ro')
x1=[xi,x1];y1=[yi,y1]; %将选的点存于向量x1,y1
end
x1=[b,x1];y1=[gx(n),y1];
xx=linspace(a,b,n); %定义自变量xx
%计算不同的插值函数：x1,y1为节点，xx为输入自变量
ynearest=interp1(x1,y1,xx,'nearest');
ylinear=interp1(x1,y1,xx,'linear');
yspline=interp1(x1,y1,xx,'spline');
%多项式拟合指令[p,s] = polyfit(x,y,n)，n为拟合多项式次数，x,y为被拟合数据，
%p为拟合多项式的系数，s是用来做误差 估计和预测的数据结构。
[p,c]=polyfit(x1,y1,4);
ypolyfit=polyval(p,xx); %用polyval(p,x)计算系数为p的多项式在标量或向量x处的值
subplot(2,2,1),h=plot(xx,ynearest,'r-');set(h,'linewidth',2) %画图
subplot(2,2,2),h=plot(xx,ylinear,'r-');set(h,'linewidth',2);
subplot(2,2,3),h=plot(xx,yspline,'r-');set(h,'linewidth',2)
subplot(2,2,4),h=plot(xx,ypolyfit,'r-');set(h,'linewidth',2) 】
3.3应用、思考和练习
3.3.1若干函数的插值和拟合练习
3.3.2几个应用问题
1. 机床加工和水深流速问题
2. 内燃机转角与升程的关系
3. 随高度变化的大气压强
4. 交通事故的调查
3.3.4多元函数的插值
plot(fi,x,'o')
fi=[91 105 110 115 120 124 128];
x=[0 1.0869 1.9710 2.7555 3.3986 4.9073 8.3409];
plot(fi,x,'o')
x =
0 1.0869 1.9710 2.7555 3.3986 4.9073 8.3409
1. 矩形温箱的温度
2. 航行区域的警示线
3.3.5 Fourier级数和周期函数的经验公式
zxy3_2.m
【 clf,clear,
n=10;m=3;x=1:1:12;
y=[3.1 3.8 6.9 12.7 16.8 20.5 24.5 25.9 22.0 16.1 10.7 5.4];
z=(pi/6)*x;plot(z(1:12),y(1:12),'o');hold on
k=1:m; %计算数据矩阵。
for i=1:n
X(i,:)=[1 cos(k*z(i)) sin(k*z(i))];
end
c=inv(X'*X)*X'*y(1:n)', %求解。
z1=linspace(0,2*pi,30);
s=[]; %开始计算F－级数的部分和。
for i=1:30;
sd=[1 cos(k*z1(i)) sin(k*z1(i))]'; %注意k是向量。
sd=c.*sd; s=[s,sum(sd)];
end
plot(z1,s,'r-'),hold on, xlabel('月份'),ylabel('平均气温') 】
3.3.6由实验到理论:从开普勒定律和牛顿万有引力定律
x=[0 3 5 7 9 11 12 13 14 15];
y=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6];
x1=[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9];
y1=[3.1950, 3.2299, 3.2532, 3.2611, 3.2516, 3.2282, 3.1870, 3.1266, 3.0594, 2.9757];
clf,subplot(1,2,1),plot(x,y,'k.'),axis([-0.5 16 -0.5 4]),subplot(1,2,2),plot(x1,y1,'k.') T=[0.241 0.615 1.000 1.881 11.862 29.457];
R=[0.387 0.723 1.000 1.524 5.203 9.539];
T=T.^2, R=R.^3 plot(T,R,'o')
―――――――――――――――――――――――――
运行结果（图形略，显示为直线）
T = 0.0581 0.3782 1.0000 3.5382 140.7070 867.7148
R = 0.0580 0.3779 1.0000 3.5396 140.8515 867.9777
实验４. 微分、积分和微分方程
4.1. 知识要点和背景：微积分学基本定理
对于定义在[a,b]区间上的函数y=f(x),将[a,b]区间n等分，步长h=(b-a)/n.设分点为a=x0<x1<x2<…<xn=b,y0=f(x0),y1=f(x1),y2=f(x2),…,yn=f(xn).则
（2.2.1）
(2.2.2)
4.2 实验与观察（Ⅰ）：数值微积分
4.2.1实验：积分定义、微分方程和微积分基本定理的联系
◆步骤1.
【 n=4;n1=n+1;x=linspace(0,pi,n1);f=myfun2_2(x);[0:n;x;f] 】
◆步骤2.
【 y(1)=0; y(2)=y(1)+f(1)*(x(2)-x(1));
P_intial=[x(1),y(1)],P_final=[x(2),y(2)], 】
◆步骤3.
【 y(3)=y(2)+f(2)*(x(3)-x(2));
P_intial=[x(2),y(3)], P_final=[x(3),y(3)] 】
◆ 步骤4。
【 Dy=y(3)-y(2) 】
【 DS=f(2)*(x(3)-x(2)) 】
◆步骤5 .
4.2.2 求解数值积分的Matlab积分命令
◆观察cumsum指令的功能。
【 x=[1 1 1 1 1]；I=cumsum(x) 】
对于矩阵X，cumsum(X,dim)返回一与X 相同维数的矩阵，cumsum(X)＝cumsum(X,1)返回X的按行累积和(第１维)．cumsum(X,2)返回X的按列累积和（第２维）．
>>x=[1 1 1 ;1 1 1;1 1 1],I1=cumsum(x,1),I2=cumsum(x,2)
x =
1 1 1
1 1 1
1 1 1
I1 =
1 1 1
2 2 2
3 3 3
I2 =
1 2 3
1 2 3
1 2 3
◆观察：用累积和命令cumsum计算积分。
【 x=linspace(0,pi,50); y=sin(x);
T=cumsum(y)*pi/(50-1);I=T(50) 】
◆观察梯形公式计算的效果。
【 x=linspace(0,pi,50); y=sin(x);T=trapz(x,y) 】
◆ 观察辛普森积分公式的计算效果。
【 I1=quad('sin',0,pi), I2=quad8('sin',0,pi), 】
◆ 观察：用解常微分方程的方法计算积分
。
【 y0=0;[x,y]=ode23('myfun2_2',[0,pi],y0);s=length(y);y(s) 】
【 y0=0;[x,y]=ode45('myfun2_2',[0,pi],y0);s=length(y);y(s) 】
4.2.3 观察程序及其说明
zxy4_1.m （观察方程的解曲线族，图4.1）
【 clf,clear,a=0;b=pi;c=0;d=2.2;n=20;
[X,Y]=meshgrid(linspace(a,b,n),linspace(c,d,n));
z=X.*Y; Fx=cos(atan(sin(X)));Fy=sin(atan(sin(X)));
quiver(X,Y,Fx,Fy,0.5),hold on,axis([a,b,c,d])
[x,y]=ode45('zxy4_1f',[0,pi],0);
[x1,y1]=ode45('zxy4_1f',[0,pi],0.2);
[x2,y2]=ode45('zxy4_1f',[0,pi],0.4);
[x3,y3]=ode45('zxy4_1f',[0,pi],0.6);
plot(x,y,'r.-',x1,y1,'r.-',x2,y2,'r.-',x3,y3,'r.-'), xlabel('x') ,ylabel('y') 】
zxy4_1f.m
【 function dy=zxy4_1f(x,y)
dy=sin(x); 】
. 观察 程序zxy4_2.m （图4.1～4.3）
【 clf, a=0;b=4*pi;n=31;
x=linspace(a,b,n); y=zeros(1,n); y(1)=0; s(1)=0;
for i=1:n-1
dy=myfun2_2(x(i));
y(i+1)=y(i)+dy*(x(i+1)-x(i));
hh(i)=myfun2_2(x(i));
s(i+1)=s(i)+hh(i)*(x(i+1)-x(i));
end
a1=0.9*a;b1=1.1*b; % 设置坐标范围。
ymin=min(y');ymax=max(y');
ymin1=ymin*0.9;ymax1=ymax*1.1;
subplot(2,1,1), %在第一幅图中画垂线,和原函数。
fplot('1-cos(x)',[0,pi],'r:');axis([a1 b1 ymin1 ymax1]),hold on,
plot([x(1) x(1)],[ymin ymax]);
subplot(2,1,2), %在第二幅图中画被积函数图象。
fplot('myfun2_2',[a b],'-k');hold on,axis([a1 b1 ymin1 ymax1])
for i=1:n-1 %在第I个子区间上计算。
subplot(2,1,1),
plot([x(i+1) x(i+1)],[ymin ymax]); %在各分点画竖线。
plot([x(i) x(i+1)],[y(i),y(i)]); %画水平线。
h=plot([x(i+1) x(i+1)],[y(i),y(i+1)]); %画表示增量的垂线。
set(h,'linewidth',3,'color','b'); %设置线宽和颜色。
h1=plot([x(i) x(i+1)],[y(i) y(i+1)],'.-'); %画折线，设置图形属性 。
set(h1,'linewidth',2,'markersize',18);
subplot(2,1,2),
xx=[x(i) x(i) x(i+1) x(i+1)]; yy=[0 hh(i) hh(i) 0]; %计算矩形顶点坐标。
patch(xx,yy,[0.7 0.7 0.7]); %在第二幅图中画矩形块并填充颜色。
plot([x(i+1) x(i+1)],[ymin ymax]);
plot([x(1) x(1)],[s(i),s(i+1)],'r','linewidth',5); %在y轴上画面积增量线段。
plot(x(1),y(i+1),'r.','Markersize',18);
subplot(2,1,1),plot([x(i+1) x(i+1)],[ymin ymax]);
h=plot([x(1) x(1)],[s(i),s(i+1)]); %画相应的辅助线段。
set(h,'linestyle','-','linewidth',5,'color','red');
plot(x(1),y(i+1),'r.','Markersize',18);
plot([x(1) x(i+1)],[s(i+1) y(i+1)],'r--')
pause %暂停，n大时应该去掉。
end 】
myfun2_2.m (选择其它的函数进行实验，可修改此程序)
【 function y=myfun2_2(x)
sin(x); %y=sin(x)./(x); 】
4.3 实验与观察(Ⅱ)： Matlab符号微积分简介
4.3.1 创建符号变量
4.3.2 求符号极限limit指令
◆观察 ：求下列问题的极限：
【 syms x a
I1=limit('(sin(x)-sin(3*x))/sin(x)',x,0)
I2=limit('(tan(x)-tan(a))/(x-a)',x,a)
I3=limit('(3*x-5)/(x^3*sin(1/x^2))',x,inf) 】
4.3.3 求导指令diff
1.符号求导指令diff
◆
【 syms x y
f=sym('exp(-2*x)*cos(3*x^(1/2))')
diff(f,x)
g=sym('g(x,y)') %建立抽象函数。
f=sym('f(x,y,g(x,y))') %建立复合抽象函数。
diff(f,x)
diff(f,x,2) 】
2.数值求导指令diff
【 x=linspace(0,2*pi,50);y=sin(x);dydx=diff(y)./diff(x);
plot(x(1:49),dydx),grid 】
4.3.4 求符号积分int
【 syms x y z
I1=int(sin(x*y+z),z)
I2=int(1/(3+2*x+x^2),x,0,1)
I3=int(1/(3+2*x+x^2),x,-inf,inf) 】
4.3.5 化简、提取和代入
【 syms x a
t=(a+x)^6-(a-x)^6
t_expand=expand(t)
t_factor=factor(t_expand)
t_simplify=simplify(t) 】
◆ 观察： 将
代入表达式
中求值。
【 syms a b c x y a0 y0
f=a*b+c/x*y;
a='a0';b=1;c=4;x='x0';y=5;
t= subs(f) 】
【 syms a b c x y a0 y0
f=a*b+c/x*y;
subs(f,{a,b,c,x,y},{'a0',1,4,'x0',5}) 】
4.4应用、思考和练习
4.4.1 追击问题
1.追击问题的数值模拟
2. 追踪雷达失效的情形
3. 追线问题的解析解
【 syms y a v s0
x=sym('x(y)'); t=sym('t(y)'); %定义函数关系。
f_left=-y*diff(x,y); f_right=s0+a*t-x; %方程左、右边表达式。
r_left=diff(f_left,y), r_right=diff(f_right,y) %求导 】
【 syms y d r
xs=simplify(dsolve('D2x=-r*sqrt(1+Dx^2)/y','x(20)=0','Dx(20)=0','y')) 】
【 r=0.4; y=20:-0.01:0;
x=1/2*(y.^(1+r)*r*20^(-r)+y.^(1-r)*r*20^r-40*r-y.^(1+r)*20^(-r)+y.^(1-r)*20^r)/(-1+r^2);
shg,comet(x,y) 】
4.4.2 应用问题
1.枪支的设计
【clear,clf
x=[0.0127,0.0254,0.0381,0.0508,0.0635,0.0762,0.0889,0.1016,0.1143,0.1270,0.1397,0.1524,0.3175,0.3302,0.3429,0.3556,0.3683,0.3810,0.3937,0.4064,0.4191,0.4318,0.4445,0.4572,0.1651,0.1778,0.1905,0.2032,0.2159,0.2286,0.2413,0.2540,0.2667,0.2794,0.2921,0.3048,0.4699,0.4826,0.4953,0.5080,0.5207,0.5334,0.5461,0.5588,0.5715,0.5842,0.5969,0.6096];
p=[0.10135,0.20064,0.27303,0.31095,0.33094,0.33991,0.34474,0.33577,0.31508,0.29578,0.27717,0.26131,0.11859,0.11238,0.10687,0.10204,0.09215,0.09308,0.08894,0.08480,0.08067,0.07722,0.07377,0.07032,0.24545,0.23097,0.21718,0.20339,0.19167,0.17995,0.16823,0.15789,0.14824,0.13927,0.13238,0.12548,0.06757,0.06481,0.06205,0.05929,0.05654,0.05378,0.05102,0.04826,0.04550,0.04274,0.04067,0.03861];
[xx,k]=sort(x);
pp=p(k);plot([0,xx],[0,pp],'-'),grid on,
xlabel('枪 管 长 度 x'),
ylabel('压 强 p') 】
2.天然气井的开采量
实验５.用导数作定性分析
5.1知识要点：函数作图 —用导数定性描述函数
【 clf,x=linspace(-8,8,30);f=(x-3).^2./(4*(x-1)); plot(x,f) 】
【 fplot('(x-3)^2/(4*(x-1))',[-8,8])) 】
【 clf,x=sym('x'); f=(x-3)^2/(4*(x-1)); ezplot(f,[-8,8]) ,
title('(x-3)^2/(4*(x-1))','fontsize',11) ,
xlabel('x','fontsize',11) 】
◆按函数绘图步骤绘制完整的函数图，直接用Matlab符号演算完成必须的计算。
【 df_dx=diff('(x-3)^2/(4*(x-1))') 】
【 sym('x');factor(df_dx) 】
【 f=inline('(x-3)^2/(4*(x-1))');X1=[-1,f(-1)],X2=[3,f(3)] 】
求符号二阶导数d2y/dx2：
【 df2_dx2=diff(‘(x-3)^2/(4*(x-1))’,2) 】
二阶导数的因式分解：
【 sym('x'); factor(df2_dx2) 】
【 syms x
f_left=limit('(x-3)^2/(4*(x-1))',x,1,'left')
f_right=limit('(x-3)^2/(4*(x-1))',x,1,'right') 】
【 syms x
f_minus_inf=limit('(x-3)^2/(4*(x-1))',x,-inf,'right')
f_plus_inf=limit('(x-3)^2/(4*(x-1))',x,inf,'left') 】
③有无斜渐近线？
【 syms x, a=limit('((x-3)^2/(4*(x-1)))/x',x,inf) 】
【 b=limit('(x-3)^2/(4*(x-1))-(1/4)*x',x,inf,'left') 】
5.2实验与观察：微分方程的定性解图示
5.2.1人口增长的预测
1.Malthus模型
2.Logistic模型
【 N=dsolve('DN=r*(1-N/Nm)*N','N(t0)=N0') 】
3.微分方程解的定性分析
观察1：
◆（1）求N = N (t)的驻点和拐点。
【 syms t
dN2_dt2=diff('r*(1-N(t)/Nm)*N(t)',t) , dN2_dt2=factor(dN2_dt2) 】
4.用导数作稳定性分析
下面是绘制图5.8的参考程序。
【 clf, N=linspace(0,300,50);
dN=0.3134*(1-N/250).*N;
plot(N,dN),hold on,
plot([0 300],[0,0]),
plot([0,250/2,250],[0,0,0],'o'),
xlabel('N','fontsize',11),ylabel('dN','fontsize',11),
text(N(32),dN(32),'\leftarrow\it{d N / d t}>0,相点递增右移','fontsize',11),
text(125,dN(45),'\it{dN/dt}<0,相点递减左移 \rightarrow','fontsize',11);
h=text(251/2,1.5,'\it{N_m/2}');set(h,'fontsize',11) 】
5.观察程序及其说明
zxy5_1.m (绘制函数图象，图5.3)
【 clf, x=sym('x'); f=(x-3)^2/(4*(x-1)); g=x/4-5/4;
hold on,
h=line([-8 8],[0,0]); set(h,'color’,'red’);
h=line([0 0],[-8,8]); set(h,'color’,'red’);
line([1 1],[-8 8]);plot([-1 1 3],[-2,0,0],'o’),
ezplot(g,[-8 8]); ezplot(f,[-8,8]), %符号函数绘图
text(-1-0.5,-2-0.5,'(-1,-2)’);text(1,0-0.5,'(1,0)'); text(3,0.5,'(3,0)');
x=1.4;text(x,subs(f),'\leftarrow{(x-3)}^{2}/{4(x-1)}');
x=0.6;text(x,subs(f),'\leftarrow{(x-3)}^{2}/{4(x-1)}');
x=2.5;text(x,subs(g),'\leftarrow斜渐近线{y=x/4-5}');
text(1,-2,'\leftarrow垂直渐近线x=1');title('(x-3)^2/4(x-1)') 】
zxy5_2.m （人口预测，图5.5）
【 global p1;clf,
t1=[1790 1800 1810 1820 1830 1840 1850 1860 1870 1880 1890 1900 1910 1920 1930 1940 1950 1960 1970 1980 ];
x=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.0 72.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5];
p1(1)=3.9; p1(2)=250;p1(3)=1790;p1(4)=0.03134;
[t,N]=ode23('Logistic_fun',[1790 2100],3.9);
plot(t,N,t1,x,'o',t,250*ones(1,length(t))),axis([1790 2100 0 300]),
xlabel('t','fontsize',11), ylabel('N','fontsize',11) 】
Logistic_fun.m
【 function dN = Logistic_fun(t,N)
global p1
N0=p1(1);Nm=p1(2);t0=p1(3);r=p1(4);
dN =r*(1-N/Nm)*N; 】
zxy5_3.m （方程的解轨线和相轨线，图5.6）
【 clear,clf
global p1;
p1(1)=3.9; p1(2)=250;p1(3)=0;p1(4)=0.03134;Nm=p1(2);
tpas=linspace(0,300,1000);
plot([0 0],[0,500],':',[0 300],[Nm,Nm],':',[0 300],[Nm/2,Nm/2],':'),
axis([-50 300 0 500]),xlabel('t'),ylabel('N'),hold on
text(-30,Nm,'\it{N_m}\rightarrow'); text(-35,Nm/2,'\it{N_m/2}\rightarrow');
button=1;
while button==1
k=[];
[t0,N0,button]=ginput(1);
[t,N]=ode23('Logistic_fun',tpas,N0);
k=find(N<=Nm/2+1&N>=Nm/2-1);
ts=tpas(k);Ns=N(k); text(-35,N0,'\it{N_0}\rightarrow');
plot(t0,N0,'o',ts,Ns,'square',t,N,':'),hold on
comet(t,N),pause,comet(0*ones(length(t),1),N)
end 】
这一程序是不难读懂的。
5.3 应用、思考和练习
5.3.1.函数作图
◆（2） 下面的绘图较复杂一些，是一个很好的练习。
zxy5_4.m
【 clf,n=2000;a=-4;b=6;c=-8;d=8;
t=linspace(a,b,n);
x=(t.^2)./(t-1);y=t./(t.^2-1);
kx=find(abs(x)>=d);x(kx)=NaN;
ky=find(abs(y)>=d);x(ky)=NaN;
plot(t,[x;y],'.','markersize',3),
hold on,plot([a b],[0,0],'r',[0 0],[c,d],':'),axis([a b c d]),
xlabel('t'),ylabel('x and y') ,
text(-3.8,7,'put any key to show x=x(t)');pause,comet(t,x),
text(-3.8,6,'put any key to show y=y(t)');pause,comet(t,y) 】
5.3.2.平衡点的分类
5.3.3定性分析的应用
1.捕鱼业持续的收获
画定性分析图的程序zxy5_4.m
【 clf,clear,N=50; r=4.4;E=[0.5 2.2 6.5];x=linspace(0,N,30);
f1=r*x.*(1-x/N);plot(x,r.*x,':','linewidth',2),axis([0 50 0 80]),hold on
text(x(10),r*x(10),['\leftarrow y=rx, r= ',num2str(r)])
for i=1:3
f2(i,:)=E(i)*x;
text(x(5),f2(i,5),['\leftarrow y=',num2str(E(i)),'x'])
end
plot(x,f1,x,f2),hold on,
text(x(22),f1(22),['\leftarrow dx/dt=rx(1-x/N)']),xlabel('x'),ylabel('dx/dt') 】
2. 蚜虫生长和跃变
实验６．最优化实验
6.1知识要点与背景
6.1.1 由简入繁： 最佳水槽断面问题的推广
6.1.2 微分法求最大和最小
◆问题4 （1）求受检点：
，
zxy6_1.m
【 syms x1 x2 %定义符号变量。
f=x1^3-x2^3+3*x1^2+3*x2^2-9*x1; % 函数z。
v=[x1 x2];df=jacobian(f,v) %计算雅可比
。
[X,Y]=solve(df(1),df(2));[X,Y] % 用指令solve求驻点。 】
zxy6_2.m
【 clf,xmin=-5;xmax=3.5;ymin=-3;ymax=5;
x1=linspace(xmin,xmax,30);x2=linspace(ymin,ymax,30);
[X1,X2]=meshgrid(x1,x2);Z= X1.^3.-X2.^3+3*X1.^2+3*X2.^2-9*X1;
contour(X1,X2,Z,60);,hold on, xp=[-3,1,-3,1];yp=[0 0 2 2];
plot(xp,yp,'ro'),axis([xmin xmax ymin ymax]),colorbar
xlabel('x_1'),ylabel('x_2'),
for i=1:length(xp)
text(xp(i),yp(i),['\leftarrow (',num2str(xp(i)),',',num2str(yp(i)),')'] )
end 】
6.2 实验与观察(Ⅰ)：模拟盲人下山的迭代寻优法
zxy6_3.m（盲人下山的模拟）
【 clf, a=-2;b=4;
xmin=a;xmax=b;ymin=a;ymax=b; %设置变量范围和坐标轴显示范围。
x1=linspace(xmin,xmax,100);x2=linspace(ymin,ymax,100);
[X1,X2]=meshgrid(x1,x2);
[Z,DZ1,DZ2]=zxy6_3f(X1,X2); %计算函数和梯度向量。
contour(X1,X2,Z,30), %画等值线图。
axis([xmin xmax ymin ymax]),hold on,
axis equal, %该命令将使横轴、纵轴具有相同比例，避免失真。
plot([1.46808510638298],[1.148936170212776],'o'), %标注最优点。
axis([xmin xmax ymin ymax])
x=[];y=[]; %开始用鼠标选点，按左键选点，按右键中止选点过程。
disp('Select a point by put on mouse left-key')
%disp指令,在命令窗口显示文字。
disp('Stop selecting point by put on mouse right-key')
button=1; %button和ginput命令结合使用可用鼠标选点， 按左键时button＝1。
x=[];y=[];
while button==1
[xi,yi,button]=ginput(1);
%ginput(n)用鼠标选n个点，xi,yi分别为点的横坐标和纵坐标。
plot([xi],[yi],'r.','MarkerSize',10)，hold on, %画所选的点。
[zi,dz1,dz2]=zxy6_3f(xi,yi); %计算函数值和梯度方向。
v=zi;
contour(X1,X2,Z,[v v],'-'), %在点所在的高度画一条等高线。
axis([xmin xmax ymin ymax]),
x=[x,xi];y=[y,yi];
H_line2=plot(x,y); %画已走的路径连线。
set(H_line2,'color','red','linewidth',2); %设置颜色和线宽。
xt=xi-dz1;yt=yi-dz2;
H_line=plot([xi xt],[yi yt],'k:','linewidth',1); %画最速下降方向路径。
end %若按左键button=1，继续循环。若按右键，button~=1,循环终止 。 】
zxy6_3f.m（模拟山谷的二次函数程序）
【 function [f,df1,df2]=zxy6_4f(x1,x2)
f=8*x1.*x1+9*x2.*x2-10*x1.*x2-12*x1-6*x2; %计算函数值。
if nargout > 1
df1=2*8*x1-10*x2-12*ones(size(x1)); %计算梯度向量。
df2=2*9*x2-10*x1-6*ones(size(x2));
end 】
6.3 .实验与观察（Ⅱ）：Matlab优化工具箱简介
6.3.1多元函数无约束优化指令fminunc和fminsearch
1. 观察：运行香蕉函数的优化程序bandemo.m
2. 使用fminunc和fminsearch指令
◆ 观察：用inline生成函数。
【 f=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2'),
x=[2,2],y=f(x), %代入一个点计算看看效果 。 】
3. bandemo.m的简化和剖析
zxy6_4.m
【 clf; clear
%以下程序段是画香蕉函数图形。
xx = [-2:0.125:2]'; yy = [-1:0.125:3]'; [x,y]=meshgrid(xx',yy') ;
meshd = 100.*(y-x.*x).^2 + (1-x).^2; conts = exp(3:20);
xlabel('x1'),ylabel('x2'),title('Minimization of the Banana function')
contour(xx,yy,meshd,conts), hold on
plot(-1.9,2,'ro'),text(-1.9,2,'Start Point')
plot(1,1,'ro'),text(1,1,'Solution')
%优化程序段开始。
x0=[-1.9,2]; %赋初值。
l=1;
while l %while 语句是可以重复运行下面的程序段，直至l=0退出循环。
clc %清除命令窗口的全体内容 。
%以下程序段是在命令窗口显示相应的文字内容 。
disp(' ')
disp(' Choose any of the following methods to minimize the …
banana function')
disp('')
disp(' UNCONSTRAINED: 1) BFG direction ')
disp(' 2) DFP direction')
disp(' 3) Steepest Descent direction')
disp(' 4) Simplex Search')
disp(' 0) Quit')
method=input('Select method : '); % input 从键盘输入控制变量method数据。
switch method %Switch体开始。
case 0 %当method＝0，终止程序。
hold off
disp('End of demo')
break %break指令：中止程序。
case 1 %当method＝1，采用BFGS法。
clf,hold on %每一个case中重新画等值线图,下面的程序段是重新画图。
xlabel('x1'),ylabel('x2'),
title('Minimization of the Banana function')
contour(xx,yy,meshd,conts)
plot(-1.9,2,'ro'), text(-1.9,2,'Start Point')
plot(1,1,'ro'), text(1,1,'Solution')
% 这里是学习的重点： OPTIONS是控制fminunc和fminsearch指令的重要参数，
%用optimset('属性'，'属性值',…)指令改变设置，可以容易地控制算法。
OPTIONS=optimset('LargeScale','off');
%fminunc默认的大规模算法是“信赖域方法”，这是一种有效的算法；
%将LargeScale的属性设置为off时,fminunc的默认中等规模的算法就是BFGS方法。
OPTIONS = optimset(OPTIONS,'gradobj','on'); %使用解析梯度。
%定义梯度函数和画图函数banplot6_4。
GRAD=inline('[100*(4*x(1)^3-4*x(1)*x(2))+2*x(1)-2;…
100*(2*x(2)-2*x(1)^2); banplot6_4(x)]');
f=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2'); %定义目标函数。
disp('[x,fval,exitflag,output] = fminunc({f,GRAD},x0,OPTIONS);');
%（调用fminunc指令，输出x,fval分别为最优点和最优函数值，exitflag和output
% 提供算法的一些信息,读者可在程序结束后，键入output或exitflag查看这些信息）
[x,fval,exitflag,output] = fminunc({f,GRAD},x0,OPTIONS);
hold off
disp(' ')
disp('Strike any key for menu')
pause
case 2 %当method＝2，采用DFP法。
clf, xlabel('x1'),ylabel('x2'),
title('Minimization of the Banana function')
contour(xx,yy,meshd,conts), hold on
plot(-1.9,2,'ro'), text(-1.9,2,'Start Point')
plot(1,1,'ro'), text(1,1,'Solution')
OPTIONS=optimset('LargeScale','off');
OPTIONS = optimset(OPTIONS,'gradobj','on');
OPTIONS=optimset(OPTIONS,'HessUpdate','dfp');
% 将HessUpdate属性设置为dfp就使fminunc指令采用DFP法。
GRAD=inline('[100*(4*x(1)^3-4*x(1)*x(2))+2*x(1)-2;…
100*(2*x(2)-2*x(1)^2); banplot6_4(x)]');
f=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2');
disp('[x,fval,exitflag,output] = fminunc({f,GRAD},x0,OPTIONS);');
[x,fval,exitflag,output] = fminunc({f,GRAD},x0,OPTIONS);
hold off
disp(' ')
disp('Strike any key for menu')
pause
case 3 %当method＝3，采用最速下降法。
clf, xlabel('x1'),ylabel('x2'),
title('Minimization of the Banana function')
contour(xx,yy,meshd,conts)
hold on
plot(-1.9,2,'ro'), text(-1.9,2,'Start Point')
plot(1,1,'ro'), text(1,1,'Solution')
OPTIONS=optimset('LargeScale','off');
OPTIONS = optimset(OPTIONS,'gradobj','on');
OPTIONS=optimset(OPTIONS,'HessUpdate','steepdesc');
%将HessUpdate属性设置为steepdesc就使fminunc指令采用最速下降法。
GRAD=inline('[100*(4*x(1)^3-4*x(1)*x(2))+2*x(1)-2;…
100*(2*x(2)-2*x(1)^2); banplot6_4(x)]');
f=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2');
disp('[x,fval,exitflag,output] = fminunc({f,GRAD},x0,OPTIONS);');
[x,fval,exitflag,output] = fminunc({f,GRAD},x0,OPTIONS);
hold off
disp(' ')
disp('Strike any key for menu')
pause
case 4 %当method＝4，采用单纯形方法。
clf,hold on, xlabel('x1'),ylabel('x2'),
title('Minimization of the Banana function')
contour(xx,yy,meshd,conts),
plot(-1.9,2,'ro'), text(-1.9,2,'Start Point')
plot(1,1,'ro'), text(1,1,'Solution')
OPTIONS=optimset('LargeScale','off');
OPTIONS = optimset(OPTIONS,'gradobj','off');
%该方法不使用导数，所以要设置gradobj属性为off。
f=inline('[100*(x(2)-x(1)^2)^2+(1-x(1))^2; banplot6_4(x)]');
%如果要画迭代过程的中间图，就要编制一个画图程序 banplot6_4，
% 套用本程序的格式定义目标函数。
disp('[x,fval,exitflag,output] = fminsearch(f,x0,OPTIONS);');
[x,fval,exitflag,output] = fminsearch(f,x0,OPTIONS);
%fminsearch 是多变量函数寻优的单纯形法指令，用法和fminunc是类似的。
hold off
disp(' ')
disp('Strike any key for menu')
pause
end
end 】
banplot6_4.m
【function out = banplot6_4(x)
plot(x(1),x(2),'O','Erasemode','none')
drawnow; % Draws current graph now
out = []; 】
6.3.2其它的优化算法指令
1.多变量约束优化指令fmincon
2. 线性规划linprog指令
【 f = [-5; -4; -6]; A = [1 -1 1;3 2 4;3 2 0];b = [20; 42; 30];
lb = zeros(3,1); x0=[1,1,0];
options=optimset('Diagnostics','on', 'largescale','off');
%查看诊断信息并采用单纯形算法
[x,fval,exitflag,output,lambda]= …linprog(f,A,b,[],[],lb,[],x0,options)
%没有等式约束且变量无上界,故需置Aeq=[ ];Aeb=[ ]；ub=[ ]; 】
3. 二次规划quadprog指令
4. 一元函数寻优fminbnd指令
5.非线性最小二乘指令lsqnonlin和非线性数据拟合指令lsqcurvefit
程序如下（zxy6_5.m）
【 clf;x=1:10; y=2+2*x; %选直线上的10个点。
a0=[0.1,0.4]; y1=exp(x*a0(1))+exp(x*a0(2)); %计算一条曲线。
[a,resnorm,residual] = lsqnonlin('zxy6_5f',a0); % 求最优解初始点a0。
disp('a='), disp(a);
disp('resnorm='), disp(resnorm)
y2=exp(x*a(1))+exp(x*a(1));
plot( －x,y,x,y1,'r:',x,y2,'o-',x,residual,'.-'),grid on
legend('直线','猜测的曲线','解曲线','残差') 】
函数子程序为（zxy6_5f.m）
【 function F = zxy6_5f(a)
x = 1:10;
F = 2 + 2*x-exp(a(1)*x)-exp(a(2)*x); 】
{ Optimization terminated successfully:
Norm of the current step is less than OPTIONS.TolX
a= 0.2578 0.2578
6.4 应用、思考与练习
6.4.1 .计算最佳水槽断面面积
zxy6_6S.m
【 function s=zxy6_6S(x)
l=24;a(1)=x(3);a(2)=x(4); xs0=(0.5*l-x(1)-x(2));
xs1=xs0+x(1)*cos(a(1)); xs2=xs1+x(2)*cos(a(2));
h1=x(1)*cos(a(1)); h2=x(2)*cos(a(2));
s=(xs0+xs1)*h1+(xs1+xs2)*h2;s=-s; 】
zxy6_6.m
【 clf,A=[1,1,0,0;-1,-1,0.0,0.0];b=[12,0]';
lb=[0,0,0,0]';ub=[100,100,pi/2,pi/2]';x0=[4,4,pi/3,pi/3]';
[s,fval] = fmincon('zxy6_6S',x0,A,b,[],[],lb,ub) %最优化计算
%以下是绘制最优断面的图形，首先计算坐标点。将底边放在x轴上，并让断面关于
%y轴对称。逆时针计算坐标点，使之成为一个封闭的图形)
x(1)=(24-2*s(1)-2*s(2))/2; y(1)=0;
x(2)=x(1)+s(1)*cos(s(3)); y(2)=s(1)*sin(s(3));
x(3)=x(2)+s(2)*cos(s(4)); y(3)=y(2)+s(2)*sin(s(4));
x(4)=-x(3); y(4)=y(3);x(5)=-x(2); y(5)=y(2);x(6)=-x(1); y(6)=y(1);
x(7)=x(1);y(7)=y(1); %首尾相接。
%plot(x,y),axis equal %用这命令可画出封闭图形。
patch(x,y,'y'); axis equal, %用patch命令画块对象并填充颜色。 】
{ s = 4.8000 4.8000 0.6283 1.2566
fval = -88.6373 }
6.4.2 对约束优化的讨论
6.4.3.工程优化问题的计算
1. 啤酒配方问题: 线性规划
2. 储能飞轮的设计
3. 齿轮减速器设计
实验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时分别对应函数
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. 混沌对初值的敏感性
for r=[3.4]
x(1)=0.2;
for j=2:30
x(j)=r*x(j-1)*(1-x(j-1));
end
end
clf,plot(x,'ko-')
4. 其它迭代函数
――――――――――――――――――――――――――――――――――――――――
练习与思考
（7）相变――作为转折点的不动点
市场上有一种产品 A，想知道它在市场上占有的份额将如何变化，是增加、减少还是保持不变？为此可以对市场进行调查。
假定产品A在当前市场中占有的份额的百分比是p 。做多次抽样调查，调查的方法是，每次找三个顾客了解他们愿意购买产品A的意向。由于人一般会具有从众心理，所以按多数原则决定每一组的购买意向：如果在一个被调查的三人组中多数人表示愿意购买产品 A，也就是三人组中至少有两个人愿意购买A，则认为该组是愿意购买产品A的。根据这种随机抽样调查，例如调查了500组，有300 组愿意购买产品A。那么根据统计的观点，可以认为
7)供需曲线的均衡
8)相变问题[7]151
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) 】
另外可用solve指令对方程组（7.15）进行符号求解。程序如下
【 syms x1 x2 x3 %用syms 对每个符号变量进行说明
f1='sin(x1)+x2+x3^2*exp(x1)-4'; %象这样定义各个方程
f2='x1+x2*x3';
f3='x1*x2*x3';
[x1,x2,x3]=solve(f1,f2,f3); %用solve指令求解
solution=[x1,x2,x3]' 】
{ solution =
[ 0, 0, 0]
[ 0, 0, 4]
[ 2, -2, 0] }
也就是方程组的解有三组（0，0，2）、（0 0 －2）和（0 4 0）。在数值求解中已经求得了第一组解。下面来验证一下，将上面三组值代入方程计算之
fval1=zxy7_7f([0,0,2])',fval2=zxy7_7f([0,0,-2])',fval3=zxy7_7f([0,4,0])'
结果为
fval1 =
0 0 0
fval2 =
0 0 0
fval3 =
0 0 0
所以是解。
【 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',14)
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做线性代数题（矩阵的符号演算）
◆ 试题 设
，求 a的值使得 1. a4不能由a1、a2、a3线性表示。 2. a4可以由a1、a2、a3线性表示，并写出表达式。
【 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;
%A=[cos(a) -sin(a);sin(a) cos(a)];
%break
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. 随机变量的模拟：逆概率法
【 clf
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,mu,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<l); %求整根报废钢材的下标。
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_chenpin);
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)','fontsize',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'; %代入sx的具体数值。
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;
xx=[];
for i=1:length(cc);
c=cc(i);
[t,x]=ode45('zxy10_2f',tspan,x0); xx=[x(:,1),xx];
end
animinit('one'); %逐条观察振动图形。
for i=1:length(cc)
c=cc(i)
plot3(t,c*ones(length(t),1),xx(:,i),'r:'),hold on,
view(30,60); %适当选择观察角度（请查阅该指令的帮助，了解用法）。
comet3(t,c*ones(length(t),1),xx(:,i)),axis([ 0 4 -0.2 1.2 -1.1 1.1])
end 】
10.2.2振动弹簧的实时动画
zxy10_2.m
【 animinit('onecart1 Animation')
axis([-2 6 -10 10]); hold on; u=2;
xy=[ 0 0 0 0 u u u+1 u+1 u u;
-1.2 0 1.2 0 0 1.2 1.2 -1.2 -1.2 0];
x=xy(1,:);y=xy(2,:);
% Draw the floor under the sliding masses
plot([-10 20],[-1.4 -1.4],'b-','LineWidth',2);
hndl=plot(x,y,'b-','EraseMode','XOR','LineWidth',2);
set(gca,'UserData',hndl);
for t=1:0.025:100;
u=2+exp(-0.00*t)*cos(t);
x=[0 0 0 0 u u u+1 u+1 u u];
hndl=get(gca,'UserData');
set(hndl,'XData',x);
drawnow
end 】
10.3.3物理问题的数值模拟
下面列举两个物理模拟的例子，用Matlab模拟它们是十分有趣的。
1.多普勒效应的模拟
【 x0=500;v=60;y0=30;c=330;w=1000;t=0:0.001:30;
r=sqrt((x0-v*t).^2+y0.^2);t1=t-r/c;
u=sin(w*t)+sin(1.1*w*t);u1=sin(w*t1)+sin(1.1*w*t1);
sound (u);pause (5),sound (u1) 】
2. 用image指令模拟 两点（双缝）光干涉图案
◆观察：读取Matlab中的mpgcover.jpg图形文件并画出。
【 clf,A=imread('mpgcover','jpg');image(A) 】
双缝干涉参考程序
zxy10_6.m
【 Lamda=0.0000006;d=2;z=1000; ymax=5*Lamda*z/d;n=1000;
x=[0,4];y=[-ymax,ymax]; %设定屏幕。
ys=linspace(-ymax, ymax, n);L1=sqrt((ys-d/2).^2+z^2);
L2=sqrt((ys+d/2).^2+z^2);Phi=2*pi*(L2-L1)/Lamda;B=4*cos(Phi/2).^2;
clf;figure(gcf);NCLevels=255;Br=(B/4.0)*NCLevels;Br=Br';
subplot (1,2,1),image(x,y, Br);colormap(gray(NCLevels));
%colormap('default');
%colormap('hot');
Subplot(1,2,2),plot(B,ys) 】
实验11. 傅氏分析与小波分析
11.1 知识要点 — 傅氏分析与小波分析
11.1.1傅氏分析
11.1.2小波分析
11.2 观察与实验
11.2.1 信号频谱分析
xsf11_1.m (信号生成与显示程序)
【 dalt=0.002; %采样间隔
t=0:0.002:1.2;
rn=randn(1,length(t));rn(1:300)=0; %产生随机序列
s=sin(2*pi*10*t)+sin(2*pi*50*t)+rn; %生成模拟信号
save singal1 dalt s;
clear;
load singal1;
t=[0:length(s)-1]*dalt;plot(t,s,'k');Ylabel('幅值');Xlabel('时间');
title('模拟信号'); 】
xsf11_2.m (模拟信号singal1的频谱分析)
【 clear;load singal1; t=[0:length(s)-1]*dalt;
subplot(211);plot(t,s); Ylabel('幅值'); Xlabel('时间');title('原始信号');
fs=fft(s,512); %快速傅氏变换
pp=fs.*conj(fs)/512; %计算功率谱
ff=(0:255)/512/dalt; %计算各点对应的频率值
subplot(212);plot(ff,pp(1:256)); Ylabel('功率谱密度');Xlabel('频率');
title('信号功率谱'); 】
11.2.2 如何得到小波函数
◆用下面的程序产生并显示与小波db4相关联的四个滤波器组。
xsf11_3.m
【 wname='db4';
[LD,HD,LR,HR]=wfilters(wname);
subplot(221);stem(LD);title('低通分解滤波器');grid;
subplot(223);stem(LR);title('低通重构滤波器');grid;
subplot(222);stem(HD);title('高通分解滤波器');grid;
subplot(224);stem(HR);title('高通重构滤波器');grid; 】
运行结果见图11.3。
◆用下面的程序观察小波函数db4的迭代生成过程。
xsf11_4.m
【 for I=1:10
[phi,psi,xval]=wavefun('sym4',I);
plot(xval,psi+2*I);title('小波函数的生成过程');
hold on
end
hold off 】
结果见图11.4。
11.2.3单尺度一维离散小波分解与重构
xsf11_5.m（信号singal1的单尺度离散小波分解与重构）
【 clear;
load singal1; %装载模拟信号
ls=length(s);
t=[0:ls-1]*dalt;
subplot(511);plot(t,s);Ylabel('S'); %显示原始信号
[C,D]=dwt(s,'db4'); %用小波db4对信号进行单尺度分解
subplot(523);plot(C);Ylabel('C'); %显示低频分解系数
subplot(524);plot(D);Ylabel('D'); %显示高频分解系数
SCR=upcoef('a',C,'db4',1,ls); %用低频分解系数重构
SDR=upcoef('d',D,'db4',1,ls); %用高频分解系数重构
subplot(513);plot(t,SCR);Ylabel('SCR'); %显示低频重构信号
subplot(514);plot(t,SDR);Ylabel('SDR'); %显示高频重构信号
SR=idwt(C,D,'db4',ls); %对信号进行完全重构
subplot(515);plot(t,SR);Ylabel('SR'); %显示完全重构后的信号 】
11.2.4多尺度分解与重构
xsf11_
【 clear;
load singal1; %装载模拟信号
ls=length(s);
t=[0:ls-1]*dalt;
subplot(711);plot(t,s);Ylabel('S'); %显示原始信号
[C,L]=wavedec(s,3,'db4'); %用小波db4对信号进行多尺度分解（三层）
C3=appcoef(C,L,'db4',3); %提取最低频分解系数
D3=detcoef(C,L,3); %提取次低频分解系数
D2=detcoef(C,L,2); %提取次高频分解系数
D1=detcoef(C,L,1); %提取最高频分解系数
subplot(789);plot(C3);Ylabel('C');
subplot(7,8,10);plot(D3);
subplot(746);plot(D2);
subplot(724);plot(D1);
SRC3=wrcoef('a',C,L,'db4',3); %用最低频分解系数重构
SRD3=wrcoef('d',C,L,'db4',3); %用次低频分解系数重构
SRD2=wrcoef('d',C,L,'db4',2); %用次高频分解系数重构
SRD1=wrcoef('d',C,L,'db4',1); %用最高频分解系数重构
subplot(713);plot(t,SRC3);Ylabel('SRC3');
subplot(714);plot(t,SRD3);Ylabel('SRD3');
subplot(715);plot(t,SRD2);Ylabel('SRD2');
subplot(716);plot(t,SRD1);Ylabel('SRD1');
SR=waverec(C,L,'db4'); %对信号进行完全重构
subplot(717);plot(t,SR);Ylabel('SR'); 】
11.3 应用、 思考与练习
11.3.1信号的奇异性检测
11.3.2信号去噪
11.3.3信号的压缩
11.3.4练习
实验12．金融分析实验
12.1知识要点和背景：最优投资组合及其计算
12.1.1包含无风险证券的投资组合
12.1.2无风险证券投资组合的计算
mzy1.m
【 ra=0.12;rb=0.08;rf=0.06;rp=0.1;
Aeq=[ra,rb,rf;1 1 1];beq=[rp 1]';lb=[0 0 -100]';
x0=[1 1 1]'/3;
[x,fval]=fmincon('mzy1f',x0,[],[],Aeq,beq,lb,[],[]),
pause;x0=x
[x,fval]=fmincon('mzy1f',x0,[],[],Aeq,beq,lb,[],[]),
x0=[6/19,20/19,-7/19] 】
mzy1f.m
【 function f=mzy1f(x)
v(1,1)=10;v(1,2)=0;v(1,3)=0;v(2,1)=0;v(2,2)=1;v(2,3)=0;v(3,:)=0;
f=x'*v*x; 】
12.2 机会的价值
12.2.1 简单二项式模型机会价值
以上计算方法是从期末开始向期初倒推，称为后向式动态规划。Maltab程序为：
mzy2.m % 一期二项式方法计算Call options
【 clear;
S0=input('输入当前股票价格:');
E=input('输入协定执行价格:');
u=input('输入上升比例u:');
d=input('输入下降比例d:');
r=input('输入无风险利率r:');
if any(S0 <= 0 | E<=0 | r<0 | u<=0 | d<0 | u<1+r | 1+r<d )
error(sprintf('输入 S0, E, u, d 和 r > 0. 输入 u>(1+r)>d.'))
break
end
fu=max(u*S0-E,0)；fd=max(d*S0-E,0)；
p=(1+r-d)/(u-d)；f=(p*fu+(1-p)*fd)/(1+r);
disp('所求的合约价值为:');f 】
12.2.2 两期二项式模型机会价值
mzy3.m (一期二项式方法计算Call options)
【 clear;
SO=input('输入当前股票价格:');
E=input('输入协定执行价格:');
u=input('输入上升比例u:');
d=input('输入下降比例d:');
r=input('输入无风险利率r:');
if any(SO <= 0 | E<=0 | r<0 | u<=0 | d<0 | u<1+r | 1+r<d )
error(sprintf('输入 SO, E, u, d 和 r > 0. 输入 u>(1+r)>d.'))
break
end
fuu=max(u*u*SO-E,0)；fud=max(u*d*SO-E,0)；
fdd=max(d*d*SO-E,0)；p=(1+r-d)/(u-d)；
fu=(p*fuu+(1-p)*fud)/(1+r)；fd=(p*fud+(1-p)*fdd)/(1+r)；
f=(p*fu+(1-p)*fd)/(1+r);
disp('所求的合约价值为:'); f 】