解線性代數(shù)方程組

1、求解線性方程組的直接解法5.3特殊矩陣的三角分解 實對稱矩陣的ldlt分解設(shè)A是實對稱陣，且A的所有順序主子式均不為零，則LDR分解中 R=LT,故可用以作LDLT分解.這就是說，當A的對角元素非零時，我們可 以作LU分解,也就得到LDLT分解丄相同，是單位上三角陣,U的對角元素 構(gòu)成D.不過沒有利用對稱性，存儲量運算量都未能節(jié)省一預(yù)計是一半。試 用n=3的計算表格說明如何實現(xiàn)節(jié)省。d1=u 11=anU12=a12l21=U12/d1U13=a13l31=U13/d1d2=U22=a22-l21U12U23=a23-l21U13l32=U23/d2U33= a33-l31 U13-l 32U

2、23這樣，可用上半部元素逐列計算 D,LT。也可用下半部元素逐行計算L,D引進輔助量ti, t2代替Uij,U2j,并利用對稱性得到：d1=ant1 =a21l21= t1/d1d2= a22-t1121t1 =a31I31=t1/d1t2=a32-t1l21I32=t2/d2d3=a33-t1l31-t2l32據(jù)此不難寫出LDLT分解A=LDLT的計算公式和程序(逐行計算L,D).di=aiifor i=2: nfor j=1:i-1tj = aij-ljltl-lj2t2-lj,j-ltj-1lij=tj/djenddi=aii-liltl-li2t2-li,i-lti-1end存儲約n(

3、n+1)/2單元,乘加運算各約n3/6.利用LDLt分解解Ax= b分四步：1 .分解 A=LDLT2 解 Lg=b求 g3 .解 Dy=g求 y4 .解 L Tx=y 求 x 實對稱正定矩陣的LLt分解A實對稱正定時順序主子式皆正,可作LDLT,D的對角元素皆正，有正的平方根。因此有LLT分解A=LLT,L下三角陣,對角元素皆正,是LDLT 中的LD1/2.乃可用上半部元素逐列計算 LT. 1/2ln=anl21= a12/lnI31=a13/|11l22=(a22-l212)1/2I32=(a23-l21l31)/l2222l33=a33-l31 -l32也可用下半部元素逐行計算L.計算表

4、格和算法安排如下I1/2ln=anl21= a21/ln2 1/2|22=22-|21 )l31= a31/lnI32=(a32-l31l21)/l22I33=(a33-l312-l322) "2lii=aii1/2fori=2: nfor j=1:i-1lij=(aij-lillj 1-li2lj2-li,j-llj,j-1”djjend-l2J22 2lii = (aii - li1 - li2end存儲量，運算量同LDLt分解，但要n次求平方根. 利用LLT分解解Ax=b分三步：1 .分解 A=LLT2 解 Lg=b 求 g3 .解 L Tx=g 求 x 三對角方程組的追趕法消

5、去法或LU分解用于三對角方程組有特殊形式，即稱追趕法設(shè)Ax=f:b1X1+ C1X2=f1aiXi-1 + biXi+ CiXi+ 1 =fii= 2,3, n-1anXn-什 bnXn=fnA是三對角陣，則L,U同樣結(jié)構(gòu).L的對角元素為a 2, a 3,，a n, U的對角 元素為B 1, B 2,B n,上對角元素同A.1 .分解 A=LU: B 1= b1, a i =ai/ B i-1, B i= b-a i Ci-1,i=2,3，,n2. 解 Lg=f 求 g: g1=f1,gi =fi- a ifi-1, i=2,3，,n3. 解 UX=g 求 x: Xn=gn/ B n,x i=

6、(gi-CiXi+1)/ B i,i=n-1,n-2,1編程時，A可用三個一維數(shù)組，f用一個一維數(shù)組 丄,U存入Ao g, x 存入f o還有一種計算格式，消去時用主元素除主行元素，即分解A為下三角 矩陣和單位上三角矩陣之積，相當于對AT作LU分解._biCifj-Pi化)gi _a2b2 c2f2a2%億)g2+ + +aT+ + +m+ +cn _1a+ +g)m-anbnfn 一1 1anPngn 一括號中是單位上三角矩陣的上對角元素計算步驟：1 分解 A=LU: B i=bi，丫 i=cB 1，B i=bi-ai y i-i, y i =G/B i, i=2,3,n 2解 Lg=f 求

