黄河小浪底调水调沙综合项目工程数学实验实验报告.doc

《数学实验》实验报告

题目：黄河小浪底调水调沙工程

姓名： 胡 迪
学号： 14622
专业：信息与计算科学

黄河小浪底调水调沙问题

6月至7月黄河进行了第三次调水调沙实验，特别是初次由小浪底、三门峡和万家寨三大水库联合调度，采用接力式防洪预泄放水，形成人造洪峰进行调沙实验获得成功。整个实验期为20多天，小浪底从6月19日开始预泄放水，至到7月13日恢复正常供水结束。小浪底水利工程按设计拦沙量为75.5亿m3,在这之前，小浪底共积泥沙达14.15亿t。这次调水调沙实验一种重要目就是由小浪底上游三门峡和万家寨水库泄洪，在小浪底形成人造洪峰，冲刷小浪底库区沉积泥沙，在小浪底水库开闸泄洪后来，从6月27日开始三门峡水库和万家寨水库陆续开闸放水，人造洪峰于29日先后到达小浪底，7月3日达到最大流量2700 ，使小浪底水库排沙量也不断地增长。表1是由小浪底观测站从6月29日到7月10日检测到实验数据。
表1 实验观测数据 （ 单位：水流为，含沙量为）
日期
6.29
6.30
7.1
7.2
7.3
7.4
时间
8：00
20:00
8：00
20:00
8：00
20:00
8：00
20:00
8：00
20:00
8：00
20:00
水流量
1800
1900
2100
2200
2300
2400
2500
2600
2650
2700
2720
2650
含沙量
32
60
75
85
90
98
100
102
108
112
115
116

日期
7.5
7.6
7.7
7.8
7.9
7.10
时间
8：00
20:00
8：00
20:00
8：00
20:00
8：00
20:00
8：00
20:00
8：00
20:00
水流量
2600
2500
2300
2200

1850
1820
1800
1750
1500
1000
900
含沙量
118
120
118
105
80
60
50
30
26
20
8
5
注：以上数据重要是依照媒体公开报道成果整顿而成。
当前，依照实验数据建立数学模型研究下面问题：
（1）给出估算任意时刻排沙量及总排沙量办法；
（2）拟定排沙量与水流量变化关系。

核心词：拟合，SAS，Matlab，线性回归，调水调沙实验

问题分析：
1、对于问题一，所给数据中水流量x和含沙量h乘积即为该时刻排沙量y即：y=hx。
2、对于问题二，研究排沙量与排水量关系，从实验数据中可以看出，开始排沙量随水量增长而增长，而后随水流量增长而减少，显然变化关系并非线性关系，为此，把问题分为两某些，从水流量增长到最大值为第一阶段，从水流量最大值到结束为第二阶段，分别来研究水流量与排沙量之间函数关系。
模型假设：
1、水流量和排沙量都是持续，不考虑上游泄洪所带来含沙量和外界带来含沙量。
2、时间是持续变化，所取时间点依次为1,2,3，…,24,单位时间为12h。
模型建立与求解：
<一>对于问题一，由于排沙量与时间散点图基本符合正态曲线，如图二所示。因此，排沙量对数与时间函数关系就应当符合二次函数关系，因而排沙量取对数后，再与时间t进行二次回归，排沙量取自然后数据见表2.
假设排沙量与时间函数关系数学模型是

两边取对数得
Lny=at^2+bt+c
先由表二做出排沙量自然对数lny与时间t散点图见图一，并运用SAS软件进行拟合，得到排沙量自然对数与时间回归方程为：
Lny=-0.0209t^2+0.4298t+10.6321
由回归拟合参数表可知回归方程是明显，由于有关系数人R^2=0.9629,误差均方S^2=0.0543，阐明回归曲线拟合效果较好。
因此排沙量与时间之间函数关系式为

图二：排沙量对时间曲线图

时间点

1
1800
32
57600
10.96128
2
1900
60
114000
11.64395
3
2100
75
157500
11.96718
4
2200
85
187000
12.13886
5
2300
90
207000
12.24047
6
2400
98
235200
12.36819
7
2500
100
250000
12.42922
8
2600
102
265200
12.48824
9
2650
108
286200
12.56445
10
2700
112
302400
12.61951
11
2720
115
312400
12.65332
12
2650
116
307400
12.63591
13
2600
118
306800
12.63395
14
2500
120
300000
12.61154
15
2300
118
271400
12.51135
16
2200
105
231000
12.35017
17

80
160000
11.98293
18
1850
60
111000
11.61729
19
1820
50
91000
11.41864
20
1800
30
54000
10.89674
21
1750
26
45500
10.72547
22
1500
20
30000
10.30895
23
1000
8
8000
8.987197
24
900
5
4500
8.411833

最后对所求出函数关系在区间[0,24]之间进行积分

成果为总排沙量1.93962亿吨，此与媒体报道排沙量几乎同样。

<二>对于第二个问题，两个阶段数据如表三、表四所示

表三：第一阶段实验数据
序号
1
2
3
4
5
6
7
8
水流量x
1800
1900
2100
2200
2300
2400
2500
2600
含沙量h
32
60
75
85
90
98
100
102

表四：第二阶段实验观测数据
序号
1
2
3
4
5
6
7
8
水流量x
2650
2600
2500
2300
2200

1850
1820
含沙量h
116
118
120
118
105
80
60
50

对于第一阶段，有表四用MATLAB作图（如图三）可以看出其变化趋势，咱们用多项式做最小二乘拟合。

设三次拟合函数关系h=a0+a1x+a2x^2+a3x^3
其中a0,a1,a2,a3,为待定系数。
四次拟合函数关系h= a0+a1x+a2x^2+a3x^3+a4x^4
其中a0,a1,a2,a3,a4为待定系数。

