典型数值分析模型

数值分析模型

第一节插值法在生产和科学研究中,经常出现这样的问题:由实验或测量得到的某一函数在一系列点处的值,需要构造一个简单函数作为函数的近似表达式:,使得

这类问题称为插值问题.称为被插值函数,称为插值函数称为插值节点;式(6-1)称为插值条件.常用的插值函数是多项式与样条函数.拉格朗日(lagrange)插值取n次多项式pn(x)作为插值函数

pn(x)=a0+a1x+a2x2+…+anxn

利用插值条件有：其系数行列式为n+1阶范德蒙行列式，因插值节点互不相同，所以方程组的解存在且唯一。

其系数行列式为范德蒙(Vandermonde)行列式:由于插值节点互不相同,所以上述行列式不等于0,故由克莱姆(Cramer)法则知,方程组(6-3)的解存在而且是唯一的.实际上比较简单的方法不是解方程组(6-3),而是构造一组插值基函数.为此,首先求满足条件

的n次多项式因为式（6-4）表明，除点以外，其他所有的节点都是n次多项式的零点，故设

其中A为待定常数。由可得所以

称之为拉格朗日插值基函数。利用插值基函数（6-5），可以构造多项式就是满足插值条件的拉格郎日插值问题的解，称式（6-6）为拉格朗日插值多项式。特别地，当n=1时称为线性插值，其插值多项式为：满足从几何上看,为过两点的直线。当n=2时，称为抛物线插值，其插值多项式为：满足从几何上看为过点和的一条抛物线。埃尔米特插值许多插值问题不但要求在节点上函数值相等，而且还要求对应的导数值也相等，甚至要求高阶导数也相等，满足这种要求的插值多项式被称为埃尔米特（Hermite）插值多项式设在节点上，要求插值多项式，满足条件由于（6-7）式给出了2n+2个条件，因此可以唯一确定一个次数不超过2n+1的多项式，其形式为根据（6-7）式来确定显然非常复杂。仿照拉格朗日插值多项式的基函数方法，可先求插值基函数，及共用2n+2个，每一个基函数都是2n-1次多项式，且满足条件于是满足条件（6-7）的插值多项式可写成由条件（6-8）式显然有

利用拉格朗日插值基函数令其中为（6-5）式所表示的基函数。由条件（6-8）式可得整理得:

解出对两边取对数求导可得于是

同理仿照拉格朗日插值余项的证明方法，若在内的2n+2阶导数存在，则其插值余项为其中三次样条插值分段线性插值，具有良好的稳定性和收敛性，但光滑性较差。在数学上若函数（曲线）的k阶导数存在且连续，则称该曲线具有k阶光滑性。易见，分段线性插值不光滑，这影响了它在某些工程技术实际问题中的应用。例如：在船体、飞机等外形曲线的设计中，不仅要求曲线连续而且还要求曲线的曲率连续，这就要求插值函数具有连续的二阶导数。为解决这一类问题，就产生了三次样条插值。所谓样条（Spline），本来是指一种绘图工具，它是一种富有弹性的细长木条，在飞机或轮船制造过程中，被用于描绘光滑的外形曲线。使用时，用压铁将其固定在一些给定的节点上，在其他地方任其自然弯曲，然后依样画下的光滑曲线，就称为样条曲线。它实际上是由分段三次曲线拼接而成，在连续点即节点上，不仅函数自身是连续的，而且它的一阶和二阶导数也是连续的。从数学上加以概括，可得到样条函数的定义如下：三次样条函数记作S（x），a≤x≤b，满足：①在每个小区间[xi，xi+1]（i=0,1,…,n-1）是三次多项式②在每个内节点xi（i=1,2,…,n-1）上具有二次连续导数③S（xi）=yi，i=0,1,2,…,n由三次样条函数中的条件①知，S（x）有4n个待定系数。由条件②知，S（x）在n-1个内节点上具有二阶连续导数，即满足条件：

共有3（n-1）个条件。由条件③，知S（xi）=yi（i=0,1,2,…,n），共有n+1个条件。因此，要确定一个三次样条，还需要外加4n-3（n-1）-（n+1）=2个条件，最常用的三次样条函数S（x）的边界条件有两类：第一类边界条件：第二类边界条件：特别地，,称为自然边界条件第三类边界条件：

