数学实验与Matlab

实 验 一 矩阵运算与Matlab命令

日常矩阵及其运算 矩阵应用实例: 榄球防护用品的生产管理

应用问题 一个工厂生产三种橄榄球用品 ： 防护帽、 垫肩、臀垫。 需要不同数量的：硬塑料 、 泡沫塑料 尼龙线 、 劳动力。 需要不同数量的：硬塑料 、 泡沫塑料 尼龙线 、 劳动力。 为监控生产，管理者对它们之间的关系十分关心。 为把握这些量的关系，他列出下面的表

原料产品关系表

订单 管理者接到四份订单如上表所示。 问应该如何计算每份订单所需的原材料，以便组织生产？

将表格写成矩阵形式

计 算 输入下面Matlab指令 A=[4 2 3;1 3 2;1 3 3;3 2 2], 计 算 输入下面Matlab指令 A=[4 2 3;1 3 2;1 3 3;3 2 2], B=[35 20 60 45;10 15 50 40;20 12 45 20] C=A*B 请自行计算观看结果 A=[4 2 3;1 3 2;1 3 3;3 2 2], B=[35 20 60 45;10 15 50 40;20 12 45 20] C=A*B

Matlab基本指令 向量的创建和运算

1. 直接输入向量 》x1=[1 2 4],x2=[1,2,1],x3=x1’ 运行结果 x1 = 1 2 4 x2 = 1 2 1

2.冒号创建向量 x1=3.4:6.7, x2=3.4:2:6.7 x3=2.6:-0.8:0 运算结果 x1 = 3.4000 4.4000 5.4000 6.4000 x2 = 3.4000 5.4000 x3 = 2.6000 1.8000 1.0000 0.2000

3.生成线性等分向量 指令x=linspace(a,b,n) 在[a,b]区间产生 n 个等分点(包括端点) 结果 x = 0 0.2500 0.5000 0.7500 1.0000

工作空间 在Matlab窗口创建向量后并运行后，向量就存在于工作空间，可以被调用。 x1=3.4:6.7, x2=3.4:2:6.7, x3=2.6:-0.8:0; x=linspace(0,1,5), whos

向量的运算 设x=[x1 x2 x3]; y=[y1 y2 y3];为两个三维向量，a,b为标量。 向量的数乘：a*x=[a*x1 a*x2 a*x3] 向量的平移： x+b=[x1+b x2+b x3+b] 向量和： x+y=[x1+y1 x2+y2 x3+y3] 向量差： x-y=[x1-y1 x2-y2 x3-y3] 数的乘幂： 如 a^2 x=[1 2 3]; y=[2 2 2]; a=10; a*x, x+a, x+y, x-y, a^2

元素群运算(四则运算） x.*y=[x1*y1 x2*y2 x3*y3] (元素群乘积) x./y=[x1/y1 x2/y2 x3/y3] (元素群右除，右边的y做分母) x.\y=[y1/x1 y2/x2 y3/x3] (元素群左除，左边的x做分母) x.^5=[x1^5 x2^5 x3^5] (元素群乘幂) 2.^x=[2^x1 2^x2 2^x3] (元素群乘幂) x.^y=[x1^y1 x2^y2 x3^y3] (元素群乘幂) x.*y ,x./y,x.\y,x.^5, 2.^x,x.^y

元素群运算（函数计算） Matlab有许多内部函数，可直接作用于向量产生一个同维的函数向量。 x=linspace(0,4*pi,100);（产生100维向量x） y=sin(x); (y也自动为100维向量) y1=sin(x).^2; y2=exp(-x).*sin(x); 观察结果 x=linspace(0,4*pi,100); y=sin(x); y1=sin(x).^2; y2=exp(-x).*sin(x); plot(x,y,'-',x,y1,'-',x,y2,'.-')

创建矩阵（数值矩阵的创建） 直接输入法创建简单矩阵。 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] 观察运行结果 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]

创建矩阵（符号矩阵的创建） 用指令“syms”说明符号变量。 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] 运行 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]

矩阵的运算(矩阵的加减、数乘、乘积) C=A1+B1 D=A1-B1 syms c, cA=c*A1 A2=A1(:,1:3), B1 G=A2*B1 C=A1+B1 D=A1-B1 syms c, cA=c*A1 A2=A1(:,1:3), B1 G=A2*B1

矩阵的运算(矩阵的加减、数乘、乘积) 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, 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 ----------------------

矩阵的运算(左除和右除) 左除“ \ ”： 求矩阵方程AX=B的解；（ A 、B的行要保持一致） 右除“ / ”： 解为 X=A\B； 当A为方阵且可逆时有X=A\B=inv(A)*B； 右除“ / ”： 求矩阵方程XA=B的解 （A 、B的列要保持一致） 解为 X=B/A ， 当A为方阵且可逆时有X=B/A=B*inv(A)

矩阵的运算(左除和右除) 求矩阵方程： 解：有(A-2I)B＝A 程序 ： 观察结果： 设A、B满足关系式：AB＝2B+A,求B。 B=inv(A-2*eye(3))*A, B=(A-2*eye(3))\A 观察结果： A=[3 0 1; 1 1 0;0 1 4]; B=inv(A-2*eye(3))*A, B=(A-2*eye(3))\A

分块矩阵（矩阵的标识） 1.矩阵元素的标识 ： A(i,j)表示矩阵A 的第 i 行 j 列的元素； 2.向量标识方式 A(vr,vc)： vr=[i1,i2,…,ik]、vc=[j1,j2,…,ju]分别是含有矩阵A的行号和列号的单调向量。 A(vr,vc)是取出矩阵A的第i1,i2,…,ik行与j1,j2,…,ju列交叉处的元素所构成新矩阵。

分块矩阵（矩阵的标识） 取出A的1、3行和1、3列的交叉处元素 构成新矩阵A1 程序 观察结果 3 0 5 1 0;2 3 1 2 1], vr=[1, 3];vc=[1, 3]; A1=A(vr, vc) 观察结果 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分为四块，并把它们赋值到矩阵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] 结果 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，将它的第1行变为0。 程序： 删除上面矩阵A的第1、3行。 结果 3 0 5 1 0;2 3 1 2 1] A(1,:)=[0 0 0 0 0]; A 删除上面矩阵A的第1、3行。 A([1,3],:)=[ ] 结果 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 --------------------------------------- A([1,3],:)=[ ]

生成特殊矩阵 全1阵 全零阵： 单位阵： 随机阵： ones(n), ones(m,n), ones(size(A)) zeros(n) ，zeros(m,n), zeros(size(A)) 常常用于对某个矩阵或向量赋0初值 单位阵： eye(n)，eye(m,n) 随机阵： rand(m,n)， rand(n)=rand(n,n)用于随机模拟，常和rand('seed',k)配合使用。

生成特殊矩阵 将rand指令运行多次，观察结果。 程序： 结果 y1=rand(1,5), y2=rand(1,5), rand('seed',3), x1=rand(1,5), rand('seed',3), x2=rand(1,5) 结果 y1=rand(1,5), y2=rand(1,5), rand('seed',3), x1=rand(1,5), rand('seed',3), x2=rand(1,5)

常用矩阵函数 det(A) ： 方阵的行列式； rank(A)： 矩阵的秩； eig(A)： 方阵的特征值和特征向量； trace(A)： 矩阵的迹； rref(A)： 初等变换阶梯化矩阵A svd(A)： 矩阵奇异值分解。 cond(A)： 矩阵的条件数；

数据的简单分析 1.当数据为行向量或列向量时，函数对整个向量进行计算. 2.当数据为矩阵时，命令对列进行计算，即把每一列数据当成同一变量的不同观察值。 max(求最大)、min(求最小)、mean(求平均值)、sum(求和)、std(求标准差)、cumsum(求累积和)、median(求中值)、diff(差分)、sort(升序排列)、sortrows(行升序排列)等等。

