西安交通大学计算方法B完整版

第一章绪论

1.1数值计算

现代科学的发展，已导致科学与技术的研究从定性前进到定量，尤其是现代数字计算机的出现及迅速发展，为复杂数学问题的定量研究与解决，提供了强有力的基础。

通常我们面对的理论与技术问题，绝大多数都可以从其物理模型中抽象出数学模型，因此，求解这些数学模型已成为我们面临的重要任务。

一、本课程的任务：

寻求解决各种数学问题的数值方法——如何将高等数学的问题回归到初等数学（算术）的方法求解——了解计算的基础方法，基本结构（否则只须知道数值软件）——并研究其性质。

立足点：

面向数学——解决数学问题

面向计算机——利用计算机作为工具

充分发挥计算机的功能，设计算法，解决数学问题

例如：迭代法、并行算法

二、问题的类型

1、离散问题：例如，求解线性方程组b

Ax=——从离散数据：矩阵A和向量b，求解离散数据x；

2、连续问题的离散化处理：例如，数值积分、数值微分、微分方程数值解；

3、离散问题的连续化处理：例如，数据近似，统计分析计算；

1.2数值方法的分析

在本章中我们不具体讨论算法，首先讨论算法分析的基础——误差。

一般来讲，误差主要有两类、三种（对科学计算）：

1）公式误差——“截断误差”，数学?计算，算法形成——主观（人为）：数学问题-数值方法的转换，用离散公式近似连续的数学函数进行计算时，一般都会发生误差，通常称之为“截断误差”；——以后讨论

2）舍入误差及输出入误差——计算机，算法执行——客观（机器）：由于计算机的存储器、运算器的字长有限，在运算和存储中必然会发生最末若干位数字的舍入，形成舍入误差；在人机数据交换过程中，十进制数和二进制数的转换也会导致误差发生，这就是输入误差。这两种误差主要是由于计算机的字长有限，采用浮点数系所致。

首先介绍浮点数系

一、计算机上的运算——浮点运算

面向计算机设计的算法，则先要讨论在计算机上数的表示。科学记数法——浮点数：约定尾数中小数点之前的数全为零，小数点后第一个数不能为零。 目前，一般计算机都采用浮点数系，一个存储单元分成首数和尾数：

首数 尾数(位)

