《数学模型》实验指导书通识版Word格式.docx

1、 b=regress（y,X） b,bint,r,rint,stats= regress（y,X） b,bint,r,rint,stats= regress（y,X,alpha）说明：b=regress（y,X）返回X处y的最小二乘拟合值。该函数求解线性模型：y=X+是p1的参数向量；是服从标准正态分布的随机干扰的n1的向量；y为n1的向量；X为np矩阵。bint返回的95%的置信区间。r中为形状残差，rint中返回每一个残差的95%置信区间。Stats向量包含R2统计量、回归的F值和p值。例1：设y的值为给定的x的线性函数加服从标准正态分布的随机干扰值得到。即y=10+x+ ；求线性拟合方程

2、系数。程序： x=ones（10,1） （1:10） y=x*10;1+normrnd（0,0.1,10,1） b,bint=regress（y,x,0.05）结果： x = 1 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 1 9 1 10y = 10.9567 11.8334 13.0125 14.0288 14.8854 16.1191 17.1189 17.9962 19.0327 20.0175b = 9.9213 1.0143bint = 9.7889 10.0537 0.9930 1.0357即回归方程为：y=9.9213+1.0143x1.2多项式曲线拟合函数：p

3、olyfit（ ） p=polyfit（x,y,n） p,s= polyfit（x,y,n）x,y为数据点，n为多项式阶数，返回p为幂次从高到低的多项式系数向量p。矩阵s用于生成预测值的误差估计。（见下一函数polyval）例2：由离散数据x.1.2.3.4.5.6.7.8.9y1.41.61.91.5拟合出多项式。 x=0:.1:1; y=.3 .5 1 1.4 1.6 1.9 .6 .4 .8 1.5 2 n=3; p=polyfit（x,y,n） xi=linspace（0,1,100）; z=polyval（p,xi）; %多项式求值 plot（x,y,o,xi,z,k:,x,y,b）

4、 legend（原始数据,3阶曲线）p = 16.7832 -25.7459 10.9802 -0.0035多项式为：16.7832x3-25.7459x2+10.9802x-0.0035曲线拟合图形：也可由函数给出数据。例3：x=1:20,y=x+3*sin（x） x=1:20; y=x+3*sin（x）; p=polyfit（x,y,6） xi=1inspace（1,20,100）; z=polyval（p,xi）; %多项式求值函数 plot（x,y,o,xi,z,k: legend（原始数据,6阶曲线）0.0000 -0.0021 0.0505 -0.5971 3.6472 -9.72

5、95 11.3304再用10阶多项式拟合 程序：y=x+3*sin（x）;p=polyfit（x,y,10）xi=linspace（1,20,100）;z=polyval（p,xi）;plot（x,y,o,xi,z,k:,x,y,b）legend（原始数据,10阶多项式 Columns 1 through 7 0.0000 -0.0000 0.0004 -0.0114 0.1814 -1.8065 11.2360 Columns 8 through 11 -42.0861 88.5907 -92.8155 40.2671可用不同阶的多项式来拟合数据，但也不是阶数越高拟合的越好。1.3 多项式曲

6、线求值函数：polyval（ ） y=polyval（p,x） y,DELTA=polyval（p,x,s）y=polyval（p,x）为返回对应自变量x在给定系数P的多项式的值。y,DELTA=polyval（p,x,s） 使用polyfit函数的选项输出s得出误差估计YDELTA。它假设polyfit函数数据输入的误差是独立正态的，并且方差为常数。则YDELTA将至少包含50%的预测值。1.4 多项式曲线拟合的评价和置信区间函数：polyconf（ ） Y,DELTA=polyconf（p,x,s） Y,DELTA=polyconf（p,x,s,alpha）Y,DELTA=polyconf