数据的简单分析 观察：生成一个3×6的随机数矩阵，并 将其各列排序、求各列的最大值与各列 元素之和。 程序 结果 rand('seed',1);A=rand(3,6), Asort=sort(A), Amax=max(A), Asum=sum(A) 结果 rand('seed',1);A=rand(3,6), Asort=sort(A), Amax=max(A), Asum=sum(A)

实验二 函数可视化与Matlab作图

函数的可视化 令 f (x), g (x)是周期函数吗？观察它们的图象。 程序 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) 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)

绘制平面曲线(plot指令) plot(x,y)： plot(y)： 以x为横坐标、y为纵坐标绘制二维图形 x,y是同维数的向量； 相当于x=[1,2,…,length(y)]时情形。

绘制平面曲线（绘制多个图形) 1. plot(x,[y1;y2;…]), 2. plot(x,y1), hold on, plot(x,y2), hold off 3. plot(x,y1,x,y2,…) 4.plotyy 两个坐标系，用于绘制不同尺度的函数。

绘制平面曲线（线型、点形和颜色的控制） plot(x,y,‘颜色＋线型＋点形’) plot(x,y,‘颜色＋线型＋点形’,x,y,‘颜 色＋线型＋点形’,… ) 句柄图形和set命令改变属性值，可套用： h=plot(x,y), set(h,‘属性’,属性值,‘属性’,属性值,…) 也可用plot(x,y,'属性','属性值')设置图形对象的属性。

绘制平面曲线（属性变量和属性值） 线宽：LineWidth 点的大小： MarkerSize 线型：LineStyle 颜色：color

绘制平面曲线（例） 观察： 程序 观察结果 改变绘图的线型和颜色。 用grid on 指令为图形窗口加上 网格线，并改变网格的线 型和字体的大小。 程序 h=plot([0:0.1:2*pi],sin([0:0.1:2*pi])); set(h,'LineWidth',5,'color','red'); grid on set(gca,'GridLineStyle','-','fontsize',16) 观察结果 h=plot([0:0.1:2*pi],sin([0:0.1:2*pi])); -------------------------------------- set(h,'LineWidth',5,'color','red'); ------------------------------------ grid on set(gca,'GridLineStyle','-','fontsize',16)

绘制平面曲线（坐标轴的控制） axis指令 axis([xmin xmax ymin ymax])： axis([xmin xmax ymin ymax zmin ymax])： 设定三维图形的坐标范围 ； 其中xmin<x<xmax， ymin<y<ymax ，zmin<z<zmax。

绘制平面曲线（gca属性控制） 改变当前轴对象句柄gca属性 用set(gca,‘属性’,属性值,…)可改变字体大小、坐标刻度等轴对象的内容。例如： set(gca,'ytick',[-1 -0.5 0 0.5 1]) 将 y 坐标按向量[-1 -0.5 0 0.5 1]将刻度分成4格； set(gca,'yticklabel','a|b|c|d|e') 改变y坐标刻度的说明。

绘制平面曲线（gca属性控制,例） 设置y坐标的刻度并加以说明，并改变字体的大小。 程序 运行结果 plot([0:0.1:2*pi],sin([0:0.1:2*pi]),'k.-',);grid on, axis([0 6.3 -1.1 1.1])， set(gca,'ytick',[-1 -0.5 0 0.5 1]), set(gca,'yticklabel','a|b|c|d|e'), set(gca,'fontsize',20) get(gca) 运行结果 plot([0:0.1:2*pi],sin([0:0.1:2*pi]),'k.-','linwidth',3);grid on,axis([0 6.3 -1.1 1.1]) ---------------------------------------- set(gca,'ytick',[-1 -0.5 0 0.5 1]), ------------------------------------- set(gca,'yticklabel','a|b|c|d|e'), ----------------------------------- set(gca,'fontsize',20) ------------------------------------ get(gca)

绘制平面曲线（文字标注） xlabel（‘x轴名称’）；ylabel（‘y轴名称’）；zlabel（‘z轴名称’）； title(‘图形标题’)； xlabel（‘x轴名称’）；ylabel（‘y轴名称’）；zlabel（‘z轴名称’）； text(‘说明文字’)：创建说明文字； gtext（'说明文字'）：用鼠标在特定位置输入文字。 文字标注常用符号： \pi （π）；\alpha （α）；\beta （β）； \leftarrow （左箭头） \rightarrow （右箭头）； \bullet （点号）

绘制平面曲线（程序讲解，exp2_1.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}'); % 'HorizontalAlignment','right' 表示箭头所指的曲线对象在 文字的右边。

绘制平面曲线（程序讲解，exp2_1.m） text(6,sin(6),'\fontsize{15}The Value \it{sin(t)} at {\itt}=6\rightarrow\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 ']); legend('\itsin(t)','{\itAe}^{-\alphat}') 注1： num2str： ['string1' ,num2str,'string2']，用方括号 注2： legend 请结合图形观察此命令的使用

图形窗口的创建和分割 subplot(m,n,k)命令。 例：将一个图形分为9个子图，在第k个子图画 sin(kx) 的图象. 程序： 在图形区域中显示多个图形窗口。 m为上下分割数，n为左右分割数，k为第k子图编号。 例：将一个图形分为9个子图，在第k个子图画 sin(kx) 的图象. 程序： 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 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

若干有用的指令 clf：清除图形窗口已有的内容. shg：显示图形窗口。 clear、 clear x：清除工作空间的已有变量。 figure(n): 打开第n个图形窗口 help: … : 续行号

绘制二元函数 基本步骤： 1.生成二维网格点 2. 计算函数在网格点上的值 3. 绘制函数图形

三维绘图（ 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,'+') 观察结果 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,'+') －－－－－－－－－－－－－－ whos

三维绘图（计算函数值，定义域裁减） 程序： for i=1:n for j=1:n if (1-X(i,j))<eps1|X(i,j)-Y(i,j)<eps1 z(i,j)=NaN; else z(i,j)=1000*sqrt(1-X(i,j))^-1.*log(X(i,j)-Y(i,j)); end 运行exp2_1.m

三维绘图（绘图指令） mesh(X,Y,z) ： meshz(X,Y,z)： meshc(X,Y,z)： surf(X,Y,z)： 除了具有mesh的功能外，还在曲面的下方画出函数z=f(x,y)的等值线图， surf(X,Y,z)： 也是三维绘图指令，与mesh的区别在于mesh绘出彩色的线，surf绘出彩色的面， 运行exp2_1,观察效果

三维绘图（等值线指令） colorbar： 表现二维函数的图形的另一种方式是绘制等 值线图。 contour(X,Y,z,n)： n条等高线，n可缺省； contourf(X,Y,z,n)： 等值线间用不同的颜色填满，有更好的视觉效果； contour3(X,Y,z,n)： 在三维空间画出等值线图 colorbar： 将颜色与函数值对应起来显示在图中。

三维绘图（等值线指令，继续exp2_2显示效果） clf,contour(X,Y,z,40),colorbar contourf(X,Y,z,40),colorbar contour3(X,Y,z,40),colormap([0,0,0]) 为等值线标上函数值: 可套用下面程序的格式. [cs,h]=contour(X,Y,z,15); clabel(cs,h,'labelspacing',244) labelspace是数值标记之间相隔的宽度，默认值为144, 这里取了244， clf,contour(X,Y,z,40),colorbar ---------------------------------- contourf(X,Y,z,40),colorbar contour3(X,Y,z,40),colormap([0,0,0]) [cs,h]=contour(X,Y,z,15); clabel(cs,h, 'labelspacing',244)

空间曲线和运动方向的表现 一条空间曲线可以用矢量函数表示为 它的速度矢量表现为曲线的切矢量:

空间曲线和运动方向的表现 很显然飞行曲线方程为：

绘制空间曲线（指令） plot3(x,y,z)： quiver(X,Y,u,v)：绘制二维矢量， quiver3（X,Y,Z,u,v,w）： 在坐标矩阵点[X,Y]处绘制矢量[u,v], 其中u为矢量的x坐标，v为矢量的y 坐标，其维数不小于2。 quiver3（X,Y,Z,u,v,w）： 绘制三维矢量，用法与quiver类似。 Gradient： [Fx,Fy,Fz]=gradient(F)为函数F数值梯度

