Ch3-線性代數

1、83-90第3章線性代數 應用MATLAB函數，可以求解數學問題。在這一章中，我們主要介紹MATLAB在線性代數中的應用，其內容包括：行列式的求值，矩陣的基本計算，矩陣的初等變換和矩陣的秩，矩陣的分解，求解線性方程組，向量的內積和正交，方陣的特征值和特征向量，矩陣的對角化，求解二次型的標準形，以及判別二次型的正定性。3-1 行列式的求值 在線性代數中，行列式是一個基本工具，其應用比較廣泛。在MATLAB中我們只需要借助函數det 就可以求出行列式的值，其格式為：det(A)其中A為n階方陣。范例3-1求矩陣A=的行列式的值。程序設計：>> clear>> A=1 0 2

2、 1;-1 2 2 3;2 3 3 1;0 1 2 1;>> det(A)運行結果：ans = 14程序說明：1 clear命令的作用是清除內存中的變量。2 矩陣的輸入可以有兩種格式，除程序中的輸入方式外，還可以如下輸入：A=1，0，2，1；-1，2，2，3；2，3，3，1；0，1，2，1例題分析： 在線性代數中，我們可以用det(A)的值，來判斷矩陣是否可逆，即矩陣是否奇異（可逆矩陣即為非奇異矩陣）。但是，在MATLAB中，由于計算機運算存在舍入誤差，用det(A)=0判斷A是否可逆，并不是完全有效；如用行列式的絕對值abs(det(A)<tol判斷時，tol的值又不好確定

3、。所以，我們常根據矩陣的條件數cond(A)判斷矩陣是否奇異（關于矩陣的條件數可以參看3-2-11節(jié)）。范例3-2計算行列式程序設計：>> clear>> syms a b c d %聲明變量>> A=a 1 0 0;-1 b 1 0;0 -1 c 1;0 0 -1 d; %生成符號矩陣>> DA=det(A)運行結果：DA =a*b*c*d+a*b+a*d+c*d+1程序說明：函數det也可以用于計算含有變量的行列式。3-2 矩陣的基本運算 矩陣的基本運算主要介紹矩陣的加、減，數與矩陣相乘，矩陣與矩陣相乘，矩陣相除，矩陣的冪運算，矩陣的轉置和共

4、軛，矩陣的逆和偽逆，矩陣的跡，矩陣和向量的范數以及矩陣的條件數。矩陣的基本運算命令如表3-1所示。表3-1矩陣的基本運算命令列表命令意義命令意義+矩陣相加conj矩陣的共軛-矩陣相減inv矩陣的逆*矩陣相乘pinv矩陣的偽逆矩陣左除trace矩陣的跡/矩陣右除norm矩陣的范數矩陣的冪norm向量的范數矩陣的轉置cond矩陣的條件數說明：表中僅列出具有代表性的部分命令，關于該部分內容的命令的格式及使用在下面的章節(jié)中都有介紹。3-2-1 矩陣的加、減 進行加減運算的矩陣，要求維數相同，即行數和列數都分別相同。進行加、減運算時，矩陣相應位置的元素相加、減。范例3-3求矩陣A=與矩陣B=的和與差。程

5、序設計：>> clear>> A=1 2 3;2 1 2;3 3 1;>> B=3 2 4;2 5 3;2 3 1;>> C=A+B;>> D=A-B;>> C,D運行結果：C = 4 4 7 4 6 5 5 6 2D = -2 0 -1 0 -4 -1 1 0 0例題分析：1. 進行加減運算的矩陣必須是同型的。2. 在進行矩陣相加的運算時，A+B和B+A的值相同，滿足加法交換率。3-2-2數與矩陣相乘數與矩陣相乘，是數與矩陣中的每個元素相乘。范例3-4求矩陣A=與5的乘積。程序設計：>> clear>&

6、gt; A=1 0 1;2 1 1;1 2 1;>> B=5*A>> C=A*5運行結果：B = 5 0 5 10 5 5 5 10 5C = 5 0 5 10 5 5 5 10 5程序說明：5*A與A*5的值相同。3-2-3 矩陣與矩陣相乘設是A一個m×n的矩陣，B是一個n×s矩陣，A= B=則A與B的乘積C（）是m×n的矩陣，且 在線性代數中，我們就是按照上面的公式求得矩陣與矩陣的乘積，同時，要注意：兩矩陣相乘時，第一個矩陣（左矩陣）的列數必須等于第二個矩陣（右矩陣）的行數。范例3-5求矩陣A=與矩陣B=的乘積。程序設計：>>

7、; clear>> A=1 2 3;2 1 2;3 3 1;>> B=3 2 4;2 5 3;2 3 1;>> C=A*B>> D=B*A運行結果：C = 13 21 13 12 15 13 17 24 22D = 19 20 17 21 18 19 11 10 13程序說明：比較C和D，可以看出A*B和B*A的結果完全不同。3-2-4 矩陣相除 在MATLAB中，矩陣相除可以運用運算符“”（左除）和“”（右除），而在線性代數中并沒有定義矩陣的除法，下面我們通過例子看看矩陣的左除和右除。范例3-6求矩陣A=與矩陣B=相除。程序設計：>>

8、; clear>> A=1 2 3;4 2 1;2 1 3;>> B=2 1 2;1 2 1;3 2 1;>> C=AB %矩陣左除，相當于inv(A)*B,inv(A)為矩陣A的逆C = 0.3333 0.6000 -0.2000 -0.6667 -0.4000 0.8000 1.0000 0.4000 0.2000>> D=A/B %矩陣右除，相當于A*inv(B)D = 1.3333 1.3333 -1.0000 0 -0.5000 1.5000 1.6667 0.1667 -0.5000 程序說明：1. 矩陣的左除和右除概念完全不同，要注

9、意區(qū)分。2. 可以利用矩陣的左除求解線性方程組AX=b，其中X= Ab。3. 可以利用矩陣的右除求解線性方程組XA=b，其中X=A/b。錯，應為：X=b/A范例3-7求矩陣A=與矩陣B=相除。程序設計：>> A=1 2 3;4 5 6;7 8 9;>> B=9 8 7;6 5 4;3 2 1;>> C=ABWarning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.541976e-018.(Type "warning off M

10、ATLAB:nearlySingularMatrix" to suppress this warning.)C = -27.0000 -26.0000 -17.0000 42.0000 41.0000 24.0000 -16.0000 -16.0000 -8.0000>> D=A/BWarning: Matrix is singular to working precision.(Type "warning off MATLAB:singularMatrix" to suppress this warning.)D = Inf Inf Inf Inf

11、Inf Inf Inf Inf Inf程序說明：在線性代數中,我們可以用矩陣行列式的值來判斷矩陣是否奇異，現在我們手算一下A和B的行列式值，|A|=|B|=0，在MATLAB中計算A和B的行列式值:>> detA=det(A)detA = 0>> detB=det(B)detB = 0由計算結果可以看出兩種方法的計算結果相同，A和B都是奇異矩陣。但我們不建議使用矩陣的行列式值來判斷矩陣是否奇異，因為當行列式的值很小時，由于計算機的舍入誤差可能求出行列式的值為零，而無法得出正確的結論。3-2-5 矩陣的冪運算進行冪運算的矩陣必須是方陣，當冪指數p為正整數時，冪運算即為矩陣

12、A的自乘運算，p為矩陣自乘的次數；當冪指數p為負整數時，冪運算即為的自乘運算，-p為矩陣自乘的次數。當冪指數p為非整數時，AP=V*。其中V為方陣A的特征向量，D=為方陣A的特征值對角矩陣。 標量的冪運算定義為：pA= V*,其中V，D由矩陣A的特征值分解而得，V為矩陣A的特征向量矩陣，D為矩陣A的特征值矩陣。范例3-8求矩陣A=的三次方。程序設計：>> clear>> A=1 0 1;2 1 1;1 2 1;>> B=A3運行結果：B = 8 6 6 15 11 12 18 12 14程序說明：也可以利用矩陣的自乘求解本題，程序如下：>> A=

13、1 0 1;2 1 1;1 2 1;91-98>>C=A*A*AC=8 6 615 11 1218 12 14兩種方法運算的結果相同。范例 3-9 求5 的 次冪。程序設計： >> clear>> A=1 2;2 1;>> B=5AB = 62.6000 62.4000 62.4000 62.6000>> %現在用特征值分解求5A 關于矩陣的特征值分解可以參看3-4-5節(jié) %求矩陣A的特征向量和特征值矩陣>> V,D=eig(A)V = -0.7071 0.7071 0.7071 0.7071D = -1 0 0 3>

14、;> C=V*(5D)*V(-1)C = 62.6000 62.4000 62.4000 62.6000程序說明：1 函數eig（A）所求得矩陣V即為特征向量矩陣，D為特征值矩陣，特征值-1對應的特征向量為-0.7071 0.7071，特征值3對應的特征向量為0.7071 0.7071，特征值對應的特征向量為矩陣V的對應列向量。2 比較B和C可以發(fā)現，它們的值相等。3-2-6 求矩陣的轉置矩陣 矩陣的轉置是指把矩陣的換成同序數的列從而得到新矩陣。它也算一種基本運算，可以用運算符“”來求矩陣的轉置，通過判斷矩陣與其轉置矩陣的關系可以判斷矩陣的對稱性。 范例 3-10求矩陣A=的轉置矩陣A。

15、程序設計：>> clear>> A=1,0,1,1;2,1,0,1;2,1,2,1;1,0,2,1;>> A'運行結果：ans =1 2 2 1 0 1 1 0 1 0 2 2 1 1 1 1范例 3-11判斷矩陣A=是否對稱。程序設計：>>clear>>A=1 2;2 1;>>B=A;>>if(A= =B)Fprintf(A是對成矩陣)else if(A= =(-B)fprintf(A是反對稱矩陣)elseprintf(A既不是對稱矩陣也不是反對稱矩陣)endend運行結果:A是對稱矩陣例題分析:1

16、若A=A,則矩陣為對稱矩陣;若A+A=0,則矩陣為反對稱矩陣.2 也可以利用函數isequal來判斷矩陣是否對稱,若isequal(A,A)=1,則矩陣為對稱矩陣, 若isequal(-A,A)=1,則矩陣A為反對稱矩陣。范例 3-12求復數矩陣A=的轉置矩陣。程序設計：>> clear>> A=1+2*i 2;1 i;>> B=A'運行結果：B =1.0000 - 2.0000i 1.0000 2.0000 0 - 1.0000i程序說明：1 從運行結果可以看出，復數轉置矩陣是由原矩陣轉置后的共軛復數組成。2 如果僅要求將復數矩陣轉置，可以命令B=

17、A. >> A=1+2*i 2;1 i;>> B=A.'B =1.0000 +2.0000i 1.0000 2.0000 0 +1.0000i3-2-7 求矩陣的共軛矩陣 復數矩陣的共軛與復數的共軛類似,復數矩陣的共軛矩陣與復數矩陣的實部相同，虛部相反，可運用函數conj（A）求得。 范例 3-13 求矩陣A=和B=的共軛矩陣。程序設計：>> clear>> A=1 2;3 4;>> B=1+2i 2+2i;1-i,3+i;>> C=conj(A);>> D=conj(B);>> C,D運行

18、結果：C =1 2 3 4D =1.0000 - 2.0000i 2.0000 - 2.0000i 1.0000 + 1.0000i 3.0000 - 1.0000i3-2-8 求矩陣的逆和偽逆 如果矩陣A是方陣且是非奇異的，我們可以用函數inv（A）求得逆矩陣。如果A不是方陣或A是奇異陣，則可以用函數pinv（A）求得矩陣A的偽逆。偽逆只是有逆的某些性質，與逆不同，它們的格式如下： X=inv（A） X為A的逆 X=pinv（A） X為A的偽逆，其求解是建立在奇異值分解和把小于默認誤差的奇異值當作0的基礎上計算的，且滿足A*X*A=A，X*A*X=X且A*X和X*A都是Hermitian矩陣

19、 X=pinv（A，tol） tol為指定的誤差范例 3-14求矩陣A=和矩陣B=的逆矩陣。程序設計：>> clear>> A=1 2 1;2 4 2;2 1 1;>> B=3 1 2;1 2 2;3 1 4;運行結果：>> C=inv(A)Warning: Matrix is singular to working precision.(Type "warning off MATLAB:singularMatrix" to suppress this warning.)C =Inf Inf Inf Inf Inf Inf I

20、nf Inf Inf>> D=inv(B)0.6000 -0.2000 -0.2000 0.2000 0.6000 -0.4000 -0.5000 0 0.5000程序說明：1 如果矩陣不可逆，則運行結果會給出警告信息，說明矩陣是奇異的。2 矩陣的逆也可以用表達式A(-1)求得.3 如果矩陣A可逆,則可以利用系數矩陣A的逆求解線形方程組AX=b,其中 X=inv(A)*b。范例 3-15求上例中矩陣A的偽逆。程序設計：>> A=1 2 1;2 4 2;2 1 1;>> C=pinv(A)運行結果：C =-0.0727 -0.1455 0.6364 0.127

21、3 0.2545 -0.36360.0182 0.0364 0.0909 程序說明： 如果A是奇異陣或不是方陣，則可以利用矩陣A的偽逆求解線形方程組AX=b，其中X=pinv（A）*b，如果方程組有精確解，則求得的是方程組的精確解；如果方程組的解不唯一，則求得的是方程組的最小二乘意義上的解。3-2-9 矩陣的跡 矩陣的跡是指矩陣的對角線元素的和，也等于矩陣的特征值的和，用函數trace來求。范例3-16 求矩陣A=的對角線元素之和。程序設計：>> clear>> A=1 2 1;2 1 2;1 1 2;>> a=trace(A)a = 4>> E

22、=eig(A)E = 4.3028 -1.0000 0.6972>> b=sum(E)b =4.0000程序說明：1E=eig（A）是求矩陣A的特征值列向量2函數sum（E）是求矩陣E的列向量的元素之和，也可以通過sum（A，dim）確定矩陣是行元素相加還是列元素相加，1為列元素相加，2為行元素相加，默認值為1，如： >>sum(A,1) ans= 4 4 5 >>sum(A,2) Ans= 4 5 43 矩陣的對角線元素之和等于矩陣的特征值之和,a=b。3-2-10 矩陣和向量的范數 范數是與距離相類似的概念，不過范數是建立在線形空間上的，它利用了線形空間

23、的代數結構，本質是描述線形空間中元素的“長度”和“大小“，它也能表示元素之間的距離。對于矩陣而言，其范數為=。在MATLAB中，函數norm可以求矩陣的范數，其格式如下：norm（A） 矩陣A的最大奇異值， max（svd（A）norm（A，2） 矩陣A的2-范數，等同于norm（A）norm（A，1） 矩陣A的1-范數（列范數），其值為max（sum（abs（A）norm（A，inf） 矩陣A的無窮范數（行范數），其值為max（sum（abs（A）norm（A，fro） 矩陣A的 F范數，其值為sqrt（sum（diag（A*A）向量的范數可以表示為，我們可以用函數norm求出向量的范數，其

24、格式如下：norm（V，P） 向量的P范數，值為sum（abs（V）.P）(1/P)norm（V） 向量的2-范數，值為sqrt(sum(abs(V).2)norm（V，1） 向量的1-范數，值為sum（abs（V）norm（V，inf） 向量的+范數，值為max(abs(V)norm（V，-inf） 向量的-范數，值為min(abs(V)范例 3-17 求矩陣A=的列范數，2-范數，,行范數, F范數。程序設計：>> clear>> A=2 3 1;-4 2 2;3 4 3;>> norm(A,1) %矩陣A的列范數ans = 9>> norm

25、(A) %矩陣A的2-范數ans = 6.8834>> norm(A,inf) %矩陣A的行范數ans = 10>> norm(A,'fro') %矩陣A的F范數ans =8.4853程序說明： 在求矩陣的2-范數時，也可以用norm（A，2）范例 3-18 求向量V=1 2 -3的1-范數，2-范數和無窮范數。程序設計：>> V=1 2 -3;>> norm(V) %向量V的2-范數ans =3.7417>> norm(V,1) %向量的1-范數ans =6>> norm(V,inf) %向量的+范數an

26、s =3>> norm(V,-inf) %向量V的-范數ans =13-2-11 矩陣的條件數 在MATLAB中我們可以用函數cond 、rcond、 condest、 condeig來求出矩陣的條件數。根據矩陣的條件數的值，我們可以判斷矩陣是否“病態(tài)“或”奇異“。病態(tài)矩陣和奇異矩陣起概念不同：例如，在求解AX=b時，如果A的微小變化會引起解的劇烈變化，則矩陣A即為病態(tài)矩陣。而奇異矩陣是指行列式為0的矩陣，等價于可逆矩陣。函數具體格式有以下幾種形式：c=cond（A） 返回矩陣的2-范數條件數，值為矩陣A的最大奇異值與最小奇異值的比，值越大，說明矩陣越接近奇異c=cond（A，p）

27、 返回矩陣的p-條件范數，值為norm（A，p）*norm（inv（A），p），其中p取值1、2、inf或froc=rcond（A） 返回矩陣的條件數的倒數值，如果矩陣A是良態(tài)的，返回值接近于1；如果矩陣A是病態(tài)的, 返回值接近于0c=condest（A） 返回方陣的1-范數的條件數下界值c，v= condest（A） 返回的向量v滿足norm（A*v，1）=norm（A，1）*norm（v，1）/cc=condest（A，t） t為矩陣的列數的正整數，t的默認值為2condeig（A） 返回矩陣A的特征值的條件數向量V，D，s= condeig（A） V為矩陣A的特征向量，D為A的特征值，s

28、為矩陣A的特征值的條件數范例 3-19求解矩陣A=的條件數。程序設計：>>clear>>A=2 -1 1;1 2 2;2 3 1;99-106>> cond(A,1) %矩陣A的1-范數條件數ans = 6>> cond(A,2) % 矩陣A的2-范數條件數ans = 4.3735>> cond(A,inf) %矩陣A的無窮范數條件數ans = 7.0000>> cond(A,'fro') %矩陣A的F范數條件數ans = 5.6051>> rcond(A) %矩陣A條件數的倒數值ans = 0

29、.1667>> condest(A) %矩陣A的1-范數條件數的下界值ans = 6>> V,D,s=condeig(A) %矩陣A的特征值條件數V = -0.0000 -0.6581 -0.3737 -0.7071 0.6226 -0.3950 -0.7071 0.4233 0.8392D = 4.0000 0 0 0 2.3028 0 0 0 -1.3028s = 1.4907 1.51401.0734程序說明：1. 由計算出的cond(A)的值,可以判斷出矩陣A是非奇異的2. 有計算的rcond(A)的值,可以判斷出矩陣是良態(tài)的。下面通過解方程組證明矩陣是良態(tài)的,

30、設b=2 5 6,求解AX=b; >> A=2 -1 1;1 2 2;2 3 1;>> b=2 5 6'>> X1=AbX1 = 1 1 1>> %現將矩陣A的元素都改變0.01,再比較運行結果>>A1=2.01 -0.99 1.01;1.01 2.01 2.01;2.01 3.01 1.01;>> X1=A1bX1 = 0.9901 1.00000.9901比較X和X1的值可以看出,它們看出A的微小變化僅引起解的更微小變化,矩陣A是良態(tài)的。范例3-20判斷方程組的性態(tài)。程序設計:>> A=1 1/2

31、1/3;1/2 1/3 1/4;1/3 1/4 1/5;>> rcond(A)ans = 0.0013>> b=11/6 13/12 47/60'>> X=AbX = 1.0000 1.0000 1.0000 %現將方程組中的A和b都舍入成兩位有效數字變成A1和b1，再求解A1*X=b1。>> A1=1 0.50 0.33;0.50 0.33 0.25;0.33 0.25 0.20;>> b1=1.8 1.1 0.78'>> X1=A1b1X1 = -6.2222 38.2540 -33.6508程序說明:

32、1. 由計算出的rcond(A)的值都可以判斷出矩陣是病態(tài)的。2. 在解方程組AX=b時，如果A或b的微小變化，只引起解的更微小變化，則方程組是良態(tài)的，如果A或b的微小變化導致解的劇烈變化，則方程組是病態(tài)的，在MATLAB中可以根據矩陣A的條件數判斷方程組病態(tài)與否，條件數越大，說明方程組越接近與病態(tài)。3. 當X1和b有微小的變化時，所求得的解X1就變化太大了，因此，該方程組嚴重病態(tài)的。3-3 矩陣的初等變換和矩陣的秩3-3-1 矩陣的初等行變換矩陣的初等行變換包含以下三種變換：1. 對調兩行；2. 以數k0乘某一行中的所有元素；3. 把某一行所有元素的k倍加到另一行對應元素上去。在線性代數中，

33、我們常把矩陣通過行變換化為行最簡形，即非零行向量的第一個元素為1，且含這些元素的列的其他元素為0。利用矩陣的行最簡形，可以求出矩陣的秩、矩陣的逆、向量組的最大無關組等，在MATLAB中使用函數rref或rrefmovie就可以把矩陣化為行最簡形，其格式為：R=rref(A) 給出矩陣A的行最簡形。R，jb=rref(A) jb是一個向量，r=length(jb)是矩陣A的秩，A(：jb)為矩陣A的列向量基，jb表示列向量基的所在列數。R,jb=rref(A,tol) tol為指定的精度。Rrefmovie(A) 給出每一步過程。范列3-21將矩稱A=化為行最簡形。程序設計：>> c

34、lear>> A=1 1 1;1 2 -5;2 3 -4;>> B=rref(A)運行結果：B = 1 0 7 0 1 -6 0 0 0范列3-22利用矩陣的初等行變換求矩陣A=的逆。程序設計：>> B=1 2 3 1 0 0;2 2 1 0 1 0;3 2 2 0 0 1;>> format rat>> C=rref(B)C = 1 0 0 -1/3 -1/3 2/3 0 1 0 1/6 7/6 -5/6 0 0 1 1/3 -2/3 1/3 >> D=C(:,4:6) %取矩陣C的4到6列，D即為矩陣A的逆矩陣D =

35、-1/3 -1/3 2/3 1/6 7/6 -5/6 1/3 -2/3 1/3 例題說明： 由線性代數的知識可知，矩陣A和其同型的單位矩陣E組成增廣矩陣B，對B進行初等行變換，當矩陣A變?yōu)閱挝魂嚂r，單位矩陣E變?yōu)锳的逆。范例3-23求向量組的最大無關組。程序設計：>> clear>> a1=1 2 -2 1'>> a2=2 -3 2 1'>> a3=3 -1 0 2'>> a4=3 2 1 2'>> A=a1 a2 a3 a4A = 1 2 3 3 2 -3 -1 2 -2 2 0 1 1

36、1 2 2 >> R,jb=rref(A)R = 1 0 1 0 0 1 1 0 0 0 0 1 0 0 0 0 jb = 1 2 4 >> A(:,jb)ans = 1 2 3 2 -3 2 -2 2 1 1 1 2 程序說明：由jb=1 2 4，說明矩陣A的第1、2、4列為一個最大不相關組，即a1、a2、a4是向量組的一個最大無關組。3-3-2 矩陣的秩矩陣的秩是指矩陣的行或列向量不相關的個數，可以用函數rank(A)求得。其格式如下：rank(A) 返回矩陣A中比誤差tol大的奇異值個數，tol默認為tol=max(size(A)*norm(A)*epsrank(

37、A,tol)返回矩陣A中比誤差tol大的奇異值個數 范例 3-24求矩陣A=的秩。程序設計：>> clear>> A=2 1 1 2;1 2 2 1;1 2 1 2;2 2 1 1;>> rank(A)運行結果：ans = 4 程序說明：由rank(A)=4可以得知，矩陣A的行向量或列向量完全無關。3-4 矩陣的分解在求解線性方程組、矩陣的特征值和特征向量等過程中都要用到矩陣的分解，在MATLAB中有專門的矩陣分解函數，以方便我們對矩陣進行分解。這里我們主要介紹Cholesky分解、LU分解、奇異值分解、特征值分解、Hessenberg分解和Schur分解。

38、3-4-1 對稱正定矩陣的Cholesky分解 如果矩陣A是對稱正定的，則A可以進行Cholesky分解，即A=RR，其中R為上三角矩陣。在MATLAB中我們用函數chol來求得對稱正定矩陣的Cholesky分解，其格式如下：R=chol(A) 如果矩陣A是對稱正定的，我們就可以得到上三角矩陣R，且R*R=A，在計算求R時，只應用矩陣A的對角和上三角元素；如果不是對稱正定的，那么產生出錯信息。R，q=chol(A) 有兩個輸出，不產生出錯信息。如果A是對稱正定的，p的值為0，R為上三角矩陣；如果A不是對稱正定的，則p的值為正整數，當A是滿陣時，R是階數為q=p-1的上三角矩陣，且R*R=A(1

39、:q,1:q)，當A是稀疏矩陣時，R是q×n的上三角矩陣，且R*R所得矩陣的前q行q列的元素與稀疏矩陣A的對應位置元素相等。范例 3-25 求5階pascal矩陣的Cholesky分解。程序設計：>> clear>> A=pascal(5)A = 1 1 1 1 1 1 2 3 4 5 1 3 6 10 15 1 4 10 20 35 1 5 15 35 70 >> R,p=chol(A)R = 1 1 1 1 1 0 1 2 3 4 0 0 1 3 6 0 0 0 1 4 0 0 0 0 1 p = 0 程序說明：1. pascal(5)是生成5

40、階的pascal矩陣。2. 由p的值可以看出A是對稱正定的。3. 計算R*R的值，可以看出R*R為5階pascal矩陣>> R'*Rans = 1 1 1 1 1 1 2 3 4 5 1 3 6 10 15 1 4 10 20 35 1 5 15 35 70 3-4-2 LU分解 矩陣A的LU分解也稱為三角分解，A=LU，其中，L為下三角矩陣，U為上三角矩陣。當L為單位下上三角矩陣時，稱為Doolittle分解；當U為單位上三角矩陣時，稱為Crout分解。在MATLAB中，我們通過函數lu求得矩陣的LU分解，但是與線性代數中的L、U形式也有些許的區(qū)別，其格式如下：L，U=l

41、u(A) U為上三角矩陣，L為下三角矩陣或為下三角矩陣的變換形式，滿足L*U=AL，U，P=lu(A) L為下三角矩陣，L為上三角矩陣，P為單位矩陣的行變換矩陣，滿足P*A=L*U 在MATLAB中，行列式和矩陣的逆求解就是建立在LU分解基礎之上的： det(A)=det(L)*det(U) inv(A)=inv(U)*inv(L)范例3-26求矩陣A=的LU分解。程序設計：>> clear>> A=2 1 1 2;1 2 3 2;2 4 1 1;3 1 2 3;>> L,U=lu(A)L = 0.6667 0.1000 -0.1200 1.0000 0.3

42、333 0.5000 1.0000 0 0.6667 1.0000 0 0 1.0000 0 0 0U = 3.0000 1.0000 2.0000 3.0000 0 3.3333 -0.3333 -1.0000 0 0 2.5000 1.5000 0 0 0 0.2800>> l,u,p=lu(A)l = 1.0000 0 0 0 0.6667 1.0000 0 0 0.3333 0.5000 1.0000 0 0.6667 0.1000 -0.1200 1.0000u = 3.0000 1.0000 2.0000 3.0000 0 3.3333 -0.3333 -1.0000

43、0 0 2.5000 1.5000 0 0 0 0.2800p = 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 107-114程序說明 ：運用兩種格式分解進行LU分解，所得到的上三角矩陣U和u完全相同，而L和l不同，它們可以通過行變換變?yōu)橄嗤木仃?，關系為L=inv（p）*l。3-4-3 QR分解由線性代數可知，任何實的非奇異矩陣A都可以分解為正交矩陣Q和上三角矩陣R的乘積 ，即A=Q*R，我們可以用函數qr對 矩陣A進行QR分解（也叫正交直角分解），其格式如下：Q，R=qr（A）R為與A同階的上三角矩陣，Q為 正交矩陣，滿足A=Q*RQ，R，E=qr（A）R為對角元素的

44、絕對值遞減的 上三角矩陣，Q為 正交矩陣，E為單位矩陣的變換矩陣，滿足A*E=Q*RQ，R=qr（A，0）產生矩陣A的“經濟型 ”分解，如獲 A 為m×n（m>n）的矩陣，只計算出n列的正交矩陣QQ，R，E=qr（A，0）計算出的E使R的對角元素的 絕對值遞減，滿足Q*R=A（：E）qr（A）得到LAPACK的DEGEQRF或ZGEQRF程序的輸出Triu（qr（A）得到上三角矩陣R范例3-27求矩陣的QR分解。程序設計：>> clear>> A=1 2 3 2;2 1 1 3;3 1 2 1;2 2 1 4;>> Q,R=qr(A)Q =

45、-0.2357 0.7671 -0.5914 -0.0783 -0.4714 -0.1227 0.1428 -0.8616 -0.7071 -0.4603 -0.3671 0.3916 -0.4714 0.4296 0.7036 0.3133R = -4.2426 -2.5927 -3.0641 -4.4783 0 1.8105 1.6877 2.4242 0 0 -1.6621 1.6927 0 0 0 -1.0966>> Q1,R1,E=qr(A)Q1 = -0.3651 0.7303 0.5765 -0.0320 -0.5477 -0.1826 -0.1601 -0.8006 -0.1826 0.5477 -0.8006 0.1601 -0.7303 -0.3651 0.0320 0.5765R1 = -5.4772 -2.7386 -3.4689 -2.9212 0 2.7386 1.2780 1.0954 0 0 -2.0817 0.2562 0 0 0 0