称为周期边界条件。

三次样条插值不仅光滑性好，而且稳定性和收敛性都有保证，具有良好的逼近性质。样条插值函数的建立构造满足条件的三次样条插值函数的表达式可以有多种方法。下面我们利用的二阶导数值表达，由于在区间上是三次多项式，故在上是线性函数，可表示为其中对积分两次并利用及，可定出积分常数，于是得三次样条表达式上式中是未知的，为确定对求导得

由此可得类似地可求出在区间上的表达式，从而得利用可得

其中对第一类边界条件，可导出两个方程

如果令则式（6-14）及其（6-16）可写出矩阵通过求解上述三对角矩阵可求得对于第二类边界条件，直接得端点方程

如果令，则式（6-14）及式（6-18）也可以写成矩阵（6-17）的形式对于第三类边界条件，可得其中则式（6-14）及式（6-19）可以写成矩阵形式求解上述矩阵可得曲线拟合通过实验等方法观测到反映某个函数y=f(x)的数据(xi，yi)（i=1,2,…,n），要求利用这些数据构造出y=f（x）的近似表达式y=P（x），上面介绍的插值法就是寻求近似函数的方法之一。但由于实验观测数据不可避免地带有误差，甚至是较大的误差，所以使用插值法是不合适的，它会保留数据的误差。因此，不必要求近似函数y=P（x）满足

P（xi）=yi（i=1,…,n），而只要求偏差按某种标准最小，以反映所给数据的总体趋势，消除局部波动的影响，这就是曲线拟合问题。这样的函数P（x）称为拟合函数。拟合的准则：衡量一个函数P（x）同所给数据（xi，yi）（i=1,…,n）的偏差

的大小，常用的准则有如下三种：(1)最小二乘准则：使偏差的平方和最小，即(2)最小一乘准则：使偏差的绝对值之和为最小，即

(3)极小极大准则：使偏差的最大绝对值最小，即计算的方法(1)最小二乘准则下的计算方法

设为m个线性无关的函数，对给定的数据（xi，yi）（i=1,2,…,n），求

使

最小。利用极值的必要条件得到关于的线性方程组

则方程组可表示为，其中

由于线性无关，所以G是列满秩，GTG是可逆矩阵，方程组的解存在且唯一，并且常用的拟合曲线:(a)取，得直线拟合(b)取，得多项式拟合

(c)取，得多元线性拟合(d)取，则得曲线拟合

还有许多非线性拟合，例如，S曲线

可通过变量代换，令，，化为对

a1，a2的线性函数。一般地，已知一组数据，先画出数据的散点图，在直观判断的基础上，选几种曲线分别作拟合，选择偏差平方和Q最小的曲线。

(2)最小一乘准则下的计算方法

设为m个线性无关的函数，对给定的数据（xi，yi）（i=1,2,…,n），求

使

达到最小。

易见式中目标函数含有绝对值，为了去掉上式中的绝对值，令

则Ui≥0，Vi≥0，且此时（6-20）式可改写为

相应地（6-21）式可改写为综合以上几个式子，求得的问题可归结为如下线性规划问题：为任意常数(3)极小极大准则下的计算方法

设为m个线性无关的函数，对给定的数据（xi，yi）（i=1,2,…,n），求

使

达到最小。上式可改写为记，则求解的问题可归结为如下线性规划问题：为任意常数第二节非线性方程求根

