插值和拟合在水流流量计算中的运用
摘要
在解决实际问题的生产（或工程）实践和科学实验过程中，通常需要通过研究某些变量之间的函数关系来帮助我们认识事物的内在规律和本质属性，而这些变量之间的未知函数关系又常常隐含在从验、观测得到的一组数据之中。因此本文根据实际中一个水流的流量计算问题, 采用插值和拟合的方法对原始数据进行探究,最后得到拟合函数（系数）方程。我们先分别采用Lagrange插值、分段线性插值和三次样条插值等一维插值法分析比较后，得出三次样条插值最合适估计水流的流量；再对结果进行曲线拟合求解，根据最小二乘原理，分别用3次和5次拟合多项式逼近，仔细比较后，3次拟合效果最佳。整个模型求解程序均是在MATLAB7.1中运行的。
关键词: 插值拟合最小二乘法方差分析
目录
一.插值分析法 (1)
1.1拉格朗日多项式插值 (1)
1.2 分段线性插值法 (2)
1.3 三次样条插值 (2)
1.4一维插值法总结 (5)
二. 曲线拟合 (5)
2.1 最小二乘法拟合法 (5)
三．数值实验 (6)
3.1实例测试 (6)
3.2题设假设： (6)
3.3实例解答 (6)
四、模型检验 (13)
五、总结 (14)
参考文献 (14)
附件 (15)
数学软件期末论文
09统计一班 20093370 马宏发
1 一.插值分析法
插值的定义：求过已知有限个数据点的近似函数的方法。
常用的插值法：拉格朗日多项式插值、牛顿插值、分段线性插值、Hermite 插值和三次样条插值。本文主要分析拉格朗日多项式插值、分段线性插值、三次样条插值。
1.1拉格朗日多项式插值：
已知函数()f x 在区间[],a b 上的n+1个结点01,,,n x x x 处的函数为i f(x )(i=0,1,…n)，若存在简单函数()P x ,使得(x )=()i i P f x (i=0,1,…n)成立,就()P x ()f x 关于节点01,,...,n x x x 的插值函数,点01,,...,n x x x 称为差值节点,包含差值节点的区间[a,b]称为插值区间，而()f x 称为被，求插值函数插()P x 的方法称为插值法。若()P x 是一个次数至多为n 次的项式,即
2012()n n n P x a a x a x a x =++++ ,i a 为实数，就称()P x 为插值多项式，相应的插值法称为多项
式插值。若()P x 为分段的多项式,就是分段差值。
n 次多项式有n +1个待定系数，由插值条件(x )=()i i P f x 恰好给出n +1个方程
2010200020112111212......................................................n n n n n n n
n n n n a a x a x a x y a a x a x a x y a a x a x a x y ??++++=??++++=????????++++=?? (1) 先构造一组基函数0110011()...()()..()(),(0,1,....)()...()()..()n j i i n i j i i i i i i n i j
j i x x x x x x x x x x l x i n x x x x x x x x x x -+=-+≠-----===-----∏ （2） 显然()i l x 是n 次多项式且满：0,(),(0,1,....)1,i j j i l x i n j i
≠?==?=?称n 次多项式称n 次多项式()i l x 为节点01,,...,n x x x 上的n 次插值基函数。 设0,01,12,2,000()()()()()()()()()()()(),(0,1,....)n n n
j n n n n n n n i i i i i j i j j i x x P x l x f x l x f x l x f x l x f x f x l x f x i n x x ===≠??- ?=+++=== ?- ???∑∑∏ （3） 多项式（3）称为n 次拉格朗日(lagrange)插值多项式，函数,(),0,1,,k n l x k n = 称为拉格朗日插值基函数。特别地，当n ＝1，2时，n 次拉格朗日(lagrange)插值多项式即为线性差值多项式和抛物插值多项式。
1.2 分段线性插值法
假设区间[a,b]的连续函数g(x)在n+1个节点01...n a x x x b =<<<=上的函数值
(),0,1,...,j i g x y j n ==。则得到xy 平面上的n+1个数据点(,)j j x y 。连接相邻数据点11(,)j j x y --，
(,)j j x y 得到n 条线段，它们组成一条折线。把区间[a,b]上这n 条折线段表示的函数称为被插值函数g(x)关于这n+1个数据点的分段线性插值函数，记作I(x)具有如下性质：1）I(x)可以分段表示，在每个小区间1[,]j j x x -上，它是线性函数，即
1）11
1;11
(),j j j j
j j j j
j j x x x x I x y y x x x x x x x -------=+≤≤--
2）(),0,1,...;j I x y j n ==
3）()I x 在[a,b]上连续。若构造插值函数1
111
11,[,(0)],(),[,()],0i i i i i i i i i i i x x x x x i x x x x l x x x x i n x x 时舍去时舍去，其它---+++-?∈=?-??-=∈=?
-????
因0
()()()()n n
i i i i i i I x f x l x y l x ====∑∑ 则
当g(x)在[a,b]上连续时，分段线性插值函数I(x)具有良好的收敛性，即而且当g(x)在[a,b]上二阶导数连续时，对于任意[,]x a b ∈有
{}{}2
2121()()(),max ,max ()8
j j j n a x b M R x g x I x h h x x M g x -≤≤≤≤''=-≤
=-=
其中用()I x 计算x 点的插值时，只用到x 左右的两个节点，计算量与节点个数n+1无关。但n 越大，分段越多，插值误差越小。实际上用数据点作插值计算时，分段线性插值就足够了，如
数学、物理中用的特殊函数表，数理统计中用的概率分布表等。
1.3 三次样条插值
分段线性插值不光滑，这影响了它在某些工程技术实际问题中的应用。例如：在船体、飞机等外形曲线的设计中，不仅要求曲线连续而且还要求曲线的曲率连续，这就要求插值函数具有连续的二阶导数。为解决这一类问题，就产生了三次样条插值。 从数学上加以概括，可得到样条函数的定义如下： 三次样条函数记作()S x ，a x b ≤≤，满足：
①在每个小区间1[,]i i x x +(0,1,,1)i n =-是三次多项式。 ②在每个内节点i x (0,1,,1)i n =-上具有二次连续导数。
③(),0,1,,i i S x y i n ==
由三次样条函数中的条件①知，()S x 有4n 个待定系数。由条件②知，()S x 在1n -个内节点上具有二阶连续导数，即满足条件：
(0)(0)'(0)'(0),(0,1,2,,1)(0)"(0)
i i i i i i S x S x S x S x i n S x S x -=+??
-=+=-??''-=+?
共有3(1)n -个条件。由条件③，知(),0,1,,i i S x y i n ==，共有1n +个条件。因此，要确定一个三次样条，还需要外加43(1)(1)2n n n ---+=个条件，最常用的三次样条函数()S x 的边界条件有两类：第一类边界条件：00'(),'()n n S x f S x f ''==
第二类边界条件：00''(),''()n n S x f S x f ''''==特别地，0''()0,''()0n S x S x ==,称为自然边界条件。 第三类边界条件： 00()(),'()'(),n n S x S x S x S x ==0()()n S x S x ''''= 称为周期边界条件。
三次样条插值不仅光滑性好，而且稳定性和收敛性都有保证，具有良好的逼近性质。样条插值函数的建立。构造满足条件的三次样条插值函数()S x 的表达式可以有多种方法。
下面我们利用()S x 的二阶导数值"()(0,1,,)j j S x M j n ==表达()S x ，由于()S x 在区间
1[,]j j x x + 上是三次多项式，故"()S x 在1[,]j j x x +上是线性函数，可表示为 11
"()j j j
j j
j
x x x x S x M M h h ++--=+ （5）
其中1j j j h x x +=-对"()S x 积分两次并利用()j j S x y =及11()j j S x y ++=，可定出积分常数，于是得三次样条表达式
3
3221111
1()()()()()
6666j j j j j j j j
j
j j j j
j
j j
x x x x M h x x M h x x S x M M y y h h h h +++++----=++-+- (0,1,,1)j n =- （6） 上式中(0,1,,)j M j n =是未知的，为确定(0,1,,)j M j n =，对()S x 求导得 2
2
1111
()()'()226
j j j j
j j
j
j j j
j
j
x x x x y y M M S x M M h h h h ++++----=-++
-
（7）
由此可得11'(0)3
6
6
j j j j
j j j h h y y S x M M ++-+=--
+
。
同样求出()S x 在区间1[,]j j x x -上的表达式，从而得11111
'(0)6
3
j j j j j j j j h h y y S x M M h -------=
+
+
利用 '(0)'(0)j j S x S x +=-可得 112(1,2,1)j j j j j j M M M d j n μλ-+++==- （8） 其中 111,,0,1,,j j j j j j
j j
h h j n h h h h μλ---=
=
=++；11111[,][,]
6
6[,,]j j j j j j j j j j
f x x f x x d f x x x h h +--+--==+
11()()
[,]j j j j j
f x f x f x x h ++-=
；11111[,][,]
[,,]j j j j j j j j
f x x f x x f x x x h -+--+-=
（9）
对第一类边界条件00'()','()'n n S x f S x f ==，可导出两个方程
010100
1116
2([,]')6
2('[,])n n n n n n M M f x x f h M M f f x x h ---?+=
-??
??
+=-??
（10）
如果令00010061,([,]'),1n d f x x f h λμ==
-=，11
6('[,])n n n n n d f f x x h --=-则式（8）及其（10）可写出矩阵 000
111
1
1
1112
22
2n n n n n
n n M d M d M d M d
λμλμλμ----????
?? ? ? ? ? ?
? ? ?
?=
? ? ? ? ?
? ? ? ??
?????
（11） 通过求解上述三对角矩阵可求得01,,,n M M M 。
对于第二类边界条件，直接得端点方程00","n n M f M f == （12） 如果令00000,2",2"n n d f d f λμ====，则式（8）及式（12）也可以写成矩阵（11）的形式。 对于第三类边界条件，可得 0n M M = 112n n n n n M M M d λμ-++= （13） 其中 011010,1n n n n n n h h h h h h λμλ---=
=-=
++， 01101
[,][,]
6n n n n f x x f x x d h h ---=+ 则式（8）及式（13）可以写成矩阵形式
1111222211112
222n n n n n n n n M d M d M d M d λμμλμλλμ----?????? ??? ? ??? ? ??? ?= ??? ? ??? ? ??? ??????? 求解上述矩阵可得 121,,,,n n M M M M -。
1.4一维插值法总结
拉格朗日插值函数在整个插值区间上有统一的解析表达式，其形式关于节点对称，光滑性好。但缺点同样明显，这主要体现在高次插值收敛性差（龙格现象）；增加节点时前期计算作废，导致计算量大；一个节点函数值的微小变化（观测误差存在）将导致整个区间上插值函数都发生改变，因而稳定性差等几个方面。因此拉格朗日插值法多用于理论分析，在采用拉格朗日插值方法进行插值计算时通常选取n<7。
分段线性插值函数（仅连续）与三次样条插值函数（二阶导数连续）虽然光滑性差，但他们都克服了拉格朗日插值函数的缺点，不仅收敛性、稳定性强，而且方法简单实用，计算量小。因而应用十分广泛。分段线性插值，具有良好的稳定性和收敛性，但光滑性较差。在数学上若函数（曲线）的k 阶导数存在且连续，则称该曲线具有k 阶光滑性。易见，分段线性插值不光滑，这影响了它在某些工程技术实际问题中的应用。
二. 曲线拟合
已知一组数据，即平面上的n 个点(,),=1,2,3,....,i i x y i n ，i x 互不相同，寻求一个函数
（曲线）y = f (x)，使f (x)在某种准则下与所有数据点最为接近，即曲线拟合得最好。
2.1 最小二乘法拟合法
设12(),(),,()m x x x ???为m 个线性无关的函数，对给定的数据(,),(1,2,,)i i x y i n =， 求1122()()()()m m P x a x a x a x ???=+++使22121
1
(,,,)[()]n
n
m i i i i i Q a a a y P x ε====-∑∑最小。
利用极值的必要条件
0(1,2,,)k
Q
k m a ?==?。 得到关于12,,,m a a a 的线性方程组 11
1
211
1
1()(())0
()(())0
()(())0
n
m
i i k k i i k n
m
i i k k i i k n
m
m i i k k i i k x y a x x y a x x y a x ??????======?-=???-=????
?-=??∑∑∑∑∑∑则方程组可表示为T T G Ga G Y =， 其中12(,,,)T m a a a a =，12(,,,)T n Y y y y =，112111222212()
()()()()()()
()()m m n n m n n m
x x x x x x G x x x
????????????
?
?
=
?
?
??
由于1{(),,()}m x x ??线性无关，所以G 是列满秩，T G G 是可逆矩阵，方程组的解存在且唯一，并且1()T T a G G G Y -=。取12()1,(),,()m m x x x x x ???===， 得多项式拟合方程12()m m P x a a x a x =+++。
三．数值实验
为了更准确明白插值分析与拟合方法的实用性，采用以下示例进行测试.
3.1实例测试
许多社区没有测量流入或流出水塔的水量装置，他们只能代之以每小时测量水塔中的水位，其误差不超过0.5%。更重要的是，当水塔中的水位下降到最低水位L时水泵就启动向水塔输水，直到最高水位H，但也不能测量水泵的供水量。因此，当水泵在输水时不容易建立水塔中水位和水泵工作时用水量之间的关系。水泵每天输水一次或两次，每次约二小时。试估计任意时刻（包括水泵在输水工作的时候）从水塔流出的流量f（t），并估计一天的总水量。表1给出了从第一次测量开始的以秒为单位的时刻，以及该时刻的高度单位为百分之一英尺的水位测量值。
（1）影响水箱流量的唯一因素是该区公众对水的普通需要；
（2）水泵的灌水速度为常数；
（3）从水箱中流出水的最大流速小于水泵的灌水速度；
（4）每天的用水量分布都是相似的；
（5）水箱的流水速度可用光滑曲线来近似；
（6）当水箱的水容量达到514×103g 时，开始泵水；达到677.6×103g 时，便停止泵水3.3实例解答
首先依照表1所给数据，用MATLAB作出时间—水位散点图（图1）。
09统计一班 20093370 马宏发
7
图1
下面来计算水箱流量与时间的关系。
根据图1一种简单的处理方法为将表1中的数据分为三段，然后对每一段的数据做如下处理： 设某段数据 ，相邻数据中点的平均流速用下面的公式（流速=（右端点的水位－左端点的水位）/区间长度）：11
1(
)2i i i
i i i
x x y y v x x ++++-=- 每段数据首尾点的流速用下面的公式计算：001220()(34)()v x y y y x x =-+- 122()(34)()n n n n n n v x y y y x x ---=-+--
用以上公式求得时间与流速之间的数据如表2。
0011{(,),(,),,(,)}n n x y x y x y
8
由表2作出时间—流速散点图如图2。
图2
3.3.1插值法分析
由表2，对水泵不工作时段
1，2采取插值方法，可以得到任意时刻的流速，从而可以知道任意时刻的流量。我们分别采取拉格朗日插值法，分段线性插值法及三次样条插值法；对于水泵工作时段1应用前后时期的流速进行插值，由于最后一段水泵不工作时段数据太少，我们将它忽略，只对水泵工作时段2进行插值处理。我们总共需要对四段数据（第1，2未供水时段，第1供水时段，混合时段）进行插值处理，下面以第1未供水时段数据为例分别用三种方法算出流量函数和用水量（用水高度）。 调用附件程序1实现图3结果
第一未供水时间-流速图
图3
运行结果：lglrjf = 145.6250 fdxxjf =147.1469 sancytjf =145.6885
图中曲线lglr 、fdxx 和scyt 分别表示用拉格朗日插值法，分段线性插值法及三次样条插值法得到的曲线。
由表1知，第1未供水时段的总用水高度为146(=968－822)，可见上述三种插值方法计算的结果与实际值（146）相比都比较接近。考虑到三次样条插值方法具有更加良好的性质，建议采取该方法。其他三段的处理方法与第1未供水时段的处理方法类似，这里不再详细叙述，只给出数值结果和函数图像.
第一供水段时间—流速示意图(程序代码见附件程序2)如下
09统计一班 20093370 马宏发
9
16
18202224262830323436第一供水段时间流速图
图4
运行结果：lglrjf =56.4426 fdxxjf =49.6051 sancytjf =53.5903 第2未供水段时间—流速示意图（程序代码见附件程序3）如下
第2未供水段时间—流速示
意图
图5
运行结果：lglrjf = 258.8664 fdxxjf = 258.9697 sancytjf = 258.6547 第2供水段时间—流速示意图（程序代码见附件程序4）如下
1015
20
25
3035
40
45
第2供水段时间—流速示意图
图6
运行结果：lglrjf =104.1526 fdxxjf =73.7929 sancytjf =79.8172
下图是用分段线性及三次样条插值方法得到的整个过程的时间—流速函数示意图。（程序代码见附件程序5）
总时间—流速函数示意图
图7
运行结果：fdxxjf = 534.4311 sancytjf = 541.8148
3
表4是对一天中任取的4个时刻分别用3种方法得到的水塔水流量近似值.
10
09统计一班 20093370 马宏发
11
3.3.2拟合法分析
（1）拟合水位—时间函数
从表1中的测量记录看，一天有两次供水时段和三次未供水时段，分别对第1，2未供水时段的测量数据直接作多项式拟合，可得到水位函数（注意，根据多项式拟合的特点，此处拟合多项式的次数不宜过高，一般以3～6次为宜）。对第3未供水时段来说，数据过少不能得到很好的拟合。
设t ，h 分别为已输入的时刻和水位测量记录（由表1提供，水泵启动的4个时刻不输入），这样第1未供水时段各时刻的水位可由MATLAB 程序完成（程序代码见附件程序6）：
8.2
8.48.68.899.29.49.69.8第1未供水时段的时间—水位拟合函数图
图8
（2）确定流量—时间函数
对于第1，2
未供水时段的流量可直接对水位函数求导（程序代码见附件程序7）：
一未供水
22
24
2628
3032
34
二未供水
图9
运行结果： wgsysl1 =145.6657 wgsysl2 =260.6561 下图为5次多项式拟合，相比较显然较三次拟合的效果好。（程序代码见附件程序8
）：
14
151617
1819202122
23
24一未供水
二未供水
图10
而第1供水时段的流量则用前后时期的流量进行拟合得到。为使流量函数在t=9和t=11连续，我们只取4个点，用三次多项式拟合得到第1供水时段的时间—流量图，可以看到与总时间—流速函数示意图(图7)中的相应部分比较吻合。
1416
18
20
22
24
26
28
30
32
第一供水段流量
图11（程序代码见附件程序9）
运行结果： gsysl1 =49.8215。
在第2供水时段之前取t=19.96,20.84两点的流量，用第3未供水时段的3个记录做差分得到两个流量数据21.62，18.48，然后用这4个数据做三次多项式拟合得到第2供水时段与第3未供水时段的时间—流量图(图7)，可以看到与总时间—流速函数示意图中的相应部分也比较吻合。
1617
18
19
20
21
22
23
24
25
第2供水时段与第3未供水时段时间流速图
图12（程序代码见附件程序10）
运行结果：gsysl2 =73.9635
（3）一天总用水量的估计
分别对供水的两个时段和不供水的两个时段积分（流量对时间）并求和得到一天的总用水量约为530.1068（此数据是总用水高度，单位为cm ）。表6列出各段用水量，与插值法算得的表4相比，二者较为吻合。
（4）结果分析
由表3可以看出，使用三次样条插值法得到的结果(145.6885，258.6547)与表1中记录的下降高度146cm， 260cm相差不大，说明插值结果与原始数据比较吻合。
由表6可以得全天的用水量约为526.7148*233.3475*10=1229075.82升
（5）流量及总用水量的检验
计算出各时刻的流量可用水位记录的数值微分来检验，各时段的用水高度可以用实际记录的水位下降高度来检验。例如，算得第1未供水段的用水量高度是145.67cm，而实际记录的水位下降高度为968－822=146cm，两者是吻合的；算得第2未供水段的用水量高度是260.66cm，而实际记录的水位下降高度为1082－822=260cm，两者也是吻合的。从算法设计和分析可知，计算结果与各时段所用的拟合多项式的次数有关。表7给出的是对第1，2未供水时段用五多项式拟合后得到的用水量结果。
表7
四、模型检验
1、分三段（水泵未工作）的实际用水量与由模型推算的用水量的差异很小，见表5，最
大的不超过4%，三段总计不超过0.5%。其中由模型推算的用水量为水流量函数f (t)的积分；实际用水量为水塔内水的体积之差。
表5 分三段的实际用水量与模型计算出的24小时用水量比较（单位：立方米）
实际用水量模型用水量绝对误差相对误差
第一段[0，8.967] 345.3682 338.0511 7.3171 2.1%
第二段[10.954，20.839] 616.3161 620.8341 4.5180 0.7%
第三段[22.958，25.908] 151.7308 156.6000 4.8692 3.2%
三段总计 1113.4151 1115.4852 2.0701 0.2%
2、两次充水期间，水泵注入量的差异大约为3个立方米，不到0.5%。
水泵充水量=充水后的水量+充水期间的流出量—充水前的水量。
第一次水泵充水量=2564+117.4463-1948=733.4463立方米；
第二次水泵充水量=2564+120.7052-1948=736.7052立方米。
由于水泵功率=充水量/充水时间，由上述计算知，水泵的功率大约为367立方米/小时，而且两次冲水期间计算出的功率大致相等，这与实际问题是一致的。
3、用水高峰的比较
实际与模型之间相差无几。实际用水高峰可近似地用差商最大值点表示为t=11，即上午11点钟左右；模型得到的用水高峰可由模型得到的水流量函数依次在单位区间上积分或者找出水流量函数的最大值点得出，它们都在11点左右。
五、总结
1、优点
（1）模型灵活性好、稳定性强，可用于那些拥有地方性的竖直圆柱型水塔的小城镇和乡镇。模型中的输入数据可以是任何近似均匀的时间间隔时的水位，时间间隔大约2小时。
（2）模型中的数学概念简单，并且容易理解。只用到数值计算知识。
（3）模型容易实现，且给出了一天里水流速度和用水量的精确估计。
2、缺点
（1）本模型受水塔的几何形状限制。
（2）光滑曲线的逼近方法不能模拟真实流动中流速的微小变化，实际流动中流速可能会有一种高程度的“噪音”，即激波出现。
3、改进与推广
（1）我们可以在模型中用一个参数来限定不同几何形状的水塔。
（2）可以通过对流速数据进行回归分析等一系列处理，以便得到一些随机变化的特征。
参考文献
[1] 亨塞尔曼.美.精通Matlab.清华大学出版社.2006-05
[2] 马修斯.美数值方法(MATLAB版)(第4版).电子工业出版社.2005-12
[3] 司守奎.数学建模算法与程序.海军航空工程学院
[4] 陆元鸿.数理统计方法[M] .上海:华东理工大学出版社 .2005
[5] 杨桂元等.数学模型应用实例.合肥工业大学.2007-06
[6]陈桂明，戚红雨，潘伟编著，Matlab 数理统计（6.X），北京：科学出版社，2002
[7]谢云荪，张志让.数学实验.北京：科学出版社.2000
[8]飞思科技产品研发中心编著.MATLAB6.5 辅助优化计算与设计.北京：电子工
业出版社，2003。
附件
程序1：
t = [0,0.46,1.38,2.395,3.41,4.425,5.44,6.45,7.465,8.45,8.97];
v = [29.89,21.74,18.48,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.29];
t0 = 0 : 0.1 : 8.97;
lglr = lglrcz(t, v, t0);
lglrjf = 0.1 * trapz(lglr)
fdxx = interp1(t, v, t0);
fdxxjf = 0.1 * trapz(fdxx)
scyt = interp1(t, v, t0, 'spline');
sancytjf = 0.1 * trapz(scyt)
plot(t, v, '*', t0, lglr, ':', t0, fdxx ,'-.', t0, scyt, 'b')
gtext('lglr')
gtext('fdxx')
gtext('scyt')
程序2：
t = [8.45,8.97,10.95,11.49,12.49];
v = [16.35,19.29,33.50,29.63,31.52];
t0 = 8.97 : 0.1 : 10.95;
lglr = lglrcz(t, v, t0);
lglrjf = 0.1 * trapz(lglr)
fdxx = interp1(t, v, t0);
fdxxjf = 0.1 * trapz(fdxx)
scyt = interp1(t, v, t0, 'spline');
sancytjf = 0.1 * trapz(scyt)
plot(t, v, '*', t0, lglr, ':', t0, fdxx ,'-.', t0, scyt, 'b')
gtext('lglr')
gtext('fdxx')
gtext('scyt')
程序3：
t=[10.95,11.49,12.49,13.42,14.43,15.44,16.37,17.38,18.49,19.50,20.40,20.84]; v=[33.50,29.63,31.52,29.03,26.36,26.09,24.73,23.64,23.42,25.00,23.86,22.17]; t0 = 10.95 : 0.1 : 20.84;
lglr = lglrcz(t, v, t0);
lglrjf = 0.1 * trapz(lglr)
fdxx = interp1(t, v, t0);
fdxxjf = 0.1 * trapz(fdxx)
scyt = interp1(t, v, t0, 'spline');
sancytjf = 0.1 * trapz(scyt)
plot(t, v, '*', t0, lglr, ':', t0, fdxx ,'-.', t0, scyt, 'b')
gtext('lglr')
gtext('fdxx')
gtext('scyt')
程序4：
t = [20.40,20.84,23.88,24.43,25.45,25.91];
v = [23.86,22.17,27.09,21.62,18.48,13.3];
t0 = 20.84 : 0.1 : 23.88;
lglr = lglrcz(t, v, t0);
lglrjf = 0.1 * trapz(lglr)
fdxx = interp1(t, v, t0);
fdxxjf = 0.1 * trapz(fdxx)
scyt = interp1(t, v, t0, 'spline');
sancytjf = 0.1 * trapz(scyt)
plot(t, v, '*', t0, lglr, ':', t0, fdxx ,'-.', t0, scyt, 'b')
gtext('lglr')
gtext('fdxx')
gtext('scyt')
程序5：
t =
[0,0.46,1.38,2.395,3.41,4.425,5.44,6.45,7.465,8.45,8.97,10.95,11.49,12.49,13.42 ,14.43,15.44,16.37,17.38,18.49,19.50,20.40,20.84,23.88,24.43,25.45,25.91];
v =
[29.89,21.74,18.48,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.29,33.50,29.63, 31.52,29.03,26.36,26.09,24.73,23.64,23.42,25.00,23.86,22.17,27.09,21.62,18.48,1 3.30];
t0 = 0:0.1:23.88;
fdxx = interp1(t, v, t0);
fdxxjf = 0.1 * trapz(fdxx)
scyt = interp1(t, v, t0, 'spline');
sancytjf = 0.1 * trapz(scyt)
plot(t, v, '*', t0, fdxx ,'-.', t0, scyt, 'b')
gtext('fdxx')
gtext('scyt')
程序6：
t=[0,0.92,1.84,2.95,3.87,4.98,5.90,7.00,7.93,8.97];
h=[9.68,9.45,9.31,9.13,8.98,8.81,8.69,8.52,8.39,8.22];
c1=polyfit(t(1:10),h(1:10),3);
tp1=0:0.1:8.97;
x1=polyval(c1,tp1);
plot(tp1,x1)
数学软件期末论文
程序7：
t=[0,0.92,1.84,2.95,3.87,4.98,5.90,7.00,7.93,8.97,10.95,12.03,12.95,13.88,14.98 ,15.90,16.83,17.93,19.04,19.96,20.84,23.88,24.99,25.66];
h=[9.68,9.48,9.31,9.13,8.98,8.81,8.69,8.52,8.39,8.22,10.82,10.50,10.21,9.94,9.6 5,9.41,9.18,8.92,8.66,8.43,8.22,10.59,10.35,10.18];
c1=polyfit(t(1:10),h(1:10),3);
c2=polyfit(t(11:21),h(11:21),3);
a1=polyder(c1);
a2=polyder(c2);
tp1=0:0.01:8.97;
tp2=10.95:0.01:20.84;
x13=-polyval(a1,tp1);
x113=-polyval(a1,[0:0.01:8.97]);
wgsysl1=100*trapz(tp1,x113)% */计算第1未供水时段的总用水量/*
x14=-polyval(a1,[7.93,8.97]); % */为下面的程序准备数据/*
x23=-polyval(a2,tp2);
x114=-polyval(a2,[10.95:0.01:20.84]);
wgsysl2=100*trapz(tp2,x114) % */计算第2未供水时段的总用水量/*
x24=-polyval(a2,[10.95,12.03]);% */为下面的程序准备数据/*
x25=-polyval(a2,[19.96,20.84]);% */为下面的程序准备数据/*
subplot(1,2,1)
plot(tp1,x13*100) % */与第1，2未供水时段的时间—流量图单位保持一致/* subplot(1,2,2)
plot(tp2,x23*100) % */与第1，2未供水时段的时间—流量图单位保持一致/*
程序8：
t=[0,0.92,1.84,2.95,3.87,4.98,5.90,7.00,7.93,8.97,10.95,12.03,12.95,13.88,14.98 ,15.90,16.83,17.93,19.04,19.96,20.84,23.88,24.99,25.66];
h=[9.68,9.48,9.31,9.13,8.98,8.81,8.69,8.52,8.39,8.22,10.82,10.50,10.21,9.94,9.6 5,9.41,9.18,8.92,8.66,8.43,8.22,10.59,10.35,10.18];
c1=polyfit(t(1:10),h(1:10),5);
c2=polyfit(t(11:21),h(11:21),5);
a1=polyder(c1);
a2=polyder(c2);
tp1=0:0.01:8.97;
tp2=10.95:0.01:20.84;
x13=-polyval(a1,tp1);
x113=-polyval(a1,[0:0.01:8.97]);
wgsysl1=100*trapz(tp1,x113) % */计算第1未供水时段的总用水量/*
x14=-polyval(a1,[7.93,8.97]); % */为下面的程序准备数据/*
x23=-polyval(a2,tp2);
x114=-polyval(a2,[10.95:0.01:20.84]);
wgsysl2=100*trapz(tp2,x114)% */计算第2未供水时段的总用水量/*
x24=-polyval(a2,[10.95,12.03]);% */为下面的程序准备数据/*
x25=-polyval(a2,[19.96,20.84]);% */为下面的程序准备数据/*
subplot(1,2,1)
09统计一班20093370 马宏发17
plot(tp1,x13*100) % */与第1，2未供水时段的时间—流量图单位保持一致/* subplot(1,2,2)
plot(tp2,x23*100) % */与第1，2未供水时段的时间—流量图单位保持一致/*
程序9：
dygsdsy=[7.93,8.97,10.95,12.03];
dygsdls=[x14,x24];
nhjg=polyfit(dygsdsy,dygsdls,3);
nhsj=7.93:0.1:12.03;
nhlsjg=polyval(nhjg,nhsj);
gssj1=8.97:0.01:10.95;
gs1=polyval(nhjg,[8.97:0.01:10.95]);
gsysl1=100*trapz(gssj1,gs1) % */该语句计算第1供水时段的总用水量/*
plot(nhsj,100*nhlsjg)
程序10：
t3=[19.96,20.84,t(22),t(23)];
ls3=[x25*100,21.62,18.48];
nhhddxsxs=polyfit(t3,ls3,3);
tp3=19.96:0.01:25.91;
xx3=polyval(nhhddxsxs,tp3);
gssj2=20.84:0.01:24;
gs2=polyval(nhhddxsxs,[20.84:0.01:24]);
gsysl2=trapz(gssj2,gs2) % */该语句计算第2供水时段的总用水量/*
plot(tp3,xx3)
18
请教：已知管道直径D，管道内压力P，能否求管道中流体的流速和流量？怎么求 已知管道直径D，管道内压力P，还不能求管道中流体的流速和流量。你设想管道末端有一阀门，并关闭的管内有压力P，可管内流量为零。管内流量不是由管内压力决定，而是由管内沿途压力下降坡度决定的。所以一定要说明管道的长度和管道两端的压力差是多少才能求管道的流速和流量。 对于有压管流，计算步骤如下： 1、计算管道的比阻S，如果是旧铸铁管或旧钢管，可用舍维列夫公式计算管道比阻s=0.001736/d^5.3 或用s=10.3n2/d^5.33计算，或查有关表格； 2、确定管道两端的作用水头差H=P/(ρg)，)，H 以m为单位；P为管道两端的压强差(不是某一断面的压强），P以Pa为单位； 3、计算流量Q：Q = (H/sL)^(1/2) 4、流速V=4Q/(3.1416d^2) 式中：Q―― 流量，以m^3/s为单位；H――管道起端与末端的水头差，以m^为单位；L――管道起端至末端的长度，以m为单位。 管道中流量与压力的关系 管道中流速、流量与压力的关系 流速:V=C√(RJ)=C√[PR/(ρgL)] 流量：Q=CA√(RJ)=√[P/(ρgSL)] 式中：C――管道的谢才系数；L――管道长度；P――管道两端的压力差；R――管道的水力半径；ρ――液体密度；g――重力加速度；S――管道的摩阻。 管道的内径和压力流量的关系 似呼题目表达的意思是：压力损失与管道内径、流量之间的关系，如果是这个问题，则正确的答案应该是：压力损失与流量的平方成正比，与内径5.33方成反比，即流量越大压力损失越大，管径越大压力损失越小，其定量关系可用下式表示： 压力损失（水头损失）公式（阻力平方区） h=10.3*n^2 * L* Q^2/d^5.33 上式严格说是水头损失公式，水头损失乘以流体重度后才是压力损失。式中n――管内壁粗糙度；L――管长；Q――流量；d――管内径 在已知水管:管道压力0.3Mp、管道长度330、管道口径200、怎么算出流速与每小时流量？ 管道压力0.3Mp、如把阀门关了，水流速与流量均为零。（应提允许压力降） 管道长度330、管道口径200、缺小单位，管道长度330米？管道内径200为毫米？其中有无阀门与弯头，包括其形状与形式。 水管道是钢是铸铁等其他材料，其内壁光滑程度不一样。 所以无法计算。 如果是工程上大概数，则工程中水平均流速大约在0.5－－1米/秒左右，则每小时的流量为：0.2×0.2×0.785×1（米/秒，设定值）×3600＝113（立方/小时） 管道每米的压力降可按下式计算：
插值法和拟合实验报告 一、 实验目的 1.通过进行不同类型的插值，比较各种插值的效果，明确各种插值的优越性； 2.通过比较不同次数的多项式拟合效果，了解多项式拟合的原理； 3.利用matlab 编程，学会matlab 命令； 4.掌握拉格朗日插值法； 5.掌握多项式拟合的特点和方法。 二、 实验题目 1.、插值法实验 将区间[-5,5]10等分，对下列函数分别计算插值节点 k x 的值，进行不同类型 的插值，作出插值函数的图形并与)(x f y =的图形进行比较： ;11)(2x x f += ;a r c t a n )(x x f = .1)(42 x x x f += (1) 做拉格朗日插值； (2) 做分段线性插值； (3) 做三次样条插值. 2、拟合实验 给定数据点如下表所示： 分别对上述数据作三次多项式和五次多项式拟合，并求平方误差，作出离散函数 ),(i i y x 和拟合函数的图形。 三、 实验原理 1.、插值法实验
∏∑∏∏∏∑∑≠==≠=≠=≠=+-==--= =-= ==-=-=----==++==j i j j i i i i i n i i n n j i j j n j i j j i i n j i j j n i i i n i i n n n o i n i i n x x x x x y x l x L x x c n i x x c x x x c x x x x x x x x c y x l x L y x l y x l y x l x L ,00 ,0,0,01100 00 )(l )()() (1 ,1,0, 1)()(l ) ()())(()()()()()()()(， 故， 得 再由，设 2、拟合实验
学生实验报告
了解插值与拟合的基本原理和方法；掌握用MATLAB计算插值与作最小二乘多项式拟合和曲线拟合的方法；通过范例展现求解实际问题的初步建模过程； 通过动手作实验学习如何用插值与拟合方法解决实际问题，提高探索和解决问题的能力。这对于学生深入理解数学概念，掌握数学的思维方法，熟悉处理大量的工程计算问题的方法具有十分重要的意义。 二、实验仪器、设备或软件：电脑,MATLAB软件 三、实验内容 1．编写插值方法的函数M文件； 2．用MATLAB中的函数作函数的拟合图形； 3．针对实际问题，试建立数学模型，并求解。 四、实验步骤 1．开启软件平台——MATLAB，开启MATLAB编辑窗口； 2．根据各种数值解法步骤编写M文件； 3．保存文件并运行； 4．观察运行结果(数值或图形)； 5．写出实验报告，并浅谈学习心得体会。 五、实验要求与任务 根据实验内容和步骤，完成以下具体实验，要求写出实验报告（实验目的→问题→数学模型→算法与编程→计算结果→分析、检验和结论→心得体会）。 1．天文学家在1914年8月的7次观测中，测得地球与金星之间距离（单位：米），并取得常用对数值，与日期的一组历史数据如下表： 由此推断何时金星与地球的距离（米）的对数值为9.93518？ 解：输入命令
days=[18 20 22 24 26 28 30]; distancelogs=[9.96177 9.95436 9.94681 9.93910 9.93122 9.92319 9.91499]; t1=interp1(distancelogs,days,9.93518) %线性插值 t2=interp1(distancelogs,days,9.93518,'nearest') %最近邻点插值 t3=interp1(distancelogs,days,9.93518,'spline') %三次样条插值 t4=interp1(distancelogs,days,9.93518,'cubic') %三次插值 计算结果： t1 = 24.9949 t2 = 24 t3 = 25.0000 t4 =
流量与管径、压力、流速的一般关系2007年03月16日星期五13:21 一般工程上计算时，水管路，压力常见为0.1--0.6MPa，水在水管中流速在1--3米/秒，常取1.5米/秒。流量=管截面积X流速=0.002827X管内径的平方X流速(立方米/小时)。其中，管内径单位：mm ，流速单位：米/秒，饱和蒸汽的公式与水相同，只是流速一般取20--40米/秒。水头损失计算Chezy 公式 Chezy 这里：Q ——断面水流量（m3/s）C ——Chezy糙率系数（m1/2 /s）A ——断面面积（m2）R ——水力半径（m）S ——水力坡度（m/m）根据需要也可以变换为其它表示方法: Darcy-Weisbach公式 由于这里：hf ——沿程水头损失（mm3/s） f ——Darcy-Weisbach水头损失系数（无量纲）l ——管道长度（m）d ——管道内径（mm）v ——管道流速（m/s）g ——重力加速度（m/s2） 2 水力计算是输配水管道设计的核心，其实质就是在保证用户水量、水压安全的条件下，通过水力计算优化设计方案，选择合适的管材和确经济管径。输配水管道水力计算包含沿程水头损失和局部水头损失，而局部水头损失一般仅为沿程水头损失的5~10%，因此本文主要研究、探讨管道沿程水头损失的计算方法。 1.1 管道常用沿程水头损失计算公式及适用条件管道沿程水头损失是水流摩阻做功消耗的能量，不同的水流流态，遵循不同的规律，计算方法也不一样。输配水管道水流流态都处在紊流区，紊流区水流的阻力是水的粘滞力及水流速度与压强脉动的结果。紊流又根据阻力特征划分为水力光滑区、过渡区、粗糙区。管道沿程水头损失计算公式都有适用范围和条件，一般都以水流阻力特征区划分。水流阻力特征区的判别方法，工程设计宜采用数值做为判别式，目前国内管道经常采用的沿程水头损失水力计算公式及相应的摩阻力系数，按照水流阻力特征区划分如表1。 沿程水头损失水力计算公式和摩阻系数表1 阻力特征区适用条件水力公式、摩阻系数符号意义水力光滑区>10 雷诺数h：管道沿程水头损失v：平均流速d：管道内径γ：水的运动粘滞系数λ：沿程摩阻系数Δ：管道当量粗糙度q：管道流量Ch：海曾-威廉系数C：谢才系数R：水力半径n：粗糙系数i：水力坡降l：管道计算长度紊流过渡区10< <500 （1）（2）紊流粗糙区>500 达西公式是管道沿程水力计算基本公式，是一个半理论半经验的计算通式，它适用于流态的不同区间，其中摩阻系数λ可采用柯列布鲁克公式计算，克列布鲁克公式考虑的因素多，适用范围广泛，被认为紊流区λ的综合计算公式。利用达西公式和柯列布鲁克公式组合进行管道沿程水头损失计算精度高，但计算方法麻烦，习惯上多用在紊流的阻力过渡区。海曾—威廉公式适用紊流过渡区，其中水头损失与流速的 1.852次方成比例（过渡区水头损失h∝V1.75~2.0）。该式计算方法简捷，在美国做为给水系统配水管道水力计算的标准式，在欧洲与日本广泛应用，近几年我国也普遍用做配水管网的水力计算。谢才公式也应是管道沿程水头损失通式，且在我国应用时间久、范围广，积累了较多的工程资料。但由于谢才系数C采用巴甫洛夫公式或曼宁公式计算确定，而这两个公式只适用于紊流的阻力粗糙区，因此谢才公式也仅用在阻力粗糙区。另外舍维列夫公式，前一段时期也广泛的用做给水管道水力计算，但该公式是由旧钢管和旧铸铁管 3 管材试验资料确定的。而现在国内采用的金属管道已普遍采用水泥砂浆和涂料做内衬，条件已发生变化，因此舍维列夫公式也基本不再采用。 1.2 输配水管道沿程水头损计算的实用公式输配水管道沿程水头计算时，先采用判别水流的阻力特征用，再选择相应的公式计算，科学合理，但操作麻烦，特别在流速是待求的未知数时，需要采用试算的方法确定雷诺数（Re）很不方便。为了使输配水管道水力计算能满足工程设计的需要，又可以方便的选择计算公式和进行简捷的计算，根据多年来管道水力计算的经验，《室外给水设计规
利用matlab实现插值与拟合实验 张体强1026222 张影 晁亚敏 [摘要]：在测绘学中，无论是图形处理，还是地形图处理等，大多离不开插值与拟合的应用，根据插值与拟合原理，构造出插值和拟合函数，理解其原理，并在matlab平台下，实现一维插值，二维插值运算，实现多项式拟合，非线性拟合等，并在此基础上，联系自己所学专业，分析其生活中特殊例子，提出问题，建立模型，编写程序，以至于深刻理解插值与拟合的作用。 [关键字]: 测绘学插值多项式拟合非线性拟合 [ Abstract]: in surveying and mapping, whether the graphics processing, or topographic map processing and so on, are inseparable from the interpolation and fitting application, according to the interpolation and fitting theory, construct the fitting and interpolation function, understanding its principle, and MATLAB platform, achieve one-dimensional interpolation, two-dimensional interpolation, polynomial fitting, non-linear fitting, and on this basis, to contact their studies, analysis of their living in a special example, put forward the question, modeling, programming, so that a deep understanding of interpolation and fitting function. [ Key words]: Surveying and mapping interpolation polynomial fitting nonlinear
本节概要 本节讨论喷管内流量、流速的计算。工程上通常依据已知工质初态参数和背压，即喷管出口截面处的工作压力，并在给定的流量等条件下进行喷管设计计算，以选择喷管的外形及确定其几何尺寸；有时也需就已有的喷管进行校核计算，此时喷管的外形和尺寸已定，须计算在不同条件下喷管的出口流速及流量。在喷管的计算中要注意到背压对确定喷管出口截面上压力的作用。 本节内容 4.8.1 流速计算及其分析 4.8.2 临界压力比 4.8.3 流量计算及分析 4.8.4 例题 本节习题 4-24、4-25、4-26、4-27、4-29 下一节 流速计算及其分析 1.喷管出口截面的流速计算 2.压力比对流速的影响 …喷管出口截面的流速计算 据能量方程,气体在喷管中绝热流动时任一截面上的流速可由下式计算： (4-28) 因此，出口截面上流速： (4-28a) 或(4-28b)
在入口速度较小时，上式中可忽略不计，于是： (4-28c) (4-28)各式表明，气流的出口流速取决于气流在喷管中的绝热焓降。值得注意的是，上述各式中焓的单位是J/kg。 如果理想气体可逆绝热流经喷管，可据初态参数(p1,T1)及速度求取滞止参数， 然后结合出口截面参数如p2按可逆绝热过程方程式求出T2从而计算h2再求得；对水蒸汽 可逆绝热流经喷管，可以利用h-s图，根据进口蒸汽的状态查得初态点1，通过点1作垂线与喷管出口截面上压力p2相交，得出状态点2，从点1和2可查出h1和h2，代入式(4-28)即可求出出口流速。 ☆ 式子对理想气体和实际气体均适用；与过程是否可逆无关，但不可逆绝 热流动，若用可逆的关系求出h2在求得的需修正，若h2是不可逆过程终态的焓，则求出的不需修正。 式的适用范围是什么？是否与过程的可逆与否有关？与工质的性质有关？ 返回
实验三 插值法与拟合实验 一、实验目的 1. 通过本实验学会利用程序画出插值函数，并和原图形相比较 2. 通过本实验学会拟合函数图形的画法，并会求平方误差 二、实验题目 1. 插值效果的比较 实验题目：区间[]5,5-10等分，对下列函数分别计算插值节点k x 的值，进行不同类型的插值，作出插值函数的图形并与)(x f y =的图形进行比较： 2 11)(x x f +=； x x f arctan )(=； 4 41)(x x x f += （1） 做拉格朗日插值； （2） 做三次样条插值. 2. 拟合多项式实验 实验题目：给定数据点如下表所示： 分别对上述数据作三次多项式和五次多项式拟合，并求平方误差，作出离散函数),(i i y x 和拟合函数的图形. 三、实验原理 本实验应用了拉格朗日插值程序、三次样条插值程序、多项式拟合程序等实验原理. 四、实验内容 1（1） figure x=-5:0.2:5; y=1./(1+x.^2); plot(x,y,'r'); hold on %拉格朗日插值 x1=-5:1:5; y1=1./(1+x1.^2); xx=-4.5:0.5:4.5; yy=malagr(x1,y1,xx); plot(xx,yy,'+') %三次样条插值 dy0=1./(1+25); dyn=1./(1+25);
m=maspline(x1,y1,dy0,dyn,xx); plot(xx,m,'ok') 1（2） x=-5:0.2:5; y=atan(x); plot(x,y,'r'); hold on %拉格朗日插值 x1=-5:1:5; y1=atan(x1); xx=-4.5:0.5:4.5; yy=malagr(x1,y1,xx); plot(xx,yy,'+') %三次样条插值 dy0=1./(1+25); dyn=1./(1+25); m=maspline(x1,y1,dy0,dyn,xx); plot(xx,m,'ok') 1（3） x=-5:0.2:5; y=x.^2./(1+x.^4); plot(x,y,'r'); hold on %拉格朗日插值 x1=-5:1:5; y1=x1.^2./(1+x1.^4); xx=-4.5:0.5:4.5; yy=malagr(x1,y1,xx); plot(xx,yy,'+') %三次样条插值 dy0=1./(1+25); dyn=1./(1+25); m=maspline(x1,y1,dy0,dyn,xx); plot(xx,m,'ok') 2. x=[-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5]'; y=[-4.45 -0.45 0.55 0.05 -0.44 0.54 4.55]'; plot(x,y,'or'); hold on %三次多项式拟合 p1=mafit(x,y,3);
MATLAB中的插值、拟合与查表 插值法是实用的数值方法，是函数逼近的重要方法。在生产和科学实验中，自变量x与因变量y的函数y = f(x)的关系式有时不能直接写出表达式，而只能得到函数在若干个点的函数值或导数值。当要求知道观测点之外的函数值时，需要估计函数值在该点的值。 如何根据观测点的值，构造一个比较简单的函数y=φ(x)，使函数在观测点的值等于已知的数值或导数值。用简单函数y=φ(x)在点x处的值来估计未知函数y=f(x)在x点的值。寻找这样的函数φ(x)，办法是很多的。φ(x)可以是一个代数多项式，或是三角多项式，也可以是有理分式；φ(x)可以是任意光滑（任意阶导数连续）的函数或是分段函数。函数类的不同，自然地有不同的逼近效果。在许多应用中，通常要用一个解析函数（一、二元函数）来描述观测数据。 根据测量数据的类型： 1．测量值是准确的，没有误差。 2．测量值与真实值有误差。 这时对应地有两种处理观测数据方法： 1．插值或曲线拟合。 2．回归分析（假定数据测量是精确时，一般用插值法，否则用曲线拟合）。 MATLAB中提供了众多的数据处理命令。有插值命令，有拟合命令，有查表命令。 2.2.1 插值命令 命令1 interp1 功能一维数据插值（表格查找）。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。各个参量之间的关系示意图为图2-14。 格式 yi = interp1(x,Y,xi) %返回插值向量yi，每一元素对应于参量xi，同时由向量x 与Y的内插值决定。参量x指定数据Y的点。若Y为一矩阵，则按Y的每列计算。yi是阶数为length(xi)*size(Y,2)的输出矩阵。 yi = interp1(Y,xi) %假定x=1:N，其中N为向量Y的长度，或者为矩阵Y的行数。 yi = interp1(x,Y,xi,method) %用指定的算法计算插值： ’nearest’：最近邻点插值，直接完成计算；
第十章 插值与拟合方法建模 在生产实际中，常常要处理由实验或测量所得到的一批离散数据，插值与拟合方法就是要通过这些数据去确定某一类已经函数的参数，或寻求某个近似函数使之与已知数据有较高的拟合精度。插值与拟合的方法很多，这里主要介绍线性插值方法、多项式插值方法和样条插值方法，以及最小二乘拟合方法在实际问题中的应用。相应的理论和算法是数值分析的内容，这里不作详细介绍，请参阅有关的书籍。 §1 数据插值方法及应用 在生产实践和科学研究中，常常有这样的问题：由实验或测量得到变量间的一批离散样点，要求由此建立变量之间的函数关系或得到样点之外的数据。与此有关的一类问题是当原始数据 ),(,),,(),,(1100n n y x y x y x 精度较高，要求确定一个初等函数)(x P y =（一般用多项式或分段 多项式函数）通过已知各数据点（节点），即n i x P y i i ,,1,0,)( ==，或要求得函数在另外一些点（插值点）处的数值，这便是插值问题。 1、分段线性插值 这是最通俗的一种方法，直观上就是将各数据点用折线连接起来。如果 b x x x a n =<<<= 10 那么分段线性插值公式为 n i x x x y x x x x y x x x x x P i i i i i i i i i i ,,2,1,,)(11 1 11 =≤<--+--= ----- 可以证明，当分点足够细时，分段线性插值是收敛的。其缺点是不能形成一条光滑曲线。 例1、已知欧洲一个国家的地图，为了算出它的国土面积，对地图作了如下测量：以由西向东方向为x 轴，由南向北方向为y 轴，选择方便的原点，并将从最西边界点到最东边界点在x 轴上的区间适当的分为若干段，在每个分点的y 方向测出南边界点和北边界点的y 坐标y1和y2，这样就得到下表的数据（单位：mm ）。
流量与管径、压力、流速的一般关系 一般工程上计算时，水管路，压力常见为0.1--0.6MPa，水在水管中流速在1--3米/秒，常取1.5米/秒。 流量=管截面积X流速=0.002827X管内径的平方X流速(立方米/小时)。 其中，管内径单位：mm ，流速单位：米/秒，饱和蒸汽的公式与水相同，只是流速一般取20--40米/秒。 水头损失计算Chezy 公式 这里： Q ——断面水流量（m3/s） C ——Chezy糙率系数（m1/2/s） A ——断面面积（m2） R ——水力半径（m） S ——水力坡度（m/m） 根据需要也可以变换为其它表示方法: Darcy-Weisbach公式
由于 这里： h f——沿程水头损失（mm3/s） f ——Darcy-Weisbach水头损失系数（无量纲） l ——管道长度（m） d ——管道内径（mm） v ——管道流速（m/s） g ——重力加速度（m/s2） 水力计算是输配水管道设计的核心，其实质就是在保证用户水量、水压安全的条件下，通过水力计算优化设计方案，选择合适的管材和确经济管径。输配水管道水力计算包含沿程水头损失和局部水头损失，而局部水头损失一般仅为沿程水头损失的5~10%，因此本文主要研究、探讨管道沿程水头损失的计算方法。 1.1 管道常用沿程水头损失计算公式及适用条件 管道沿程水头损失是水流摩阻做功消耗的能量，不同的水流流态，遵循不同的规律，计算方法也不一样。输配水管
道水流流态都处在紊流区，紊流区水流的阻力是水的粘滞力及水流速度与压强脉动的结果。紊流又根据阻力特征划分为水力光滑区、过渡区、粗糙区。管道沿程水头损失计算公式都有适用范围和条件，一般都以水流阻力特征区划分。 水流阻力特征区的判别方法，工程设计宜采用
压力和流速的关系如何计算 两管道之间的压差=a*l*p*u*u/2d 单位为pa a 为管道的摩擦系数，与管道的新旧和材质有关系。 l为你所取两点之间的距离单位为米 p为流体的密度 kg/m3 u为管内流体的流速，单位为米/秒 d为管子的管径，单位为米 请教：已知管道直径D，管道内压力P，能否求管道中流体的流速和流量？怎么求 已知管道直径D，管道内压力P，还不能求管道中流体的流速和流量。你设想管道末端有一阀门，并关闭的管内有压力P，可管内流量为零。管内流量不是由管内压力决定，而是由管内沿途压力下降坡度决定的。所以一定要说明管道的长度和管道两端的压力差是多少才能求管道的流速和流量。 对于有压管流，计算步骤如下： 1、计算管道的比阻S，如果是旧铸铁管或旧钢管，可用舍维列夫公式计算管道比阻s=0.001736/d^5.3 或用s=10.3n2/d^5.33计算，或查有关表格； 2、确定管道两端的作用水头差H=P/(ρg)，)，H 以m为单位；P为管道两端的压强差(不是某一断面的压强），P以Pa为单位； 3、计算流量Q：Q = (H/sL)^(1/2) 4、流速V=4Q/(3.1416d^2) 式中：Q―― 流量，以m^3/s为单位；H――管道起端与末端的水头差，以m^为单位；L――管道起端至末端的长度，以m为单位。
管道中流量与压力的关系 管道中流速、流量与压力的关系 流速:V=C√(RJ)=C√[PR/(ρgL)] 流量：Q=CA√(RJ)=√[P/(ρgSL)] 式中：C――管道的谢才系数；L――管道长度；P――管道两端的压力差；R――管道的水力半径；ρ――液体密度；g――重力加速度；S――管道的摩阻。 管道的内径和压力流量的关系 似呼题目表达的意思是：压力损失与管道内径、流量之间的关系，如果是这个问题，则正确的答案应该是：压力损失与流量的平方成正比，与内径5.33方成反比，即流量越大压力损失越大，管径越大压力损失越小，其定量关系可用下式表示： 压力损失（水头损失）公式（阻力平方区） h=10.3*n^2 * L* Q^2/d^5.33 上式严格说是水头损失公式，水头损失乘以流体重度后才是压力损失。式中n――管内壁粗糙度；L――管长；Q――流量；d――管内径
实验报告三 一、实验目的 通过本实验的学习，各种插值法的效果，如多项式插值法，牛顿插值法，样条插值法，最小二乘法拟合（即拟合插值），了解它们各自的优缺点及插值。 二、实验题目 1、 插值效果比较 实验题目：将区间[]5,5-10等份，对下列函数分别计算插值节点k x 的值，进行不同类型的插值，作出插值函数的图形并与)(x f y =的图形进行比较： 211)(x x f +=；x x f arctan )(=；4 2 1)(x x x f +=。 （1） 做拉格朗日插值； （2） 做三次样条插值。 2、 拟合多项式实验 实验题目：给定数据点如下表所示： i x -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 i y -4.45 -0.45 0.55 0.05 -0.44 0.54 4.55 分别对上述数据做三次多项式和五次多项式拟合，并求平方误差，作出离散函数()i i y x ,和拟合函数的图形。 三、实验原理 n 阶拉格朗日插值 设已知n x x x ,,,10 及()()()x L n i x f y n i i ,,,1,0 ==为不超过n 次的多项式，且满足 插值条件()().,,1,0n i y x L i i n ==由对()x L 2的构造经验，可设 ()()()()(),11000 n n n i i i n y x l y x l y x l y x l x L +++==∑= 其中，()()n i x L i ,,1,0 =均为n 次多项式且满足() .,,1,0,, ,0, ,1n j i j i j i x l j i =?? ?≠==不难验 证，这样构造出的()x L n 满足插值条件。因此问题归结为求()()n i x l i ,,1,0 =的表达式。因 ()i j x i ≠是n 次多项式()x l i 的n 个根，故可设
第七讲插值方法与数据拟合 § 7.1 引言 在工程和科学实验中，常常需要从一组实验观测数据(x i , y i ) (i= 1, 2, …, n) 揭示自变量x与因变量y 之间的关系，一般可以用一个近似的函数关系式y = f (x) 来表示。函数f (x) 的产生办法因观测数据与要求的不同而异，通常可采用两种方法：插值与数据拟合。 § 7.1.1 插值方法 1．引例1 已经测得在北纬32.3?海洋不同深度处的温度如下表： 根据这些数据，我们希望能合理地估计出其它深度（如500米、600米、1000米…）处的水温。 解决这个问题，可以通过构造一个与给定数据相适应的函数来解决，这是一个被称为插值的问题。 2．插值问题的基本提法 对于给定的函数表 其中f (x) 在区间[a, b] 上连续，x0，x1，…，x n为[a, b] 上n + 1个互不相同的点，要求在一个性质优良、便于计算的函数类{P(x)} 中，选出一个使 P(x i ) = y i，i= 0, 1, …, n(7.1.1) 成立的函数P(x) 作为 f (x) 的近似，这就是最基本的插值问题（见图7.1.1）。 为便于叙述，通常称区间[a, b] 为插值区间，称点x0，x1，…，x n为插值节点，称函数类{P(x)} 为插值函数类，称式(7.1.1) 为插值条件，称函数P(x) 为插值函数，称f (x) 为被插函数。求插值函数P(x) 的方法称为插值法。 § 7.1.2 数据拟合 1．引例2 在某化学反应中，已知生成物的浓度与时间有关。今测得一组数据如下： 根据这些数据，我们希望寻找一个y = f (t) 的近似表达式（如建立浓度y与时间t之间的经验公式等）。从几何上看，就是希望根据给定的一组点（1, 4.00），…,（16, 10.60），求函数y = f (t) 的图象的一条拟合曲
实验 2 插值与拟合 系班姓名学号 【实验目的】 1、掌握用MATLAB计算拉格朗日、分段线性、三次样条三种插值的方法，改变节点的数目， 对三种插值结果进行初步分析。 2、掌握用MATLAB作线性最小二乘的方法。 3、通过实例学习如何用插值方法与拟合方法解决实际问题，注意二者的联系和区别。 【实验内容】 预备：编制计算拉格朗日插值的M文件： 以下是拉格朗日插值的名为y_lagrl的M文件： function y=y_lagr1(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k)+s; end y(i)=s; end 第1题（d） 选择函数y=exp(-x2) (-2≤x≤2),在n个节点上（n不要太大，如5~11）用拉格朗日、分段线性、三次样条三种插值方法，计算m个插值点的函数值（m要适中，如50~100）。通过数值和图形输出，将三种插值结果与精确值进行比较。适当增加n，在作比较，由此作初步分析。 运行如下程序： n=7;m=61;x=-2:4/(m-1):2; y=exp(-x.^2); z=0*x; x0=-2:4/(n-1):2; y0=exp(-x0.^2); y1=y_lagr1(x0,y0,x); y2=interp1(x0,y0,x); y3=interp1(x0,y0,x,'spline');
[x'y'y1'y2'y3'] plot(x,z,'w',x,y,'r--',x,y1,'b:',x,y2,'m',x,y3,'b') gtext('y=exp(-x^2)'),gtext('Lagr.'),gtext('Piece.-linear.'),gtext ('Spline'), 将三种插值结果y1,y2,y3与精确值y 项比较，显然y1在节点处不光滑,拉格朗日插值出现较大的振荡,样条插值得结果是最好的.增加n 值（使n=11），再运行以上程序，得到的图形如右图所示，比较这两个图可发现，节点增加后，三种插值方法结果的准确度均有所提高，因此可近似地认为：增加节点个数可以提高插值结果的准确程度。 第3题 用给定的多项式，如y=x 3-6x 2+5x-3，产生一组数据（x i ,y i ,i=1,2,…,n ），再在yi 上添加随机干扰（可用rand 产生（0，1）均匀分布随机数，或用randn 产生N （0，1）分布随机数），然后用x i 和添加了随机干扰的y i 作3次多项式拟合，与原系数比较。如果作2或4次多项式拟合，结果如何？ 解：2 编制y_2_3.m 文件 n=15; x=0:8/(n-1):8; y=x.^3-6*x.^2+5*x-3; z=0*x; y0=y+rand(1,15); f=polyfit(x,y0,m); r=polyval(f,x) pl2ot(x,z,'k',x,y,'r:'r,'b') 程序及运行结果如下：m=2 ,y_2_3 f = 5.9888 -31.9916 17.6679 m=3 ,y_2_3
源压力4公斤，用25管径接出，不同管长流量是不等的呀！分别用1米、100米、1000米、100米长的管接出，流量悬殊很大！ 举例说明：200mm管道，水管起端与未端的压力差为0.8Mpa，当管道长度不同时，流量是不等的： 假设管道长度L=100m,水管起端与未端的压力差是0.8Mpa，即管道两端的水头差H=80m,则流量可计算如下： 200mm管道的摩阻S=9.029 流量 流速V=4Q/3.1416D^2=9.49 m/s 假设管道长度L=1000m,水管起端与未端的压力差是0.8Mpa，即管道两端的水头差H=80m,则流量可计算如下： 200mm管道的摩阻S=9.029 流量 流速V=4Q/3.1416D^2= 2.99m/s 由以上计算知，管道长度不同，其他条件一样，流量是不等的。水头差与管道长度的比值（H/L）称为水力坡度，有水头（或有压力）不一定有流量，流量大小取决与水力坡度（或压力坡度），而不是取决于水头（或压力）！ 安的梦|六级 DN 15、DN 25、DN50管径的截面积分别为： DN15：152平方毫米，合0.0177平方分米。
DN25：252平方毫米，合0.0491平方分米。 DN50：502平方毫米，合0.1963平方分米。 设管道流速为V=4米/秒，即V=40分米/秒，且1升=1立方分米，则管道的流量分别为（截面积乘以流速）： DN15管道： 流量Q=0.0177*40=0.708升/秒， 合2.55立方米/小时。 DN25管道： 流量Q=0.0491*40=1.964升/秒， 合7.07立方米/小时。DN50管道： 流量Q=0.1963*40=7.852升/秒， 合28.27立方米/小时。 注： 必须给定流速才能计算流量，上述是按照4米/秒计算的。 是600立方米/小时吧！ 流量Q=600立方米/小时= 0.167m^3/s 流速V=4Q/(3.1416D^2)= 4*0.167/(3.1416*0.090^2)= 26.25 m/s 流速非常大！可能是600立方米/日，则： 流量Q=600立方米/日= 0.00694 m^3/s 流速V=4Q/(3.1416D^2)= 4*0.00694/(3.1416*0.090^2)= 1.09 m/s 流速属于正常。
1、流速计算： 按照伯努利方程，假设条件为水平管，管口为大气压。 则p1+ρ1gz1+(1/2)*ρ1v1^2=p2+ρ2gz2+(1/2)*ρ2v2^2 由于ρ1gz1=ρ2gz2；v1=0；p2=0.1MPa；ρ2为水的密度=1000kg/m3； p1=1.1MPa(管道内的绝对压力）； 公式化简为：p1=p2+(1/2)*ρ2v2^2 按照已知条件计算得出v2=44.72m/s 这是管道敞口端的计算流速，实际中不会有这么高，因为管道敞口端压力不一定是大气压。 2、流量计算： Q=ρ.s.v2=1000*3.14/4*0.2*0.2*44.72=1404 kg/s 每小时的出水量=1404*3600/1000=5054（吨） 这个计算值明显偏大，但是计算结果是这样，我无奈。 根据我实际中见到的自来水管道的水量估算，压力为4公斤，管径为DN40，每小时最大的流量大概16吨。按照这个比例折下来你的管子每小时流量大概为1000吨。 DN15、DN25、DN50管径的截面积分别为： DN15：152*3.14/4=176.625平方毫米，合0.0177平方分米。 DN25：252*3.14/4=490.625平方毫米，合0.0491平方分米。 DN50：502*3.14/4=1962.5平方毫米，合0.1963平方分米。 设管道流速为V=4米/秒，即V=40分米/秒，且1升=1立方分米，则管道的流量分别为（截面积乘以流速）： DN15管道：流量Q=0.0177*40=0.708升/秒， 合2.55立方米/小时。 DN25管道：流量Q=0.0491*40=1.964升/秒， 合7.07立方米/小时。 DN50管道：流量Q=0.1963*40=7.852升/秒， 合28.27立方米/小时。 注：必须给定流速才能计算流量，上述是按照4米/秒计算的。 任何气体流量的计算都可以用密度乘速度乘面积来计算，你给的条件中面积已经知道了，密度可以通过压力和温度来计算(用理想气体公式或者查表)，速度虽然计算不出来，但是可以用两个公式解方程得到。根据流量连续(即按阀前截面计算流量应等于按阀后截面计算的流量)，设两个未知数，阀前速度和阀后速度，得到一个方程。在根据能量守恒原理，阀前压力势能+阀前动能= 阀后的压力势能+阀后的动能。两个方程，两个未知数，可以解出速度，这样就可以计算出流量。这只是理论上的方法，可能和实际的流量有微小的差别 题目分析：流量为1小时10吨，这是质量流量，应先计算出体积流量，再由体积流量计算出管径，再根据管径的大小选用合适的管材，并确定管子规格。
计算方法实验 题目： 班级： 学号： 姓名：
目录 计算方法实验 (1) 1 实验目的 (3) 2 实验步骤 (3) 2.1环境配置： (3) 2.2添加头文件 (3) 2.3主要模块 (3) 3 代码 (4) 3.1主程序部分 (4) 3.2多项式方程部分 (4) 3.3核心算法部分 (8) 3.4数据结构部分 (13) 4运行结果 (19) 4.1拉格朗日插值法运行结果 (19) 4.2牛顿插值法运行结果 (20) 4.3多项式拟合运行结果 (20) 5总结 (21) 拉格朗日插值法 (21) 牛顿插值法 (21) 多项式拟合 (21) 6参考资料 (22)
1 实验目的 1.通过编程对拉格朗日插值法、牛顿插值法以及多项式拟合数据的理解 2.观察上述方法的计算稳定性和求解精度并比较各种方法利弊 2 实验步骤 2.1环境配置： VS2013，C++控制台程序 2.2添加头文件 #include "stdio.h" #include "stdlib.h" #include "stdafx.h" 2.3主要模块 程序一共分成三层，最底层是数据结构部分，负责存储数据，第二层是交互部分，即多项式方程部分，负责输入输出获得数据，最上层是核心的算法部分，负责处理已获得的数据。具体功能如下： ●数据结构部分 数据结构部分是整个程序的最底层，负责存储部分。因方程系数作为数据元素插入和删除操作较少，而顺序表空间利用率大且查看方便，故此程序选用顺序表保存系数。数据结构文件中写的是有关顺序表的所有基本操作以供其他文件调用。本次实验使用列主元高斯消元法作为求解方程组的方法，所以也用了二维顺序表存储数组。综上，数据结构部分文件是前两个试验的文件内容和，稍作修改。 ●常系数微分方程部分 多项式方程部分是程序的第二层，内容主要是常系数微分方程导数的计算和显示菜单部分。 ●算法部分 算法部分分为两个文件，一个是插值部分，一个是拟合部分。 插值部分文件负责有关插值的核心算法，处于整个程序最上层部分，负责拉格朗日插值法和牛顿插值法的具体实现过程。调用方程文件的函数，将获得的数据进行处理运算，将结果返回给方程主函数和输出的第二层。每种方法有两个函数，一个为仅仅实现一次插值的算法，另一个是和方程部分联系的
Matlab中插值拟合函数汇总和使用说明 命令1 interp1 功能一维数据插值（表格查找）。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。 x：原始数据点 Y：原始数据点 xi：插值点 Yi：插值点 格式 (1)yi = interp1(x,Y,xi) 返回插值向量yi，每一元素对应于参量xi，同时由向量x 与Y 的内插值决定。参量x 指定数据Y 的点。 若Y 为一矩阵，则按Y 的每列计算。yi 是阶数为 length(xi)*size(Y,2)的输出矩阵。 (2)yi = interp1(Y,xi) 假定x=1:N，其中N 为向量Y 的长度，或者为矩阵Y 的行数。(3)yi = interp1(x,Y,xi,method) 用指定的算法计算插值： ’nearest’：最近邻点插值，直接完成计算； ’linear’：线性插值（缺省方式），直接完成计算； ’spline’：三次样条函数插值。对于该方法，命令interp1 调用
函数spline、ppval、mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函数。命令spline 用它们执行三次样条函数插值； ’pchip’：分段三次Hermite 插值。对于该方法，命令interp1 调用函数pchip，用于对向量x 与y 执行分段三次内插值。该方法保留单调性与数据的外形； ’cubic’：与’pchip’操作相同； ’v5cubic’：在MATLAB 5.0 中的三次插值。 对于超出x 范围的xi 的分量，使用方 法’nearest’、’linear’、’v5cubic’的插值算法，相应地将返回NaN。对其他的方法，interp1 将对超出的分量执行外插值算法。 (4)yi = interp1(x,Y,xi,method,'extrap') 对于超出x 范围的xi 中的分量将执行特殊的外插值法extrap。(5)yi = interp1(x,Y,xi,method,extrapval) 确定超出x 范围的xi 中的分量的外插值extrapval，其值通常取NaN 或0。 例1 1. 2.>>x = 0:10; y = x.*sin(x); 3.>>xx = 0:.25:10; yy = interp1(x,y,xx); 4.>>plot(x,y,'kd',xx,yy) 例2 1.
曲线拟合的数值计算方法实验 郑发进 2012042020022 【摘要】实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合（curve fitting）是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。曲线直线化是曲线拟合的重要手段之一。对于某些非线性的资料可以通过简单的变量变换使之直线化，这样就可以按最小二乘法原理求出变换后变量的直线方程，在实际工作中常利用此直线方程绘制资料的标准工作曲线，同时根据需要可将此直线方程还原为曲线方程，实现对资料的曲线拟合。常用的曲线拟合有最小二乘法拟合、幂函数拟合、对数函数拟合、线性插值、三次样条插值、端点约束。 关键词曲线拟合、最小二乘法拟合、幂函数拟合、对数函数拟合、线性插值、三次样条插值、端点约束 一、实验目的 1.掌握曲线拟合方式及其常用函数指数函数、幂函数、对数函数的拟合。 2.掌握最小二乘法、线性插值、三次样条插值、端点约束等。 3.掌握实现曲线拟合的编程技巧。 二、实验原理 1.曲线拟合 曲线拟合是平面上离散点组所表示的坐标之间的函数关系的一种数据处理方法。用解析表达式逼近离散数据的一种方法。在科学实验或社会活动中，通过实验或观测得到量x与y的一组数据对（X i，Y i)（i=1，2，...m），其中各X i 是彼此不同的。人们希望用一类与数据的背景材料规律相适应的解析表达式，y=f(x，c）来反映量x与y之间的依赖关系，即在一定意义下“最佳”地逼近或拟合已知数据。f(x，c)常称作拟合模型，式中c=(c1，c2，…c n)是一些待定参数。