其中首数存放数的指数（或“阶”）部分，尾数存放有效数字。对于β进制，尾数字长为t 位的浮点数系),,,(U L t F β中的（浮点）数，可以用以下形式表示：

t

j j

d d l t

t d d

d x fl ,,3,2,

01

1)221()(ΛΛ=<≤<≤?++±=ββββββ

此处，指数l （称为阶）限制在U l L ≤≤范围内。

以下记实数系中的实数为R x ∈，它在浮点数系),,,(U L t F β中对应的浮点数记为),,,()(U L t F x fl β∈——β进制，t 尾数位数，U L ,阶的范围。

几乎所有近代计算机都采用“二进制”（即2=β）：位、字节和字分别是指位数不同的二进制数。例如

位是一个二进制数（即0或1）；字节是8个二进制数字；上表的最后一列是字节。

单精度浮点数（single precision ）按32位存储，双精度浮点数（double precision ）按64位存储。精度用于指明每个浮点数保留多少位以及尾数和阶数各分配多少位。单精度浮点数的尾数为23位、阶数为8位；双精度浮点数的尾数为53位（包含符号位）、阶数为11位（包含符号位）。双精度浮点数的等价二进制数如下所示：

{4444448444444764434421Λ43421Λ位位尾数

符号位

位指数（含符号位）645211ddddd ddd

f bb bbbb 浮点数的特点：

1、 实数转换到浮点数——浮点化，〈缺点：〉总会产生误差（除极个别的情

况：Λ,2,1,0,2±±==l x l ）

按 四舍五入，绝对误差：t l x fl x -≤-β2

1

)(（举例），

〈优点：〉浮点化产生的相对误差有界（与数字本身的数量级无关）

t x x fl x x -≤-=

12

1

)()(βδ 注：设实数R x ∈，则按β进制可表达为：

Λ

ΛΛΛ1,,,3,2,01

1)1

1221(+=<≤<≤?++++++±=t t j j

d d l

t t d t t d d

d x ββ

βββββ按四舍五入的原则，当它进入浮点数系),,,(U L t F β时，若β2

1

1<

+t d ，则 l t

t d d

d x fl ββββ?++±=)221(

)(Λ

若β211

≥+t d ，则 l t

t d d d x fl ββββ?+++±=)1221()(Λ

对第一种情况：

t l l

t

l t t d x fl x -++=?≤

?+=-βββββ21)2

1(1)(

)(1

1

Λ 对第二种情况：

t l l

t

l t t d x fl x -++=?≤?--=

-βββ

βββ21)21(1)(11Λ 就是说总有： t

l x fl x -≤

-β2

1)( 另一方面，由浮点数要求的 β<≤11d ， 有l x ββ

1

≥

，将此两者相除，便得

t x x fl x -≤-12

1

)(β

2、每一个浮点数系的数字有限： 1)1()1(21++---L U t ββ

3、浮点数系中的运算非自封闭,（因为数字有限等）

例：在)5,5,4,10(-F 中，32102001.,105420.?=?=-y x ，运算y x *和y x /,的结果显然已不在此浮点数系内，而y x +或y x -也不在此浮点数系内，需

)(y x fl ο结果才在此浮点数系内。

浮点运算应注意：

1）避免产生大结果的运算，尤其是避免小数作为除数参加运算； 2）避免“大”“小”数相加减；

3）避免相近数相减，防止大量有效数字损失； 4）尽可能简化运算步骤，减少运算次数。 原因：

记t

x x fl x x -=-==12

1)(max

)(max βδ?，由x x fl x ?≤-)(，可得： y x y x fl y x οοο?≤-)()( （“。”表示任意一种四则运算）

此处 ? 是由机器字长（实质上是尾数字长t 的大小）确定的常数（它反映了实际运算的精度）。

显然，若要求运算的舍入误差小，应使运算结果（如：y x y x y x ,,?±）较小。 尤其是小分母运算：

y

y x y x δ

δ=-+，小y ?大误差。 其次，当浮点数系中两个数量级相差较大的数相加（或减），注意到由于浮点数系中数字字长的有限性，可能导致“大数吃小数”。例如，)5,5,4,10(-F 中

13102317.,108231.-?=?=y x ，则

x y x =?±?=?±?=±-3313100000.108231.102317.108231.

似乎)(x y y <<没有参加运算。

第三，同样，由于浮点数系中数字字长的有限性，当两个相近数相减时：例如，在)5,5,8,10(-F 中，331082317832.,1082317844.?=?=y x ，两数相减：

31012000000.-?=-y x ，计算结果仅剩2位有效数字，而原来参加运算的数字有8

位有效数字，这将严重影响最终计算结果的精度。

二、 算法分析

作为一个可用的算法，必须考虑其效率和可靠性，

定义：计算机完成一个乘法和一个加法，称为一个浮点运算（记为flop ）； 注：由于计算机在运算时，加（减）法所耗时间远少于乘（除）法，所以通常只须计算乘法的次数，因此也有“一个算法需要多少个‘乘法’”这一提法。

1、计算效率——可计算性（计算复杂性——空间、时间） 例：解线性方程组b Ax =的Grammar 方法：A

A i

i x =

，其中A 是方程

组系数矩阵A 对应的行列式，而i A 则是以右端向量b 替代A 的第i 列所得矩阵对应的行列式。由线性代数知识可知，若)(ij α=A ，则

∑-=

n

n

ni

i i i i i J αααΛΛ2

1

2

121),()1(A ，

其中),,,(21n i i i J Λ是由{n ,,2,1Λ}变换到{n i i i Λ,,21}所需的置换次数。可见每计算一个行列式，需要!)1(n n ?-个浮点运算；因此，按Grammar 方法解方程组约需!)1(2n n N ?-= 个浮点运算。当20=n 时 201070728.9?≈N ，用一个运算速度为秒/108的计算机进行求解，约需510078.3?年（日前报道我国计算机已达到/1038408?秒，这仍需近10年）。而n=20的方程组应该说是一个小型的方程组。因此Grammar 方法是一个不能接受的算法，它缺乏可计算性。第二章将介绍的Gauss 消去法和迭代法就有较高的效率，具有很好的可计算性。

2、计算可靠性

作为算法，除了考虑其效率外，必须重视可靠性，它包括两方面： 问题的性态 和 方法的稳定性 问题的性态

所计算的问题当原始数据发生小扰动时，问题的解一般也发生扰动。问题的性态——问题的解对原始数据发生变化的敏感性。

原始数据小扰动?问题解 ???—问题是病态的—大变化—问题是良态的

—小扰动

例：线性方程组：

?????????=++=++=++60475141311213413121

6113121321321321x x x x x x x x x 的解是：????? ??=111X

若将方程组的系数改写成具有2位有效数字的小数：

?????=++=++=++78

.020.025.033.01.125.033.050.08.133.050.000.132132121x x x x x x x x x 的解则变成：???? ??--=65.3325.3822.6~X ；

这是一个典型的病态方程组。

一般：由原始数据 ?x 计算结果)(x f ，

扰动后的数据 ?x ~计算结果)~(x f ， 若对问题f 存在常数m ，满足关系式：

x

x

x m

x f x f x f ~)()~()(-≤- 或 x

x

x m ~ f(x)

)

x ~f(-f(x)sup -= 则称（相对误差之比的上界）m 为该问题的条件数，记作 )(f cond m =；

由微积分中值定理知识容易计算出，当x x ~≈时，)

()

()(x f x f x

f cond '≈ 。 稍后我们在第二章将对线性方程组的性态作进一步的分析。 算法的数值稳定性：

例：计算8,,1,0,1

1Λ==?-n dx e x e I x n n ；

解：由微积分知识，可得计算公式，① 11--=n n nI I ，②)1(1

1n n I n

I -=-，我们将准确值与计算结果列表如下：

由上表可见，方法①中，原始步的误差，随着计算步数的增加被严重地放大，特别是8I 竟变成负数（注意：被积函数是非负函数），而方法②则相反；这是因

为方法①中，若前步有误差δ：δ+=--11~

k k I I ，则

δδk I k kI I k I k k k k -=--=-=--111~

1~， 说明1-k I 的误差δ，经一步计算后被扩大了k 倍，随着k 的增大，误差将被大大地扩大；而通过同样的分析可知：方法②中k I 的误差则被缩小k 倍。

算法的数值稳定性：算法对初始误差导致的最终误差的可控性，如果最终误差被有效地控制，则称算法是稳定的，否则就是不稳定的。

第二章 线性方程组求解

线性方程组：

???????=+++=+++=+++?=)1.2(22112

22221

2111212111n

n nn n n n n n n x x x x x x x x x βαααβαααβαααΛΛΛΛΛΛΛΛΛb Ax

其中

.,,212

12

122221

11211?????

?

??????=????????????=??

?

????

?????=n n nn n n n n b x x x x A βββαααααααααM M Λ

ΛΛΛΛΛΛ

2.1 Gauss 消去法

2.1.1 消去法

设计方法原则：面向计算机（事先未知元素，编制程序），例：

1123233

203212221212

1=?

??

??

??????=?-=-=+?=-=+x x x x x x x x x

基本思想：降维——N 维问题转化为N-1维问题——逐次降维，依次进行 消去过程——对方程组（2.1）由上而下逐步消去对角元)1,,2,1(-=n k kk Λα 以下的 ),,1(n k i ik Λ+=α，使之转化为如下等价形式，达到目标：

)2.2(222221

1212111??

?

?

??

?==++=+++n n nn n n n n x x x x x x βαβααβαααΛΛΛΛΛΛ

从而，可进入回代过程，再由下而上，逐步解得 )1,2,,1,(Λ-=n n k x k ：

这儿的)1,,2,1(-=n k kk Λα——主元

对问题（2.1）设011≠α，目标：将第2～n 方程的1x 的系数121,,n ααΛ转化

为0；方法：“第k 个方程”-?11

1αα

k “第1个方程”（),,3,2n k Λ=，得到

??

?

???

?=++=++=+++)1()1(2)1(2)

