MATLAB计算方法3解线性方程组计算解法

1、解线性方程组解线性方程组 bAx x解向量解向量列列迭代法：格式，无穷序迭代法：格式，无穷序误差，有限步，精确解误差，有限步，精确解直接法：理论，无舍入直接法：理论，无舍入第三章线性方程组数值解法第三章线性方程组数值解法3.3.1 1 直接法直接法一、一、 GaussGauss消去法消去法设设 有有1,22111,2211)1.3(1,222221211, 11212111nnnnnnnnininiinnnnnnaxaxaxaaxaxaxaaxaxaxaaxaxaxa行。行。行，加到第行，加到第第第把把化为零；化为零；将将用用 1),2(111111iaaniaaii 以后各步类似。以后各步类

2、似。或或题：题：问问?001111 aa消消 元：元： 用用MatlabMatlab实现顺序实现顺序GaussGauss消去法消去法在在MatlabMatlab程序编辑器中输入：程序编辑器中输入：function x=nagauss(a,b,flag) %function x=nagauss(a,b,flag) %解线形方程组解线形方程组ax=bax=b，a a为系数矩阵，为系数矩阵，b b为右端列向量，为右端列向量，flagflag若为若为0 0，则显示中间，则显示中间过程，否则不显示，默认为过程，否则不显示，默认为0 0，x x为解向量为解向量if nargin3,flag=0;endif

3、 nargin3,flag=0;endn=length(b); a=a,b;n=length(b); a=a,b;% % 消元消元for k=1:(n-1)for k=1:(n-1) a(k+1):n,(k+1):(n+1)=a(k+1):n,(k+1):(n+1)-a(k+1):n,(k+1):(n+1)=a(k+1):n,(k+1):(n+1)-a(k+1):n,k)/a(k,k)a(k+1):n,k)/a(k,k)* *a(k,(k+1):(n+1);a(k,(k+1):(n+1); a(k+1):n,k)=zeros(n-k,1); a(k+1):n,k)=zeros(n-k,1);

4、if flag=0,a,end if flag=0,a,endend end % % 回代回代x=zeros(n,1);x=zeros(n,1);x(n)=a(n,n+1)/a(n,n);x(n)=a(n,n+1)/a(n,n);for k=n-1:-1:1for k=n-1:-1:1 x(k,:)=(a(k,n+1)- x(k,:)=(a(k,n+1)-a(k,(k+1):n)a(k,(k+1):n)* *x(k+1):n)/a(k,k);x(k+1):n)/a(k,k);endend程序运行程序运行结果：结果：二二 列主元素列主元素GaussGauss消去法消去法-计算结果可靠计算结果可靠

5、行：行：，对调，对调使使找行号找行号111111max)1(1raarinir ）；（个元素成为个元素成为行第行第第第行行加到第加到第行行第第为为消消消元：用消元：用1, 2, 1, 3, 2,1:01111111111 njniaaaaajiiaaaaijijijii到此原方程组化为到此原方程组化为1,22 1,22 1,222221, 11212111 nnnnnnninininnnnnnaxaxaaxaxaaxaxaaxaxaxa .2max)2(222222行行对调对调，使，使找找raarinir ），；，（行，则行，则第第行行第第消为消为把把消元：用消元：用1,3 , 2,4 , 3

6、 2:), 4 , 3( 02222222222 njniaaaaaiaaniaaijijijii( (上三角方程组上三角方程组) ) （3.2)3.2)1, 1,222221, 11212111 nnnnnnnnnnnaxaaxaxaaxaxaxa直到直到(n-1) (n-1) 原方程组化为原方程组化为以上为消元过程。以上为消元过程。)1,.,2, 1(111,1, nnkxaaaxaaxnkjjkjnkkkknnnnn （3.3)3.3)( (n) n) 回代求解公式回代求解公式 Gauss Gauss 列主元消去法：列主元消去法： 优点优点 - - 计算结果更可靠；计算结果更可靠； 缺点