绘制空间曲线（程序讲解exp2_3） 运行程序 exp2_3.m clf,t=linspace(0,1.5,20); x=t.^2;y=(2/3)*t.^3;z=(6/4)*t.^4-(1/3)*t.^3; plot3(x,y,z,'r.-,'linewidth',1,'markersize',10),hold on Vx=gradient(x);Vy=gradient(y);Vz=gradient(z); h=quiver3(x,y,z,Vx,Vy,Vz),set(h,'linewidth',1),grid on axis([0 1.5 0 1.5 0 40]) xlabel('x'),ylabel('y'),zlabel('z')，box on 运行程序 clf, t=linspace(0,1.5,20);%: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.-','linewidth',1,'markersize',10),hold on %Vx=gradient(x);Vy=gradient(y);Vz=gradient(z); h=quiver3(x,y,z,Vx,Vy,Vz),set(h,'linewidth',1),grid on axis([0 1.5 0 1.5 0 40]) xlabel('x'),ylabel('y'),zlabel('z') box on

应用、思考和练习 用平行截面法讨论由曲面z=x^2-y^2构成的马鞍面形状。 对于二重积分，积分指示线方法是很有用的，当然你首先得了解一下什么是积分指示线法，请查阅高等数学相关的内容，然后设计一个数学实验，然后用Matlab的绘图工具表现这一方法。

应用、思考和练习 绘制微分方程 y/dx=xy, y(0)=0.4的斜率场， 并将解曲线画在图中，观察斜率场和解曲线的关系。

应用、思考和练习 地球表面的气温差异很大，而且随时间变化，要绘制气温分布图绝不是一件容易的事情。但是赤道温度最高，而在两级最冷，在中间地带则是过渡带。所以可粗略将这种气温分布情况用图形表现出来，试绘制地球表面的气温分布示意图。

应用、思考和练习

应用、思考和练习(若干特殊图形 ) 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) 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)

应用、思考和练习(周期函数的推广）

应用、思考和练习(线性P周期函数）

应用、思考和练习(线性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)

应用、思考和练习(线性P周期函数） p周期函数是一种特殊的线性p周期函数， 为M=0的情形。 如何证明这一点？

应用、思考和练习（循环比赛的名次） 有若干支球队参加循环比赛，他们两两交锋，每场比赛只计胜负，不允许平局，循环赛结束后要根据他们的比赛成绩排列名次。 一种方法是计算得分，得分是每支球队获胜的场次，根据各队的得分排出名次，决定冠军队。

应用、思考和练习（循环比赛的名次 虽然计算各队的得分很容易，但有时按得分排名的方法并不一定合理。 假定有4支球队, 记做v1~v4。在一次循环赛中，v1得分为2，v2得分为2，v3得分为1，v4得分为1，可以把得分写成4维向量的形式：s = [2 2 1 1]T， 在这种情况下应该如何决定[v1, v2,v3,v4]的名次呢？

实验三 函数式－直接确定型模型

从系统分析的观点理解函数 y = f（x） 函数不是枯燥的数学符号 函数是系统 函数是数学模型 是描述自然现象的有力工具 x：自变量，y：因变量，f： 映射规则 函数不是枯燥的数学符号 函数是系统 函数是数学模型 是描述自然现象的有力工具

黑箱模型和经验函数 白箱：映射规则f 已知； 灰箱：映射规则f 部分已知； 黑箱：映射规则f 未知。 对于黑箱模型，只知道输入输出的数据，

经验函数(机床加工问题） 用程控铣床加工机翼断面的下轮廓线时 每一刀只能沿x方向和y方向走非常小的一步。 表3-1给出了下轮廓线上的部分数据 试完成加工所需的数据，画出曲线.

交通事故的调查（司机有责任吗？） 一辆汽车在拐弯时急刹车，结果冲进路边的沟里． 警察闻讯赶到现场，对汽车留在路上的刹车痕迹进行了测量，确定刹车痕迹近似为一圆弧。 讯问司机时，他说，当车进入弯道时刹车失灵，且进入弯道时汽车时速为40英里／小时。 此速度系该路的速度上限。 通过验车证实该车制动器在事故发生时确实失灵， 但司机所说的车速是否真实呢？

交通事故的调查（司机有责任吗？） 通常，作一条基准线来测量刹车痕迹. （水平）距离x沿基准线测得，（垂直）距离y与x垂直. 表3-6给出了刹车痕迹的测量有关值. 警察测量了路的坡度，发现这段路是平的． 请给出一个使警察可以核对速度的计算办法．

航行区域的警示线 某海域上频繁地有各种吨位的船只经过。 为保证船只的航行安全，有关机构在低潮时对水深进行了测量,表3-8是他们提供的测量数据： 表3-8. 水道水深的测量数据 x 129.0 140.0 103.5 88.0 185.5 195.0 105.5 y 7.5 141.5 23.0 147.0 22.5 137.5 85.5 z 4 8 6 8 6 8 8 x 157.5 107.5 77.0 81.0 162.0 162.0 117.5 y -6.5 -81.0 3.0 56.5 -66.5 84.0 -33.5 z 9 9 8 8 9 4 9

航行区域的警示线 其中（x, y）为测量点，z为（x, y）处的水深(英尺)。 船的吨位可以用其吃水深度来反映， 分为 4英尺、4.5英尺、5英尺和 5.5英尺 4 档。 航运部门要在矩形海域（75，200）×（－50，150）上为不同吨位的航船设置警示标记。 请根据测量的数据描述该海域的地貌，并绘制不同吨位的警示线，供航运部门使用。 水深z是区域坐标（x, y）的函数z= z (x, y）， 测量数据只是它的部分取值。 可绘制函数图象和等值线图，将不同吃水线标记图上

插值与拟合（基本原理和区别） 已知有n +1个节点（xj,yj），j = 0,1,…, n 其中xj互不相同 节点(xj, yj)可看成由某个函数 y= f（x）产生 f 的解析表达式可能十分复杂 或不存在封闭形式, 也可以是未知的

插值与拟合（基本原理和区别） 插值： 插值指令 yi=interp1(x1,y1,xi,'method') 构造一个相对简单的函数 y=g(x) 使g通过全部节点 即使g (xj) = yj，j=0,1,…, n 用g (x)作为函数f ( x )的近似。 插值指令 yi=interp1(x1,y1,xi,'method') 对应于插值函数yi=g(xi) 其中x1,y1为节点向量 method＝四个选项： ‘nearest’ 为近邻插值；‘linear’为线性插值； ‘spline’ 为样条插值； 'cubic'为立方函数插值。

插值与拟合（基本原理和区别） 多项式拟合 拟合指令polyfit、polyval 对给定的数据（xj,yj），j = 0,1,…, n 选取适当阶数的多项式(也可采用其它形式的函数) 例如二次多项式g(x)=ax^2+bx+c 使g(x)尽可能逼近（拟合）这些数据 拟合指令polyfit、polyval 用p=polyfit(x1,y1,m)做 m 次多项式拟合 拟合数据向量为x1,y1 多项式系数为p=[p(1),…,p (m),p (m+1)] 即g(x)=p(1)x^m+…p (m)x+p (m+1) 用y = polyval(p,x)计算在x 处 多项式的值 y

观察插值、拟合的效果 运行观察程序exp3_1.m 选取一个已知函数作为参考，并将这一函数的图象用虚线显示在图中。 观察程序允许用鼠标选取节点 按鼠标左键选点，按右键选最后一个点 观察不同的选点方式对各种插值和拟合效果的影响 clf figure(1) a=-1;b=1;n=100; %g=inline('1/(1+x^2)'); g=inline('x^2-x^4'); %g=inline('cos(x)^10'); xx=linspace(a,b,n); for i=1:n gx(i)=g(xx(i)); end ymin=min(gx)*0.8;ymax=max(gx)*1.2; a1=0.5 subplot(2,2,1),plot(xx,gx,'k'),grid,hold on,axis([a b ymin ymax]),title('近邻插值') subplot(2,2,2),plot(xx,gx,'k'),grid,hold on,axis([a b ymin ymax]),title('线性插值') subplot(2,2,3),plot(xx,gx,'k'),grid,hold on,axis([a b ymin ymax]),title('样条插值') subplot(2,2,4),plot(xx,gx,'k'),grid,hold on,axis([a b ymin ymax]),title('多项式拟合') button=1; x1=[a];y1=[gx(1)]; while button==1 [xi,yi,button]=ginput(1); subplot(2,2,1),h=plot(xi,yi,'ko') subplot(2,2,2),h=plot(xi,yi,'ko') subplot(2,2,3),h=plot(xi,yi,'ko') subplot(2,2,4),h=plot(xi,yi,'ko') x1=[xi,x1];y1=[yi,y1]; x1=[b,x1];y1=[gx(n),y1]; ynearest=interp1(x1,y1,xx,'nearest'); %ynearest=interp1(x1,y1,xx,'cubic'); ylinear=interp1(x1,y1,xx,'linear'); yspline=interp1(x1,y1,xx,'spline'); [p,c]=polyfit(x1,y1,4); ypolyfit=polyval(p,xx); subplot(2,2,1),h=plot(xx,ynearest,'k-');set(h,'linewidth',2) subplot(2,2,2),h=plot(xx,ylinear,'k-');set(h,'linewidth',2); subplot(2,2,3),h=plot(xx,yspline,'k-');set(h,'linewidth',2) subplot(2,2,4),h=plot(xx,ypolyfit,'k-');set(h,'linewidth',2)

程序注解（inline指令） 定义内联函数：inline指令 g=inline('x^2-x^4'); 程序

程序注解（ginput） [x,y,button] = ginput(n) 用鼠标在屏幕选n个点，返回这n个点，存于x,y中。 1为左键、2为中间键、 3为右键。

程序注解（插值拟合） xx=linspace(a,b,n); %定义自变量xx ynearest=interp1(x1,y1,xx,'nearest'); ylinear=interp1(x1,y1,xx,'linear'); yspline=interp1(x1,y1,xx,'spline'); [p,c]=polyfit(x1,y1,4); ypolyfit=polyval(p,xx);

程序注解（插值拟合） 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)

插值拟合效果观察 沿曲线选取３个节点，保持等间隔。当 节点较少时，插值的效果如何？ 加密节点，共８个等距节点，观察插值 的效果，如果去掉中间的一个节点，插 值效果又会如何？ 有意偏离原来的曲线，如果误差较大， 将会怎样呢？

实验四 微分、积分和微分方程

定积分--连续求和

定积分--连续求和

三种方法计算数值积分 （1）定义法，取近似和的极限。 （3）解微分方程计算定积分 高等数学中不是重点内容 但数值积分的各种算法却是基于定义建立的 （2）用不定积分计算定积分。 不定积分是求导的逆运算， 而定积分是连续变量的求和（曲边梯形的面积） 表面上看是两个完全不同的概念， 通过牛顿－莱布尼兹公式联系在一起， （3）解微分方程计算定积分

微积分学基本定理 特别，F(b)-F(a) 就是所需的定积分. 在高等数学中总是期望求出不定积分的封闭解. 但数值积分是更有用的工具。 牛顿－莱布尼兹公式不愧为微积分的“基本定理”。

基本定理的推广（解微分方程计算定积分）

基本定理的推广（解微分方程计算定积分）

解微分方程的 Eular折线法

解微分方程的 Eular折线法 将区间n = 4等分（共有5个分点)；计算分点和 相应的函数值 （x(1),x(2), x(3) x(4) x(5)) （f(1),f(2) ,f(3) ,f(4) ,f(5)) 在第一个子区间[x(1),x(2)]上，画出折线段 y(2)=y(1)+f(1)*(x-x(1))代替解曲线段y(x)， 这里y(1)=y0=0 折线段的起点为[x(1), y(1)]，终点为[x(2),y(2)]. 运行exp4_1.m，观察第二、三、四子区间的情况。 n=4;n1=n+1;x=linspace(0,pi,n1);f=myfun2_2(x);[0:n;x;f] y(2)=y(1)+f(1)*(x(2)-x(1)) P_intial=[x(1),y(1)],P_final=[x(2),y(2)], ```````````````````````````````````````````````` clf, a=0;b=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]) 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); 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) x(i+1)],[s(i+1) y(i+1)],'r--') pause %ÔÝÍ££¬n´óʱӦ¸ÃÈ¥µô¡£

