用matlab解决插值问题PPT学习教案

1、会计学1用用matlab解决插值问题解决插值问题 MATLAB是一种交互式的以矩阵为基础的系统计算平台,它用于科学和工程的计算与可视化。 Matlab的含义是矩阵实验室(Matrix Laboratory),是美国Math Work公司于1982年推出的一套高性能的数值计算和可视化软件,它集数值分析、矩阵计算、信号数值分析、矩阵计算、信号处理和图形显示处理和图形显示于一体,已发展成为国际上最优秀的科技应用软件之一。 矩阵矩阵是是MATLAB的核心的核心. 第1页/共87页第2页/共87页第3页/共87页第4页/共87页 x1=3.4 x1=3.4: :6.7,6.7, x2=3.4x2=3.4

2、:2:2:6.7,6.7, x3=2.6x3=2.6:-0.8:-0.8:0 0第5页/共87页第6页/共87页第7页/共87页第8页/共87页nx x . . y = x 1 y 1 x 2 y 2 y = x 1 y 1 x 2 y 2 x3y3x3y3第9页/共87页n观察结果观察结果第10页/共87页第11页/共87页第12页/共87页第13页/共87页第14页/共87页第15页/共87页第16页/共87页n第17页/共87页第18页/共87页第19页/共87页第20页/共87页第21页/共87页观察结果观察结果第22页/共87页第23页/共87页第24页/共87页)sin(cos(s

3、in()(xxxxf)sin(cos(sin(2 . 0)(xxxxxg问：问：f (x), g (x)f (x), g (x)是周期函数吗？观察它们的图象。是周期函数吗？观察它们的图象。解：程序解：程序clf, x=linspace(0,8clf, 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(x);x+sin(x+cos(x+sin(x);plot(x,y1,k:,x,y2,k-) plot(x,y1,k:,x,y

4、2,k-) legend(sin(x+cos(x+sin(x),0.2x+sin(legend(sin(x+cos(x+sin(x),0.2x+sin(x+cos(x+sin(x),2)x+cos(x+sin(x),2)例例1 1 令令第25页/共87页度的函数。度的函数。第26页/共87页n或或plot(x,y,plot(x,y,属性属性,属性值属性值) )设置图形对象的属性。设置图形对象的属性。第27页/共87页n线宽：线宽：LineWidthLineWidthn点的大小：点的大小： MarkerSizeMarkerSizen线型：线型：LineStyleLineStylen颜色：颜色：C

5、olorColor第28页/共87页第29页/共87页第30页/共87页第31页/共87页第32页/共87页第33页/共87页第34页/共87页第35页/共87页mesh(X,Y,z) mesh(X,Y,z) ：n在三维空间中绘出由在三维空间中绘出由(X,Y,z)(X,Y,z)表示的曲面表示的曲面；meshz(X,Y,z)meshz(X,Y,z)：n除了具有除了具有meshmesh的功能外，还画出上下高度线，的功能外，还画出上下高度线，meshc(X,Y,z)meshc(X,Y,z)：n除了具有除了具有meshmesh的功能外，还在曲面的下方画出函数的功能外，还在曲面的下方画出函数z=f(x,

6、y)z=f(x,y)的等值线图，的等值线图，surf(X,Y,z)surf(X,Y,z)：n也是三维绘图指令，与也是三维绘图指令，与meshmesh的区别在于的区别在于meshmesh绘出彩色的线，绘出彩色的线，surfsurf绘出彩色的面，绘出彩色的面，第36页/共87页第37页/共87页, , )(),(),()( :ttztytxtrn它的速度矢量表现为曲线的切矢量:)(,)(,)()(dttdzdttdydttdxdtdrtv第38页/共87页第39页/共87页)3/1 (3 ,)3/2( ,()6 ,2 ,2()(34320232ttttdttttttrtn很显然飞行曲线方程为： 程

7、序：程序：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, hold off第40页/共8

8、7页第41页/共87页 表表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.6第42页/共87页近似。近似。第43页/共87页yi=interp1(x，y，xi，method)插值方法插值方法被插值点被插值点插值节点插值节点xixi处的插处的插值结果值结果nearest ：最邻近插值：最邻近插值linear ： 线性插值；线性插值；spline ： 三次样条插值三次样条插值；cubic ： 立方插值。立方插值。缺省时：缺省时： 线性插值。线性插值。V4:V4: 注

9、意：所有的插值方法都要求注意：所有的插值方法都要求x x是单调的，并且是单调的，并且xi不不能够超过能够超过x的范围。的范围。用用MATLABMATLAB作一维插值计算作一维插值计算第44页/共87页第45页/共87页 要求要求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)被插值点插值方法插值节点被插值点的函数值nearestnearest 最邻近插值最邻近插值line

10、arlinear 双线性插值双线性插值cubiccubic 双三次插值双三次插值缺省时缺省时, , 双线性插值双线性插值用用MATLABMATLAB作网格节点数据的插值作网格节点数据的插值( (二维二维) )第46页/共87页 cz =griddata（x，y，z，cx，cy，method） 要求要求cxcx取行向量，取行向量，cycy取为列向量取为列向量。被插值点插值方法插值节点被插值点的函数值nearestnearest 最邻近插值最邻近插值linearlinear 双线性插值双线性插值cubiccubic 双三次插值双三次插值v4 Matlab提供的插值方法提供的插值方法缺省时缺省时,

11、, 双线性插值双线性插值用用MATLABMATLAB作散点数据的插值计算作散点数据的插值计算第47页/共87页第48页/共87页第49页/共87页第50页/共87页 已知一组（二维）数据，即平面上已知一组（二维）数据，即平面上 n个点个点（xi,yi) i=1,n, 寻求一个函数（曲线）寻求一个函数（曲线）y=f(x), 使使 f(x) 在某种准则下与所有数据点最为接近，即曲线拟合得最好。在某种准则下与所有数据点最为接近，即曲线拟合得最好。 +xyy=f(x)(xi,yi)i i 为点为点（xi,yi) 与与曲线曲线 y=f(x) 的距离的距离最常用的方法是线性最小二乘拟和最常用的方法是线性最

12、小二乘拟和第51页/共87页第52页/共87页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 （数组）数组）输入同长度输入同长度数组数组X，Y拟合多项式拟合多项式次数次数多项式拟合指令多项式拟合指令第53页/共87页即要求出二次多项式即要求出二次多项式:3221)(axaxaxf中中 的的123(,)Aa a a 使得使得:1121() iii

13、f xy最小例例 对下面一组数据作二次多项式拟合对下面一组数据作二次多项式拟合xi00.10.20.30.40.50.60.70.80.91yi-0.4471.9783.286.166.167.347.669.589.489.3011.2第54页/共87页2）计算结果：）计算结果： = -9.8108, 20.1293, -0.0317解：解：用多项式拟合的命令用多项式拟合的命令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

14、,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) %作出数据点和拟合曲线的图形作出数据点和拟合曲线的图形第55页/共87页一般非线性最小二乘拟和一般非线性最小二乘拟和,实际上是无实际上是无约束最优化问题约束最优化问题命令：命令：lsqcurvefit、lsqnonlin等等第56页/共87页第57页/共87页)x(x)f(limf(x)dxi1i1n0ii0ba第58页/共87页第59页/共87页第60页/共87页第61页/共87页第62页/共87页 1、在解、在解n个未知函数的方程组时，个未

15、知函数的方程组时，x0和和x均为均为n维向量，维向量，m-文件中的待解方程组文件中的待解方程组应以应以x的分量形式写成的分量形式写成. 2、使用、使用Matlab软件求数值解时，高软件求数值解时，高阶微分方程必须等价地变换成一阶微分阶微分方程必须等价地变换成一阶微分方程组方程组.注意注意: :第63页/共87页 设位于坐标原点的甲舰向位于设位于坐标原点的甲舰向位于x轴上点轴上点A(1, 0)处的乙舰发射导弹，导弹头始终对准乙舰。处的乙舰发射导弹，导弹头始终对准乙舰。如果乙舰以最大的速度如果乙舰以最大的速度v0(是常数是常数)沿平行于沿平行于y轴的直线行驶，轴的直线行驶，导弹的速度是导弹的速度是

16、5v0，（1）画出导弹运行的曲线方程）画出导弹运行的曲线方程.（2）乙舰行驶多远时，）乙舰行驶多远时， 导弹将它击中？导弹将它击中？导弹追击问题第64页/共87页由由(1),(2)消去消去t整理得模型整理得模型:(3) 151)1 (2yyx第65页/共87页1.建立建立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)/(1-x); 2. 取取x0=0，xf=0.9999，建立主程序，建立主程序ff6.m如下如下： x0=0，xf=0.9999 x,y=ode23(eq1,x0

17、 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）化为一阶微分方程组。第66页/共87页第67页/共87页 许多生产计划与管理问题都可以归纳为最优许多生产计划与管理问题都可以归纳为最优化问题化问题, , 最优化模型是数学建模中应用最广泛的最优化模型是数学建模中应用最广泛的模型之一模型之一, ,其内容包括线性规划、整数线性规划其内容包括线性规划

18、、整数线性规划、非线性规划、动态规划、变分法、最优控制等、非线性规划、动态规划、变分法、最优控制等. . 近几年来的全国大学生数学建模竞赛中，几近几年来的全国大学生数学建模竞赛中，几乎每次都有一道题要用到此方法乎每次都有一道题要用到此方法. .目标函数 约束条件 可行解域 xpjxhmixgtsxfji, 2 , 10)(, 2 , 1, 0)(. .)(min第68页/共87页xpjxhmixgtsxfji, 2 , 10)(, 2 , 1, 0)(. .)(min若)() ,() ,(xhxgxfii都是线性 函 数 ， 则 称 为 线 性 规 划 问题 ， 简 称L P问 题 第69页/

19、共87页用用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 第70页/共87页3、模型：min z=cX bAXts. .beqXAeq VLBXVUB命令：1 x=linprog（c，A，b，Aeq, beq, VLB，VUB） 2 x=linprog（c，A，b，Aeq, beq, VLB，VUB, X0） 注意：若没有等式约束: ,

20、则令Aeq= , beq= ; 其中X0表示初始点 beqXAeq4、命令：x,fval=linprog()返回最优解及处的目标函数值返回最优解及处的目标函数值fval.第71页/共87页解解: 编写编写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,120111 .321321xxxxxxts第72页/共87页 某车间有甲、乙两台机床，可用

21、于加工三种工件。假定这两台车床的可用台时数分别为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 例例（任务分配问题（任务分配问题)第73页/共87页解解 设在甲车床上加工

22、工件1、2、3的数量分别为x1、x2、x3，在乙车床上加工工件1、2、3的数量分别为x4、x5、x6。可建立以下线性规划模型：6543218121110913minxxxxxxz 6 , 2 , 1, 09003 . 12 . 15 . 08001 . 14 . 0500600400 x . .654321635241ixxxxxxxxxxxxtsi第74页/共87页S.t.Xz8121110913min 9008003 . 12 . 15 . 000000011 . 14 . 0X500600400100100010010001001X ，0654321xxxxxxX改写为：第75页/共87

23、页编写编写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)第76页/共87页结果结果: x = 0.0000 600.0000 0.0000 400.0000 0.0000 500.0000 fval =1.3800e+004 即

24、在甲机床上加工即在甲机床上加工600个工件个工件2,在乙机床上在乙机床上加工加工400个工件个工件1、500个工件个工件3，可在满足条件，可在满足条件的情况下使总加工费最小为的情况下使总加工费最小为13800。第77页/共87页 1. 首先建立首先建立M文件文件fun.m,定义目标函数定义目标函数F（X）:function f=fun(X);f=F(X);非线性规划非线性规划 其中其中X为为n维变元向量，维变元向量，G(X)与与Ceq(X)均为非线性函数组均为非线性函数组 成的向量成的向量. 用用Matlab求解上述问题，基本步骤分三步：求解上述问题，基本步骤分三步：第78页/共87页 3.

25、建立主程序建立主程序.非线性规划求解的函数是非线性规划求解的函数是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) (6) x,fval= fmincon(.) 输出极值点M文件迭代的初值参数说明变量上下限注意：注意：fmincon函数可能会给出局部最优解，函数可能会给出局部最优解， 这与初值这与初值X0的选取有关。的选取有关。第79页/共87页 1、写成标准形式写成标准形式： s.t. 005