7、（p,x,s）使用polyfit函数的选项输出s给出Y的95%置信区间Y1-alpha为置信度。例4：给出上面例1的预测值及置信度为90%的置信区间。 x=0: y=.3 .5 1 1.4 1.6 1.9 .6 .4 .8 1.5 2 n=3; p,s=polyfit（x,y,n） alpha=0.05;Y,DELTA=polyconf（p,x,s,alpha） 结果： p = 16.7832 -25.7459 10.9802 -0.0035s = R: 4x4 double df: 7 normr: 1.1406Y = -0.0035 0.8538 1.2970 1.4266 1.3434

8、1.1480 0.9413 0.8238 0.8963 1.2594 2.0140DELTA = 1.3639 1.1563 1.1563 1.1589 1.1352 1.1202 1.13521.1589 1.1563 1.1563 1.36391.5.稳健回归函数：robust（ ）稳健回归是指此回归方法相对于其他回归方法而言，受异常值的影响较小。 b=robustfit（x,y） b,stats=robustfit（x,y） b,stats=robustfit（x,y,wfun,tune,const）b返回系数估计向量；stats返回各种参数估计；wfun指定一个加权函数；tune为调协

9、常数；const的值为on（默认值）时添加一个常数项；为off 时忽略常数项。例5：演示一个异常数据点如何影响最小二乘拟合值与稳健拟合。首先利用函数y=10-2x加上一些随机干扰的项生成数据集，然后改变一个y的值形成异常值。调用不同的拟合函数，通过图形观查影响程度。x=（1:10）;y=10-2*x+randn（10,1）;y（10）=0;bls=regress（y,ones（10,1） x） %线性拟合brob=robustfit（x,y） %稳健拟合scatter（x,y）hold onplot（x,bls（1）+bls（2）*x,:）plot（x,brob（1）+brob（2）*x,r）

10、结果 ： bls = 8.4452 -1.4784brob = 10.2934 -2.0006分析：稳健拟合（实线）对数据的拟合程度好些，忽略了异常值。最小二乘拟合（点线）则受到异常值的影响，向异常值偏移。1.6 自定义函数拟合对于给定的数据，根据经验拟合为带有待定常数的自定义函数。所用函数：nlinfit（ ） beta,r,J=nlinfit（X,y,fun,betao）beta返回函数fun中的待定常数；r表示残差；J表示雅可比矩阵。X,y为数据；fun自定义函数；beta0待定常数初值。例6：在化工生产中获得的氯气的级分y随生产时间x下降，假定在x8时，y与x之间有如下形式的非线性模型

11、： 现收集了44组数据，利用该数据通过拟合确定非线性模型中的待定常数。x y x y x y8 0.49 16 0.43 28 0.418 0.49 18 0.46 28 0.4010 0.48 18 0.45 30 0.4010 0.47 20 0.42 30 0.4010 0.48 20 0.42 30 0.3810 0.47 20 0.43 32 0.4112 0.46 20 0.41 32 0.4012 0.46 22 0.41 34 0.4012 0.45 22 0.40 36 0.4112 0.43 24 0.42 36 0.3614 0.45 24 0.40 38 0.4014

12、0.43 24 0.40 38 0.4014 0.43 26 0.41 40 0.3616 0.44 26 0.40 42 0.3916 0.43 26 0.41 首先定义非线性函数的m文件：model.mfunction yy=model（beta0,x） a=beta0（1）; b=beta0（2）; yy=a+（0.49-a）*exp（-b*（x-8）;x=8.00 8.00 10.00 10.00 10.00 10.00 12.00 12.00 12.00 14.00 14.00 14.00. 16.00 16.00 16.00 18.00 18.00 20.00 20.00 20.0

13、0 20.00 22.00 22.00 24.00. 24.00 24.00 26.00 26.00 26.00 28.00 28.00 30.00 30.00 30.00 32.00 32.00. 34.00 36.00 36.00 38.00 38.00 40.00 42.00; y=0.49 0.49 0.48 0.47 0.48 0.47 0.46 0.46 0.45 0.43 0.45 0.43 0.43 0.44 0.43. 0.43 0.46 0.42 0.42 0.43 0.41 0.41 0.40 0.42 0.40 0.40 0.41 0.40 0.41 0.41. 0.4