7、 g: gi=fi/B i,gi=(fi-aigi-i)/B i, i=2,3，,n3.解 Ux=g 求 x: xn=gn,x i=gi-丫匚Xi+i,i=n-i,n-2,i三對角矩陣是帶形矩陣的特例所謂帶形矩陣是那些主對角線附近幾 條對角線以外元素皆零的矩陣，即aij工0,僅當-mivj-i<m2.帶形矩陣的LU 分解也保持結(jié)構(gòu)5.4 向量和矩陣的范數(shù)引入實數(shù)的絕對值和復(fù)數(shù)的模(也稱絕對值)來表示實數(shù)和復(fù)數(shù)的”大小”從而帶來許多用處.例如,數(shù)列收斂的概念就是通過絕對值來表示的.范 數(shù)這個概念就是這些表示”大小”的數(shù)值普遍化.它在研究數(shù)值計算方法的 收斂性和穩(wěn)定性中有著重要的應(yīng)用. 向量

8、的范數(shù)定義1.如果向量xe Rn(或Cn)的某個實值函數(shù)N(x)=|XI，滿定理1. Cn中任意兩種向量范數(shù)I足條件：1.正定性:I x I> 0, I x I =0iffx=02.齊次性:I cx I = Ic I I x ICJ3.三角不等式：I x+y I <IxI +I y|X - y1 勻 x-y|則稱Cn中定義了向量范數(shù)丨x丨為向量x的范數(shù)。 可見向量范數(shù)是向量的一種具有特殊性質(zhì)的實值函數(shù)。常用向量范數(shù)有：(令x=( xi,x2，,xn)T )1-范數(shù)：1x I 1= I Xi I + I X2 I + + I Xn I2-范數(shù)：|x I2=( IXiI 2+IX2I

9、2+IXnI2)1/2s -范數(shù)：1x I o0=max( I xi I , I X2 I ,I Xn I )易得Ix I 比<丨 x I 2<I x I 1 < n1/2 I x I 2< n I x I 旳P-范數(shù)：nXp =(遲 Xi P)1P,其中 PJ1嚴).iTx| ：., | x| 是等價的，即m,M>0使m I x I :. <1 x I < M I x I可根據(jù)范數(shù)的連續(xù)性來證明它.由定理1可得定理2.lim x(k)二x”u lim._x(k)x” =0，其中*為向量的任一種范數(shù)。此時 稱x(k)收斂于x,記作x(k) x(kx),

10、或 lim x(k) = x。 矩陣的范數(shù)定義2. 設(shè) ： X Cn nX R,滿足1. 正定性：丨 X | > 0, | X | =0 iff X=02. 齊次性:| cX |= | c| X |，-c C3. 三角不等式：|X+Y | <| X |+| Y |4. 相容性：| XY| <| X | | Y|則稱cn n中定義了矩陣范數(shù),| X |為矩陣X的范數(shù).注意:矩陣X可視為n2維向量，故有前三條性質(zhì).因此定理1,2中向量的等價性 和向量序列收斂的概念與性質(zhì)等也適合于矩陣.第四條，是考慮到矩陣乘法關(guān)系而設(shè)| Ax | <| A | | x |所謂由向量范數(shù)導(dǎo)出的

11、矩陣范數(shù)與該向量范數(shù)就是相容的定理3.設(shè)A是n x n矩陣,|是n維向量范數(shù)則| A | =max | Ax | / | x | =1 = max | Ax | / | x | , x 0是一種矩陣范數(shù)，稱為由該向量范數(shù)導(dǎo)出的矩陣范數(shù)或算子范數(shù)，它們具有相容性或者說是相容的。單位矩陣的算子范數(shù)為1?？梢宰C明任一種矩陣范數(shù)總有與之相容的向量范數(shù).例如定義：| x| = | X | ,X=(XXx)常用的三種向量范數(shù)導(dǎo)出的矩陣范數(shù)是n1-范數(shù):| A | 1= max | Ax | 1 | x | 1=1= max 瓦 ajj1芻勺72-范數(shù):| A | 2=max | Ax | 2/ II x

12、| 2=1=1，刀是 ATA 的最大特征值.ng -范數(shù)：| A | o0=max | Ax | / | x | =1= mgxZ 兇 1 一蟲 jm此外還有Frobenius范數(shù):Afn=(為.aiji,j弓2f2 .它與向量2-范數(shù)相容14 / 11 矩陣譜半徑定義3.設(shè)A是nx n矩陣，入是其特征值，i=1,2，,n.稱：(a)二m丄護為A的譜半徑譜半徑是矩陣的函數(shù)，但非矩陣范數(shù).對任一矩陣范數(shù)有如下關(guān)系：PA)W | A I因為任一特征對入x,Ax= 2x,令X=(xxx),可得AX = X兩邊取范數(shù)，由矩 陣范數(shù)的相容性和齊次性就導(dǎo)出結(jié)果.定理3.矩陣序列I ,A,A2,Ak,收斂于

