精通MATLAB科学计算（第3版）(王正林)09-3r

《精通MATLAB科学计算（第3版）(王正林)09-3r》由会员分享，可在线阅读，更多相关《精通MATLAB科学计算（第3版）(王正林)09-3r（39页珍藏版）》请在金锄头文库上搜索。

1、第 章线性方程组求解在自然科学和工程技术中， 很多问题可以归结为求解线性方程组。采用MATLAB ,不仅可以利用其提供的相关函数直接解决一些简单的线性方程组， 而且可以通过简洁的编程来解决一些复杂的线性方程组。通过本章的介绍， 读者既能应用MATLAB中相应的函数求解线性方程组， 又能通过编 程 ， 灵活使用迭代法和其他的特殊解法来求解线性方程组。9 .1 求逆法MATLAB中求解线性方程最直接的方法是矩阵求逆法， 它适用于系数矩阵的数据无规律且系数矩阵的阶数比较小的情况。如果系数矩阵的阶数太大的话， 系数矩阵求逆就需要花费很长时间。在 MATLAB中用求逆法求解， 线性方程可以直接用左除法“

2、 ” 求 解 ， 也可以用求逆函数 inv()求解。【 例 9-1】 左除法和求逆法求解线性方程组实例。用左除法和求逆法求解如下线性方程组，X + 2%2 + 3x3 = 1 x=inv(A )*b输出计算结果为： X =0.3333-1.66671.3333由计算结果可知, 方程组的解为 的42户3 = 03333,-1.6667,1.3333。9.2分解法矩阵分解法是指根据一定的原理用某种算法将系数矩阵分解成若干个矩阵的代数运算 ， 常用的分解是乘积分解， 而分解后的新矩阵一般是某种特殊矩阵。常用的分解有LU分解、Q R分解、Cholesky分解、Schur分解、Hessenberg分解和

