西安交通大学计算方法B.docx

《西安交通大学计算方法B.docx》由会员分享，可在线阅读，更多相关《西安交通大学计算方法B.docx（145页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第一章绪论1.1 数值计算现代科学的发展，已导致科学与技术的研究从定性前进到定量，尤其是现代数字计算机的出现及迅速发展，为复杂数学问题的定量研究与解决，提供了强有力的基础。通常我们面对的理论与技术问题，绝大多数都可以从其物理模型中抽象出数学模型，因此，求解这些数学模型已成为我们面临的重要任务。一、 本课程的任务：寻求解决各种数学问题的数值方法一如何将高等数学的问题回归到初等数学（算术）的方法求解了解计算的基础方法，基本结构（否则只须知道数值软件）并研究其性质。立足点：面向数学解决数学问题面向计算机利用计算机作为工具充分发挥计算机的功能，设计算法，解决数学问题例如：迭代法、并行算法二、 问题的类

2、型1、离散问题：例如，求解线性方程组Ax = b从离散数据：矩阵A和向量b,求解离散数据x；三、 连续问题的离散化处理：例如，数值积分、数值微分、微分方程数值解;四、 离散问题的连续化处理：例如，数据近似，统计分析计算；1.2 数值方法的分析在本章中我们不具体讨论算法，首先讨论算法分析的基础一误差。一般来讲，误差主要有两类、三种（对科学计算）：1）公式误差“截断误差”，数学c计算，算法形成主观（人为）：数学问题-数值方法的转换，用离散公式近似连续的数学函数进行计算时，一般都会发生误差，通常称之为“截断误差”；以后讨论2）舍入误差及输出入误差计算机，算法执行一客观（机器）：由于计算机的存储器、运

3、算器的字长有限，在运算和存储中必然会发生最末若干位数字的舍入，形成舍入误差；在人机数据交换过程中，十进制数和二进制数的转换也会导致误差发生，这就是输入误差。这两种误差主要是由于计算机的字长有限，采用浮点数系所致。首先介绍浮点数系一、计算机上的运算浮点运算面向计算机设计的算法，则先要讨论在计算机上数的表示。科学记数法浮点数：约定尾数中小数点之前的数全为零，小数点后第一个数不能为零。目前，一般计算机都采用浮点数系，一个存储单元分成首数和尾数：XXXXXXX首数/尾数（/位）其中首数存放数的指数（或“阶”）部分，尾数存放有效数字。对于。进制，尾数字长为t位的浮点数系尸（6,/,。）中的（浮点）数，可

4、以用以下形式表示：d. d d J d. Pflx）=（+-）xR . n a .P p2,0djJ3, j =2,3,/此处，指数/（称为阶）限制在范围内。以下记实数系中的实数为xe R ,它在浮点数系尸（4/心。）中对应的浮点数记为川E）eF（0,t,L,U）进制，尾数位数，阶的范围。儿乎所有近代计算机都采用“二进制（即夕=2）：位、字节和字分别是指位数不同的二进制数。例如十进制转换二进制11=2。0000000122=20000001044=220000010088=230000100099=2+2000010011010=23+2000010102727=16+8+2+1=24+23+

5、21+2000110111字节位是一个二进制数（即0或1）；字节是8个二进制数字；上表的最后一列是字节。单精度浮点数（single precision）按32位存储，双精度浮点数（double precision）按64位存储。精度用于指明每个浮点数保留多少位以及尾数和阶数各分配多少位。单精度浮点数的尾数为23位、阶数为8位；双精度浮点数的尾数为53位（包含符号位）、阶数为11位（包含符号位）。双精度浮点数的等价二进制数如下所示：64位bbbb bb f ddd ddddd 11位指数（含符号位）符旗位52位尾数浮点数的特点：1、实数转换到浮点数浮点化，缺点：总会产生误差（除极个别的情况：x

6、=2/,/=0,l,2,-）按四舍五入，绝对误差：（举例），优点：浮点化产生的相对误差有界（与数字本身的数量级无关）I (5(x)|= X二胆)T注：设实数xeH,则按进制可表达为： d d d.14d 4P 夕2pt pt 4-1。 dj /，J =2,3,JJ +1按四舍五入的原则，当它进入浮点数系尸（,/,U）时，若4+1；夕，则1 fI7（x）=（-+-|+）xp/P外衣1 4 d 2+1/若4+i N 彳夕，贝U 力。）=（大+-57+-）x2 夕/2,对第一种情况：卜-/。）|=（翳+-”4（g）x =；，-对第二种情况：l-v-./M=r-一x号（；）x/#就是说总有：另一方面，

