北科大数学模型第二次拟合及插值实验.doc
实验一 拟合和插值教学目的 .了解最小二乘法的原理. .通过实例的学习,懂得如何用拟合和插值的方法解决实际的问题,并能注意它们的联系与区别,会用 Matlab来求解 教学内容 .拟合与插值的原理及简单分类 .相应问题的实例建模及用软件求解的实现 .练习与上机实验的内容 插值：求过已知有限个数据点的近似函数。拟合：已知有限个数据点，求近似函数，不要求过已知数据点，只要求在某种意义下它在这些点上的总偏差最小。插值和拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的许多工程技术中提出的计算问题对插值函数的光滑性有较高要求，如飞机的机翼外形，内燃机的进、排气门的凸轮曲线，都要求曲线具有较高的光滑程度，不仅要连续，而且要有连续的曲率，这就导致了样条插值的产生。所谓样条（Spline）本来是工程设计中使用的一种绘图工具，它是富有弹性的细木条或细金属条。绘图员利用它把一些已知点连接成一条光滑曲线（称为样条曲线），并使连接点处有连续的曲率。9拟合：Zj2.my97 =233.4286y98 =253.9286 x = -260.2912 264.3161 0.0006课堂练习与作业：1， 所有例题上机实现.见例题下方有答案输出2. x0=0:pi/10:pi/2;y0=sin(x0);f=polyfit(x0,y0,5);poly2str(f,'x')ans =0.0057725 x5 + 0.0059562 x4 - 0.17215 x3 + 0.0022051 x2 + 0.99969 x + 4.1234e-16function f=fun2(x,tdata);f=x(1)*exp(x(2)*tdata);xd=1:1:8;yd=15.3 20.5 27.4 36.6 49.1 65.6 87.87 117.6;x0=15.3 0;X=lsqcurvefit(fun2,x0,xd,yd)X =11.4251 0.29143. 用下列数据拟合函数中的参数。数据序号 y/kg x1/cm2 x2 x3 115.0223.735.491.211415.9423.525.181.98212.6222.344.321.351514.3321.864.861.59314.8628.845.041.921615.1128.955.181.37413.9827.674.721.491713.8124.534.881.39515.9120.835.351.561815.5827.655.021.66612.4722.274.271.501915.8527.295.551.70715.8027.575.251.852015.2829.075.261.82814.3228.014.621.512116.4032.475.181.75913.7624.794.421.462215.0229.655.081.701015.1828.965.301.66 2315.7322.114.901.811114.2025.774.871.642414.7522.434.651.821217.0723.175.801.902514.3520.045.081.531315.4028.575.221.66建立fun3.m文件function f=fun3(x,t)f=exp(-x(1).*t(:,1).*sin(x(2).*t(:,2)+t(:,3).*t(:,3);导入数据变量A、yy,yy为y值列矢量，A为x1、x2、x3值矩阵主程序x=lsqcurvefit(fun3,x0, A,yy);k1=x(1)k1=x(2)结果4.数学建模实例：人口预报问题 指数增长模型（马尔萨斯人口模型） 参数进行估计列表与画图分析； 参数进行估计列表与画图分析以美国1790-1990人口为例指数增长：建立maer.m文件：function f= maer(x,t)f=x(1)*exp(x(2)*t);主程序t=1:21;y=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.7 179.3 204.0 226.5 251.4;x0=0.05 0.05;x=lsqcurvefit(maer,x0, t,y);x0=x(1)r=x(2)y2000=x0*exp(r*22)y2015=x0*exp(r*37)求出x0，r画图对比预测2000年人数为309.2196百万人而2015年2839.6百万人,指数增长过快阻滞增长模型建立zuzhi.m文件function f=zuzhi(x,t)f=x(1)./(1+(x(1)./x(2)-1).*exp(-x(3).*t);主程序clear all t=1:21; y=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.7 179.3 204.0 226.5 251.4; x0=400 5.3 0.22;x=lsqcurvefit(zuzhi,x0,t,y);xm=x(1)x0=x(2)r=x(3)yy=xm./(1+(xm./x0-1).*exp(-r*t);plot(t,yy,'r-',t,y,'*')y2000=251.4+r*251.4*(1-251.4/xm)预测2000年人数为272.2099百万人2015年人数为392.87百万人与321.62百万人差距略大5. 2005年长江水质的评价与预测做第3问，预测未来十年可饮用水和第二类水的比例。x=xlsread('a.xls');for i=1:90 y(i)=x(i,3)/(x(i,1)+x(i,3)+x(i,5);endy=y'x1=174,179,183,189,207,234,220,256,270,285'x2=9205,9513,9171,13127,9513,9924,8892,10210,9980,9405'for i=1:10 y1(i)=y(9*i-8);endy1=y1'for i=1:10 y2(i)=y(9*i-7);endy2=y2'for i=1:10 y3(i)=y(9*i-6);endy3=y3'for i=1:10 y4(i)=y(9*i-5);endy4=y4'for i=1:10 y5(i)=y(9*i-4);endy5=y5'for i=1:10 y6(i)=y(9*i-3);endy6=y6'for i=1:10 y7(i)=y(9*i-2);endy7=y7'for i=1:10 y8(i)=y(9*i-1);endy8=y8'for i=1:10 y9(i)=y(9*i);endy9=y9'X=x1 x2 ones(10,1); b1,bint1,r1,rint1,stats1=regress(y1,X) b2,bint2,r2,rint2,stats2=regress(y2,X) b3,bint3,r3,rint3,stats3=regress(y3,X) b4,bint4,r4,rint4,stats4=regress(y4,X) b5,bint5,r5,rint5,stats5=regress(y5,X) b6,bint6,r6,rint6,stats6=regress(y6,X) b7,bint7,r7,rint7,stats7=regress(y7,X) b8,bint8,r8,rint8,stats8=regress(y8,X) b9,bint9,r9,rint9,stats9=regress(y9,X) x3=294.4 304.7 315.5 326.8 338.5 350.8 363.6 376.9 390.6 404.9' x4=9894.1 9894.1 9894.1 9894.1 9894.1 9894.1 9894.1 9894.1 9894.1 9894.1' for i=1:10 y10(i)=b1(1)*x3(i)+b1(2)*x4(i)+b1(3); y11(i)=b2(1)*x3(i)+b2(2)*x4(i)+b2(3); y12(i)=b3(1)*x3(i)+b3(2)*x4(i)+b3(3); y13(i)=b4(1)*x3(i)+b4(2)*x4(i)+b4(3); y14(i)=b5(1)*x3(i)+b5(2)*x4(i)+b5(3); y15(i)=b6(1)*x3(i)+b6(2)*x4(i)+b6(3); y16(i)=b7(1)*x3(i)+b7(2)*x4(i)+b7(3); y17(i)=b8(1)*x3(i)+b8(2)*x4(i)+b8(3); y18(i)=b9(1)*x3(i)+b9(2)*x4(i)+b9(3); end结果：y10 = 0.4724 0.4801 0.4881 0.4965 0.5052 0.5143 0.5238 0.5337 0.5439 0.5545y11 = 0.4012 0.4044 0.4078 0.4113 0.4150 0.4188 0.4228 0.4270 0.4313 0.4358y12 = 0.4679 0.4730 0.4785 0.4841 0.4900 0.4962 0.5026 0.5092 0.5161 0.5233y13 = 0.4429 0.4530 0.4636 0.4748 0.4863 0.4984 0.5110 0.5241 0.5376 0.5516y14 = 0.4569 0.4672 0.4779 0.4891 0.5007 0.5129 0.5256 0.5388 0.5524 0.5666y15 = 0.4429 0.4534 0.4644 0.4759 0.4878 0.5004 0.5134 0.5270 0.5409 0.5555y16 = 0.5339 0.5500 0.5669 0.5846 0.6029 0.6222 0.6422 0.6631 0.6845 0.7069y17 = 0.3983 0.3998 0.4014 0.4031 0.4049 0.4067 0.4086 0.4106 0.4127 0.4148y18 =0.5687 0.5889 0.6100 0.6321 0.6550 0.6790 0.7041 0.7301 0.7569 0.78496、2011赛题对高度与浓度进行插值；x1=xlsread('ha.xls','¸½¼þ1');x=x1(:,2);y=x1(:,3);z=x1(:,4);x2,y2=meshgrid(0:300:30000,0:200:20000);z1=griddata(x,y,z,x2,y2,'v4');subplot(1,2,1)contourf(x2,y2,z1,15);subplot(1,2,2)mesh(x2,y2,z1)x1=xlsread('ha.xls','¸½¼þ1');x2=xlsread('ha.xls','¸½¼þ2');x=x1(:,2);y=x1(:,3);z=x2(:,2);x2,y2=meshgrid(0:300:30000,0:200:20000);z1=griddata(x,y,z,x2,y2,'v4');subplot(1,2,1)contourf(x2,y2,z1,15);subplot(1,2,2)mesh(x2,y2,z1)x1=xlsread('ha.xls','¸½¼þ1');x2=xlsread('ha.xls','¸½¼þ2');x=x1(:,2);y=x1(:,3);z=x2(:,3);x2,y2=meshgrid(0:300:30000,0:200:20000);z1=griddata(x,y,z,x2,y2,'v4');subplot(1,2,1)contourf(x2,y2,z1,15);subplot(1,2,2)mesh(x2,y2,z1)x1=xlsread('ha.xls','¸½¼þ1');x2=xlsread('ha.xls','¸½¼þ2');x=x1(:,2);y=x1(:,3);z=x2(:,4);x2,y2=meshgrid(0:300:30000,0:200:20000);z1=griddata(x,y,z,x2,y2,'v4');subplot(1,2,1)contourf(x2,y2,z1,15);subplot(1,2,2)mesh(x2,y2,z1)7、2011赛题对各指标进行主成分分析。各成分分析如下：Cux1=xlsread('ha.xls','sheet1');x2=xlsread('ha.xls','sheet2');x0=x1(:,2:4);y0=x2(:,5);r=corrcoef(x0);xb=zscore(x0);yb=zscore(y0);c,s,t=princomp(xb)contr=cumsum(t)/sum(t)c = 0.5873 -0.2855 0.7573 0.5639 0.8155 -0.12990.5805 -0.5034 -0.6400contr = 0.6501 0.8363 1.0000Hgx1=xlsread('ha.xls','sheet1');x2=xlsread('ha.xls','sheet2');x0=x1(:,2:4);y0=x2(:,6);r=corrcoef(x0);xb=zscore(x0);yb=zscore(y0);c,s,t=princomp(xb)contr=cumsum(t)/sum(t)c = 0.5873 -0.2855 0.7573 0.5639 0.8155 -0.12990.5805 -0.5034 -0.6400contr = 0.6501 0.8363 1.0000Nix1=xlsread('ha.xls','sheet1');x2=xlsread('ha.xls','sheet2');x0=x1(:,2:4);y0=x2(:,7);r=corrcoef(x0);xb=zscore(x0);yb=zscore(y0);c,s,t=princomp(xb)contr=cumsum(t)/sum(t)c = 0.5873 -0.2855 0.7573 0.5639 0.8155 -0.12990.5805 -0.5034 -0.6400contr = 0.6501 0.8363 1.0000Pbx1=xlsread('ha.xls','sheet1');x2=xlsread('ha.xls','sheet2');x0=x1(:,2:4);y0=x2(:,8);r=corrcoef(x0);xb=zscore(x0);yb=zscore(y0);c,s,t=princomp(xb)contr=cumsum(t)/sum(t)c = 0.5873 -0.2855 0.7573 0.5639 0.8155 -0.12990.5805 -0.5034 -0.6400contr = 0.6501 0.8363 1.0000Znx1=xlsread('ha.xls','sheet1');x2=xlsread('ha.xls','sheet2');x0=x1(:,2:4);y0=x2(:,9);r=corrcoef(x0);xb=zscore(x0);yb=zscore(y0);c,s,t=princomp(xb)contr=cumsum(t)/sum(t)c = 0.5873 -0.2855 0.7573 0.5639 0.8155 -0.12990.5805 -0.5034 -0.6400contr = 0.6501 0.8363 1.0000Asx1=xlsread('ha.xls','sheet1');x2=xlsread('ha.xls','sheet2');x0=x1(:,2:4);y0=x2(:,2);r=corrcoef(x0);xb=zscore(x0);yb=zscore(y0);c,s,t=princomp(xb)contr=cumsum(t)/sum(t)c = 0.5873 -0.2855 0.7573 0.5639 0.8155 -0.12990.5805 -0.5034 -0.6400contr = 0.6501 0.8363 1.0000Cdx1=xlsread('ha.xls','sheet1');x2=xlsread('ha.xls','sheet2');x0=x1(:,2:4);y0=x2(:,3);r=corrcoef(x0);xb=zscore(x0);yb=zscore(y0);c,s,t=princomp(xb)contr=cumsum(t)/sum(t)c = 0.5873 -0.2855 0.7573 0.5639 0.8155 -0.1299 0.5805 -0.5034 -0.6400contr = 0.6501 0.8363 1.0000Crx1=xlsread('ha.xls','sheet1');x2=xlsread('ha.xls','sheet2');x0=x1(:,2:4);y0=x2(:,4);r=corrcoef(x0);xb=zscore(x0);yb=zscore(y0);c,s,t=princomp(xb)contr=cumsum(t)/sum(t)c = 0.5873 -0.2855 0.7573 0.5639 0.8155 -0.12990.5805 -0.5034 -0.6400contr = 0.6501 0.8363 1.0000课外拓展：7、对2011赛题对高度与浓度进行克立金插值或者shepard插值，二选一；8、2005年长江水质的评价与预测做第4问，如何控制？
