三章解線性方程組的直接方法

1、用二早用牛雙I土刀 住羽陽旦依力 壯在這章中我們要學(xué)習(xí)線性方程組的直接法，特別是適合用數(shù)學(xué)軟件在計算機上求解 的方法.方程組的逆矩陣解法及其MATLA跑序3.1.3線性方程組有解的判定條件及其MATLA醒序判定線性方程組 A m0 X =b是否有解的MATLA典序function RA,RB,n=jiepb(A,b)B=A b;n=length(b); RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica0,disp(請注意：因為RA=RB,所以此方程組無解.) returnendif RA=RBif RA=ndisp(請注意：因為RA=RB=n ,所以此方程組

2、有唯一解 .) elsedisp(請注意：因為RA=RB A=2 3 -1 5;3 1 2 -7;4 1 -3 6;1 -2 4 -7;b= 0; 0; 0; 0; RA,RB,n=jiepb(A,b)運行后輸出結(jié)果為請注意：因為 RA=RB=n ,所以此方程組有唯一解RA = 4,RB =4,n =4在MATLAB：作窗口輸入X=Ab,運行后輸出結(jié)果為 X =(0 0 0 0).(2)在MATLABT作窗口輸入程序 A=3 4 -5 7;2 -3 3 -2;4 11 -13 16;7 -2 1 3;b= 0; 0; 0; 0;RA,RB,n=jiepb(A,b)運行后輸出結(jié)果24 / 13請

3、注意：因為RA=RB A=4 2 -1;3 -1 2;113 0; b=2;10;8; RA,RB,n=jiepb(A,B)運行后輸出結(jié)果請注意：因為RA=RB,所以此方程組無解.RA =2,RB =3,n =3(4)在MATLABT作窗口輸入程序 A=2 1 -1 1;4 2 -2 1;2 1 -1 -1;b=1; 2; 1; RA,RB,n=jiepb(A,b)運行后輸出結(jié)果請注意：因為RA=RB0,disp(請注意：因為RA=RB,所以此方程組無解.,) return end if RA=RBif RA=ndisp(請注意：因為RA=RB=n ,所以此方程組有唯一解 .,)X=zeros

4、(n,1); X(n)=b(n)/A(n,n);for k=n-1:-1:1X(k)=(b(k)-sum(A(k,k+1:n)*X(k+1:n)/A(k,k); end elsedisp(,請注意：因為RA=RBA=5 -1 2 3;0 -2 7 -4;0 0 6 5;0 0 0 3;b=20; -7; 4;6;RA,RB,n,X=shangsan(A,b)運行后輸出結(jié)果請注意：因為RA=RB=n,所以此方程組有唯一解.RA = RB =4,4, 4,25 / 13X =2.4 -4.0 -1.0 2.03.3 高斯（Gaus9消元法和列主元消元法及其 MATLA能序高斯消元法及其MATLAB

5、?序用高斯消元法解線性方程組AX =b的MATLA翼序function RA,RB,n,X=gaus(A,b)B=A b; n=length(b); RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica0,disp(請注意：因為RA=RB,所以此方程組無解.) returnendif RA=RBif RA=ndisp(請注意：因為RA=RB=n ,所以此方程組有唯一解.)X=zeros(n,1); C=zeros(1,n+1);for p= 1:n-1for k=p+1:nm= B(k,p)/ B(p,p);B(k,p:n+1)= B(k,p:n+1)-m* B

6、(p,p:n+1);end endb=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); endelsedisp(請注意：因為RA=RBn ,所以此方程組有無窮多解.)endend例3.3.2用高斯消元法和 MATLA叁序求解下面的非齊次線性方程組，并且用逆矩陣解方程組的方法驗證.x1 x2 + x3 3x4 =1, A=1 -1 1 -3; 0 -1 -1 1;2 -2 -4 6;1 -2 -4 1; b=1;0; -1;-1; RA,RB,n,X

7、 =gaus (A,b)運行后輸出結(jié)果RA =4 RB =4 n =4請注意：因為RA=RB=n ,所以此方程組有唯一解X =0-0.50000.50000列主元消元法及其MATLAB?序用列主元消元法解線性方程組AX= b的MATLA翼序function RA,RB,n,X=liezhu(A,b)B=A b; n=length(b); RA=rank(A);RB=rank(B);zhica=RB-RA;26 /13if zhica0,disp(請注意：因為RA=RB,所以此方程組無解.,) returnendif RA=RBif RA=ndisp(請注意：因為RA=RB=n ,所以此方程組有

8、唯一解.,)X=zeros(n,1); C=zeros(1,n+1);for p= 1:n-1Y,j=max(abs(B(p:n,p); C=B(p,:);B(p,:)= B(j+p-1,:); B(j+p-1,:)=C;for k=p+1:nm= B(k,p)/ B(p,p);B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);end endb=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); endelsedisp(請注意

9、：因為RA=RB A=0 -1-1 1;1 -1 1 -3;2 -2 -4 6;1 -2 -4 1;b=0;1;-1;-1; RA,RB,n,X=liezhu(A,b)運行后輸出結(jié)果請注意：因為RA=RB=n ,所以此方程組有唯一解.RA = 4 , RB = 4 , n = 4 , X =0 -0.5 0.5 03.4 LU分解法及其MATLA能序判斷矩陣LU分解的充要條件及其MATLAB?#判斷矩陣A能否進行LU分解的MATLA翼序function hl=pdLUfj(A)n n =size(A); RA=rank(A);if RA=ndisp(請注意：因為A白n階行列式hl等于零，所以A

10、不能進行LU分解.A的秩 RA如下：),RA,hl=det(A);returnendif RA=nfor p=1:n,h(p)=det(A(1:p, 1:p);,endhl=h(1:n);for i=1:nif h(1,i)=0disp(請注意：因為A白r階主子式等于零，所以 A不能進行LU分解.A的秩RA和各階順序主子式值 hl依次如下：,),hl;RA, returnend27 /13endifh(1,i)=0disp(請注意：因為A的各階主子式都不等于零，所以 A能進行LU分 解.A的秩RA和各階順序主子式值 hl依次如下：)hl;RA endend123、12 3、1(1)1127；1

11、2 7；(3)1d 5 6Jd 5 6;d例 3.4.1判斷下列矩陣能否進行 LU分解，并求矩陣的秩解 (1)在MATLAB作窗口輸入程序2 3 2 35 6 , A=1 2 3;1 12 7;4 5 6;hl=pdLUfj(A)運行后輸出結(jié)果為請注意：因為A的各階主子式都不等于零，所以 A能進行LU分解.A的秩RA和 各階順序主子式值hl依次如下：RA = 3 , hl = 1 10 -48(2)在MATLAS：作窗口輸入程序 A=1 2 3;1 2 7;4 5 6;hl=pdLUfj(A)運行后輸出結(jié)果為請注意：因為A的r階主子式等于零，所以 A不能進行LU分解.A的秩RA和各階 順序主子

12、式值hl依次如下：RA = 3 , hl =1 0 12(3)在MATLAS：作窗口輸入程序 A=1 2 3;1 2 3;4 5 6;hl=pdLUfj(A)運行后輸出結(jié)果為請注意：因為A的n階行列式hl等于零，所以A不能進行LU分解.A的秩RA如下RA = 2 , hl = 0直接LU分解法及其MATLA醒序?qū)⒕仃嘇進行直接LU分解的MATLA叁序function hl=zhjLU(A)n n =size(A); RA=rank(A);if RA=ndisp(請注意：因為A白n階行列式hl等于零，所以A不能進行LU分解.A 的秩 RA如下:),RA,hl=det(A);returnendif

13、 RA=nfor p=1:nh(p)=det(A(1:p, 1:p);endhl=h(1:n);for i=1:nif h(1,i)=0disp(請注意：因為A的r階主子式等于零，所以A不能進行LU分解.A 的秩RA和各階順序主子式值hl依次如下：),hl;RAreturnendendifh(1,i)=0disp(請注意：因為A的各階主子式都不等于零，所以A能進行LU分解.A的秩RA和各階順序主子式值hl依次如下：)28 / 13for j=1:nU(1,j尸A(1,j);endfor k=2:nfor i=2:nfor j=2:nL(1,1)=1;L(i,i)=1; if ijL(1,1)=

14、1;L(2,1)=A(2,1)/U(1,1);L(i,1)=A(i,1)/U(1,1);L(i,k)=(A(i,k)- L(i,1:k-1)*U(1:k-1,k)/U(k,k); elseU(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j); end end end end hl;RA,U,L end end例 3.4.3用矩陣進行直接LU分解的MATLAB程序分解矩陣I_ 0_ 1012120014 30 3解 在MATLAB：作窗口輸入程序 A=1 0 2 0;0 1 0 1;1 2 4 3;0 1 0 3; hl=zhjLU(A)運行后輸出結(jié)果請注意：因為A的各階主子式

15、都不等于零，所以A能進行LU分解.A的秩RA和各階順序主子式值hl依次如下：RA = 4L = 1 0 0 0 TOC o 1-5 h z U = 1 0 2 001000 1 0 112100 0 2 10101hl = 1 1 2 4判斷正定對稱矩陣的方法及其 MATLAB?序 判斷矩陣A是否是正定對稱矩陣的 MATLA叁序function hl=zddc(A)n n =size(A);for p=1:nh(p)=det(A(1:p, 1:p);endhl=h(1:n);zA=A; for i=1:nif h(1,i)0disp(請注意：因為A的各階順序主子式hl都大于零，所以A是正定的.

16、A的轉(zhuǎn)置矩陣zA和各階順序主子式值hl依次如下：,)hl;zA end例3.4.5判斷下列矩陣是否是正定對稱矩陣：29 / 13(1)0.11122213-313844419(2)1121130與1、11100-43；(3)點 11 了-600190011不飛11002209-6C21u1-60解 (1)在MATLAB作窗口輸入程序 A=0.1 2 3 4;-1 2 -3 4;11 21 13 41;5 7 8 9;hl=zddc (A)運行后輸出結(jié)果請注意：A不是對稱矩陣請注意：因為A的各階順序主子式 hl不全大于零，所以A不是正定的.A的轉(zhuǎn) 置矩陣zA和各階順序主子式值 hl依次如下：=1

17、/10-111222173-3138444191/1011/5-1601/1053696/5zAhl =因此，A即不是正定矩陣，也不是對稱矩陣.(2)在MATLAB：作窗口輸入程序 A=1 -1 2 1;-1 3 0 -3;2 0 9 -6;1 -3 -6 19,hl=zddc(A)運行后輸出結(jié)果A = 1-121-130-3209-61-3-619請注意：A是對稱矩陣請注意：因為A的各階順序主子式hl都大于零，所以A是正定的.A的轉(zhuǎn)置矩陣zA和各階順序主子式值hl依次如下：zA = 1-121-130-3209-61-3-619hl = 12624(3)在MATLAB：作窗口輸入程序 A=1

18、/sqrt(2) -1/sqrt(2) 0 0; -1/sqrt(2) 1/sqrt(2) 0 0; 00 1/sqrt(2) -1/sqrt(2); 0 0 -1/sqrt(2) 1/sqrt(2), hl=zddc (A)運行后輸出結(jié)果A= 985/1393 -985/139300-985/13939851393 -985/139300-985/1393 985/1393請注意：A是對稱矩陣請注意：因為A的各階順序主子式 hl不全大于零，所以A不是正定的.A的轉(zhuǎn)置 矩陣zA和各階順序主子式值 hl依次如下： TOC o 1-5 h z zA = 985/1393

19、-985/139300-985/13939851393 -985/139300-985/1393985/1393hl = 985/1393000可見，A不是正定矩陣，是半正定矩陣；因為 A = A T因此，A是對稱矩陣.(4)在MATLAB：作窗口輸入程序 A=-2 1 1;1 -6 0;1 0 -4;hl=zddc (A)運行后輸出結(jié)果A = -2111 -6030 / 1310 -4請注意：A是對稱矩陣請注意：因為A的各階順序主子式 hl不全大于零，所以A不是正定的.A的轉(zhuǎn)置 矩陣zA和各階順序主子式值 hl依次如下：zA = -211 hl = -2 11 -3

20、81 -6010 -4可見A不是正定矩陣，是負(fù)定矩陣；因為 A = A T因此，A是對稱矩陣.3.5 求解線性方程組的LU方法及其MATLA程序解線性方程組的楚列斯基(Cholesky )分解法及其MATLAB?序例3.5.1 先將矩陣A進行楚列斯基分解，然后解矩陣方程 AX = b,并用其他方法 驗證.1-121、,1、-130-33A =209-6,b =5-3-619I7解在工作窗口輸入A=1 -1 2 1;-1 3 0 -3; 2 0 9 -6;1 -3 -6 19;b1=1:2:7; b=b1; R=chol(A);C=A-R*R,R1=inv(R);R2=R1;x=R1*R2*b,

21、Rx=Ab-x運行后輸出方程組的解和驗證結(jié)果3.5.2 解線性方程組的直接LU分解法及其MATLAB?序例3.5.2 首先將矩陣A直接進行LU分解，然后解矩陣方程 AX = b101001210、1、1r2b =3-13 A=10A直接進行2 0;01LU分解.在MATLAB作窗口輸入程序0 1;12 4 3;01 0 3;b=1;2;-1;5;hl=zhjLU(A),A-L*U 運行卡輸出LU分解請注意：和各階順序主子式值RA = 4U = 1000因為A的各階主子式都不等于零, hl依次如下：所以A能進行LU分解.A的秩RA010020200112L = 1010hl = 10121100

22、10200014x = Rx = 1.0e-014 * C = 1.0e-015 *-8.0000-0.710500000.3333-0.08330-0.4441003.66670.222000002.00000.13320000U的積A = LU .A分解為一個單位下三角形矩陣L和一個上三角形矩陣(2)在工作窗口輸入 U=1 0 2 0;0 1 0 1;0 0 2 1;0 0 0 2; L=1 0 0 0;0 1 0 0;1 2 1 0;0 1 0 1;31 / 13b=1;2;-1;5;U1=inv(U); L1=inv(L); X=U1*L1*b,x=Ab運行后輸出方程組的解X = 8.

23、500000000000000.50000000000000-3.750000000000001.500000000000003.5.3 解線性方程組的選主元的 LU方法及其MATLA勰序例3.5.3先將矩陣A進彳T LU分解，然后解矩陣方程 AX =b其中0.1-1A =11I5解方法1根據(jù)(3.28 )式編寫34、1-342b =1341-189 322217MATLAB序，然后在工作窗口輸入 A=0.1 2 3 4;-1 2 -3 4;11 21 13 41;5 7 8 9; b=1;2;-1;5;L U P=LU(A), U1=inv(U); L1=inv(L); X=U1*P = 0

24、010010010000001X =-1.20133.36770.0536-1.4440 TOC o 1-5 h z L1*P*b 運行后輸出結(jié)果 L = 1.0000000-0.0909 1.0000000.0091 0.4628 1.000000.4545 -0.6512 0.2436 1.0000 U =11.0000 21.0000 13.0000 41.0000 03.9091-1.8182 7.7273003.7233 0.0512000 -4.6171方法2根據(jù)(3.29)式編寫MATLA程序，然后在工作窗口輸入 A=0.1 2 3 4;-1 2 -3 4;11 21 13 41

25、;5 7 8 9;b=1;2;-1;5; F U=LU(A), U1=inv(U); F1=inv(F); X=U1*F1*bF=0.0091 0.4628 1.00000-0.0909 1.0000001.00000000.4545 -0.6512 0.2436 1.0000X =-1.2013 3.3677 0.0536 -運行后輸出結(jié)果U=11.0000 21.0000 13.0000 41.00000 3.9091 -1.8182 7.727300 3.7233 0.0512000 -4.61711.4440用LU分解法解線T方程組 A m必X = b的MATLA翼序function

26、RA,RB,n,X,Y=LUjfcz(A,b)n n =size(A);B=A b; RA=rank(A); RB=rank(B);for p=1:nh(p尸det(A(1:p, 1:p);endhl=h(1:n);for i=1:nif h(1,i)=0disp(請注意：因為A的r階主子式等于零，所以A不能進行LU分解.A 的秩R解口各階順序主子式值hl依次如下：)hl;RA returnend endif h(1,i)=0disp(請注意：因為A的各階主子式都不等于零，所以A能進行LU分解.A的秩RA和各階順序主子式值hl依次如下：)X=zeros(n,1); Y=zeros(n,1);

27、C=zeros(1,n);r=1:1;32 /13for p=1:n-1max1,j=max(abs(A(p:n,p); C=A(p,:);A(p,:)= A(j+P1,:); C= A(j+P1,:);g=r(p); r(p)= r(j+P1); r(j+P1)=g;for k=p+1:nH= A(k,p)/A(p,p); A(k,p) = H; A(k,p+1:n)=A(k,p+1:n)- H* A(p,p+1:n);endend丫(1)=B(r(1);for k=2:nY(k)= B(r(k)- A(k,1:k-1)* Y(1:k-1);endX(n)= Y(n)/ A(n,n);for

28、 i=n-1:-1:1X(i)= (Y(i)- A(i, i+1:n) * X (i+1:n)/ A(i,i);end end RA,RB,n,X,Y；3.6 誤差分析及其兩種MATLA能序用MATLA歆件作誤差分析例3.6.2解下列矩陣方程 AX =b,并比較方程（1）和（2）有何區(qū)別，它們的解有何變化.其中11 / 21 / 31 /41 / 51 / 61 /711/21 /31/41 /51 / 61 /71 / 821/ 31 /41/51 /61 / 71 /81 /92(1) A =1/41 /51/61 /71 /81 /91 /10,b 二21/ 51 / 61 / 71 /

29、81 / 91/101/1121/ 61 / 71 / 81 /91 /101 / 111 /122J/ 71 /81 / 91/101/111/121 / 13,a1I .0011 / 21 / 31 / 41 / 51 / 61 / 7J1/21 / 31 /41/51 / 61 / 71 /821 / 31 / 41 /51 / 61 / 71 /81 / 92(2)A =1 / 41 / 51 / 61 / 71 / 81 / 91 / 10,b =21 / 51 / 61 /71 / 81 / 91 / 101 / 1121 / 61 / 71 /81 / 91 /101/111 /

30、 122、1/71 / 8 1 / 91 / 101/111 / 121 /13，I2解(1)矩陣方程AX =b的系數(shù)矩陣A為7階希爾伯特(Hilbert )矩陣，我們可以用下列命令計算n階希爾伯特矩陣h=hilb(n)% 輸出 h 為n 階 Hilbert 矩陣在MATLABL作窗口輸入程序 A=hilb(7);b=1;2;2;2;2;2;2;X=Ab運行后輸出 AX=b 的解為 X = (-35 504 -1260-420020790 -2772012012 )T .(2)在MATLAB：作窗口輸入程序 B =0.001,zeros(1,6);zeros(6,1),zeros(6,6);A

31、=(B+hilb(7); b=1;2;2;2;2;2;2;X=Ab運行后輸出方程的解為X= (-33 465 -966 -5181 22409 -29015 12413)T.在MATLABL作窗口輸入程序 X =-33, 465,-966,-5181,22409,-29015,12413;33 /13X1 =-35,504,-1260,-4200,20790,-27720,12012,; wu=X1- X運行后輸出方程(1)和(2)的解的誤差為X =(-2 39 -294 981 -1619 1295 - 401 )T. 一 i6.001萬程(1)和(2)的系數(shù)矩陣的差為 A=，常數(shù)向量相同，

32、則Ax = b。6貨。6源J的解的差為 舉=(2 39 -294 981 -1619 1295 -401 )T. A的微小變化， 引起X的很大變化，即 X對A的擾動是敏感的.求P條件數(shù)和討論AX =b解的性態(tài)的MATLAB序求P條件數(shù)和討論AX =b解的性態(tài)的MATLA翼序function Acp=zpjxpb(A)Acw = cond (A, inf);Ac1= cond (A,1);Ac2= cond (A,2);Acf = cond (A,fro );dA=det(A);if (Acw50)&(dA Acp =zpjxpb(hilb(7); Acp,det(hilb(7)運行后輸出結(jié)果請

33、注意：AX=b是病態(tài)的，A的8條件數(shù),1條件數(shù),2條件數(shù)，弗羅貝尼烏斯條 件數(shù)和A的行列式的值依次如下：ans = 1.0e+008 * 9.8519 9.8519 4.7537 4.8175 0.0000ans = 4.8358e-025(2)在 MATLAB：作窗口輸入程序 A=2 3 -1 5;3 1 2 -7;4 1 -3 6;1 -2 4 -7;Acp=zpjxpb(A); Acp運行后輸出結(jié)果AX=b是良態(tài)的，A的8條件數(shù),1條件數(shù),2條件數(shù)，弗羅貝尼烏斯條件數(shù)和A的行列式的值依次如下：ans =14.1713 19.4954 8.2085 11.4203 327.0000用P范數(shù)

34、討論AX = b解和A的性態(tài)的MATLAB?序用P范數(shù)討論 AX =b解和A的性態(tài)的MATLA翼序function Acp=zpjwc(A,jA,b,jb,P)Acp = cond (A,P);dA=det(A); X=Ab;dertaA=A-jA;PndA=norm(dertaA, P);dertab=b-jb;Pndb=norm(dertab,P);34 / 13if Pndb0 jX=Ajb; Pnb= norm(b, P);PnjX = norm(jX,P); dertaX=X-jX; PnjdX= norm(dertaX, P);jxX= PnjdX/PnjX; PnjX =norm

35、(jX,P);PnX = norm(X,P); jxX= PnjdX/PnjX; xX= PnjdX/PnX;Pndb=norm(dertab,P);xAb=Pndb/Pnb;Pnbj=norm(jb,P); xAbj=Pndb/Pnbj;Xgxx= Acp*xAb; end if PndA0jX=jAb; dertaX=X-jX;PnX = norm(X,P); PnjdX= norm(dertaX, P);PnjX = norm(jX,P); jxX= PnjdX/PnjX;xX= PnjdX/PnX;PnjA=norm(jA,P); PnA=norm(A,P);PndA=norm(der

36、taA,P);xAbj= PndA/PnjA;xAb= PndA/PnA; Xgxx= Acp*xAb; end if (Acp 50)&(dA jA =1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769; A=hilb(7);b=1;1/3;4;2;2;2;2;jb=1;0.3333;4;2;2;2;2; Acp=zpjwc(A,jA,b,jb,inf)運行后輸疝結(jié)果請注意：AX=b是病態(tài)的，