本节主要讨论单变量非线性方程的求根问题，这里。在科学与工程计算中有大量方程求根问题，其中一类特殊的问题是多项式方程其中系数为实数。方程的根，又称为函数的零点，它使，若可分解为其中为正整数，且。当时，则称单根，若称为重根，或为的重零点。若是的重零点，且充分光滑，则当为代数多项式（6-25）时，根据代数基本定理可知，n次方程在复数域中有且只有n个根（含复根，m重根为m个根），n=1，2时方程的根是大家熟悉的，n=3，4时虽有求根公式但比较复杂，可在数学手册中查到，但已不适合于数值计算，而时就不能用公式表示方程的根。因此，通常对的多项式方程求根与一般连续函数方程（6-20）一样都可采用迭代法求根。迭代发要求先给出根的一个近似，若且，根据连续函数性质可知在内至少有一个实根，这时称为方程（6-20）的有根区间。通常可通过逐次搜索法求得方程的有根区间。例1.求方程的有根区间。解：f(0)<0,f(1)<0,f(2)>0,f(3)>0,f(4)<0,f(5)<0,f(6)>0;由此可知，方程的有根区间为[1,2],[3,4],[5,6];练习:求方程的有根区间。

1．不动点迭代法将方程（6-24）改写成等价的形式若要求满足，则；反之亦然，称为函数的一个不动点。求的零点就等价于求的不动点，选择一个初始近似值，将它代入（6-26）右端，即可求得可以如此反复迭代计算

称为迭代函数。如果对任何，由（6-27）得到的序列有极限则称迭代方程（6-27）收敛，且为的不动点，故称（6-27）为不动点迭代法。上述迭代是一种逐次逼近法，其基本思想是将隐式方程（6-26）归结为一组显式的计算公式（6-27），就是说，迭代过程实质上是一个逐步显式化的过程。方程的求根问题在xy平面上就是在确定曲线与直线y=x交点，对于的某个近似值，在曲线上可确定一点,它以为横坐标，而纵坐标则等于。过引平行x轴的直线，设此直线交直线y=x于点，然后过再作平行于y轴的直线，它与曲线的交点记作，则点的横坐标为，纵坐标则等于，在曲线上得到点列，其横坐标分别为依公式求得的迭代值。如果点列趋向于点，则相应的迭代值收敛到所求的根。

2．不动点的存在性与迭代法的收敛性首先考察在上不动点的存在唯一性。定理1设满足以下两个条件：（1）对任意的有。（2）存在正常，使对任意都有则在上存在唯一的不动点。在的不动点存在唯一的情况下，可得到迭代法（6-27）收敛的一个充分条件。定理2设满足定理1中的两个条件，则对任意，由（6-27）得到的迭代序列收敛到的不动点，并有误差估计迭代过程是个极限过程。在用迭代法进行实际计算时，必须按精度要求控制迭代次数。误差估计式（6-29）原则上可用于确定迭代次数，但它由于含有信息L而不便于实际应用。根据式(6-30)，对任意的正整数p有在上式中令知由此可见，只要相邻两次计算结果的偏差足够小即可保证近似值具有足够的精度。对定理1和定理2中的条件（2），在使用时如果且对任意有则由中值定理可知对有它表明定理中的条件（2）可用（6-30）代替。

3．局部收敛性与收敛阶上面给出了迭代序列在区间上的收敛性，通常称为全局收敛性。有时不易检验定理的条件，实际应用时通常只在不动点的邻近考察其收敛性，即局部收敛性。 定义1设有不动点，如果存在的某个邻域，对任意，迭代（6-26）产生的序列。且收敛到，则称迭代法（6-26）局部收敛。定理3设为的不动点，在的某个邻域连续，且，则迭代法（6-26）局部收敛。为衡量迭代法收敛速度的快慢，我们给出下列定义。定义2设迭代过程收敛于方程的根，如果迭代误差当时成立下列渐进关系式

则称该迭代过程是P阶收敛的，特别地，P=1时称线性收敛，P>1时称超线性收敛，P=2时称平方收敛。定理4对于迭代过程，如果在所求根的邻近连续，并且则该迭代过程在点邻近是p阶收敛的。上述定理告诉我们，迭代过程的收敛速度依赖于迭代函数的选取，如果当时，则该迭代过程只可能是线性收敛。第三节牛顿法及其收敛性

1.牛顿法对于方程如果是线性函数，则它的求根是容易的，牛顿法实质上是一种线性化方法，其基本思想是将非线性方程逐步归结为某种线性方程来求解。设已知方程有近似根（假定），将函数在点展开，有于是方程可近似地表示为

