用matlab解决插值问题

《用matlab解决插值问题》由会员分享，可在线阅读，更多相关《用matlab解决插值问题（87页珍藏版）》请在金锄头文库上搜索。

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

2、图;实验三：插值和拟和实验三：插值和拟和;实验四：微分、积分和常微分方程实验四：微分、积分和常微分方程;实验五：最优化方法实验五：最优化方法;如何撰写数学建模论文如何撰写数学建模论文.实实 验验 一一 矩阵运算与矩阵运算与Matlab命令命令Matlab基本指令基本指令向量的创建和运算向量的创建和运算1. 直接输入向量直接输入向量x1=1 2 4, x2=1,2,1, x3=x12. 冒号创建向量冒号创建向量 x1=3.4:6.7, x2=3.4:2:6.7, x3=2.6:-0.8:0 3.生成线性等分向量生成线性等分向量指令指令x=linspace(a,b,n) 在在a,b区间产区间产生生

3、 n 个等分点个等分点(包括端点包括端点) 例如：例如：x=linspace(0,1,5)结果结果 x = 0 0.2500 0.5000 0.7500 1.0000工作空间工作空间在在Matlab窗口创建向量后并运行后，窗口创建向量后并运行后，向量就存在于工作空间向量就存在于工作空间（Workspace），可以被调用。），可以被调用。 向量的运算向量的运算 设设三三维维向向量量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=x

4、1+y1 x2+y2 x3+y3向量向量差差： x-y=x1-y1 x2-y2 x3-y3 向量对应元素的运算向量对应元素的运算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=x15 x25 x35 (乘幂乘幂) 2.x=2x1 2x2 2x3x.y=x1y1 x2y2 x3y3函数计算函数计算Matlab有有许许多多内内部部函函数数，可可直直接接作作用用于于向量产生一个同维的函数向量。向量产生一个同维的函数

5、向量。如如：x=linspace(0,4*pi,100);（产产生生100维向量维向量x） y=sin(x); (y也自动为也自动为100维向量维向量) y1=sin(x).2; y2=exp(-x).*sin(x); 观察结果观察结果创建矩阵（数值矩阵的创建）创建矩阵（数值矩阵的创建） 直接输入法创建简单矩阵。直接输入法创建简单矩阵。 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.3000 1.

6、7321 2.4000 -0.9589 7.3891 6.0000矩阵的运算矩阵的运算(矩阵的加减、数乘、乘积等矩阵的加减、数乘、乘积等)A, A_trans=AH=1 2 3 ; 2 1 0 ; 1 2 3 ,K=1 2 3 ; 2 1 0 ; 2 3 1h_det=det(H), k_det=det(K),H_inv=inv(H), K_inv=K-1矩阵的运算矩阵的运算(左除和右除左除和右除)左除左除“ ”： 求矩阵方程求矩阵方程AX=B的解；（的解；（ A 、B的行要保持一致）的行要保持一致） 解为解为 X=AB； 当当A为方阵且可逆时有为方阵且可逆时有X=AB=inv(A)*B；右除

