黄河小浪底调水调沙.doc

数值分析课程论文

论文题目：
指导老师：
学 院：
专 业：
姓 名：
学 号：

【实验课题】 黄河小浪底调水调沙问题

【实验目标】
（1）加深对插值及数据拟合知识的理解；
（2）学会利用拟合实现计算有关数值方法；
（3）验证插值拟合所预言的数值现象；
（4）改进曲线拟合既有算法；
（5）掌握最小二乘法的基本原理，并会通过计算机解决实际问题。
【理论概述与算法描述】
为了确定排沙量与时间，排沙量与水流量的函数关系，我们需要对数据进行曲线拟合，所以通过Matlab对数据进行插值拟合，提高精确度，使图像变得光滑，然后利用多项式进行拟合。当多项式次数越高拟合也越准确，但是数据受到的影响较多，所以这里的数据也不是准确值，因此我们只取三次进行拟合，也方便了后续的计算。
符号说明
t: 时间或时间点
v: 水流量
S: 含沙量
V: 排沙量
【实验问题】
在小浪底水库蓄水后，黄河水利委员会进行了多次试验，特别是2004年6月至7月进行的黄河第3次调水调沙试验具有典型意义．这次试验首次由小浪底、三门峡和万家寨三大水库联合调度，进行接力式防洪预泄放水，形成人造洪峰进行调沙试验获得成功．这次调水调沙试验的一个重要目的就是由小浪底上游的三门峡和万家寨水库泄洪，在小浪底形成人造洪峰，冲刷小浪底库区沉积的泥沙．在小浪底水库开闸泄洪以后，从6月27日开始三门峡水库和万家寨水库陆续开闸放水，人造洪峰于29日先后到达小浪底，7月3日达到最大流量2720 m3／s，使小浪底水库的排沙量也不断地增加．表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
2000
1850
1820
1800
1750
1500
1000
900
含沙量
118
120
118
105
80
60
50
30
26
20
8
5
试根据实验数据建立模型解决下面问题：
(1) 给出估算任意时刻的排沙量及总排沙量的方法；
(2) 确定排沙量与水流量的变化关系。
【实验过程与结果】
(1) 给出估算任意时刻的排沙量及总排沙量的方法

通过Matlab工具将所知道的数据显示为直观的图像，如下所示，具体程序见附录的%tuxing.m。

通过观察图像，看出其变化并不光滑，而且也没有特定的表现出服从某种分布的趋势。
于是分别对含沙量和水流量进行插值拟合，便可以得到下面图像和结果，具体程序见附录%hansha.m和%liuliang.m

所得到的拟合函数为： y = 0.014*x^{3} - 1.3*x^{2} + 21*x + 16
即含沙量与时间的关系式为：S=0.014*t^3-1.3*t^2+21*t+16

所得到的拟合函数为：y = 0.13*x^{3} - 14*x^{2} + 2.4e+002*x + 1.5e+00即水流量与时间的关系式为：v=0.13*t^3-14*t^2+2.4e+002*t+1.5e+003
因为某一时刻的排沙量V=v(t)S(t)，所以将所拟合出来的多项式带入上式，通过Matlab进行计算可以得到下面答案，程序见附录%jisuan.m。
ans=91/50000*t^6-73/200*t^5+2429/100*t^4-14573/25*t^3+2866*t^2+35340*t+24000
即排沙量与时间的关系为：
V=0.0018*t^6-0.365*t^5+24.29*t^4-582.92*t^3+2866*t^2+35340*t+24000
由于这里的多项式次数过高，所以对其再进行一次拟合，结果如下，程序见附录%paisha.m。

拟合后的函数为V= 95*t^3-5.5e+003*t^2+7.7e+004*t-3.2e+004，通过图像可以看出排沙量与时间服从正态分布所以化成形式e的指数形式进行拟合。
得到了拟合函数，下面就可以计算出这几天的总排沙量，通过Matlab编程可以计算出定积分，结果如下，程序详见附录%jisuan.m。

ans =170366976000

即总含沙量为1.704亿吨。

(2) 确定排沙量与水流量的变化关系。
先利用Matlab将排沙量和水流量的相关数据反映到图像中。具体程序 见附录%paishui.m。