1(2

)1(22)1(2211212111n n nn n n n n n x x x x x x x βααβααβαααΛΛΛΛΛΛΛΛ

现在只须关心由第2～n 方程形成的n-1维方程组，以后可仿此进行。

消去：第k 步)1,,2,1(-=n k Λ：设0≠kk α，以第k 行为基准，消去以下各 行中k 列kk α以下的),,1(n k i ik Λ+=α，令 ,kk

ik

ik l αα=施行：第k 行—?ki l 第i 行 ?新的第k 行，

元素变化：（0?kk α），k ik i i kj ik ij ij l l βββααα-=-=~,~，

经过1-n 步消去（注意：nn α以下无元可消），得到)2.2(式。〈注意：每计算1个i ij ik l βα,,仅需1个浮点运算〉

回代：第一步 nn

n

n x αβ= ， 第二步 Λ,][1,1,11

-----=

n n n n n n x n x ααβ，

第1+-k n 步 ;1,2,,1)]

([11ΛΛ-=++-=++n k x x x kk

n kn k k k k k αααβ

运算量：

消元：n 元方程组：第1步消元：对第),,3,2(n i i Λ=行，共n-1行； 每1行：计算),,3,2(1n i l i Λ=，1个乘除法（或“浮点运算”），

计算新的11)1(11)1();,,,3,2(βββαααi i i j i ij ij l n j l -==-=Λ共 1)1(+-n 个

乘除法（或“浮点运算”），

第1次元共1]1)1(1)[1(2-=+-+-n n n 个乘除法，因此，消元的运算量

n n n n k k k n

k n k n k 6523)1()1(2312

1222-+=-=-=-∑∑∑===； 回代：第1步，求n x 需要1个乘除法（或“浮点运算”），即：第2步求1

-n x 用2个乘除法（或“浮点运算”），一般，第k 步求1+-k n x 用k 个浮点运算，因此，回代的运算量

)1(21

+=∑

=n n

k n

k ； Gauss 消去法 的 总运算量：3

323n

n n T -+= 。

例2-1：解线性方程组

???

??=--=++=++2333220221

321321x x x x x x x x 解：利用增广矩阵（因为线性方程组的解仅与系数与右端常数项相关）

????

? ??-=??????

? ?

?-???

?? ??--????? ??--=→-=→=111212

131201

21211312012120313322

01212112323121x l l l 例： ???

?? ??---450433562112 ；解?????

? ??-=329419x

???

??

??

??----??????? ??--=???????

??------551214712230

12120113112154212132225201144230

12

；解??????? ??=1111x ：???????

?

?-??????? ??-=

???????

??--105424

1711343

1

02112311112160185669463225713443

1

2；解 ??

??

???

??-=20

51x ??

?

?

?

??

????????? ??=???????

?

?242418246812622432116711311111902568116144642781101694124321；解??????

? ??--=1111x

??????? ??-------??????? ??=??????? ??41711203213043211134123121577554435443405432304321

；解????

?

?

?

??=43

21x

???????

?

?---??????? ??--=??????? ??--------11253

312170

2

21130112011121411214772010143170221；????

???

??-=1121x

2.1.3 （列）主元消去法

算法中，若第k 步：0)=kk i α

或 n k i ii ik

kk ,,1)Λ+=<<αα

