数学建模——水塔流量问题.doc

实验十四　水塔流量问题
【实验目得】
１．了解有关数据处理得基本概念与原理。
2。初步了解处理数据插值与拟合得基本方法，如样条插值、分段插值等。
３。学习掌握用MATLAB命令处理数据插值与拟合问题。
【实验内容】
某居民区有一供居民用水得圆形水塔，一般可以通过测量其水位来估计水得流量.但面临得困难就是,当水塔水位下降到设定得最低水位时,水泵自动启动向水塔供水,到设定得最高水位时停止供水,这段时间就是无法测量水塔得水位与水泵得供水量。通常水泵每天供水一两次,每次约两小时。水塔就是一个高1２、2米、直径17、4米得正圆柱。按照设计,水塔水位降到约8、2米时,水泵自动启动，水位升到约1０、8米时水泵停止工作.
某一天得水位测量记录如表１所示,试估计任何时刻（包括水泵正供水时）从水塔流出得水流量，及一天得总用水量。
表１　水位测量启示录(//表示水泵启动）
时刻（h）
水位（ｃm）
0
968
０、92
９48
1、８4
931
2、95
9１3
3、８7
89８
4、98
８８1
5、90
869
7、01
852
7、９3
8３9
8、97
8２2
时刻（h)
水位（ｃｍ)
9、９8
//
１０、9２
//
10、９5
1082
12、03
10５0
12、95
１０21
13、88
994
14、９8
965
1５、90
941
1６、８３
9１8
17、9３
8９2
时刻(h）
水位(cｍ)
１9、０4
866
１9、9６
84３
20、８４
8２2
22、01
//
２2、９6
/／
23、8８
10５9
２４、99
103５
2５、91
1018

【实验准备】
　　在生产实践与科学研究中，常常遇到这样得问题:由实验或测量得到得一批离散样点,需要确定满足特定要求得曲线或曲面（即变量之间得函数关系或预测样点之外得数据)。如果要求曲线(面）通过所给得所有数据点(即确定一个初等函数通过已知各数据，一般用多项式或分段多项式），这就就是数据插值。在数据较少得情况下，这样做能够取得好得效果.但就是,如果数据较多，那么插值函数就是一个次数很高得函数，比较复杂。如果不要求曲线（面)通过所有得数据点，而就是要求它反映对象整体得变化趋势,可得到更简单实用得近似函数,这就就是数据拟合.函数插值与曲线拟合都就是要根据一组数据构造一个函数作为近似，由于近似得要求不同，二者在数学方法上就是完全不同得.
　1。数据插值得基本方法
　拉格朗日插值
若知道函数＝在互异得两个点与处得函数值与,而想估计该函数在另一点处得函数值，最自然得想法就是作过点（，）与点（,)得直线=，用作为准确值得近似值,如果得到得结果误差太大,还可增加一点得函数值,即已知＝在互异得三个点，与处得函数值，与，可以构造过这三点得二次曲线=，用作为准确值得近似值。
一般得，若已知=在互异得＋１个点,,…，处得函数值，,…，，则可以考虑构造一个过这+１个点得次数不超过得多项式
　 　＝＋+…+＋　　 　 (1）
通过所有+１个点,即满足
　　 　 　　 =,=0,1，…， 　　　　　　 　　　 （2）
然后用作为准确值得近似值。这样构造出来得多项式称为得次拉格朗日插值多项式或插值函数。
　　分段插值
　 多项式历来都被认为就是最好得逼近工具之一，它插值光滑，但不具有收敛性,会随着节点数目增多而次数升高，一般不宜采用高次多项式（如＞７)插值，否则逼近得效果往往就是不理想得,甚至发生龙格振荡（当节点数目不断增大时,在区间中部趋于,但对于区间两端得,并不趋于,也称龙格现象）。
　 在插值范围较小,用低次插值往往就能奏效。最直观得办法就就是将各数据点用折线连接起来，这种增加节点，用分段低次多项式插值得化整为零得处理方法称作分段插值法,即不去寻求整个插值区间上得一个高次多项式，而就是把区间划分为若干个小区间。如果
　 ＝＜<…<= 　 　 　 　 (3）
那么分段线性插值公式为
＝＋，〈≤,＝0，１,…，　 （4）
分段线性插值通常有较好得收敛性与稳定性,算法简单，克服了龙格现象,其缺点就是不如拉格朗日插值多项式光滑。
　 样条插值
分段线性插值函数在节点得一阶导数一般不存在，且不光滑,这就导致了样条插值函数得提出。在机械制造、航海、航空工业中，经常需要解决下列问题:已知一些数据点（,），（,),（,）,如何全部通过这些数据点作一条比较光滑得曲线呢？
绘图员解决了这一问题，首先把数据描绘在平面上，再把一根富有弹性得细直条（称为样条)弯曲，使其一边通过这些数据点,用压铁固定其形状，沿样条边绘出一条光滑得曲线，往往要用几根样条，分段完成上述工作,同时也应让连接点处保持光滑.对绘图员用样条画出得曲线,进行数学模拟，就导出了样条函数得概念。如今已经成为了一个应用极为广泛得数学分支。现在数学上所说得样条,实质上指分段多项式得光滑连接。
　 设有区间[,]得一个划分如(３)式，称分段函数为次样条函数,若它有：
(1）在每个小区间上得次数不超过多项式；
（２)=
　(3）在区间［，］上有—１阶连续得导数；
用样条函数作出得插值称为样条插值,工程上广泛采用三次样条插值。
2.曲线拟合得基本方法
　 曲线拟合问题就是指:已知平面上个点(，），=0，1,…,,互不相同，寻求函数=，使在某种准则下与所有数据点最为接近,即曲线拟合得最好。
　　线性最小二乘法就是解决曲线拟合最常用得方法，其基本思路就是，令
　 　　　 ＝+＋…＋ 　　 　 　(5）
其中就是事先选定得一组函数,系数（＝０,１,…，,＜）待定。寻求，使得残差平方与　　　　 　= 　 　 　 　 　 (6）
达到最小。这里得建模原理实质上与实验七中得回归分析就是一致得。
　 ３。数据插值与拟合得MＡTLAＢ命令
　　对于多项式插值与拟合，有一个方便得方法
p = polyｆit( ｘ , y , m　） 用m次多项式拟合向量数据(x,ｙ),返回多项式(１）得降冥系数.当m≥n-1时,poｌyfｉt实现多项式插值，p返回多项式得系数向量；
y ＝ ｐolyval( p , x )　求polｙfit所得得多项式在x处得插值y，它就是准确值ｆ(x)得近似值；
　 对于一维与二维插值,其命令格式如下
yi　= inteｒp1( ｘ ,　ｙ ， xi ， methｏd ) 求解一维插值问题,x，y分别表示数据点得横、纵坐标向量，且x要单调。xi为插值点,它不能走出x得范围.ｍetｈoｄ为可选参数,对应四种插值方法：
最近邻点插值:'nｅarest’；　 线性插值：＇lｉnear’；（就是缺省值）
三次样条插值：'ｓｐliｎe'； 　　 分段三次插值:＇cuｂic'
ZI =　iｎteｒp2( Ｘ　， Y , Z , XＩ ， YI　, metｈod ）　求解二维插值，X　，Y分别就是m、ｎ维族自变量，均单调递增;Ｚ就是ｍ×n维矩阵，标明相应于所给数据得函数值.向量XI　, YＩ给定网格点得横、纵坐标,ZI返回函数在网格坐标（XI，YI）处得函数值。XI,YI应就是方向不同得向量,即一个就是行向量,一个就是列向量。mｅｔhoｄ得定义与iｎterｐ1命令中得定义就是一致得；
ZI　= gｒidｄaｔa(　x , y , ｚ ,　XＩ ， YI , ｍeｔhod )　插值基点为散乱得节点,各参数定义与二维插值中一致，只不过向量x，y散乱无序,同时mｅｔｈod方法中多了一种MATLAＢ提供得网格插值方法V４;
有关上述命令得详细内容可在MATLAB帮助文件中查阅.
对于线性最小二乘拟合,用得较多得就是多项式拟合,使用前面所讲得pｏｌyｆit命令；若要寻求f（ｘ）任意得非线性函数，则称为非线性最小二乘拟合，MATＬAB提供了两个求解命令：cuｒfit与leａstsq。二者都要事先定义Ｍ-函数文件,但定义方式稍有不同:
c = curvefit(　'fun’ ，　c0 ,　xdatａ ， yｄａtａ，optｉons)　求含参量非线性函数ｆun中得参变量c，使残差平方与（６)最小，xdata,ydatａ为数据点得横、纵坐标向量,ｃ0为参变量得迭代初始值,optｉonｓ见实验一表1（它可以缺省)；非线性函数fun得M－函数文件定义方式为：ｆun（ c ，　xdaｔa）;
c　= leastsq(　’fｕn' , c0 ，　optioｎs） 　求含参量非线性函数fuｎ中得参变量c，使得各数据点函数值fｕｎ得平方与最小,命令中各参数定义与curvefｉt命令中一致,非线性函数fｕn得Ｍ—函数文件定义方式为：fuｎ（　c , xdata　,　yｄata），这里得ｆun已经与数据点得函数值向量ｙdata作差；
【实验方法与步骤】
　　1．引例问题得分析
　　流量就是单位时间流出得水得体积，由于水塔就是圆柱形,横截面积就是常数，在水泵不工作时段,流量很容易从水位对时间得变化率算出,问题就是如何估计水泵供水时段得流量。
　 水泵供水时段得流量只能靠供水时段前后得流量拟合得到,作为拟合得原始数据,我们希望水泵不工作时段得流量越准确越好。我们可以考虑先用表中数据拟合水位～时间函数,然后对之求导即可得到各时段得流量。有了任意时刻得流量,就不难计算一天得总用水量.其实，水泵不工作时段得用水量可以由测量记录直接得到,如由某一时段水位下降量乘以水塔得截面积（水塔截面积就是常数＝(17、4/２)2＝２37、8（m２))就得到这一时段得用水量。这个数值可以还可以用来检验拟合效果。
　 流量就是时间得连续函数，只取决于水位差,与水位本身无关，与水泵就是否工作无关。按照Ｔorricellｉ定律从小孔流出得液体得速度正比于水面高度得平方根,题目给出水塔得最高与最低水位分别为1０、８米与8、2米(设出水口得水位为０），因为＝１、１5，可以忽略水位对流速得影响。简单起见,计算中将流量定义为单位时间流出得水得高度,即水位对时间变化率得绝对值(水位就是下降得）。
　　水泵第1次供水时段为＝9、0到＝１1、0(小时),第2次供水时段为=20、８到=23、０（小时)。这就是根据最高与最低水位分别为10、8米与８、2米,及表1得水位测量记录作出得假设，其中前3个时刻直接取自实测数据（精确到0、１小时）,最后１个时刻来自每次供水约两小时得已知条件（从记录瞧，第2次供水时段应在记录得2２、9６小时之后不久结束)。水泵工作时单位时间得供水量大致为常数，这个常数应该大于单位时间得平均流量。
　　首先考虑拟合水位~时间函数，从表1测量记录瞧,一天有两个供水时段（以下称第1供水时段与第２供水时段)，与三个水泵不工作时段（简称第1时段=0到=8、97，第２时段=10、95到=20、8４，第３时段=23以后）。对第1、２时段得测量数据可直接分别作多项式拟合，得到水位函数。为使拟合曲线比较光滑，多项式次数不要太高，一般在3～6次。由于第3时段只有3个测量记录,无法对这一时段得水位作出较好得拟合.
　 接着确定流量~时间函数,对于第1、２时段只需将水位函数求导数即可，对于两个供水时段得流量，则用供水时段前后(水泵不工作时段）得流量拟合得到,并将拟合得到得第2供水时段流量外推，将第３时段流量包含在第2供水时段内。
　最后一天总用水量等于两个水泵不工作时段与两个供水时段（将第3时段包含在第2供水时段内）用水量之与，它们都可以由流量对时间得积分再乘以水塔截面积得到.
2.ＭATLAＢ命令求解
　拟合第1、2时段得水位,并导出流量，t，h为时刻与水位测量记录（水泵启动得4个时刻不输入），程度代码如下：
>＞　ｔ＝[0　0、9２ １、８4 2、9５　３、8７ 4、9８ 5、90 7、０1 ７、93 8、97 １0、95　1２、０3 12、95　13、８8 14、９8 15、90　１6、83 １7、９3　1９、０4 １9、96 ２0、８4 2３、８８　24、99 25、91];
＞〉　ｈ＝［968 948 931 ９13 8９8 881 869 ８５2　839 ８２2 10８2 1０5０ １02１ ９９4 965 941　918 892 ８6６ ８43 ８22 1059 １035 1018]；
〉〉 c1=poｌyfｉt（t（1：10）,h（1：１0)，3）;％　用3次多项式拟合第1时段得水位
〉〉 a１=poｌyder(c1);%　对拟合得多项式求导数得到第１时段流量
〉＞ ｔp1＝０:0、1：９；％　对第1时段得时刻进行划分
>> x1=abs(pｏlyvaｌ(a1，ｔp１）);% 计算第1时段各时刻得流量
类似地，可计算第2时段各时刻得流量.
>> c2＝pｏlｙfit(t（11:２1），ｈ(11:21）,3）;
〉> a2=pｏlyder(c２）；
〉> tp２=11：0、１:2０、８；
>>　x2=ａbs(polｙvａl（a2,tp２））；
　在第1供水时段（t=９～11）之前(即第1时段)与之后（即第2时段)各取几点,其流量已经得到,用它们拟合水泵第1供水时段得流量。为使流量函数在t=９与t=11连续，我们简单地只取４个点，拟合3次多项式(即曲线必过这4个点)，实现如下：
>> xx1=ａbs（pｏlyval(a１,[8 ９］）)；
〉〉 xx2=ａbs（polyｖaｌ（ａ2，［１1 12］)）；
〉〉　xｘ12＝[xx1,ｘｘ２]；
〉>　c1２=polyfiｔ（［8　9 1１ 12]，xx12,3);%　拟合水泵供水时段得流量函数
＞> tｐ1２＝９：0、1:11；
>〉 x1２＝poｌyvａl（ｃ1２,tp12）； % 计算第1供水时段各时刻得流量
在第２供水时段之前取t＝２0,20、８两点得流水量，第3时段仅有３个水位记录，我们用差分得到流量，然后用这4个数值拟合第２供水时段得流量：
〉〉 dt３=ｄiff（t（２2：２４)）;%　最后3个时刻得两两之差
＞〉 ｄh3=diff(ｈ（2２：2４)）；%　最后3个水位得两两之差
>> dｈｔ3=－dh3、/dt3；％% 用差分计算ｔ（22)与t（２3）得流量
＞> t3=［2０ 20、8 t(２2） ｔ（23)］；％ 取第2时段20，２0、8两点与第3时段23、88,24、99两点
>〉 ｘｘ3=［abs(polyｖａl（a2,ｔ３（1:2)））,ｄht3]; 取第2时段２0,2０、8两点与第３时段23、88，24、９9两点得流量
>> c3＝ｐolyfit（t3,ｘx3，3）%　拟合出第２水泵供水时段得流量函数
>〉 ｔp３=20、8:0、1:24；
>＞ x3=pｏlｙｖａl(ｃ3，tp3);%　输出第２供水时段（外推到t＝２4）各时刻得流量
　 求第1、２时段与第１、2供水时段流量得积分之与，就就是一天总用水量.虽然诸时段得流量已表示为多项式函数，积分可以解析地算出,这里仍可用数值积分计算:
〉>　y１=0、1＊trａpz(x１)%　第1时段用水量,0、1为积分步长
ｙ1 =
　146、18１5
>〉　ｙ２=０、1＊trapz（ｘ2） %　第2时段用水量
y2 ＝
　　 　２５8、0441
〉> y12=０、１＊trapz(x12) %　第１水泵供水时段用水量
y12　=
５0、3９90
>〉 y3＝0、1＊trａpz（x3） %　第2水泵供水时段用水量
ｙ3 =
74、9138
>＞ y＝(y１+y２+y12＋y3）*237、８＊0、01％ 总用水量为水位差乘以水塔截面积,0、01就是因为流量单位为厘米
y =
　 1、2５92e+0０3
【结果分析】
计算出来得各时段用水量可以用测量记录来检验,y1可用第1时段水位测量下降高度为9６8－822＝146来检验,类似地,y2用１0８2－8２2＝260来检验.
供水时段流量得一种检验方法如下：供水时段用水量加上水位上升值260就是该时段泵入得水量，除以时间长度得到水泵得功率(单位时间泵入水量)，而两个供水时段得功率应大致相等。第1、2时段水泵得功率计算如下：
>〉 p1＝(y12+260)／２
ｐ1 =
　 １55、1995
>〉 tp2＝2０、8:0、1：２3；
〉> xp2=poｌyval(c3，ｔp2);
＞〉 ｐ２=(0、1＊ｔｒapz（x3）+260）／２、2
p2 =
152、2335
可以瞧到，两次水泵泵水得功率差别不大。下面就是水塔一天得流量曲线图:

图14、1　当取三次多项式拟合得流量曲线图
由图１4、1我们可以瞧到，流量曲线与原始记录基本上相吻合,但在第1时段与第1泵水时段得交接处曲线不太光滑，这说明我们采用3次曲线通过4点得做法不够好,应该多取几点进行拟合。0点到１0点很流量很低，10点到下午３点即中午时间段就是用水高峰期.
【练习与思考】
1．假定某天得气温变化见下表，试找出这一天得气温变化规律：
时刻(h）

0
1１
12
温度℃（ｔ）

５
１6

8
时刻（h）

8
19
20
２1
２2
23
2４

温度℃（t)
3

２４
２２
20
18
17
1６

2.在化工生产中常常需要知道丙烷在各种温度与压力下得导热系数。下面就是实验得到得一组数据:
（℃）
68
68
８７
87
1

（103）
９、7９８1
13、324
9、０0７8
１3、３5５
9、7９18
14、２７7
9、6563
12、4６３

0、08４8
0、0８9７
0、０7６2
０、0８07
0、0696
０、0753
0、0611
０、06５1
试求＝９９（℃）与＝10、3（1０3）下得导热系数。
３.下表给出了某一海域以码为单位得直角坐标为、得水面一点处以英尺为单位得水深，水深数据就是在低潮时测得得。船得吃水深度为5英尺,问在矩形区域（８52０0）×（－4015０）里,哪些地方船要避免进入。
水道测量数据——在低潮测得得水深

1２9
140
10３、５
88
185、5
195
１05、5
1５7、5
10７、5
７7
８1
16２
162
117、5

7、５
１４1、5
23
1４７
22、5
137、5
８５、5
－6、５
－81
3
56、5
－６6、5
84
－３３、5

4
8
６
8
6
8
8
９
９
8
8
9
４
9