14、0 0.40 0.40 0.38 0.41 0.40 0.40 0.41 0.38 0.40 0.40 0.39 0.39 beta0=0.30 0.02;betafit = nlinfit（x,y,model,beta0）betafit = 0.38960.1011 即：a=0.3896 ，b=0.1011 拟合函数为：。2 插值问题在应用领域中，由有限个已知数据点，构造一个解析表达式，由此计算数据点之间的函数值，称之为插值。2.1 一元插值一元插值是对一元数据点（xi,yi）进行插值。线性插值：由已知数据点连成一条折线，认为相临两个数据点之间的函数值就在这两点之间的连线上。一般来说，数据点

15、数越多，线性插值就越精确。yi=interp1（x,y,xi,linear） %线性插值zi=interp1（x,y,xi,spline） %三次样条插值wi=interp1（x,y,xi,cubic） %三次多项式插值yi、zi、wi为对应xi的不同类型的插值。x、y为已知数据点。已知数据：求当xi=0.25时的yi的值。x=0:y=.3 .5 1 1.4 1.6 1 .6 .4 .8 1.5 2;yi0=interp1（x,y,0.025,linearxi=0:.02:yi=interp1（x,y,xi,）;zi=interp1（x,y,xi,splinewi=interp1（x,y,xi

16、,cubic,xi,yi,r+,xi,zi,g*,xi,wi,k.-原始点线性点三次样条三次多项式yi0 = 0.3500要得到给定的几个点的对应函数值，可用：xi = 0.2500 0.3500 0.4500yi =1.2088 1.5802 1.3454 2.2二元插值二元插值与一元插值的基本思想一致，对原始数据点（x,y,z）构造见世面函数求出插值点数据（xi,yi,zi）。单调节点插值函数，即x,y向量是单调的。调用格式1：zi=interp2（x,y,z,xi,yi,linear） liner 是双线性插值 （缺省）调用格式2：zi=interp2（x,y,z,xi,yi,neare

17、st） nearest 是最近邻域插值 调用格式3：zi=interp2（x,y,z,xi,yi,spline） spline是三次样条插值这里x和y是两个独立的向量，它们必须是单调的。z是矩阵，是由x和y确定的点上的值。z和x,y之间的关系是z（i,:）=f（x,y（i） z（:,j）=f（x（j）,y） 即：当x变化时，z的第i行与y的第i个元素相关，当y变化时z的第j列与x的第j个元素相关。如果没有对x,y赋值，则默认x=1:n, y=1:m。n和m分别是矩阵z的行数和列数。已知某处山区地形选点测量坐标数据为：x=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5y=0 0.

18、5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6海拔高度数据为：z=89 90 87 85 92 91 96 93 90 87 82 92 96 98 99 95 91 89 86 84 82 84 96 98 95 92 90 88 85 84 83 81 85 80 81 82 89 95 96 93 92 89 86 86 82 85 87 98 99 96 97 88 85 82 83 82 85 89 94 95 93 92 91 86 84 88 88 92 93 94 95 89 87 86 83 81 92 92 96 97 98 96 93 95 84 82

19、 81 84 85 85 81 82 80 80 81 85 90 93 95 84 86 81 98 99 98 97 96 95 84 87 80 81 85 82 83 84 87 90 95 86 88 80 82 81 84 85 86 83 82 81 80 82 87 88 89 98 99 97 96 98 94 92 87其地貌图为：对数据插值加密形成地貌图。.5:5;y=0:6;z=89 90 87 85 92 91 96 93 90 87 82 87 88 89 98 99 97 96 98 94 92 87;mesh（x,y,z） %绘原始数据图xi=linspace（

