本节书摘来自华章出版社《数值分析（原书第2版）》一 书中的第2章，第2.6节，作者：（美）Timothy Sauer，更多章节内容可以访问云栖社区“华章计算机”公众号查看。
2.6 用于对称正定矩阵的方法
对称矩阵由于它们的特殊结构，和一般的矩阵相比，它们只有一半数量的独立元素，在线性方程组求解中占据一个有利的位置.这就提出了一个问题，形如LU的分解是否可以用一半的计算代价实现，并且仅仅使用一半的内存.对于对称正定矩阵，可以使用楚列斯基（Cholesky）分解实现这个目的.
对称正定矩阵还允许以一个非常不同的方式求解Ax=b，该方式并不依赖于矩阵分解.这种新方法，被称为共轭梯度方法，对于大规模的稀疏矩阵非常有用，该方法也属于迭代方法.
本节首先定义对称正定矩阵的概念，然后显示所有的对称正定矩阵A，使用楚列斯基分解方法，可分解为A=RTR，其中R是一个上三角矩阵.所以，问题Ax=b可以使用两步回代进行求解，和非对称情况下的LU分解类似.在本节的最后我们介绍共轭梯度方法以及预条件.
2.6.1 对称正定矩阵
定义2.12 如果AT=A，则n×n矩阵A是对称矩阵.如果对于所有向量x≠0，xTAx>0，则矩阵A是正定矩阵.
例2.26 证明矩阵A=22
25是对称正定矩阵.
显然A是对称矩阵.为了证明它是正定矩阵，我们使用定义:xTAx=［x1 x2］22
25x1
x2
=2x21+4x1x2+5x22
=2(x1+x2)2+3x22该式子总是非负，而且不可能为0，除非x2=0以及x1+x2=0，而这两个条件放在一起意味着x=0.
例2.27 证明对称矩阵A=24
45不是正定矩阵.
用配方法计算xTAx:
117xTAx=［x1 x2］24
45x1
x2
=2x21+8x1x2+5x22
=2(x21+4x1x2)+5x22
=2(x1+2x2)2-8x22+5x22
=2(x1+2x2)2-3x22例如，令x1=-2，x2=1，使得结果小于0，和正定矩阵的定义矛盾.
注意到对称正定矩阵是非奇异矩阵，因为不存在非零向量x满足Ax=0.此外，对于这类特殊的矩阵还有三个重要的性质.
性质1 如果n×n矩阵A是对称矩阵，则A是正定矩阵，当且仅当所有特征值是正数.
证明 定理A.5指出，一组单位特征向量是正交向量，并张出Rn空间.如果A是正定矩阵，对于非零向量v，Av=λv，则00.另一个方面，如果A的所有特征值是正数，则写出任何一个非0的向量x=c1v1+…+cnvn，其中vi是正交的单位向量，ci不全为0.则xTAx=(c1v1+…+cnvn)T(λ1c1v1+…+λncnvn)=λ1c21+…+λnc2n>0,因而A是正定矩阵.
在例2.26中A的特征值是6和1.例2.27中A的特征值近似是7.77和-0.77.
性质2 如果A是n×n对称正定矩阵，X是一个满秩n×m矩阵,n≥m，则XTAX是m×m对称正定矩阵.
证明 由于(XTAX)T=XTAX，所以矩阵对称.为了证明正定，考虑一个非零的m维向量v.注意到vT(XTAX)v=(Xv)TA(Xv)≥0，由于矩阵A是正定矩阵，仅当Xv=0时等号成立.由于X满秩，它的列向量线性无关，因而Xv=0意味着v=0.
定义2.13 方阵A的主子矩阵是一个方的子矩阵，其对角线元素是矩阵A的对角线元素.
性质3 任何对称正定矩阵的主子矩阵也是对称正定矩阵.
证明 参见习题12.
例如，如果a11a12a13a14
a21a22a23a24
a31a32a33a34
a41a42a43a44是对称正定矩阵，则a22a23
a32a33也是对称正定矩阵.118
2.6.2 楚列斯基分解
为了验证主要想法，我们从一个2×2的例子开始.在这个例子中包含对所有重要问题的讨论，从简单例子到一般规模的扩展仅仅需要额外的记录.
考虑对称正定矩阵ab
bc从对称正定矩阵的性质3我们知道a>0.并且我们知道矩阵A的值ac-b2是正的，这是由于矩阵的值为对应特征值的乘积，从性质1知道所有特征值都是正数，所以乘积也是正数.使用上三角矩阵将矩阵A写做A=RTR意味着ab
bc=a0
uvau
0v=aua
uau2+v2我们希望检查如上的分解是否可能.从我们对于矩阵的值的了解，左右两侧进行比较得到u=b/a以及v2=c-u2.注意，我们知道v2=c-(b/a)2=c-b2/a＞0.这也验证了v是一个实数，楚列斯基分解如下A=ab
bc=a0
bac-b2/aaba
0c-b2/a=RTR对于2×2对称正定矩阵成立.楚列斯基分解不唯一，很显然，我们同样可以选择c-b2/a负的平方根作为v.
下面的结论保证，相同的思想对于n×n情况也成立.
定理2.14（楚列斯基分解定理） 如果A是n×n对称正定矩阵，则存在上三角n×n矩阵R满足A=RTR.
证明 我们通过对于矩阵大小n使用归纳法构造R.n=2的情况前面已经解决.考虑矩阵A分割为A=a|bT
bC
其中，b是一个（n-1）维的向量，C是一个（n-1）×（n-1）子矩阵.我们将使用分块乘法（见附录A.2节）进行简化.如2×2的情况下一样，设u=b/a.令A1=C-uuT，定义可逆矩阵S=a|uT
0
I
0
119得到ST1|0 … 0
0
A1
0S=a|0 … 0
uI1|0 … 0
0
A1
0a|uT
0
I
0
=a|bT
buuT+A1
=A 注意到A1对称正定.这是由于，根据性质2，如下矩阵也是对称正定的1|0 … 0
0
A1
0=(ST)-1AS-1因而根据性质3，（n-1）×（n-1）主子矩阵A1也是对称正定的.通过归纳假设，A1=VTV，其中V是上三角矩阵.最后，定义上三角矩阵R=a|uT
0
V
0并检查RTR=a|0 … 0
uVT
a|uT
0
V
0=a|bT
buuT+VTV
=A完成证明.
显式构造的定理证明，对应于楚列斯基分解的标准算法.矩阵R从外向里构造.首先我们找到r11=a11，并令R的最上面一行的其他元素为uT=bT/r11.然后从（n-1）×（n-1）下主子矩阵中减去uuT，然后进行相同的步骤，填上矩阵的第2行.重复进行这些步骤，直到确定R的所有行.根据定理，在构造矩阵的每个阶段，新生成的主子矩阵是正定矩阵，因而根据性质3，矩阵的左上角元素是正的，可以进行平方根操作.该方法可以直接放到后面的算法里.我们使用“冒号”表达子矩阵.120
楚列斯基分解
结果中矩阵R是上三角矩阵，满足A=RTR.
例2.28 找出如下矩阵的楚列斯基分解:4-22
-22-4
2-411R的最上面一行是R11=a11=2，其他元素是R1,2∶3=［-2,2］/R11=［-1,1］:R=2|-1 1
从A的2×2下主子矩阵中减去外积uuT=-1
1-11得到|[-]
2-4
-411-|[-]
1-1
-11=|[-]
1-3
-310现在我们在2×2的子矩阵中重复这个过程，找到R22=1以及R23=-3/1=-3:R= 2[|]-1[-] 1
1[|]-3[-]
下面的1×1主子矩阵是10-（-3）（-3）=1，因而R33=1.矩阵A的楚列斯基因子如下R=2-11
01-3
001 对于对称正定矩阵A，求解Ax=b，和LU分解的方式相同.既然A=RTR是两个三角矩阵的乘积，我们需要求解下三角矩阵系统RTc=b和上三角系统Rx=c来求解x.
2.6.3 共轭梯度方法
共轭梯度方法的引入（Hestenes和Steifel，1952）将稀疏矩阵问题的迭代求解方法带入了一个新的时代.尽管这种方法没有很快流行起来，但是一旦引入有效的预条件，其他方式难以处理的问题现在都可以解决了.共轭梯度方法的成功应用带来了更多的技术进步以及迭代方法的一个新时代.121正交 本书中我们第一次应用正交是以一种间接的方式求解一个和正交没有明显关系的问题.共轭梯度方法在求解正定的n×n线性方程组中，持续地找出和消去n个正交的误差成分.通过使用两两正交的余项确定的方向，可以降低算法的复杂度.我们在第4章中将进一步讨论这个问题，在GMRES方法中的culminating方法，是共轭梯度法在非对称问题中的对应方法.
共轭梯度的思路依赖于内积思想的推广.因为（v，w）=（w，v），以及对于标量α和β，有（αv+βw，u）=α（v，u）+β（w，u），所以欧几里得内积(v,w)=vTw对称并对于输入v和w线性.欧几里得内积也是正定的，当v≠0时，（v，v）>0.
定义2.15 令A是对称正定的n×n矩阵.对于两个n维向量v和w，定义A内积(v,w)A=vTAw当（v，w）A=0时，向量v和w为A共轭.
注意到新的内积定义继承了矩阵A的对称、线性、正定的性质.由于A是对称矩阵，因而A内积也是对称的:(v,w)A=vTAw=(vTAw)T=wTAv=(w,v)A.如果A正定，则A内积也是线性、对称正定矩阵，并且(v,v)A=vTAv>0其中v≠0.
严格来讲，共轭梯度方法是一个直接方法，使用下面的有限步，就可以得到对称正定系统Ax=b的解x:
共轭梯度方法
首先对该迭代进行一种非正式的描述，然后是对定理2.16重要事实的证明.共轭梯度方法在每一步中更新三个向量.向量xk是第k步122时的近似解.向量rk表示近似解xk的余项.从定义看，r0显然是这样的余项，在迭代过程中，注意到Axk+1+rk+1=A(xk+αkdk)+rk-αkAdk
=Axk+rk因而对k进行归纳得到rk=b-Axk.最后，变量dk表示用于更新xk得到改进的xk+1时所使用的新的搜索方向.
该方法能够成功的原因在于所有的余项都和前面的余项正交.如果能做到所有的余项正交，该方法搜索所有的正交方向，在经过至多n步就可以得到余项为零的正确解.实现所有余项的正交的关键在于选择搜索方向dk，并使之两两共轭.共轭的概念推广了正交的概念，并据此在算法的名字中也包含“共轭”.
现在我们解释对于αk和βk的选择.从先前余项向量所张的空间中选择方向向量dk，这可以从伪码中的最后一行看到.为了保证下一余项向量和前面所有的余项向量都正交，需要精确选择αk使得新的余项rk+1和方向dk正交:xk+1=xk+αkdk
b-Axk+1=b-Axk-αkAdk
rk+1=rk-αkAdk
0=dTkrk+1=dTkrk-αkdTkAdk
αk=dTkrkdTkAdk这和算法中所写的αk不同，但是注意到dk-1和rk正交，我们有dk-rk=βk-1dk-1
rTkdk-rTkrk=0使得重写rTkdk=rTkrk成立.第二，选择系数βk，以保证dk两两A共轭:dk+1=rk+1+βkdk
0=dTkAdk+1=dTkArk+1+βkdTkAdk
βk=-dTkArk+1dTkAdk如下面的（2.47）所示，对于βk，可以重写为算法中的简单形式.
下面的定理2.16证实所有共轭梯度迭代所生成的rk彼此正交.由于它们都是n维向量，至多有n个rk可以两两正交，因而rn或者前面的rk必须为0，也就求解了Ax=b.因而至多n步后，共轭梯度方法得到一个解.理论上，该方法是直接方法，而不是迭代方法.
在转向共轭梯度方法保证成功的定理之前，使用精确算术做一个例子，这将会具有启发意义.
例2.29 使用共轭梯度方法求解22
25u
v=6
3
123使用上面的算法x0=0
0,r0=d0=6
3
α0=6
3T6
36
3T22
256
3=456·18+3·27=521
x1=0
0+5216
3=10/7
5/7
r1=6
3-52118
27=121/7
-2/7
β0=rT1r1rT0r0=144·5/4936+9=1649
d1=121/7
-2/7+16496
3=180/49
-120/49
α1=12/7
-24/7T12/7
-24/7180/49
-120/49T22
25180/49
-120/49=710
x2=10/7
5/7+710180/49
-120/49=4
-1
r2=121/7
-2/7-71022
25180/49
-120/49=0
0由于r2=b-Ax2=0，解是x2=［4,-1］.
定理2.16 令A为对称正定的n×n矩阵，b≠0是一个向量.在共轭梯度方法中，假设rk≠0，其中k（a） 后面Rn的三个子空间等价:〈x1,…,xk〉=〈r0,…,rk-1〉=〈d0,…,dk-1〉（b） 余项rk两两正交rTkrj=0,其中j（c） 方向dk两两A共轭dTkAdj=0,其中j证明 （a） 对于k=1，注意到〈x1〉=〈d0〉=〈r0〉，因为x0=0.由定义知xk=xk-1+αk-1dk-1.这意味着通过归纳，〈x1,…,xk〉=〈d0,…,dk-1〉.相似地，利用dk=rk+βk-1dk-1,表明〈r0,…,rk-1〉等于〈d0,…,dk-1〉.
对于（b）和（c），使用归纳.但k=0时不需要证明.假设（b）和（c）对于k成立，我们将证明（b）和（c）对于k+1成立.对rk+1的定义在左侧乘上rTj:rTjrk+1=rTjrk-rTkrkdTkAdkrTjAdk(2.46)如果j≤k-1，则使用归纳假设（b），可得到rTjrk=0.由于rj可以表达为d0，…，dj的124组合，从归纳假设（c）得到rTjAdk=0，并且（b）成立.另一方面，如果j=k，则由于使用归纳假设（c），dTkAdk=rTkAdk+βk-1dTk-1Adk=rTk，由（2.46）知rTkrk+1=0，这证明了（b）.
既然rTkrk+1=0，（2.46）指出,对于j=k+1，rTk+1rk+1rTkrk=-rTk+1AdkdTkAdk(2.47)这与在dk+1的定义的左侧乘上dTjA一起得到dTjAdk+1=dTjArk+1-rTk+1AdkdTkAdkdTjAdk(2.48)如果j=k，则利用矩阵A的对称性，从（2.48）知dTkAdk+1=0.如果j≤k-1，则Adj=(rj-rj+1)/αj（从rk+1的定义可得到）与rk+1正交，表明（2.48）右侧的第一项是零，使用归纳假设，第二项也是零，这就完成了（c）的论证.
在例2.29中，注意到r1和r0正交，这与定理2.16所声明的一致.这种正交性是共轭梯度法成功的关键：每个新的余项ri和所有前面的ri正交.如果一个ri是零，则Axi=b，xi就是解.否则，在n步循环后，rn和n个两两正交的向量r0，…，rn-1所张成的空间正交.而r0，…，rn-1这n个向量是Rn空间中的所有正交向量.因而rn必须是零向量，所以Axn=b.
共轭梯度法在某种方式上比高斯消去法简单.例如，写出的代码看起来更简单，其中不需要担心高斯消去法中的行交换.两种方法都是直接方法，在有限步骤后可以得到理论正确的解.这就带来了两个问题:为什么共轭梯度法比高斯消去法好，以及为什么共轭梯度法常常被看做是迭代方法？
为了回答这两个问题，首先计算操作的次数.在循环中走一轮需要一次矩阵向量乘法Adn-1以及一些内积计算.矩阵向量乘法本身在每一步中需要n2次的乘法（以及相同数量的加法），在所有n步计算后需要n3次的乘法.和高斯消去的n3/3的计算次数相比，计算次数乘上了3倍，代价更大.
如果A是稀疏矩阵，这个问题就不同了.假设对于n3/3次高斯消去操作，n大得难以再进行高斯消去.虽然高斯消去必须在所有步骤都完成才能得到解x，但是共轭梯度方法在每一步中都给出一个解xi.
后向误差和余项的欧几里得长度，随着每一步下降，因而至少对于这种度量方式，Axi在每一步变得和b越来越接近.因而通过检测ri，可能得到一个足够好的xi，并不必做完n步.在这种情况下，共轭梯度方法和迭代方法没有区别.
当A是一个病态矩阵时，由于对舍入误差累计的敏感性，这种方式的优势很快就消失了.实际上，共轭梯度法在病态矩阵上性能比部分主元的高斯消去法还要差.在当前，这种问题通过预条件得到缓解，主要是将问题转化为良态矩阵系统，然后再实施共轭梯度法.我们将在下一节中讨论预条件的共轭梯度法.
共轭梯度法的名字来自于共轭梯度法实际进行的操作:在n维空间的二次抛物125面的斜面上滑下来.名字中“梯度”表示通过微积分寻找最速下降的方向，“共轭”不是表示单个步骤和其他步正交，而是表示余项之间的正交性质.该方法的几何细节以及动因很有趣.对此，Hestenes和Steifel［1952］最初的论文给出一个完整的描述.
例2.30 对系统（2.45）使用共轭梯度法，n=100000.
共轭梯度法20步后，计算解x和精确解（1，…，1）之间的差异用向量的无穷范数度量，该误差小于10-9.在PC上所需要的运行时间小于1秒.
2.6.4 预条件
通过使用预条件技术，可以使得诸如共轭梯度方法的迭代方法收敛速度加快.迭代方法的收敛率通常直接或者间接依赖于系数矩阵A的条件数.预条件方法的思想是降低问题中的条件数.
n×n的线性方程组Ax=b的预条件形式是M-1Ax=M-1b其中M是可逆的n×n矩阵，称为预条件子.我们所要做的就是在方程两侧左乘上该矩阵.预条件矩阵试图对矩阵A逆转，从而可以有效地降低问题的条件数.概念上来讲，它试图同时做两件事:矩阵M应该（1）和矩阵A足够接近，（2）容易求逆.这两个目的常常彼此对立.
和A最接近的矩阵就是A自身.使用M=A会把问题的条件数变为1，但是一般A并不容易进行逆转，或者我们不想使用复杂的方法.最容易求逆的矩阵是单位阵M=I，但是它并不能降低条件数.完美的预条件矩阵应该居于两种极端的中间，并同时具备二者好的性质.
一种特别简单的方式是雅可比预条件子M=D，其中D是A的对角线矩阵.D的逆矩阵是D的元素的倒数.例如在一个严格对角占优矩阵中，雅可比预条件子和A相似，同时非常容易求逆.注意到，由2.6.1节中的性质3，对称正定矩阵的每个对角线元素都严格为正.因而计算倒数不是问题.
当A是n×n对称正定矩阵，我们将选择对称正定矩阵M作为预条件子.回忆2.6.3节中的M内积（v，w）M=vTMw.现在很容易描述预条件的共轭梯度方法:使用预条件方程M-1Ax=M-1b替换Ax=b.使用(v，w)M替换欧几里得内积.使用原始的共轭梯度法的原因仍然成立，这是由于矩阵M-1A相对于新的内积仍然是对称正定矩阵.
例如，(M-1Av,w)M=vTAM-1Mw=vTAw=vTMM-1Aw=(v,M-1Aw)M把2.6.3节中的算法转化为预条件版本，令zk=M-1b-M-1Axk=M-1rk是预条件系统的余项.则126αk=(zk,zk)M(dk,M-1Adk)M
xk+1=xk+αdk
zk+1=zk-αM-1Adk
βk=(zk+1,zk+1)M(zk,zk)M
dk+1=zk+1+βkdk乘上M可以进行消减(zk,zk)M=zTkMzk=zTkrk
(dk,M-1Adk)M=dTkAdk
(zk+1,zk+1)M=zTk+1Mzk+1=zTk+1rk+1通过这些简化，预条件版本的伪代码如下.
预条件共轭梯度法
k步迭代后方程Ax=b的近似解是xk.注意到并没有和M-1进行显式的相乘.由于相对简单的M，该乘法可以被适当的回代所取代.
雅可比预条件子是当前还在不断扩展和增长的可能选择的库中最简单的一个.我们将描述另外一个系列的例子，并且读者可以参考相关文献得到更加复杂的选择.
对称连续过松弛（SSOR）中，预条件子定义如下：M=(D+ωL)D-1(D+ωU)其中A=L+D+U被分为下三角部分、对角线以及上三角部分.如同SOR方法，ω是一个在0和2之间的常数.在特例中ω=1，这被称为高斯塞德尔预条件子.
如果预条件子矩阵难以求逆，一般就没什么用.注意到SSOR预条件子定义为下三角矩阵和上三角矩阵的乘积M=(I+ωLD-1)(D+ωU)，因而方程z=M-1v可以通过两次回代求解.(I+ωLD-1)c=v
(D+ωU)z=c
127对于稀疏矩阵，两次回代所花的时间和非零元素的个数成比例.换句话讲，乘上M-1并不比乘上M的计算复杂度高.
例2.31 令A表示矩阵，对角线元素Aii=i，其中i=1，…，n，Ai，i+10=Ai+10，i=cosi，i=1，…，n-10，所有其他的元素都是0.令x是一个包含n个1的向量，定义b=Ax.对于n=500，使用共轭梯度法以三种方式求解Ax=b:不使用预条件，使用雅可比预条件子，使用高斯塞德尔预条件子.
在MATLAB中矩阵定义为
图2.4显示三种不同的方式.即使对于例子中简单的矩阵，没有使用预条件，共轭梯度法收敛很慢.非常容易使用的雅可比预条件子，对求解有了显著的改善，而使用高斯塞德尔预条件子，仅仅需要大约10步就可以得到机器精度.
图2.4 对于例2.31求解使用预条件的共轭梯度法的效率.相对迭代步数，画出对应误差.
圆形:没有预条件子.方形:雅可比预条件子.菱形:高斯塞德尔预条件子2.6节习题
1.证明下面的矩阵是对称正定矩阵，把xTAx表示为平方和.
(a) 10
03 (b) 13
310 (c)100
020
003
2.通过发现满足xTAx＜0的向量x，证明下面的矩阵不是对称正定矩阵.
(a) 10
0-3 (b) 12
22
(c) 1-1
-10 (d) 100
0-20
003
128
3.使用楚列斯基分解把习题1中的矩阵表示为A=RTR.
4.证明楚列斯基分解过程对于习题2中的矩阵失效.
5.找出下面每一个矩阵的楚列斯基分解A=RTR.
(a) 12
28 (b) 4-2
-25/4 (c) 255
526 (d) 1-2
-25
6.找出下面每一个矩阵的楚列斯基分解A=RTR.
(a) 4-20
-22-3
0-310(b) 120
252
025(c) 111
122
123(d) 1-1-1
-121
-112
7.通过找出矩阵A的楚列斯基分解，以及后面的两次回代求解方程组.
(a) 1-1
-15x1
x2=3
-7(b) 4-2
-210x1
x2=10
4
8.通过找出矩阵A的楚列斯基分解，以及后面的两次回代求解方程组.
(a) 40-2
011
-213x1
x2
x3=4
2
0(b) 4-20
-22-1
0-15x1
x2
x3=0
3
-7
9.证明:如果d>4，矩阵A=12
2d正定.
10.找出所有数字d，使得A=1-2
-2d正定.
11.找出所有数字d，使得A=1-10
-121
01d正定.
12.证明:对称正定矩阵的主子矩阵对称，并且正定.（提示:考虑一个合适的X并使用性质2.）
13.通过手工进行共轭梯度法求解问题.
(a) 12
25u
v=1
1(b) 12
25u
v=1
3
14.通过手工进行共轭梯度法求解问题.
(a) 1-1
-12u
v=0
1(b) 41
14u
v=-3
3
15.在一般的标量情况下，使用共轭梯度法求解Ax=b，其中A是一个1×1矩阵.找出α1、x1，并证实r1=0以及Ax1=b.129
2.6节编程问题
1.写出共轭梯度法的MATLAB版本，并使用该程序求解如下系统：
(a) 10
02u
v=2
4(b) 12
25u
v=1
1
2.使用MATLAB版本的共轭梯度法，求解如下问题:
(a) 1-10
-121
012u
v
w=0
2
3(b) 1-10
-121
015u
v
w=3
-3
4
3.使用共轭梯度法求解系统Hx=b，其中H是n×n希尔伯特矩阵，b是全为1的向量，（a） n=4，（b） n=8.
4.使用共轭梯度法求解（2.45）的稀疏问题.（a） n=6，（b） n=12.
5.使用共轭梯度法求解（2.45），n=100，1000以及10000.报告最后余项的规模，以及所需的迭代步数.
6.令A是n×n矩阵，n=1000，元素A(i,i)=i,A(i,i+1)=A(i+1,i)=1/2,A(i,i+2)=A(i+2,i)=1/2.（a） 使用spy（A），打印非零结构.（b） 令xe是n个1组成的向量.令b=Axe，并使用没有预条件子的共轭梯度法和雅可比预条件子，以及高斯塞德尔预条件子求解.在图中比较三种方法相对于步数的运行结果误差.
7.令n=1000.使用编程问题6中的n×n矩阵A，并加上非零元素A(i,2i)=A(2i,i)=1/2,1≤i≤n/2.完成问题6中的步骤（a）和（b）.
8.令n=500，令A是n×n矩阵，对于所有i，元素A(i,i)=2,A(i,i+2)=A(i+2,i)=1/2,A(i,i+4)=A(i+4,i)=1/2，并且A(500,i)=A(i,500)=-0.1，1≤i≤495.完成问题6中的（a） 和（b）.
9.令A为编程问题8中的矩阵，但是对角线元素替换为A(i,i)=3i，完成问题8的（a）和（b）部分.
10.令C是195×195矩阵块，其中对于所有的i，C(i,i)=2,C(i,i+3)=C(i+3,i)=0.1，C(i,i+39)=C(i+39,i)=1/2,C(i,i+42)=C(i+42,i)=1/2.定义A是n×n矩阵，n=780，矩阵的对角线由4个对角矩阵块C组成，并在上三角和下三角矩阵中放置矩阵块12C.完成问题6的（a）和（b）步骤，并求解Ax=b.