图三：第一阶段水流量与排沙量之间关系图

三次多项式拟合由MATLAB拟合函数求解出a0=a1=0,a2=0.0032,a3=-2.4929.则拟合函数
h=0.0032x^2-2.4929x^3,拟合效果如图四所示
图四：三次多项式拟合效果，红线为拟合曲线

类似四次多项式拟合由MATLAB拟合函数求解出a0=a1=a2=0,a3=0.0121,a4= -7.4347则拟合函数
h=0.0121x^3-7.4347x^4,拟合效果如图五所示
图五：四次多项式拟合效果，蓝线线为拟合曲线

对于第二阶段，有表五用MATLAB作图可以看出其变化趋势，咱们用多项式做最小二乘拟合。
设三次拟合函数关系h=a0+a1x+a2x^2+a3x^3
其中a0,a1,a2,a3,为待定系数。
四次拟合函数关系h= a0+a1x+a2x^2+a3x^3+a4x^4
其中a0,a1,a2,a3,a4为待定系数。
三次多项式拟合由MATLAB拟合函数求解出a0=a1=0,a2=-0.9475,a3= 464.9601.则拟合函数
h=-0.9475x^2+464.9601x^3,拟合效果如图图六所示
图六：三次拟合函数拟合效果

类似四次多项式拟合由MATLAB拟合函数求解出a0=a1=0，a2=-0.0013，a3= 1.1219 a4=-354.5952则拟合函数
h=-0.0013x^2+1.1219x^3-354.5952x^4,拟合效果如图七所示
图七：四次拟合函数拟合效果

结论以及分析检查：
<一>用SAS软件做线性回归得到排沙量与时间函数关系式为:

再运用所求函数在区间[0,24]上进行积分得到总排沙量1.93962亿吨，这与现状基本相符。
<二>对于第一阶段三次多项式拟合由MATLAB拟合函数求解a0=a1=0,a2=0.0032,a3=-2.4929则拟合函数h=0.0032x^2-2.4929x^3
对于第一阶段四次多项式拟合由MATLAB拟合函数求解出a0=a1=a2=0,a3=0.0121,a4= -7.4347则拟合函数h=0.0121x^3-7.4347x^4
对于第二阶段三次多项式拟合由MATLAB拟合函数求解出a0=a1=0,a2=-0.9475,a3= 464.9601则拟合函数h=-0.9475x^2+464.9601x^3
对于第二阶段四次多项式拟合由MATLAB拟合函数求解出a0=a1=0，a2=-0.0013，a3= 1.1219 a4=-354.5952则拟合函数h=-0.0013x^2+1.1219x^3-354.5952x^4
讨论与推广：
1、对于第一种问题排沙量与时间不是严格正态函数关系也许与实际有些偏差，此外还可以用SAS软件进行高次多向式回归
2、对于第二个问题，由于MATLAB软件计算也许有些偏差导致拟合函数关系也许与实际有稍微偏差，此外，还可以进行高次拟合。

附录：
1、排沙量与时间关系图像MATLAB程序：
>>t=1:1:24;
>>y=[57600,114000,157500,187000,207000,235200,250000,265200,286,2400,312800,307400,306800,300000,271400,231000,160000,111000,91000,54000,45500,30000,8000,4500];
>> plot(t,y,'r')
2、对排沙量求自然对数MATLAB程序与成果：
y3=log(y)
y3 =
Columns 1 t
hrough 17
10.9613 11.6440 11.9672 12.1389 12.2405 12.3682 12.4292 12.4882 12.5644 12.6195 12.6533 12.6359 12.6340 12.6115 12.5113 12.3502 11.9829
Columns 18 through 24
11.6173 11.4186 10.8967 10.7255 10.3090 8.9872 8.4118
3、第一阶段排沙量与水流量之间关系MATLAB程序：
>>x=[1800,1900,2100,2200,2300,2400,2500,2600,2650,2700,2720];
>> h=[32,60,75,85,90,98,100,102,108,112,115];
>> x1=[2650,2600,2500,2300,2200,,1850,1820,1800,1750,1500,1000,900];
>> h1=[116,118,120,118,105,80,60,50,40,32,20,8,5];
>> plot(x,h,'r:')
4、第一阶段三次多项式拟合函数以及拟合效果程序与成果：
>> A1=polyfit(x,h,3)
>> in polyfit at 80
A1 =
0.0000 -0.0000 0.0032 -2.4929
>> z1=polyval(A1,x);
plot(x,h,'k+',x,z1,'r')
5、第一阶段四次多项式拟合函数以及拟合效果程序与成果：
>>A2=polyfit(x,h,4)
>> In polyfit at 80

A2 =
-0.0000 0.0000 -0.0000 0.0121 -7.4347
>> z2=polyval(A2,x);
>> plot(x,h,'*',x,z2,'r')
6、第二阶段三次多项式拟合函数以及拟合效果程序与成果：
>>A3=polyfit(x1,h1,3)
>> In polyfit at 80
A3 =
-0.0000 0.0006 -0.9475 464.9601
>> z3=polyval(A3,x1);
>> plot(x,h,'*',x1,z3,'b')
7、第二阶段四次多项式拟合函数以及拟合效果程序与成果：
>>A4=polyfit(x1,h1,4)
>> In polyfit at 80
A4 =
-0.0000 0.0000 -0.0013 1.1219 -354.5952
>> z4=polyval(A4,x1);
>> plot(x1,h1,'k*',x1,z4,'r:')

参照文献：
【1】姜启源，数学模型（第三版），高等教诲出版社
【2】楼顺天,MATLAB程序设计语言（第二版），西安电子科技大学出版社
【3】刘娜，在SAS中拟合ARCH/GARCH模型