7、缺点 - - 挑主元花机时更多挑主元花机时更多, , 次序次序有变动，程序复杂。有变动，程序复杂。说明：说明： (1) (1)也可采用无回代的列主元消去法也可采用无回代的列主元消去法( (叫叫Gauss-Gauss- -Jordan -Jordan消去法消去法) )，但比有回代的列主元消，但比有回代的列主元消 去法的乘除运算次数多。去法的乘除运算次数多。 (2) (2)有回代的列主元消去法所进行的乘除运算有回代的列主元消去法所进行的乘除运算 次数为次数为 ，量很小。，量很小。nnn313123 nxx,1用用MatlabMatlab实现选列主元实现选列主元GaussGauss消去法解线性方程组

8、消去法解线性方程组在在MatlabMatlab程序编辑器中输入：程序编辑器中输入：function x=nagauss2(a,b,flag) %afunction x=nagauss2(a,b,flag) %a为系数矩阵；为系数矩阵；b b为为右端列向量；右端列向量；flagflag若为若为0 0，则显示中间过程，否则不显示，则显示中间过程，否则不显示if nargin3,flag=0;endif nargink,t=a(k,:); a(k,:)=a(p,:); a(p,:)=t; endif pk,t=a(k,:); a(k,:)=a(p,:); a(p,:)=t; end% % 消元消元

9、a(k+1):n,(k+1):(n+1)=a(k+1):n,(k+1):(n+1)-a(k+1):n,(k+1):(n+1)=a(k+1):n,(k+1):(n+1)-a(k+1):n,k)/a(k,k)a(k+1):n,k)/a(k,k)* *a(k,(k+1):(n+1);a(k,(k+1):(n+1); a(k+1):n,k)=zeros(n-k,1); a(k+1):n,k)=zeros(n-k,1); if flag=0,a,endif flag=0,a,endend end % %回代回代x=zeros(n,1);x=zeros(n,1);x(n)=a(n,n+1)/a(n,n);

10、x(n)=a(n,n+1)/a(n,n);for k=n-1:-1:1for k=n-1:-1:1 x(k,:)=(a(k,n+1)- x(k,:)=(a(k,n+1)-a(k,(k+1):n)a(k,(k+1):n)* *x(k+1):n)/a(k,k);x(k+1):n)/a(k,k);endend程序运行程序运行结果结果: : uALLLLgubALLLLgubAKKkk 121121)|()|()1|()|()|(uLLLLAKK1121)( 记记11)( LLLK，则，则LUA 上上三三角角）（下下三三角角 （三角因子分解）（三角因子分解） GaussGauss消元，初等行变换，化原

11、方程组为上三消元，初等行变换，化原方程组为上三角型。角型。三矩阵三角分解法三矩阵三角分解法定义定义3.13.1 LUA 叫叫A的三角（因子）分解，其中的三角（因子）分解，其中 是是L是上三角。是上三角。下三角下三角, , L为单位下三角阵（对角元全为为单位下三角阵（对角元全为1 1），），U为上三角阵，则称为上三角阵，则称LUA 为为DoolittleDoolittle分解分解; ;L若若 是下三角，是下三角， 是单位上三角，则称是单位上三角，则称定义定义3.23.2 若若 为为CroutCrout分解。分解。UULUA 定义定义3.2 3.2 若若L 为什么要讨论三角分解？若在消元法进行前能

12、实为什么要讨论三角分解？若在消元法进行前能实 现三角分解现三角分解LUA ， 则则 bxLUbAx)(上三角方程组）上三角方程组）下三角方程组）下三角方程组）( ( yUxbLy 从而容易回代求解。从而容易回代求解。 nnnnnnaaaaaaaaa2122221112111 1直接三角分解法直接三角分解法( (DoolittleDoolittle分解为例）分解为例） 1112121nnlll nnnnuuuuuu22211211由矩阵乘法由矩阵乘法221212222221211212222121111111111111/ )( n),3,4, 2i22) 1 n),2, j22)1 )2(),

