数学建模比赛前准备的Matlab和lingo代码

1、Matlab和lingo代码Matlab1基础知识1Polyval3Polyfit3interrep13回归分析4牛顿迭代法求解非线性方程组6建模课上的代码12lingo求解部分18目标规划18第10章 数据的统计描述和分析23!7个工人，7个工作的分配问题;24案例分析25差分方程28!三阶段面试模型;30装配线平衡模型32露天矿生产的车辆安排（CMCM2003B）34Matlab基础知识相关系数矩阵的方式，通过Matlab软件进行相关性分析，得到主成分种类与重要指标的线性组合： （10）prod 连乘积for k=1:100p(k)=1-prod(365-k+1:365)/365k;end

2、plot(p)=不等于 数论函数gcd(a,b)两个整数的最大公约数lcm(a,b)两个整数的最小公倍数求整函数与截尾函数ceil(x)表示大于或等于实数x的最小整数floor(x)表示小于或等于实数x的最大整数round(x)最接近x的整数syms 表达式中包含的变量 factor(表达式) 因式分解syms 表达式中包含的变量 expand(表达式) 代数式展开 syms 表达式中包含的变量 collect(表达式，指定的变量) 合并同类项syms 表达式中包含的变量 simplify(表达式) 进行数学式化简syms 表达式和代换式中包含的所有变量 subs(表达式，要替换的变量或式子，

3、代换式) syms x Taylor(f(x), x, n, a)fplot(f(x)，xmin,xmax,ymin,ymax) syms x int(f(x), x,a,b)Polyval 计算对多项式p(x)=1+2*x+3*x2，计算在x=5,7,9的值。 p = 3 2 1; x=5,7,9; polyval(p,5 7 9)%结果为ans =86 162 262Polyfit 拟合曲线x=1,2,4,7,9,12,13,15,17; F=1.5,3.9,6.6,11.7,15.6,18.8,19.6,20.6,21.1;plot(x,F,.)%从图像上我们发现:前5个数据应与直线拟合

4、,后5个数据应与二次曲线拟合。于是键入 : a=polyfit(x(1:5),F(1:5),1); a=polyfit(x(5:9),F(5:9),2)生日概率模型for n=1:100p(n)=1-prod(365-n+1:365)/365n;endplot(p)c5=polyfit(n,p,5)c5 =-0.0000 0.0000 -0.0001 0.0023 -0.0046 -0.0020该多项式即为：在Matlab环境下继续键入下列指令： p5=polyval(c5,n); /用多项式近似计算100个概率值 plot(n,p,n,p5,.) /画出拟合多项式的图象与概率曲线作比较int

5、errep1x0=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.6plot(x0,y0) %完成第一步工作x=0:0.1:15;y=interp1(x0,y0,x); %用分段线性插值完成第二步工作plot(x,y)y=spline(x0,y0,x); plot(x,y) %用三次样条插值完成第二步工作指数模型t=1790:10:1980;x(t)=3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76.0 92.0 106.5 123.2 131.7 150.

