第二章 - 城市供水量的预测模型 - 插值与拟合算法

第二章 城市供水量的预测模型

——插值与拟合算法

2.1 城市供水量的预测问题

2.1.1 实际问题与背景

为了节约能源和水源，某供水公司需要根据日供水量记录估计未来一时间段（未来一天或一周）的用水量，以便安排未来（该时间段）的生产调度计划。现有某城市7年用水量的历史记录，记录中给出了日期、每日用水量（吨/日）。如何充分地利用这些数据建立数学模型，预测2007年1月份城市的用水量，以制定相应的供水计划和生产调度计划。 表2.1.1 某城市7年日常用水量历史记录（万吨/日）

日期 日用水量 20000101 122.1790 20000102 128.2410 …… …… 20061230 150.40168 20061231 148.2064 表2.1.2 2000-2006年1月城市的总用水量（万吨/日）

年份 2000 2001 2002 2003 2004 2005 2006 用水量 4032．41 4186．0254 4296．9866 4374．852 4435．2344 4505．4274 4517．6993 利用这些数据，可以采用时间序列、灰色预测等方法建立数学模型来预测2007年1月份该城市的用水量。如果能建立该城市的日用水量随时间变化的函数关系，则用该函数来进行预测非常方便。但是这一函数关系的解析表达式是没办法求出来的，那么能否根据历史数据求出该函数的近似函数呢？根据未知函数的已有数据信息求出其近似函数的常用方法有插值法和数据拟合。本章将介绍插值法和数据拟合，并用这两种方法对该城市的供水量进行预测。

2.2 求未知函数近似表达式的插值法

2.2.1 求函数近似表达式的必要性

一般地，在某个实际问题中，虽然可以断定所考虑的函数f(x)在区间[a,b]上存在且连续，但却难以找到它的解析表达式，只能通过实验和观测得到该函数在有限个点上的函数值（即一张函数表）。显然，要利用这张函数表来分析函数的性态，甚至直接求出其它一些点上的函数值是非常困难的。在有些情况下，虽然可以给出函数f(x)的解析表达式，但由于结构相当复杂，使用起来很不方便。面对这些情况，希望根据所得函数表（或结构复杂的解析表达式），构造某个简单函数?(x)作为未知函数f(x)的近似。插值法是解决此类问题的一种常用的经典方法，它不仅直接广泛地应用于生产实际和科学研究中，而且也是进一步学习数值计算的基础。

17

定义2.2.1 设函数y?f(x)在区间[a,b]上连续，且在n?1个不同的点

a?x0,x1,,xn?b上分别取值y0,y1,在一个性质优良、便于计算的函数类?,yn，

中，求一简单函数?(x)，使

??xi??yi?i?0,1,n? (2.2.1)

而在其它点x?xi上作为f(x)的近似。称区间[a,b]为插值区间，点x0,x1,,x为

插值节点，称（2.2.1）为f(x)的插值条件，称函数类?为插值函数类，称?(x)为函数在节点x0,x1,,xn处的插值函数。求插值函数?(x)的方法称为插值法。

插值函数类?的取法不同，所求得的插值函数?(x)逼近f(x)的效果就不同，它的选择取决于使用上的需要。常用的有代数多项式、三角多项式和有理函数等。

当选用代数多项式作为插值函数时，相应的插值问题就称为多项式插值。这里主要介绍多项式插值。

在多项式插值中，求一次数不超过n的代数多项式

Pn?x??a0?a1x?使

?anxn (2.2.2) （2.2.3） ,n?

P,,n?xi??yi?i?0,1其中a0,a1,满足插值条件(2.2.3)的多项式(2.2.2)，称为函数f(x),an为实数。

的n次插值多项式。

n?1个点n次插值多项式Pn?x?的几何意义：过曲线y?f(x)上的