7、右除“ / ”： 求矩阵方程求矩阵方程XA=B的解的解 （A 、B的列要保持一致）的列要保持一致） 解为解为 X=B/A ， 当当A为方阵且可逆时有为方阵且可逆时有X=B/A=B*inv(A)矩阵的运算矩阵的运算(左除和右除左除和右除)例例1：求矩阵方程：求矩阵方程：设设A、B满足关系式：满足关系式：AB2B+A,求求B。其中其中A=3 0 1; 1 1 0; 0 1 4。解解：有：有(A-2I)BA 程序程序 ： A=3 0 1; 1 1 0;0 1 4;B=inv(A-2*eye(3)*A,B=(A-2*eye(3)A 观察结果：观察结果：分块矩阵分块矩阵（矩阵的标识矩阵的标识）例例2 取

8、取出出A的的1、3行行和和1、3列列的的交交叉叉处处元元素构成新矩阵素构成新矩阵A1。解：程序解：程序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) 观察结果观察结果分块矩阵分块矩阵（矩阵的标识矩阵的标识）例例3 将将A分分为为四四块块，并并把把它它们们赋赋值值到到矩矩阵阵B中，观察运行后的结果。中，观察运行后的结果。解：程序解：程序nA11=A(1:2,1:2),A12=A(1:2,3:5),nA21=A(3:4,1:2),A22=A(3:4,3:5)nB=A11 A12;n A21 A22

9、 运行结果运行结果分块矩阵分块矩阵（矩阵的修改和提取）（矩阵的修改和提取）例例4 （1）修改矩阵）修改矩阵A，将它的第将它的第1行变为行变为0。解：程序：解：程序： 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 （2）删除上面矩阵）删除上面矩阵A的第的第1、3行。行。 程序：程序：n A(1,3,:)= 生成特殊矩阵生成特殊矩阵 全全1阵阵 n ones(n), ones(m,n), ones(size(A)全零阵：全零阵：n zeros(n) ，zeros(m,n), zeros(size(A) n常常用于对某

10、个矩阵或向量赋常常用于对某个矩阵或向量赋0初值初值单位阵：单位阵：neye(n)，eye(m,n) 随机阵：随机阵： nrand(m,n)， rand(n)=rand(n,n)用用于于随随机机模模拟，常和拟，常和rand(seed,k)配合使用配合使用。常用矩阵函数常用矩阵函数det(A) ： 方阵的行列式；方阵的行列式；rank(A)： 矩阵的秩；矩阵的秩；eig(A)： 方阵的特征值和特征向量；方阵的特征值和特征向量；trace(A)： 矩阵的迹；矩阵的迹；rref(A)： 初等变换阶梯化矩阵初等变换阶梯化矩阵Asvd(A)： 矩阵奇异值分解。矩阵奇异值分解。分块矩阵分块矩阵（矩阵的标识矩

11、阵的标识）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列交叉处的元素所构成新矩阵。列交叉处的元素所构成新矩阵。数据的简单分析数据的简单分析1.当当数数据据为为行行向向量量或或列列向向量量时时，函函数数对对整整个向量个向量进行计算进行计算.2.当当数数据据为为矩矩

12、阵阵时时，命命令令对对列列进进行行计计算算，即即把把每每一一列列数数据据当当成成同同一一变变量量的的不不同同观观察察值。值。常用的命令：常用的命令： max(求求最最大大)、min(求求最最小小)、mean(求求平平均均值值 )、 sum(求求 和和 )、 std(求求 标标 准准 差差 )、cumsum(求求累累积积和和)、median(求求中中值值)、diff(差差分分)、sort(升升序序排排列列)、sortrows(行行升序排列升序排列)等等等等。数据的简单分析数据的简单分析例例5 观观察察：生生成成一一个个36的的随随机机数数矩矩阵阵，并并将将其其各各列列排排序序、求求各各列列的的最

13、最大大值值与与各各列列元元素之和。素之和。解：程序解：程序sort(使每一列从小到大排序使每一列从小到大排序)nA=rand(3,6),nAsort=sort(A), nAmax=max(A), nAsum=sum(A) 观察结果观察结果实实 验验 二二函数可视化与函数可视化与Matlab作图作图绘制平面曲线绘制平面曲线(plot指令指令) plot(x,y)：n以以x为为横横坐坐标标、y为为纵纵坐标绘制二维图形坐标绘制二维图形nx,y是同维数的向量；是同维数的向量； plot(y)：n相相当当于于x=1,2,length(y)时情形。时情形。函数的可视化函数的可视化问：问：f (x), g

14、(x)是周期函数吗？观察它们的图象。是周期函数吗？观察它们的图象。解：程序解：程序clf, x=linspace(0,8*pi,100);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)例例1 令令绘制平面曲线绘制平面曲线（绘制多个图形（绘制多个图形)1. plot(x,y1;y2;),n x是是横横坐坐标标向向量量，y1;y2;是是由由若若干干函函数数的纵坐标拼成的矩阵；的纵坐标拼成的矩阵；2

15、. plot(x1,y1), hold on, plot(x2,y2), hold off3. plot(x1,y1,x2,y2,) 4. plotyyn两个坐标系，用于绘制不同尺度的函数。两个坐标系，用于绘制不同尺度的函数。绘制平面曲线绘制平面曲线（线型、点形和颜色的控制）（线型、点形和颜色的控制）（线型、点形和颜色的控制）（线型、点形和颜色的控制）plot(x,y,颜色线型点形颜色线型点形)plot(x,y,颜颜色色线线型型点点形形,x,y,颜色线型点形颜色线型点形, )句句柄柄图图形形和和set命命令令改改变变属属性性值值，可可套套用：用：nh=plot(x,y), set(h,属属性性

16、,属属性性值值,属属性性,属属性性值值,)n或或plot(x,y,属属性性,属属性性值值)设设置置图图形形对对象的属性。象的属性。绘制平面曲线绘制平面曲线（属性变量和属性值）（属性变量和属性值）（属性变量和属性值）（属性变量和属性值）n线宽：线宽：LineWidthn点的大小：点的大小： MarkerSizen线型：线型：LineStylen颜色：颜色：Color绘制平面曲线绘制平面曲线（例）（例）（例）（例）程序程序nh=plot(0:0.1:2*pi,sin(0:0.1:2*pi); nset(h,LineWidth,5,color,red); ngrid on观察结果观察结果绘制平面曲线

17、绘制平面曲线（坐标轴的控制）（坐标轴的控制）（坐标轴的控制）（坐标轴的控制）grid on 指令为图形窗口加上网格线指令为图形窗口加上网格线axis指令指令 axis(xmin xmax ymin ymax)： 设定二维图形的设定二维图形的x和和y坐标的范围；坐标的范围； axis(xmin xmax ymin ymax zmin zmax)： 设定三维图形的坐标范围设定三维图形的坐标范围 ； 其其中中xminxxmax， yminyymax ， zminzzmax。绘制平面曲线绘制平面曲线（文字标注）（文字标注）（文字标注）（文字标注）title(图形标题图形标题)；xlabel（x轴轴名名

18、称称）；ylabel（y轴轴名名称称）；）；zlabel（z轴名称轴名称）；）；text(说明文字说明文字)：创建说明文字；：创建说明文字；gtext（说说明明文文字字）：用用鼠鼠标标在在特特定定位位置置输输入入文文字。字。文字标注常用符号：文字标注常用符号： pi （）；）；alpha （）；）；beta （）；）； leftarrow（左箭头）（左箭头） rightarrow（右箭头）；（右箭头）； bullet （点号）点号）图形窗口的创建和分割图形窗口的创建和分割 subplot(m,n,k)命令命令 在在图图形形区区域域中中显显示示多多个个图图形形窗窗口口，m为为上上下下分分割数，割

19、数，n为左右分割数，为左右分割数，k为第为第k子图编号。子图编号。例例：将将一一个个图图形形分分为为4个个子子图图，在在第第k个个子子图图画画sin(kx) 的图象的图象.程序：程序： clf，b=2*pi; x=linspace(0,b,50); for k =1:4 y=sin(k * x); subplot(2,2,k),plot(x,y),axis(0,2*pi,-1,1) end若干有用的指令若干有用的指令clf：清除图形窗口已有的内容清除图形窗口已有的内容. shg：显示图形窗口。显示图形窗口。clear、 clear x：清清除除工工作作空空间间的的已已有变量。有变量。figur

20、e(n): 打开第打开第n个图形窗口个图形窗口 help: : 续行号续行号绘制二元函数绘制二元函数z=f(x,y)基本步骤：基本步骤： 1. 生成二维网格点生成二维网格点 2. 计算函数在网格点上的值计算函数在网格点上的值 3. 绘制函数图形绘制函数图形1. meshgrid指令：生成网格点指令：生成网格点观察观察meshgrid指令的效果。指令的效果。程序：程序： a=-3;b=3;c=-3;d=3;n=10; x=linspace(a,b,n); y=linspace(c,d,n); X,Y=meshgrid(x,y); plot(X,Y,+)观察结果观察结果2. 计算函数值计算函数值

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

22、（等值线指令（等值线指令）contour(X,Y,z,n)：nn条等高线，条等高线，n可缺省；可缺省；contourf(X,Y,z,n)：n等值线间用不同的颜色填满，有更好的视觉效果；等值线间用不同的颜色填满，有更好的视觉效果； contour3(X,Y,z,n)：n在三维空间画出等值线图；在三维空间画出等值线图；n画单个等高线时 contour(cx,cy,cz,I I );colorbar：n将颜色与函数值对应起来显示在图中将颜色与函数值对应起来显示在图中。空间曲线和运动方向的表现空间曲线和运动方向的表现一条空间曲线可以用矢量函数表示为一条空间曲线可以用矢量函数表示为n它的速度矢量表现为曲

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

24、现空间曲线和运动方向的表现n很显然飞行曲线方程为： 程序：程序：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(

25、z)， box on, hold off 插值和拟和插值和拟和 实实 验验 三三机床加工问题机床加工问题用程控铣床加工机翼断面的下轮廓线时用程控铣床加工机翼断面的下轮廓线时每一刀只能沿每一刀只能沿x方向和方向和y方向走非常小的一步。方向走非常小的一步。表表3-1给出了下轮廓线上的部分数据给出了下轮廓线上的部分数据但工艺要求铣床沿但工艺要求铣床沿x方向每次只能移动方向每次只能移动0.1单位单位. 这时需求出当这时需求出当x坐标每改变坐标每改变0.1单位时的单位时的y坐标。坐标。试完成加工所需的数据，画出曲线试完成加工所需的数据，画出曲线.插值插值n已知有已知有n +1个节点（个节点（xj,yj）

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

27、次样条插值；值；cubic ： 立方插值。立方插值。缺省时：缺省时： 线性插值。线性插值。V4:V4: 注意：所有的插值方法都要求注意：所有的插值方法都要求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,k+,x,y,r)grid on

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

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

30、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, y）为测量点，为测量点，z为（为（x, y）处的水深处的水深(英尺英尺)，水深，水深z是区域坐标（是区域坐标（x,

31、y）的函数）的函数z= z (x, y），），船的吨位可以用其吃水深度来反映，分为船的吨位可以用其吃水深度来反映，分为 4英英尺、尺、4.5英尺、英尺、5英尺和英尺和 5.5英尺英尺 4 档。档。 航运部门要在矩形海域（航运部门要在矩形海域（75，200）（50，150）上为不同吨位的航船设置警示标记。）上为不同吨位的航船设置警示标记。请根据测量的数据描述该海域的地貌，并绘制请根据测量的数据描述该海域的地貌，并绘制不同吨位的警示线，供航运部门使用。不同吨位的警示线，供航运部门使用。Matlab求解求解x=129 140 103.5 88 185.5 195 105.5 157.5 107.5

32、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); grid on,hold onplot(x,y,+)xla

33、bel(X),ylabel(Y)画单个等高线时 contour(cx,cy,cz,I I ); 曲线拟合曲线拟合 已知一组（二维）数据，即平面上已知一组（二维）数据，即平面上 n个点个点（xi,yi) i=1,n, 寻求一个函数（曲线）寻求一个函数（曲线）y=f(x), 使使 f(x) 在在某种准则某种准则下与所有数据点最为接近，即曲线拟合得最好。下与所有数据点最为接近，即曲线拟合得最好。 +xyy=f(x)(xi,yi)i i 为点为点（xi,yi) 与与曲线曲线 y=f(x) 的的距距离离最常用的方法是线性最小二乘拟和最常用的方法是线性最小二乘拟和多项式拟合多项式拟合n对给定的数据（对给定

