数值计算方法理论与典型例题选讲是为理工科大学本科课程“数值分析”和“计算方法”编写的教材与课外自学指导两用书,主要内容包括引言、插值法、线性方程组的直接解法与迭代法、方程求根、数据拟合与函数逼近、数值积分与数值微分、常微分方程数值解法、矩阵特征值与特征向量的计算。此外,为了兼顾学生能力的培养和考试技能的提高,并帮助其合理掌握学习**,数值计算方法理论与典型例题选讲附录包括上机实习、理工科专业“计算方法”模拟题6套、数学专业“数值分析”模拟题2套和数学专业硕士研究生入学考试“数值分析”模拟题2套,并附有答案。
数值计算方法理论与典型例题选讲可作为理工科大学数学及相关专业本科和研究生“(数值)计算方法”课程的教材,也可供从事科学与工程计算的科技工作者参考。 数值计算方法理论与典型例题选讲_雷金贵,蒋勇,陈文兵_科学出版社_

第1 章引言
作为数值计算工作的一部分, 误差估计占有重要地位. 如果数值计算得到的结果没
有进行相应的误差估计, 那么, 一般来说, 该结果就不能保证是十分可靠或可用的, 因此
一定要尽可能地进行误差估计这项工作. 那么, 什么是误差? 什么是误差估计? 怎么进行
误差估计呢? 本章将对这些看似简单的问题做一个汇总, 以方便后面章节的叙述.
本章的**是误差、有效数字与机器数系的概念与相关知识, 以及数值计算陷阱的
预防措施. 在1.3 节介绍应当掌握的典型例题.
1.1 误差、有效数字与机器数系
利用计算机进行数值计算可以按照图1.1 所示的流程进行操作.
上述流程图是用计算机解决实际问题的一般步骤, 虽然不同的书中有不同的画法, 但
本质上区别不大. 图1.1 表明, 数值计算的**步是建立数学模型, 该步骤在一般情况下
可以省略, 因为很多问题都早已建立了较为成熟的数学模型; 第二步是算法设计, 其实就
某些具体问题来说, 该步骤也不需要每个细节都由自己完成, 科技工作者只有学会\站在
巨人的肩膀上", 才能事半功倍, 但是作为初学者, 应该关注如何设计算法, 毕竟这是我们
学习这门课程的根本目的; 第三步、第四步要求将算法在计算机上实现, 严格来说这也是
数值计算的重要环节, 但本书的**不在于此, 只是在有些地方稍作论述; *后一项是结
果分析, 即分析结果的有效性和正确性, 一般来说是指误差估计以及算法的稳定性和收敛
速度的分析等.
1.1.1 误差的来源与概念
根据图1.1, 可以确定如下误差来源, 即模型误差、观测误差、截断误差和舍入误差.
那么到底什么是模型误差? 例如要研究大气运动问题, 如果研究人员研究的是有黏
性、可压缩、斜压、干气体的深厚系统, 则需要建立如下Navier-Stokes 方程组(简称N-S
方程组)
上述N-S 方程组就是天气预报问题的数学模型, 该模型考虑了气体的动力学特征, 如
压强p、速度矢量v 和密度? 以及热力学特征如温度T 、热源(或汇)Q, 还考虑了地球引
力g、摩擦力F 以及地转偏向力?2- ￡v 等因素的影响, 但没有考虑大气化学特征和大
气电场等诸多因素对气体状态的影响. 因此, 即使能够求出上述方程组的**解, 这个所
谓的\**解" 也注定不是真实大气的真实状态. 这种建立数学模型时由于舍弃化学特性
和电场等影响因素而产��的误差称为模型误差. 用数学方法解决任何实际问题都会产生
模型误差, 也就是说, 模型误差是不可避免的.
观测误差就比较好理解, 例如, 今天南京的实际气温是25℃, 但是测量人员用仪器测
量得到的温度测量值是25.3℃, 这0.3℃的差值就是今天南京大气温度的测量误差. 测量
仪器总是会产生误差, 况且观测员在读数据的时候也会产生误差, 因此测量误差也是不可
避免的. 由仪器和人为测量产生的误差通称为测量误差或观测误差.
一般来说, 舍入误差产生在存储环节, 比如用二进制表示十进制数0.01, 然后再输出
十进制表示后的结果时, 将会产生截断误差, 这种误差就是舍入误差. 事实上, 如果我们
用20 个比特存储十进制数的小数部分, 然后再转化成十进制输出, 计算机输出的结果将
是0.009 999 275 而不是0.01. 因此原来的十进制数0.01 与计算机输出结果0.009 999 275
之间的差值就是舍入误差. 舍入过程可参见图1.2.
同理, 用计算机存储数值1=3 时, 在内存中的实际值是小数, 例如, 0.333 333 33. 由于
计算机总是用有限位二进制单元来存储整数和实数, 因而舍入误差也是不可避免的.
截断误差可以这样理解, 例如N-S 方程组由于无法直接求得它的**解, 只能通过
数值方法来求解, 这种数值解一般是通过一种称为模式或差分格式的算法程序求得, 如当
前气象学界使用的WRF 模式, 则模式的\**解" 与N-S 方程组的\**解" 之间的误
差称为截断误差. 有关截断误差的讨论, 只有大家学习了微分方程数值解或数值积分、插
值法等相关章节后才能有一个比较深刻的理解, 读者在此不必深究.
总之, 在用数值方法求解实际问题的时候, 各种误差都不可避免, 并且*终的\数值
解" 与物理问题\**解" 的误差总是以上面提到的各种误差的总和形式出现. 况且如果
研究的问题很复杂, 则很难得到物理问题的\**解", 所以误差估计会较难进行. 所以,
青年学生在工作之前首先应该学会工作方法, 否则将来的工作将难以开展.
可能也有人认为既然数值计算很困难, 那为什么不能求出物理问题的**解呢?
对于这个问题, 可以用上述天气预报的例子来说明. 理论上, 现在老百姓每天从电视
台、网络或者报纸中获取的天气预报信息基本上是**气象台用数值预报模式计算出来
的N-S 方程组的数值解. 由于N-S 方程组是一个非线性的偏微分方程组, 按照人类目前
的科学技术水平并不能求出其显式的**解, 只能求其数值解. 每天, 气象台的科学家们
就是用全国各个气象观测台站的地面常规观测资料、卫星遥感资料以及雷达探空资料等
诸多数据构造的大气初始状态为初始值, 用数值预报模式计算出数值结果; 该结果经过气
象台站的软件加工处理后就是我们得到的天气预报. 天气预报有时很准确, 有时不准确,
这说明对预报结果的误差估计有时很小, 有时很大.
学习数值计算方法课程, 首先要掌握两个定义:**误差和相对误差. 所谓**误
差e, 就是指两个数的差
e = x¤ ? x (1.1)
其中, x¤ 是某个数学量或物理量(可以是函数的自变量或因变量) 的准确值, x 是x¤ 的测
量值或近似值.
两个数的相对误差定义为
er = x¤ ? x
x¤
(1.2)
由于准确值x¤ 通常无法获得, 因此常用公式
x¤ ? x
x
代替
x¤ ? x
x¤
来求近似相对误差, 有
时也称
x¤ ? x
x
是相对误差, 实际上当x 非常接近x¤ 时,
x¤ ? x
x ?
x¤ ? x
x¤ ? o 3(x¤ ? x)2′
因此当
x¤ ? x
x¤
无法计算时, 并不需要估计表达式
x¤ ? x
x
-
x¤ ? x
x¤
的大小, 我们直接就
用
x¤ ? x
x
计算相对误差er 的大小.
当存在常数±, 使得
jej = jx¤ ? xj 6 ±
时, 则称± 为**误差限; 而当存在常数° 使得
jerj =ˉˉ
x¤ ? x
x¤ ˉˉ 6 °
时, 称° 为相对误差限. **误差限可以用来估计**值, 即
x ? ± 6 x¤ 6 x + ±
通常, 要衡量一个近似值是否与其对应的**值足够接近, 不能简单看**误差是否
小, 而是要看相对误差的大小. 例如, 在测量小行星的大小时, 行星半径的**误差是几
公里, 则可以认为测量已经相当**, 但如果测量的是一个普通建筑物的高度, 相差即便
只有几米, 误差也已经很大了. 所以通常情况下, 不能仅靠**误差就断定一个近似值的
好与坏.
1.1.2 误差的传播
下面以二元函数y = f(x1; x2) 为例, 解释自变量的误差是怎样传递给因变量的.
设x1; x2 与x¤1; x¤2 分别是二元函数y = f(x1; x2) 自变量的近似值与**值, y 与y¤
是因变量的近似值与**值, 则函数值的**误差为
e(y) = y¤ ? y = f(x¤1 ; x¤2) ? f(x1; x2)
将f(x¤1 ; x¤2 ) 在(x1; x2) 处Taylor 展开到一阶导数项, 略去高阶项, 得
e(y) = y¤ ? y ? f(x1; x2) + @f(x1; x2)
@x1
(x¤1 ? x1) + @f(x1; x2)
@x2
(x¤2 ? x2) ? f(x1; x2)
= @f(x1; x2)
@x1
e(x1) + @f(x1; x2)
@x2
e(x2) (1.3)
上式说明自变量xi 的**误差e(xi) 对e(y) 的贡献是
@f(x1; x2)
@xi
e(xi), 其中, 偏导数
@f(x1; x2)
@xi
是这个方向的放大系数, i = 1; 2.
二元函数的相对误差为
er(y) = e(y)
y ?
@f(x1; x2)
@x1
e(x1)
y
+ @f(x1; x2)
@x2
e(x2)
y
= @f(x1; x2)
@x1
x1
y
er(x1) + @f(x1; x2)
@x2
x2
y
er(x2) (1.4)
利用式(1.3) 和式(1.4) 可得到一些简单二元函数, 如两个自变量x1 与x2 的和、差、
积、商的**误差和相对误差的估计式
e(x1 + x2) ? e1 + e2
e(x1 ? x2) ? e1 ? e2
e(x1x2) ? x2e1 + x1e2
e(x1=x2) ? e1=x2 ? e2x1=x2
2; x2 6= 0
er(x1 + x2) ?
x1
x1 + x2
e1r + x2
x1 + x2
e2r; x1 + x2 6= 0
er(x1 ? x2) ?
x1
x1 ? x2
e1r ?
x2
x1 ? x2
e2r; x1 ? x2 6= 0
er(x1x2) ? e1r + e2r
er(x1=x2) ? e1r ? e2r; x2 6= 0
1.1.3 有效数字
若实数x¤ 的非零近似值为x, x 的十进制规格化浮点数的标准形式为
x = §0:x1x2 … xn … xp ￠ 10m (1.5)
其中, x1 6= 0, 如果有jx¤ ? xj = je(x)j <
1
2 ￡ 10m?n, 则称x有n位有效数字.
有效数字的定义, 也可以解释为:从左边**位非零数字开始数, 到第n 位数字截止,
如果**误差小于第n 位数字的数值的0:5 倍, 则该近似数的有效数字位数为n.
关于有效数字, 需要记住的是:① 用四舍五入得到的数都是有效数字; ② 有效数字
越多, 截断误差越小, 计算结果越**.
1.1.4 机器数系
计算机内部能够处理的数据除了整型数据和字符型数据之外, *常用的是浮点数, 也
就是人们常用的实数类型.
实数类型在C 语言中称为浮点数(用f loat 定义单精度或用double 定义双精度), 在
Fortran 语言中称为实型数据(用关键字real*m 来定义实型数据, 其中, m 为一数字, 表示
实型数据的种数, 详细解释请参见Fortran 语言相关文献).
总体来看, 浮点数是计算机表示范围*广、*为常用的一种数据类型, 表示为式(1.5)
形式的数称为规格化浮点数, 式(1.5) 只是十进制形式浮点数的表示形式. 计算机能够处
理的数不仅仅是十进制的, 它有多种形式, *常见的是二进制、八进制、十六进制三种,
这与计算机用二进制存储数据有关. 1985 年, 美国电气与电子工程师学会IEEE(Institute
of Electrical and Electronics Engineers) 颁布了名为二进制浮点数标准的报告(编号:754-
1985; 2008 年, 该标准更新为IEEE 754-2008). 该报告为二进制和十进制浮点数、数据
交换格式、舍入误差算法和异常处理操作等提供了一系列的标准. 标准中还规定了单精
度、双精度以及扩展精度数据类型的格式, 这些格式被所有的计算机厂商广泛遵守, 并
做成硬件模块集成到计算机设备中. 在2008 年的新标准中, 实数用64 位二进制数表
示, 其中, 从左边起**个比特是符号位(sign indicator), 接下来的11 个比特是指数部分
(characteristic), 剩余的52 个比特表示小数部分(mantissa)[1]:
一般情况下, 一个ˉ 进制的规格化浮点数可以表示为
x = §0:a1a2 … at ￠ ˉp (1.6)
其中, a1 6= 0, 0 6 aj 6 ˉ ? 1; j = 1; 2; … ; t;L 6 p 6 U; aj 和p 都是整数. 如果把计算机
中所有浮点数的集合加上\机器零" 记为F, 则F 被如下4 个参数所描述, 分别是基数
ˉ、字长t、阶码范围[L;U]. 集合F 称为机器数系. 应该说, 机器数系F 是一个分布不均
匀的、有限个离散数据的集合, 阶码越小的地方表示的数越稠密, 阶码越大的地方表示的
数越稀疏(图1.3).
不难证明, 集合F 仅含有1 + 2(ˉ ? 1)ˉt?1(U ? L + 1) 个数.
图1.3 机器数系
A: **值*大负实数, D: **值*大正实数; B: **值*小负实数, C: **值*小正实数
图1.3 把机器数系能够表示的**值*大和**值*小的4 个数标在数轴上, 当一
个数的值大于数D 或者小于A 时, 则为上溢; 当一个非零浮点数的值介于数轴上的B 和
C 之间时, 则产生下溢, 该数被计算机当做数值0 来处理.
因此当一个实数x = §0: a1a2 … atat+1 ￠ ˉp(a1 6= 0) 被存入计算机后, 称之为机器数,
记为fl(x). 当前计算机有两种方式产生机器数fl(x), 一种是直接截取x 的前t 位数, 得
到fl(x) = §0: a1a2 … at ￠ ˉp, 称这种计算机为截断机; 还有一种计算机, 对数x 的第t+1
位数字四舍五入得到~at, 即得到fl(x) = §0: a1a2 … ~at ￠ ˉp, 称这种计算机为舍入机. 一般
来说, 用截断机处理数值计算所得结果的误差要大于等于舍入机的误差.
1.2 数值计算陷阱的防范措施
以计算机为工具进行数值计算, 具有速度快、精度高、适合重复计算的特点. 但是数
值计算有很多的特殊性, 比如会出现大数吃小数、数值溢出、自动截断等问题, 这些问题
为我们利用计算机进行数值计算设置了很多的陷阱, 必须加以防范. 下面就来详细讨论数
值计算中的这些陷阱.
1.2.1 注意防止大数吃小数
例1.1 编程计算1011 + 1. 实现上述计算的程序和结果如图1.4 和图1.5 所示.
图1.4 所示程序和图1.5(a) 及图1.5(b) 所示的程序实现两个相同数据的加法运算,
不同之处在于两个变量的数据类型不同, 因此输出的计算结果不同. 图1.5(b) 中的结果是