这是个线性方程，记其根为，则的计算公式为这就是牛顿（Newton）法。牛顿法有明显的几何解释，方程的根可解释为曲线与x轴的交点的横坐标，设是根的某个近似值，过曲线上横坐标为的点引切线，并将该切线与x轴的交点的横坐标作为的新的近似值。注意到切线方程为这样求得的值必满足（6-32）从而就是牛顿公式（6-33）的计算结果。由于这种几何背景，牛顿法亦称切线法。关于牛顿法（6-33）的收敛性，可直接由定理4得到，对（6-33）其迭代函数为由于假定是的一个单根，即，则由上式知，于是依据定理4可以断定，牛顿法在根的邻近是平方收敛的。又因，可得

2.简化牛顿法与牛顿下山法牛顿法的优点是收敛快，缺点是每步迭代要计算及，计算量较大且有时计算较困难，二是初始近似只在根附近才能保证收敛，如给的不合适可能不收敛。为克服这两个缺点，通常可用下述方法。

(1)简化牛顿法，也称平行弦法，其迭代公式为迭代函数若，即取。在根附近成立，则迭代法（6-35）局部收敛。在（6-35）中取，则称为简化牛顿法，这类方法计算量省，但只有线性收敛，其几何意义是用平行弦与x轴交点作为的近似。（2）牛顿下山法。牛顿法收敛性依赖初值的选取。如果偏离所求根较远，则牛顿法可能发散。例如，用牛顿法求解方程此方程在附近的一个根。设取迭代初值，用牛顿法公式

计算得迭代三次得到的结果有六位有效数字。但是，如果改用作为迭代初值，则依牛顿法公式（6-37）迭代一次得这个结果反而比更偏离了所求的根为了防止迭代发散，我们对迭代过程再附加一项要求，即具有单调性：满足这项要求的算法称为下山法。我们将牛顿法和下山法结合起来使用，即在下山法保证函数值稳定下降的前提下，用牛顿法加快收敛速度。为此，我们将牛顿法的计算结果与前一步的近似值适当加权平均作为新的改进值其中称为下山因子，（6-39）即为

称为牛顿下山法。选择下山因子时从开始，逐次将减半进行试算，直到能使下降条件（6-38）成立为止。若用此法解方程（6-36），当时由（6-37）求得，它不满足条件（6-38），通过逐次取半进行试算，当时可求得。此时有而，显然。由计算时，均能使条件（6-33）成立。计算结果如下：

即为的近似。一般情况只要能使条件（6-38）成立，则可得到，从而使收敛。第四节弦截法与抛物线法

用牛顿法求方程（6-24）的根，每步除计算外还要算，当函数比较复杂时，计算往往较困难，为此可以利用已求函数值来回避导数值的计算的计算，这类方法是建立插值原理基础上的，下面介绍两种常用方法。

1．弦截法设是的近似根，我们利用构造一次插值多项式，并用的根作为的新的近似根。由于

因此有这样导出的迭代公式（6-42）双点弦截法可以看作牛顿公式中的导数用差商取代的结果。现在解释这种迭代过程的几何意义。如图6-4，曲线上横坐标为的点分别记为，则弦线的斜率等于差商值

其方程是因此，按（6-42）式求得的实际上是弦线与x轴交点的横坐标。这种算法因此而称为弦截法。弦截法与切线法（牛顿法）都是线性化方法，但两者有本质的区别。切线法在计算时只用到前一步的值，而弦截法（6-42），在求时要用到前面两步的结果，因此使用这种方法必须先给出两个初始值。

2．抛物线法设已知方程的三个近似根，我们以这三点为节点构造二次插值多项式。并适当选取的一个零点作为新的近似根，这样确定的迭代过程称抛物线法。亦称密勒（Müller）法。在几何图形上，这种方法的基本思想是用抛物线与x轴交点作为所求根的近似位置（图6-5）。现在推导抛物线法的计算公式。插值多项式

有两个零点式中为了从（6-43）中定出一个值，我们需要讨论根式前正负号的取舍问题。在三个近似根中，自然假定更接近所求的根，这时，为了保证精度，我们选式（6-43）中较接近的值作为新的近似根。为此，只要取根式前的符号与的符号相同。第五节

