数值分析与算法——读书笔记（一）

chapter1

数值计算导论

1.1概述

可以将数值计算的研究目标归纳为：寻找一个能迅速完成的（迭代）算法，同时评估计算结果的准确度。
数值计算研究的核心内容是数值算法的设计与分析。
IEEE1主办的《科学与工程中的计算》杂志评选出20世纪对科学和工程的发展和实践影响最大的10个算法，它们按时间顺序依次如下：
1. Monte Carlo方法
1946年，在洛斯阿拉莫斯科学实验室工作的John von Neumann，Stan Ulam和Nick Metropolis编制了Metropolis算法，也称为Monte Carlo方法。Metropolis算法旨在通过模仿随机过程，来得到具有难以控制的大量的自由度的数值问题和具有阶乘规模的组合问题的近似解法。数字计算机是确定性问题的计算的强有力工具，但是对于随机性（不确定性）问题如何当时并不知晓，Metropolis算法可以说是最早的用来生成随机数，解决不确定性问题的算法之一。
2. 线性规划的单纯形方法
1947年，兰德公司的Grorge Dantzig创造了线性规划的单纯形方法。就其广泛的应用而言，Dantzig算法一直是最成功的算法之一。线性规划对于那些要想在经济上站住脚，同时又有赖于是否具有在预算和其他约束条件下达到最优化的能力的工业界，有着决定性的影响(当然，工业中的“实际”问题往往是非线性的；使用线性规划有时候是由于估计的预算，从而简化了模型而促成的)。单纯形法是一种能达到最优解的精细的方法。尽管理论上讲其效果是指数衰减的，但在实践中该算法是高度有效的——它本身说明了有关计算的本质的一些有趣的事情。
3. Krylov子空间叠代法
1950年，来自美国国家标准局的数值分析研究所的Magnus Hestenes, Eduard Stiefel和Cornelius Lanczos开创了Krylov子空间叠代法的研制。这些算法处理看似简单的求解形为Ax=bAx=b的方程的问题。当然隐藏的困难在于A是一个巨型的n*n 矩阵，致使代数解x=b/Ax=b/A是不容易计算的(确实，矩阵的“相除”不是一个实际上有用的概念)。
叠代法——诸如求解形为Kx(k+1)=Kx(k)+b−Ax(k)Kx(k+1)=Kx(k)+b-Ax(k)的方程，其中K 是一个理想地“接近”A 的较为简单的矩阵——导致了Krylov子空间的研究。以俄罗斯数学家Nikolai Krylov命名的Krylov子空间由作用在初始“余量”向量r(0)=b−Ax(0)r(0)=b-Ax(0)上的矩阵幂张成的。当 A是对称矩阵时，Lanczos找到了一种生成这种子空间的正交基的极好的方法。对于对称正定的方程组，Hestenes 和Stiefel提出了称为共轭梯度法的甚至更妙的方法。过去的50年中，许多研究人员改进并扩展了这些算法。当前的一套方法包括非对称方程组的求解技巧，像字首缩拼词为GMRES和Bi-CGSTAB那样的算法。(GMRES和Bi-CGSTAB分别首次出现于1986和1992 SIAM journal on Scientific and Statistical computing(美国工业与应用数学学会的科学和统计计算杂志)。
4. 矩阵计算的分解方法
1951年，橡树岭国家实验室的A1ston Householder系统阐述了矩阵计算的分解方法。研究证明能把矩阵因子分解为三角、对角、正交和其他特殊形式的矩阵是极其有用的。这种分解方法使软件研究人员能生产出灵活有效的矩阵软件包。这也促进了数值线性代数中反复出现的大问题之一的舍入误差分析问题。 (1961年伦敦国家物理实验室的James Wilkinson基于把矩阵分解为下和上三角矩阵因子的积的LU分解，在美国计算机协会(ACM)的杂志上发表了一篇题为“矩阵逆的直接方法的误差分析”的重要文章。)
5. Fortran最优编译程序
1957年，John Backus在IBM领导一个小组研制Fortran最优编译程序。Fortran的创造可能是计算机编程历史上独一无二的最重要的事件：科学家(和其他人)终于可以无需依靠像地狱那样可怕的机器代码，就可告诉计算机他们想要做什么。虽然现代编译程序的标准并不过分??Fortran I只包含23,500条汇编语言指令??早期的编译程序仍然能完成令人吃惊的复杂计算。就像Backus本人在1998年在IEEE annals of the History of computing 发表的有关Fortran I，II, III的近代历史的文章中回忆道：编译程序“所产生的如此有效的代码，使得其输出令研究它的编程人员都感到吓了一跳。”
6. 矩阵本征值计算的QR算法
1959—61年，伦敦Ferranti Ltd.的J.G. F. Francis找到了一种称为QR算法的计算本征值的稳定的方法。本征值大概是和矩阵相连在—起的最重要的数了，而且计算它们可能是最需要技巧的。把—个方阵变换为一个“几乎是”上三角的矩阵??意即在紧挨着矩阵主对角线下面的一斜列上可能有非零元素??是相对容易的，但要想不产生大量的误差就把这些非零元素消去，就不是平凡的事了。QR 算法正好是能达到这一目的的方法，基于QR 分解， A可以写成正交矩阵Q 和一个三角矩阵R 的乘积，这种方法叠代地把 A=Q(k)R(k)A=Q(k)R(k) 变成 A(k+1)==Q(k)R(k)A(k+1)==Q(k)R(k) 就加速收敛到上三角矩阵而言多少有点不能指望。20世纪60年代中期QR 算法把一度难以对付的本征值问题变成了例行程序的计算。
7. 快速分类法
1962：伦敦Elliott Brothers, Ltd.的Tony Hoare提出了快速(按大小)分类法.把n个事物按数或字母的次序排列起来，在心智上是不会有什么触动的单调平凡的事。智力的挑战在于发明一种快速完成排序的方法。Hoare的算法利用了古老的分割开和控制的递归策略来解决问题：挑一个元素作为“主元”、把其余的元素分成“大的”和“小的”两堆(当和主元比较时)、再在每一堆中重复这一过程。尽管可能要做受到严厉责备的做完全部N(N-1)/2 次的比较(特别是，如果你把主元作为早已按大小分类好的表列的第一个元素的话！)，快速分类法运行的平均次数具有O(Nlog(N)) 的有效性，其优美的简洁性使之成为计算复杂性的著名的例子。
8. 快速Fourier变换
1965年，IBM的T. J. Watson研究中心的James Cooley以及普林斯顿大学和AT＆T贝尔实验室的John Tukey向公众透露了快速Fourier变换(方法)(FFT)。应用数学中意义最深远的算法，无疑是使信号处理实现突破性进展的FFT。其基本思想要追溯到Gauss(他需要计算小行星的轨道)，但是Cooley—Tukey的论文弄清楚了Fourier变换计算起来有多容易。就像快速分类法一样，FFT有赖于用分割开和控制的策略，把表面上令人讨厌的O(N*N) 降到令人欢乐的O(Nlog(N)) 。但是不像快速分类法，其执行(初一看)是非直观的而且不那么直接。其本身就给计算机科学一种推动力去研究计算问题和算法的固有复杂性。
9. 整数关系侦查算法
1977年，BrighamYoung大学的Helaman Ferguson 和Rodney Forcade提出了整数关系侦查算法。这是一个古老的问题：给定—组实数，例如说x(1),x(2),...,x(n)x(1),x(2),...,x(n) ，是否存在整数a(1),a(2),..,a(n)a(1),a(2),..,a(n) （不全为零），使得a(1)x(1)+a(2)x(2)+...+a(n)x(n)=0a(1)x(1)+a(2)x(2)+...+a(n)x(n)=0对于n=2n=2 ，历史悠久的欧几里得算法能做这项工作、计算x(1)/x(2)x(1)/x(2) 的连分数展开中的各项。如果x(1)/x(2) 是有理数，展开会终止，在适当展开后就给出了“最小的”整数a(1)和a(2) 。欧几里得算法不终止——或者如果你只是简单地由于厌倦计算——那么展开的过程至少提供了最小整数关系的大小的下界。Ferguson和Forcade的推广更有威力，尽管这种推广更难于执行(和理解)。例如，他们的侦查算法被用来求得逻辑斯谛(logistic)映射的第三和第四个分歧点，b(3)=3.544090 和 b(4)=3.564407所满足的多项式的精确系数。(后者是120 阶的多项式；它的最大的系数是257^30 。)已证明该算法在简化量子场论中的Feynman图的计算中是有用的。
10. 快速多极算法
1987年，耶鲁大学的Leslie Greengard 和Vladimir Rokhlin发明了快速多极算法。该算法克服了N体模拟中最令人头疼的困难之一：经由引力或静电力相互作用的N个粒子运动的精确计算(想象一下银河系中的星体，或者蛋白质中的原于)看来需要O(N*N) 的计算量——比较每一对质点需要一次计算。该算法利用多极展开(净电荷或质量、偶极矩、四矩，等等)来近似遥远的一组质点对当地一组质点的影响。空间的层次分解用来确定当距离增大时，比以往任何时候都更大的质点组。快速多极算法的一个明显优点是具有严格的误差估计，这是许多算法所缺少的性质。

求解数值计算问题的一般策略是将复杂或困难的问题用解相同或相近的简单问题代替，这种近似过程通常包括以下几种：

用有限维空间代替无限维空间
用有限的过程代替无限的过程
用代数方程代替微分方程
用线性问题代替非线性问题
用低阶系统代替高阶系统
用简单函数代替复杂函数
用简单结构的矩阵代替一般的矩阵

名称	地址
ming	gams.nist.gov
NetLib	www.netlib.org
MATLAB	www.matlab.com

1.2误差

定义：设准确值xx对应的近似值为x^\hat{x}，则误差（error）

exp(x^)=x^−x

\exp(\hat{x})=\hat{x}-x
它也被称为绝对误差（absolute error）

定义：设准确值xx对应的近似值为x^\hat{x}，则相对误差（relative error）

expr(x^)=x^−xx

\exp_r(\hat{x})=\frac{\hat{x}-x}{x}
定理：设 xx与其近似值x^\hat{x}的第一位有效数字相同，均为 d0d_0，若 x^\hat{x}有 pp位正确的有效数字，则其相对误差满足：

|expr(x^)|≤1d0×10−p+1

\left | \exp_r(\hat{x}) \right |\leq \frac {1}{d_0} \times 10^{-p+1}
定理：设对 xx保留pp位有效数字后得到近似值 x^\hat{x}，则 x^\hat{x}的相对误差满足：

|expr(x^)|≤12d0×10−p+1

\left | \exp_r(\hat{x}) \right |\leq \frac {1}{2d_0} \times 10^{-p+1}
定理：设 xx的第一位有效数字为d0d_0，若近似值 x^\hat{x}的相对误差满足

|expr(x^)|≤12(d0+1)×10−p+1

\left | \exp_r(\hat{x}) \right |\leq \frac {1}{2(d_0+1)} \times 10^{-p+1}
则 x^\hat{x}具有 pp位正确地有效数字，或者在保留pp位有效数字后 x^\hat{x}和 xx的结果相等。

定理：若xx的近似值 x^\hat{x}的相对误差满足

|expr(x^)|≤12×10−p

\left | \exp_r(\hat{x}) \right |\leq \frac {1}{2} \times 10^{-p}
则 x^\hat{x}具有 pp位正确地有效数字，或者在保留pp位有效数字后 x^\hat{x}和 xx的结果相等。

近似数的精度(precision)和准确度(accuracy)两个概念含义接近，但有区别。精度与有效数字的位数有关，而准确度则与准确的有效数字位数有关。

定义：一个数值计算过程的数据传递误差指单纯由问题输入数据误差引起的计算结果的误差

定义：一个数值计算过程的计算误差指计算过程中的近似引起的计算结果的误差

定义：截断误差(truncation error)是指实际输入对应的准确结果与用给定算法精确计算得到的结果之间的差，反映了数值方法对原始数学问题的简化和近似。

定义：舍入误差(rounding error)是指给定算法经精确计算得到的结果与同样的算法经有限精度运算(即舍入运算)得到结果之间的差。

定义：问题的相对条件数：

cond=∥∥问题解的相对变化量∥∥∥∥输入数据的相对变化量∥∥

cond=\frac{\left \| 问题解的相对变化量 \right \|}{\left \| 输入数据的相对变化量 \right \|}
如果一个问题的条件数远远大于1，称这个问题是病态的。

算法的稳定性：

若计算机结果对计算过程中的舍入误差不敏感，则相应的算法为稳定的算法。
对于包含一系列步骤的计算过程，若计算机中的小扰动不被放大(传播)或放大不严重，则相应的算法是稳定的算法。

向后误差(backward error)是指假设计算过程没有近似时，将结果误差折算为初始数据上的扰动，即向后误差是等效的初始数据误差。

1.3计算机浮点数系统与舍入误差

设浮点数系统F\mathbb{F}采用β\beta进制，任意x∈Fx\in \mathbb{F}的浮点数采用类似于科学计数法的表示形式

x=±(d0+d1β+d2β2+⋯+dp−1βp−1)×βE

x=\pm(d_0+\frac{d_1}{\beta}+\frac{d_2}{\beta^2}+\cdots+\frac{d_{p-1}}{\beta^{p-1}})\times\beta^E
其中整数 did_i满足 0≤di≤β−1,i=0,⋯,p−10\leq d_i\leq \beta-1, i=0,\cdots,p-1

整数EE满足L≤E≤UL\leq E\leq U.

一个规范化的浮点数系统由四个整数完全确定：基数(或底数)β\beta、尾数长度(或精度)pp、指数下限值LL和上限值UU。

浮点数系统	β\beta pp	LL UU
IEEE单精度	2	24	-126	127
IEEE双精度	2	53	-1022	1023

1985年以后的计算机都使用IEEE标准2的浮点算术体系

IEEE定义了两个特殊值

Inf：表示无穷

NaN：表示不是数，由没有定义或不确定的操作产生。

这两个特殊值都是通过特殊的指数域取值实现的

定义：

一个浮点数系统中，最小的正数称为下溢值(underflow level,UFL)
一个浮点数系统中，最大的正数称为上溢值(overflow level,OFL)

IEEE 浮点数系统中的几个关键数值	IEEE单精度浮点数	IEEE双精度浮点数
OFL	3.403×10383.403\times10^{38} 1.797×103081.797\times10^{308} UFL 1.175×10−381.175\times10^{-38} 2.225×10−3082.225\times10^{-308} εmach\varepsilon_{mach} 5.960×10−85.960\times10^{-8} 1.110×10−161.110\times10^{-16} 定理：在规范化的浮点数系统中表示非零实数xx，其相对误差限是εmach\varepsilon_{mach}，即： ∣∣∣fl(x)−xx∣∣∣≤εmach \left \| \frac{fl(x)-x}{x} \right \| \leq \varepsilon_{mach} 1.4 保证数值计算的准确性减少舍入误差的几条建议避免中间计算结果出现上溢或下溢避免“大数吃掉小数” 指在做加减运算时，若两个操作数大小相差悬殊，较小数的信息将被较大的数吃掉避免符号相同的两相近数相减注意简化步骤，减少运算次数 the Institute of Electrical and Electronics Engineers ↩ William M. Kahan http:// www.eecs.berkeley.edu/Faculty/Homepages/kahan.html ↩

IEEE 浮点数系统中的几个关键数值

IEEE单精度浮点数

IEEE双精度浮点数

OFL

3.403×10383.403\times10^{38} 1.797×103081.797\times10^{308} UFL 1.175×10−381.175\times10^{-38} 2.225×10−3082.225\times10^{-308} εmach\varepsilon_{mach} 5.960×10−85.960\times10^{-8} 1.110×10−161.110\times10^{-16}

定理：在规范化的浮点数系统中表示非零实数xx，其相对误差限是εmach\varepsilon_{mach}，即：

∣∣∣fl(x)−xx∣∣∣≤εmach

\left | \frac{fl(x)-x}{x} \right | \leq \varepsilon_{mach}

1.4 保证数值计算的准确性

减少舍入误差的几条建议

避免中间计算结果出现上溢或下溢
避免“大数吃掉小数”

指在做加减运算时，若两个操作数大小相差悬殊，较小数的信息将被较大的数吃掉
避免符号相同的两相近数相减
注意简化步骤，减少运算次数

the Institute of Electrical and Electronics Engineers ↩
William M. Kahan http:// www.eecs.berkeley.edu/Faculty/Homepages/kahan.html ↩