13、 3 , 2(1L12) n),1,2,(j u 1. 1)1 )1(uulalaululiuLLulauauuljuLuniualauluiLauajuLuiiiiiijjjjjjiiiijjjj （列列的第的第行行的第的第列：用列：用的第的第求求（列列的第的第行行的第的第行：用行：用的第的第求求列列的第的第行行的第的第列：用列：用的第的第求求列列的第的第的第一行的第一行行：用行：用的第的第求求 .(k) 1111,2,1,2,11)0 ,0,()0 ,0,(n),1,kk,(j1kmmjkmkjkjkjkjkmmjkmkjjjjjkkkkulauauulauuullljukLku，列列的第

14、的第行行的第的第行：用行：用的第的第求求 kkkmmkimikikikkkikkmmkimikkkkkkkikiuulalaululauuulllkuiLkL 1111,2,1,1)00,()00(n),1,k(i2，即，即列列的第的第行行的第的第列：用列：用的第的第求求 ), 3 , 2(), 1( ), 1,( ), 3 , 2( ), 2 , 1( 1111111111nknkiauulalnkkjaulauniualnjauLUikkkkmmkimikikkjkmmjkmkjkjiijj分解式：分解式：例例3 31 112/ )126(23/ )06(, 111 1661230 3 2

15、2 1 251865364 5 42 1 2231323 xxxxxx所以所以用用MatlabMatlab实现实现LULU分解分解在在MatlabMatlab程序编辑器中输入：程序编辑器中输入：function L,U=nalu(a) % afunction L,U=nalu(a) % a为可逆方阵；为可逆方阵；L L返回单位下返回单位下三角矩阵；三角矩阵；U U返回上三角矩阵返回上三角矩阵n=length(a);n=length(a);U=zeros(n,n);L=eye(n,n);U=zeros(n,n);L=eye(n,n);U(1,:)=a(1,:);L(2:n,1)=a(2:n,1)

16、/U(1,1);U(1,:)=a(1,:);L(2:n,1)=a(2:n,1)/U(1,1);for k=2:nfor k=2:n U(k,k:n)=a(k,k:n)-L(k,1:k-1) U(k,k:n)=a(k,k:n)-L(k,1:k-1)* *U(1:k-1,k:n);U(1:k-1,k:n); L(k+1:n,k)=(a(k+1:n,k)-L(k+1:n,1:k-1) L(k+1:n,k)=(a(k+1:n,k)-L(k+1:n,1:k-1)* *U(1:k-U(1:k-1,k)/U(k,k);1,k)/U(k,k);endend程序运行程序运行结果结果: :2 2平方根法平方根法定

17、理定理3.13.1 设设A A对称正定，则有非奇异下三角阵对称正定，则有非奇异下三角阵L L，使使TLLA 分解方法：设分解方法：设jiijnnnnnnnnnnnnnnaallllllllllllaaaaaaaaa 其中其中 222121112122211121222211121122121222222212121222222222111111111111111111 n),3,4,(j j2L )2 22L 1) )2( n),2,(j j1L )2 )( )1 )1(2222lllalaallllLlalallLlalaallLalaljjjjjjjTTjjjjjT 列列第第行行第第列列第

18、第行行第第列列第第行行第第取正取正由矩阵乘法由矩阵乘法）（，可得：可得：nkjllalllalkkmjmkmjkkkjmkmkmkkkk, 11)(11112 及其改进分解。及其改进分解。缺点：开方运算。缺点：开方运算。主元。主元。平方根法稳定，不必先平方根法稳定，不必先受到控制受到控制所以说舍入误差的放大所以说舍入误差的放大的最大元的最大元的值不会超过的值不会超过程中程中这说明在分解过这说明在分解过故故列：列：第第行行、第、第、计算量小、计算量小优点：优点：TkmkknkkmkknkkkkmkmkmkklalaallakkLDL,A .maxmax2;1211212 大。大。可以控制舍入误差

19、的增可以控制舍入误差的增此时此时为交换前的为交换前的于是于是变。变。且元素的足码也相应改且元素的足码也相应改行行行与第行与第并交换矩阵的第并交换矩阵的第为为取取的大小的大小比较比较于是于是引进量引进量作除数作除数用小用小为了避免为了避免分解时分解时对其增广矩阵作对其增广矩阵作设方程组设方程组, 1|./),(,|)|( |max,|,), 1,(,11 ikkiiktkkkkkktinikikkkkqqkiqikikklsslsssuktussssunkkiulasuLUbAx 6121151318153312A四、解三对角方程组四、解三对角方程组追赶法追赶法）（，且按行对角占优：且按行对角占

20、优：（三对角方程组）（三对角方程组）给定给定1, 2)0( 1112112111122211 niabcacabcbddddxxxxbacbacbacbnniiiiinnnnnnnnn 1 , 2, 1/ )(/1111222111nnixcyxyxyycycyciiiiinnnnnnnn 回代过程为回代过程为 1111 1211122111122211nnnnnnnnnnbacbacbacbLU 设有分解设有分解分解：分解：利用利用列得：列得：第第行行第第，）（由矩阵乘法由矩阵乘法1, 1)(1)2(112221222222222212221111111111 iiiiicbacbacbcb

21、 00 00i i 11111 iiiii ii1 ,1 .,1 .1iiiiiiba , ,ciii 故有故有 ,111111 iiiiiiiiicbacb), 3 , 2(ni (3.1) dLUxdAx （二对角方程组）（二对角方程组）二对角方程组）二对角方程组） ( yUxdLy)1, 3 , 2( ni解解n)2,3,(i , 1111y iiiiiyddydLy设设解解得得 yUx 1)1,-n(i ,1 xyxyxiiiinn (3.1) (3.1) (3.3) (3.3) 叫追赶法，工作量小，叫追赶法，工作量小， 非常有效。非常有效。Ex:Ex:第三章第三章 4 4(3.2)(

22、3.2)（3.33.3）用用MatlabMatlab实现追赶法求解线性方程组实现追赶法求解线性方程组在在MatlabMatlab程序编辑器中输入：程序编辑器中输入：function x=pursue(a,b,flag) % afunction x=pursue(a,b,flag) % a为系数矩阵；为系数矩阵；b b为右为右端列向量；端列向量；flagflag为为0 0则显示过程则显示过程if nargin3,flag=0;end % narginif nargin3,flag=0;end % nargin为计算输入变量个为计算输入变量个数数n=length(b); a=a,b;n=lengt

23、h(b); a=a,b;for k=1:(n-1)for k=1:(n-1) a(k+1,k+1)=a(k+1,k+1)- a(k+1,k+1)=a(k+1,k+1)-a(k+1,k)/a(k,k)a(k+1,k)/a(k,k)* *a(k,k+1);a(k,k+1); a(k+1,n+1)=a(k+1,n+1)- a(k+1,n+1)=a(k+1,n+1)-a(k+1,k)/a(k,k)a(k+1,k)/a(k,k)* *a(k,n+1);a(k,n+1); a(k+1,k)=zeros(1,1); a(k+1,k)=zeros(1,1); if flag=0,a,end if flag=0

24、,a,endendend% % 回代回代x(n)=a(n,n+1)/a(n,n);x(n)=a(n,n+1)/a(n,n);for k=n-1:-1:1for k=n-1:-1:1 x(k)=(a(k,n+1)-a(k,k+1) x(k)=(a(k,n+1)-a(k,k+1)* *x(k+1)/a(k,k)x(k+1)/a(k,k)endend程序运程序运行结果行结果: : 用用MatlabMatlab实现改进平方根法求解线性方程组实现改进平方根法求解线性方程组在在MatlabMatlab程序编辑器中输入：程序编辑器中输入：function x=squareroot(a,b,flag)func

25、tion x=squareroot(a,b,flag)if nargin3,flag=0;endif nargin3,flag=0;endn=length(b);a=a,b;n=length(b);a=a,b;for k=1:(n-1)for k=1:(n-1) a(k+1,1)=a(k+1,1)/a(1,1); a(k+1,1)=a(k+1,1)/a(1,1);endend a afor k=1:(n-2)for k=1:(n-2) for j=k+1:n+1 for j=k+1:n+1 a(k+1,j)=a(k+1,j)-a(1:k,j) a(k+1,j)=a(k+1,j)-a(1:k,j

26、)* *a(k+1,1:k);a(k+1,1:k); end end for l=k+2:n for l=k+2:n a(l,k+1)=(a(l,k+1)- a(l,k+1)=(a(l,k+1)-a(1:k,k+1)a(1:k,k+1)* *a(l,1:k)/a(k+1,k+1);a(l,1:k)/a(k+1,k+1); endend a aendendfor i=1:(n-1)for i=1:(n-1) a(n,n)=a(n,n)-a(i,n) a(n,n)=a(n,n)-a(i,n)* *a(n,i);a(n,i); a(n,n+1)=a(n,n+1)-a(i,n+1) a(n,n+1)=

27、a(n,n+1)-a(i,n+1)* *a(n,i);a(n,i);endenda ax=zeros(n,1);x=zeros(n,1);x(n)=a(n,n+1)/a(n,n);x(n)=a(n,n+1)/a(n,n);for k=n-1:-1:1for k=n-1:-1:1 x(k)=(a(k,n+1)- x(k)=(a(k,n+1)-a(k,(k+1):n)a(k,(k+1):n)* *x(k+1):n)/a(k,k);x(k+1):n)/a(k,k);endend程序程序运行运行结果：结果：。向量范数向量范数的的为上为上则称则称（三角不等式性）。（三角不等式性）。有有）对任意）对任意（

28、奇次性）；（奇次性）；和任意和任意）对任意实常数）对任意实常数（非负性）（非负性）当且仅当当且仅当）对任意）对任意（个条件：个条件：如果它满足以下如果它满足以下上的实函数，上的实函数，是定义在是定义在定义：设定义：设nnnnnRyxyxRyxxccxRxcxxxRxRxxf|,3| ,2;00| ; 0| ,13|)( . |max|x| -|x| -2; |x| -1112211ininiiniixxx 范数：范数：向量的向量的范数：范数：向量的向量的范数：范数：向量的向量的数：数：最常用的是如下三种范最常用的是如下三种范.|max|;|max|;|max|A| .|max|max| ,1|

29、21|211|11|1|21 AXAAXAAXAXAAAAXRRAxxxxxnnn阵范数：阵范数：最常用的是如下三种矩最常用的是如下三种矩即即，。记作。记作矩阵范数矩阵范数的的为为称称上的任一向量范数，上的任一向量范数，是是矩阵范数：设矩阵范数：设.0| |max)( )(|max| ; |max|A|211111 BIBAAAaAanjijniniijnj 即即其中为矩阵的谱半径，其中为矩阵的谱半径，可以证明：可以证明：.| ,|A|,21312 . 321AAA及求：设例。作为矩阵特征值的上界作为矩阵特征值的上界即矩阵的任一范数均可即矩阵的任一范数均可则则为任一矩阵范数，为任一矩阵范数，：

30、设：设定理定理有有）对任意）对任意（有有对任意向量对任意向量有有对任意对任意；和任意和任意）对任意实常数）对任意实常数（当且仅当当且仅当）对任意）对任意性质：（性质：（.|)(|2 . 3.|,5|,)4(|;|,)3(| ,2; 00| ; 0| ,1AARABAABRBAxAAxRARxBABARBAAccARAcAAARAnnnnnnnnnnnnn 在在MatlabMatlab中可以用命令中可以用命令 norm(A,p)norm(A,p)来求范数。来求范数。p p为为(2,1,inf)(2,1,inf)分别对应：分别对应：2-2-范数，范数，1-1-范数和范数和- -范数。范数。举例如右

31、图：举例如右图：bAx 一一. .简单迭代法简单迭代法 1. 1.迭代法建立迭代法建立. . 考虑考虑gBxxbAx ( (矩阵矩阵B B不唯一不唯一) )对应写出对应写出 )0()()1()4 . 3( ), 2 , 1 , 0( xkgBxxkk取定初始向量取定初始向量 3.3 3.3 解线性方程组的迭代法解线性方程组的迭代法产生向量序列产生向量序列,)1()()2()1( kkxxxx若收敛若收敛, ,记记xxkk )1(lim则于则于(3.4)(3.4)两端取极限有：两端取极限有：, gxBx 上式说明：上式说明： 是解向量是解向量 , ,从而当从而当k k充分大时充分大时 )1(kx

32、注意注意: : 迭代阵迭代阵B B不唯一不唯一, ,影响收敛性。影响收敛性。 解向量解向量(3.4)(3.4)叫简单迭代法叫简单迭代法, ,B B叫迭代矩阵。叫迭代矩阵。xxx 5 . 05 . 05 . 05 . 0 3 . 32121xxxx用迭代法解方程组用迭代法解方程组例例. 1 3.4()0()()1( ）（都收敛都收敛对任意初始向量对任意初始向量简单迭代法简单迭代法充要条件）充要条件）定理定理BxgBxxkk , 2 , 1 |1|*|)4(, 2 , 1 |1|*|)3(, 2 , 1 , 0 |*|*|*,)2(.*)1(, 1|)(3 . 3)0()1()()1()()()(

33、)1()0( kxxBBxxkxxBBxxkxxBxxxRxxBkkkkkkkn，且有，且有迭代格式收敛于迭代格式收敛于对任意初始向量对任意初始向量方程组有唯一解方程组有唯一解则则给定方程组，如果给定方程组，如果充分条件判别法充分条件判别法定理定理.1:)0(都不收敛都不收敛时不能说对任意时不能说对任意）（注意注意xB 二二. .JacobiJacobi迭代法迭代法( (k=0,1,2,k=0,1,2,) )(3.7)(3.7) )()(2)(1)(2)(2)(12222)(1)(212)(1111100021 1.21 1 1knkknnnnknnkkknnkkxxaxabaxxaxxaba

34、xxaxaxbaxnn矩阵形式为：矩阵形式为： nnnkkknnnnnnnnknkkabababxxxaaaaaaaaaaaaxxx2221112212122222211111112)1()1(2)1(1 0001)()7 . 3(:)0()()1( JkJkBxJacobigxBx 收敛收敛对任意对任意迭代法迭代法收敛性收敛性简记为简记为 nknnnknknknknnkkkknnkkkgggJacobiseidelGaussxbxbxbxxbxbxbxxbxbxbx)()(22)(11)1(2)(2)(222)(121)1(21)(1)(212)(111)1(1:.迭代格式为迭代格式为设方程

35、组设方程组迭代法迭代法三三(3.8)(3.8) 称如下迭代法称如下迭代法 (3.9)(3.9) nknnnknknknknnkkkknnkkkgggxbxbxbxxbxbxbxxbxbxbx)()1(22)1(11)1(2)(2)(222)1(121)1(21)(1)(212)(111)1(1 为与为与(3.8)(3.8)对应的对应的 迭代法，其迭代矩阵迭代法，其迭代矩阵 可用可用 “代入法代入法”求得。求得。SeidelBsSeidelGauss 迭代法的收敛性迭代法的收敛性. 1)()0( sBxSeidelGauss 收敛的充要条件是收敛的充要条件是迭代法对任意迭代法对任意.)2(.)1

36、()(5 . 3迭代法收敛迭代法收敛为对称正定阵，则为对称正定阵，则若若迭代法均收敛迭代法均收敛和和迭代法迭代法为严格对角占优，则为严格对角占优，则若系数矩阵若系数矩阵：对于给定的线性方程组对于给定的线性方程组充分条件判别法充分条件判别法定理定理seidelGaussAseidelGaussJacobiA 例例3.4 3.4 分别用分别用JacobiJacobi迭代法和迭代法和SeidelSeidel迭代法求解方程组迭代法求解方程组 1015352111021210321xxxTx)0 , 0 , 0()0( 10max3)()1(31| xxkikii取初始向量取初始向量, ,要求要求 时迭代终止。时迭代终止。解解：因为系数矩阵严格对角占优，故因为系数矩阵严格对角占优，故JacobiJacobi迭代迭代法和法和SeidelSeidel方法对任意