7、由浮点数要求的144/，有国之夕，将此两者相除，便得x-/（x） w L Ix 2P2,每一个浮点数系的数字有限：2（尸-1）夕t（U-L +1）+13、浮点数系中的运算非自封闭，（因为数字有限等）例：在尸（10,4,-5,5）中，x =.5420x10-2/=.2001x103,运算x*y 和 x/y,的结果显然已不在此浮点数系内，而x + y或x-y也不在此浮点数系内，需力（x y）结果才在此浮点数系内。浮点运算应注意：1）避免产生大结果的运算，尤其是避免小数作为除数参加运算；2）避免“大”“小”数相加减；3）避免相近数相减，防止大量有效数字损失；4）尽可能简化运算步骤，减少运算次数。原因

8、：记= max,（x）|= max力（*）=,由4八国，可得：|（x y）- y）！个浮点运算；因此，按Grammar方法解方程组约需N =(2_1).！个浮点运算。当=20时 N工9.70728 x 1O20,用一个运算速度为103秒的计算机进行求解，约需3.078x1()5年(日前报道我国计算机已达到3840x1(/秒，这仍需近10年)。而n=20的方程组应该说是一个小型的方程组。因此Grammar方法是一个不能接受的算法，它缺乏可计算性。第二章将介绍的Gauss消去法和迭代法就有较高的效率，具有很好的可计算性。2、计算可靠性作为算法，除了考虑其效率外，必须重视可靠性，它包括两方面：问题的

9、性态和方法的稳定性问题的性态所计算的问题当原始数据发生小扰动时，问题的解一般也发生扰动。问题的性态问题的解对原始数据发生变化的敏感性。原始数据小扰动=问题解小扰动问题是良态的大变化问题是病态的例：线性方程组:1 111X, 4d_ 1 2 2 3 361 1 113234 121 1 1_ 47-X, + -+ Xq 345- 60的解是：X =若将方程组的系数改写成具有2位有效数字的小数:f 1 .OOx, + 0.50x2 + 0.33x = 1.8 0.50X +0.33x2 +025x3 =1.1 o.33X +0.252 + 0-20x3 = 0.78,-6.22、的解则变成：X =

10、 38.25、33.65,这是一个典型的病态方程组。一般：由原始数据x =计算结果/(x),扰动后的数据X n计算结果,/(三),若对问题/存在常数m,满足关系式:/(x)-/(T)/(x)m xf(x)-f(2Tf(X)取m = sup；二jx -XI则称(相对误差之比的上界)m为该问题的条件数，记作m = cond(f)；由微积分中值定理知识容易计算出，当x乏工时，cond(f) xZ宜 o/(x)稍后我们在第二章将对线性方程组的性态作进一步的分析。算法的数值稳定性：例：计算/“=,=0,1,8；0我们将准确值与计算结果列表如下：解：由微积分知识，可得计算公式，/“=1-/“t ,/方法/

11、()AAA17准确值.6321.3679.2642.2073.1709.1455.1268.1124.1010.6321.3679.2642.2074.1704.1480.1120.2160-.7280.6320.3680.2643.2073.1709.1455.1268.1124.1010n由上表可见，方法中，原始步的误差，随着计算步数的增加被严重地放大,特别是A竞变成负数(注意：被积函数是非负函数)，而方法则相反；这是因为方法中，若前步有误差S：7i=/i+S,则4=1一用=1一亿T 一必=4k3,说明的误差6,经一步计算后被扩大了左倍，随着左的增大，误差将被大大地扩大；而通过同样的分析可

12、知：方法中的误差则被缩小左倍。算法的数值稳定性：算法对初始误差导致的最终误差的可控性,如果最终误差被有效地控制，则称算法是稳定的，否则就是不稳定的。第二章线性方程组求解线性方程组：|X|+|22+(nXn =0、(2.1).-02M1+。22工2+=/2Ax = b u33nx2=l = X=lX|-X2=0-x2=-基本思想：降维一N维问题转化为NT维问题逐次降维，依次进行消去过程对方程组（2.1）由上而下逐步消去对角元e=1,2,，-1）以下的（i = A + l,），使之转化为如下等价形式，达到目标：|4内+anx2 h1-aynxn =0、a22x2+-+ a2nxn=2a”=瓦从而，

13、可进入回代过程，再由下而上，逐步解得xk （k = m,m- I,-,2,1）：这儿的a*/=1,2,-1）主元对问题（2.1）设工0,目标：将第2n方程的西的系数a21,a”】转化为0；方法：“第个方程”-3X 第1个方程”（无=2,3,）,得到孙QX+ a12X2+ OLXnXn = P姗+ a如=砂新的第4行，元素变化：（%*=0）,即= a厂/*%,与=Bj-likPk，经过T 步消去（注意：以下无元可消），得到（2.2）式。注意：每计算1个%,回仅需1个浮点运算回代：第一步x”=五,ann第二步X,t优一%一瓜%第-+1步=自一014+1+k = n-,-,2,1;运算量：消元：n元

14、方程组：第1步消元：对第i（i =2,3,行,共n-1行;每1行：计算/,（i =2,3,），1个乘除法（或“浮点运算”），计算新的谭=%-/“（_/=2,3,）;夕尸=4-/共（n-l）+ l个乘除法（或“浮点运算”），第1次元共（-1）口+（-1）+1=2一1个乘除法，因此，消元的运算量2 曾3扇 s（炉-1）工（炉-1）=6父-=可+可-z； k=nhlhl32。回代：第1步，求X”需要1个乘除法（或“浮点运算”），即：第2步求X”T用2个乘除法（或“浮点运算”），一般，第步求x._i用个浮点运算，因此,回代的运算量=+1）;k=l 2Gmss消去法的总运算量：T =-+ n2-。3 3

15、例21：解线性方程组Xj +2%2+工3=02%|+2x2+3x3=3-%!-3x2=2解：利用增广矩阵（因为线性方程组的解仅与系数与右端常数项相关）20134、1、20134、-1、431725213-11175;解x =2-364-91-114240、6651860,3211,、510,12/2342、fl2342、f-1)14916101126128；解X 二1182764441316241SJ 1681256190,762424,b J123430、勺、(123430)23454021-1-2-3-20234=:解*=4543321-1-1-73心55757)、43114)12-207

16、、X12-2071、-1-341-10-11-121一320-277-4021352；X =-1J2-11142;0一311121=2(o)-1一3-2-12J-2-1例: 2 -112、-6 5 -3 34054,解* =(%-3-2 1 032、 1 、-2 1-441 -102 122 5 2 -22 -13-1 31-21245 ,J 0 -22.1.3(列)主元消去法算法中，若第左步：1)%*=0或 n)akk|aik| i = k + l,-,n则按照原来的简单Gaz/ss消去法算法可能无法执行，也可能出现大误差:例：在浮点数系尸(10,5,-10,10)中计算；方程组IO-2号)

17、(H.0.250001875=准确解：x =I 23人J(2)0.499998749解：按Gauss消去法解：.10000*1()720000*10.10000*0如。2a”、20000*10.30000*10.20000*101 J-pOOOO+lO4.20000*10,10000*101、-.40000106-.20000*10%x =(.0000010.50000*10)r误差大，原因：若有误差a n %+5,瓦=片+ S则%=%-4%，P,= A -1淳瓦，则演变为%=%-4(%+6)=玛+/*，Bi =力-4(乩+6)= A +(3；这说明前步各元素的误差均被放大4倍。克服方法，将按

18、绝对值最大的元素交换到主元位置，使|4|41,使前步的误差不再被放大。.10000*107.20000*10.10000*10)由3万.20000*10.30000*10.20000*10广r .20000*10.30000*10.20000*10，小即。、loooo*io-，.20000*10.loooo*io?-.20000*1()1.30000*10.20000*10、.20000*10.10000*10JH =(.2500010.50000*10)r消元过程中，第4步消元(左=1,2,，-1)以第行为基准，消去其后的-4个方程的*(系数矩阵第上列出”以下的元素)，Gauss消去法要求主

19、元为避免出现a版=0作为主元，并使前步的误差不被放大，应使|/林区1,为此应使：|%*|=加以|七*|,|%+4-,”1,通常采用按列选主元的列主元消去法：若%/=朋1aM%/,鼠+1小.也,*|,便将第左行与第加行交换，使与交换位置，使新的第左行执行在原始Ga“ss消去法中的角色，保证将作为除数的主元i = k + l,n ,从而重复前述的Gaass消去过程。列主元消去法的效果：1 .算法稳定，即计算误差能被有效控制；2 .当矩阵力的行列式同0时，算法总可以执行完成；注：若矩阵”是对称正定或严格对角占优，则不选主元，消去法也是稳定的;矩阵严格对角占优的定义：对角元的绝对值大于该行其他元的绝对

20、值之和,n 即同同;7=1评问题：为什么系数矩阵Z严格对角占优，则不选主元的消去法也是稳定的?为保证消去法是稳定的，算法应如何进行？注意：第一步消元&产K-以=由于占优：& a2j a2y + a21awa, , Z%/ Z%/ +Z j3j=3j3aij a2I1=Z a2j + %j=3Inj3n4-2/ 丁 2111ina2 + Z 沏7=342一% 即z% =7=3 X a2j即新的第二行的对角元绝对占优，其他也可同样证明J-3例2-2：列主元消去法求解例2-122-3121223、一 1 -3 00( 23 T 1c 1行“2行 1同样得到原方程组的解x = (l -I if；2.2

21、矩阵分解2.2.1 Gauss消去法的矩阵意义矩阵的三角分解若将Gmss消去法解方程过程中产生的（i力作为一个单位下三角矩阵考虑方程的系数矩阵部分）;其对角元为1,对角线上的元素均为0 Z,的对应元素；将Gai/ss消去法消元过程最后得到的上三角矩阵对角线以下元素均为0一记作万或U （仅i、210、T2I0、lU =21213=2233L % I%C-302如例2.1,再将L与或U相乘:显然ZU相乘得到的正是原始所求解的方程的增广矩阵，而LU相乘得到的正是原始所求解的方程的系数矩阵。反过来，一个矩阵也可以通过Gmss求解的过程获得这样的“LU分解”，其中L为单位下三角矩阵，U为上三角矩阵。这样

22、将一个矩阵分解成个下三角矩阵和一个上三角矩阵的乘积形式，称为矩阵的三角（Doolittle）分解。1）消去法的矩阵意义：以例2-1说明1-1 -3 0 2J与以下矩阵乘法是一样的:2 1 0、-2 1 3T 1 2/Y 1210、223 31 11-1 -3 0 2.12 10、-2 1 3、T 1 2,可见，一般的第一步消元是：1 a aLXA,b）= lil 1 ,即即A -A : 为?： 电E : 即O : I =A A : ll，H : a a二A ;即 ob)=(Uy)其中。为上三角矩阵，y为向量：P Al 12222.2.2 矩阵分解/= lu注意到初等矩阵。具有性质:因此，我们有

23、 (4。)=匕：匚3 (U, y) = L (U, y)根据矩阵乘积的性质，有 A = LU这就是基本的矩阵三角分解Doolittle分解将矩阵分解为单位下三角矩阵L与上三角矩阵U的乘积。从矩阵乘法与行列式的关系可知，若矩阵A三角分解得到A=LU,则有:det A = det L * det U ,由于下三角矩阵或上三角矩阵的行列式是全部对角元的乘积，因此可利用三角分解求矩阵对应的行列式的值。例如：上述例2.1方程组系数矩阵对应的行列式的值是-1,下例2.3对称正定矩阵对应的行列式的值为144o当系数矩阵为的方程组可以顺利求解(不必选主元)时，解的过程显然是唯一的，因此它的Doolittle分

24、解必然也是唯的，从而可以通过待定系数的方法获得该矩阵的Doolittle分解（通常称为“直接分解”或“紧凑格式”）。2.2.3 其他三角分解除了上述的单位下三角矩阵与上三角矩阵的乘积形式以外，还有其他类型的分解，例如下三角矩阵和单位上三角矩阵的乘积，单位下三角矩阵与对角矩阵、单位上三角矩阵的乘积L0ML对于对称矩阵，可以分解成LO/7矩阵分解的形式，特别是对称正定矩阵，可以分解成GG，形式（称为Cholesky分解），其中G 为下三角矩阵由Doolittle分解的唯一性可知这些形式的分解也是唯一的。这些不同的分解都可以通过矩阵乘积的方法取得。下面例2.3以对称正定矩阵为例,说明通过矩阵乘积的方

25、法取得矩阵分解的不同形式。当然，当矩阵的Doolittle分解存在唯一时，这些不同的分解分别时唯一的，因此可以通过待定系数的方法取得，也即通过直接分解的方法取得这些分解。例 2-3：由此可知:4、)=1)。以例2-3矩阵分解为例，容易得到:由此，也很容易看到，一个对称矩阵通过Gmss消去过程得到的，U分解(其中，为单位下三角，。为上三角矩阵)，当U的对角元全部为正数时，该对称矩阵必为对称正定矩阵。由Gass消去的过程可知各次主元是”,22,，（矩阵u的对角元），可知当矩阵/的各阶顺序主子式均非零时，（原始的）Gmss消去法可以顺利完成,这是因为当各阶顺序主子式均非零时，从1，22,，”均非零，

26、即Gmss消去法可以顺利完成。进一步，当矩阵/非奇时，列主元Gmss消去法必可顺利完成：因为当矩阵/非奇时，/的第一列必不全为零，故通过选主元，可进行第一步消元,即算法L/P/4可执行,得到20,且其下方全为零，因为det（L/PE）= detZwO,所以11的代数余子式det/：也非零，（.detN = det（Z1/）=i det/；i #0），即矩阵非奇，因此下一步列主元Gmss消去可进行，由此，可完成全部消元过程。2.2.5带状矩阵的分解jI、X X XX.上带宽夕:ji + q %=0- I*下带宽p : i j + q ay =0xX X X,定理2.5若/=（DoolMe分解，则

27、L（下）带宽p,U（上）带宽g.证明：对二和三阶矩阵显然.（确定），对矩阵的阶作归纳证明:设对-1阶矩阵结论成立.令可验证（介绍分块运算）:由归纳假设=与幼，因此最后的矩阵乘法：前者是初等矩阵左乘实际上是Gauss变换矩阵相乘，后者按分块运算容易得到。特殊情况,p = l,0,|x|=0 x =0, Vx e R,2) 正齐性|ox|=|ar|x|, Va e R,xe Rn ,3) 三角不等式|x + y|x|+|y|,x,yeRn ；常用的范数1-范数|4=同+闯+%,2-范数|电=(,+,2+,/)2,8-范数 IxIL = max(|xj,|x2|,|x|),注：一般若不指某特定的范数

28、，则范数符号不作下标，记为|卜|;2、矩阵范数矩阵，与之相应的一个非负实数，记作M|,如果它满足条件：1) 非负性同，M|=0o/=0,/AeL(Rmxn),2) 正齐性 H|=|4 VaeRZAeLW),3) 三角不等式|卜+训引闻+|同，4) pB|j|B|.AeL(Rmxn),BeL(Rnxp),因为矩阵与向量相乘仍是向量，因此：特别要求一种矩阵范数与一种向量范数:5) px|同”=|心|-范数，口|卜=|人创7卜矩阵范数，其中M-非奇矩阵3、范数的等价性：IHL4M 4的叱，力嗣m2引刈，hL H2|2|x|Zr|=|/lx|2| p(/l)1),而从推导可知，(此处表示“相当于。因此

29、一般总有科普卜卜 L 5.1(1.001_2.000)卜 I上例中：八=44=3x3.001/0.003=30010.0031-0.4991.000)11 I、对于更一般的矩阵力和右端向量均发生扰动/和1的情况，有H Cond(A) M 3W囱丽M证明：恸(1-5)-,存在，且1(15|上品反证：若（1-B）奇异，即|/-网=0,则（1-3=。有非。解= Bx = x 大0又。=|0|= II1-5同|=卜-取| A 国-|忸Ml AIIMHWMI =。-|冏）国|。矛盾，所以（1-5）非奇，。-8尸存在，由恢=1（至少对常用的三种范数）:1= UII =|（/-5尸（/-砌=（/-By 8|

30、N |。-8尸卜|（/-5尸炳=|（1-）i（l-M）由此得若：（/+/Ax + Jx）=（A + Ab）, Ax = b ,.(+ AA )dx - Ab - AA - x /x =(4+ AA )-*Ab -AA-x)/(A + AA)= A(I + A1 AA)|jx|(/+ A_1 AA A-1(4b -AA x.iiii Iklhii f iiii , nin玉-/网间wj WHIhhM/肛F例IM Cond(A)(Ml 可一(徘4!旧WMil2 .4线性方程组的迭代解法迭代法：将方程组Ax = b （2.1）转化为等价的形式x = Gx + d ,（2.3）从而，将向量*伏=0,1

31、,）代入4）的右边，可计算得到一个新的值一+1： xk+l =Gxk +d ,（2.3。如此反复地进行，便得到向量序列#，4=01,;问题I）如何将方程组（2.1）的形式转化为（2.3）?II）给定初始向量x,迭代（23）产生的向量序列是否收敛？ III）向量序列?的收敛性与初始向量x的选取是否相关？ IV）如果收敛，如何估计误差？3 .4.1 Jacobi 迭代与 Gauss-Seidel 迭代首先我们回答问题I）:由方程组（2.1）： Ax = b,使等式左端仅保留向量x,其他一概放到右端,即：将方程组的第i个方程（设为#0, i = l,2,）改成：3=- ai2x2- aiMxMainxn + d，便可将原方程组改写成如下等价形式：=- ai2X2- anXn + Pa22X2=-21XIa23X3a2nXn+2-八S,（2.4）- an2x2-+ A,Jacobi迭代：将/=(xf澧,弁/代入上式右端，便可(按顺序逐行)进行计算得到?+1=(xf+,xJ+1,-,xi+1)r,这便是Jacobi迭代：x*+l =(12X2_%Xj+A)llX2+=-CC2X23X3a2nXn+2)0 q1a22，(乙)