符号微积分 用Matlab符号工具箱 （Symbolic Toolbox） 可以进行符号演算

符号微积分(创建符号变量） 创建多个符号变量； sym var syms var1 var2 … f=sym(‘符号表达式’) 创建单个符号变量； syms var1 var2 … 创建多个符号变量； f=sym(‘符号表达式’) 创建符号表达式，赋予f； equ=sym('equation') 创建符号方程 。

符号微积分(极限） limit(‘表达式’,var,a)：求当var →a，表达式的极限 例：求极限： 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) syms x a I1=limit(‘(sin(x)-sin(3*x))/sin(x)’,x,0) 运行结果

符号微积分(求导） diff(f,‘var’,n) 求 f 对变量var 的n阶导数 缺省n时为求一阶导数 缺省变量'var' 时，默认变量为x 可用来求单变量函数导数 多变量函数的偏导数 还可以求抽象函数的导数

符号微积分(求导） 例：求 syms x y f=sym('exp(-2*x) * cos(3 * x^(1/2))') diff(f,x) 运行 syms x y f=sym('exp(-2*x) * cos(3 * x^(1/2))') diff(f,x)

符号微积分(求导） syms x y 运行 例：求 g=sym('g(x,y)') f=sym('f(x,y,g(x,y))') diff(f,x) diff(f,x,2) 运行 syms x y g=sym('g(x,y)') f=sym('f(x,y,g(x,y))') diff(f,x) diff(f,x,2)

符号微积分(积分） int(f,var)：求函数f的不定积分； int(f, var, 积分下限，积分上限)： syms x y z 例：求不定积分 syms x y z I1=int(sin(x*y+z),z) syms x y z I1=int(sin(x*y+z),z)

符号微积分(积分） syms x y z I2=int(1/(3+2*x+x^2),x,0,1) I3=int(1/(3+2*x+x^2),x,-inf,inf) syms x y z I2=int(1/(3+2*x+x^2),x,0,1) I3=int(1/(3+2*x+x^2),x,-inf,inf)

符号微积分(化简、提取和代入) 符号运算的结果比较繁琐，使用化简指令可对其进行化简。 但是不能指望机器可以完成一切，人的推理往往必须的。 常用的化简指令如下 展开指令：expand(表达式)； 因式分解：factor(表达式) 降幂排列：collect（表达式，var） ； 一般化简：simplify(A)；

符号微积分(化简、提取和代入) 观察： 将展开(a+x)^6-(a-x)^6，然后作因式分解。 t_expand=expand(t) t_factor=factor(t_expand) t_simplify=simplify(t) 观察结果 syms x a t=(a+x)^6-(a-x)^6 t_expand=expand(t) t_factor=factor(t_expand) t_simplify=simplify(t) Taylor(sin(x))

数值微积分（梯形公式和辛普森公式） trapz(x,y)，按梯形公式计算近似积分； 其中步长x=[x0 x1 … xn]和函数值y=[f0 f1 …fn]为同维向量， q = quad('fun',a,b,tol,trace,P1,P2,...) （低阶方法，辛普森自适应递归法求积） q = quad8('fun',a,b,tol,trace,P1,P2,...)（高阶方法，自适应法Cotes求积） 在同样的精度下高阶方法quad8要求的节点较少。 x=linspace(0,pi,50); y=sin(x);T=trapz(x,y) --------------------------------------- I1=quad('sin',0,pi), I2=quad8('sin',0,pi),

［x,y］=ode23('fun',tspan,y0,option) （ 低阶龙格－库塔函数） （高阶龙格－库塔函数）

应用、思考和练习(追击问题) 我缉私雷达发现，距离d处有一走私船正以匀速a沿直线行驶，缉私舰立即以最大速度（匀速v）追赶。 若用雷达进行跟踪，保持船的瞬时速度方向始终指向走私船， 缉私舰的运动轨迹是怎样的？是否能够追上走私船？ 如果能追上，需要用多长时间？

应用、思考和练习(追击问题)

应用、思考和练习(追击问题) xs=simplify(xs1) r = dsolve(‘eq1,eq2,…’, ‘cond1,cond2,…’, ‘v’) 方程的符号解 syms y d r xs1= dsolve('D2x=-r*sqrt(1+Dx^2)/y','x(20)=0','Dx(20)=0','y') xs=simplify(xs1) 运行结果，画彗星图 syms y d r xs1=dsolve('D2x=-r*sqrt(1+Dx^2)/y','x(20)=0','Dx(20)=0','y') xs=simplify(xs1) ---------------------------------- 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)

应用、思考和练习(追击问题) xs=simplify(xs1) r = dsolve(‘eq1,eq2,…’, ‘cond1,cond2,…’, ‘v’) 方程的符号解 syms y d r xs1= dsolve('D2x=-r*sqrt(1+Dx^2)/y','x(20)=0','Dx(20)=0','y') xs=simplify(xs1) 运行结果，画彗星图 syms y d r xs1=dsolve('D2x=-r*sqrt(1+Dx^2)/y','x(20)=0','Dx(20)=0','y') xs=simplify(xs1) ---------------------------------- 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)

应用、思考和练习(追击问题，如果雷达失效) 当缉私舰雷达发现d处有一走私船后，雷 达突然损坏 若假定走私船作匀速直线运动（但不知 方向），且缉私舰艇速度v大于走私船速 度a， 则缉私舰应采用什么样的航行路线，不 管走私船从哪个方向逃跑，都能追捕上 它？

实时动画制作（见实验10） 观察：模拟弹簧振动 讨论最简单的情形，一弹簧系统作横向运动，其位移由u=2+cos(t) 所决定， 仿真弹簧的振动

实时动画制作(初始化、见实验10） x=xy(1,:);y=xy(2,:); 程序讲解 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,:); plot([-10 20],[-1.4 -1.4],'k-','LineWidth',2); hndl=plot(x,y,'k-','EraseMode','XOR','LineWidth',2) 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,:); plot([-10 20],[-1.4 -1.4],'k-','LineWidth',2); hndl=plot(x,y,'k-','EraseMode','XOR','LineWidth',2)

实时动画制作(初始化、见实验10）zxy10-2 hndl=get(gca,'UserData'); set(gca,'UserData',hndl); for t=1:0.025:1000; 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,'YData',y); drawnow end set(gca,'UserData',hndl); for t=1:0.025:1000; 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,'YData',y); drawnow end _________________ 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)

电影动画制作（zxy7_3) moviein、 getframe、movie指令 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次