6、7 179.3 204.0 226.5 ;y=log(x(t);a=polyfit(t,y,1)r=a(1),x0=exp(a(2)x1=x0.*exp(r.*t);plot(t,x(t),r,t,x1,b)%阻滞增长模型（或 Logistic 模型）%建立函数文件curvefit_fun2.mfunction f=curvefit_fun2 (a,t)f=a(1)./(1+(a(1)/3.9-1)*exp(-a(2)*(t-1790);在命令文件main.m中调用函数文件curvefit_fun2.m % 定义向量（数组）x=1790:10:1990;y=3.9 5.3 7.2 9.6 12

7、.9 17.1 23.2 31.4 38.6 50.2 62.9 76 . 92 106.5 123.2 131.7 150.7 179.3 204 226.5 251.4;plot(x,y,*,x,y); % 画点，并且画一直线把各点连起来hold on; a0=0.001,1; % 初值% 最重要的函数，第1个参数是函数名（一个同名的m文件定义），第2个参数是初值，第3、4个参数是已知数据点a=lsqcurvefit(curvefit_fun2,a0,x,y); disp(a= num2str(a); % 显示结果% 画图检验结果xi=1790:5:2020;yi=curvefit_fun

8、2(a,xi);plot(xi,yi,r);% 预测2010年的数据x1=2010;y1=curvefit_fun2(a,x1)hold off回归分析设回归模型为 y=0+1x,在MATLAB命令窗口中键入下列命令进行回归分析(px_reg11.m)x=0.1:0.01:0.18;x=x,0.2,0.21,0.23;y=42,41.5,45,45.5,45,47.5,49,55,50,55,55.5,60.5;X=ones(12,1),x; %一元回归b,bint,r,rint,stats=regress(y,X,0.05);b,bint,stats,rcoplot(r,rint)得结果和图

9、b = 27.0269 140.6194bint = 22.3226 31.7313 111.7842 169.4546stats = 0.9219 118.0670 0.0000 3.1095结果含义为0=27.0269 1=140.61940的置信区间是 22.3226,31.73131的置信区间是 111.7842,169.4546R2=0.9219 F=118.0670, p10-4.R是衡量y与x的相关程度的指标，称为相关系数。R越大，x与y关系越密切。通常R大于0.9才认为相关关系成立。F是一统计指标p是与F对应的概率，当 p0.05时，回归模型成立。此例中 p=0 10-40.0

10、5，所以，所得回归模型成立。牛顿迭代法求解非线性方程组1实验原理对于非线性方程在处按照多元函数的泰勒展开，并取线性项得到其中2 数据来源计算非线性方程组给定初值，要求求解精度达到0.000013 实验步骤步骤一：首先建立函数,方程编程如下，将F.m保存到工作路径中：function f=F(x)f(1)=x(1)2-10*x(1)+x(2)2+8;f(2)=x(1)*x(2)2+x(1)-10*x(2)+8;f=f(1),f(2) ;步骤二：建立函数,用于求方程的jacobi矩阵,将DF.m保存到工作路径中：function df=DF(x)df=2*x(1)-10,2*x(2);x(2)2+

11、1,2*x(1)*x(2)-10; %jacobi矩阵是一阶偏导数以一定方式排列成的矩阵。步骤三：编程牛顿迭代法解非线性方程组，将newton.m保存在工作路径中：clear,clc;x=0,0;%求第二组解x=4,4f=F(x);df=DF(x);fprintf(%d %.7f %.7fn,0,x(1),x(2);N=4;for i=1:Ny=dff;x=x-y;f=F(x);df=DF(x);fprintf(%d %.7f %.7fn,i,x(1),x(2);if norm(y)0.break;elseendendh=ezplot(x2-10*x+y2+8,-6,6,-6,6);set(h

12、,color,r);hold onezplot(x*y2+x-10*y+8,-6,6,-6,6);grid on运行结果如下：第一组解：0 0. 0.1 0. 0.2 0. 0.3 0. 0.4 1. 1.第二组解：0 4. 4.1 2. 3.2 2. 3.3 2. 3.4 2. 3.function x2=newtonmethod(f,x0,tol) %牛顿迭代法% f is the symbol varsyms x;k=0; maxlter=1000; x1=x0; df=diff(f);while kmaxlter k=k+1 x2=x1-double(subs(f,x,x1)/subs

13、(df,x,x1) if abs(x2-x1)tol disp(found!); break; end if abs(subs(f,x,x2)tol disp(found!); break; end if k=maxlter-1 error(reach maxlter!); break; end x1=x2;endclear,clcsyms x;p1=2; p2=3; h1=5; h2=6; s=20; x0=5;dc=-3*p1*h1*x/sqrt(h12+x2)5)+3*p2*h2*(s-x)/sqrt(h22+(s-x)2)5;x=newtonmethod(dc,x0,1e-6)func

14、tion xk=fast_down(f,x0,tol) %最速下降法% f is the symbol functionsyms x1 x2 t;k=0; maxlter=1000; xk=x0; p=-1*diff(f,x1); diff(f,x2);while kmaxlter k=k+1; pk=double(subs(p,x1,x2,xk(1,k),xk(2,k); g=subs(f,x1,x2,xk(1,k)+t*pk(1),xk(2,k)+t*pk(2); kk=0; t1=0; while kkmaxlter kk=kk+1; dg=diff(g,t); ddg=diff(dg,

15、t); t2=t1-subs(dg,t,t1)/subs(ddg,t,t1); if abs(t2-t1)tol break; end if kk=maxlter-1 error(found t reach maxlter!); end t1=t2; end xk(:,k+1)=xk(:,k)+t2*pk; if norm(xk(:,k+1)-xk(:,k)0.98plot(j:j+pipei_width,i,r);plot(j:j+pipei_width,i+pipei_height,r);plot(j,i:i+pipei_height,r);plot(j+pipei_width,i:i+p

16、ipei_height,r);endendendtocfigure;imshow(pipeitu2gray);二元线性回归方程求红葡萄质量与红葡萄酒理化指标对红葡萄酒质量的影响x1=1.731 1.048 2.094 0.727 0.883 0.586 0.62 1.88 2.13 1.096 0.758 0.493 0.679 1.468 0.673 0.646 0.958 0.578 0.838 0.921 0.659 1.227 0.969 0.692 1.538 0.666；x2=1. 1. 1. 0. 1. 0. 0. 1. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0.

17、 1. 1. 1. 1. 1. 0. 0. 0. ；y=8.585 9.7769.488 8.794 8.68 8.607 8.5669.3639.338 8.404 8.642 7.741 9.257 8.696 8.054 9.282 9.0978.39 9.662 9.0789.3439.5259.829.0768.725 9.085； X=ones(28,1),x1,x2； b,bint,r,rint,s=regress(y,X)syms a b l h ; syms x; f = 2 * a *sqrt(1-(x/b-1)2) * l; v = int(f,x,0,h) %对表达式求

18、积分jinyou1 = xlsread(问题A附件1：实验采集数据表.xls,无变位进油,C2:D79); chuyou1 = xlsread(问题A附件1：实验采集数据表.xls,无变位出油,C2:D75); g1=jinyou1(:,2)*10(-2); %高度MM g1 乘以0.01使得x轴为个位数g2=chuyou1(:,2)*10(-2); %高度MM g2syms a b h l a=0.89*10;b=0.6*10;l=2.45*10;V=(a*l/b).*(h-b).*(2.*h.*b-h.2).0.5)+b.2.*asin(h-b)./b)+0.5.*b.2*pi);V1=s

19、ubs(V,h,g1) ; %理论上进油后储油量V1V2=subs(V,h,g2) ; %理论上进油后储油量V2V3=jinyou1(:,1)+262; %实际总油V3V4=3706.91+262-chuyou1(:,1); %实际剩下油量Lc1=V1-V3; %理论减掉实际的误差c2=V2-V4;x1=abs(c1)./V1; %相对误差x2=abs(c2)./V2; figure(1)plot(g1,V1,.r,g1,V3,.c);xlabel(进油后储油量); %h1ylabel(罐内油位高度); %V1title(罐内油位高度随进油后储油量的变化情况)legend(储油理论量,储油实际

20、量)grid on;附录三（三次样条插值得到的罐体实际储油量L）data= xlsread(问题A附件2：实际采集数据表.xls,实际储油罐的采集数据,E2:F604); h=data(:,1);v=data(:,2);hi=300:.001:2700;vi=interp1(h,v,hi,spline);plot(h,v,hi,vi,r);hi=400:100:2700;vi=interp1(h,v,hi,spline)DVD在线租赁根据上述模型，编写Lingo程序如下：model:sets:dvd/1.100/:total;!total为网站现有的数量;huiyuan/1.1000/; !自

21、定义0-1向量;pianai(huiyuan,dvd):data0;bianliang(huiyuan,dvd):data1;Endsets!目标函数（总体满意度最大）;max=sum(pianai(i,j):data1(i,j)*data0(i,j);! x (i,j) ,t (i)为0-1变量;for(bianliang:bin(data1); ! 约束条件每人可租数Yi=3或0;for(huiyuan(j):sum(bianliang(j,i):data1(j,i)=3); ! 所租总数小于网站现有的数量;for(dvd(i):sum(bianliang(i,j):data1(i,j)t

22、ol*abs(x(k) x(k+1)=x(k)-feval(fv,x(k)/feval(df,x(k); b=x(k+1)-x(k); k=k+1; if(kn) error(Error: Reached maximum iteration times); endendy=x(k-1);if nargout1 z=k-1;endclear allx0=1;y0=1;n=10;tol=1e-6;x(1)=x0;y(1)=y0;i=1;u=1 1;k(1)=1;while(norm(u)tol*norm(x(i),y(i) A=2*x(i),y(i);x(i),-y(i); b=4-x(i)2-y

23、(i)2,1-x(i)2+y(i)2; u=Ab; x(i+1)=x(i)+u(1); y(i+1)=y(i)+u(2); i=i+1;k(i)=i; if(in)error(n is full); endendk,x,yA=0.15 0.1 0.2;0.3 0.05 0.3;0.2 0.3 0;d=50 150 100;B=eye(3)-A;x=Bdc=inv(B)v=20:20:140/3.6; v2=v.2;x=v;v2; d=6.5,17.8,33.6,57.1,83.4,118.0,153.5; % 输出最小二乘解 左除a=xd;k1=a(1), k2=a(2)dd=x*a %辛普森

24、公式z1=quad(1./(1-sin(x),0,pi/4) %与精确值比较dz1=z1-sqrt(2)%Gauss-Lobatto公式z2=quadl(1./(1-sin(x),0,pi/4)%与精确值比较dz2=z2-sqrt(2)%梯形公式 x=0:pi/400:pi/4; y=1./(1-sin(x);z3=trapz(y)*pi/400z4=trapz(x,y)%与精确值比较dz3=z3-sqrt(2)dz4=z4-sqrt(2) % 按照表1输入原始数据x=0:0.2:5, 4.8:-0.2:0;y=5 4.71 4.31 3.68 3.05 2.5 2.05 1.69 1.4 1

25、.18 1 0.86 0.74 0.64 0.57 0.5 . 0.44 0.4 0.36 0.32 0.29 0.26 0.24 0.2 0.15 0 -1.4 -1.96 -2.37 -2.71 . -3 -3.25 -3.47 -3.67 -3.84 -4 -4.14 -4.27 -4.39 -4.49 -4.58 -4.66 . -4.74 -4.8 -4.85 -4.9 -4.94 -4.96 -4.98 -4.99 -5; % 逆时针方向转90度，节点（x, y）变为（u, v）v0=x; u0=-y; % 按0.05的间隔在u方向产生插值点u=-5:0.05:5; % 在v方向计

26、算分段线性插值v1=interp1(u0,v0,u); % 在v方向计算三次样条插值v2=spline(u0,v0,u); % 在（x, y）坐标系输出结果v1 v2 -u subplot(1,3,1),plot(x,y),axis(0 5 -5 5) gtext(原轮廓线,FontSize,12)subplot(1,3,2),plot(v1,-u),axis(0 5 -5 5) gtext(分段线性插值,FontSize,12)subplot(1,3,3),plot(v2,-u),axis(0 5 -5 5) gtext(三次样条插值,FontSize,12)lingo求解部分6发点8收点运

27、输问题model:!6发点8收点运输问题;sets: warehouses/wh1.wh6/: capacity; vendors/v1.v8/: demand; links(warehouses,vendors): cost, volume;endsets!目标函数; min=sum(links: cost*volume);!需求约束; for(vendors(J): sum(warehouses(I): volume(I,J)=demand(J);!产量约束; for(warehouses(I): sum(vendors(J): volume(I,J)= required(J);end计算

28、的部分结果为Global optimal solution found at iteration: 0 Objective value: 22.00000 Variable Value Reduced Cost REQUIRED( MON) 20.00000 0. REQUIRED( TUE) 16.00000 0. REQUIRED( WED) 13.00000 0. REQUIRED( THU) 16.00000 0. REQUIRED( FRI) 19.00000 0. REQUIRED( SAT) 14.00000 0. REQUIRED( SUN) 12.00000 0. START

29、( MON) 8. 0. START( TUE) 2. 0. START( WED) 0. 0. START( THU) 6. 0. START( FRI) 3. 0. START( SAT) 3. 0. START( SUN) 0. 0.从而解决方案是：每周最少需要22个职员，周一安排8人，周二安排2人，周三无需安排人，周四安排6人，周五和周六都安排3人，周日无需安排人。例4.18 求解最优化问题其LINGO代码如下：model: min=fx+fy; fx=if(x #gt# 0, 100,0)+2*x; fy=if(y #gt# 0,60,0)+3*y; x+y=30;end1 4 灵敏

30、性分析（Range，Ctrl+R）用该命令产生当前模型的灵敏性分析报告：研究当目标函数的费用系数和约束右端项在什么范围（此时假定其它系数不变）时，最优基保持不变。灵敏性分析是在求解模型时作出的，因此在求解模型时灵敏性分析是激活状态，但是默认是不激活的。为了激活灵敏性分析，运行LINGO|Options，选择General Solver Tab， 在Dual Computations列表框中，选择Prices and Ranges选项。灵敏性分析耗费相当多的求解时间，因此当速度很关键时，就没有必要激活它。下面我们看一个简单的具体例子。例5.1某家具公司制造书桌、餐桌和椅子，所用的资源有三种：木料

31、、木工和漆工。生产数据如下表所示：每个书桌每个餐桌每个椅子现有资源总数木料8单位6单位1单位48单位漆工4单位2单位1.5单位20单位木工2单位1.5单位0.5单位8单位成品单价60单位30单位20单位若要求桌子的生产量不超过5件，如何安排三种产品的生产可使利润最大？用DESKS、TABLES和CHAIRS分别表示三种产品的生产量，建立LP模型。max=60*desks+30*tables+20*chairs;8*desks+6*tables+chairs=48;4*desks+2*tables+1.5*chairs=20;2*desks+1.5*tables+.5*chairs=8;tabl

32、es=5;求解这个模型，并激活灵敏性分析。这时，查看报告窗口（Reports Window）,可以看到如下结果。Global optimal solution found at iteration: 3 Objective value: 280.0000 Variable Value Reduced Cost DESKS 2. 0. TABLES 0. 5. CHAIRS 8. 0. Row Slack or Surplus Dual Price 1 280.0000 1. 2 24.00000 0. 3 0. 10.00000 4 0. 10.00000 5 5. 0.“Global optimal solution found at iteration: 3”表示3次迭代