数学实验与Matlab

1、数学实验与数学实验与Matlab 05101520250102030-10-50510 MATLAB是一种交互式的以矩阵为基础的系统计算平台,它用于科学和工程的计算与可视化。 Matlab的含义是矩阵实验室(Matrix Laboratory),是美国Math Work公司于1982年推出的一套高性能的数值计算和可视化软件,它集数值分析、矩阵计算、信号处数值分析、矩阵计算、信号处理和图形显示理和图形显示于一体,已发展成为国际上最优秀的科技应用软件之一。 矩阵矩阵是是MATLAB的核心的核心. 目目 录录 MatlabMatlab基本指令基本指令向量的创建和运算向量的创建和运算x1=1 2 4,

2、 x2=1,2,1, x3=x1x1=1 2 4, x2=1,2,1, x3=x1 x1=3.4 x1=3.4: :6.7,6.7, x2=3.4x2=3.4:2:2:6.7,6.7, x3=2.6x3=2.6:-0.8:-0.8:0 0指令指令x=linspace(a,b,nx=linspace(a,b,n) ) 在在a,ba,b区间产区间产生生 n n 个等分点个等分点( (包括端点包括端点) ) 例如：例如：x=x=linspacelinspace(0,1,5)(0,1,5)结果结果 x = x = 0 0.2500 0.5000 0.7500 1.0000 0 0.2500 0.500

3、0 0.7500 1.0000在在MatlabMatlab窗口创建向量后并运行后，窗口创建向量后并运行后，向量就存在于工作空间向量就存在于工作空间（WorkspaceWorkspace），可以被调用。），可以被调用。 设三维向量设三维向量x=x1 x2 x3; y=y1 y2 y3; x=x1 x2 x3; y=y1 y2 y3; ，a, ba, b为标量。为标量。向量的向量的数乘数乘：a a* *x=ax=a* *x1 ax1 a* *x2 ax2 a* *x3x3向量的向量的平移平移： x+b=x1+b x2+b x3+bx+b=x1+b x2+b x3+b向量向量和和： x+y=x1+y

4、1 x2+y2 x3+y3x+y=x1+y1 x2+y2 x3+y3向量向量差差： x-y=x1-y1 x2-y2 x3-y3 x-y=x1-y1 x2-y2 x3-y3 x x. .* *y=x1y=x1* *y1 x2y1 x2* *y2 x3y2 x3* *y3 (y3 (乘积乘积) )x x./ ./y=x1/y1 x2/y2 x3/y3 (y=x1/y1 x2/y2 x3/y3 (右除，右右除，右边的边的y y做分母做分母) )x x. .y=y1/x1 y2/x2 y3/x3 (y=y1/x1 y2/x2 y3/x3 (左除，左左除，左边的边的x x做分母做分母) )x x.5=x

5、15 x25 x35 (5=x15 x25 x35 (乘幂乘幂) ) 2 2.x=2x1 2x2 2x3x=2x1 2x2 2x3x x.y=x1y1 x2y2 x3y3y=x1y1 x2y2 x3y3MatlabMatlab有许多内部函数，可直接作用于有许多内部函数，可直接作用于向量产生一个同维的函数向量。向量产生一个同维的函数向量。如：如：x=linspace(0,4x=linspace(0,4* *pi,100);pi,100);（产生（产生100100维向量维向量x x） y=sin(xy=sin(x); ); (y(y也自动为也自动为100100维向量维向量) ) y1=sin(x)

6、.2;y1=sin(x).2; y2=exp(-x). y2=exp(-x).* *sin(x); sin(x); 观察结果观察结果 直接输入法创建简单矩阵。直接输入法创建简单矩阵。 A=1 2 3 4; 5 6 7 8; 9 10 11 12A=1 2 3 4; 5 6 7 8; 9 10 11 12 B=-1.3, B=-1.3,sqrt(3)sqrt(3);(1+2);(1+2)* *4/5,4/5,sin(5)sin(5); ;exp(2exp(2),6),6 观察运行结果观察运行结果 A = 1 2 3 4 5 6 7 8 9 10 11 12 B = -1.3000 1.7321