通过观察可以看出，其关系是分段的，所以按时间进行分段拟合，拟合原理同问题（1）相同，于是可以得到分段前后的拟合多项

y = - 7.5e-005*x^{3} + 0.43*x^{2} - 5.2e+002*x + 3.6e+004

y = 2.3e-005*x^{3} - 0.066*x^{2} + 1.9e+002*x - 1.9e+005

综上，就可以得到排沙量与水流量的关系式为
- 7.5e-5*v^3+0.43*v^2-5.2e+2*v+3.6e+4 0<=t<9
V=
2.3e-5*v^3-0.066*v^2+1.9e+2*v-1.9e+5 9<=t<=24

【结果分析、讨论与结论】
⑴ 第一个问题，排沙量与时间不是严格的正态函数关系可能与实际有些偏差。
⑵ 第二个问题，由于MATLAB软件的计算可能有些偏差导致拟合的函数关系可能与实际有稍微偏差。
⑶ 另外使用多项式拟合虽然简便实用，但是也存在一定的缺陷，当所用的拟合多项式次数较低以及样本点较少时，拟合程度差，精确度较低。而且选用较大次数多项式进行拟合时，计算量和工作量较差，为了减少误差，我们可以通过增加模型中的拟合多项式的次数，并且加入误差评估，来对模型进行完善。

【附程序】
%tuxing.m
T=1:24;
S=[32 60 75 85 90 98 100 102 108 112 115 116 118 120 118 105 80 60 50 30 26 20 8 5 ];
W=[1800 1900 2100 2200 2300 2400 2500 2600 2650 2700 2720 2650 2600 2500 2300 2200 2000 1850 1820 1800 1750 1500 1000 900];
subplot(2,1,1);
plot(T,S);
hold on;
plot(T,S,'.');
title('时间与含沙量关系');
xlabel('时间t/12h');ylabel('含沙量/公斤每立方米');
subplot(2,1,2);
plot(T,W);
hold on;
plot(T,W,'.');
title('时间与水流量关系');
xlabel('时间t/12h');ylabel('水流量/立方米每秒');

%hansha.m
T=1:24;
S=[32 60 75 85 90 98 100 102 108 112 115 116 118 120 118 105 80 60 50 30 26 20 8 5 ];
x=1:0.1:24;
y=interp1(T,S,x,'spline');
plot(T,S,'.',x,y);
title('时间与含沙量关系拟合图');
xlabel('时间t/12h');ylabel('含沙量/公斤每立方米');

%liuliang.m
T=1:24;
W=[1800 1900 2100 2200 2300 2400 2500 2600 2650 2700 2720 2650 2600 2500 2300 2200 2000 1850 1820 1800 1750 1500 1000 900];
x=1:0.1:24;
y=interp1(T,W,x,'spline');
plot(T,W,'.',x,y);
title('时间与水流量关系拟合图');
xlabel('时间t/12h');ylabel('水流量/立方米每秒');

%jisuan.m
syms t;
S=0.014*t^3-1.3*t^2+21*t+16;
v=0.13*t^3-14*t^2+2.4e+002*t+1.5e+003;
V=v*S;
simple(V);
syms t;
V=95*t^3-5.5e+003*t^2+7.7e+004*t-3.2e+004;
int(12*60*60*V,t,0,24)

%paisha.m
t=1:24;
V=0.0018*t.^6-0.365*t.^5+24.29*t.^4-582.92*t.^3+2866*t.^2+35340*t+24000;
plot(t,V);
title('时间与排沙量关系图')

%paishui.m
t=1:24;
v=0.13*t.^3-14*t.^2+2.4e+002*t+1.5e+003;
V= 95*t.^3-5.5e+003*t.^2+7.7e+004*t-3.2e+004;
plot(v,V,'.');
title('整理图')
figure;
t=1:9;
v=0.13*t.^3-14*t.^2+2.4e+002*t+1.5e+003;
V= 95*t.^3-5.5e+003*t.^2+7.7e+004*t-3.2e+004;
plot(v,V,'.');
title('前半段图')
figure;
t=10:24;
v=0.13*t.^3-14*t.^2+2.4e+002*t+1.5e+003;
V= 95*t.^3-5.5e+003*t.^2+7.7e+004*t-3.2e+004;
plot(v,V,'.');
title('后半段图')