孩子成长和学生考试成绩问题ti（从11岁起年龄）00.81.42.02.43.24.0增长高度hi（cm）00.742.255.258.2515.0021.38ti（从11岁起年龄）4.85.46.07.08.010.0增长高度hi（cm）26.2528.8830.6032.253335孩子成长问题一个男孩在11岁长到21岁过程中，身高的变化如表5-1所示，试找一个最佳的函数曲线来表示这个男孩的成长过程。

表5-1智商（Iq）105110120116122130114102复习时间（t）10126131682015考试成绩（g）7579688591799876学生考试成绩问题：对8个学生测量其智商Iq和课后复习某门课时间t

及该门课考试成绩g，得表5-2，试研究该门课考试成绩与智商和课后复习时间之间的关系。

表5-2

Mathematica和MATLAB求解

（1）Mathematica命令

利用Mathematica可计算曲线拟合，命令输入格式为

例如，可以利用Mathematica中的Fit命令计算本节中的实际问题。

孩子成长问题的计算：首先利用表5—2中的离散数据，画出散点图。

孩子成长问题ch531

文件名：ch531.ma

d1={{0,0},{0.8,0.74},{1.4,2.25},{2.0,5.25},{2.4,8.25},{3.2,15.00},{4.0,21.38},{4.8,26.25},{5.4,28.88},{6.0,30.60},{7.0,32.25},{8.0,33},{10.0,35}};

gp=ListPlot[d1,PlotStyle->{PointSize[0.01]}]

从图5—3可见，取正弦级数为拟合曲线较为合适，为此令

用Mathematica计算

中的参数a1,a2,a3。

输入命令：

f=Fit[d1,{Sin[Pi*t/20],Sin[3*Pi*t/20],Sin[5*Pi*t/20]},t]

fp=Plot[f,{t,0,10}]

Show[gp,fp]

运行后显示

学生考试成绩问题的计算：最小二乘准则

利用表5-2中的离散数据，在Mathematica中输入：文件名：ch532

d1={{105,10,75},{110,12,79},{120,6,68},{116,13,85},{122,16,91},{130,8,79},{114,20,98},{102,15,76}};

g=Fit[d1,{1,iq,t},{iq,t}]

执行后输出

0.736555+0.473084iq+2.10344t

最小一乘准则

取则线性规划问题为

在Mathematica中输入：

学生成绩问题ch533

文件名：ch533.ma

c=u1+v1+u2+v2+u3+v3+u4+v4+u5+v5+u6+v6+u7+v7+u8+v8;

m={a1+105*a2+10*a3-u1+v1==75,

a1+110*a2+12*a3-u2+v2==79,

a1+120*a2+6*a3-u3+v3==68,

a1+116*a2+13*a3-u4+v4==85,

a1+122*a2+16*a3-u5+v5==91,

a1+130*a2+8*a3-u6+v6==79,

a1+114*a2+20*a3-u7+v7==98,

a1+102*a2+15*a3-u8+v8==76};

ConstrainedMin[c,m,{u1,v1,u2,v2,u3,v3,u4,v4,u5,v5,u6,v6,u7,v7,u8,v8,a1,a2,a3}]//N

执行后输出

{13.9318,{u1->0,v1->2.47727,u2->0,v2->0,u3->2.36364,

v3->0,u4->0,v4->1.25,u5->1.81818,v5->0,u6->0,

v6->0,u7->0,v7->0,u8->6.02273,v8->0,a1->5.59091,

a2->0.431818,a3->2.15909}}

从输出结果易见

g=5.59091+0.431818iq+2.15909t

极小极大准则

取则线性规划问题为：

在Mathematica中输入：

学生成绩问题ch534

文件名：ch534.ma

c=R;

m={a1+105*a2+10*a3-u1+v1==75,

a1+110*a2+12*a3-u2+v2==79,

a1+120*a2+6*a3-u3+v3==68,

a1+116*a2+13*a3-u4+v4==85,

a1+122*a2+16*a3-u5+v5==91,

a1+130*a2+8*a3-u6+v6==79,

a1+114*a2+20*a3-u7+v7==98,

a1+102*a2+15*a3-u8+v8==76,

R-u1-v1>=0,R-u2-v2>=0,

R-u3-v3>=0,R-u4-v4>=0,

R-u5-v5>=0,R-u6-v6>=0,

R-u7-v7>=0,R-u8-v8>=0};ConstrainedMin[c,m,{u1,v1,u2,v2,u3,v3,u4,v4,u5,v5,u6,v6,u7,v7,u8,v8,a1,a2,a3,R}]//N