20、0,5,50）; %加密横坐标数据到50个yi=linspace（0,6,80）; %加密纵坐标数据到60个xii,yii=meshgrid（xi,yi）; %生成网格数据zii=interp2（x,y,z,xii,yii, %插值mesh（xii,yii,zii） %加密后的地貌图hold on % 保持图形xx,yy=meshgrid（x,y）;plot3（xx,yy,z+0.1,ob） %原始数据用O绘出2.3二元非等距插值zi=griddata（x,y,z,xi,yi,指定插值方法）插值方法有： linear % 线性插值 （默认） bilinear % 双线性插值 cubic % 三

21、次插值 bicubic % 双三次插值 nearest % 最近邻域插值例：用随机数据生成地貌图再进行插值x=rand（100,1）*4-2;y=rand（100,1）*4-2;z=x.*exp（-x.2-y.2）;ti=-2:.25:2;xi,yi=meshgrid（ti,ti）; % 加密数据zi=griddata（x,y,z,xi,yi）;% 线性插值mesh（xi,yi,zi）plot3（x,y,z,该例中使用的数据是随机形成的，故函数griddata可以处理无规则的数据。实验三 数学规划与Lingo程序设计奶制品的生产与销售一、实验问题一加工厂用牛奶生产A1,A2两种奶制品，1桶牛奶

22、可以在甲设备上用12小时生产成3公斤A1，或者在乙设备上用8小时加工成4公斤A2。根据市场需求，生产的A1,A2全部可以销售出去。每公斤A1可获利24元，每公斤A2获利16元。现在加工厂每天可以有50桶牛奶，每天正式工人总劳动时间为480小时，且甲类设备每天至多加工100公斤A1，乙类设备的加工能力没有上限。为增加工厂的获利，工厂开发了奶制品的深加工技术：用2小时和3元加工费，可将1公斤A1加工成0.8公斤高级奶制品B1，也可将1公斤A2加工成0.75公斤高级奶制品B2；每公斤B1能获利44元，每公斤B2能获利32元试为该厂制订一个生产销售计划，使每天的净利润最大，并讨论以下问题：1）若投资3

23、0元可以增加供应1桶牛奶，投资3元可以增加1小时劳动时间，应否作这些投资?若每天投资150元，可赚回多少?2）每公斤高级奶制品B1，B2的获利经常有10的波动，对制订的生产销售计划有无影响?若每公斤B1的获利下降10，计划应该变化吗?二、实验目的： 1）学会lingo软件的基本使用方法。2）加强学生对数学规划的认识和理解；培养学生建模的能力。三、LINGO预备知识1. LINGO入门当你在windows下运行LINGO系统,窗口外层是主框架窗口，包含了所有菜单命令和工具条，其它所有的窗口将被包含在主窗口之下。在主窗口内的标题为LINGO Model LINGO1的窗口是LINGO的默认模型窗口

24、，建立的模型都都要在该窗口内编码实现。例1.1 如何在LINGO中求解如下的LP问题：在模型窗口中输入如下代码：model:max = 72*x1+64*x2;milk x1 + x250;time 12*x1+8*x2480;cpct 3*x1100;end 注意：模型输入lingo软件时，乘号*不能省略，每句结束时必须加；。默认所有变量为非负，大于号和大于等于、小于号和小于等于等价。最后以end结束。2 基本运算符及函数2.1.算术运算符算术运算符是针对数值进行操作的。LINGO 提供了5 种二元运算符： 乘方 乘 除 加 减LINGO 唯一的一元算术运算符是取反函数“”。这些运算符的优先级由高到底为：高 （取反）低 运算符的运算次序为从左到右按优先级高低来执行。运算的次序可以用圆括号“（）”来改变。2.2 变量界定函数变量界定函数实现对变量取值范围的附加限制，共4 种：bin（x） 限制x 为0 或1bnd（L,x,U） 限制LxUf