34、的数据（xj,yj），），j = 0,1, n；n选取适当阶数的多项式，如二次多项式选取适当阶数的多项式，如二次多项式g(x)=ax2+bx+c；n使使g(x)尽可能逼近（拟合）这些数据，但尽可能逼近（拟合）这些数据，但是不要求经过给定的数据（是不要求经过给定的数据（xj,yj）；）；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）输出拟合多项式系数输出拟合多项式系数a=a1,am,am+1 （数组）数组

35、）输入同长度输入同长度数组数组X，Y拟合多项式拟合多项式次数次数多项式拟合指令多项式拟合指令即要求出二次多项式即要求出二次多项式:中中 的的使得使得:例例 对下面一组数据作二次多项式拟合对下面一组数据作二次多项式拟合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解：解：用多项式拟合的命令用多项式拟合的命令1）输入命令：）输入命令：x=0:0.1:1; y=-0.447,1.978,3.28,6.16,

36、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) %作出数据点和拟合曲线的图形作出数据点和拟合曲线的图形可化为多项式的非线性拟和可化为多项式的非线性拟和一般非线性最小二乘拟和一般非线性最小二乘拟和,实际上是无约实际上是无约束最优化问题束最优化问题命令：命令：lsqcurvefit、lsqnonlin等等 微微分分、积积分分和和常常微微分分方程方程实验四实验四数值积分常用的方法数值积分常用的方法（1）用不定积分计算定积分：）用不定积分计算定积分： 牛牛顿顿莱莱布布尼尼兹兹公公式