执行后输出

{3.42529,{u1->0,v1->3.42529,u2->1.07854,v2->2.34674,u3->3.42529,v3->0,u4->0.469349,v4->2.95594,

u5->1.72222,v5->1.70307,u6->2.22031,v6->1.20498,

u7->0,v7->3.42529,u8->3.42529,v8->0,a1->1.86207,a2->0.48659,a3->1.86207,R->3.42529}}

从输出结果易见

g=1.86207+0.48659iq+1.86207t

以上三种准则的计算误差可列表如下：

准则误差平方和误差绝对值之和最大误差最小二乘47.6616.27

4.54最小一乘52.84

13.94

6.02极小极大

57.90

18.49

3.43

第六节估计水塔的水流量

某居民区的民用自来水是由一个圆柱形的水塔提供，一般可以通过测量其水位来估计水的流量。但问题是，当水塔水位下降到设定的最低水位时水泵自动启动向水塔供水，到设定的最高水位时停止供水，在水泵自动加水期间无法测量水塔的水位和水泵的供水量。通常水泵每天供水两次，每次约两小时。水塔是一个高12.2米，直径17.4米的正圆柱，当水塔的水位降至最低水位，约8.2米时，水泵自动启动供水当水塔的水位升高到一个最高水位，约10.8米时，水泵停止工作。表5-4是某一天的水位测量记录数据，测量了28个时刻，但是由于其中有4个时刻遇到水泵正在向水塔供水，而无水位记录，表5-4中用符号∥表示。试估计任何时刻（包括水泵正供水时）的用水率，及一天的总用水量。

模型假设及符号说明模型假设假设水塔中流出的水流量只受社区的日常生活需要的影响，水的消耗每天大致差不多。由Torricelli定律知，从水塔流出的最大流速正比于水位高度的平方根，题目中给出水塔的最高和最低水位分别为10.82米和8.22米，所以对于这两种高度，最大水流速度的比约为

这表明我们可以假设水塔中水位对水流速度影响忽略不计。水泵工作时单位时间的供水量大致为常数，这个常数大于单位时间内从水塔中流出的水流的最大流速，这是因为居民区内一直需要用水，不允许水塔中的水用光。水塔中水流量是时间的连续光滑函数，与水泵工作与否无关。这是因为虽然就个别用户而言可能用水量有较大的变化，但由于个人的用水量与整个居民区用水量相比是非常小的，从统计意义上来讲，不太可能同时整个社区的用水量增长或减少。水泵工作与否完全取决于水塔内水位的高度，且每次加水的工作时间为2小时，根据表6-1中的数据可知，水泵第一次供水时间段为[8.967,10.954]，第二次供水时间段为

[20.839,23.880]符号说明

t：测量的时刻；

h：水位的高度；

v：水塔中水的体积；

f（t）：水塔中水流速度，即流量；

yi（i=1~3）：第i时段用水量，即没有供水时间段用水量；

y12：第1和第2供水时间段用水量之和；

y：一天中总用水量。

问题分析

问题要求任意时刻的用水率，即求单位时间流出的水的体积，一般称为水流速度或流量。由于水塔是一个圆柱体，体积可以很容易地通过水位高度h计算出来，这样在水泵不工作的时间段，水流速度就可以从体积对时间的导数计算出来，由于没有水的体积关于时间的函数表达式，而只能利用问题中给定的原始数据表6-1，通过公式

，计算出离散的在测量时刻的

