2013数学建模论文(高级数值分析)

1、2013上数学建模论文题目：城市空气质量评估及预测学院：信息学院班级：10网络工程2学号：201030520118姓名：梁智扬2013年6月21日一、问题的提出我国是一个人口大国，城市众多，但由于工业的发展，我们的很多城市都受到了不同程度的污染尤其是空气污染，直接对我们造成伤害。近年来，虽然我国大气污染防治工作取得了很大的成效，但由于各种原因，我国大气环境面临的形势仍然非常严峻。要实现大气环境与社会经济协调发展，开展城市环境空气质量影响因素及其影响机制的研究是十分必要的，也是改善城市环境质量迫切需要解决的实际问题。1）从中国的实际情况的特点出发，利用附件中的相关数据（以及相关文献和有关补充数据

2、），建立数学模型给出十个城市空气污染严重程度的科学排名；2）建立模型对广州市11月的空气质量状况进行预测；3）收集必要的数据，建立模型分析影响城市空气污染程度的主要因素；4）评价所建立的数学模型中的优点和存在的不足。二、模型的假设1）广州周边风速适中，扩散一般，忽略周围城市污染物的扩散对广州市的影响;2）在预测月份中，气候未出现大的异常；3）与相似年份的数据基本吻合，未出现大规模变化，污染潜势和气象参数基本一样；4）相似年份的数据来源可靠，具有使用价值；5）广州市的车辆未出现大幅增加或减少，大范围工厂兴建和拆毁；6）空气质量的原始数据并不一定完全相同，但为了简化问题，计算所得出来的数据只要是同

3、一级别，如35和47一样，无区别；7）预测广州11月份的空气质量时，忽略特殊年份的特殊数据，假设其对结果无影响；8）广州在相似年份并没有出现中度、重度污染，故忽略这两种情况，得出数据无效。四、问题的分析与模型的建立由于空气质量的高低一般和年份、临近月份存在关系，而且会因为季节的不同而不同，比如2005、2006等年份的空气质量大体相同但也存在联系，再比如有些地区的冬季和春季的空气质量也存在差异。研究某地区未来空气质量可以借助以前已有的数据，利用最小二乘法对数据进行拟合，最后用Matlab软件求解拟合曲线方程。其次在分析数据时为了避免数据太少对预测结果造成偶然误差，在预测某一月的空气质量时可以将

4、其相似月份一起分析。在分析往年数据时发现规律，有如下数据统计：编号污染指数平均数IIIIII1III2178.11067131279.01069111379.21269100469.1147250581.3156592其中：第一列编号1,2,3,4,5分别代表2005,2006,2007,2008,2009；第二列数值表示九月十月十一月污染指数的平均值；第三列数值表示空气质量优的天数；第三列数值表示空气质量良的天数；第四列表示空气轻微污染的天数；第五列数值表示空气轻度污染的天数。污染指数平均数从这些散点图可以看出除了第四个点的分布比较特殊外，其它四个点的分布还是比较有规律的，可以看出这四个点基

5、本呈线性分布。一等级出现次数从一等级出现次数我们可以看出这些散点图除第一个点之外，其余四点的分布还是很有规律的，这四个点呈抛物线态分布。二等级出现次数从二等级出现次数的散点图来看，除了第四个点之外，其它四点的分布还是很有规律的，这四点呈开口向下的抛物线分布。1412108642三1等级出现次数*A*系列1*1111100123456从三1等级出现次数的散点图可以看出除了第四个点外，其余四点的分布有规律，基本呈开口向上的抛物线。再来单独考虑空气质量等级为三2的情况，由于广州以往的数据中这三个月之中出现三2等级的次数很少，所以可以先不考虑。4232曲线拟合的最小二乘法方法介绍多项式数据拟合【3】由

6、于预测广州十一月份的空气质量是借助于以往的数据，分析反映空气质量的各种指标的分布规律，又因为往年数据有限，所以给预测带来了困难，这不仅是发现其规律有困难，还会使预测结果产生很大的误差，一旦方法选择不好，误差就可能严重偏离实际。为此，类如插值法求解方程，待定系数法求解方程等方法就不适和本题，插值法对外插的值收敛性要求比较高，而且外插的值误差也较大；待定系数法必须先发现公式然后求解，显然不符合本模型。所以对于预测广州十一月份的空气质量还是选用最小二乘法，有两个优点，其一，最小二乘法求出的结果误差小，可以最大程度上拟合这些点，其二，最小二乘法的求解可以利用Matlab等软件求解。下面为最小二乘法的原