37、式,但但是是有有时时得得不不到原函数；到原函数；（2）定义法，取近似和的极限：）定义法，取近似和的极限：n高等数学中不是重点内容高等数学中不是重点内容n但数值积分的各种算法却是基于定义建立的但数值积分的各种算法却是基于定义建立的 数值微积分数值微积分（梯形公式和辛普森公式梯形公式和辛普森公式）trapz(x,y)，按梯形公式计算近似积分；，按梯形公式计算近似积分； 其中步长x=x0 x1 xn和函数值y=f0 f1 fn为同维向量；q = quad(fun,a,b,tol,trace,P1,P2,.) （低阶方法，辛普森自适应递归法求积分）低阶方法，辛普森自适应递归法求积分）q = quad8

38、(fun,a,b,tol,trace,P1,P2,.)（高高阶方法，自适应法阶方法，自适应法Cotes求积分）求积分） 在在同同样样的的精精度度下下高高阶阶方方法法quad8要要求求的的节节点点较较少。少。例：数值积分例：数值积分X = 0:pi/100:pi;Y = sin(X);Z1 = trapz(X,Y)Z2=quad(sin,0,pi)Z1 = 1.998z2 = 2.0000求解常微分方程求解常微分方程一阶常微分方程数值解法一阶常微分方程数值解法x,y=ode23(fun,tspan,y0,option) （低阶龙格库塔函数）（低阶龙格库塔函数）x,y=ode45(fun,tspa

39、n,y0,option) （高阶龙格库塔函数）（高阶龙格库塔函数）其中，(1) tspan=t0，tf，t0、tf为自变量的初值和终值 (2) option用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令为：options=odeset （reltol,rt,abstol,at）, rt，at：分别为设定的相对误差和绝对误差. 1、在解、在解n个未知函数的方程组时，个未知函数的方程组时，x0和和x均为均为n维向量，维向量，m-文件中的待解方程文件中的待解方程组应以组应以x的分量形式写成的分量形式写成. 2、使用、使用Matlab软件求数值解时，高软件求数值解时，高阶微分