前言
第1章 引言
1.1 误差、有效数字与机器数系
1.2 数值计算陷阱的防范措施
1.3 典型例题分析
第2章 插值法
2.1 插值问题
2.2 Lagrange(拉格朗日)插值法
2.3 Newton插值多项式与差商
2.4 差分与等距节点插值
2.5 Hermite(埃尔米特)插值
2.6 分段插值法
2.7 三次样条插值函数
第3章 线性方程组的直接解法
3.1 问题提出
3.2 Gauss(高斯)消去法
3.3 追赶法
3.4 矩阵的三角分解
3.5 向量范数和矩阵范数
3.6 摄动理论与误差分析初步
第4章 解线性方程组的迭代法
4.1 Jacobi迭代法与Gae-seidel迭代法的构造
4.2 迭代法的收敛性
4.3 SOR(逐次超松弛)迭代法
第5章 方程求根
5.1 方程根的存在性、**性与二分法
5.2 迭代法的基本概念与收敛性
5.3 加速方法
5.4 Newton—Raphson迭代法
5.5 割线法
5.6 代数方程求根
第6章 数据拟合与函数逼近
6.1 矩阵的广义逆
6.2 方程组的*小二乘解
6.3 矩阵的正交分解与方程组的*小二乘解
6.4 正交多项式
6.5 数据拟合
6.6 函数逼近初步
第7章 数值积分与数值微分
7.1 数值积分的基本思想与代数精度
7.2 Newton-cotes(牛顿科茨)型求积公式
7.3 区间逐次二分法与Romberg算法
7.4 Gauss(高斯)型积分公式
7.5 数值微分简介
第8章 常微分方程数值解法
8.1 常微分方程初值问题
8.2 Runge-Kutta(龙格-库塔)法
8.3 单步法的收敛性和稳定性
8.4 线性多步法
8.5 常微分方程组与边值问题的数值解法
第9章 矩阵特征值与特征向量的计算
9.1 幂法与反幂法
9.2 Jacobi(雅**)方法
9.3 QR算法
附录1 上机实习
附录2 理工科专业期末考试模拟题
附录3 数学专业考试模拟题
参考文献

雷金贵和蒋勇等编著的《数值计算方法理论与典型例题选讲》利用计算机进行数值计算工作的核心是:①研究建立数学模型;②研究数值算法的构造方法,并探索算法的特点与规律,如计算速度、收敛性和稳定性等问题;③理论联系实际,研究如何编制相应的算法软件,以及如何降低软件的计算复杂度、时间复杂度,减少内存的占用,降低算法的结构复杂性,增加算法的稳定性等问题这些也是计算方法的研究内容。