计算方法第1章预备知识与误差分析

计算方法北京工业大学应用数理学院杨中华第一章预备知识与误差分析1.1绪论

随着计算机的问世，出现了为计算机提供算法的学科—计算数学，计算方法即是计算数学中的重要课程。

概括起来计算机所涉及的算法可分为两类：数值算法和非数值算法。非数值算法在“数据结构”课程讲授，而计算方法这门课程所讲授内容是最为常用的数值算法，更确切的说计算方法讲授我们在高等数学、线性代数、常微分方程课程中常见的一些数学问题的数值求解方法。计算方法也可称为数值分析，计算方法和数值分析略有区别，前者着重于各算法的实现，后者兼论算法的理论分析。

由于本课程的授课对象主要是计算机学院和软件学院的本科学生，在算法的理论基础方面不追究精巧的数学推导，力图在算法的编程实现方面给出一些作者的体会，在设计算法的各个步骤时，尽量给出便于改写成程序语言的计算公式。

本课程在讲授过程中还将给出主要算法的C程序，这些程序仅用于辅助算法过程的学习和数值试验，不一定适用用于商业软件开发，所以程序中没有过多的考虑数值计算中的异常情况，如果将这些程序用于商业软件开发，还需在程序中添加异常情况的保护机制，以增强程序的健壮性。绪论1.本课程主要内容

包括计算数学中经典的数学问题：误差分析、非线性方程的数值解法、线性方程组的数值解法、插值方法与数据拟合、数值积分与数值微分、矩阵特征值与特征向量、常微分方程的数值解法。绪论2.学习计算方法的重要性和必要性

对于任何数学问题的求解，最为理想的解决方案当然是能够给出其解的解析表达式，但是有时这是难以实现的，如下述代数方程的求解

当n>4时其解析形式的解是不存在。还有很多数学问题的解都无法用解析形式的表达式表示；另一方面，某些数学问题的数据并不是用解析表达式给出的而是用表格形式给出的离散数据，当然其解也没有解析表达式；再者某些数学问题的解即使可以用解析方式表达出来，按这类表达式计算，时间代价也是不可接受的。绪论例1.1考虑一个求解n阶线性方程组的问题：绪论由线性代数知识我们知道，当其系数矩阵对应的行列式不等于零时，即|A|≠0，该线性方程组有唯一一组解，根据克莱姆法则，有：我们考虑用克莱姆法则求解此问题的时间代价。

通常度量一个算法的时间复杂度是统计该算法中乘除法的运算次数，我们知道每计算一个n阶行列式，需要乘法次数为n!，克莱姆法则求解线性方程组需要计算n+1个行列式，所需总的乘法次数为：

绪论如果求一个20阶线性方程组，把n=20代入，得乘法次数为：若以运算速度为每秒一亿次(108)浮点数乘除法运算的计算机昼夜不停的连续计算，其耗时约为：

这个耗时数还不包括求解过程中的加减运算以及更耗时的读写内存数据操作所需要的时间。绪论

用Gauss消去法在同样速度的计算机上，每秒可以求解出32679个20阶的线性方程组，由此可以看出求解线性方程组用Gauss消去法非常有效，因此对于稍微大一点规模的线性方程组没有任何理由选择克莱姆法则解决此类问题。对程序员的忠告：千万不要以为计算机的速度不是问题，选择数学方法不当计算机的速度太慢将是一个很大的问题！

但是如果用Gauss消去法求解此线性方程组，其乘除法次数约仅为：

我们再看一个实例，从中可以发现，有时直接使用高等数学中给出的很简单明了的数学表达式进行计算并不一定能够得到我们预期的结果。

例1.2考虑导数的近似计算问题，根据导数的定义绪论再根据极限的定义可知，当Δx足够小时可以用上式右端的差商近似此导数，如果取分别用单精度和双精度进行此近似计算，所得结果列表如下：表1.1差商近似导数的单精度计算结果绪论变量差差商误差1.00e-0022.73191118e+0001.363e-0021.00e-0042.71953058e+0001.249e-0031.00e-0062.54005599e+0001.782e-0011.00e-008-8.25483990e+0001.097e+0011.00e-010-8.25484009e+0028.282e+0021.00e-012-8.25483984e+0048.255e+0041.00e-014-8.25484000e+0068.255e+006精确值为：表1.2差商近似导数的双精度计算结果绪论变量差差商误差1.00e-0022.73191866e+0001.364e-0021.00e-0042.71841775e+0001.359e-0041.00e-0062.71828319e+0001.359e-0061.00e-0082.71828181e+0002.107e-0081.00e-0102.71828193e+0001.014e-0071.00e-0122.71856951e+0002.877e-0041.00e-0142.69448092e+0002.380e-002精确值为：对程序员的忠告：千万不要以为按照一个优美的数学公式和数学算法就可以得到正确的计算结果，误差的传播与积累可能使最后的计算结果与真正的解相差万里！