则按照原来的简单Gauss 消去法算法可能无法执行，也可能出现大误差：

例：在浮点数系)10,10,5,10(-F 中计算；方程组

???? ??=???? ?????? ??-2132210215x x 准确解：??????=*499998749.025*******.0x 解：按Gauss 消去法解：

????→????

?

?

?*******=-6

1002000021111114

1020000.1030000.1020000.1010000.1020000.1010000.l ???

?

?

?*-*-***-66

114

1020000.1040000.1010000.1020000.1010000. T )1050000.1000000.(~00

*?=x

误差大，原因：若有误差 δββδαα+?+?k k kj kj , 则

k ik i i kj ik ij ij l l βββααα-=-=~,~，则演变为

δβδβββδαδαααik

i k ik i i ik

ij kj ik ij ij l l l l +=+-=+=+-=~

)(?,~)(?； 这说明前步各元素的误差均被放大ik l 倍。克服方法，将按绝对值最大的元素交换到主元位置，使1≤ik l ，使前步的误差不再被放大。

???? ?

?*****????→????

?

?

?******??→????

?

?

?******-?*=--111111050000.21114111

111114

1010000.1020000.1020000.1030000.1020000.1010000.1020000.1010000.1020000.1030000.1020000.1020000.1030000.1020000.1010000.1020000.1010000.5

21l 行

行 T )1050000.1025000.(~

00*?=x

消元过程中，第k 步消元)1,,2,1(-=n k Λ以第k 行为基准，消去其后的k n -个方程的ik α（系数矩阵第k 列kk α以下的元素）

，Gauss 消去法要求主元0≠kk α。为避免出现0=kk α作为主元，并使前步的误差不被放大，应使1≤ik l ，为此应使：},,{,1k n k k kk kk Max ααααΛ+=，通常采用按列选主元的列主元消去法：若},,{1k n k k kk k m Max ααααΛ+=，便将第k 行与第m 行交换，使k m α与kk α交换位置，使新的第k 行执行在原始Gauss 消去法中的角色，保证将作为除数的主元n k i ik

kk ,,1Λ+=≤αα，从而1≤ik l ，重复前述的Gauss 消去过程。

列主元消去法的效果：

1. 算法稳定，即计算误差能被有效控制；

2. 当矩阵A 的行列式0≠A 时，算法总可以执行完成；

注：若矩阵A 是对称正定或严格对角占优，则不选主元，消去法也是稳定的；