体积V，因此可以考虑用差商代替微商，也就是用离散代替连续的思想。为提高计算精度，采用二阶差商，即由于所有数据被水泵两次供水分割成三组数据对每组数据的中间数据采用中心差商，前后两个数据不能采用中心差商，改用向前差商、向后差商或用中点公式进行差商。中心差商公式：向前差商公式：向后差商公式：中点公式：以上分析了水泵不工作的时段，用水率的计算。对于水泵供水时段的用水率，计算难度较大，我们只好用供水时间段前后的用水率进行插值或拟合而得到。有了任何时刻的用水率，可以采用数值积分计算一天的总用水量。

模型建立及求解通过以上对问题的分析，现在的问题已转化为根据某一天已测量的时刻水塔中水的流速，产生在整个区间（24小时）上的函数或函数值，一般来说插值和拟合是两种最常用的方法。（1）计算水流速度并画出散点图

MATLAB程序如下：

%估计水塔流量ch541

%计算流速，并画流速散点图

%文件名：ch541.m

t=[0,0.921,1.843,2.949,3.871,4.978,5.900,7.006,7.928,8.967,...

10.954,12.032,12.954,13.875,14.982,15.903,16.826,17.931,...

19.037,19.959,20.839,23.880,24.986,25.908];

h=[9.677,9.479,9.308,9.125,8.982,8.814,8.686,8.525,8.388,...8.220,10.820,10.500,10.210,9.936,9.653,9.409,9.180,...

8.921,8.662,8.433,8.220,10.591,10.354,10.180];

%计算水塔中水的体积

v=(pi*17.4^2/4)*h;

%对每组数据的中间数据计算中心差商

%对每组数据不能计算中心差商的计算向前或向后差商

fori=1:2

f(i)=-(-v(i+2)+4*v(i+1)-3*v(i))/(2*(t(i+1)-t(i)));%计算向前差商

end

fori=3:8

f(i)=-(-v(i+2)+8*v(i+1)-8*v(i-1)+v(i-2))/(12*(t(i+1)-t(i)));%计算中心差商

end

fori=9:10

f(i)=-(3*v(i)-4*v(i-1)+v(i-2))/(2*(t(i)-t(i-1)));%计算向后差商

end

fori=11:12f(i)=-(-v(i+2)+4*v(i+1)-3*v(i))/(2*(t(i+1)-t(i)));%计算向前差商

end

fori=13:19

f(i)=-(-v(i+2)+8*v(i+1)-8*v(i-1)+v(i-2))/(12*(t(i+1)-t(i)));%计算中心差商

end

fori=20:21

f(i)=-(3*v(i)-4*v(i-1)+v(i-2))/(2*(t(i)-t(i-1)));%计算向后差商

end

i=22;f(i)=-(-v(i+2)+4*v(i+1)-3*v(i))/(2*(t(i+1)-t(i)));%计算向前差商

i=23;f(i)=-(v(i+1)-v(i-1))/((t(i+1)-t(i-1)));%用中点方法计算差商

i=24;f(i)=-(3*v(i)-4*v(i-1)+v(i-2))/(2*(t(i)-t(i-1)));%计算向后差商

disp(‘水塔流速’)

f

plot(t,f,‘b*’)

title(‘流速散点图’);xlabel(‘时间（小时）’);ylabel(‘流速（立方米/小时）’);

文件ch541.m执行后，可得水塔中水的流速（经整理成表6—5）以及流速散点图6—8。

表6-5图6-8

（2）模型及计算

通过对不同插值方法的比较，结合假设，考虑到流速应该是时间的连续光滑函数，所以采用三次样条插值模型计算用水率函数

f（x）。

①计算用水率

首先用三次样条插值计算用水率函数f（x），MATLAB程序如下：

估计水塔流量ch542

%文件名：ch542.m

%用三次样条插值计算用水率函数f(t)，并画出流速图

yt=0:1/3600:24;

ys=interp1(t,f,yt,‘spline’);%计算样条插值

plot(yt,ys)

title(‘样条插值下的流速图’);xlabel(‘时间（小时）’);ylabel(‘流速（立方米/小时）’);

执行文件ch542.m后，可得样条插值下的流速图5—6。

②计算一天的总用水量

用三次样条插值模型得到的函数f（x）在时间区间[0，24]上做数值积分可得一天的总用水量。

%用数值积分计算一天总用水量

quat=trapz(ys)*(1/3600);

di