3、奇异分解等。9.2.1 LU分解法矩阵的LU分解也叫三角分解， 还有的书称之为Doolittle分 解 ， 它是将矩阵分解为一个单位下三角矩阵与上三角矩阵的乘积。只要矩阵非奇异， 这种分解总是可以进行的。MATLAB提供了 lu 函数来求矩阵A 的 LU分 解 ， 常用格式为：也, S= lu(Z )： 产生一个上三角矩阵U 和一个变换形式的下三角矩阵Z ( 进行了行交换 ) ， 使得A=LU ；,(/,P=lu(/1)： 产生一个上三角矩阵U 和一个下三角矩阵乙及置换矩阵P , 使得PA=LU ；分解以后， 方 程 组 = 的解可写为x=LA(ZA ) 或x=(A(ZAP* ) 。【 例9-

4、21 LU分解法求解线性方程组实例。用 LU分解法求解以下线性方程组，1.5xi + 3x2 - 0.83 + 414 = 42xi + 9% 3 + 10%4 = 0 lx + 4.8x2 0.6x3 + = 114%| + 12.3%2 4M + 5X4 = 2 解 ：用 LU分解法求解， 在 MATLAB命令窗口中输入求解程序： A =1.5 3 -0 .8 4;2 0 9 1 0 ;-7 4.8 -0 .6 1;14 12.3 -4 5 ; b=4 0 1 - 2 ; Lz U =lu (A);1 187精通MATLAB科学计算( 第 2 忡 - - - - - - - - - - -

5、 - - - - - - - - - - - - - - - - - x= U (Lb)输出计算结果为： x =-0.3721-0.8291-1.53361.4547由计算结果可知, 方程组的解 为 即 X2, X3, X4 = -0.3721, -0.8291, -1.5336, 1.4547。9.2.2 QR分解法矩阵的Q R 分解就是把矩阵分解为一个正交矩阵和一个上三角矩阵的乘积。 M A T L A B中对矩阵A进行Q R 分解的函数调用格式如下：Q, K= qr(/) ： 产生一个正交矩阵。和一个上三角矩阵A , 使得；Q, K, E= qr(/) ： 产生一个正交矩阵Q 、一个上三

6、角矩阵R以及一个置换矩阵E , 使得AE=QR ；实现分解以后， 方程组4r=6 的解可写为x= K(Q功) 或x= E(K(QS) ) 。【 例 9-31 Q R 分解法求解线性方程组实例。用 Q R 分解法求解以下线性方程组，x + 0.5x2 + 0.3333%3 + 0.25%4 = 10.5%1 + 0.3333%2 + 0.25x3 + 0.2%4 = 20.3333x1 + 0.25x2 + 0.2%3 + 0.1667x4 = 20.25xi + 0.2必 + 0.1667x3 + 0.1429x4 =1 。解 ：用 Q R 分解法求解， 在 M A T L A B 命令窗口中

7、输入求解程序： A=1.00000.50000.33330.25000.50000.33330.25000.20000.33330.25000.20000.16670.2500;0.2000;0.1667;0.1429 ; b = l 2 2 1 , ; Q ,R=qr (A) x=R(Q b)输出。、A的计算结果为：Q =-0.83810.5226-0.1540-0.0263-0.4191-0.44170.72780.3157-0.2794-0.5288-0.1395-0.7892-0.2095-0.5021-0.65360.5261-1.1932-0.6705-0.4749-0.36980

8、-0.1185-0.1257-0 .1175188 第 9 章 线性方程组求解0 0 -0.0062 -0.0096000 0.0002输出解的计算结果为：x =1.0e+003 *0.1160-1.44003.6000-2.3800由计算结果可知, 方程组的解为 的广2,孙孙 = 口的, -1440,3600,-2380。9.2.3 Cholesky 分解法当系数矩阵N正定且对称时， 它 能 分 解 为 以 下 的 形 式 : 其 中 K 为上三角矩阵，肝 为R的转置， 是下三角矩阵， 这种分解称为Cholesky分解。MATLAB中与Cholesky分解有关的函数如下：R=chol(Z):

9、 对 进行 Cholesky 分 解 ， 使得 A = RXR ;K,p=chol(4): 对 A 进行 Cholesky 分 解 ， 使得 A=RrRa如果4 对称正定， 则 p=0 ； 否 则 ， p 为一正整数。如果4 满 秩 ， / ? 是为2 -1 阶的上三角矩阵， 且有相/ ? = 4 1 :仍- 1) , l:(p -l)； 实现分解后， 方 程 组 的 解 写 成 上 大 ( 如 财的形式。【 例 9-4 Cholesky分解法求解线性方程组实例。用 Cholesky分解法求解以下线性方程组，9xi - 36切 + 3013 = 1* -36xi + 1922 一 180工 3

10、 = 130%1-180%2+180%3=1 0解：用 Cholesky分解法求解， 在 MATLAB命令窗口中输入求解程序： A=9 -36 30;-36 192 -180;30 -180 180); b=ones (3,1); R=chol(A) x=R(Rzb)输出/ ? 的计算结果为：R =3.0000 -12.0000 10.00000 6.9282 -8.66030 0 2.2361输出解的计算结果为： 189精通MATLAB科学计算( 第2版i -X =1.83331.08330.7833由计算结果可知， 方程组的解为 莺/ 2,必 = 1.8333,1.0833,0.7833

11、。起 合 在 MATLAB中 ， 当N正定但不对称时， 用 K=chol(N)也能得出一个矩阵，MATLAB并不报错， 但是只要验证一下Rt*R就会看出R *R并不等于A , 因此用公式x=R(*M)就会得出错误的结果。190 第9章 线性方程组求解9 .2 .4其他分解法MATLAB中还提供了矩阵的奇异值分解、Hessenberg分解和Schur分解等函数， 可用来求解线性方程组， 如表9-1所示。表9 - 1其他的矩阵分解函数函数及常用语法说 明USM =svd(j)奇异值分解， 使得P.H =hess()Hessenberg分解. 使得A = P*H*/， 其中P为酉矩阵U,T =sch

12、ur(/l)Schui分解， 使得力。 ， 其中U为正交矩阵，T为Schur矩阵【 例 9-5】 奇异值分解法求解线性方程组实例。用奇异值分解法求解以下线性方程组 ，X) + 0.5x2 + 0.3333%3 + 0.25x4 + 0.2x5 = 10.5%1 +0.3333x2 +0.25%3 +0.2x4 +0.1667x5 = 0UzSrV=svd(A)b = l 0 1 0 1 zrx=V *inv(S)*UI *b输出分解的I/、S、H的计算结果为：U =-0.76790.6019-0.21420.04720.0062-0.4458-0.27590.7241-0.4327-0.116

13、7-0 .3 2 1 6-0.42490.12050.66740.5062-0.2534-0.4439-0.30960.2330-0.7672-0.2098-0.4290-0.5652-0 .55760.3762S=1.5671000000.2085000000.011400000 0.000300000 0.,0000V=-0.76790.6019-0.21420.04720.0062 191精通MATLAB科学计算（ 第2蚪-0.4458-0.3216-0.2534-0.2098-0.2759-0.4249-0.4439-0.42900.72410.1205-0.3096-0.5652-0

14、.43270.66740.2330-0.5576-0.11670.5062-0.76720.3762输出的计算结果为：x =1.0e+004 *0.0865-1.24674.4659-5.84042.5426由计算结果可知， 方程组的解为，x1,X2,x3,X4,x5 = 104*0.0865,-1.2467,4.4659,-5.8404,2.5426。选W奇异值分解很有用， 将系数矩阵进行奇异值分解以后，A的奇异值都在S的对角线上， 这样就可以粗略估计N的条件数， 看看N是否是病态矩阵，以此决定相应的求解方法。如果系数矩阵对称， 那么奇异值分解后的。和/ 都是正交矩阵， 而 S 是对角矩阵，

15、 这样求解方程就十分方便。【 例 9-6 Hessenberg分解法求解线性方程组实例。用 Hessenberg分解法求解以下线性方程组，闲+工2 +工3 +工4 = 1X + 2X2 + 3汹 + 414 = 4X + 3X2 + 6x3 + 1014 = 7X + 4X2 +1013 + 204 = -2 o解 ：用 Hessenberg分解法求解， 在 MATLAB命令窗口中输入求解程序： A = 1 1 1 1;1 2 3 4;1 3 6 10;1 4 10 20; P,H=hess(A) b=l 4 7 -21; x=P*inv(H)*P1*b输出分解的P 、的计算结果为：P =0.

16、7754 -0.6247 -0.0925 0-0.6092 -0.7015 -0.3698 0192 第9章 线性方程组求解0.166200.34310-0.92450 1.0, 00000.2147-0.265400-0.26541.08441.0802001.08027.7009-10.816700-10.816720.0000输出的计算结果为：10.0000-33.000036.0000-12.0000由计算结果可知, 方程组的解为,x i,x2, X 3 = 10 ,-3 3,36 ,-12 o起 合 Hessenberg矩阵指的是第一子对角线以下的元素为0的矩阵。 Hessenber

17、g分解法分解后产生正交矩阵和Hessenberg矩 阵 ， 正交矩阵的逆就是其转置； 如果系数矩阵对称的话， 产生的Hessenberg矩阵就是三对角矩 阵 ， 这种矩阵的求逆速度也很快， 因此用这种方法求解线性方程组也很快。【 例9-7】 Schur分解法求解线性方程组实例。Schur分解法求解以下线性方程组，汨+ 工2 + 工3 + 工4 = 1X + 2X2 + 3冷 + 414 = 4X 1 + 3x2 + 6x3 + 10x4 = 7X + 42 +10工3 + 20%4 = -2 o解 ：用Schur分解法求解， 在M A T L A B命令窗口中输入求解程序: A = 1 11

18、1；1 23 4;1 36 10;1 410 20; b= l 4 7 -2 UzT= schur (A )x= U* inv(T) * Uf* b输出分解的 U、7的计算结果为：U =0.3087-0.78730.53040.0602-0.72310.16320.64030.20120.59460.53210.39180.4581-0.1684-0.2654-0.39390.8638T =0.038000000.4538001 1939.3精 通MATLAB科 学 计 算 （ 第2版i -0 0 2.2034 0000 26.3047输出的计算结果为：x =10.0000-33.000036

19、.0000-12.0000由计算结果可知, 方程组的解为XI,X2,X3,X4 = 10,-33,36,-12。通 多 Schur矩阵是上三角矩阵， 且其对角元素为被分解矩阵的特征值。 Schur分解法在分解后产生正交矩阵和Schur矩 阵 ， 如果系数矩阵对称的话，产生的Schur矩阵就是对角阵， 显然对角阵的逆矩阵非常容易得到， 因此用这种方法求解线性方程组也很快。迭代法在实际应用中（ 例如在运筹学、图论等领域中） , 往往会出现这样一种情况： 系数矩阵阶数很高（ 例 如 ICT ） , 但系数矩阵含零元素相对较多， 如果用前面的几种分解法去求解的 话 ， 非零元素反而增多了， 这就有点得

20、不偿失了。本节介绍另一类常用的求解方法迭代法。迭代法是将求一组解转换为求一个近似解序列的过程， 并用最终的近似解来逼近真实解。迭代法需要考虑以下3 个重要的问题。（ 1 ） 迭代的初始值考虑初始值的选取是否有范围限制， 不同的初始值对最终的迭代结果是否有影响。（ 2 ） 迭代算法考虑怎样由当前的迭代结果得出下次迭代的初始值； 由于不同的算法会带来不同的误差 ， 所以应考虑算法使误差放大还是减少。（ 3 ） 迭代的收敛性考虑迭代是否收敛， 收敛的过程是快还是慢。以上问题的具体分析请参考数值分析的教材， 下面讲述采用MATLAB编程实现的几种常见迭代法。9 .3 .1逐次逼近法对于n阶线性方程组，

21、Ax=b的系数矩阵N （ 假设N为非奇异） 进行如下分解A=Q-C194 第 9 章 线性方程组求解其中。为非奇异。则方程组可变换为x=Bx+r其中B=Q C , r=Q： b ,那么迭代过程可以写成xk+l=Bxk+r这种迭代过程称为逐次逼近法。取不同的。和 C , 就得到不同的迭代算法。通常情况下 ， 逐次逼近法收敛的充分条件是迭代矩阵的谱半径小于lo9 .3 .2里查森迭代法里查森迭代法可以说是最简单的迭代法了， 它采用的迭代公式如下：xk+i=(I-A)xk+b在 MATLAB中编程实现的里查森迭代法函数为： richasono功 能 ： 用 里 查 森 迭 代 法 求 线 性 方 程

22、 组 的 解 。调用格式： x,n=richason(, b,xO, eps, M)其 中 ，A为线性方程组的系数矩阵；分为线性方程组中的常数向量；xO为迭代初始向量；eps为解的精度控制( 此参数可选) ；M 为迭代步数控制( 此参数可选) ；x 为线性方程组的解；为求出所需精度的解实际的迭代步数。理查森迭代法的MATLAB程序代码如下：function x,n=richason(A,b,xO,eps,M)% 采用里查森迭代法求线性方程组Ax=b的解W 线性方程组的系数矩阵：A告线性方程组中的常数向量：b% 迭代初始向量：xO3解的精度控制：eps为迭代步数控制：M年线性方程组的解：x若求出

23、所需精度的解实际的迭代步数：nif(nargin = 3)eps = 1. Oe-6; %eps表示迭代精度1 195精通MATLAB科学计算( 第2忡 - - - - - - - - - - - - - - - - - - - - - - - - - - - - -M = 200; *M表示迭代步数的限制值elseif(nargin = 4)M = 200;endI =eye(size(A);xl=x0;x=(I-A)*x0+b;n=l;告迭代过程while(norm(x-xl)eps)xl=x;x=(I-A)*xl+b;n = n + 1; 会n为最终求出解时的迭代步数if(n=M)dis

24、p (Warning:迭代次数太多， 可能不收敛！ , ) ;return;endend【 例 98 】 理查森迭代法求解线性方程组实例。用理查森迭代法求解以下线性方程组 ， 其中取初始值为 0,0,0 。1.0170汨 一 0.0092M-0.0095x3 = 1, -0.0092X + 0.9903x2 + 0.0136x3 = 0-0.0095xi +0.0136x2 + 0.9898x3 =1 0解 ：用理查森迭代法求解， 在 MATLAB命令窗口中输入求解程序： A = 1.0170 -0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136

25、 0.9898; b=l 0 1; x0 = 0 0 0/_; x,n=richason(A,b,xO)输出的计算结果为：X =0.9738-0.00471.0010输出的迭代次数为：n =g身可 见 ， 经过5 步迭代， 理查森迭代法求出了方程组的解为：196 第9章 线性方程组求解阳,xz,X3 = 0.9738； 6.0047,1.0010 l巧， 巧， 演 =0.9738, -0 .0 0 4 7 , 1 .0 0 1 。对上述迭代计算结果进行验证， 在 MATLAB命令窗口中输入如下程序： A*x输出结果为：ans =1.00000.00001.0000经过验证， 也确定了解的正确性

26、。 一般情况下， 理查森迭代法难以收敛， 因为要使/ - /的谱半径小于1 , 通常要求力是严格对角占优的矩阵。 197精 通MATLAB科 学 计 算 ( 第2版i -9.3.3 Jacobi 迭代法如果系数矩阵”的主对角元全不为0 , 在上节4 的分解中取Q=DC=D-A其中O 是由/ 的主对角元素组成的对角阵， 则有B=I-D A , r=D b , 迭代公式为：x=(I-D Axk+D b这种迭代方法称为Jacobi迭代法。在 MATLAB中编程实现的Jacobi迭代法函数为： jacobio功 能 ： 用 Jacobi迭代法求线性方程组4r 的解调用格式： x,w= jacobi (

27、A, h,xO, eps, varargin)其 中 ，A为线性方程组的系数矩阵；8 为线性方程组中的常数向量；xO为迭代初始向量；eps为解的精度控制( 此参数可选) ；varargin为迭代步数控制( 此参数可选) ；x 为线性方程组的解；为求出所需精度的解实际的迭代步数。Jacobi迭代法的MATLAB程序代码如下：function x,n=jacobi(A,b,xO,eps,varargin)为采用Jacobi迭代法求线性方程组Ax=b的解电线性方程组的系数矩阵：A年线性方程组中的常数向量：b卷迭代初始向量：xO优解的精度控制：eps带迭代步数控制：varargin为线性方程组的解：x

28、生求出所需精度的解实际的迭代步数：nif nargin=3eps= 1. Oe-6;M = 200;elseif nargin第9章 线性方程组求解M = varargin1;endD=diag (diag (A); 务求A的对角矩阵L=-tril (Az-1) ; % 求A的下三角阵U=-triu (A, 1); 彩求A的上三角矩阵B=D(L+U);f=Db;x=B*xO+f;n=l; % 迭代次数多迭代过程while norm(x-xO)=epsxO=x;x =B*xO+f;n=n+l;if(n=M)disp (，Warning:迭代次数太多， 可能不收敛！,) ；return;enden

29、d【 例 9-9 Jacobi迭代法求解线性方程组实例。用 Jacobi迭代法求解以下线性方程组 ， 其中取初始值为1,1,1。0.9889xj -0.0005x2 -0.0002x3 = 1 x0 = ones (3Z1); xz n=jacobi(A,b,xO)输出的计算结果为：X =1.0114-0.00311.0062输出的迭代次数为：n =4Y 199精通MATLAB科学计算（ 第2版i -可 见 ， 经过4 步迭代， Jacobi法求出了方程组的解为：XI,X2,X3 = 1,0114,-0.0031,1.0062。 Jacobi迭代法对任意初始向量M） 收敛的充要条件为B的谱半径

30、小于1 , 即 8的所有特征值的绝对值都小于lo200 第 9 章 线性方程组求解9.3.4 Gauss-Seidel 迭代法仍然设系数矩阵N的主对角元全不为零， 如果对N作以下分解：A=(D-L)-U其中。的意义同Jacobi迭代法， ，为下三角矩阵， 。为上三角矩阵， 它的迭代公式为:xk+i=(D-)- Uxk+(D-L)-b在 MATLAB中编程实现的Gauss-Seidel迭代法函数为： gauseidelo功 能 ： 用 Gauss-Seidel迭 代 法 求 线 性 方 程 组 的 解 。调用格式：gauseidel (4 4 的, eps,防 。其 中 ， A 为线性方程组的系

31、数矩阵；b 为线性方程组中的常数向量；xO为迭代初始向量；eps为解的精度控制( 此参数可选) ；M 为迭代步数控制( 此参数可选) ；x 为线性方程组的解；为求出所需精度的解实际的迭代步数。Gauss-Seidel的 MATLAB程序代码如下：function x,n=gauseidel(A,b,xO,eps,M)会采用Gauss-Seidel迭代法求线性方程组Ax=b的解年线性方程组的系数矩阵：A% 线性方程组中的常数向量：b他迭代初始向量：xO软解的精度控制：eps多迭代步数控制：M会线性方程组的解：x务求出所需精度的解实际的迭代步数：nif nargin=3eps= 1. Oe-6;M

32、 = 200;elseif nargin = 4M = 200;elseif nargin=epsxO=x;x=G*xO+f;n=n+l;if(n=M)disp (Warning:迭代次数太多， 可能不收敛！， ) ；return;endend【 例 9-10】 GaussSeidel迭代法求解线性方程组实例。用 Gauss Seidel迭代法求解以下线性方程组， 其中取初始值为1,1,1。1.4449xi +0.7948x2 +0.880 卜3 =1, 0.6946XI +1.9568x2 +0.1730x3 =00.6213为 + 0.5226x2 + 1.9797x3 = 1 o解 ：用

33、 Gauss-Seidel迭代法求解， 在 MATLAB命令窗口中输入求解程序： A = 1.4449 0.7948 0.8801;0.6946 1.9568 0.1730;0.6213 0.5226 1.9797; b=l 0 1已 ； x0=zeros (3,1); x,n=gauseidel(A,b,xO)输出的计算结果为：X =0.5929-0.24430.3835输出的迭代次数为：n =n.a可 见 ， 经 过 11步迭代， Gauss-Seidel迭代法求出方程组的解为：Xi,x2,x3 = 0.5929,- 0.2443,0.3835。9 .3 .5超松弛迭代法202 - - -

34、 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 第 9章 线性方程组求解如果对系数矩阵”作以下分解：A=(D-a)L)-(l-a)D+a)U)其中。、L 、。的意义与Gauss-Seidel迭代法相同。 是一个事先选好的常数， 称为松弛因子， 当 时 叫 超 松 弛 法 ， 也叫SOR迭 代 法 ；当。1时叫低松弛法。其迭代公式为：x*+1 = ( -尸 ( 1 - co)D+a)UXk+(o(D-(oL)b关于超松弛迭代法的

35、收敛性有如下结论：若系数矩阵/ 对称正定， 当 0 。 2 , SOR迭代法收敛。在 MATLAB中编程实现的超松弛迭代法函数为： SORo功 能 ： 用 超 松 弛 迭 代 法 求 线 性 方 程 组 的 解 。调用格式： x, n=SOR (A, b, xO, w, eps, M)其 中 ，A为线性方程组的系数矩阵；6 为线性方程组中的常数向量；xO为迭代初始向量；w 为松弛因子；eps为解的精度控制( 此参数可选) ；M 为迭代步数控制( 此参数可选) ；x 为线性方程组的解；n为求出所需精度的解实际的迭代步数。超松弛迭代法的MATLAB程序代码如下：function x,n=SOR(A

36、,b,xO,w,eps,M)带采用超松弛迭代法求线性方程组Ax=b的解传线性方程组的系数矩阵： A省线性方程组中的常数向量： b与迭代初始向量： xO% 松弛因子： w为解的精度控制： eps超迭代步数控制： M% 线性方程组的解： x省求出所需精度的解实际的迭代步数： nif nargin=4eps= 1.Oe-6;M = 200;elseif nargin 4error 203精通MATLAB科学计算( 第2版i -returnelseif nargin =5M = 200;endif (w=2) 与 收敛条件要求error;return;endD=diag (diag (A) ) ; %

37、 求A的对角矩阵L=-tril (Az-1); 考求A的下三角矩阵U=-triu(Azl) ; % 求A的上三角矩阵B=inv(D-L*w)* (1-w)*D+w*U);f=w*inv(D-L*w)*b;x=B*xO+f;n=l; 恭迭代次数带迭代过程while norm(x-xO)=epsx0=x;x =B*xO+f;n=n+l;if(n=M)disp (1 Warning:迭代次数太多， 可能不收敛！， ) ；return;endend【 例 9-11 超松弛迭代法求解线性方程组实例。用超松弛迭代法求解以下线性方程组 ， 其中取初始值为0,0,0。4xi + 3x2 = 24 xz n=S

38、OR(Azb,xOz1.25)输出的计算结果为：X =3.00004.0000-5.0000输出的迭代次数为：204 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 第 9章 线性方程组求解n =14可 见 ， 经 过 14步迭代， SOR迭代法求出了方程组的解为R,X2,X3 = 3,4,-5。超松弛迭代法还有一种改进形式， 叫做对称逐次超松弛迭代法(SSOR), 它采用的是两步迭代公式：(D-a)L) x*+i/

39、2=0 (fJxk+b)+( 1 -a)Dxi；(D-coU) Xk+i=(o (Lxk+b)+( 1 -v)Z)x*.+i/2在 MATLAB中编程实现的对称逐次超松弛迭代法函数为： SSORo功 能 ： 对 称 逐 次 超 松 弛 迭 代 法 求 线 性 方 程 组 的 解 。调用格式： x, =SSOR (4 b, xO, iv, eps, M)o其 中 ，A为线性方程组的系数矩阵；b 为线性方程组中的常数向量；xO为迭代初始向量；w 为松弛因子；eps为解的精度控制( 此参数可选) ；M 为迭代步数控制( 此参数可选) ；x 为线性方程组的解；为求出所需精度的解实际的迭代步数。SSOR

40、算法用MATLAB实现如下：function xz n=SSOR(A,b,xO,w,eps,M)务采用对称逐次超松弛迭代法求线性方程组Ax=b的解金线性方程组的系数矩阵： A告线性方程组中的常数向量： b超迭代初始向量： xO% 松弛因子： w三解的精度控制： eps告迭代步数控制： M与线性方程组的解： x% 求出所需精度的解实际的迭代步数： nif nargin=4eps= 1.Oe-6;M = 200;elseif nargin4errorreturnelseif nargin =5M = 200; 205精通MATLAB科学计算（ 第2蝌务求A的对角矩阵% 求A的下三角矩阵求A的上三角

41、矩阵endif(w=2)error;return;endD=diag(diag(A);L=-tril(A,-l);U=-triu(A,1);Bl=inv(D-L*w)*(1-w)*D+w*U);B2=inv(D-U*w)*(1-w)*D+w*L);fl=w*inv(D-L*w)*b;f2=w*inv(D-U*w)*b;xl2=Bl*xO+fl;x =B2*xl2+f2;n=l;% 迭代过程while norm(x-xO)=epsxO=x;xl2=Bl*xO+fl;x =B2*xl2+f2;楮第一步的迭代矩阵带第二步的迭代矩阵3第一步的迭代得到的中间值称迭代次数n=n+l;if(n=M)disp

42、 (Warning:迭代次数太多， 可能不收敛！ D;return;endend【 例 9-12 对称逐次超松弛迭代法求解线性方程组实例。用 SSOR迭代法求解以下线性方程组， 其中取初始值为 0Q0 。4为 + 3% 2 = 24第 9 章 线性方程组求解输出的迭代次数为：n =34%可 见 , 经过34步迭代, SSOR迭代法求出方程组的解为M,X2,X3 = 3,4,-5。上面的两个例子采用的是同一个线性方程组， 且采用的松弛系数都是1.25 , 从求解所需的步骤来看, SOR迭代法（14步迭代） 比 SSOR迭代法（ 34步迭代） 要快一些。9 .3 .6两步迭代法两步迭代法的迭代公式

43、如下：(D-L)Xk+Mk Uxk+b(D-U)xk+ 1 =Lxk+i.2+b其中D、4。的定义和前面一样， 即D为矩阵/ 分解的对角矩阵，L为矩阵4 分解的下三角矩阵， 。为矩阵4 分解的上三角矩阵。在 MATLAB中编程实现的两步迭代法函数为： twostepo功 能 ： 用 两 步 迭 代 法 求 线 性 方 程 组 的 解 。调用格式： x,/?= twostep （A, b, xO, eps, varargin）o其 中 ，A为线性方程组的系数矩阵；b 为线性方程组中的常数 向 量 ；xO为迭代初始 向 量 ；eps为解的精度控制（ 此参数可选） ；varargin为迭代步数控制（

44、 此参数可选） ；x 为线性方程组的解；n为求出所需精度的解实际的迭代步数。两步迭代法用MATLAB实现如下：function x,n=twostep（ A,b,xO,eps,varargin）为采用两步迭代法求线性方程组Ax=b的解软线性方程组的系数矩阵： A告线性方程组中的常数向量： b告迭代初始向量： xO志解的精度控制： eps为迭代步数控制：varargin% 线性方程组的解： x% 求出所需精度的解实际的迭代步数： nif nargin=3/ 207精通MATLAB科学计算( 第2忡 - - - - - - - - - - - - - - - - - - - - - - - - -

45、 - - - -eps= 1. Oe-6;M = 200;elseif nargin=epsxO =x;xl2=Bl*xO+fl;x =B2*xl2+f2;n=n+l;if (n=M)disp ( Warning:迭代次数太多， 可能不收敛！， ) ；return;endend【 例 9-13 两步迭代法求解线性方程组实例。 采用两步迭代法求解以下线性方程组,其中取初始值为0,0,0。4xi + 3% 2 = 24 xz n=twostep(A,b,xO)输出的计算结果为：X =3.00004.0000208 第9章 线性方程组求解-5.0000输出的迭代次数为：n =30如可 见 ， 经过3

46、0步迭代， 两步迭代法求出了方程组的解为5 户2,均 = 3 ,4 , 。从上面两个例子可以看出， 与两步迭代法比， SSOR迭代法稍慢一些。9 .3 .7梯度法梯度法相对于前面的几种对系数矩阵的简单分解而构造的迭代法而言， 它不需考虑收敛的问题， 而且求解的速度也相对较快。对于系数矩阵对称正定的方程组， 梯度法有最速下降法、共羯梯度法和预处理的共拆梯度法这三种常见的求法， 下面分别进行讲述。1 . 最速下降法最速下降法的基本思路是将求线性方程组：Ax=b的解转化为求二次泛函的极小值问题。( p(x) =xAx-2bYx通常的做法是： 先任意给定一个初始向量， 然后确定一个搜索方向和搜索步长，

47、 如此循环直到找到极小值。最速下降法采取的搜索方向是当前点的负梯度方向， 每次的搜索步长都使得泛函取极小值。在 MATLAB中编程实现的最速下降法函数为： fastdowno功 能 ： 用最速下降法求线性方程组Nx=b的解。调用格式： x,n= fastdown A, h, xO, eps)o其 中 ， A 为线性方程组的系数矩阵；b 为线性方程组中的常数向量；xO为初始向量；eps为解的精度控制( 此参数可选) ；x 为线性方程组的解；n 为求出所需精度解的实际的搜索步数。最速下降法用MATLAB实现如下：function x,n=fastdown(A,b,xO,eps)会采用最速下降法求线

48、性方程组Ax=b的解 209精通MATLAB科学计算（ 第2蝌常线性方程组的系数矩阵：A% 线性方程组中的常数向量：b% 初始向量：xO故解的精度控制：eps传线性方程组的解：x% 求出所需精度的解实际的搜索步数：nif(nargin = 3)eps = 1 .Oe-6;endr = b-A*xO;d = dot(r,r)/dot(A*rz r);x = xO+d*r;n=l;与迭代过程while(norm(x-xO)eps)xO = x;r = b-A*xO;d = dot(rz r)/dot(A*rz r);x = xO+d*r;n = n + 1;end【 例 9-14】最速下降法求解线

49、性方程组实例。采用最速下降法求解以下线性方程组，且分别取迭代初始值1,1,1,1和2,10,-1,2 , 分析迭代初值对迭代步数的影响。X + 0.5-2 + 0.3333x3 + 0.25-4 = 10.5xi + 0.3333工 2 + 0.25x3 + 0.2x4 = 00.3333x1 + 0.25%2 + 0.2%3 + 0.1667x4 = 10.25%) +0.2-2 +0.1667-3 + 0.1429x4 = 0解：用最速下降法求解， 在 MATLAB命令窗口中输入求解程序：A = 1.00000.50000.33330.2500b= 1 0 1 01;0.50000.333

50、30.25000.20000.33330.25000.20000.16670.2500;0.2000;0.1667;0.1429;x0=l 1 1 1 1;告迭代初值x, n=fastdown(A,b,xO)输出的计算结果为：x =1.0e+003 *0.3288-3.60888.5860210 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 第9章 线性方程组求解-5.5404输出的迭代次数为：n 101383可

51、见 ， 当迭代初始值为1,1,1,1时 ， 经 过 101383步迭代， 最速下降法求出了方程组的解为XI,X2,X3,X4 = 328.8, -3608.8, 8586, -5540.4 o如果改变迭代的初始值， 那么迭代步数会怎样变化呢？在 MATLAB命令窗口中接着输入： x0=2 10 -1 2 , ; 3新的迭代初值 xz n =fastdow n(Az b,x0)输出的计算结果为：X =1.0e+003 *0.3288-3.60888.5860-5.5404输出的迭代次数为：n =8136581365可 见 ， 当迭代初始值为2,10,-1,2时 ， 经 过 81365步迭代， 最

52、速下降法求出了方程组的解为xi,X2,X3,X4 = 328.8, -3608.8, 8586, -5540.4 0从迭代结果可以看出， 选用新的迭代初值， 所需步数大大降低了， 说明对最速下降法而 言 ， 初始值的选取很重要， 选得好的话， 能大大降低求解所需的迭代步数。2 . 共聊梯度法共轨梯度法是从整体来寻找最佳的搜索方向的。它的思路是第一步仍然取负梯度方向作为搜索方向， 对于接下来的各步， 在过当前点由负梯度向量和上一步的搜索向量组成的平面内寻找最佳的搜索方向， 因此它比最速下降法快。在 MATLAB中编程实现的共辗梯度法函数为： conjgrado功 能 ： 共 辗 梯 度 法 求

53、线 性 方 程 组 的 解 。调用格式： X,H= conjgrad (A, b, x0)o其 中 ， A 为线性方程组的系数矩阵；6 为线性方程组中的常数向量；x 0 为初始向量；Y 211精通MATLAB科学计算( 第2版i -X为线性方程组的解；n为求出所需精度的解实际的搜索步数。共轨梯度法用MATLAB实现如下：function x,n=conjgrad(A,b,xO)他采用共轨梯度法求线性方程组Ax=b的解名线性方程组的系数矩阵：A阳线性方程组中的常数向量：b% 初始向量：xO名线性方程组的解：x超求出所需精度的解实际的搜索步数：nif(nargin = 3)eps = 1.Oe-6

54、;endrl = b-A*xO;pl = rl;d = dot(rl,rl)/dot(plzA*pl);x = xO+d*pl;r2 = rl-d*A*pl;f = dot(r2,r2)/dot(rlz rl);p2 = r2+f*pl;n = 1;for(i=l:(rank(A)-1)xO = x;pl = p2;rl = r2;d = dot(rl,rl)/dot(pl,A*pl);x = xO+d*pl;r2 = rl-d*A*pl;f = dot(r2,r2)/dot(rlz rl);p2 = r2+f*pl;n = n + 1;endd = dot(r2,r2)/dot(p2zA*p

55、2);x = x+d*p2;n = n + 1;【 例 9-15】共辗梯度法求解线性方程组实例。 采用共拆梯度法求解以下线性方程组,其中取初始值为0,0,0,0,0。25xi - 300x2 +1050x3 - 1400x4 + 630x5 = 5-300x1 + 4800x2 -1 8900%3 + 26880x4 -1 2600工 5 = 3- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 第9章 线性方程组求解X

56、 =解: 用共轨梯度法求解， 在 MATLAB命令窗口中输入求解程序：A= 25 -3001050-1400630;-300 4800-1890026880-12600;1050 -1890079380-11760056700;-1400 26880-117600179200-88200;630 -12600 b= 5 3 - 1 0 -2 匕 ； x0=0 0 0 0 0 ; x ,n = co n jg rad (A ,bz xO)输出的计算结果为：56700-882004 4 1 0 0 ; ;5.76672.91671.93101.43331.1349输出的迭代次数为：n =6可 见

57、， 当迭代初始值为 0,0, 0,0,0 时 ， 经过6 步迭代， 共聊梯度法求出了方程组的解为 西， 必， 冷， M,X5 = 5.7667,2.9167,1.9310,1.4333,1.1349 。输入下面的表达式以验证结果， 在 MATLAB命令窗口中输入： A*x输出的计算结果为：ans =5.00003.0000-1.0000-0.0000-2.0000共聊梯度法求解的速度非常快， 通常比最速下降法提高了好几个数量级。3 . 预处理的共拆梯度法当 为 病 态 方 程 组 时 ， 共辗梯度法会收敛得很慢。预处理技术是在用共轨梯度法求解之前对系数矩阵做一些变换后再求解的技术。在方程组两边

58、都乘以一个非奇异矩阵F , 令A=FAFTx=PTxh=FbY 213精通MATLAB科学计算( 第2忡 - - - - - - - - - - - - - - - - - - - - - - - - - - - - -对 方 程 用 共 朝 梯 度 法 求 解 ， 则严x 为原方程的解。理想情况下是取尸= 4 、 但对病态矩阵， 求逆太耗费时间， 实用的取法步骤为：步骤 1 , 取产=diag(l/au , 1/做2l/nn)；步骤2 , 取系数矩阵”的三条对角线构成的矩阵， 再对其进行Cholesky分 解 ， 将尸取为分解后的上三角矩阵。另外尸还有一些其他的取法， 具体可参考数值计算的书

59、籍。MATLAB中有一个对应的处理函数pcg(), 它的用法如下：X, flag = peg (A, bz tol, maxit,Ml,M2, xO, pl,p2,.)；x,flag,relres = peg(A,b,tol,maxit,Ml,M2,xO,pl,p2,)；x, flag, relres, iter = peg (A,b, tol,maxit,Ml,M2, xO,plzp2,.)；x, flag, relres, iterz resvec = peg (A, b, tol, maxit, Ml, M2Z xOz pl, p2,.)； 其 中 ， tol为迭代的精度， maxit为

60、迭代的最大步数， M l和M 2就是预处理矩阵； flag的取值为： 0 表示在指定迭代次数之内按要求精度收敛； 1表示在指定迭代次数内不收敛； 2 表示 为坏条件的预处理因子； 3 表示两次连续迭代完全相同； 4 表示标量参数太小或太大； relres 为相对误差： norm(6-*x)/norm(6)； iter为计算x 的迭代次数； resvec为每次迭代的残差： normS-N*xO)。【 例 9-16】预处理的共轨梯度法求解线性方程组实例。采用预处理的共轨梯度法求解以下线性方程组， 其中取初始值为0,0,0,0,0。25xi 300x2 + 1050x3 - 1400x4 + 630

61、xs = 5 300xi + 4800x2 18900x3 + 26880x4 - 12600x5 = 3- 1050xi -18900x2 +79380x3 -117600x4 +56700x5 =-1-1400X1 + 26880x2 -117600x3 +179200x4 - 88200x5 = 0630xi - 12600%2 + 56700x3 - 88200x4 +44100x5 = 2解 ：用预处理的共轨梯度法求解， 在 MATLAB命令窗口中输入求解程序：A= 25-3001050-1400630;-3004800-1890026880-12600;1050-189007938

62、0-11760056700;-140026880-117600179200-88200;630-1260056700-8820044100;b= 5 3 -10 -2x0=0 0 0 00M=pascal(5); 告预处理矩阵x,flag,rezit=pcg(Azbz1. e- 8, 1000, M,M,x0)214 第9章 线性方程组求解9 ；输出的预处理矩阵M =111111 12 33 64 105 151 14 510 1520 3535 70输出的计算结果为：5.76672.91671.93101.43331.1349输出的flag取值为0 , 表示在指定迭代次数之内按要求精度收敛f

63、lag =0输出的相对误差：re =3 .6125e-0105 . 4 9 9 0 .0 0 9输出的计算X 的迭代次数：it = 109可 见 ， 当迭代初始值为 0,0,0。0 时 ， 经 过 皿 步 迭 代 ， 预处理的共森梯度法求出了方程组的解为 的, X2,X3,X4,X 5 = 5.7667, 2.9167,1.9310,1.4333,1.1349 。9 .3 .8其他迭代法MATLAB还有一些迭代法， 其用法和参数与预处理共轨梯度法基本上一样， 如表9-2所示。表 9 - 2 其他的线性方程组迭代法求解函数函数及语法说 明x = symmlq(4,b)线性方程组的LQ解法x =

64、bicg(46)线性方程组的双共规梯度解法x = bicgstab(J,d)线性方程组的稳定双共较梯度解法x = lsqr(4b)线性方程组的共聊梯度的LSQR解法Y 215精通MATLAB科 学 计 算 ( 第2蚪【 例 9-17】最小残差法求解线性方程组实例。 采用最小残差法求解以下线性方程组，x = gmres(J)线性方程组的广义最小残差解法x = minres(/4,6)线性方程组的最小残差解法x = qmr(4b)线性方程组的准最小残差解法3% 2 + 8工3 = 1 即 + 4X2 + 7工3 = 1-2xi + 6x2 + 9x3 = 1解 ：用最小残差法求解， 在 MATLA

65、B命令窗口中输入求解程序： A = 0 3 8;1 4 7; -2 6 9; b= 1 1 1 1; xz f la g ,r e , it= m in res(A ,b )输出的计算结果为：X =0.07540.08070.0826输出的flag取值为1 , 表示在指定迭代次数之内不收敛：fla g =1输出的相对误差：re =0.0519输出的计算X 的迭代次数i t =3由计算结果可知， 经过迭代3 步 ， flag= l,表明默认情况下不收敛， 因此需要选择其他的解法来求解。216 第 9 章 线性方程组求解9.4特殊解法9 .4 .1三对角矩阵的追赶法对于系数矩阵为对角占优的三对角矩

66、阵， 有一种特殊的解法叫追赶法。追赶法首先对系数矩阵N做如下三角分解：A=LU其中为下三角的两对角矩阵， 。为上三角的两对角矩阵， 且主对角线元素全为lo求解时先求解0 = 6 的解y , 再 求 解 但 的 解 x。在 MATLAB中编程实现的追赶法函数为： followupo功 能 ： 用 追 赶 法 求 线 性 方 程 组 的 解 。调用格式： x= followup (A,b)o其 中 ， A 为线性方程组的系数矩阵；b 为线性方程组中的常数向量；x 为线性方程组的解。追赶法用MATLAB实现如下：function x=followup(Azb)巷采用追赶法求线性方程组Ax=b的解常线

67、性方程组的系数矩阵：A带线性方程组中的常数向量：b为线性方程组的解：xn = rank(A);for(i=l:n)if(A(i,i)=0)disp ( Error:对角有元素为0 !， ) ；return;endend;d = ones(n,1);a = ones(n-1,1);c = ones (n-1);for(i=l:n-1)a(i,l)=A(i+lzi);c(i, 1)=A(i,i+1);d(i,l)=A(i,i);endd(n, 1) = A(n,n);g求解Ly=b的解y ,解保存在b中 ，1 217精通MATLAB科学计算( 第2忡 - - - - - - - - - - - -

68、 - - - - - - - - - - - - -f o r( i= 2 :n)d ( izl) = d ( i ,l) - ( a ( i - 1 ,1 )/ d ( i - lz1 ) ) * c ( i - lr1 );b ( izl) = b ( i ,l) - ( a ( i - lzl ) / d ( i - l , l ) ) * b ( i - l , l ) ;end务求解Ux=y的解x ,x (n ,1) = b ( n ,1 )/d ( n ,1 );f o r ( i= ( n - 1 ) : - 1 : 1)x ( i , 1) = ( b ( i ,1 )- c

69、( iz1 )* x (i + 1 ,1 ) / d ( i , l ) ;end【 例 9-18】追赶法求解线性方程组实例。采用追赶法求解以下线性方程组，2.5即 + 工2 = 1X + 1.5X2 + 工 3 = 1X2 + 0.5%3 + 必 =1冷 + 0.5%4 + 工 5 = 1% 4 +1.5工 5 + 冗6 = 1落X 5 + 2.56=1解：用追赶法求解， 在 MATLAB命令窗口中输入求解程序： b = o n es(6,1) x= follow up(A ,b)A = 2 . 50001.00000000；1.00001.50001.0000000;01.00000.50

70、001.000000;001.00000.50001.00000;0001.00001.50001.0000;00001.00002.5000输出的计算结果为：X =0.4615-0.15380.76920.7692-0.15380.4615可 见 ， 追赶法求出了方程组的解为：xi,X2,x3,X4,X5,X6 = 0.4615, -0.1538,0.7692,0.7692, -0.1538,0.4615 o最后可以输入下面的表达式以验证结果的正确性， 在 MATLAB命令窗口中输入： A*x输出的计算结果为：ans =218 第 9 章 线性方程组求解1.00001.00001.00001

71、.00001.00001.0000可 见 ， 追赶法求解的结果是正确的。9 .4 .2快速求解法MATLAB中有一个通用的求解线性方程组的函数， 它的用法如下：x = linsolve(/4, b, options)其意义为快速求解方程， 其中”的结构由options决 定 ， 内容如表93 所示。表 9-3 options取值表options说 明LT上三角矩阵UT下三角阵UHESS上三角Hessenberg矩阵SYM实对称阵POSDEF正定矩阵RECT一般矩阵TRANSA指出求解的方程是Ax=b还是Ax=b , ， 指A的共枕转置【 例 9-19】快速求解法求解线性方程组实例。采用快速解法

72、求解以下线性方程组，4xi - 42 + 2x3 = 1 optt.TRANSA=false; x=linsolve(A,broptt)输出的计算结果为：X =0.0426-0.1915词 x=lsqnonneg (A, b) % 用最小二乘法求解输出的计算结果为：000.0826从上述结果可以看出， 采用不同的解法， 得到的结果不完全相同。9 .5 .2有无穷组解的线性方程组的解法对于此类方程的解法， 可分为以下3 个步骤：步 骤 1-士 求 的 一 个 特 解 ；步骤2七 求 / 1 4 0 的通解；步骤3七 将特解与通解组合成最终的解。【 例 9-21 有无穷组解的线性方程组求解实例。求

73、解以下线性方程组，4X 1 + 8X2 - 6X3 + 2%4 = 1X + 3% 2 + 9X3 + 5X4 = 25% 1 +11X2+ 3冷 + 7X4 = 33xj + 52 15工 3 3x4 = - 1解：在 MATLAB命令窗口中输入求解程序： A=4 8 -6 2;1 3 9 5;5 11 3 7;3 5 -15 -3; b=l 2 3 -1 1; D=A b S=rref (D) % 采用增广矩阵法求解输出的计算结果为：S =1.00000 -22.5000 -8.5000 -3.25000 1 .00000 00 010.5000 4.50000 00 01.750000上

74、面例子是采用增广矩阵法求解的， 从输出的S 结果可以看出， 此方程的一个特解为x0= -3.25 1.75 0 0T , 其基础解系有两个基向量：-22.5 10.5 1 0T , -8.5 4.5 0 1T222 则最终的解为x=x（） + / ： i f i + A ： 2$。第9章 线性方程组求解 223精通MATLAB科 学 计 算 ( 第2蚪9.6小结本章详细讲述了在MATLAB中求解线性方程组的常用解法， 包括求逆法、分解法、迭代法、特殊解法， 以及非齐次线性方程组的解法。每种求解方法也各有其优缺点， 适用于不同的方程组。选择线性方程组解法的主要原则可总结如下：(1 ) 首先判断系数矩阵是否为方阵，是方阵就转下面的步骤，否则此方程要么无解,要么有无穷解， 此时用非齐次线性方程组的解法求解。(2 )观察系数矩阵的特征， 是否为对称矩阵或正定矩阵， 它们适合用矩阵分解法求解。(3 ) 稀疏矩阵和病态矩阵可考虑用迭代法求解。(4 )特殊解法适合特殊方程。224