从以上两个表格可以看出，用差商近似导数并非Δx愈小愈精确，但是按高等数学中导数的定义Δx愈小愈接近导数值

，实际计算结果似乎有悖于高等数学中所学的知识！

但是如果用本课程中的误差分析的理论则很容易解释此类现象为何发生，因此具有这些误差分析的知识对编程人员来说将是很重要的也是必须的。绪论1.2误差分析

所谓误差就是真实值或精确值与近似值之差。用计算机解决数学问题的过程中会产生多种误差，而引起误差的原因是多方面的，下面我们通过考察解决实际问题的过程来具体说明误差的来源。

1.误差的来源及误差类型

一般使用计算机解决实际问题须经过如下几个过程：实际问题数学模型数值方法程序设计计算结果

(1)模型误差、根据实际问题建立数学模型的过程中通常会忽略某些次要因素而对问题进行简化，由此产生的误差称为模型误差。例如自由落体运动规律的公式

是忽略了空气阻力的影响，因此产生的误差就是模型误差。

误差分析(1.6)

(2)参数误差、很多数学模型都含有一些参数，而有些参数往往又是观测得到的近似值，如此取得的近似参数与真实参数值之间的误差称为参数误差或观测误差。

在自由落体运动规律的公式中重力加速度通常取g=9.8，实际上参数g的值与落体所在位置的地球纬度和海拔高度有关，此近似参数将导致出现参数误差。

(3)截断误差、通常在选定数学模型及参数后还要选择数值算法，很多数值算法往往是将一个无限过程截断为有限过程，此类误差称为截断误差或方法误差。例如前述的差商近似导数的计算方法实质上是在如下展开式中“截断”上式中

的二阶以上无穷小项得到的计算公式，由此产生的误差就是截断误差。误差分析

(4)舍入误差、因为计算机表示浮点数的字长有限(通常4字节或8字节)，除极少量能够准确表示的数据外绝大部分是在超过界定的某位4舍5入所得的近似值，比如圆周率等存入计算机内存时均需要在某位进行4舍5入，此种误差称为舍入误差。

虽然舍入误差在单步运算时或许微不足道，但是在一个较为复杂且连贯的数值算法中，舍入误差可能会积累、传播，所以我们必须加以关注。误差分析

2.绝对误差

定义1.1

设x是准确值，x*是它的一个近似值，称

为x的绝对误差，简称误差，误差的绝对值的上限ε*

称为x*的绝对误差限，简称误差限。

例如在使用圆周率时通常取为π*=3.14，其绝对误差及误差限为

在工程技术或商品规格中通常用如下方式表示产品的误差限：误差分析(1.8)(1.7)

3.相对误差

度量误差的另一种方式是相对误差，例如从10±0.1和1000±1中我们肯定不会认为绝对误差限小的比绝对误差限大的更精确。

定义1.2设x为准确值,x*是它的一个近似值，称比值

为近似值x*的相对误差，相对误差的绝对值的上限记为

称为x*的相对误差限，即有

误差分析(1.9)(1.10)

对于x=10±0.1和1000±1，我们可以计算出它们的相对误差分别为：

从绝对误差看，后者误差限更大些，但从相对误差的角度看，后者更为精确。

在(1.9)式中，在可以确定准确值x的前提下也可以取为x作为分母。误差分析

4.有效数字

定义1.3设x为准确值，其近似值为

这里

为整数，且

则称近似值x*有p位有效数字。

例1.3圆周率

，如果取3.1415和3.14159作为其近似值，问分别有几位有效数字？

解：(1)取

则3.1415有4为有效数字，最后一位数字5是无效数字。误差分析(1.11) (2)取

则取3.14159有6为有效数字。

我们在近似计算中四舍五入保留小数点以后n位数字时有效数字的个数，如同样保留小数点后4位数字 1.414213… 近似值1.4142有5位有效数字 0.012353… 近似值0.0124有3位有效数字 28.08918… 近似值28.0892有6位有效数字

总之：对精确值进行四舍五入所得结果后所有非零数字均是有效数字。

误差分析

定理1.1设近似值x*具有(1.11)式的形式，如果有n位有效数字或是由其准确值4舍5入后得到的近似值，则x*的相对误差限为：

证明：如果

有n位有效数字，则其绝对误差限为

相对误差限为

如果x*是由其准确值4舍5入后得到的近似值，则此近似值必然有n位有效数字，从而结论也必然成立。误差分析

5.近似计算中减少误差的几个策略 1)避免两个相近的数相减

我们知道一个近似值其有效数字愈多愈为精确，如果两个相近的数相减势必减少其有效数字的个数，使得相对误差变大。

例1.4已知

x=1.5846,y=1.5839，求x-y的值。