7、理：首先可以在Matlab或其它画图软件中描出散点图，观察散点图的分布特征选择合适的多项式，如设用一个m次多项式mPm(=)=a0-ax-a2x2amxm=丐农来拟合一组数据（二，T：）（i=l,2,，n）。但拟合出来的直线存在一定的偏差，由散点图可以看出影响空气质量的几种参数随时间的分布并不规律，而且由于某些特殊的点还要舍去，所以分析偏差很重要。因而可以如下进行偏差的分析：在节点二处，汎：二）的偏差设为最少二乘法的基本思想是对所有数据点，拟合函数式的偏差二的平方和nnnm=吩X（亏）一h=丐珀2取最小值。由于（:八?：）为已知，故将e视作H二）的系数扎（j=0,l,，m）的函数。不同的拟合多

8、项式，有不同的一组系数、：（j=0,l,，m）,因而有不同的e值即于是上述拟合问题便归结为求求多元函数的极值问题，欲使0=0日0，归，.-日口）取最小值，贝丘J匚，必须满足条件30=0（k=l,2,，m）经过上面的方法处理之后可以得出使偏差最小的表达式，令=工严tk=Syixfi则式可以表达成为:#includevstdlib.h#includevmalloc.h#includeSmooth(double*x,double*y,double*a,intn,intm,double*dtl,double*dt2,double*dt3);voidmain()inti,n,m;double*x,*y,

9、*a,dt1,dt2,dt3,b;n=12;/12个样点m=4;3次多项式拟合b=0;/x的初值为0/*分别为x,y,a分配存贮空间*/x=(double*)calloc(n,sizeof(double);if(x=NULL)printf(”内存分配失败n);exit(0);y=(double*)calloc(n,sizeof(double);if(y=NULL)printf(内存分配失败n);exit(0);a=(double*)calloc(n,sizeof(double);if(a=NULL)printf(内存分配失败n);exit(0);for(i=1;iv=n;i+)xi-1=b+(

10、i-1)*5;/*每隔5取一个点，这样连续取12个点*/y0=0;y1=1.27;y2=2.16；y3=2.86;y4=3.44;y5=3.87;y6=4.15;y7=4.37;y8=4.51;y=4.58;y10=4.02;y11=4.64;/*xi-1点对应的y值是拟合已知值*/Smooth(x,y,a,n,m,&dt1,&dt2,&dt3);/*调用拟合函数*/for(i=1;iv=m;i+)printf(a%d=%.10fn,(i-1),ai-1);printf(拟合多项式与数据点偏差的平方和为：n);printf(%.10en,dt1);printf(拟合多项式与数据点偏差的绝对值之

11、和为：n);printf(%.10en,dt2);printf(拟合多项式与数据点偏差的绝对值最大值为：n);printf(%.10en,dt3);free(x);/*释放存储空间*/free(y);/*释放存储空间*/free(a);/*释放存储空间*/Smooth(double*x,double*y,double*a,intn,intm,double*dt1,double*dt2,double*dt3)/(x,y,a,n,m,dt1,dt2,dt3)/double*x;/*实型一维数组，输入参数，存放节点的xi值*/double*y;/*实型一维数组，输入参数，存放节点的yi值*/doub

12、le*a;/*双精度实型一维数组，长度为m。返回m1次拟合多项式的m个系数*/intn;/*整型变量，输入参数，给定数据点的个数*/intm;/*整型变量，输入参数，拟合多项式的项数*/double*dt1;/*实型变量，输出参数，拟合多项式与数据点偏差的平方和*/double*dt2;/*实型变量，输出参数，拟合多项式与数据点偏差的绝对值之和*/double*dt3;/*实型变量，输出参数，拟合多项式与数据点偏差的绝对值最大值*/inti,j,k;double*s,*t,*b,z,d1,p,c,d2,g,q,dt;/*分别为s,t,b分配存贮空间*/s=(double*)calloc(n,s

13、izeof(double);if(s=NULL)printf(”内存分配失败n);exit(0);t=(double*)calloc(n,sizeof(double);if(t=NULL)printf(内存分配失败n);exit(0);b=(double*)calloc(n,sizeof(double);if(b=NULL)printf(内存分配失败n);exit(0);z=0;for(i=l;iv=n;i+)z=z+xi-l/n;/*z为各个x的平均值*/b0=1;d1=n;p=0;c=0;for(i=1;i1)t1=1;t0=-p;d2=0;c=0;g=0;for(i=1;i=4)for(k=j-2;k=2;k-)sk-1=-p*tk-1+tk-2-q*bk-1;sO=-p*tO-q*bO;d2=0;c=0;g=0;for(i=1;i=1;k-)q=q*(xi-1-z)+sk-1;d2=d2+q*q;c=yi-1*q+c;g=(xi-1-z)*q*q+g;c=c/d2;p=g/d2;q=d2/d1;d1=d2;aj-1=c*sj-1;tj-1=sj-1;for(k=j-1;k=