应用、思考和练习(枪支的设计) 枪支发火后，气体压强随子弹在膛内的运动 而变化。枪管长度x的单位为m。 压强p是 距离x的函数，通过实测得到了的一批数据， 子弹射出枪管时的出口速度是设计者关心的问题，如果一只枪管长0.6096m，其膛孔面积4.56×10-5m2，子弹重量0.956N，试决定这种型号枪支的出口速度； 更一般的，确定出口速度和枪管长度的关系曲线，绘制这一曲线，并作出适当的标记。这样的问题和你在高等数学中处理的积分有什么区别吗？

应用、思考和练习(天然气井的开采量) 东方天然气公司钻了一口新的气井，他们希望研究一下将这口井于供气管路联接的经济性 计算此井的压强随时间的变化曲线，由此得到流量Q与时间t的关系，以此估计此井的总开采量。

实验六 最优化实验

最佳水槽断面问题 （矩形断面） 用宽 l = 24 cm的长方铁板 折成一个断面为矩形 的水槽，问怎样的折 法可使水槽的断面面 积达到最大

最佳水槽断面问题 （梯形断面） 将问题１推广等腰梯形的水槽，问怎样折法可使水槽断面面积达到最大?

最佳水槽断面问题 （对称五边形断面）将铁板折成如图所示的对称五边形，问怎样的折法可使水槽的断面面积达到最大？

最佳水槽断面问题（五边）

最佳水槽断面问题（五边） 运行zxy6_6 [s,fval] = fmincon('zxy6_6S',x0,A,b,[],[],lb,ub) 求解

最佳水槽断面问题（多边和无限边） 优化变量数与最大断面面积的关系 断面形状 优化变量数 最大断面积 cm2 矩形断面 1 72 矩形断面 1 72 梯形断面 2 83.14 对称五边形 4 88.637 将铁板折成对称7边形、9边形，一般为对称2n+1边形 可以期望最大断面面积得到进一步的增加 随之而来是计算代价也随之增加。

最佳水槽断面问题（无限边）

最佳水槽断面问题（无限边）

最佳水槽断面问题（无限边）

微分法求最大和最小（高等数学）

微分法求最大和最小（高等数学）

微分法求最大和最小（高等数学） 运行zxy6_1.m syms x1 x2 %定义符号变量。 运行zxy6_2画等值线图并将点标注在图上 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)) % 用指令solve求驻点。 运行zxy6_2画等值线图并将点标注在图上

微分法求最大和最小（高等数学） jacobian(f,v)： 计算函数的符号梯度、雅可比矩阵 如：若f(v),v=[v1 v2] 则df=[df/dv1 df/dv2] 如：若f(v)=[f1(v) f2(v)] ,v=[v1 v2] 则df=[df1/dv1 df1/dv1 df2/dv1 df2/dv2]

微分法求最大和最小（高等数学） solve指令： solve('eqn1','eqn2',…,'eqnn') 求n个方程eqn1，eqn2，… ，eqnn所构成的方程组的根（符号解）

盲人下山与迭代寻优 一个盲人处于山上的某一点x0 要下到谷底，他应如何做？

盲人下山与迭代寻优

多元函数无约束优化指令fminunc、fminsearch的剖析 Matlab优化工具箱简介 多元函数无约束优化指令fminunc、fminsearch的剖析

Matlab优化工具箱简介 观察： 在命令窗口键入bandemo 选择不同方法观察对香蕉函数的优化结果和过程。

Matlab优化工具箱简介 [x,fval,exitflag,output,grad,hessian]= fminunc(fun,x0,options,P1,…) 其中有些项可以缺省，如 exitflag,output,grad,hessian,options,P1,P 2,… 等等。 x0是初始点； fun是目标函数，可以用inline指令或建立M文件的方法生成目标函数；

Matlab优化工具箱简介 参数的传递： 输出： 使用变量P1,P2,…可在目标函数和主程序之间需要传递某些参数 也可使用全局变量Gobal说明来进行传递。 输出： x,fval,exitflag,output,grad,hessian为输出信息， 最优点、最优函数值、算法结束的状态 （exitflag > 0 算法收敛；=0达到最大步骤而停止；<0算法不收敛）、算法结束后的某些信息（如迭代次数、所使用的优化方法等，可在命令窗口查看output的具体内容）、梯度值和海森矩阵，除x之外均可缺省。

Matlab优化工具箱简介 OPTIONS（控制参数） optimset(‘属性’，‘属性值’,…) OPTIONS是一个数组，有多项内容 修改OPTIONS默认值，如默认属性‘LargeScale’＝‘on’， ( 使用“信赖域算法”)。 如果要使用其它方法，就要修改此项设置。