40、方程必须等价地变换成一阶微分阶微分方程必须等价地变换成一阶微分方程组方程组.注意注意: : 设位于坐标原点的甲舰向位于设位于坐标原点的甲舰向位于x轴上点轴上点A(1, 0)处的乙舰发射导弹，导弹头始终对准处的乙舰发射导弹，导弹头始终对准乙舰。如果乙舰以最大的速度乙舰。如果乙舰以最大的速度v0(是常数是常数)沿沿平行于平行于y轴的直线行驶，轴的直线行驶，导弹的速度是导弹的速度是5v0，（1）画出导弹运行的曲线方程）画出导弹运行的曲线方程.（2）乙舰行驶多远时，）乙舰行驶多远时， 导弹将它击中？导弹将它击中？导弹追击问题由由(1),(2)消去消去t整理得模型整理得模型:1.建立建立m-文件文件eq

41、1.m function dy=eq1(x,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1/5*sqrt(1+y(1)2)/(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）处击中乙舰）处击中乙舰令y1=y, y2=y1，将方程（3）化为一阶微分方程组。 最优化方法实验五最优

42、化方法最优化方法 许多生产计划与管理问题都可以归纳为最优许多生产计划与管理问题都可以归纳为最优化问题化问题, , 最优化模型是数学建模中应用最广泛的最优化模型是数学建模中应用最广泛的模型之一模型之一, ,其内容包括线性规划、整数线性规划、其内容包括线性规划、整数线性规划、非线性规划、动态规划、变分法、最优控制等非线性规划、动态规划、变分法、最优控制等. . 近几年来的全国大学生数学建模竞赛中，几近几年来的全国大学生数学建模竞赛中，几乎每次都有一道题要用到此方法乎每次都有一道题要用到此方法. .目标函数 约束条件 可行解域 线性规划线性规划用用MATLAB优化工具箱解线性规划优化工具箱解线性规划

43、min z=cX 1、模型：命令：x=linprog（c，A，b） 2、模型：min z=cX 命令：x=linprog（c，A，b，Aeq, beq）注意：若没有不等式： 存在，则令A= ，b= .3、模型：min z=cX VLBXVUB命令：1 x=linprog（c，A，b，Aeq, beq, VLB，VUB） 2 x=linprog（c，A，b，Aeq, beq, VLB，VUB, X0） 注意：若没有等式约束: , 则令Aeq= , beq= ; 其中X0表示初始点 4、命令：x,fval=linprog()返回最优解及处的目标函数值返回最优解及处的目标函数值fval.解解: 编写

44、编写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) 某车间有甲、乙两台机床，可用于加工三种工件。假定这两台车床的可用台时数分别为800和900，三种工件的数量分别为400、600和500，且已知用三种不同车床加工单位数量不同工件所需的台时数和加工费用如下表。问怎样分配车床的加工任务，才能既满足加工工件的要求，又使加工费用最低？ 例例（任务分配问题（任务分配问题)解解 设在甲车床上加工工件1、2、3的数量分