解：x–y=1.5846–1.5839=0.0007，显然x，y均有5位有效数字，而

x–y仅有1位有效数字，使得其有效数字大大的减少。误差分析

例1.5已知

x=18.496，y=18.493，取4位有效数字，计算x-y的近似值，并估计其相对误差。

解：先取4为有效数字，x*=18.50，y*=18.49，再计算其差

x*–y*=18.50–18.49=0.01

实际上按准确值计算 x–y=18.496–18.493=0.003

相对误差为：误差分析

在编程序时，可采取以下措施减少误差尽量使用双精度定义变量当两个接近的量相减时，最好做等价变换避免损失有效数字，如当x的值较大时令：当x与y非常接近时令：误差分析

2)避免绝对值太小的数作为分母

如

时，计算比值

会使绝对误差变大，通常也做等价变换，如

时：

3)防止大数“吃掉”小数

例1.6计算定积分

，其中N是一个非常大的正数。

解：根据Newton-Leibniz公式得到误差分析

但是如果按此式进行计算，由于N+1时1可能被N吃掉，所以实际计算时可能是N+1=N，因此计算结果可能是

但是实际上

结果严重失真，正确的计算方法应该是：误差分析

另一方面，根据积分中值定理有

我们对于

编程实际计算，结果如下：

因为

，断定后两个结果的精度优于第一个计算结果。误差分析1.3向量和矩阵的范数

如同误差分析是计算方法的预备知识一样，有关范数的概念也是学习计算方法所必须的基础知识。

通常用绝对值来度量实数的大小，用模来度量复数的大小，用范数来度量向量、矩阵的大小。

范数是论述、证明某些数值方法的收敛性和稳定性所必须的工具，本节将简单介绍范数的有关概念。 1.向量的范数

定义1.4如果定义在

上的一个实值函数

，对于任意的

都有 1)非负性：

且仅当

时等号成立 2)齐次性： 3)三角不等式：

则称该实值函数

为

上的一个向量范数。

满足定义1.4中三个条件的实值函数均是向量范数，但最常用的向量范数有： 1-范数 2-范数 ∞-范数向量和矩阵的范数(1.12)(1.13)(1.14)

对于2-范数与内积的关系，有著名的Cauchy-Schwarz(柯西—许瓦尔兹)不等式：

我们可以证明这三种向量范数均满足定义1.4中的三个条件，仅以2-范数为例加以证明，另两个留作习题。

证明：条件1)非负性是显然的，对于条件2)的齐次性

对于条件3)向量和矩阵的范数(1.16)

以上三种向量范数可以统一于p-范数定义之下，

p-范数定义如下：

当p=1、2时(1.17)式与(1.12)、(1.13)显然是一致的，对于p=∞情况，考虑下不等式向量和矩阵的范数

即有

对此式取p→∞时的极限，则有(1.17)例1.7求向量

前述的三种向量范数。

解：

如前所述向量范数可以度量向量的“大小”，而范数又不是唯一的，虽然对应同一个向量不同的范数值是不一样的，但是我们可以证明任意两种向量范数是等价的，首先给出向量范数等价性的定义。向量和矩阵的范数

定义1.5设有两种范数

，如果存在常数

使

成立，则称两种范数是等价的。

按此等价性的定义，常用的三种范数是等价的。事实上很容易证明以下三式成立

三种范数关系的综合不等式为：向量和矩阵的范数(1.18)

2.矩阵的范数

定义1.6如果定义在

上的一个实值函数

，对于任意的

都有 1)非负性：

且仅当

时等号成立 2)齐次性： 3)三角不等式：

4)相容性：

则称该实值函数

为

上的一个矩阵范数。

一个最常用而且便于计算的矩阵范数为Frobenius(弗罗贝纽斯)范数，定义如下：

该矩阵范数也简称为矩阵的F-范数，如果把矩阵按相邻行的首尾连接构成向量，则F-范数正是此向量的2-范数。向量和矩阵的范数(1.19)

以下讨论矩阵范数与向量范数之间的关系与相容性问题。

定义1.7对于给定

上的向量范数

和

上的矩阵范数

如果有

则称矩阵范数

与向量范数

是相容的。

我们也可以按相容性条件定义从属于向量范数的矩阵范数，即对于指定的向量范数

来定义矩阵的范数：

如此定义的矩阵范数

称为向量范数的从属范数，可以证明此从属范数满足矩阵范数定义的4个条件。向量和矩阵的范数(1.20)

从属于向量范数的矩阵1-范数、2-范数、

范数的表达式：

其中

是矩阵

的最大特征值。

对于以上3种从属范数我们均可以给出证明，此处我们仅就(1.21)式的1-范数给出证明向量和矩阵的范数(1.21)(1.22)(1.23)向量和矩阵的范数

证明：将矩阵A按列分块记为

假设第k列