Matlab优化工具箱简介 内联函数inline f=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2'), x=[2,2],y=f(x), 备注 f=inline('100*(x(2)-x(1)^2)^2+(1-x(1))^2'), x=[2,2],y=f(x),

Matlab优化工具箱简介 用M文件生成目标函数（套用如下格式） myfun.m function [f,g,H] = myfun(x) ％ 关键词function不可省，函数myfun和M文件同名。 f = … % 计算目标函数x 。 if nargout > 1 %如有两个输出量（[目标函数，梯度]）。 g = … % 计算g为函数x点的解析梯度（可省）。 if nargout > 2 %如有三个输出量（[目标函数，梯度，海森阵]）。 H = … % H为函数在x点的海森阵，（可省）。 end

Matlab优化工具箱简介（ zxy6_4 讲解运行） bandemo.m的简化和剖析 程序zxy6_4.m是对bandemo.m的简化 基本结构为： （1） 绘制香蕉函数的等值线图，并将Start Point和Solution标在图形上。 （2）用Switch语句结构，允许程序选用BFGS、DFP、最速下降法和单纯形法等四种优化方法。

Matlab优化工具箱简介 多变量约束优化指令fmincon 上面的命令等价于 [x,fval,exitflag,output,lambda,grad,hessian]= fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options,P1,P2, …) 上面的命令等价于

Matlab优化工具箱简介 线性规划linprog指令 算法选择 ： [x,fval,exitflag,output,lambda]=linprog(f,A,b,Aeq,beq,lb,ub,x0,options) 算法选择 ： options=optimset(‘largescale’,‘off’)，单纯形方法； options=optimset('largescale','on')，内点法(默认)。

Matlab优化工具箱简介 一元函数寻优fminbnd指令 此时 x, x1, x2 是标量，f (x) 为标量函数。 [x,fval,exitflag,output] = fminbnd(fun,x1,x2,options,P1,P2,...) 此时 x, x1, x2 是标量，f (x) 为标量函数。

Matlab优化工具箱简介 Quadprog：解二次规划 lsqnonlin： 解非线性最小二乘 lsqcurvefit：非线性数据拟合 lsqnonneg：非负系数的最小二乘法。 lsqlin： 约束最小二乘

应用思考与练习（计算最佳水槽断面面积） 试推导对称2n+1边形面积的一般公式 选择一系列的 n 值，仿照zxy6_6.m计算它们的最大断面面积，观察计算结果的规律性。 在工程实践中并不能保证每一次计算都能够成功，但是本问题即使不成功，你是否也能洞察结果？

应用思考与练习（盲人约束下山） 对盲人下山问题，引入一个有界约束区域，试用图形表现函数在区域边界上的图象。 可以用等值线或用曲顶柱体曲面显示函数在区域变化的情况。 不过建议单独用二维绘图指令plot画出它们的曲线图，观察函数在边界的极值情况。

应用思考与练习（盲人约束下山） 结合高等数学知识， 如果可能，用 Matlab符号演算指令求出函数在不同约 束下的极值点和最值点（例如可用 Lagrange函数方法解决这些问题）。 你也可以在盲人下山模拟中对有约束的 情况进行讨论，这时盲人应该如何前进 呢？

应用思考与练习（啤酒配方问题） 某啤酒厂希望用原料掺水的办法生产一种复合标准的低成本啤酒。其标准要求为：酒精含量为3.1％；发酵前平均比重在1.034～1.040之间；颜色在8～10EBC单位之间；每升混合物中，蛇麻子脂的含量在20～25mg之间。 请根据相关数据算出最优配方。

应用思考与练习（ 储能飞轮的设计） 下面的表达式用于设计储能用的飞轮。 准则是储藏的能量最大。 用约束条件限定了重量、直径、转速和厚度， 试计算最优解。你能确定算出的解是最优的吗？

应用思考与练习（齿轮减速器设计） 抽去各变量的物理意义，齿轮减速器最优设计模型如下： 这是一个具有7个变量、23个约束的优化问题。试对其进行计算。 你可能会遇到很大的困难，你能想办法解决这些困难吗？

应用思考与练习（齿轮减速器设计）

实验七 方程求根、不动点和迭代

隐函数的存在定理的可视化

隐函数的存在定理的可视化 选择特殊的例子 运行zxy7_1.m， 画出曲面z=F(x,y)、x-y平面的图像和它们的交线。 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-2.m [x,fval,exitflag,output]=fzero(fun,x0,options) zxy7-2.m

蛛网图与不动点迭代 称满足方程 f(x)=x的点x为函数f的不动点 求函数f的不动点。可以从一个初始点x0出发，以格式 xn+1=f(xn)进行迭代； 得到x0,x1,x2,……,xn,….. 如果该数列是收敛的，则

蛛网图与不动点迭代

蛛网图与不动点迭代 运行观察程序zxy7_3， 理解蛛网图的原理

简单和复杂：二次函数迭代和混沌 观察 对二次函数 f(x)=rx(1-x) 进行的迭代，其中0 < r < 1是一个可变参数。 1） 固定若干个不同的的值，观察迭代序列的的极限； 迭代N次，略去前n个迭代值，并将后N – n个迭代值画在r-x坐标系中（zxy7_4) 2）用蛛网图观察三种不同类型的迭代。(zxy7_5) 3）加密r的取值，得到加密Feigenbaum图。 （zxy7_4改变参数)

实验八 线性代数实验

向量组的线性关系

向量组的线性关系（产生相关向量，运行zxy8_1） 产生向量：产生m个n维向量，且各向量的分量均在[a , b]之间。 clear n=3; m=2; a=-10; b=10; rand('seed',32), A = unifrnd(a,b,[n,m]) 组合向量 ：产生 m = 2个组合系数, 将组合得到的新向量并入矩阵 A中： x = unifrnd(-1,1,[1,m]), A(:,3)=x(1)*A(:,1)+x(2)*A(:,2) 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)

向量组的线性关系（产生相关向量，运行zxy8_1) A = -1.7383 -9.1707 0.0256 8.7064 6.6219 -9.0842 4.7096 7.5495 -0.0246 -8.9245 -3.5331 8.3272 -5.7470 3.9105 -0.0038 -0.2352 -6.6197 2.1934

向量组的线性关系（产生相关性的判别）

Gauss消元法（运行rref(A),rrefmovie(A)） rref(A) 将矩阵A做行初等变换阶梯化。 B=rref(A) B = 1.0000 0 -0.0011 -0.5360 0.5851 0.2590 0 1.0000 -0.0026 -0.8478 -0.8330 0.9415 0 0 0 0 0 0 rrefmovie(A)：观察到行初等变换的过程 rref(A) rrefmovie(A)