如图(xi,yi)i(?0,1,n作一条,n次代数曲线y?Pn(x)作为曲线y?f(x)的近似，2.2.1所示。

Yy?f?x?y?Pn?x?y0y1x1yn0ax0xnbX

图2.2.1 插值多项式的几何直观图

18

2.2.2 插值多项式的存在唯一性

由插值条件（2.2.3）知，P,n?x?的系数ai?i?0,1n?满足线性方程组：

n?a0?a1x0?anx0?y0?n?a0?a1x1??anx1?y1 （2.2.4） ??n??a0?a1xn??anxn?yn由线性代数知，线性方程组的系数行列式是n?1阶范德蒙（Vandermonde）行列式，且

1x0x1xn2x0nx0 V?11x122xnx1nnxn????xi?xj?

i?1j?0ni?1因x0,x1,xn是区间[a,b]上的不同点，上式右端乘积中的每一个因子xi?xj?0，

于是系数行列式不等于0，即方程组（2.2.4）的解存在且唯一。从而得出插值多项式的存在唯一性定理。 定理2.2.1 若插值节点x0,x1,xn互不相同，则满足插值条件（2.2.3）的n次插值多项式（2.2.2）存在且唯一。

2.3 求插值多项式的Lagrange（拉格朗日）法

在上一节里，插值多项式存在唯一性的证明过程不仅指出了满足差值条件的n次插值多项式是存在唯一性的，而且也提供了插值多项式的一种求法，即通过解线性方程组（2.2.4）来确定其系数ai。但是，当未知数个数多时，这种做法的计算工作量大，不便于实际应用。Lagrange基于用简单插值问题的插值函数表示一般的插值函数思想，给出一种求插值函数的简便方法——Lagrange插值法。

2.3.1 Lagrange插值基函数

先考虑简单的插值问题：对节点xi(i?0,1,n)中任意一点xk?0?k?n?做

一n次多项式lk(x)使它在该点上取值为1,而在其余点xi?i?0,1,k?1,k?1,上取值为零, 即满足插值条件

,n??1i?k （2.3.1） lk(xi)??0i?k?（2.3.1）表明n个点xi?i?0,1,k?1,k?1,n?都是n次多项式lk(x)的零点，故可

19

设lk(x)?Ak(x?x0)(x?x1)条件lk(x)?1可得

Ak?(x?xk?1)(x?xk?1)(x?xn)，其中Ak为待定系数，由

(xk?x0)1(xk?xk?1)(xk?xk?1)(xk?xn)

故

lk(x)?(x?x0)(x?xk?1)(x?xk?1)(x?xn) （2.3.2）

(xk?x0)(xk?xk?1)(xk?xk?1)(xk?xn)对应于每一节点xk?0?k?n?都能求出一个满足插值条件（2.3.1）的n次插值多项式（2.3.2），这样，由（2.3.2）式可以求出n?1个n次插值多项式

l0(x),l1(x),容易看出，这组多项式仅与节点的取法有关，称它们为在n?1,ln(x)。

个节点上的n次基本插值多项式或n次插值基函数，即Lagrange插值基函数。

2.3.2 Lagrange（拉格朗日）插值多项式

利用Lagrange插值基函数立即可以写出满足插值条件（2.2.3）的n次插值多项式

y0l0(x)?y1l1(x)? ?ynln(x) （2.3.3）

事实上，由于每个插值基函数lk(x)(k?0,1,,n)都是n次多项式，故其线性

组合（2.3.3）必是不高于n次的多项式，同时，根据条件（2.2.1）容易验证多项式（2.3.3）在节点xi处的值yi(i?0,1,,n)为Pn(xi)，因此，它就是待求的n次

插值多项式Pn(x)。形如(2.3.3)的插值多项式称为Lagrange插值多项式，记为

Ln?x??y0l0(x)?y1l1(x)? ??ykk?0n?ynln(x)

(x?x0)(x?xk?1)(x?xk?1)(x?xn) （2.3.4）

(xk?x0)(xk?xk?1)(xk?xk?1)(xk?xn)x?x0x?x1?y1 （2.3.5） x0?x1x1?x0令n?1，由（2.3.4）即得两点插值公式

L1(x)?y0即

L1(x)?y0?y1?y0(x?x0) （2.3.6） x1?x0这是一个线性函数，用线性函数L1(x)近似代替函数f(x)，在几何上就是通过曲

20

线y?f(x)上两点(x0,y0),(x1,y1)做一直线y?L1(x)近似代替曲线y?f(x)，如见图2.3.1所示，故两点插值又称线性插值。

y?f?x?y(xy0,0)(x,y)11y?L1?x?0x0x1x

图2.3.1 线性插值的几何直观图

令n?2，由（2.3.4）可得常用的三点插值公式：

L2?x??y0(x?x0)(x?x2)(x?x0)(x?x1)(x?x1)(x?x2) （2.3.7） ?y1?y2(x0?x1)(x0?x2)(x1?x0)(x1?x2)(x2?x0)(x2?x1)这是一个二次函数，用二次函数L2(x)近似代替函数f(x)，在几何上就是通过曲线y?f(x)上的三点(x0,y0),(x1,y1),(x2,y2)作一抛物线y?L2(x)，近似地代替曲线y?f(x)，如图2.3.1所示，故称为三点插值或二次插值，也称为抛物线插值。

y(x,y)00y?L2?x?(x,y)11(x,y)22y?f?x?0x0x1x

图2.3.2二次插值的几何直观图

例2.3.1 已知100?10,121分别用线性插值和抛物线插值求?11,144?12115的值。

解 因为115在100和121之间，故取节点x0?100，x1?121，相应地有y0?10，

y1?11，于是，由线性插值公式（2.3.5）可得

L1(x)?10?x?121x?100?11?

100?121121?100故用线性插值求得的近似值为：

115?121115?100?11??10.714

100?121121?100同理，用抛物线插值公式（2.3.7）所求得的近似值为：

115?L1(115)?10?

21

115?L2?115??10??115?121??115?144??11??115?100??115?144??100?121??144?121??121?100??121?144?115?100??115?121???12? ?144?100??144?121??10.732将所得结果与115的精确值10.7328相比较，可以看出抛物线插值的精确度较

好。为了便于上机计算，我们常将拉格朗日插值多项式（2.3.4）改写成公式

（2.3.8）的对称形式

?n?x?xj?Ln(x)??yk???j?0?xk?xjk?0?j?k??n??? （2.3.8） ??????可用二重循环来完成Ln(x)值的计算，先通过内循环，即先固定k，令j从0到

?x?xjn(j?k)累乘求得 lk(x)????j?0?xk?xjj?kn?，然后再通过外循环，即令k从0到n，累???加得出插值结果Ln(x)。

2.3.3 插值余项

在插值区间[a,b]上用插值多项式Pn(x)近似代替f(x)，除了在插值节点xi上没有误差以外，在其他点上一般有误差。若记

Rn(x)?f(x)?Pn(x)

则Rn(x)就是用Pn(x)近似代替f(x)时所产生的截断误差，称Rn(x)为插值多项式

Pn(x)的余项。关于插值多项式的误差有定理2.3.1的估计式。 定理2.3.1 设f(x)在区间[a,b]上有直到n?1阶导数,x0,x1,n?1个互异的节点，Pn(x)为满足条件:P,n(xi)?f(xi)(i?0,1,xn为区间[a,b]上,n) 的n次插值多

项式，则对于任何x??a,b?有

f(n?1)(?)（2.3.9） Rn(x)??n?1(x)

(n?1)!其中?n?1(x)??(x?xi) ,??(a,b)且依赖于x。

i?0n 22

证明 由插值条件Pn(xi)?f(xi)知Rn(xi)?0(i?0,1,的零点,故可设

,n)，即插值节点都是Rn(x)Rn(x)?K(x)?n?1(x) (2.3.10)

其中K(x)为待定函数。下面求K(x)，对区间[a,b]上异于xi的任意一点x?xi作辅助函数：

F(t)?f(t)?Pn(t)?K(x)?n?1(t)

不难看出F(t)具有如下特点 （1）F(x)?F(xi)?0(i?0,1,,n) (2.3.11)

（2）在[a,b]上有直到n?1阶导数，且

（2.3.12） F(n?1)(t)?f(n?1)(t)?K(x)(n?1)!

等式（2.3.11）表明F(t)在[a,b]上至少有n?2个互异的零点，根据罗尔(Rolle)定理，在F(t)的两个零点之间，F'(t)至少有一个零点，因此，F'(t)在(a,b)内至少有n?1个互异的零点，对F'(t)再应用罗尔定理，推得F''(t)在(a,b)内至少有n个互异的零点。继续上述讨论，可推得F(n?1)(t)在(a,b)内至少有一个零点，若记为?，则F(n?1)(?)?0，于是由（2.3.12）式得

f(n?1)(?)?K(x)(n?1)!?0

f(n?1)(?) K(x)?(n?1)!将它代入(2.3.10)即得（2.3.9）。对于x?xi，（2.3.9）显然成立。

例2.3.2 在例2.3.1中分别用线性插值和抛物插值计算了115近似值，试估计它们的截断误差。

解 用线性插值求f(x)?x的近似值，其截断误差由插值余项公式（2.3.9）知

R1(x)?1f''(?)?2(x)21????3/2(x?x0)(x?x1)???x0,x1?8

23

现在x0?100，x1?121，x?115，故

1R1(115)??(115?100)(115?121)max??3/2??[100,121]8

1??15?6?10?3?0.011258当用抛物插值求f(x)?x的近似值时，其截断误差为

R2(x)???521f'''(?)?3(x)3!1?(x?x0)(x?x1)(x?x2),???x0,x2?16

将x0?100,x1?121,x2?144,x?115代入，即得

R2(115)?1(115?100)(115?121)(115?144)?10?5?0.0017 162.3.4 插值误差的事后估计法

在许多情况下，要直接应用余项公式（2.3.9）来估计误差是很困难的，下面将以线性插值为例，介绍另一种估计误差的方法。

设x0?x?x1?x2且f(xi)(i?0,1,2)为已知，若将用x0,x1两点做线性插值求得y?f(x)的近似值为y1，用x0,x2两点作线性插值所求得y?f(x)的近似值记为

y2，则由余项公式（2.3.9）知：

1f''(?1)(x?x0)(x?x1),?1?[x0,x1]2

1y?y2?f''(?2)(x?x0)(x?x2),?2?[x0,x2]2y?y1?假设f\x)在区间?x0,x2?中变化不大，将上面两式相除，即得近似式

y?y1x?x1? y?y2x?x2即

y?y1?x?x1(y2?y1) （2.3.13） x2?x1近似式（2.3.13）表明，可以通过两个结果的偏差y2?y1来估计插值误差，这种

24

直接利用计算结果来估计误差的方法，称为事后估计法。

在例2.3.1中，用x0?100,x1?121做节点，算得的115近似值为y1?10.714，同样，用x0?100,x2?144做节点，可算得115的另一近似值y2?10.682。

通过（2.3.13）可以估计出插值y1的误差为：

115?y1?115?121(10.682?10.714)?0.00835

144?1212.4 求插值多项式的Newton法

拉格朗日插值法通过求出拉格朗日基函数而方便地表示插值多项式，但是当增加节点时，所有的基函数都需重新计算，因此其计算过程不具继承性。由线性代数可知，任何一个不高于n次的多项式，都可表示成函数

1,x?x0,(x?x0)(x?x1),,(x?x0)(x?x1)(x?xn?1)

的线性组合，即可将满足插值条件P(xi)?yi(i?0,1,,n) 的n次多项式写成形式

(x?xn?1)

a0?a1(x?x0)?a2(x?x0)(x?x1)?其中ak?k?0,1,?an(x?x0)(x?x1)这种形式的插值多项式称为牛顿﹙Newton﹚插,n?为待定系数。

值多项式，把它记成Nn(x),即

Nn(x)?a0?a1(x?x0)?a2(x?xo)(x?x1)??an(x?x0)(x?x1) (x?xn?1) （2.4.1）

因此，牛顿插值多项式Nn(x)是插值多项式Pn(x)的另一种表示形式。在牛顿插值多项式中要用到的差分与差商等概念，与数值计算的其它方面有着密切的关系。

2.4.1 向前差分与Newton向前插值公式

设函数f(x)在等距节点xk?x0?kh(k?0,1,,n)处的函数值f(xk)?yk为已

知，其中h是正常数，称为步长，称两个相邻点xk和xk?1处函数值之差 yk?1?yk为函数f(x)在点xk处以h为步长的一阶向前差分﹙简称一阶差分﹚，记?yk，即

?yk?yk?1?yk

于是，函数f(x)在各节点处的一阶差分依次为

?y0?y1?y0,?y1?y2?y1,,?yn?1?yn?yn?1

又称一阶差分的差分?2yk????yk???yk?1??yk为二阶差分。

25

一般地，定义函数f(x)在点xk处的m阶差分为：

?myk??m?1yk?1??m?1yk

为了便于计算与应用，通常采用表格形式计算差分，如表2.4.1所示。

表2.4.1 差分计算表

xkyyy1k?y?y?y0k?y2k?y3k?y4kxxxxx001?y20?y3102yy2?y33?y2?y4012?y31?y224y?y34

在等距节点xk?x0?kh(k?0,1,,n)情况下，可以利用差分表示牛顿插值多项式

﹙2.4.1﹚ 的系数，并将所得公式加以简化。

事实上，由插值条件Nn(x0)?y0即可得a0?y0。 再由插值条件Nn(x1)?y1可得：

a1?y1?y0?y0 ?x1?x0h由插值条件Nn(x2)?y2可得：

y2?y0??y0(x2?x0)y2?2y1?y0?2y0h??

(x2?x0)(x2?x1)2hh2!?h2a2?一般地，由插值条件Nn(xk)?yk可得

?kyoak?(k?1,2,kk!?h,n)

于是，满足插值条件Nn(xi)?yi的插值多项式为：

26

?y0?2y0Nn(x)?y0?(x?x0)?(x?x0)(x?x1)?h2!?h2?ny0?(x?x0)(x?x1)n!?hn(x?xn?1)

令x?x0?th(t?0)，并注意到xk?x0?kh则牛顿插值多项式可简化为

t(t?1)2t(t?1)(t?n?1)n?y0???y0 （2.4.2） 2!n!这个用向前差分表示的插值多项式，称为牛顿向前插值公式，简称前插公式。它

Nn(x0?th)?y0?t?y0?适用于计算x0附近的函数值。

由插值余项公式（2.3.9），可得前插公式的余项为：

Rn(x0?th)?t(t?1)(t?n)n?1n?1hf(?),??(x0,xn) （2.4.3）

(n?1)!例2.4.1 从给定的正弦函数表出发计算sin(0.12),并估计截断误差。

表2.4.2 sinx的函数及差分表 x 0.1 0.2 0.3 0.4 0.5 0.6 sinx 0.09983 0.19867 0.29552 0.38942 0.47943 0.56464 ?y 0.09884 0.09685 0.09390 0.09001 0.08521 ?2y -0.00199 -0.00295 -0.00389 -0.00480 ?3y -0.00094 -0.00094 -0.00091 解 因为0.12介于0.1与0.2之间，故取x0?0.1，此时

t?x?x00.12?0.1??0.2， h0.1为求?y0,?2y0,?3y0，构造差分表2.4.2。

若用线性插值求sin(0.12)的近似值，则由前插公式（2.4.2）立即可得

sin(0.12)?N1(0.12)?0.09983?0.2?0.09984?0.11960

用二次插值得：

sin(0.12)?N2(0.12)?0.09983?0.2?0.09884??N1(0.12)?0.00016?0.119760.2?(0.2?1)?(?0.00199) 2用三次插值得：

27

sin(0.12)?N3(0.12)?N2(0.12)??0.119710.2?(0.2?1)?(0.2?2)?(?0.00096)

6因N3(0.12)与N2(0.12)很接近，且由差分表2.4.2可以看出，三阶差分接近于常数（即?4y0接近于零），故取N3(0.12)?0.11971作为sin(0.12)的近似值，此时由余项公式（2.4.3）可知其截断误差为：

R3(0.12)?0.2?(0.2?1)?(0.2?2)?(0.2?3)?(0.1)4?sin(0.4)?0.000002

242.4.2 向后差分与Newton（牛顿）向后插值公式

在等距节点xk?x0?kh(k?0,1,,n)下，除了向前差分外，还可引入向后差

分和中心差分，其定义和记号分别如下：

y?f(x)在点xk处以h为步长的一阶向后差分和m阶向后差分分别为：

?yk?yk?yk?1?yk??mm?1yk??m?1yk?1(m?2,3,)

y?f(x)在xk点处以h为步长的一阶中心差分和m阶中心差分分别为：

?yk?yk?12?yk?12?myk??m?1y其中yk?12??m?1yk?12(m?2,3,)

k?12h?h????f?xk??,y1?f?xk??，各阶向后差分与中心差分的计算，可通

2?k?22???过构造向后差分表与中心差分表来完成。

利用向后差分，可简化牛顿插值多项式（2.4.1），导出与牛顿前插公式（2.4.2）类似的公式。

事实上，若将节点的排列次序看作xn,xn?1x0，那么?2.4.1)可写成：

Nn(x)?b0?b1(x?xn)?b2(x?xn)(x?xn?1)?根据插值条件Nn(xi)?yi(i?n,n?1,Nn(xn?th)?yn?t?yn??bn(x?xn)(x?xn?1)(x?x1)

,1,0)得到用向后差分表示的插值多项式：

t(t?1)2t(t?1)(t?n?1)n?yn???yn （2.4.4） 2!n!其中t<0，插值多项式(2.4.4)称为牛顿向后插值公式，简称后插公式。它适

28

用于计算xn附近的函数值。由插值余项公式（2.3.9），可写出后插公式的余项为：

Rn(xn?th)?t(t?1)(t?n)n?1(n?1)hf(?),(??(x0,xn)) （2.4.5）

(n?1)!例2.4.2 已知函数表如表2.4.2，计算sin(0.58)，并估算其截断误差。 解 因为0.58位于表尾x5?0.6附近，故用后插公式（2.4.4）计算sin(0.58)的近似值。

一般的，为了计算函数在x5处的各阶向后差分，应构造向后差分表。但由向前差分与向后差分的定义可以看出，对同一函数表来说，构造出来的向后差分表与向前差分表在数据上完全相同。因此，表2.4.2中用下划线标出的各数依次给出了sinx在x5?0.6处的函数值和向后差分值。因三阶向后差分接近于常数，故

x?x50.58?0.6???0.2， h0.1于是由牛顿向后插值公式（2.4.4）得：

用三次插值进行计算，且t?sin(0.58)?N3(0.58)?0.56464?(?0.2)?0.08521?(?0.2)?(?0.2?1)?(?0.00480)2(?0.2)?(?0.2?1)?(?0.2?2) ??(?0.00091)6?0.54802因为在整个计算中，只用到四个点x?0.6,0.5,0.4,0.3上的函数值，固由余项公式（2.4.5）知其截断误差为：

R3(0.58)??0.2?(?0.2?1)?(?0.2?2)?(?0.2?3)?(0.1)4?sin(0.6)?0.000002

242.4.3 差商与牛顿基本插值多项式

当插值节点非等距分布时，就不能用差分来简化牛顿插值多项式，此时可用差商这个新概念来解决。

设函数f(x)在一串互异的点xi0,xi1,xi2,上的值依次为

f(xi0)、f(xi1)、f(xi2)称函数值之差 f(xi1)?f(xi0)与自变量之差xi1?xi0的比值

f(xi1)?f(xi0)xi1?xi0为

函数f(x)关于 xi1,xi0点的一阶差商，记作f[xi0,xi1]

29

例2.4.3 f[x0,x1]?f(x1)?f(x0)f(x2)?f(x1),f[x1,x2]?,x1?x0x2?x1。

称一阶差商的差

f[xi1,xi2]?f[xi0,xi1]xi2?xi0为函数f(x)关于点xi0、xi1、xi2的二

阶差商，记作

f[xi0,xi1,xi2]。

f[x1,x2]?f[x0,x1]。

x2?x0例2.4.4 f[x0,x1,x2]?一般地，可通过函数f(x)的m?1阶差商定义的m阶差商如下：

f[xi0,xi1,xim]?f[xi1,,xim]?f[xi0,xim?xi0,xim?1]

差商计算也可采用表格形式（称为差商表）进行，如表2.4.3所示。

表2.4.3 差商计算表

xk f(xk) 一阶差商 二阶差商 三阶差商 x0 f(x0) x1 f(x1) x2 f(x2) f[x0,x1] f[x1,x2] f[x2,x3] f[x0,x1,x2] f[x1,x2,x3] f[x0,x1,x2,x3] x3 f(x3) 差商具有下列重要性质（证明略）： （1）函数f(x)的m阶差商f[x0,x1,线性组合表示，且

,xm]可由函数值f?x0?、f?x1?、、f?xm?的

f[x0,x1,,xm]??i?0m(xi?x0)f?xi?(xi?xi?1)(xi?xi?1)(xi?xm).

（2）差商具有对称性，即任意调换节点的次序，不影响差商的值。例如

f[x0,x1,x2]?f[x1,x0,x2]?f[x1,x2,x0]?（3）当f?m??x?在包含节点xij?j?0,1,在xi0,xi1,,m?的某个区间上存在时，

xim

30

之间必有一点?，使

f[xi，xi1，,xim]?0f?m?m!???

（4）在等距节点 xk?x0?kh?k?0,1,与差商，且有下面关系：

,n?情况下，可同时引入m?m?n?阶差分

?my0f[x0,x1,xm]?m!?hm m?ynf[xn,xn?1,,xn?m]?m!?hm引入差商的概念后，可利用差商表示牛顿插值多项式（2.4.1）的系数。事实上，从插值条件出发，可以像确定前插公式中的系数那样，逐步地确定（2.4.1）中的系数

??a0?f(x0)???ak?f[x0,x1,,xk](k?1,2,,n)

故满足插值条件Nn?xi??yi?i?0,1,,n?的n次插值多项式为：

（2.4.6）

Nn?x??f?x0??f?x0,x1??x?x0??f?x0,x1,x2??x?x0??x?x1???f[x0,x1,,xn](x?x0)(x?x1)(x?xn?1)（2.4.6）称为牛顿基本插值多项式，常用来计算非等距节点上的函数值。 例2.4.5 试用牛顿基本差值多项式按例1要求重新计算115的近似值。 解 先构造差商表

x 100 121 144 x 一阶差商 二阶差商 10 11 12 0.043478 0.047619 -0.000094 由上表可以看出牛顿基本插值多项式（2.4.6）中各系数依次为：

f(x0)?10 f[x0,x1]?0.047619 f[x0,x1,x2]??0.000094

故用线性插值所得的近似值为：

31

115?N1(115)?10?0.047619?(115?100)?10.7143

用抛物插值所求得的近似值为：

115?N2(115)?N1(115)?(?0.000094)?(115?100)?(115?121)?10.7228 所得结果与例2.4.1一致。比较例2.4.1和例2.4.5的计算过程可以看出，与拉格朗日插值多项式相比较，牛顿插值多项式在增加节点时不需要重新计算前面的各项，只需在后面增加相应的项，即牛顿插值法具有很好的继承性。

由差值多项式的存在唯一性定理知，满足同一组插值条件的拉格朗日多项式（2.3.4）与牛顿基本插值多项式（2.4.6）是同一多项式。因此，余项公式（2.3.9）也适用于牛顿插值。但是在实际计算中，有时也用差商表示的余项公式：

Rn(x)?f[x0,x1,来估计截断误差（证明略） 注意 上式中的n?1阶差商f[x0,x1,出f[x0,x1, ,xn,x]?n?1(x) （2.4.7）

故不能准确地计算,xn,x]与f(x)的值有关，

,xn,x]的精确值，只能对它做估计。

例2.4.6 当四阶差商变化不大时，可用f[x近0,x1,x2,x3,x4]似代替

f[x0,x1,x2,x3, x]2.5 求插值多项式的改进算法

2.5.1 分段低次插值

例2.3.2、例2.4.1表明适当地提高插值多项式的次数，有可能提高插值的

精度。但是决不可由此得出结论，认为插值多项式的次数越高越好。 例2.5.1 对函数

1f(x)?(?1?x?1) 21?25x2先以xi??1?i(i?0,1,5)为节点作五次插值多项式P5(x)，再以

51xi??1?i(i?0,1,1为0)节点作十次插值多项式P10(x)，并将曲线51f(x)?,y?5P(x),?y10P(x)?,(x?描绘在同一坐标系中，如图[1,1])2.5.121?2x5所示。可以看出，虽然在局部范围中，例如在区间[?0.2,0.2]中，P5(x)较10(x)比P好地逼近f(x),但从整体上看，P5(x)较好地逼近，尤其是在10(x)并非处处都比P区间[?1,1]的端点附近。进一步的分析表明，当n增大时，该函数在等距节点下

]1?,的高次插值多项式 ，在[1

两端会发生激烈的振荡。这种现象称为龙格（Runge)

32

现象。这表明，在大范围内使用高次插值，逼近的效果可能不理想。另一方面，插值误差除来自截断误差外，还来自初始数据的误差和计算过程中的舍入误差。插值次数越高，计算工作越大，积累误差也可能越大。

图2.5.1 多次插值的比较图

因此，在实际计算中，常用分段低次插值进行计算，即把整个插值区间分成若干小区间，在每个小区间上进行低次插值。 例2.5.2 当给定n?1个点x0?x1??xn上的函数值y0?y1??yn后，若要计

算点x?[xi?1,xi]处函数f(x)的近似值，可先选取两个节点xi?1,xi使x?[xi?1,xi]，然后在小区间[xi?1,xi]上作线性插值，即得

f(x)?P1(x)?yi?1x?xix?xi?1 （2.5.1） ?yixi?1?xixi?xi?1这种分段低次插值叫分段线性插值。在几何上就是用折线代替曲线，如图2.5.2

所示。故分段线性插值又称折线插值。

图2.5.2 折线插值几何直观图

类似地，为求y?f(x)的近似值，也可选取距点x最近的三个节点xi?1,xi,xi?1进行二次插值，即取

33

f(x)?P2(x)?k?i?1?[y?(xkj?i?1j?ki?1i?1x?xjk?xj)] （2.5.2）

这种分段低次插值叫分段二次插值。在几何上就是用分段抛物线代替曲线，故分段二次插值又称分段抛物插值。为了保证xi?1,xi,xi?1是距点x最近的三个节点，（2.5.2）中的i可通过下面方法确定：

1?1x?x??x1?x2?0?2?11?i??jx?x?x?xj?xj?1???j?1j?22?1?n?1?xn?2?xn?1??x?xn?2?

2.5.2 三次样条插值

分段低次插值虽然具有计算简单、稳定性好、收敛性有保证且易在计算机上

实现等优点，但它只能保证各小段曲线在连接点上的连续性，却不能保证整条曲线的光滑性，这就不能满足某些工程技术上的要求。

自六十年代开始，由于航空、造船等工程设计的需要而发展起来的所谓样条（Spline)的插值方法，既保留了分段低次插值多项式的各种优点，又提高了插值函数的光滑性。今天，样条插值方法已成为数值逼近的一个极其重要的分支，在许多领域里得到越来越广泛的应用。本节介绍应用最广泛且具有二阶连续导数的三次样条插值函数。 三次样条插值函数的定义

对于给定的函数表 x f(x) x0 y0 x1 y1 ? ? xn yn 其中a?x0?x1??xn?b，若函数S(x)满足：

,n)上是不高于三次的多项式；

（1）在每个子区间[xi?1,xi](i?1,2,（2）S(x),S'(x),S\x)在[a,b]上连续； （3）满足插值条件 S(xi)?yi(i?0,1,则称S(x)为函数f(x)关于节点x0,x1,,n)；

,xn的三次样条插值函数。

边界条件问题的提出与类型

单靠一张函数表是不能完全确定一个三次样条插值函数的。事实上，由条件

34

（1）知，三次样条插值函数S(x)是一个分段三次多项式，若用Si(x)表示它在第

i个子区间[xi?1,xi]上的表达式，则Si(x)形如：

Si(x)?ai0?ai1x?ai2x2?ai3x3,x?[xi?1,xi]这里有四个待定系数aij(j?0,1,2,3)。子区间共有n个，确定S(x)需要确定4n个待定系数。

另一方面，要求分段三次多项式S(x)及其导数S'(x),S\x)在整个插值区间只要在各子区间的端点xi(i?1,2,[a,b]上连续，

（3）可得待定系数应满足的4n?2个方程为:

故由条件（2），,n?1)连续即可。

?S(xi?0)?S(xi?0)(i?1,2,...,n?1)?S'(x?0)?S'(x?0)(i?1,2,...,n?1)ii??S''(x?0)?S''(x?0)(i?1,2,...,n?1)ii??S(xi)?yi(i?0,1,...,n)? (2.5.3)

由此可以看出，要确定个待定系数还缺少两个条件，这两个条件通常在插值区间

[a,b]的边界点a,b处给出，称为边界条件。边界条件的类型很多，常见的有：

''（1）给定一阶导数值S'(x0)?y0，称为第一种边界条件； ,S'(xn)?yn\\（2）给定二阶导数值S''(x0)?y0，称为第二种边界条件，特别地，,S''(xn)?yn满足自然边界条件的三次样条插值函数称S''(x0)?S''(xn)?0称为自然边界条件，为自然样条插值函数。

（3）当f(x)是周期为b?a的函数时，则要求S(x)及其导数都是以b?a为周期的函数，相应的边界条件称为第三种边界条件或周期边界条件：

S'(x0?0)?S'(xn?0),S''(x0?0)?S''(xn?0)三次样条插值函数的求法

。

虽然可以利用方程组（2.5.3）和边界条件求出所有待定系数aij，从而得到三次样条插值函数S(x)在各个子区间[xi?1,xi]的表达式Si(x)。但是，这种做法的计算工作量大，不便于实际应用。下面介绍一种简便的方法。

设在节点xi处S(x)的二阶导数为

S\xi)?Mi(i?0,1,

35

,n)

因为在子区间[xi?1,xi]上S(x)?Si(x)是不高于三次的多项式，其二阶导数必是线性函数（或常数）。于是，有

Si''(x)?Mi?1x?xix?xi?1?Mi,x?[xi?1,xi]

xi?1?xixi?xi?1记hi?xi?xi?1则有

Si''(x)?Mi?1xi?xx?xi?1 ?Mihihi连续积分两次得：

(xi?x)3(x?xi?1)3Si(x)?Mi?1?Mi?Ai?x?xi?1??Bi (2.5.3)

6hi6hi其中Ai,Bi为积分常数。利用插值条件

Si(xi?1)?yi?1,Si(xi)?yi

易得

1yi?yi?11B?y?Mi?1hi2 , Ai??(Mi?Mi?1)ii?16hi6将它们代入（2.5.3），整理得

(xi?x)3(x?xi?1)3Mx?xMx?xi?1Si(x)?Mi?1?Mi?(yi?1?i?1hi2)i?(yi?ihi2) 6hi6hi6hi6hix?[xi?1,xi],i?1,2,,n （2.5.4）

,n)，这n?1个值，就可定出三次

综合以上讨论可知，只要确定Mi(i?0,1,样条插值函数。

为了求出Mi(i?0,1,利用一阶导函数在子区间连接点上连续的条件即 ,n)，

S'(xi?0)?S'(xi?0)Si'(xi?0)?S'i?1(xi?0)得

Si??x???Mi?1 （2.5.5）

?xi?x?2h2?Mi?x?xi?1?2hi2 （2.5.6）

故

Si??xi?0??yi?yi?1hih?Mi?1?iMi （2.5.7） hi63 36

将（2.5.6）中的i改为i?1，即得S'(x)在子区间[xi,xi?1]上的表达式Si'?1(x)，并由此得

Si??1?xi?0??yi?1?yihi?1h?Mi?i?1Mi?1 (2.5.8) hi?136将（2.5.7）、（2.5.8）代入(2.5.5)整理后得

hih?hhy?yy?yi?1 Mi?1?ii?1Mi?i?1Mi?1?i?1i?i636hi?1hi两边同乘以

6，即得方程组

hi?hi?1,n?1?

hihi?16?yi?1?yiyi?yi?1?Mi?1?2Mi?Mi?1?????i?1,2,hi?hi?1hi?hi?1hi?hi?1?hi?1hi?若记

hihi?1???,???1??ii?ih?hh?hii?1ii?1??6 (2.5.9) ?g??f?x,x?fx,x??iii?1i?1i??hi?hi?1?????i?1,2,,n?1?则所得方程组可简写成

???iMi?1?2Mi??iMi?1?gi?i?1,2,,n?1?

即

?1M0?2M1??1M2?g1???2M1?2M2??2M3?g2? （2.5.10） ?????n?1`Mn?2?2Mn?1??n?1Mn?gn?1这是一个含有n?1个未知数、n?1个方程的线性方程组。要确定Mi的值，还需用到边界条件。

?和S??xn??yn?已知，在第(1)种边界条件下，由于S??x0??yo可以得到包含Mi另外两个线性方程。由(2.5.6)知，S(x)在子区间??x0,x1?上的导数为：

S1??x???M0?x1?x?2h12?M1?x?x0?2h12?y1?y0h1??M1?M0? h16?立即可得 故由条件S??x0??y0 37

???M0y0h1y1?y0h1???M1?M0? 2h16即

2M0?M1??6?y1?y0?（2.5.11） ?y?0?

h1?h1??，可得 同理，由条件S??xn??ynMn?1?2Mn?yn?yn?1?6??（2.5.12） y??n? hn?hn?,Mn的线性方

将(2.5.10)、(2.5.11)、(2.5.12)合在一起，即得确定M0,M1,程组：

?2???1?????12?1?n?121??M0??g0???M??g???1??1????? ??（2.5.13） ??????n?1??Mn?1??gn?1?2????Mn????gn??其中

6??g??x0,x1??y0?0hf??1 （2.5.14） ??g?6y??f?x,x?nn?n?1n?hn?????\??，Mn?S''?xn??yn在第(2)种边界条件下，由M0?S???x0??y0已知，在方

程组(2.5.14)中实际上只包含有n?1个未知数M1,M2,?2???2?????并且可以改写成： Mn?1，

?12?2?n?22?n?1?????M1??g1??1y0??M???g22??????????? （2.5.15） ??????n?2??Mn?2??gn?2????2????Mn?1????gn?1??n?1yn?在第(3)种边界条件下，由S??x0?0??S??xn?0?，直接可得

（2.5.16） M0?Mn

由条件S???x0?0??S???xn?0?可得

38

?M0hy?yn?1hnh1y1?y0h1???M`1?M0??Mnn?n??Mn?Mn?1? 2h162hn6注意到y0?yn和M0?Mn，上式整理后得：

hh16?y1?y0yn?yn?1?M1?nMn?1?2Mn???? h1?hnh1?hnh1?hn?h1hn?若记

hnh16,?n??1??n,gn??f?x1,xn??f?xn?1,xn?? h1?hnh1?hnh1?hn?n?则所得方程可简写成：

（2.5.17） ?nM1??nMn?1?2Mn?gn

将(2.5.10)(2.5.16)(2.5.17)合在一起，即得确定M1,M2, Mn的线形方程组：

?1??M1??g1??2?1????M??g?2?22???2??2???????? （2.5.18） ???????2?Mgn?1n?1n?1n?1????????n2???n???Mn????gn??利用线性代数知识，可以证明方程组(2.5.13)、(2.5.15)及(2.5.18)的系数矩阵都是非奇异的，从而都有唯一确定的解。

针对不同的边界条件，解相应的方程组(2.5.13)、(2.5.15)或(2.5.18)，求出M0,M2,，就可以得到S(x)在各子区间上的Mn的值，将它们代入（2.5.4）

表达式。

综上分析，有如下定理 定理2.5.1 对于给定的函数表 x y?f(x) x0 x1 ? xn y0 y1 ? yn (a?x0?x1?...?xn?b)，满足第一、第二和第三种边界条件的三次样条插

值函数是存在且唯一的。

三次样条插值函数S(x)的求解过程在下面的例子中给出了详细的说明。 例2.5.3 已知函数y?f(x)的函数值如下

x -1.5 0 1 2 y 0.125 -1 1 9

39

在区间[?1.5,2]上求三次样条插值函数S(x)，使它满足边界条件：

S'(?1.5)?0.75,S'(2)?14

解 先根据给定数据和边界条件算出?i,?i,gi，写出确定Mi的线性方程组。在本例中，给出的是第一种边界条件，确定Mi(i?0,1,2,3)的线性方程组，形如（2.5.13）。

由所给函数表知：

h1?1.5h2?1h3?1f[x0,x1]??0.75f[x1,x2]?2f[x2,x3]?8

于是由?i,?i,gi(i?1,2,,n?1)的算式（2.5.9）知：

?1?0.6,?2?0.5,?1?0.4,?2?0.5,g1?6.6,g2?18

由第一种边界条件下g0与gn的计算公式（2.5.14）知：

g0?6???6f?x0,x1??y0?h1??g3?6??f?y3x2,x3??36 ?h3??故确定M0,M1,M2与M3的方程组为：

100??M0???6??2?0.620.40??M??6.6????1???? （2.5.19） ?00.520.5??M2??18???????0012???M3??36?解此方程组得到S''(x)在各节点xi上的值Mi为：

M0??5,M1?4,M2?4,M3?16，

将Mi代入（2.5.4），即得S(x)在各子区间上的表达式Si(x)(i?1,2,由（2.5.4）知，S(x)在[x0,x1]上的表达式为：

S1?x??M0,n)。

?x1?x?6h13?M1?x?x0?6h13MM??x?x??x?x0??y0?0h12?1??y1?1h12?

6h6h????11在本例中，将

x0??1.5,代入，整理后得

x1?0,y0?0.125,y1??1,M0??5,M1?4

S1(x)?x3?2x2?1x?[?1.5,0]，

40

同理可得：

S2(x)?2x2?1x?[0,1] S3(x)?2x3?4x2?6x?3x?[1,2]

故所求三次样条插值函数为：

?x3?2x2?1(?1.5?x?0)???S?x???2x2?1(0?x?1)

?322x?4x?6x?3(1?x?2)???上述求三次样条插值函数的方法，其基本思路和特点是: 先利用一阶导数S?(x)在内节点xi(i?1,2,确定二阶导数Mi?S''?xi??i?0,1,.n?1)上的连续性以及边界条件，列出

,n?的线性方程组（在力学上称为三弯矩方程

组），并由此解出Mi，然后用Mi来表达S?x?。

通过别的途径也可求三次样条插值函数。例如，可以先利用二阶导数在内节

点上的连续性以及边界条件，列出确定一阶导数

mi?S??xi??i?0,1,,n?

的线性方程组（在力学上称为三转角方程组），并由此解出mi，然后用mi来表达

S?xi?。在有些情况下，这种表达方法与前者相比较，使用起来更方便。读者可以参阅相关书籍，在此不赘述。

2.6 求函数近似表达式的拟合法

在科学实验和生产实践中，经常要从一组实验数据(xi,yi)(i?1,2,,m)出发，

寻求函数y?f(x)的一个近似表达式y??(x)（称为经验公式）。从几何上，就是希望根据给出的m个点(xi,yi)，求曲线y?f(x)的一条近似曲线y??(x)。这就是曲线拟合的问题。

多项式插值虽然在一定程度上解决了由函数表求函数的近似表达式问题，但用它来解决这里提出的问题，有明显缺陷。

首先，实验提供的数据通常带有测试误差。如要求近似曲线y??(x)严格地通过所给的每个数据点(xi,yi)，就会使曲线保持原有的测试误差。当个别数据的误差较大时，插值效果显然是不理想的。

41

其次，由实验提供的数据往往较多（即m较大），用插值法得到的近似表达式，明显地缺乏实用价值。

因此，怎样从给定的一组数据出发，在某个函数类?中寻求一个“最好”的函数?(x)来拟合这组数据，是一个值得讨论的问题。

随着拟合效果“好”、“坏”标准的不同，解决此类问题的方法也不同。这里介绍一种最常用的曲线拟合方法，即最小二乘法。

2.6.1 曲线拟合的最小二乘法

如前所述，在一般情况下，我们不能要求近似曲线y?f(x)严格地通过所有数据点(xi,yi)，亦不能要求所有拟合曲线函数在xi处的偏差（亦称残差）

?i??(xi)?yi(i?1,2,,m)都严格地趋于零。但是，为了使近似曲线尽量反映所

给数据点的变化趋势，要求?i都较小还是必要的。达到这一目标的途径很多，常见的有：

(1)选取?(x)，使偏差绝对值之和最小，即

?????(x)?yiii?1i?1mmi（2.6.1） ?min

(2)选取?(x)，使偏差最大，绝对值最小，即

max?i?max?(xi)?yi?min （2.6.2）

1?i?m1?i?m(3)选取?(x)，使偏差平方和最小，即

??i?1m2i ??[?(xi)?yi]2?min （2.6.3）

i?1m为了方便计算、分析与应用，我们较多地根据“偏差平方和最小”的原则（称为最小二乘原则）来选取拟合曲线y??(x)。按最小二乘原则选择拟合曲线的方法，称为最小二乘法。

本章要着重讨论的线性最小二乘问题，其基本提法是：对于给定数据表

x x0 x1 ? xm y y0 y1 ? ym 要求在某个函数类??{?0(x),?1(x),,?n(x)}（其中n?m）中寻求一个函数

* ?an?n(x) （2.6.4）

?*??(x)?ao?(x)?a1?1(x)? 42

使?*(x)满足条件

?[?(x)?y]*iii?1m2（2.6.5） ?min?[?(xi)?yi]2

?(x)??i?1m式中，?(x)?a0?0(x)?a1?1(x)??an?n(x)是函数类?中任一函数。

满足关系式（2.6.5）的函数?*(x)，称为上述最小二乘问题的最小二乘解 。 由此可知，用最小二乘法解决实际问题包含两个基本环节：先根据所给数据点的变化趋势与问题的实际背景确定函数类?,即确定?(x)所具有的形式；然后按最小二乘法原则（2.6.3）求取最小二乘解?*(x)，即确定其系数

ak*(k?0,1,,n)。

由最小二乘解（2.6.4）应满足条件（2.6.5）知，点a0*,a1*,数

2,an*是多元函

S(a0,a1,**的极小值点，从而a0,a1,?n?,an)????ak?k(xi)?yi?

i?1?k?o?m*满足方程组： ,an?S?0，(k?0,1,2,?ak,n)

即

??(x)?a?(x)?a?(x)?...?a?(x)?y??0，(k?0,1,2,ki00i11inniii?1m,n)

亦即

a0??k(xi)?0(xi)?a1??k(xi)?1(xi)?...?an??k(xi)?n(xi)???k(xi)yi

i?1i?1i?1i?1mmmm若对任意的函数h(x)和g(x)，引入记号

(h,g)??h(xi)g(xi) （2.6.6）

i?1m则上述方程组可以表示成

a0(?k,?0)?a1(?k,?1)?...?an(?k,?n)?(?k,f)(k?0,1,写成矩阵形式即

43

,n)

?(?0,?0)(?0,?1)?(?,?)(?,?)11?10???(?,?)(?,?)n1?n0(?0,?n)??a??(?,f)?00?(?1,?n)??a??(?,f)??1???1? （2.6.7） ??????????a(?,f)?nn????(?n,?n)?方程组(2.6.7)称为法方程组。

事实上，最小二乘法的法方程可以用下面的方法形成。

在?(x)?a0?0(x)?a1?1(x)??an?n(x)中，当x?xi,i?1,2,,m时，令

?(xi)?yi，即得方程组

?a0?0(x1)?a1?1(x1)??an?n(x1)?y1?a?(x)?a?(x)??a?(x)?y?002112nn22 ????a0?0(xm)?a1?1(xm)??an?n(xm)?ym将其写成矩阵形式

??0(x1)?1(x1)???0(x2)?1(x2)????0(xm)?1(xm)??0(x1)?1(x1)??0(x2)?1(x2)?令 ????0(xm)?1(xm)???????????n(xm)??an??ym??n(x1)??a0??y0?????n(x2)??a1?????y1?

?a0??y0?????ay1?1????Am?(n?1)，??，?y，则方程组可写

???????????n(xm)?a?ym??n??n(x1)??n(x2)??为

Am?(n?1)??y

将方程两边同时乘以AT，则可得到

ATAm?(n?1)??ATy

这就是最小二乘法的法方程（2.6.7）。

当?0(x),?1(x),...,?n(x)线性无关时，可以证明它有唯一解

*** a0?a0,a1?a1,...,an?an并且相应的函数(2.6.4)就是满足条件(2.6.5)的最小二乘解。 定理2.6.1 对任意给定的一组实验数据(xi,yi)(i?1,2,数类

44

,m)(其中xi互异),在函

????0(x),?1(x),...?n(x)?(n?m)(?0,?1,...,?n线性无关)

中,存在唯一的函数

***?*?a0?0(x)?a1?1(x)?...?an?n(x)

使得关系式(2.6.5)成立，并且其系数ai*(i?0,1,...,n)可以通过解方程组(2.6.7)得到。

作为曲线拟合的一种常用的情况，若讨论的是代数多项式拟合，即取

?0(x)?1,?1(x)?x,...,?n(x)?xn，

则由(2.6.6)知：

(?j,?k)??xix??xij?k?j,k?0,1,jkii?1i?1mm,n?

(?k,f)??xikyi(k?0,1,i?1m,n)

故相应的法方程组为：

??m??m??xi?i?1?...?m?xni??i?1???xi?1mmi?xi?12i....?xi?1mn?1i??a??m?...?x?0y?i???i?1i?1??????mm???...?xin?1??a1???xiyi??i?1??? （2.6.8）??i?1 ???......??????m?m???xny???...?xi2n?an?????ii?i?1?????i?1?mni下面通过两个具体的例子来说明用最小二乘法解决实际的问题的具体步骤与某些技巧。

例2.6.1 某种铝合金的含铝量为x%，其熔解温度为yoc，由实验测得x与y的数据如表2.6.1左边三列。使用最小二乘法建立x与y之间的经验公式。 解 根据前面的讨论，解决问题的过程如下： （1）将表中给出的数据点(xi,yi)(i?1,2,如图2.6.1所示。 ,6)描绘在坐标纸上，

图2.6.1 数据的散点图

45

（2）确定拟合曲线的形式。由图2.6.1可以看出，六个点位于一条直线的附近，故可以选用线性函数（直线）来拟合这组实验数据，即令

?(x)?a?bx （2.6.9）

其中a,b为待定常数。

（3）建立法方程组。由于问题归结为一次多项式拟合问题，故由(2.6.8)知，相

??6应的法方程组形如：?6???xi?i?1??6?xy??ii??a??i?1?i?1?????? 66??2??b?xxy?i???ii?i?1??i?1?6经过计算（表2.6.1）即得确定待定系数a,b的法方程组：

6a?396.9b?1458? （2.6.10） ??396.6a?28365.28b?101176.3（4）解法方程组（2.6.10）得：

a?95.3524,b?2.2337

代入(2.6.9)即得经验公式：

y?95.3524?2.2337x （2.6.11）

表2.6.1 多项式拟合法方程系数计算表

i 1 2 3 4 5 6 xi 36.9 46.7 63.7 77.8 84.0 87.5 396.6 yi 181 197 235 270 283 292 1458 xi2 1361.61 2180.89 4057.69 6052.84 7056.00 7656.25 28365.28 xiyi 6678.9 9199.9 14969.5 21006.0 23772.0 25550.0 101176.3 ?所得经验公式能否较好地反映客观规律,还需通过实践来检验。由(2.6.11)式算出的函数值(称为拟合值)

yi?95.3524?2.2337xi(i?1,2,与实际值有一定的偏差，见表2.6.2。

由表2.6.2可以看出，偏差的平方和

,m)

??i?162i?26.6704

46

其平方根（称为均方误差）

??2i?5.164

在一定程度上反映了所得经验公式的好坏。

同时，由表2.6.2还可以看出，最大偏差为 max?i?3.22。

1?i?6表2.6.2 拟合偏差分析表 i 1 36.9 177.78 181 -3.22 10.37 2 46.7 199.67 197 2.67 7.13 3 63.7 237.64 235 2.64 6.97 4 77.8 269.13 270 -0.87 0.76 5 84.0 282.98 283 -0.02 0.0004 6 87.5 290.80 292 -1.20 1.44 xi yi yi ?i?yi?yi ?i2 ??2i26.6704 如果认为这样的误差是允许的话，就可以用经验公式(2.6.11)来计算含铝量在36.9%~87.5%之间的溶解度。否则，就要用改变函数类型或者增加实验数据等方法来建立新的经验公式。

例2.6.2 在某化学反应里，测得生成物的浓度y%与时间t的数据见表2.6.3，试用最小二乘法建立t与y的经验公式 。

表2.6.3 相应时间的浓度数据表 i 1 4.00 9 10.00 2 6.40 10 10.20 3 8.00 11 10.32 4 8.80 12 10.42 5 9.22 13 10.50 6 9.50 14 10.55 7 9.70 15 10.58 8 9.86 16 10.60 y t y 解 将已知数据点(ti,yi)(i?1,2,,16)描述在坐标纸上，得散点图2.6.2。

图2.6.2 数据的散点图

47

由图2.6.2及问题的物理背景可以看出，拟合曲线y??(x)应具下列特点： （1）曲线随着t的增加而上升，但上升速度由快到慢。

（2）当t?0时，反应还未开始，即y?0；当t??时，y趋于某一常数。故曲线应通过原点(或者当t?0时以原点为极限点)，且有一条水平渐近线。 具有上述特点的曲线很多。选用不同的数学模型，可以获得不同的拟合曲线与经验公式。

下面提供两种方案：

方案1 设想y??(t)是双曲线型的，并且具有形式

t （2.6.12） at?b此时，若直接按最小二乘法原则去确定参数a和b, 则问题归结为求二元函数

y?S(a,b)??(i?116ti（2.6.13） ?yi)2

ati?b的极小点，这将导致求解非线性方程组：

16??Sti2ti??2(?yi)?0??2?a(at?b)at?b?i?1ii ?16tt?Sii???2(?yi)?0?2?ati?bi?1(ati?b)??a给计算带来了麻烦。可通过变量替换来将它转化为关于待定参数的线性形函数。为此，将(2.6.12)改写成

1b?a? yt于是，若引入新变量

y(1)?1(1)1,t? yt则(2.6.12)式就是

y(1)?a?bt(1)

同时，由题中所给数据表2.6.3可以算出新的数据表2.6.4。这样，问题就归结为：根据数据表2.6.4，求形如y(1)?a?bt(1)的最小二乘解。参照例2.6.1的做法，解方程组

??16??16(1)??ti?i?1

??16(1)?t????yi?a??i?1?????i?1?，

1616?(1)2??b?(1)(1)?(t)t?i???iyi?i?1??i?1?(1)i1648

即得

a?80.6621,b?161.6822。

代入(2.6.12) ,即得经验公式

y?表2.6.4 算出数据表 i 1 2 0.50000 3 0.33333 … … 16 0.06250 t （2.6.14）

80.6621t?161.6822ti(1)?1 ti1 yi1.00000 yi(1)?0.25000 0.15625 0.12500 … 0.09434 方案2 设想y??(t)具有指数形式

y?aeb/t(a?0,b?0) （2.6.15）

为了求参数a和b时，避免求解一个非线形方程组，对上式两边取对数

blny?lna?，

t此时，若引入新变量

1y(2)?lny,t(2)?

t并记A?lna,B?b，则上式就是

y(2)?A?Bt(2)

又由表2.6.3可算出新的数据表2.6.5。

表2.6.5 算出数据表 i 1 2 0.50000 3 0.33333 … … 16 0.06250 ti(2)?1 ti1.00000 yi(2)?lnyi 1.38629 1.85630 2.07944 … 2.36085 于是将问题归为：根据数据表2.6.5，求形如y(2)?A?Bt(2)的最小二乘解。

参照方案1，写出相应的法方程组并解之，即得

A??4.4807,B??1.0567

A?0.01，13b2?5B?? 于是 a?e 49

故得另一个经验公式

y?0.011325e?1.0567t （2.6.16）

将两个不同的经验公式(2.6.14)和(2.6.16) 的均方误差和最大偏差进行比

较，见表2.6.6。从均方误差与最大偏差这两个不同角度看，后者均优于前者。因此，在解决实际问题时，常常要经过反复分析，多次选择，计算与比较，才能获得好的数学模型。

表 2.6.6 不同经验公式的误差比较表 经验公式 （3.9）式 （3.11) 式 均方误差 最大偏差 1.19?10?3 0.34?10?3 0.568?10?3 0.277?10?3 下面以常用的多项式拟合为例，说明最小二乘法在计算机上实现的步骤。 设有一组实验数据(xi,yi)(i?1,2,...,m)，今要用最小二乘法求一n(n?m)次多项式曲线

?n(x)?a0?a1x?...?anxn

来拟合这组数据。显然，求?n(x)的实质就是要确定其系数ai(i?1,2,...,n)。 由前面讨论可知，问题可归结为建立和求解法方程组(2.6.8)。为了便于编制程序并减少工作量，引入矩阵

?1x1x12x1n??y1???y?2n?1xxx222?,Y??2? C?????????2n?1xxx???ym?mmm??则法方程组(2.6.8)的系数矩阵（用A表示）和右端向量（用b表示）分别为:

A?CTC b?CTY。

2.6.2 加权最小二乘法

在实际问题中测得的所有实验数据，并不是总是等精度、等地位的。显然，对于精度较高或地位较重要（这应根据具体情况来判定）的那些数据(xi,yi)，应当给予较大的权。在这种情况下，求给定数据的拟合曲线，就要采用加权最小二乘法。

用加权最小二乘法进行曲线拟合的要求与原则：对于给定的一组实验数据

(xi,yi)(i?1,2,求一个函数

,m)，要求在某个函数类??{?0(x),?1(x),,?n(x)}(n?m)中，寻

**?*(x)?ao?o(x)?a1?1(x)?*?an?n(x)

50

使

?W[?(x)?y]*iiii?1m2?min?Wi[?(xi)?yi] 2?(x)??i?1m其中 ?(x)?ao?o(x)?a1?1(x)?为函数类?中任一函数；Wi(i?1,2,?an?n(x)

,m)是一列正数，称为权，它的大小反映了

数据(xi,yi)地位的强弱。显然，求?*(x)的问题可归结为求多元函数：

S(ao,a1,**的极小点(ao,a1,,an)??Wi[?ak?k(xi)?yi]2

i?1k?0mn*，但,an)。采用最小二乘解的求法，仍可得法方程组（2.6.7）

其中

(?k,?j)??Wi?k(xi)?j(xi)(k,j?0,1,i?1m,n)，

(?k,f)??Wi?k(xi)yi(k?0,1,i?1m,n)。

作为特例，如果选用的拟合曲线为?n(x)?ao?a1x?组为

?m??Wi?i?1?m??Wixi?i?1??m?Wxn??ii?i?1mm?anxn，那么相应的法方程

?Wxi?1mii?Wxi?12ii?Wxi?1mn?1ii??m?a??WxWy??o??ii???i?1??i?1???mm??n?1???aWxWxy?ii?1???iii? （2.6.17） i?1????i?1??????????m?m??n????Wixi2n????an???Wixiyi?i?1??i?1?nii例2.6.3 已知一组实验数据(xi,yi)及权Wi如表2.6.7。若x与y之间有线性关系，试用最小二乘法确定系数a和b。

表2.6.7 实验数据表

i 1 14 2 2 2 27 4 11 3 12 6 28 4 1 8 40 Wi Xi Yi

51

解 因为拟合曲线为一次多项式曲线（直线）

?1(x)?a?bx,

故相应的法方程组如（2.6.17）。将表中各已知数据代入，即得方程组

?54a?216b?701 ?216a?984b?3580?解之得

a??12.885,b?6.467

2.6.3 利用正交函数作最小二乘法拟合

在前几节，虽然从原则上解决了最小二乘意义下的曲线拟合问题，但在实际

计算中，由于当n较大，例如n?7，法方程组往往是病态的，因而给求解工作带来了一定困难。近年来，产生了许多解决这一困难的新方法。本节将简要介绍利用正交函数作最小二乘拟合的基本原理，以及利用正交多项式拟合的一种行之有效的方法。

对于{xi}和权{Wi}(i?1,2,,m)，若一组函数

?o(x),?1(x),,?n(x)(n?m)

满足条件：

k?j?0(?k,?j)??Wi?k(xi)?j(xi)??(k,j?0,1,i?1?Ak?0k?jm,n) （2.6.18）

则称?o(x),?1(x),,?n(x)是关于点集{xi}带权{Wi}的正交函数族。

n)都是多项式时，就称?o(x),?1(x),,?n(x)是关于点集

特别，当?k(x)(k?0,1,{xi}带权{Wi}的一组正交多项式。

如果在提到正交函数或正交多项式时，没有提到权Wi，就意味着权都是1。 若所考虑的函数类

??{?o(x),?1(x),,?n(x)}

,m)的正交函数族，则由条件

中的基函数是关于给定点集{xi}和权{Wi}(i?1,2,(2.6.18)知，在法方程组(2.6.7)的系数矩阵中，非对角线上元素

(?k,?j)?0(k?j)，

此时法方程组简化为：

52

?(?o,?o)?(?1,?1)??????ao??(?o,f)???a??(?,f)???1???1? （2.6.19） ??????????(?n,?n)??an??(?n,f)?只要由此解出

*ak?ak(k?0,1,,n)

就可得到最小二乘法解：

**?*(x)?ao?o(x)?a1?1(x)?*（2.6.20） ?an?n(x)

由条件(2.6.18) ，知

(?k,?k)?0(k?0,1,故易解方程组（2.6.19)，且有：

,n)，

(?,f)a?k?(?k,?k)*k?W?(x)yikii?1mikii?1mi(k?0,1, ,n) （2.6.21）

?W[?(x)]2这样就避免了求解病态方程组。

若函数类?的基函数?o(x),?1(x),,?n(x)是关于给定点集{xi}和权{Wi}的正

*交函数族，则可以直接由(2.6.21)式算出待定参数ak，进而写出最小二乘解

(2.6.20)。因此，问题归结为对给定的函数类?，寻求一组由正交函数族组成的基函数的问题。

构造正交函数组的方法很多。下面以多项式为例，介绍一种具体的方法。这种做法是以下述定理为基础的。

定理2.6.2 对于给定的点集{xi}和权{Wi}(i?1,2,,m)利用递推公式：

?o(x)?1???1(x)?x??1 （2.6.22） ???(x)?(x??)?(x)???(x)k?1kkk?1?k?1?k?1??Wx[?iii?1mmk(xi)]2 （2.6.23）

(xi)]2?W[?ii?1k?k??W[?(x)]ikim2?W[?ii?1i?1m(k?1,2,（2.6.24） ,n?1;n?m)

k?1(xi)]253

构造的函数族?o(x),?1(x),且?k(x)(k?1,2,,?n(x)是关于点集{xi}和权{Wi}的一组正交多项式,

,n)是k次多项式，其最高次项xk的系数为1。

例2.6.4 已知一组实验数据如表2.6.8，试用最小二乘法求一条二次拟合曲线。

表2.6.8 实验数据表

i 1 0.0 0.0 2 0.9 10.0 3 1.9 30.0 4 3.0 50.0 5 3.9 80.0 6 5.0 110.0 xi yi *解 采用边构造正交多项式?k(x)边求最小二乘解系数ak的顺序来求拟合曲线。

由公式(2.6.22)及(2.6.21)立即可得：

?o(x)?1，a*oy??6i?280?46.667 6由公式(2.6.23)、(2.6.22)及（2.6.21)依次可得：

a1??x6i?14.7?2.45； 6?1(x)?x?2.45；

a*1?(x)y???[?(x)]1i1ii2?392?22.254

17.615由公式(2.6.23), (2.6.22)、 (2.6.24)及(2.6.21)依次可得：

a2?x[?(x)]??[?(x)]i1i21i2?2244.359?2.5183；

17.61517.615?2.9358 6?1[?(x)]???[?(x)]1ioi??2(x)?(x?2.5183)(x?2.45)?2.9385?x2?4.9683x?3.2340

*a2???(x)y?[?(x)]2i2ii2?82.8926?2.247

36.8916故所求拟合曲线为：

***y?ao?o(x)?a1(x)?a2?2(x)?2.247x2?110.9x?0.5888

54

2.7 城市供水量预测的简单方法

2.7.1 供水量增长率估计与数值微分

供水量的增长率就是供水量函数的导数，因此，要估计供水量的增长率就需要求出供水量函数的导数。用供水量近似函数的导数作为其导数的近似估计是自然的想法。作为多项式插值的应用，本节介绍两种求函数导数近似值的方法。

2.7.2 利用插值多项式求导数

若函数f(x)在节点xi?i?0,1,就可作f(x)的n次插值,n?处的函数值已知，

多项式Pn(x)，并用Pn(x)近似代替f(x)，即

f(x)?Pn(x) (2.7.1)

由于Pn(x)是多项式，容易求导数，故对应于f?x?的每一个插值多项式Pn(x)，就易建立一个数值微分公式

f??x??P'n?x?，

这样建立起来的数值微分公式，统称为插值型微分公式。

必须注意，即使Pn(x)与f?x?的近似程度非常好，导数f'?x?与Pn'(x) 在某些点上的差别仍旧可能很大，因而，在应用数值微分公式时，要重视对误差

的分析。

由插值余项公式（2.3.9）知

f?n?1?????n?1?x?d?f?n?1????? （2.7.2）???f?x??pn?x???n?1?x?? ??n?1?!?n?1?!dx???x?作出估计。由于式中?是x未知函数，故x?xi时，无法利用上式误差f??x??pn但是，如果我们限定求某个节点xi处的导数值，那么（2.7.2）右端第二项之值应为零，此时有

f?n?1??????xi????1?xi? f??xi??pn?n?n?1?!若将它写成带余项的数值微分公式，即

f?n?1??????xi????1?xi? （2.7.3）f??xi??pn?n

n?1!??其中?在x0,x1组成。

该式右端由两部分，即导数的近似值和相应的截断误差xn之间。

55

由（2.7.3），当n?1时，插值节点为x0,x1，记h?x1?x0，得带余式的两点公式

f(x1)?f(x0)h??f(x)??f??(?)0??h2 ??[x0,x1] (2.7.4) ?f(x)?f(x)h10?f?(x)??f??(?)1??h2前一公式的实质是用f(x)在x0处的向前差商作为f?(x0)的近似值，后一公式则是用f(x)在x1处的向后差商作为f?(x1)的近似值。

当n?2且节点为xk?x0?kh(k?0,1,2)时,由(2.7.3)可得带余项的三点公式：

?1h2?f?(x0)?2h??3f(x0)?4f(x1)?f(x2)??3f???(?)?1h2?f?(x1)???f(x0)?f(x2)??f???(?)??[x0,x2] (2.7.5) ?2h6??1h2f???(?)?f?(x2)??f(x0)?4f(x1)?3f(x2)??2h3?中间一个公式的实质是用f(x)在x1处的中心差商作为f?(x1)的近似值，它与前后两公式相比较，其优越性是显然的。

用插值多项式Pn(x)作为f(x)的近似函数，还可建立高阶的数值微分公式。 例2.7.1 带余式的二阶三点公式

???1h2(4)?????f(x)?f(x)?2f(x)?f(x)??hf(?)?f(?)???001212??2h6????1h2(4)? (2.7.6) f??(x1)?2?f(x0)?2f(x1)?f(x2)??f(?2)?h12????1h2(4)?????f(x)?f(x)?2f(x)?f(x)?hf(?)?f(?)???201212??2h6????2.7.3 利用三次样条插值函数求导

由三次样条插值函数知，对于给定函数表 x x0 x1 ? xn y?f(x) y0 y1 ? yn (a?x0?x1?...?xn?b)和适当的边界条件，可以写出三次样条插值公式S(x)，

56

并用S(x)近似代替f(x),即

f(x)?S(x)，x?[a,b]

由于S(x)是一个分段三次多项式，在各子区间[xi?1,xi](i?1,2,导数，故可建立数值微分公式

f?(x)?Si'(x)?,n)上容易求出

hi122??M(x?x)?M(x?x)?f[x,x]?(Mi?1?Mi) (2.7.7) ii?1i?1ii?1i?2hi?61[Mi(x?xi?1)?Mi?1(xi?x)] x?[xi?1,xi],i?1,2,,n (2.7.8) hi1在节点xi??1?0.1i(i?0,1,1?25x2f??(x)?Si\(x)?利用函数f(x)?上的函数值和边,20)界条件S'(?1)?0.0740,S'(1)??0.0740，构造三次样条插值公式S(x)，并用它来计算f(x)和f?(x)在下列点：

xk??1?0.02k(k?0,1,2,处的近似值。计算结果见表2.7.1。

表2.7.1 近似值与准确值比较表 近似值 ,100) 准确值 x -1.00 -0.92 -0.84 -0.76 -0.68 -0.60 -0.52 -0.44 -0.36 -0.28 -0.20 -0.12 -0.04 S(x) 0.03846 0.04513 0.05365 0.06476 0.07961 0.1000 0.1289 0.1711 0.2359 0.3375 0.5000 0.7372 0.9594 S'(x) 0.074 0.09369 0.1209 0.1594 0.2125 0.3000 0.4319 0.6457 1.003 1.579 2.563 3.157 1.885 f(x) 0.03846 0.04513 0.05365 0.06477 0.07962 0.1000 0.1289 0.1712 0.2358 0.3378 0.5000 0.1353 0.9615 f'(x) 0.07639 0.09367 0.1209 0.1594 0.2155 0.3000 0.4318 0.6451 1.001 1.598 2.500 3.244 1.849 由表2.7.1可以看出，利用三次样条插值函数S(x)及其导数来逼近被插值函数f(x)及其导数，其效果很好。

57

2.7.4 城市供水量预测

现在来解决本章2.1.1节提出的城市供水量预测问题。 用插值方法预测07年1月份城市的用水量

预测2007年1月份城市的用水量有两种办法：一是先预测1月份每天的日用水量，求和后即得到1月用水总量；二是直接用2000-2006每年1月份的总用水量来预测。

（1）用06年每天的日用水量预测

以预测07年1月1日为例，这里仅用06年全部365天的日用水量作为插值节点(xi,yi)(i?1,2,...,365)，其中xi?i，表示第i天，yi即为第i天的日用水量（万吨/日）。其散点图如图2.7.1所示。设插值多项式为f(x)，则07年1月1日的日用水量的预测值即为f(x)在x?366处的函数值。

170165160日用水量(万吨/日)155150145140135050100150200250日期(从06年第一天开始)300350400

图2.7.1 06年城市日用水量散点图

分别用拉格朗日与牛顿插值方法得到364次插值多项式，并计算得f(366)分别为1.8781e+109、-9.3842e+167，结果的误差显然太大。实际中，当节点数较多时，一般用分段低次插值，现在用三次样条插值法求解，其样条曲线如图2.7.2所示。对应的插值结果为133.6641。

但若用这个三次样条函数来估计f(x),x?367,...,396，即该月后30天的日用水量，结果为98.1477、33.0292、…、-4.8837，误差逐渐增大。由此可见，仅用06年一年数据预测07年1月份每天的日用水量时，即便使用低次插值，误差也会很大。

（2）用2000-2006年1月份每天的日用水量预测

现用2000-2006年1月每天的日用水量预测07年1月对应天的日用水量。同样，以预测07年1月1日为例。以2000-2006每年1月1日的日用水量作为插值节点(xi,yi)(i?1,2,...,7)，则f(8)即为所要求的07年1月1日用水量的估计值。同理，其它30天的日用水量也对应地求出。

58

170165160日用水量(万吨/日)155150145140135130050100150200250日期(从06年第一天开始)300350400

图2.7.2 三次样条曲线

采用拉格郎日、牛顿和三次样条插值函数进行计算，三种方法计算的结果基本相同，求得07年1月每天日用水量分别为143.6962，139.7505,?,146.5117。对它们求和得07年1月份总用水量的预测值为4517.6993（万吨/日）。其计算结果如图2.7.3所示（”+”号所示为217个已知的前七年日用水量节点，”*”号所示为预测的31天的用水量值）。

拉格朗日插值155150145155150145牛顿插值155150145样条插值日用水量(万吨/日)日用水量(万吨/日)140135130125120140135130125120日用水量(万吨/日)1401351301251200200400日期(前七年1月每日)0200400日期(前七年1月每日)0200400日期(前七年1月每日)

图2.7.3 三种插值函数曲线

（4）用2000-2006年1月城市的总用水量预测

由表2.1.2得到7个插值节点(xi,yi)，其中xi?i,i?1,2,...,7，其散点图如2.7.4所示。用三次样条插值得f(8)即所求的07年1月用水总量的估计值为4378.1390（万吨/日）。同理得到的07年1月用水总量估计值如下表2.7.2所示。

表2.7.2 插值法得到07年1月用水总量估计值（万吨/日） 插值节点 预测结果 06年日用水量 误差大 00-06年1月每天日用水量 4517.6993（万吨/日） 00-06每年1月用水总量 4378.1390（万吨/日） 59

46004500月用水总量(万吨/日)4400430042004100400012345年份(从2000年开始)67

图2.7.4 散点图

用数据拟合方法预测2007年1月份城市的用水量 （1）用06年每天的日用水量预测

由06年全部365天的日用水量散点图2.7.1可知，这些点并不简单地呈现线性或二次关系，但有很强的聚集性。我们试图用几个多项式进行拟合。用Matlab工具箱得到如表2.7.3拟合结果。

表2.7.3 拟合次数与残差的关系 均方根误差 线性 5.366 二次 4.689 三次 4.55 四次 4.4 五次 3.993 六次 3.897 七次 3.781 八次 3.731 九次 3.701 拟合结果显示，九次拟合的均方根误差是最小的，但四次曲线拟合的效果已经与它很接近了，拟合曲线如图2.7.5所示。

图2.7.5 四次拟合曲线

拟合函数的表达式为：

f(x)?1.389e-008x4 -1.142e-005x3+0.00282x2-0.18x+149.7。

计算f(x)在x?366,,396处的函数值分别得到07年1月31天的日用水量为

150.8524，150.8718，?，152.9422。对它们求和得到1月份总用水量的估计值为4700.4297（万吨/日）。

（2）用2000-2006年1月份每天的日用水量预测

由2000-2006年1月1日的日用水量得到7个插值节点，利用拟合函数计算

f(8)即07年1月1日日用水量估计值。

60

因插值节点只有7个，若用6次拟合曲线刚好可过所有插值节点。图2.7.6显示的是六次拟合结果。

图2.7.6六次拟合曲线

但是计算f(8)得07年1月1日日用水量为5.74456，产生了较大的误差。

通过1到5次拟合曲线的均方根误差指标值可知，用线性拟合结果最佳。拟合曲线表达式为f(x)=3.5088x +117.6164，计算f(8)，得07年1月1日日用水量估计值为145.6869。拟合曲线如图2.7.7所示

图2.7.7线性拟合曲线

同理用线性拟合得到第2到31天日用水量估计值为143.7937、?、148.3227。相加得到1月份用水总量的估计值为4654.8538（万吨/日）。 （3）用2000-2006年1月城市的总用水量预测

同理，尽管六次拟合曲线能完全拟合已知节点，但将对计算点产生较大误差，选用三次拟合，其拟合曲线如图2.7.8所示。

图2.7.8三次拟合曲线

61

拟合函数为：f(x)? 0.7678x3 -20.47 x2+ 201.3x +3854，计算得f(8)的近似值为4547.0724（万吨/日）。

上述拟合方法得到的07年1月用水总量估计值如表2.7.4所示：

表2.7.4 拟合法得到07年1月用水总量估计值（万吨/日） 数据 结果 06年日用水量 4700.4297（万吨/日） 00-06年1月每天日用水量 4654.8538（万吨/日） 00-06每年1月用水总量 4547.0724（万吨/日） 习 题2

1.求经过下列已知点的最低次代数多项式

x 0 1.5 5.1 y -1 4.25 35.21 2.已知函数表如下

x 10 11 12 13 lnx 2.3026 2.3979 2.4849 2.5649 试分别用线性插值与二次插值计算ln(11.75)的近似值，并估计截断误差 3.设x0,x1,,xn为任意给定的n?1个互不相同的节点，证明：

(1)若f(x)为不高于n次的多项式，则f(x)关于这组节点的n次插值多项式就是它自己；

（2）若lk(x)(k?0,1,,n)是关于这组节点的n次基本插值多项式，则有恒等式

mmxl(x)?x,(m?0,1,?kkk?0n,n)

4.已知函数表如下

x 0.0 0.2 0.4 0.6 0.8 ex 1.0000 1.2214 1.4918 1.8221 2.2255 (1)分别构造向前差分表与向后差分表；

(2)分别用三点与四点前插公式计算e0.13的近似值，并估计误差； (3)分别用三点与四点前插公式计算e0.72的近似值，并估计误差；

(4)构造差商表，并分别用三点与四点牛顿基本插值公式计算e0.12的近似值．

62

5.设f(x)为n次多项式，试证明当k?n时差商f[x,x0,x1,互异)为n?k次多项式，而当k?n时其值恒为零．

,xk](其中x0,x1,,xk6.今要在区间[?4,4]上构造f(x)?ex在等距节点下的函数表．问怎样选取函数表的步长，才能保证用二次插值求ex的近似值时，其截断误差不超过106 7.对于给定的插值条件

x 0 1 2 3 y 0 0 0 0

试分别求满足下列边界条件的三次样条插值函数： (1)S''(0)?1, S''(3)?0 (2)S'(0)?1, S'(3)?0 8.已知一组实验数据如下：

x 2 4 6 8 y 2 11 28 40 试用最小二乘法求一次和二次拟合多项式，并分别算出均方误差与最大偏差。 9.试根据下面五组测试数据，用最小二乘法求出一个经验公式，并计算均方误差。

xi 293 313 343 363 383 yi 28.98 30.2 32.9 35.9 38.8 10.用最小二乘法求一形如W?Ct?的经验公式(其中C和?为待定系数)，使之与下列数据相拟合。

ti 1 2 4 8 16 32 64 Wi 4.22 4.02 3.85 3.59 3.44 3.02 2.59 11.用最小二乘法求一形如y?a?bx2的多项式，使之与下列数据相拟合。

Xi 19 25 31 38 44 Yi 19.0 32.3 49.0 73.3 97.8

63

12.利用正交多项式，对第9题给出的实验数据，拟合一条二次曲线。

利用SPSS软件对以上数据进行二次最优拟合，拟合优度分别达到99.8%、99.1%，两个厂的拟合图形分别如图2.1.1和图2.1.2.根据以上处理后的数据，利用MATLAB软件产生了矩阵，然后进行二次拟合，得到了精度为10的两水厂一月份供水量g1,g2（吨）与时间x（年）的变化函数表达式： 一水厂月供水量g1函数关系式：

g1??6.9173.41536x2?1040417.918x?23982938.23

二水厂月供水量g2函数关系式：

g2??44109.84679x2?664616.9561x?14816050.48

64

本文来源：https://www.bwwdw.com/article/foor.html