13、零的充分必要條件是PA)<1.5.5 誤差分析 病態(tài)現(xiàn)象例3給出一個方程組順序消去法解的誤差很大，主元素法解的誤差很小該方程組數(shù)據(jù)有微小變化時解的變化也小但有些方程組不是這樣的，數(shù)據(jù)有微小變化時解的變化大換句話說后一種方程組對數(shù)據(jù)變化敏 感，前一種方程組對數(shù)據(jù)變化不敏感，這兩種方程組(和相應(yīng)的矩陣)分別 稱為病態(tài)的和良態(tài)的.例5. 病態(tài)方程組11J 1.000111.00012爲卜1例6.病態(tài)矩陣-' 11/21/31/41'16-120240-140 1H4 =1/21/31/41/5-1201200-270016801/31/41/51/6240-27006480-4

14、200-1/41/51/61/7 一1 1401680-42002800 一H4取五位有效數(shù)字，其逆誤差在前面第二、三位上_ 16.248-122.72246.49-144.20-122.721229.9-2771.31726.1246.49-2771.36650.1-4310.0.-144.201726.1-4310.02871.1 一 擾動分析與矩陣條件數(shù)現(xiàn)在考慮系數(shù)、右端項有擾動時解的變化，也就是數(shù)據(jù)有誤差時解的誤差. 設(shè)Ax=b,右端項有擾動 A(x+x)=b+b,A可逆解皆存在惟一，其差6 x=A-1 6 b, | 6 x |<| A-1 II § b I ,I 6

15、x| / | x| < ( | A-1 | | A | ) | 6 b| / | b|再考慮系數(shù)有擾動(A+A) (x+x)=b.首先，當A可逆,I A-1 II'： A I <1時 A A 可逆.因為此時 PA-1A) < | A-1A | < | A-1 I IA I <1,I + A-1 A可逆，從而A+A=A(I +A-V A)可逆.原方程與擾動方程解皆 存在惟一,二方程相減有rrrr1 rrA x= - A( x+ x),x= - AA(x+ x)1兩邊取范數(shù)可得Ix |<| A- II A I ( I x I + Ix I )從而有XA_

16、AA x t A|、A| |A類似的方法不難導(dǎo)出一般情況下，即系數(shù)、右端項都有擾動時的估計rM+H麗冋丿注意到估計式表明:I A-1 I I A I不大對解的影響也不大；丨A-1 I I A丨越 大，擾動對解的影響也越大.這就是說該向量是方程組敏感性以及病態(tài)或 良態(tài)的度量，稱為矩陣的 條件數(shù)，記為Cond(A)、= I A-1 I、I AL.它有如 下性質(zhì)：1. Co nd(A) > 12. Co nd(cA)=C on d(A),c 03. Cond(A)2= I A-1 I 21 A I 2=心幾稱為譜條件數(shù)。汕為分別 是AhA的最大和最小特征值故正交矩陣,酉矩陣的譜條件數(shù)為1.在例

17、 1 中有 Cond (A)1=2.0001 x 10 .例 2 中 Co nd (H 4)1=28000.另外, 計算機計算解可歸結(jié)為數(shù)據(jù)有一定擾動的準確解，因而可據(jù)以事先估計計算解的誤差(向后誤差分析). 事后誤差分析計算解的誤差還可根據(jù)下列定理用計算解的剩余量估計* * 定理4.設(shè)x和x分別為非奇異方程組 Ax= b(工0)的準確解和近似解，r為x的剩余 量 r = b- Ax* 貝 U蘭 Co nd(A)目b*11因為 I bI = I Ax I <I A I I xI , I x -x I = I A- r I<I A- II r I . 由此可見對病態(tài)方程組剩余量小時誤

18、差還可能很大 .例7.解方程組0.780x 什0.563x2=0.2170. 913x1+0.659x2=0.254T *tT *T解x=(1,-1) ,x =(0.341,-0.087) ,r=(-0.000 001,0) ,x-x =(-0.659,0.913)二實驗部分本章實驗內(nèi)容:實驗題目：Gauss消元法，追趕法,范數(shù)。實驗內(nèi)容：編制用Gauss消元法求解線性方程組 Ax=f的程序。 編制用追趕法求解線性方程組Ax=f的程序。 編制向量和矩陣的范數(shù)程序。實驗?zāi)康模毫私釭auss消元法原理及實現(xiàn)條件，熟練掌握Gauss消元法解方程組的算法，并能計算行列式的值。 掌握追趕法，能利用追趕法