Gauss消元法(同解方程） Rref(A)

Gauss消元法(解） 向量形式

Gauss消元法(基础解系） Ax = 0 的基础解系可由下式计算 ： X=[-B(1:r,r+1:m+s);eye(l)] 其中r = rank(A)，l = m+s-r为基础解系的个数。 r=2;m=2;s=4; B=[-B(1:2,r+1:m+s);eye(m+s-r)] 问题：如何用Matlab解一般的非线性齐次方程组，如A(:,1:4)X=A(:,7)？ rref(A) r=2;m=2;s=4; B=[-B(1:2,r+1:m+s);eye(m+s-r)]

应用练习与思考(平面四杆机构设计)

应用练习与思考(平面四杆机构设计) 某操纵装置采用四杆铰链机构。已知两连架杆（L1,L3）输入角和输出角满足下表数据所示的对应关系 ，机架长度L4= 50mm，试确定其余三杆的长度。

应用练习与思考(平面四杆机构设计)

应用练习与思考(平面四杆机构设计) 确定四杆的长度，并用Matlab绘图指令用图形 表示你的设计结果。你需要设计一种表现方案， 使人可以很明白的看出你的设计结果是对的。 如果只利用表前两组对应角度的值，设计方案 还是唯一的吗？计算一下结果。同样给出直观 表示。体会到自由变量含义？ 如果表中值为4组对应角度的值，你就遇到超定方程了。它没有精确解，只有近似解。你愿意用Matlab去解它吗？试一试。

应用练习与思考(用Matlab做线性代数题)

应用练习与思考(用Matlab做线性代数题) 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); 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); A

矩阵的相似化简

矩阵的相似化简 A=rand(3,3) [pc,lamda]=eig(A) Bl=A*pc Br=pc*lamda Bl=Br

矩阵的相似化简 选择方阵A，如 二阶方阵 A=[1/5,99/100;1,0]; 选择一个初始点(二维列向量)，按下面的 公式进行迭代: xk+1=Axk 观察这些迭代点位置和趋向

矩阵的相似化简(程序zxy8_2.m迭代部分) Clear,clf, n=40;a=-20*100;b=-a;c=a;d=b;p=0.1; A=[1/5,99/100;1,0];axis([a b c d]),grid,hold on button=1 while button==1 [xi,yi,button]=ginput(1);plot(xi,yi,'ko'), hold on,X0=[xi;yi];X=X0; for i=1:n X=[A*X,X0]; h=plot(X(1,1),X(2,1),'k.',X(1,1:2),X(2,1:2),'k:'); set(h,'MarkerSize',6), grid,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) end clear clf a=-20*100 a=-20*100;b=-a;c=a;d=b;p=0.1; n=40; A=[1/5,99/100;1,0]; %A=[0.9,0.2;0.1,0.8]; %A=[2,1;0,1]; %A=[1/2,0;0,1/3]; %A=[ 0 1/2; -1/2 0] axis([a b c d]), grid,hold on button=1 while button==1 [xi,yi,button]=ginput(1); plot(xi,yi,'ko'),hold on, X0=[xi;yi]; X=X0; for i=1:n X=[A*X,X0]; h=plot(X(1,1),X(2,1),'g.',X(1,1:2),X(2,1:2),'k:');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, %pause end %axis([0 5 0 5]) p=1000; 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; plot(x,z1,'linewidth',2), hold on h=quiver([500,501]',[-1000,-1001]',[pc(1,1),0]',[pc(2,1),0]',p) set(h,'linewidth',2,'color','red'),

矩阵的相似化简(程序zxy8_2.m画直线) p=60; 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; plot(x,z1,'linewidth',2), h=quiver([500,501]',[-1000,-1001]',[pc(1,1),0]',[pc(2,1),0]',p) set(h,'linewidth',2,'color','red'),

矩阵的相似化简(高维线性离散动力系统) 动力系统理论的基本目的是了解迭代过程的最终或渐进性态， 如果迭代过程是一时间为自变量的微分方程，则该理论试图预言微分方程的解在遥远的将来或遥远的过去的最终性态。 如果过程是像函数迭代的离散过程，则这种理论希望了解随着n变大，迭代点的最终性态。 就是说，动力系统提出了一个听起来非数学的问题：这些点跑到哪里去？当它们到达那里又在干些什么？

矩阵的相似化简(高维线性离散动力系统) 记线性动力系统L(x)，它是 L(x)=Lx 其中 L 是适当维数的方阵， 我们只考虑2～3维的情形。

矩阵的相似化简(高维线性离散动力系统)

矩阵的相似化简(高维线性离散动力系统) 考察前面在2 维情况的例子中的标准形 A=[2, 0; 0,1/2] 它具有混合特征值 运行zxy8_3.m 进行观察迭代情形。 观察B＝PAP^-1迭代的情况，比较异同。

主成分分析和线性变换 气象分析预报中需要分析很多变量指标。 何抓住主要特点，用较少的指标代替原来较多的指标，又能综合反映原来较多的指标信息，就是实际工作者所关心的问题。（降维） 主成分分析方法为此提供了一种有效的手段。

主成分分析和线性变换 设有两个变量指标： X1：代表某地10月副高强度指数； X2：代表该地10月的副高面积指数。

主成分分析和线性变换 运行zxy8_4.m 画出散点图。 如何找到适当的坐标轴，使信息损失尽可能小？

主成分分析和线性变换

主成分分析和线性变换 求协方差矩阵 求正交矩阵P，满足 运行观察 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) [P,latent,explained] = pcacov(C) 运行观察

主成分分析和线性变换(运行zxy8_4)

线性变换 运行zxy8_5

数学实验与Matlab 制作人 周 晓 阳 华中科技大学数学系 zxyhust@netease.com 周 晓 阳 华中科技大学数学系 zxyhust@netease.com http://math.hust.edu.cn/~zhou/index.html

Galton钉板和二项分布 分布列的意义

Galton钉板模拟 英国生物统计学家Galton设计了Galton板 右边是一个5层Galton钉板

Galton钉板模拟(原理) 在一板上钉有n排钉子 自顶端扔进一小球任其自由下落 在下落过程中小球碰到钉子，左右落下的机会相等 最后小球落入底板中的某一个格子 图中用0 1 2 3 4 5表示这6个格子

