科学计算与数学建模课件.ppt

科学计算与数学建模,中南大学数学科学与计算技术学院,城市供水量的预测模型,第2章 城市供水量的预测模型 插值与拟合算法,2.1 城市供水量的预测问题,2.1.1 实际问题与背景 为了节约能源和水源，某供水公司需要根据日供水量记录估计未来一时间段（未来一天或一周）的用水量，以便安排未来（该时间段）的生产调度计划。现有某城市7年用水量的历史记录，记录中给出了日期、每日用水量（吨/日）。如何充分地利用这些数据建立数学模型，预测2007年1月份城市的用水量，以制定相应的供水计划和生产调度计划。,表2.1.1某城市7年日常用水量历史记录（万吨/日）,表2.1.22000-2006年1月城市的总用水量（万吨/日）,利用这些数据，可以采用时间序列、灰色预测等方法建立数学模型来预测2007年1月份该城市的用水量。如果能建立该城市的日用水量随时间变化的函数关系，则用该函数来进行预测非常方便。但是这一函数关系的解析表达式是没办法求出来的，那么能否根据历史数据求出该函数的近似函数呢？根据未知函数的已有数据信息求出其近似函数的常用方法有插值法和数据拟合。本章将介绍插值法和数据拟合，并用这两种方法给出该城市供水量进行预测。,2.2 求未知函数近似表达式的插值法,2.2.1 求函数近似表达式的必要性 一般地，在某个实际问题中，虽然可以断定所考虑的函数在区间上存在且连续，但却难以找到它的解析表达式，只能通过实验和观测得到该函数在有限个点上的函数值（即一张函数表）。显然，要利用这张函数表来分析函数的性态，甚至直接求出其它一些点上的函数值是非常困难的。在有些情况下，虽然可以写出函数的解析表达式，但由于结构相当复杂，使用起来很不方便。面对这些情况，总希望根据所得函数表（或结构复杂的解析表达式），构造某个简单函数作为的近似。插值法是解决此类问题的一种比较古老的、然而却是目前常用的方法，它不仅直接广泛地应用于生产实际和科学研究中，而且也是进一步学习数值计算方法的基础。,定义2.2.1 设函数 在区间 上连续，且在 个不同的点 上分别取值，在一个性质优良、便于计算的函数类 中，求一简单函数，使(2.2.1)而在其它点 上作为 的近似。称区间 为插值区间，点 为插值节点，称（2.2.1）为 的插值条件，称函数类 为插值函数类，称 为函数在节点 处的插值函数。求插值函数 的方法称为插值法。插值函数 类 的取法不同，所求得的插值函数 逼近 的效果就不同，它的选择取决于使用上的需要。常用的有代数多项式、三角多项式和有理函数等。当选用代数多项式作为插值函数时，相应的插值问题就称为多项式插值。,在多项式插值中，求一次数不超过 的代数多项式(2.2.2)使(2.2.3)其中 为实数。满足插值条件(2.2.3)的多项式(2.2.2)，称为函数 的 次插值值多项式。次插值多项式 的几何意义：过曲线 上的 个点 作一条 次代数曲线 作为曲线 的近似，如图2.2.1所示。,图2.2.1,2.2.2插值多项式存在唯一性与求插值多项式的解方程组方法 由插值条件（2.2.3）知，的系数 满足线性方程组,由线性代数知，线性方程组的系数行列式是 阶范德蒙（Vandermonde）行列式，且,（2.2.4）,因 是区间 上的不同点，上式右端乘积中的每一个因子，于是系数行列式不等于0，即方程组（2.2.4）的解存在且唯一。从而得出下面的结论：定理2.2.1 若节点 互不相同，则满足插值条件（2.2.3）的次插值多项式（2.2.2）存在且唯一。,在上一节里，不仅指出了插值多项式的存在唯一性，而且也提供了它的一种求法，即通过解线性方程组（2.2.4）来确定其系数。但是，当未知数个数多时，这种做法的计算工作量大，不便于实际应用，Lagrange插值法就是一种简便的求法。,2.3 求插值多项式的Lagrange法,2.3.1 Lagrange插值基函数 先考虑一下简单的插值问题：对节点 中任意一点 做一 次多项式 使它在该点上取值为1,而在其余点 上取值为零,即（2.3.1）表明 个点 都是 次多项式 的零点，故可设 其中 为待定系数，由条件 可得,对应于每一节点 都能求出一个满足插值条件（2.3.1）的 次插值多项式（2.3.2），这样，由（2.3.2）式可以求出 个 次插插多项式。容易看出，这组多项式仅与节点的取法有关，称它们为在 个节点上的 次基本插值多项式或 次插值基函数。,（2.3.2）,2.3.2 Lagrange（拉格朗日）插值多项式 利用插值基函数立即可以写出满足插值条件（2.2.3）的 次插值多项式（2.3.3）事实上，由于每个插值基函数 都是 次多项式，故其线性组合（2.3.3）必是不高于 次的多项式，同时，根据条件（2.2.1）容易验证多项式（2.3.3）在节点 处的值 为，因此，它就是待求的 次插值多项式。形如(2.3.3)的插值多项式称为拉格朗日插值多项式，记为,（2.3.4）,令，由（2.3.4）即得两点插值公式（2.3.5）即（2.3.6）这是一个线性函数，用线性函数 近似代替函数，在几何上就是通过曲线 上两点 做一直线 近似代替曲线（见图2.3.1)，故两点插值又名线性插值。,令，由（2.3.4）又可得常用的三点插值公式：这是一个二次函数，用二次函数近似代替函数，在几何上就是通过曲线 上的三点 作一抛物线，近似地代替曲线（图2.3.1），故称为三点插值(二次插值)。,（2.3.7）,图2.3.1,例2.3.1 已知 分别用线性插值和抛物插值求 的值。解 因为115在100和121之间，故取节点，相应地有，于是，由线性插值公式（2.3.5）可得 故用线性插值求得的近似值为：,图2.3.2,仿上，用抛物插值公式（2.3.7）所求得的近似值为：将所得结果与 的精确值 相比较，可以看出抛物插值的精确度较好。为了便于上机计算，我们常将拉格朗日插值多项式（2.3.4）改写成公式（2.3.8）的对称形式,（2.3.8）,编程框图如图2.3.3，可用二重循环来完成 值的计算，先通过内循环，即先固定，令 从0到 累乘求得 然后再通过外循环，即令 从0到，累加得出插值结果,图2.3.3,2.3.3 插值余项 在插值区间 上用插值多项式 近似代替，除了在插值节点 上没有误差以外，在其他点上一般有存在误差的。若记 则 就是用 近似代替 时所产生的截断误差，称为插值多项式 的余项。关于误差有如下定理2.3.1中的估计式。定理2.3.1 设 在区间 上有直到 阶导数,为区间 上 个互异的节点，为满足条件:的 次插值多项式，则对于任何 有,（2.3.9）,其中 且依赖于证明 由插值条件 知，即插值节点都是 的零点,故可设(2.3.10)其中 为待定函数。下面求。对区间 上异于 的任意一点 作辅助函数：不难看出具有如下特点（1）（2）在 上有直到 阶导数，且,(2.3.11),（2.3.12）,等式（2.3.11）表明 在 上至少有 个互异的零点，根据罗尔(Rolle)定理，在 的两个零点之间，至少有一个零点，因此，在 内至少有 个互异的零点，对 再应用罗尔定理，推得 在 内至少有 个互异的零点。继续上述讨论，可推得 在 内至少有一个零点，若记为，则，于是由（2.3.12）式得将它代入(2.3.10)即得（2.3.9）。对于，（2.3.9）显然成立。,例2.3.2 在例2.3.2中分别用线性插值和抛物插值计算了 近似值，试估计它们的截断误差。解 用线性插值求 的近似值，其截断误差由插值余项公式（2.3.9）知 现在，故,当用抛物插值求 的近似值时，其截断误差为 现在，代入，即得,2.3.4 插值误差的事后估计法 在许多情况下，要直接应用余项公式（2.3.9）来估计误差是很困难的，下面将以线性插值为例，介绍另一种估计误差的方法。设 且 为已知，若将用 两点做线性插值求得 的近似值为，用 两点作线性插值所求得 的近似值记为，则由余项公式（2.3.9）知：假设 在区间 中变化不大，将上面两式相除，即得近似式,即 近似式（2.3.13）表明，可以通过两个结果的偏差 来估计插值，这种直接利用计算结果来估计误差的方法，称为事后估计法。,（2.3.13）,在例2.3.1中，用 做节点，算得的 近似值为，同样，用 做节点，可算得 的另一近似值。通过（2.3.13）可以估计出插值 的误差为：,2.4 求插值多项式的Newton法,由线性代数可知，任何一个不高于 次的多项式，都可表示成函数 的线性组合，即可将满足插值条件 的 次多项式写成形式 其中 为待定系数。这种形式的插值多项式称为牛顿Newton插值多项式，我们把它记成,即 因此，牛顿插值多项式 是插值多项式 的另一种表示形式,与拉格朗日插值多项式相比较，不仅克服了“增加一个节点时整个计算机工作必须重新开始”见例2.3.1的缺点，而且可以节省乘除法运算次数。同时，在牛顿插值多项式中用到的差分与差商等概念，又与数值计算的其它方面有着密切的关系.,（2.4.1）,2.4.1 向前差分与Newton（牛顿）向前插值公式 设函数 在等距节点 处的函数值 为已知，其中 是正常数，称为步长，称两个相邻点 和 处函数值之差 为函数 在点 处以 为步长的一阶向前差分简称一阶差分，记，即 于是，函数 在各节点处的一阶差分依次为 又称一阶差分的差分 为二阶差分。一般地，定义函数 在点 处的 阶差分为：为了便于计算与应用，通常采用表格形式计算差分，如表2.4.1所示。,在等距节点 情况下，可以利用差分表示牛顿插值多项式2.4.1 的系数，并将所得公式加以简化。事实上，由插值条件 即可得。,表2.4.1,再由插值条件 可得：由插值条件 可得：一般地，由插值条件 可得:于是，满足插值条件的插值多项式为：,令 并注意到 则可简化为 这个用向前差分表示的插值多项式，称为牛顿向前插值公式，简称前插公式。它适用于计算表头 附近的函数值。由插值余项公式（2.3.9），可得前插公式的余项为：,（2.4.2）,（2.4.3）,例2.4.1 从给定的正弦函数表表2.4.2左边两列出发计算,并估计截断误差。,表2.4.2,解 因为0.12介于0.1与0.2之间，故取，此时 为求 构造差分表2.4.2，表中长方形框中各数依次为 在 处的函数值和各阶差分。若用线性插值求 的近似值，则由前插公式（2.4.2）立即可得 用二次插值得：,用三次插值得：因 与 很接近，且由差分表2.4.2可以看出，三阶差分接近于常数（即 接近于零），故取 作为 的近似值，此时由余项公式（2.4.3）可知其截断误差为：,2.4.2 向后差分与Newton（牛顿）向后插值公式 在等距节点 下，除了向前差分外，还可引入向后差分和中心差分，其定义和记号分别如下：在点 处以 为步长的一阶向后差分和 阶向后差分分别为：在点 处以 为步长的一阶中心差分和 阶中心差分分别为：,其中，各阶向后差分与中心差分的计算，可通过构造向后差分表与中心差分表来完成（参见表2.4.2）。利用向后差分，可简化牛顿插值多项式（2.4.1），导出与牛顿前插公式（2.4.2）类似的公式，即：若将节点的排列次序看作，那么2.4.1)可写成：根据插值条件 得到一个用向后差分表示的插值多项式：,其中t0，插值多项式(2.4.4)称为牛顿向后插值公式，简称后插公式。它适用于计算表尾 附近的函数值。由插值余项公式（2.3.9），可写出后插公式的余项为：例2.4.2 已知函数表同例2.4.1，计算，并估算其截断误差。,（2.4.4）,（2.4.5）,解 因为0.58位于表尾 附近，故用后插公式（2.4.4）计算 的近似值。一般的，为了计算函数在 处的各阶向后差分，应构造向后差分表。但由向前差分与向后差分的定义可以看出，对同一函数表来说，构造出来的向后差分表与向前差分表在数据上完全相同。因此，表（2.4.2）用“”线标出的各数依次给出了 在 处的函数值和向后差分值。因三阶向后差分接近于常数，故用三次插值进行计算，且 于是由向后差分公式（2.4.4）得：,因为在整个计算中，只用到四个点 上的函数值，固由余项公式（2.4.5）知其截断误差为：,2.4.3 差商与牛顿基本插值多项式 当插值节点非等距分布时，就不能引入差分来简化牛顿插值多项式，此时可用差商这个新概念来解决。设函数 在一串互异的点 上的值依次为 我们称函数值之差 与自变量之差 的比值 为函数 关于 点的一阶差商，记作,例2.4.3 称一阶差商的差商 为函数 关于点 的二阶差商（简称二阶差商），记作，例2.4.4 一般地，可通过函数 的 阶差商定义的 阶差商如下：差商计算也可采用表格形式（称为差商表），如表2.4.3所示：,二阶差商,表2.4.3,三阶差商,差商具有下列重要性质（证明略）：,（1）函数 的 阶差商 可由函数值 的线性组合表示，且（2）差商具有对称性，即任意调换节点的次序，不影响差商的值。例如（）当 在包含节点 的某个区间上存在时，在 之间必有一点，使（4）在等距节点 情况下，可同时引入 阶差分与差商，且有下面关系：,引入差商的概念后，可利用差商表示牛顿插值多项式（2.4.1）的系数。事实上，从插值条件出发，可以像确定前插公式中的系数那样，逐步地确定（2.4.1）中的系数 故满足插值条件 的 次插值多项式为：（2.4.6）（2.4.6）称为牛顿基本插值多项式，常用来计算非等距节点上的函数值。,例2.4.5 试用牛顿基本差值多项式按例1要求重新计算的近似值。解 先构造差商表 由上表可以看出牛顿基本插值多项式（2.4.6）中各系数依次为：,故用线性插值所得的近似值为：用抛物插值所求得的近似值为：所得结果与例2.4.1一致。比较例2.4.1和例2.4.5的计算过程可以看出，与拉格朗日插值多项式相比较，牛顿差值多项式的优点是明显的。由差值多项式的存在唯一性定理知，满足同一组差值条件的拉格朗日多项式（2.3.4）与牛顿基本差值多项式（2.4.6）是同一多项式。因此，余项公式（2.3.9）也适用于牛顿插值。但是在实际计算中，有时也用差商表示的余项公式：（2.4.7）来估计截断误差（证明略）注意 上式中的 阶差商 与 的值有关，故不能准确地计算出 的精确值，只能对它做一种估计。,例2.4.6 当四阶差商变化不大时，可用近似代替,2.5 求插值多项式的改进算法,2.5.1 分段低次插值例2.3.2、例2.4.1表明适当地提高插值多项式的次数，有可能提高计算结果的准确程度。但是决不可由此得出结论，认为插值多项式的次数越高越好。例2.5.1 对函数先以 为节点作五次插值多项式，再以 为节点作十次插值多项式，并将曲线,描绘在同一坐标系中。如图2.5.1所示：,虽然在局部范围中，例如在 区间中，比 较好地逼近,但从整体上看，并非处处都比 较好地逼近，尤其是在 区间的端点附近。进一步的分析表明，当 增大时，该函数在等距节点下的高次插值多项式，在 两端会发生激烈的振荡。这种现象称为龙格（Runge)现象。这表明，在大范围内使用高次插值，逼近的效果可能不理想。另一方面，插值误差除来自截断误差外，还来自初始数据的误差和计算过程中的舍入误差。插值次数越高，计算工作越大，积累误差也可能越大。因此，在实际计算中，常用分段低次插值进行计算，即把整个插值区间分成若干小区间，在每个小区间上进行低次插值。,例2.5.2 当给定 个点 上的函数值 后，若要计算点 处函数值 的近似值，可先选取两个节点 使，然后在小区间上作线性插值，即得 这种分段低次插值叫分段线性插值。在几何上就是用折线代替曲线，如图2.5.2所示。故分段线性插值又称折线插值.,类似地，为求 的近似值.也可选取距点最近的三个节点 进行二次插值，即取 这种分段低次插值叫分段二次插值。在几何上就是用分段抛物线代替曲线，故分段二次插值又称分段抛物插值。为了保证 是距点最近的三个节点，（2.5.2）中 的可通过下面方法确定：用图表示为：,2.5.2 三次样条插值 分段低次插值虽然具有计算简单、稳定性好、收敛性有保证且易在电子计算机上实现等优点，但它只能保证各小段曲线在连接点上的连续性，却不能保证整条曲线的光滑性(如图2.5.2中的折线)，这就不能满足某些工程技术上的要求。从六十年代开始，首先由于航空、造船等工程设计的需要而发展起来的所谓样条（Spline)的插值方法，既保留了分段低次插值多项式的各种优点，又提高了插值函数的光滑性。今天，样条插值方法已成为数值逼近的一个极其重要的分支，在许多领域里得到越来越广泛的应用。本节介绍应用最广泛且具有二阶连续导数的三次样条插值函数。三次样条插值函数的定义 对于给定的函数表,其中，若函数 满足：（1）在每个子区间 上是不高于三次的多项式；（2）在 上连续；（3）满足插值条件 则称 为函数关于节点 的三次样条插值。边界条件问题的提出与类型 注：单靠一张函数表是不能完全确定一个三次样条插值函数的。事实上，由条件（1）知，三次样条插值函数 是一个分段三次多项式，若用 表示它在第 个子区间 上的表达式，则形如：这里有四个待定系数。子区间共有 个，确定 需要确定 个待定系数。,另一方面，要求分段三次多项式 及其导数 在整个插值区间 上连续，只要在各子区间的端点 连续即可。故由条件（2），（3）可得待定系数应满足的 个方程为:(2.5.3)由此可以看出，要确定个待定系数还缺少两个条件，这两个条件通常在插值区间 的边界点 处给出，称为边界条件。边界条件的类型很多，常见的有：,（1）给定一阶导数值（2）给定二阶导数值 特别地，称为自然边界条件，满足自然边界条件的三次样条插值函数称为自然样条插值函数。（3）当 是周期为 的函数时，则要求 及其导数都是以 为周期的函数，相应的边界条件为三次样条插值函数的求法 虽然可以利用方程组（2.5.3）和边界条件求出所有待定系数，从而得到三次样条插值函数 在各个子区间 的表达式。但是，这种做法的计算工作量大，不便于实际应用。下面介绍一种简便的方法。,设在节点 处 的二阶导数为 因为在子区间 上 是不高于三次的多项式，其二阶导数必是线性函数（或常数）。于是，有 记 则有 连续积分两次得：(2.5.3),其中 为积分常数。利用插值条件易得：将它们代入（2.5.3），整理得（2.5.4）综合以上讨论可知，只要确定，这 个值，就可定出三次样条插值函数。,为了求出，利用一阶函数在子区间连接点上连续的条件即（2.5.5）得（2.5.6）故（2.5.7）将（2.5.6）中的 改为，即得 在子区间 上的表式，并由此得：(2.5.8),将（2.5.7）、（2.5.8）代入(2.5.5)整理后得 两边同乘以，即得方程组 若记(2.5.9),可得则所得方程组可简写成 即（2.5.10）这是一个含有 个未知数、个方程的线性方程组。要确定 的值，还需用到边界条件。在第(1)种边界条件下，由于 和 已知，可以得到包含 另外两个线性方程。由(2.5.6)知，在子区间 上的导数为：,故由条件 立即可得 即（2.5.11）同理，由条件，可得（2.5.12）将(2.5.10)、(2.5.11)、(2.5.12)合在一起，即得确定 的线性方程组：,（2.5.13）其中（2.5.14）在第(2)种边界条件下，由，已知，在方程组(2.5.14)中实际上只包含有 个未知数，并且可以改写成：,在第(3)种边界条件下由，直接可得（2.5.16）由条件 可得 注意到 和，上式整理后得：,（2.5.15）,若记 则所得方程可简写成：（2.5.17）将(2.5.10)(2.5.16)(2.5.17)合在一起，即得确定 的线形方程组：（2.5.18）,利用线性代数知识，可以证明方程组(2.5.13)、(2.5.15)及(2.5.18)的系数矩阵都是非奇异的，从而都有唯一确定的解。针对不同的边界条件，解相应的方程组(2.5.13)、(2.5.15)或(2.5.18)，求出的值，将它们代入（2.5.4），就可以得到在各子区间上的表达式。综上分析，有 定理2.5.1 对于给定的函数表 满足第（1）或第（2）或第（3）种边界条件的三次样条插值函数是存在且唯一的。三次样条插值函数的具体求解过程在下面的例子中给出了详细的说明。,例2.5.3 已知函数的函数值如下 在区间 上求三次样条插值函数，使它满足边界条件：解 先根据给定数据和边界条件算出，写出确定 的线性方程组。在本例中，给出的是第（1）种边界条件，确定 的线性方程组，形如（2.5.13）。由所给函数表知：于是由 的算式（2.5.9）知：,由第（1）边界条件下 与 的计算公式（2.5.14）知：故确定 与 的方程组为：（2.5.19）然后解所得方程组，得到 在各节点 上的值。在本例中，解（2.5.19）得：,最后将所得 代入（2.5.4），即得 在各子区间上的表达式。由（2.5.4）知，在 上的表达式为：在本例中，将 代入，整理后得 同理可得：,故所求三次样条插值函数为：第（1）边界条件下计算三次样条插值函数 在 处函数值的程序框图如图2.5.4,上述求三次样条插值函数的方法，其基本思路和特点是:先利用一阶导数 在内节点 上的连续性以及边界条件，列出确定二阶导数 的线性方程组（在力学上称为三弯矩方程组），并由此出，然后用 来表达。通过别的途径也可求三次样条插值函数。例如，可以先利用二阶导数在内节点上的连续性以及边界条件，列出确定一阶导数 的线性方程组（在力学上称为三转角方程组），并由此出，然后用 来表达。在有些情况下，这种表达方法与前者相比较，使用起来更方便。,2.6 求函数近似表达式的拟合法,在科学实验和生产实践中，经常要从一组实验数据 出发，寻求函数 的一个近似表达式（称为经验公式）。从几上，就是希望根据给出的 个点，求曲线 的一条近似曲线。因此，这是一个曲线拟合的问题。,多项式插值虽然在一定程度上解决了由函数表求函数的近似表达式问题，但用它来解决这里提出的问题，有明显缺陷。首先，实验提供的数据通常带有测试误差。如要求近似曲线 严格地通过所给的每个数据点，就会使曲线保持原有的测试误差。当个别数据的误差较大时，插值效果显然是不理想的。,其次，由实验提供的数据往往较多（即 较大），用插值法得到的近似表达式，明显地缺乏实用价值。因此，怎样从给定的一组数据出发，在某个函数类 中寻求一个“最好”的函数 来拟合这组数据，是一个值得讨论的问题。随着拟合效果“好”、“坏”标准的不同，解决此类问题的方法也不同。这里介绍一种最常用的曲线拟合方法，即最小二乘法。,2.6.1 曲线拟合的最小二乘法 如前所述，在一般情况下，我们不能要求近似曲线 严格地通过所有数据点，亦不能要求所有拟合曲线函数在 处的偏差（亦称残差）。都严格地趋于零。但是，为了使近似曲线尽量反映所给数据点的变化趋势，要求都较小还是必要的。达到这一目标的途径很多，常见的有：,(1)选取，使偏差绝对值之和最小，即(2)选取，使偏差最大，绝对值最小，即(3)选取，使偏差平方和最小，即,为了方便计算、分析与应用，我们较多地根据“偏差平方和最小”的原则（称为最小二乘原则）来选取拟合曲线 按最小二乘原则选择拟合曲线的方法，称为最小二乘法。本章要着重讨论的线性最小二乘问题，其基本提法是：对于给定数据表,要求在某个函数类（其中）中寻求一个函数 使 满足条件 式中 是函数类 中任一函数。,满足关系式（2.6.5）的函数，称为上述最小二乘问题的最小二乘解。由上可知，用最小二乘法解决实际问题包含两个基本环节：先根据 所给数据点的变化趋势与问题的实际背景确定函数类,即确定 所具有的形式；然后按最小二乘法原则（2.6.3）求取最小二乘解，即确定其系数。,由最小二乘解（2.6.4）应满足条件（2.6.5）知，点 是多元函数 的极小点，从而 满足方程组：即,亦即若对任意的函数 和，引入记号则上述方程组可以表示成,写成矩阵形式即（2.6.7）方程组(2.6.7)称为法方程组。,事实上，最小二乘法的法方程可以用下面的方法形成。在中，当 时，令，即得方程组：将其写成矩阵形式：,令则方程组可写为：将方程两边同时乘以，则可以得到 这就是最小二乘法的法方程（2.6.7）。当 线性无关时，可以证明它有唯一解，并且相应的函数（2.6.4）就是满足条件（2.6.5）的最小二乘解。,当,并且相应的函数(2.6.4)就是满足条件(2.6.5)的最小二乘解。综上分析可得定理2.6.1 对任意给定的一组实验数据(其中 互异),在函数类(线性无关)中,存在唯一的函数 使得关系式(2.6.5)成立，并且其系数,可以通过解方程组(2.6.7)得到。作为曲线拟合的一种常用的情况，若讨论的是代数多项式拟合，即取 则由(2.6.6)知：,故相应的法方程组为：（2.6.8）,下面通过两个具体的例子来说明用最小二乘法解决实际的问题的具体步骤与某些技巧。例2.6.1 某种铝合金的含铝量为，其熔解温度为，由实验测得 与 的数据如表2.6.1左边三列。使用最小二乘法建立与之间的经验公式。解 根据前面的讨论，解决问题的过程如下：（1）将表中给出的数据点 描绘在坐标纸上，如图2.6.1所示,图 2.6.1（2）确定拟合曲线的形式。由图2.6.1可以看出，六个点位于一条直线的附近，故可以选用线性函数（直线）来拟合这组实验数据，即令（2.6.9）其中 为待定常数。,（3）建立法方程组。由于问题归结为一次多项式拟合问题，故由(2.6.8)知，相应的法方程组形如 经过计算（表2.6.1）即得确定待定系数 的法方程组：,（2.6.10）,（4）解法方程组（2.6.10）得代入(2.6.9)即得经验公式：表2.6.1,（2.6.11）,所得经验公式能否较好地反映客观规律,还需通过实践来检验。由(2.6.11)式算出的函数值(称为拟合值)与实际值有一定的偏差。,表2.6.2,由表2.6.2可以看出，偏差的平方和 其平方根（称为均方误差）在一定程度上反映了所得经验公式的好坏。同时，由表2.6.2还可以看出，最大偏差为 如果认为这样的误差都允许的话，就可以用经验公式(2.6.11)来计算含铝量在 之间的溶解度。否则，就要用改变函数类型或者增加实验数据等方法来建立新的经验公式。,例2.6.2 在某化学反应里，测得生成物的浓度 与时间 的数据表见表2.6.3，是用最小二乘法建立 与 的经验公式。解 将已知数据点 描述在坐标纸上，见图2.6.2。由图2.6.2及问题的物理背景可以看出，拟合曲线 应具下列特点：,表2.6.3,（1）曲线随着 的增加而上升，但上升速度由快到慢。（2）当 时，反应还未开始，即;当 时，趋于某一常数。故曲线应通过原点(或者当 时以原点为极限点)，且有一条水平渐近线。具有上述特点的曲线很多。选用不同的数学模型，可以获得不同的拟合曲线与经验公式。下面提供两种方案:,图2.6.2,方案1 设想 是双曲线型的，并且具有形式 此时，若直接按最小二乘法原则去确定参数 和,则问题归结为求二元函数 的极小点，这将导致求解非线性方程组：,（2.6.13）,（2.6.12）,给计算带来了麻烦。可通过变量替换来将它转化为关于待定参数的线性形函数。为此，将(2.6.12)改写成 于是，若引入新变量则(2.6.12)式就是,同时，由题中所给数据表2.6.3可以算出新的数据表2.6.4这样，问题就归结为：根据数据表2.6.4，求形如 的最小二乘解 参照例2.6.1的做法，解方程组,表2.6.4,方案2 设想 具有指数形式（2.6.15）为了求参数 和 时，避免求解一个非线形方程组，对上式两边取对数 此时，若引入新变量 并记，则上式就是,又由表2.6.3可算出新的数据表2.6.5。于是将问题归为：根据数据表2.6.5，求形如 的最小二乘解。参照方案1，写出相应的法方程组并解之，即得,表2.6.5,于是 故得另一个经验公式（2.6.16）将两个不同的经验公式(2.6.14)和(2.6.16)进行比较（表2.6.6），从均方误差与最大偏差这两个不同角度看，后者均优于前者。因此，在解决实际问题时，常常要经过反复分析，多次选择，计算与比较，才能获得好的数学模型。,下面以常用的多项式拟合为例，说明最小二乘法在电子计算机上实现的步骤。设有一组实验数据，今要用最小二乘法求一 次多项式曲线来拟合这组数据。显然，求 的实质就是要确定其系数。,表 2.6.6,由前面讨论可知，问题可归结为建立和求解法方程组(2.6.8)。为了便于编制程序并减少工作量，引入矩阵 则法方程组(2.6.8)的系数矩阵（用 表示）和右端向量（用 表示）分别为:相应的程序框图如图2.6.3所示,其中中间矩阵 的生成方法见图2.6.4。,2.6.2 加权最小二乘法 在实际问题中测得的所有实验数据，并不是总是等精度、等地位的。显然，对于精度较高或地位较重要（这应根据具体情况来判定）的那些数据，应当给予较大的权。在这种情况下，求给定数据的拟合曲线，就要采用加权最小二乘法。用加权最小二乘法进行曲线拟合的要求与原则：对于给定的一组实验数据,要求在某个函数类中，寻求一个函数使其中,为函数类 中任一函数；是一列正数，称为权，它的大小反映了数据地位的强弱。显然，求 的问题可归结为求多元函数：的极小点。采用最小二乘解的求法，仍可得法方程组（2.6.7），但其中,作为特例，如果选用的拟合曲线为 那么相应的法方程组为,（2.6.17）,解 因为拟合曲线为一次多项式曲线（直线）故相应的法方程组如（2.6.17）。将表中各已知数据代入，即得方程组 解之得,2.6.3 利用正交函数作最小二乘法拟合 在前几节，虽然从原则上解决了最小二乘意义下的曲线拟合问题，但在实际计算中，由于当 较大，例如，法方程组往往是病态的，因而给求解工作带来了一定困难。近年来，产生了许多解决这一困难的新方法。本节将简要介绍利用正交函数作最小二乘拟合的基本原理，以及利用正交多项式拟合的一种行之有效的方法。,对于 和权，若一组函数 满足条件：则称 是关于点集 带权 的正交函数族。特别，当 都是多项式时，就称 是关于点集 带权 的一组正交多项式。,（2.6.18）,如果在提到正交函数或正交多项式时，没有提到权，就意味着权都是1。若所考虑的函数类 中的基函数是关于给定点集 和权 的正交函数族，则由条件(2.6.18)知，在法方程组(2.6.7)的系数矩阵中，非对角线上元素,此时法方程组简化为：只要由此解出 就可得到最小二乘法解：,（2.6.19）,（2.6.20）,由条件(2.6.18)，知 故易解方程组（2.6.19)，且有：这样就避免了求解一个病态方程组。,（2.6.21）,若函数类 的基函数 是关于给定点集 和权 的正交函数族，则可以直接由(2.6.21)式算出待定参数，进而写出最小二乘解(2.6.20)。因此，问题归结为对给定的函数类，寻求一组由正交函数族组成的基函数的问题。构造正交函数组的方法很多。下面以多项式为例，介绍一种具体的方法。这种做法是以下述定理为基础的：,构造的函数族 是关于点集 和权 的一组正交多项式,且 是 次多项式，其最高次项 的系数为1,（2.6.24）,例2.6.4 已知一组实验数据如表2.6.8，试用最小二乘法求一条二次拟合曲线。解 采用边构造正交多项式 边求最小二乘解系数 的顺序来求拟合曲线。由公式(2.6.22)及(2.6.21)立即可得：,表2.6.8,由公式(2.6.23)、(2.6.22)及（2.6.21)依次可得：,，,由公式(2.6.23),(2.6.22)、(2.6.24)及(2.6.21)依次可得：故所求拟合曲线为：,2.7 城市供水量预测的简单方法,2.7.1 供水量增长率估计与数值微分 作为多项式插值的应用，本节介绍两种求函数导数的近似值的方法。2.7.2 利用插值多项式求导数 若函数 在节点 处的函数值已知，就可作 的 次插值多项式，并用 近似代替，即(2.7.1)由于 是多项式，容易求导数，故对应于 的每一个插值多项式，就易建立一个数值微分公式 这样建立起来的数值微分公式，统称为插值型微分公式。,必须注意，即使 与 的近似程度非常好，导数 与 在某些点上的差别仍旧可能很大，因而，在应用数值微分公式时，要重视对误差的分析。由插值余项公式（2.3.9）知（2.7.2）由于式中 是 未知函数，故 时，无法利用上式误差 作出估计。但是，如果我们限定求某个节点 处的导数值，那么（2.7.2）右端第二项之值应为零，此时有,若将它写成带余项的数值微分公式，即 其中 在 之间。该式右端由两部分，即导数的近似值和相应的截断误差组成。由（2.7.3）,当 时，插值节点为，记 得带余式的两点公式 前一公式的实质是用 在 处的向前差商（分子是向前差分的差商）作为 的近似值，后一公式则是用 在 处的向后差商作为 的近似值。,（2.7.3）,(2.7.4),当 且节点为 时,由(2.7.3)可得带余项的三点公式：中间一个公式的实质是用 在 处的中心差商作为 的近似值，它与前后两公式相比较，其优越性是显然的。用插值多项式 作为 的近似函数，还可用来建立高阶的的数值微分公式。,(2.7.5),例2.7.1 带余式的二阶三点公式2.7.3 利用三次样条插值函数求导 由三次样条插值函数知，对于给定函数表,(2.7.6),和适当的边界条件，可以写出三次样条插值公式，并用 近似代替,即由于 是一个分段三次多项式，在各子区间 上容易求出导数，故可建立数值微分公式 利用函数在节点上的函数值和边界条件,(2.7.7),(2.7.8),构造三次样条插值公式，并用它来计算 和 在下列点 处的近似值。计算结果如表2.7.1。,表2.7.1,由表2.7.1可以看出，利用三次样条插值函数 及其导数来逼近被插值函数 及其导数，其效果是相当好的。,2.7.4 城市供水量预测 现在回过头来解决本章2.1.1提出的城市供水量预测问题。用数值方法进行预测时，一般是根据数据的散点图特征采用插值或拟合来实现。2.7.4.1用插值方法预测2007年1月份城市的用水量 预测2007年1月份城市的用水量有两种办法：一是先预测1月份每天的日用水量，求和后即得到1月用水总量；二是直接用2000-2006每年1月份的总用水量来预测。,用06年每天的日用水量预测 以预测07年1月1日为例，由于数据量过大，07年1月的用水量与前一年相关性大，而与2000-2005年用水量的相关性不大。故仅用06年全部365天的日用水量作为插值节点，其中，表示第i天，即为第i天的日用水量（万吨/日）。其散点图如图2.7.1所示。设插值多项式为，则07年1月1日的日用水量的预测值即为 在 处的函数值。,图2.7.1 06年城市日用水量散点图,下面分别用拉格朗日与牛顿插值方法得到365次插值多项式，并计算得 分别为1.8781e+109、-9.3842e+167。这些结果的误差大得让人无法接受。说明在大区间上做高次插值不可行，实际中，当节点数较多时，一般用分段低次插值。这里用三次样条插值法求解，其样条曲线如图2.7.2所示。对应的插值结果为133.6641。,图2.7.2 三次样条曲线,但若用这个三次样条函数来估计，即该月后30天的日用水量，结果为98.1477、33.0292、-4.8837，误差逐渐增大。且从第四天开始的估计值就开始为负值-70.3188。结果说明，在用插值多项式进行外插时，当计算点远离插值区间时，误差增大。因此，外插时一般要求计算点靠近插值区间。综上所述，仅用06年一年数据预测07年1月份每天的日用水量时，即便使用低次插值，误差也会很大。因此用365个数据点做插值不合理。,2.用2000-2006年1月份每天的日用水量预测 采用2000-2006年1月每天的日用水量预测07年1月对应天的日用水量。同样，以预测07年1月1日为例。以2000-2006每年1月1日的日用水量作为插值节点，则 即为所要求的07年1月1日用水量的估计值。同理，其它30天的日用水量也对应地求出。此时，待插节点只有1个，且与插值区间距离为1。因插值节点只有7个，我们采用拉格郎日、牛顿和三次样条插值函数进行计算，三种方法计算的结果基本相同，求得07年1月每天日用水量分别为143.6962，139.