7、2.4000 -0.9589 7.3891 6.0000矩阵的运算矩阵的运算( (矩阵的加减、数乘、乘积等矩阵的加减、数乘、乘积等) )A, A, A_transA_trans=A=A H=1 2 3 ; 2 1 0 ; 1 2 3 ,H=1 2 3 ; 2 1 0 ; 1 2 3 ,K=1 2 3 ; 2 1 0 ; 2 3 1K=1 2 3 ; 2 1 0 ; 2 3 1h_deth_det= =detdet(H), k_det(H), k_det= =detdet(K(K), ),H_invH_inv= =invinv(H), K_inv(H), K_inv=K=K-1-1矩阵的运算矩阵

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

9、右除) )例例1 1：求矩阵方程：求矩阵方程：设设A A、B B满足关系式：满足关系式：ABAB2B+A,2B+A,求求B B。其中其中A=3 0 1; 1 1 0; 0 1 4A=3 0 1; 1 1 0; 0 1 4。解解：有：有(A-2I)B(A-2I)BA A 程序程序 ： A=3 0 1; 1 1 0;0 1 4;A=3 0 1; 1 1 0;0 1 4;B=B=invinv(A-2(A-2* *eyeeye(3)(3)* *A,A,B=(A-2B=(A-2* *eyeeye(3)A(3)A 观察结果：观察结果：例例2 2 取出取出A A的的1 1、3 3行和行和1 1、3 3列的交

10、叉处元列的交叉处元素构成新矩阵素构成新矩阵A1A1。解：程序解：程序A=1 0 1 1 2;0 1 -1 2 3;A=1 0 1 1 2;0 1 -1 2 3; 3 0 5 1 0;2 3 1 2 1, 3 0 5 1 0;2 3 1 2 1, vr=1, 3; vc vr=1, 3; vc=1, 3;=1, 3;A1=A(vr, vcA1=A(vr, vc) ) 观察结果观察结果）例例3 3 将将A A分为四块，并把它们赋值到矩阵分为四块，并把它们赋值到矩阵B B中，观察运行后的结果。中，观察运行后的结果。解：程序解：程序nA11=A(1:2,1:2),A12=A(1:2,3:5),A11=

11、A(1:2,1:2),A12=A(1:2,3:5),nA21=A(3:4,1:2),A22=A(3:4,3:5)A21=A(3:4,1:2),A22=A(3:4,3:5)nB=A11 A12;B=A11 A12;n A21 A22 A21 A22 运行结果运行结果例例4 4 （1 1）修改矩阵）修改矩阵A A，将它的第，将它的第1 1行变为行变为0 0。解：程序：解：程序： A=1 0 1 1 2;0 1 -1 2 3;A=1 0 1 1 2;0 1 -1 2 3; 3 0 5 1 0;2 3 1 2 1 3 0 5 1 0;2 3 1 2 1 A(1,:)=0 0 0 0 0; A A(1,

12、:)=0 0 0 0 0; A （2 2）删除上面矩阵）删除上面矩阵A A的第的第1 1、3 3行。行。 程序：程序：n A(1,3,:)= A(1,3,:)= 全全1 1阵阵 n ones(n), ones(m,n), ones(size(A)ones(n), ones(m,n), ones(size(A)全零阵：全零阵：n zeros(n) zeros(n) ，zeros(m,n), zeros(size(A) zeros(m,n), zeros(size(A) n常常用于对某个矩阵或向量赋常常用于对某个矩阵或向量赋0 0初值初值单位阵：单位阵：neye(n)eye(n)，eye(m,n)

13、 eye(m,n) 随机阵：随机阵： nrand(m,n)rand(m,n)， rand(n)=rand(n,nrand(n)=rand(n,n) )用于随机模用于随机模拟，常和拟，常和rand(seed,krand(seed,k) )配合使用配合使用。det(Adet(A) ) ： 方阵的行列式；方阵的行列式；rank(A)rank(A)： 矩阵的秩；矩阵的秩；eig(Aeig(A) )： 方阵的特征值和特征向量；方阵的特征值和特征向量；trace(A)trace(A)： 矩阵的迹；矩阵的迹；rref(Arref(A) )： 初等变换阶梯化矩阵初等变换阶梯化矩阵A Asvd(Asvd(A)

14、)： 矩阵奇异值分解。矩阵奇异值分解。1. 1. 矩阵元素的标识矩阵元素的标识 ： A(i,j)A(i,j)表示矩阵表示矩阵A A 的第的第 i i 行行 j j 列的元素；列的元素； 2. 2. 向量标识方式向量标识方式 A(vr,vcA(vr,vc) )： vr=i1,i2,ikvr=i1,i2,ik、vc=j1,j2,juvc=j1,j2,ju 分别分别是含有矩阵是含有矩阵A A的行号和列号的单调向量。的行号和列号的单调向量。 A(vr,vcA(vr,vc) )是取出矩阵是取出矩阵A A的第的第i1,i2,iki1,i2,ik行行与与j1,j2,juj1,j2,ju列交叉处的元素所构成新

15、矩阵。列交叉处的元素所构成新矩阵。1.1.当数据为当数据为行向量或列向量行向量或列向量时，函数对时，函数对整整个向量个向量进行计算进行计算. .2.2.当数据为当数据为矩阵矩阵时，命令对时，命令对列列进行计算，进行计算，即把每一即把每一列列数据当成同一变量的不同观察数据当成同一变量的不同观察值。值。常用的命令：常用的命令： maxmax( (求最大求最大) )、minmin( (求最小求最小) )、meanmean( (求平均求平均值值 ) ) 、 s u ms u m ( ( 求 和求 和 ) ) 、 s t ds t d ( ( 求 标 准 差求 标 准 差 ) ) 、cumsumcums

16、um( (求累积和求累积和) )、medianmedian( (求中值求中值) )、diffdiff( (差分差分) )、sortsort( (升序排列升序排列) )、sortrowssortrows( (行行升序排列升序排列) )等等等等。例例5 5 观察：生成一个观察：生成一个3 36 6的随机数矩阵，的随机数矩阵，并将其各列排序、求各列的最大值与各列并将其各列排序、求各列的最大值与各列元素之和。元素之和。解：程序解：程序nA=A=randrand(3,6),(3,6),nAsortAsort= =sortsort(A), (A), nAmax=Amax=maxmax(A(A), ), n

17、AsumAsum= =sumsum(A(A) ) 观察结果观察结果 plot(x,y)plot(x,y)：n以以x x为横坐标、为横坐标、y y为纵坐标绘制二维图形为纵坐标绘制二维图形nx,yx,y是同维数的向量；是同维数的向量； plot(yplot(y) )：n相当于相当于x=1,2,length(y)x=1,2,length(y)时情形。时情形。)sin(cos(sin()(xxxxf)sin(cos(sin(2 . 0)(xxxxxg问：问：f (x), g (x)f (x), g (x)是周期函数吗？观察它们的图象。是周期函数吗？观察它们的图象。解：程序解：程序clfclf, x=l

18、inspace(0,8, x=linspace(0,8* *pi,100);pi,100);y1=sin(x+cos(x+sin(x); y1=sin(x+cos(x+sin(x); y2=0.2y2=0.2* *x+sin(x+cos(x+sin(xx+sin(x+cos(x+sin(x););plot(x,y1,k:,x,y2,k-) plot(x,y1,k:,x,y2,k-) legend(sin(x+cos(x+sin(x),0.2x+sin(xlegend(sin(x+cos(x+sin(x),0.2x+sin(x+cos(x+sin(x),2)+cos(x+sin(x),2)例例1

19、 1 令令1. plot(x,y1;y2;),1. plot(x,y1;y2;),n x x是横坐标向量，是横坐标向量，y1;y2;y1;y2;是由若干函数是由若干函数的纵坐标拼成的矩阵；的纵坐标拼成的矩阵；2. plot(x1,y1), 2. plot(x1,y1), hold onhold on, , plot(x2,y2), plot(x2,y2), hold offhold off3. plot(x1,y1,x2,y2,) 3. plot(x1,y1,x2,y2,) 4. plotyy4. plotyyn两个坐标系，用于绘制不同尺度的函数。两个坐标系，用于绘制不同尺度的函数。plot(

20、x,y,plot(x,y,颜色线型点形颜色线型点形) )plot(x,y,plot(x,y,颜色线型点形颜色线型点形,x,y,x,y,颜颜色线型点形色线型点形, ), )句柄图形和句柄图形和setset命令改变属性值，可套命令改变属性值，可套用：用：nh=plot(x,y),h=plot(x,y), set(h, set(h,属性属性, ,属性值属性值, ,属性属性, ,属性值属性值,),)n或或plot(x,yplot(x,y, ,属性属性,属性值属性值) )设置图形对象设置图形对象的属性。的属性。n线宽：线宽：LineWidthLineWidthn点的大小：点的大小： MarkerSize

21、MarkerSizen线型：线型：LineStyleLineStylen颜色：颜色：ColorColor程序程序nh=plot(0:0.1:2h=plot(0:0.1:2* *pi,sin(0:0.1:2pi,sin(0:0.1:2* *pi); pi); nset(h,LineWidth,5,color,red); set(h,LineWidth,5,color,red); ngrid ongrid on观察结果观察结果grid on grid on 指令为图形窗口加上网格线指令为图形窗口加上网格线axisaxis指令指令 axis(xmin xmax ymin ymaxaxis(xmin

22、xmax ymin ymax) )： 设定二维图形的设定二维图形的x x和和y y坐标的范围；坐标的范围； axis(xmin xmax ymin ymax zmin axis(xmin xmax ymin ymax zmin zmaxzmax) )： 设定三维图形的坐标范围设定三维图形的坐标范围 ； 其中其中xminxxmaxxminxxmax， yminyymaxyminyymax ， zminzzmaxzminzzmax。title(title(图形标题图形标题) )；xlabelxlabel（x x轴名称轴名称）；）；ylabelylabel（y y轴名轴名称称）；）；zlabelzl

23、abel（z z轴名称轴名称）；）；text(text(说明文字说明文字) )：创建说明文字；：创建说明文字；gtextgtext（ 说明文字说明文字 ）：用鼠标在特定位置输入文）：用鼠标在特定位置输入文字。字。文字标注常用符号：文字标注常用符号： pi pi （ ）；）；alpha alpha （ ）；）；beta beta （ ）；）； leftarrowleftarrow（左箭头）（左箭头） rightarrowrightarrow（右箭头）；（右箭头）； bullet bullet （点号）（点号） subplot(m,n,k) subplot(m,n,k)命令命令 在图形区域中显示多

24、个图形窗口，在图形区域中显示多个图形窗口，mm为上下分为上下分割数，割数，n n为左右分割数，为左右分割数，k k为第为第k k子图编号。子图编号。例例：将一个图形分为：将一个图形分为4 4个子图，在第个子图，在第k k个子图画个子图画sin(kxsin(kx) ) 的图象的图象. .程序：程序： clfclf，b=2b=2* *pi; x=linspace(0,b,50);pi; x=linspace(0,b,50); for k =1:4 for k =1:4 y=sin(k y=sin(k * * x);x); subplotsubplot(2,2,k),plot(x,y),axis(0

25、,2(2,2,k),plot(x,y),axis(0,2* *pi,-1,1)pi,-1,1) end endclfclf：清除图形窗口已有的内容：清除图形窗口已有的内容. . shgshg：显示图形窗口。：显示图形窗口。clearclear、 clear xclear x：清除工作空间的已：清除工作空间的已有变量。有变量。figure(n): figure(n): 打开第打开第n n个图形窗口个图形窗口 help:help: : : 续行号续行号基本步骤：基本步骤： 1. 1. 生成二维网格点生成二维网格点 2. 2. 计算函数在网格点上的值计算函数在网格点上的值 3. 3. 绘制函数图形绘

26、制函数图形观察观察meshgridmeshgrid指令的效果。指令的效果。程序：程序： a=-3;b=3;c=-3;d=3;n=10;a=-3;b=3;c=-3;d=3;n=10; x=linspace(a,b,n x=linspace(a,b,n); ); y=linspace(c,d,n y=linspace(c,d,n); ); X,Y= X,Y=meshgridmeshgrid(x,y(x,y); ); plot(X,Y plot(X,Y,+),+)观察结果观察结果 如，如，z=peaks(X,Yz=peaks(X,Y); ); mesh(X,Y,z) mesh(X,Y,z) ：n在三

27、维空间中绘出由在三维空间中绘出由(X,Y,z)(X,Y,z)表示的曲面表示的曲面；meshz(X,Y,zmeshz(X,Y,z) )：n除了具有除了具有meshmesh的功能外，还画出上下高度线，的功能外，还画出上下高度线，meshc(X,Y,zmeshc(X,Y,z) )：n除了具有除了具有meshmesh的功能外，还在曲面的下方画出函数的功能外，还在曲面的下方画出函数z=f(x,y)z=f(x,y)的等值线图，的等值线图，surf(X,Y,z)surf(X,Y,z)：n也是三维绘图指令，与也是三维绘图指令，与meshmesh的区别在于的区别在于meshmesh绘出彩绘出彩色的线，色的线，s

28、urfsurf绘出彩色的面，绘出彩色的面，）contour(X,Y,z,n)contour(X,Y,z,n)：nn n条等高线，条等高线，n n可缺省；可缺省；contourf(X,Y,z,ncontourf(X,Y,z,n) )：n等值线间用不同的颜色填满，有更好的视觉效果；等值线间用不同的颜色填满，有更好的视觉效果； contour3(X,Y,z,n)contour3(X,Y,z,n)：n在三维空间画出等值线图；在三维空间画出等值线图；colorbarcolorbar：n将颜色与函数值对应起来显示在图中将颜色与函数值对应起来显示在图中。一条空间曲线可以用矢量函数表示为一条空间曲线可以用矢量

29、函数表示为, , )(),(),()( :ttztytxtrn它的速度矢量表现为曲线的切矢量:)(,)(,)()(dttdzdttdydttdxdtdrtv plot3(x,y,z)plot3(x,y,z)：n绘制三维空间曲线，用法和绘制三维空间曲线，用法和plotplot类似。类似。quiver(X,Y,u,v)quiver(X,Y,u,v)：绘制二维矢量，：绘制二维矢量，n在坐标矩阵点在坐标矩阵点X,YX,Y处绘制矢量处绘制矢量u,v, u,v, 其中其中u u为矢量的为矢量的x x坐标，坐标，v v为矢量的为矢量的y y 坐标，其维数不坐标，其维数不小于小于2 2。quiver3quiv

30、er3（X,Y,Z,u,v,wX,Y,Z,u,v,w）：）：n绘制三维矢量，用法与绘制三维矢量，用法与quiverquiver类似。类似。gradientgradient：n Fx,Fy,FzFx,Fy,Fz=gradient(F)=gradient(F)为函数为函数F F数值梯度数值梯度)3/1 (3 ,)3/2( ,()6 ,2 ,2()(34320232ttttdttttttrtn很显然飞行曲线方程为： 程序：程序：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.-,linewi

31、dth,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, hold off 用程控铣床加工机翼断面的下轮廓线时用程控铣床加工机翼断面的下轮廓线时每一刀只能沿每一刀只能沿x x方向和方向和y y方向走非常小的一步。方向走非常小的一步。表表3-13-1给出了下轮廓线上的部分数

32、据给出了下轮廓线上的部分数据但工艺要求铣床沿但工艺要求铣床沿x x方向每次只能移动方向每次只能移动0.10.1单位单位. . 这时需求出当这时需求出当x x坐标每改变坐标每改变0.10.1单位时的单位时的y y坐标。坐标。试完成加工所需的数据，画出曲线试完成加工所需的数据，画出曲线. . 表表3-1 3-1 机翼断面的下轮廓线上的数据机翼断面的下轮廓线上的数据x0357911 12 13 14 15y01.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6n已知有已知有n +1n +1个节点（个节点（x xj j,y,yj j），），j = 0,1, nj = 0,1, n，其

33、中其中x xj j互不相同，节点互不相同，节点(x(xj j, y, yj j) )可看成由某个函可看成由某个函数数 y= fy= f（x x）产生；）产生；n构造一个相对简单的函数构造一个相对简单的函数 y=P(x)；n使使P P通过全部节点，通过全部节点， 即即 P (xP (xk k) = y) = yk k， k k=0,1, n =0,1, n ；n用用P (x)P (x)作为作为 函数函数f ( x )f ( x )的的 近似。近似。yi=interp1(x，y，xi，method)插值方法插值方法被插值点被插值点插值节点插值节点xixi处的插处的插值结果值结果nearest ：最

34、邻近插值：最邻近插值linear ： 线性插值；线性插值；spline ： 三次样条插值；三次样条插值；cubic ： 立方插值。立方插值。缺省时：缺省时： 线性插值。线性插值。 注意：所有的插值方法都要求注意：所有的插值方法都要求x x是单调的，并且是单调的，并且xi不不能够超过能够超过x的范围。的范围。用用MATLABMATLAB作一维插值计算作一维插值计算x0=0 3 5 7 9 11 12 13 14 15 ;y0=0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 ;x=0:0.1:15;y=interp1(x0,y0,x,spline);plot(x0,y0

35、,k+,x,y,r)grid on 要求要求x0,y0 x0,y0单调；单调；x x，y y可取可取为矩阵，或为矩阵，或x x取取行向量，行向量，y y取为列向量，取为列向量，x,yx,y的值分别不能超出的值分别不能超出x0,y0 x0,y0的范围。的范围。z=interp2(x0,y0,z0,x,y,method)被插值点插值方法插值节点被插值点的函数值nearest nearest 最邻近插值最邻近插值linear linear 双线性插值双线性插值cubic cubic 双三次插值双三次插值缺省时缺省时, , 双线性插值双线性插值用用MATLABMATLAB作网格节点数据的插值作网格节点

36、数据的插值( (二维二维) ) cz =griddata（x，y，z，cx，cy，method） 要求要求cxcx取行向量，取行向量，cycy取为列向量取为列向量。被插值点插值方法插值节点被插值点的函数值nearest nearest 最邻近插值最邻近插值linear linear 双线性插值双线性插值cubic cubic 双三次插值双三次插值v4 Matlab提供的插值方法提供的插值方法缺省时缺省时, , 双线性插值双线性插值用用MATLABMATLAB作散点数据的插值计算作散点数据的插值计算某海域上频繁地有各种吨位的船只经过。为保证船只的航行安全，有关机构在低潮时对水深进行了测量,表3-

37、8是他们提供的测量数据：n表3-8. 水道水深的测量数据x 129.0 140.0 103.5 88.0 185.5 195.0 105.5y 7.5 141.5 23.0 147.0 22.5 137.5 85.5z 4 8 6 8 6 8 8x 157.5 107.5 77.0 81.0 162.0 162.0 117.5y -6.5 -81.0 3.0 56.5 -66.5 84.0 -33.5z 9 9 8 8 9 4 9其中（其中（x, yx, y）为测量点，）为测量点，z z为（为（x, yx, y）处的水深）处的水深( (英尺英尺) )，水深，水深z z是区域坐标（是区域坐标（x

38、, yx, y）的函数）的函数z= z z= z (x, y(x, y），），船的吨位可以用其吃水深度来反映，分为船的吨位可以用其吃水深度来反映，分为 4 4英英尺、尺、4.54.5英尺、英尺、5 5英尺和英尺和 5.55.5英尺英尺 4 4 档。档。 航运部门要在矩形海域（航运部门要在矩形海域（7575，200200）（5050，150150）上为不同吨位的航船设置警示标记。）上为不同吨位的航船设置警示标记。请根据测量的数据描述该海域的地貌，并绘制请根据测量的数据描述该海域的地貌，并绘制不同吨位的警示线，供航运部门使用。不同吨位的警示线，供航运部门使用。x=129 140 103.5 88

39、185.5 195 105.5 157.5 107.5 77 81 162 162 117.5;y=7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5;z=-4 -8 -6 -8 -6 -8 -8 -9 -9 -8 -8 -9 -4 -9;cx=75:0.5:200;cy=-70:0.5:150;cz=griddata(x,y,z,cx,cy,cubic);meshz(cx,cy,cz), xlabel(X),ylabel(Y),zlabel(Z)figure(2),contour(cx,cy,cz,-5 -5); g

40、rid on,hold onplot(x,y,+)xlabel(X),ylabel(Y) 已知一组（二维）数据，即平面上已知一组（二维）数据，即平面上 n个点个点（xi,yi) i=1,n, 寻求一个函数（曲线）寻求一个函数（曲线）y=f(x), 使使 f(x) 在某种准则在某种准则下与所有数据点最为接近，即曲线拟合得最好。下与所有数据点最为接近，即曲线拟合得最好。 +xyy=f(x)(xi,yi)i i 为点为点（xi,yi) 与与曲线曲线 y=f(x) 的距离的距离最常用的方法是线性最小二乘拟和最常用的方法是线性最小二乘拟和n对给定的数据（对给定的数据（x xj j,y,yj j），），j

41、 = 0,1, nj = 0,1, n；n选取适当阶数的多项式，如二次多项式选取适当阶数的多项式，如二次多项式g(x)=ax2+bx+cg(x)=ax2+bx+c；n使使g(x)g(x)尽可能逼近（拟合）这些数据，但尽可能逼近（拟合）这些数据，但是不要求经过给定的数据（是不要求经过给定的数据（x xj j,y,yj j）；）；1. 1. 多项式多项式f(x)=a1xm+ +amx+am+1拟合指令拟合指令:a=polyfit(x,y,m)2.2.多项式在多项式在x x处的值处的值y y的计算命令：的计算命令： y=y=polyvalpolyval（a,xa,x）输出拟合多项式系数输出拟合多项式

42、系数a=a1,am,am+1 （数组）数组）输入同长度输入同长度数组数组X，Y拟合多项式拟合多项式次数次数多项式拟合指令多项式拟合指令即要求出二次多项式即要求出二次多项式:3221)(axaxaxf中中 的的123(,)Aa a a 使得使得:1121() iiif xy最小例例 对下面一组数据作二次多项式拟合对下面一组数据作二次多项式拟合xi00.10.20.30.40.50.60.70.80.91yi-0.4471.9783.286.166.167.347.669.589.489.3011.22）计算结果：）计算结果： = -9.8108, 20.1293, -0.0317解：解：用多项式

43、拟合的命令用多项式拟合的命令00.20.40.60.81-20246810120317.01293.208108.9)(2xxxf1）输入命令：）输入命令：x=0:0.1:1; y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2;A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,k+,x,z,r) %作出数据点和拟合曲线的图形作出数据点和拟合曲线的图形可化为多项式的非线性拟和可化为多项式的非线性拟和一般非线性最小二乘拟和一般非线性最小二乘拟和,实际上是无约实际上是无约束最优化问题束最优化问题命令：

44、命令：lsqcurvefit、lsqnonlin等等 （1 1）用不定积分计算定积分：）用不定积分计算定积分： 牛顿莱布尼兹公式牛顿莱布尼兹公式, ,但是有时得不但是有时得不到原函数；到原函数；（2 2）定义法，取近似和的极限：）定义法，取近似和的极限：n高等数学中不是重点内容高等数学中不是重点内容n但数值积分的各种算法却是基于定义建立的但数值积分的各种算法却是基于定义建立的 )x(x)f(limf(x)dxi1i1n0ii0ba数值微积分数值微积分（梯形公式和辛普森公式梯形公式和辛普森公式）trapz(x,ytrapz(x,y) )，按梯形公式计算近似积分；，按梯形公式计算近似积分； 其中步

45、长x=x0 x1 xn和函数值y=f0 f1 fn为同维向量；q = quad(fun,a,b,tol,trace,P1,P2,.) q = quad(fun,a,b,tol,trace,P1,P2,.) （低阶方法，辛普森自适应递归法求积分）（低阶方法，辛普森自适应递归法求积分）q = quad8(fun,a,b,tol,trace,P1,P2,.)q = quad8(fun,a,b,tol,trace,P1,P2,.)（高阶方法，自适应法（高阶方法，自适应法CotesCotes求积分）求积分） 在同样的精度下高阶方法在同样的精度下高阶方法quad8quad8要求的节点要求的节点较少。较少。

46、X = 0:pi/100:pi;Y = sin(X);Z1 = trapz(X,Y)Z2=quad(sin,0,pi)Z1 = 1.998z2 = 2.0000 x,yx,y=ode23(fun,tspan,y0,option)=ode23(fun,tspan,y0,option) （低阶龙格库塔函数）（低阶龙格库塔函数）x,yx,y=ode45(fun,tspan,y0,option) =ode45(fun,tspan,y0,option) （高阶龙格库塔函数）（高阶龙格库塔函数）其中，(1) tspan=t0，tf，t0、tf为自变量的初值和终值 (2) option用于设定误差限(缺省时

47、设定相对误差10-3, 绝对误差10-6), 命令为：options=odeset （reltol,rt,abstol,at）, rt，at：分别为设定的相对误差和绝对误差. 1、在解、在解n个未知函数的方程组时，个未知函数的方程组时，x0和和x均为均为n维向量，维向量，m-文件中的待解方程文件中的待解方程组应以组应以x的分量形式写成的分量形式写成. 2、使用、使用Matlab软件求数值解时，高软件求数值解时，高阶微分方程必须等价地变换成一阶微分阶微分方程必须等价地变换成一阶微分方程组方程组.注意注意: : 设位于坐标原点的甲舰向位于设位于坐标原点的甲舰向位于x轴上点轴上点A(1, 0)处的乙

48、舰发射导弹，导弹头始终对准处的乙舰发射导弹，导弹头始终对准乙舰。如果乙舰以最大的速度乙舰。如果乙舰以最大的速度v0(是常数是常数)沿平沿平行于行于y轴的直线行驶，轴的直线行驶，导弹的速度是导弹的速度是5v0，（1）画出导弹运行的曲线方程）画出导弹运行的曲线方程.（2）乙舰行驶多远时，）乙舰行驶多远时， 导弹将它击中？导弹将它击中？导弹追击问题由由(1),(2)消去消去t整理得模型整理得模型:(3) 151)1 (2yyx1.建立建立m-文件文件eq1.m function dy=eq1(x,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1/5*sqrt(1+y(1)2

49、)/(1-x); 2. 取取x0=0，xf=0.9999，建立主程序，建立主程序ff6.m如下如下： x0=0，xf=0.9999 x,y=ode23(eq1,x0 xf,0 0); plot(x,y(:,1),b.) hold on y=0:0.01:1; plot(1,y,b*) hold off 结论结论: 导弹大致在（导弹大致在（1，0.2）处击中乙舰）处击中乙舰2151 )1 (yyx)1/(15121221xyyyy令y1=y, y2=y1，将方程（3）化为一阶微分方程组。 最优化方法最优化方法 许多生产计划与管理问题都可以归纳为最优许多生产计划与管理问题都可以归纳为最优化问题化问

50、题, , 最优化模型是数学建模中应用最广泛的最优化模型是数学建模中应用最广泛的模型之一模型之一, ,其内容包括线性规划、整数线性规划、其内容包括线性规划、整数线性规划、非线性规划、动态规划、变分法、最优控制等非线性规划、动态规划、变分法、最优控制等. . 近几年来的全国大学生数学建模竞赛中，几近几年来的全国大学生数学建模竞赛中，几乎每次都有一道题要用到此方法乎每次都有一道题要用到此方法. .目标函数 约束条件 可行解域 xpjxhmixgtsxfji, 2 , 10)(, 2 , 1, 0)(. .)(min线性规划线性规划xpjxhmixgtsxfji, 2 , 10)(, 2 , 1, 0

51、)(. .)(min若)() ,() ,(xhxgxfii都是线性 函 数 ， 则 称 为 线 性 规 划 问题 ， 简 称L P问 题 用用MATLAB优化工具箱解线性规划优化工具箱解线性规划min z=cX bAXts. .1、模型：命令：x=linprog（c，A，b） 2、模型：min z=cX bAXts. .beqXAeq命令：x=linprog（c，A，b，Aeq, beq）注意：若没有不等式： 存在，则令A= ，b= .bAX 3、模型：min z=cX bAXts. .beqXAeq VLBXVUB命令：1 x=linprog（c，A，b，Aeq, beq, VLB，VUB）

52、 2 x=linprog（c，A，b，Aeq, beq, VLB，VUB, X0） 注意：若没有等式约束: , 则令Aeq= , beq= ; 其中X0表示初始点 beqXAeq4、命令：x,fval=linprog()返回最优解及处的目标函数值返回最优解及处的目标函数值fval.解解: 编写编写M文件文件xxgh2.m如下：如下： c=6 3 4; A=0 1 0; b=50; Aeq=1 1 1; beq=120; vlb=30,0,20; vub=; x,fval=linprog(c,A,b,Aeq,beq,vlb,vub)321)436(minxxxz32120030 xxx50010

53、,120111 .321321xxxxxxts 某车间有甲、乙两台机床，可用于加工三种工件。假定这两台车床的可用台时数分别为800和900，三种工件的数量分别为400、600和500，且已知用三种不同车床加工单位数量不同工件所需的台时数和加工费用如下表。问怎样分配车床的加工任务，才能既满足加工工件的要求，又使加工费用最低？ 单 位 工 件 所 需 加 工 台 时 数 单 位 工 件 的 加 工 费 用 车 床类 型 工 件1 工 件2 工 件3 工 件1 工 件2 工 件3 可 用 台时 数 甲 0.4 1.1 1.0 13 9 10 800 乙 0.5 1.2 1.3 11 12 8 900 例例（任务分配问题（任务分配问题)解解 设在甲车床上加工工件1、2、3的数量分别为x1、x2、x3，在乙车床上加工工件1、2、3的数量分别为x4、x5、x6。可建立以下线性规划模型：6543218121110913minxxxxxxz 6 , 2 , 1, 09003 . 12 . 15 . 08001 . 14 . 0500600400 x . .654321635241ixxxxxxxxxxxxtsiS.t.Xz8121110913min 9008003 . 12 . 15 . 000000011 . 14 . 0X5006004001001000100100010