Galton钉板模拟(博彩问题) 在每一格子中放上适当价值的奖品 如依次为 10 1 0.2 0.2 1 8 (元） 如依次为 10 1 0.2 0.2 1 8 (元） 扔一次小球你要付1元给庄主 如果小球落入某个格子 你将获得相应价值的奖品 你合算吗？庄主会赚钱吗？

Galton钉板模拟(扔1万个小球) 小球落入哪一个格子是不确定的 所以要计算落入每一个格子的可能性 这些小球将堆积起来 小球的堆积形状告诉了我们什么呢？

Galton钉板模拟(程序zxy9_1.m) (1)确定钉子的位置：将钉子的横、纵坐标存储在两个矩阵X和Y之中。 (2)选取0<p<1,将[0,1]区间分成两段:[0, p) 和 [p，1]。 (3)产生随机数r=rand(1，1),如果r<p，让小球向右落下；若r>p，让小球向左落下。(见备注） (4)将这一过程重复n次，并用直线连接小球落下时所经过的点，这样就模拟了小球从顶端随机地落入某一格子的过程。 r=rand(1,1),p=0.5 r=rand(1,5)-p

Galton钉板模拟(程序zxy9_1.m) (5)模拟小球堆积的形状。 输入扔球次数 m （例如 m =100) 计算落在第 i 个格子的小球数在总球数m中所占的比例f(i) 当模拟结束时，就得到了频率: f(0),f(1),f(2),f(3),f(4) 画出它们的图形。就是小球堆积的形状

Galton钉板模拟(程序zxy9_1.m) 制作动画矩阵数据； （6）动画指令结构 moviein(n)：创建动画矩阵； getframe ：拷贝动画矩阵； movie(Mat,m)：播放动画矩阵 m 次，（zxy7_6演示、讲解，备注） 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次 %这很有趣

Galton钉板模拟(程序zxy9_1.m) 运行zxy9_1.m

一个模拟结果 扔100个小球 向右概率p=0.5 要改变参数观察一下不同的模拟结果吗？这很容易.自己动手试试吧

随机变量及其分布 当你扔小球时，你和庄家关心什么？ ？？？？？？？？？？？？ 对，是小球落入格子的编号数 X （有些绕口，但很重要） 在投球前，你不能说你的小球会落在第0个格子。 但你可以说小球将落在第X个格子 X是一个随机数 是概率论中重要的讨论对象-----随机变量！！！

随机变量及其分布 实际上, 更应该关心的是 X 的分布列 分布列是小球落在各格子里的概率： P(X=0), P(X=1), P(X=2), P(X=3), P(X=4)，P(X=5) 想一想，它是不是表现了大量投球后小球堆积的极限形状呢 备注（比较频率和概率） ballnum1/3 f= binopdf([0 1 2 3 4 5], 5, 0.5)

Bernoulli试验和二项分布 不要把Galton钉板简单地当作消遣 它是一个有用的概率模型 n重Bernoulli试验的成功次数X 服从二项分布B(n, p). 上面模拟对应于n=5, p=0.5的情形

二项分布列 随机变量X～B(n,p)，则它的分布列为

统计工具箱 用指令f= binopdf(x, n, p)可计算二项分布的分布列 用F= binocdf(x, n, p)可计算二项分布的分布函数 用R = binornd(n,p,s,m)模拟m个二项随机数

观察二项分布列 运行binopdfcompare.m 固定n , 改变p值，观察二项分布列的形状 看一看：改变向右的概率，小球的堆积形状是怎样的？ 增加钉板层数n，作进一步观察。

模拟二项分布随机变量 用R = binornd(5,0.5,1,1)模拟了一次投球的结果。多次运行它，看看你的运气。 用R = binornd(5,0.5,1,m)成批模拟了m次投球结果，看看它的堆积形状。(运行simulatingGalton.m)

数学期望和平均收益 奖品的设置 格子编号 0 1 2 3 4 5 奖品价值 5 1 0.2 0.2 1 5 格子编号 0 1 2 3 4 5 奖品价值 5 1 0.2 0.2 1 5 观察：模拟5000次抽奖过程，抽奖一次支付1元，按上表获得回报。计算总收益和一次抽奖所得的平均收益 计算理论均值 备注 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 %计算一次抽奖的平均回报 。 --------------------------------------- f=[5, 1, 0.2, 0.2, 1, 5]; x=[0:1:n]; pu=binopdf(x,n,p);EF=sum(f.*pu)

数学期望和平均收益 格子编号X 0 1 2 3 4 5 分布列p p0 p1 p2 p3 p4 p5 价值函数f f0 f1 f2 f3 f4 f5 频率 m0 / m m1/ m m2 / m m3 / m m4 / m m5/ m

应用、思考和练习（电力供应问题 ） 某车间有200台车床互相独立的工作， 由于经常需要检修、测量、调换刀具等种种原因需要停车，这使每台车床的开工率只有60％。而每台车床在开动时需耗电1kW， 显然向该车间供电200kW可以保证有足够电力供这些车床使用， 但是在电力比较紧张的情况下，给这个车间供给电力太多将造成浪费，太少又影响生产。 如何解决这一矛盾？(模拟法？）

应用、思考和练习（ 废品问题 ） 一工厂生产某种大量耗用的零件，经过统计方法估计出废品率为p=0.015。工厂将100个零件装成一盒，销售给用户。但是不久接到用户反馈意见：声称100盒产品大约只有22盒全是正品，用户希望将这个比例提高到80盒左右。 管理人员希望采取某种措施满足用户的要求。为此他们进行了技术分析，希望减少废品率，但是这样做成本太高而不现实。 一名管理人员提出了一个简单想法，他认为可以在每盒产品的100个零件之外奉送额外的若干零件，这样希望基本保证用户得到100个正品，从而满足他们提出的要求。这一方法可行吗？请用概率论知识对此进行分析。

应用、思考和练习（ 奖品的设计） 所有抽奖模型都是要赚钱的，没有人想 花费精力却一无所获，甚至亏本。 所有抽奖模型都是要赚钱的，没有人想 花费精力却一无所获，甚至亏本。 对Galton钉板抽奖模型，如何在各个 格子中摆放适当价值的奖品获取利润？ 试提出你的设计方案。

热轧机的调整（正态分布） 轧制钢材要经过两道工序，第一道是粗轧(热轧)，第二道是精扎(冷轧)。 热轧机可以调整它的额定长度值，控制轧出钢材的平均长度。 成品材具有一个规定的长度l，如果热轧出的钢材长于l,精轧时就把多余的长度切掉。 如果粗轧时轧出的钢材长度比l短，则整根钢材报废。 精轧设备精度很高，轧出的成品可以认为是完全符合规定长度要求的。 问题是如何调整热轧机的额定值，使得钢材浪费最小。

热轧机的调整（正态模型）

热轧机的调整（正态模型） 随机变量X ～N(μ,σ2)描述了热轧机轧出的钢材长度。 Y = normpdf(X,MU,SIGMA), （运行normpdfobs.m） sigma=0.2;mu=2.32 a=mu+4*sigma;b=mu+4*sigma; x=linspace(a,b,50); f=normpdf(x,mu,sigma); plot(x,f)

热轧机的调整（模拟热轧机工作） 运行观察程序zxy9_2.m 用正态分布发生器normrnd 模拟热轧的结果。 观察哪些钢材需要切割，切割多少？哪些则将整根报废。 观察正态参数的对轧钢结果的影响作为对比观测，

热轧机的调整（优化zxy9_3.m,find） 目标函数W（μ)为得到一根成品材所浪费的平均钢材长度 设成品材长度为l=2,热轧机的精度为sigma=0.2. 仿真热轧机实际操作过程：对给定的μ值，模拟热轧机轧制一批钢材，计算每得一根成品材所浪费的平均钢材长度，得到的估计值； 通过改变μ的值，绘制的变化曲线，确定使浪费达到最小的额定长度μ0，这就是热轧机应该设置的最佳额定长度的估计值。

热轧机的调整（优化）

应用练习与思考（热轧机模型） 根据上面的结果，你得到什么印象。一根成品 材长为2 米，浪费为0.44783米，这样的结果 可以接受吗？你有办法减少浪费量吗？ 图中，曲线不是光滑的，这是正常的吗？不改 变程序中的参数，多次运行程序，结果会保持 不变吗？如何解释你观察到的情况？ 你可以加密mu的分点进行计算，也可以改变一 些参数重新运行程序，观察计算结果。 你有兴趣用解析的方法解决这一问题吗？模拟的方法和解析的方法有什么联系和区别？

应用练习与思考（热轧机模型） 如何进一步减少浪费？一种想法是提高热轧机的精度，固定l的值，减少均方差重新进行计算，看看会发生什么结果？ 你可以通过绘制均方差～浪费量曲线的方法为管理者提供决策依据。 另一种方法是改变问题的提法：可以考虑将成品材分为若干级别，例如分为两个档次。长度在l1 ~ l2之间的降级使用，只有当长度小于才l1整根报废。试对这一问题进行思考，选用合适的目标函数建立优化模型，使某种意义下的浪费达到最小。并结合前面已经计算的结果，为管理者提供一个可行性建议。

应用练习与思考（如何制定胖和瘦的标准) 描述某一类人群的身高，可以使用正态随机变量。这是因为和热轧机模型一样，大多数人的身高在人群的平均身高附近波动，呈现着一种对称分布：特别高或特别矮的人是很少的，与正态分布体现的中间大两头小的特点比较吻合。 同理，也可用正态随机变量来描述这一类人群的体重。

应用练习与思考（如何制定胖和瘦的标准) simulationhigh_weight.m, hist 观察：模拟身高和体重的数据（最好你自己收集数据） 产生n个正态随机数X，代表n个人的身高； 再产生n个零均值的正态随机扰动U 用Y=rX+s+U计算n个随机数Y，代表着这n个人的体重。 simulationhigh_weight.m, hist clear,clf n=1000;x=normrnd(170,4.5,1,n); u=normrnd(0,0.5,1,n); y=0.36*x+u; 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');

应用练习与思考（如何制定胖和瘦的标准) 绘制二维直方图和二元正态分布密度函数图象。 运行程序zxy9_4.m改变参数，观察二维直方图和理论分布的图形和身高和体重的概率关系。

如何制定胖和瘦的标准(条件正态分布法)

如何制定胖和瘦的标准(回归分析) 这时将出现一图形窗口， 用鼠标水平移动十字线 可得到不同身高对应的预测值， alpha=0.01; polytool(x',y',1,alpha) 这时将出现一图形窗口， 用鼠标水平移动十字线 可得到不同身高对应的预测值， 在正上方degree输入框中可改变拟合多项式次数， 在左下角的列表框中有不同的输出选项， 可输出相应的参数，例如选择all将输出参数: beta (多项式系数)、betaci (多项式系数置 信区间)、yhat（预测值）、yci（预测值置信区 间）residuals（残差）。 alpha=0.01; polytool(x',y',1,alpha)