长江水污染状况与预测
一、 问题重述
长江作为中国第二大河，对我国的农业、经济、人民生活有着巨大的影响。
而南京 作为江苏的省会，也是长江流域的下游地区，城市的工业，农业以及生活用水大多来自长江。
然而，近年来，长江的水质不断下降，着实堪忧。
现有2003-2013年七月份长江水质pH 、DO 、CODMn 、NH3-N 等数据，通过数学建模的方法，对长江水污染状况进行预测，并提出合理的解决方案。
数据见附录。
二、 问题分析
题目中给出给出长江水质pH 、DO 、CODMn 、NH3-N 等数据，但是，这些数据比较少，并且，这四个参考量中，并不能看出哪一个或那几个参考量对长江水质有着巨大的影响，因此，需要采用灰色系统中的GM （1,1）来进行数据的预测，并进行趋势图的分析，找出主要的影响因素，在相对应的提出解决办法。
三、模型假设
1.所有数据都是真实有效的。
2.长江的水质变化满足一定函数，并且是光滑连续的。
3.长江的自净能力与纳污能力是一定的，而且沿岸的排污量满足一定的光滑连续的函数。
四、模型建立与求解
(1)数据的检验与处理
原始数据列为(0)(0)(0)(0)
(0)(1)(2)(3)(n)x =(x ,x ,x ,...,x )
,数列的级比为
(0)(k-1)(0)(k)
x λ(k)=
x ，若落在可容覆盖区间2
21
1
(,)n n X e
e
-
++=,则建立GM(1，1)模型；
若不落在可容覆盖区间内，则进行适当变换处理，如平移变换，取c ，使(0)(0)
()()k k y x c
=+的级比落在可容区间内。
对于pH:
(0)x =(7.41,7.45,7.47). DO:
(0)x =(6.37，5.88，6.95). CODMn:
(0)x =(2.5,1.8,3.8). NH3-N: (0)x =(0.1,0.22,0.26).
经程序计算后：
（2）建立GM （1,1）模型
1.原始数列经1次累加生成数列为(1)(1)(1)(1)
(1)(2)()(,,...,)n x x x x =，即(1)(0)()
()1
k
k i i x
x ==∑。
2.定义(1)
x
的灰导数为
(0)(1)(1)
()()(1)()k k k d k x x x -==-
令(1)(1)(1)(1)(2)(3)()(,,...,)n z z z z =为数列(1)x 的均值生成数列，即
(1)(1)(1)
()()(1)0.50.5k k k z x x -=+
于是定义GM(1，1)的灰微分方程模型为
(1)()()k d k az b +=
即或：
(0)(1)()()k k x az b += （1）
注：(0)
()k x 为灰导数，a 为发展系数，(1)
()k z 为白化背景值，b 为灰作用量。
将时刻k=2,3,…,n 带入（1）得
引入矩阵向量记号，
()
,T
u a b =，()
(0)(0)
(0)
(2)(3)(),,...,T
n Y x x x =，
(1)(2)(1)()
11n z B z ⎛⎫
- ⎪= ⎪ ⎪-⎝⎭
⇒Y=Bu (GM(1,1))
运行程序得
3.对应的白化模型为：
(1)()
(1)()t t dx ax b dt
+=
对应的解为：(1)
(0)(1)()
(1)
()a t t b b
x x e a a
--=-+
得到预测值：
(1)
(0)(1)
(1)
()ak k b b
x
x e a a
-+=-+，k=1,2,…,n-1.
⇒
(0)(1)(1)(1)(1)()k k k x x x ++=-,k=1,2,…,n-1.
得到趋势图，如下：
对NH3-N:(1)
0.948(1)
0.948(2)
()
0.339()t t t x e
e
--=- （t>=2）.
2
3
4
5
6
7
8
9
10
-200
02004006008001000
1200
逐年排污量
对于其它量：
对于pH:
而河水的正常pH值为6.5-8.5，故长江的pH是正常的；
对于DO:
江水的DO略有上升，但仍然在可控范围内。
C0DMn:
对于
五、分析与解决
长江中的氮污染日益严重，由图像可知，如果不采取措施，不出几年，长江的氮污染将会超出人们的控制，长江作为中国第二大河，对我国经济，农业等有着重大的意义。
一旦污染，后果不堪设想，因此我们要注重保护长江免于氮污染的危害。
有以下几种解决方案：
附录：
1.数据
2.判断
#include<stdio.h>
#include<math.h>
void main()
{
float a[3],b[2],t,c;
int i,j;
for(i=0;i<3;i++)
scanf("%f",&a[i]);
for(i=1,j=0;i<3&&j<2;i++,j++)
{
b[j]=a[i-1]/a[i];
printf("b[%d]=%f\t",j+1,b[j]);
}
printf("\n");
if(b[0]>exp(-0.5)&&b[0]<exp(0.5)&&b[1]>exp(-0.5)&&b[1]<exp(0.5)) printf("满足可容覆盖区间\n");
else
{
printf("不满足可容覆盖区间\n");
t=a[0];
while(b[0]<exp(-0.5)||b[0]>exp(0.5)||b[1]<exp(-0.5)||b[1]>exp(0.5)) {
for(i=0;i<3;i++)
a[i]+=0.001;
for(i=1,j=0;i<3&&j<2;i++,j++)
b[j]=a[i-1]/a[i];
}
c=a[0]-t;
printf("改变后\n");
for(i=0;i<3;i++)
printf("a[%d]=%f\t",i+1,a[i]);
printf("\n");
for(j=0;j<2;j++)
printf("b[%d]=%f\t",j+1,b[j]);
printf("\n");
printf("c=%f",c);
printf("\n");
}
}
3.求a,b
#include<stdio.h>
void main()
{
float a[3];
float b[2];
int i,j;
float c,d,e,f,m,n;
for(i=0;i<3;i++)
scanf("%f",&a[i]);
for(i=1,j=0;i<3&&j<2;i++&&j++)
b[j]=0.5*a[i]+0.5*a[i-1];
c=b[0]+b[1];
d=a[1]+a[2];
e=b[0]*a[1]+b[1]*a[2];
f=b[0]*b[0]+b[1]*b[1];
m=1.0*(c*d-2*e)/(2*f-c*c);
n=1.0*(d*f-c*e)/(2*f-c*c);
printf("a=%f,b=%f",m,n);
printf ("\n");
for (j=0;j<2;j++)
printf ("b[%d]=%f\t",j+1,b[j]); }。