19、求解線性方程組。 理解向量和矩陣范數(shù)定義，性質(zhì)并掌握其計算方 法.編程要求：利用Gauss消元法，追趕法解線性方程組。分析誤差。 計算算法:Gauss消元法：1. 消元過程設(shè) a0，對 k =1,2, , n1 計算rniik =aikk)/akk°pF =ajk) mikakjk)bT = b(k)5kbkk)i, j = k 1, k 2, n2. 回代過程' _(b(i) 舟 a (i) /a(i)xi 一( d l a j x j)/ aiij=t+i 二 n 1,2,1追趕法：1. 分解 Ax=f：= c1 /b1, ： j = g /(bi - ai i )(i

20、=2,3 廠，n-1)2. 解 Lg=f,求 g:g1 二什/匕® = (fi -aigt)/(bi -ai :( i =2,3, n)3解 Ux=g,求 x: Xn =yn,Xi- ： iX(i 二 n-1,n-2,2,1 )范數(shù)：常用向量范數(shù)有：(令X=( X1,x2，,Xn)T )1- 范數(shù):I x I 1= I X1 I + I X2 | + + | Xn |2- 范數(shù)：I x I 2=( I X1 | 2+| X2| 2+| Xn | 2) 1/2x -范數(shù)：| x I ： ：=max( | X1 | , | X2 | ,| Xn | ) 常用的三種向量范數(shù)導(dǎo)出的矩陣范數(shù)是

21、：1-范數(shù):I A | 1= max | Ax | 1/ I x I 1=1= max遲 a”1總應(yīng)y2-范數(shù):| A I 2=max I Ax | 2 I x | 2=1= r，入是 AA 的最大特征值.naij00 -范數(shù):I A | : =max | Ax | ：I x I .-=1= maxv 實驗例題：用Gauss消元法解方程組實驗例題：實驗例題:"2-11、-1-23X2=5J31丿N丿f丿用追趕法解二對角方程組2-1000 _-12-100A =0-12-1000-12-1000-12 一A'0.60.5"Q10.3,/b =一們0 Io |o 0JA

22、x=b,其中計算A的行列范數(shù).程序:Gauss消元法function x=Gauss(A,b)%A是線性方程組的系數(shù)矩陣,b為自由項.n=length(A)for k=1:n-1m(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1: n)-m(k+1:n,k)*A(k,k+1:n); b(k+1:n)=b(k+1: n)-m(k+1:n,k)*b(k);endx=zeros(n,1);x(n)=b( n)/A(n,n);for k=n-1:-1:1x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)/A(k,k);endx=x&

23、#39;disp(sprintf('k x(k)');for i=0:ndisp(sprintf('%d %f ',i,x(i+1); end 數(shù)值結(jié)果： x=Gauss(A,b)n =3 kx(k)01.11111110.77777822.555556程序 : 追趕法 function x,y,beta=zhuiganfa(a,b,c,f) %a,b,c 是三對角陣的對角線上的元素 ,f 是自由項 . n=length(b); beta(1)=c(1)/b(1);for i=2:nbeta(i)=c(i)/(b(i)-a(i)*beta(i-1); end y

24、(1)=f(1)/b(1);for i=2:ny(i)=(f(i)-a(i)*y(i-1)/(b(i)-a(i)*beta(i-1); end x(n)=y(n);for i=n-1:-1:1x(i)=y(i)-beta(i)*x(i+1);enddisp(sprintf('k x(k) y(k) beta(k)'); for i=0:ndisp(sprintf('%d %f ',i,x(i+1),y(i+1),beta(i+1); end 數(shù)值結(jié)果 : a=0 -1 -1 -1 -1'b=2 2 2 2 2' c=-1 -1 -1 -1 0&#

25、39; f=1 0 0 0 0'x,y,beta=zhuiganfa(a,b,c,f)kx(k) y(k)beta(k)00.833333 5.000000e-001-0.50000010.666667 3.333333e-001-0.66666720.500000 2.500000e-001-0.75000030.333333 2.000000e-001-0.80000040.166667 1.666667e-0010.000000程序：1. 列范數(shù) : function fan=lie(A) %A 為已知矩陣n=length(A) for j=1:n x(j)=0 for i=1:n x(j)=x(j)+abs(A(i,j);endendx(n)');%f',i,x(i+