45、别为x1、x2、x3，在乙车床上加工工件1、2、3的数量分别为x4、x5、x6。可建立以下线性规划模型：S.t.改写为：编写编写M文件如下文件如下: f = 13 9 10 11 12 8; A = 0.4 1.1 1 0 0 0; 0 0 0 0.5 1.2 1.3; b = 800; 900; Aeq=1 0 0 1 0 0; 0 1 0 0 1 0; 0 0 1 0 0 1; beq=400 600 500; vlb = zeros(6,1); vub=; x,fval = linprog(f,A,b,Aeq,beq,vlb,vub)结果结果: x = 0.0000 600.0000 0

46、.0000 400.0000 0.0000 500.0000 fval =1.3800e+004 即在甲机床上加工即在甲机床上加工600个工件个工件2,在乙机床上在乙机床上加工加工400个工件个工件1、500个工件个工件3，可在满足条件，可在满足条件的情况下使总加工费最小为的情况下使总加工费最小为13800。 1. 首先建立首先建立M文件文件fun.m,定义目标函数定义目标函数F（X）:function f=fun(X);f=F(X);非线性规划非线性规划 其中其中X为为n维变元向量，维变元向量，G(X)与与Ceq(X)均为非线性函数组均为非线性函数组 成的向量成的向量. 用用Matlab求解

47、上述问题，基本步骤分三步：求解上述问题，基本步骤分三步： 3. 建立主程序建立主程序.非线性规划求解的函数是非线性规划求解的函数是fmincon,命令的基命令的基本格式如下：本格式如下： (1) x=fmincon(fun,X0,A,b) (2) x=fmincon(fun,X0,A,b,Aeq,beq) (3) x=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB) (4) x=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB,nonlcon)(5)x=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB,nonlcon,options

48、) (6) x,fval= fmincon(.) 输出极值点M文件迭代的初值参数说明变量上下限注意：注意：fmincon函数可能会给出局部最优解，函数可能会给出局部最优解， 这与初值这与初值X0的选取有关。的选取有关。 1、写成标准形式写成标准形式： s.t. 2x1+3x2 6 s.t x1+4x2 5 x1,x2 0例例32、先建立M-文件 fun3.m: function f=fun3(x); f=-x(1)-2*x(2)+(1/2)*x(1)2+(1/2)*x(2)23、再建立主程序youh2.m： x0=1;1; A=2 3 ;1 4; b=6;5; Aeq=;beq=; VLB=0

49、;0; VUB=; x,fval=fmincon(fun3,x0,A,b,Aeq,beq,VLB,VUB) 4、运算结果为：运算结果为： x = 0.7647 1.0588 fval = -2.0294如何撰写数学建模论文如何撰写数学建模论文什么是数学建模什么是数学建模? 数学建模数学建模是利用数学方法解决实际问题的一种实践。即通过抽象、简化、假设、引进变量等处理过程后，将实际将实际问题用数学方式表达问题用数学方式表达，建立起数学模型，然后运用先进的数学方法及计算机技术运用先进的数学方法及计算机技术进行求解进行求解。 在实际过程中用那一种方法建模主要是根据我们对研究对象的了解程度和建模目的来决定。机理分析法建模的具体步骤大致可见右图。符合实际不符合实际交付使用，从而可产生经济、社会效益实际问题抽象、简化、假设 确定变量、参数建立数学模型并数学、数值地求解、确定参数用实际问题的实测数据等来检验该数学模型建模过程示意图 怎样撰写数学建模的论文？怎样撰写数学建模的论文？1、摘要、摘要: 问题、模型、方法、结问题、模型、方法、结果果2、问题重述、问题重述4、分析与建立模型、分析与建立模型5、模型求解、模型求解 (借助数学软件借助数学软件)6、模型检验、模型检验7、模型推广、模型推广8、参考文献、参考文献9、附录、附录3、模型假设、模型假设