解線性方程組的基本思想

1、四：基本方法基本思路將在解題的過程中得到體現(xiàn)。1（求線性方程組的唯一解或特解），這類問題的求法分為兩類：一類主要用于解低階稠密矩陣 直接法；一類是解大型稀疏矩陣 迭代法。1.1利用矩陣除法求線性方程組的特解（或一個解）方程：AX=b，解法：X=Ab，(注意此處不是/)例1-1    求方程組 的解。解:    A = ; = ；b=(1,0,0,0,1)由于>>rank(A)=5,rank( )=5    %求秩，此為R（A）=R（ ）>=n的情形，有唯一解。  &#

2、160;   >>X= Ab      %求解 X =(2.2662, -1.7218, 1.0571,-0.5940, 0.3188) 或用函數(shù)rref求解，>>sv=rref(A:b);所得sv的最后一列即為所要求的解。12 利用矩陣的LU、QR和cholesky分解求方程組的解 這三種分解，在求解大型方程組時很有用。其優(yōu)點是運算速度快、可以節(jié)省磁盤空間、節(jié)省內(nèi)存。I) LU分解又稱Gauss消去分解，可把任意方陣分解為下三角矩陣的基本變換形式（行交換）和上三角矩陣的乘積。即A=LU，L

3、為下三角陣，U為上三角陣。則：A*X=b          變成L*U*X=b所以X=U(Lb)     這樣可以大大提高運算速度。命令    L，U=lu (A)在matlab中可以編如下通用m 文件：在Matlab中建立M文件如下% exp1.mA;b;L，U=lu (A);X=U(Lb)II）Cholesky分解若A為對稱正定矩陣，則Cholesky分解可將矩陣A分解成上三角矩陣和其轉(zhuǎn)置的乘積，即：  

4、0;  其中R為上三角陣。方程    A*X=b    變成     所以   在Matlab中建立M文件如下% exp2.mA;b;R，R=chol(A);X=R(Rb)III）QR分解對于任何長方矩陣A，都可以進行QR分解，其中Q為正交矩陣，R為上三角矩陣的初等變換形式，即：A=QR方程    A*X=b    變形成     QRX=b所以 

5、    X=R(Qb)上例中    Q, R=qr(A)X=R(QB)在Matlab中建立M文件如下% exp3.mA;b;Q，R=qr(A);X=R(Qb)2求線性齊次方程組的通解(A*X=0)在Matlab中，函數(shù)null用來求解零空間，即滿足A•X=0的解空間，實際上是求出解空間的一組基（基礎(chǔ)解系）。在Matlab中建立M文件如下% exp4.mformat rat      %指定有理式格式輸出A;b=0;r=rank(A);bs=null(A，r); %一

6、組基含(n-r)個列向量% k ,k ,k % X= k *bs(:,1)+ k *bs(:,2)+ k *bs(:,n-r) 方程組的通解pretty(X)       %讓通解表達式更加精美3    求非齊次線性方程組的通解（A*X=b）非齊次線性方程組需要先判斷方程組是否有解，若有解，再去求通解。因此，步驟為：第一步：判斷AX=b是否有解，(利用基本思路的第一條)若有解則進行第二步第二步：求AX=b的一個特解第三步：求AX=0的通解第四步：AX=b的通解為： AX=0的通解加上AX=b的一個特解。在

7、Matlab中建立M文件如下% exp4.mclear allA;b;                    %輸入矩陣A,bm,n=size(A);R=rank(A);B=A b;Rr=rank(B);format rat if R=Rr&R=n            % n為未知數(shù)的個數(shù)

8、，判斷是否有唯一解x=Ab;elseif R=Rr&R<n        %判斷是否有無窮解x=Ab                   %求特解C=null(A, r )           %求AX=0的基礎(chǔ)解系,所得C為n-R列矩

9、陣，這n-R列即為對%應(yīng)的基礎(chǔ)解系                        % 這種情形方程組通解xx=k(p)*C(:,P)(p=1n-R)else X= No solution!    % 判斷是否無解end第3章 線性方程組的迭代解法3.1實驗?zāi)康?理解線性方程組計算機解法中的迭代解法的求解過程和特點，學(xué)習(xí)科學(xué)計算的方法和簡單

10、的編程技術(shù)。 3.2  概念與結(jié)論1.     n階線性方程組 如果未知量的個數(shù)為 n ,而且關(guān)于這些未知量x1,x2, ,xn 的冪次都是一次的(線性的)那末, n 個方程 a11x1+a12x2+ +a1nxn=b1 (1) an1x1+an2x2+ +annxn=bn構(gòu)成一個含n個未知量的線性方程組,稱為n階線性方程組。其中,系數(shù)a11,a1n,a21, ,a2n, ,an1, ,ann 和b1, ,bn 都是給定的常數(shù)。 方程組(1)也常用矩陣的形式表示,寫為 Ax=b其中,A是由系數(shù)按次序排列構(gòu)成的一個n階矩陣, 稱為方程組的系數(shù)矩陣

11、,x和b都是n維向量,b稱為方程組的右端向量。2. n階線性方程組的解 使方程組(1)中每一個方程都成立的一組數(shù)x1*,x2*, ,xn* 稱為式(1)的解,把它記為向量的形式,稱為解向量.3. 向量范數(shù)的三種常用范數(shù)4矩陣的四種常用范數(shù)5.譜半徑 設(shè) n´n 階矩陣A的特征值為l i(i=1,2,3n),則稱 r (A) = MAX | li| 1£i£ n為矩陣A的譜半徑. 矩陣范數(shù)與譜半徑之間的關(guān)系為: r (A) £ |A|6.嚴(yán)格(行)對角占優(yōu)陣A 如果 矩陣 A=(aij)滿足 n |aii| > S |aij| i=1,2,n, j=

12、1,j¹i 則稱方陣A是嚴(yán)格(行)對角占優(yōu)的.7.收斂定理 對任意初始向量x(0)及任意右端向量 g，由迭代x(k+1) =B x(k) +g產(chǎn)生的迭代向量序列x(k)收斂的充要條件是譜半徑r（B）<18.收斂判別條件 判別條件1： 若|B|<1, 則迭代x(k+1) =B x(k) +g 對任何初始向量x(0)都收斂. 判別條件2： 如果A為嚴(yán)格對角占優(yōu)陣，則其 Jacobi迭代和Seidel迭代對任何初始向量x(0)都收斂。 判別條件3： 如果A為對稱正定陣，則其 Seidel迭代對任何初始向量x(0)都收斂。9.迭代法的誤差估計 若|B|<1，則對迭代格式 x

13、(k+1) =B x(k) +g 有3.3 程序中Mathematica語句解釋 a*matrix 數(shù)a與矩陣matrix相乘matrix1+matrix2 矩陣matrix1和矩陣matrix2相加（注意矩陣的大小相同）matrix1.matrix2 矩陣matrix1和矩陣matrix2相乘(注意矩陣乘法的規(guī)則) Transposematrix 求矩陣matrix轉(zhuǎn)置Inversematrix 求矩陣(方陣) matrix 的逆 DiagonalMatrixlist 使用列表list中的元素生成一個對角矩陣.IdentityMatrixn 生成n階單位矩陣Maxx 求向量x中元素的最大值3

14、.4 方法、程序、實驗 解線性方程組的迭代法是將線性方程組 Ax=b 化為等價線性方程組 x=Bx+f 再由矩陣迭代格式 x(k+1)=Bx (k)+f構(gòu)造向量序列x(k)來求線性方程組解的。如果得出的向量序列x(k)收斂至某個向量x*，則可得該向量x*就是所求方程組 Ax=b 的準(zhǔn)確解.線性方程組的迭代法主要有Jocobi迭代法、Seidel迭代法和超松弛(Sor)迭代法。1. Jocobi迭代法1) Jocobi迭代法的構(gòu)造過程 假設(shè)aii¹0，依次在第i個方程解出x i , i=1,2,¼,n并令 cij = -aij /aii (i¹j) ， gi= bi

15、 /aii 就得到如下Jocobi迭代格式： x1(k+1)= c12x2(k)+c13x3(k)+×××× +c1nxn(k)+g1 x2(k+1)=c21x1(k) +c23x3(k)+×××× +c2nxn(k)+g2 。 xn(k+1)=cn1x1(k) +cn2x2(k)+×××× +cn(n-1)xn-1(k) + gn 若令則有Jocobi迭代的矩陣格式：x(k+1) = BJx(k) +gJBJ 稱為Jocobi迭代矩陣。 Jocobi迭代可以寫成如下緊湊格式

16、：在給定初始迭代向量x(0)后就可以進行Jocobi迭代求解了。2) Jacobi迭代算法1.輸入變量個數(shù)n、初值向量x（0）、迭代精度eps、系數(shù)矩陣A、常數(shù)項b 和迭代最大次數(shù)nmax2 For i=1,2,n2.1 如果|aii|<eps1,則輸出“迭代失敗”提示并終止3. Bj Ü E-D-1A4. gj Ü D-1b5.For k=1,2,nmax5.1 x ÜBj.x0+ gj 5.2 如果|x-x0|<eps ,輸出解向量x ，終止；否則x（0） Ü x6. 如果|x-x0|>eps ,輸出迭代失敗，終止。3) Jacob

17、i 迭代法程序Cleara,b,x;nmax=500;n=Input“線性方程組階數(shù)n=”;a=Input"系數(shù)矩陣A="；b=Input"常數(shù)項b="；x0=Input"輸入迭代初值向量x0"；eps1=0.000001;eps=Input"輸入精度控制eps="；DoIfAbsai,i<eps1,t1=1;Return,t1=0,i,1,n;Ift1=1, Print"Jacobi迭代法失效", d=DiagonalMatrixTableai,i,i,1,n; d1=Inversed;

18、 bj=IdentityMatrixn-d1.a; gj=d1.b; Do x=bj.x0+gj; err=MaxAbsx-x0; Print"x=",x/N," i=",i," err=",err/N; IfNerr<eps,Break,x0=x, i,1,nmax; Iferr>=eps,Print"迭代失敗 " 說明 本程序用于求線性方程組Ax=b的解。程序執(zhí)行后，先通過鍵盤輸入線性方程組階數(shù)n、系數(shù)矩陣A、常數(shù)項b、迭代初值向量x0和輸入精度控制eps，程序即可給出每次迭代的次數(shù)和對應(yīng)的迭代向量

19、序列x(k),其中最后輸出的結(jié)果即為所求的根。如果迭代超出500次還沒有求出滿足精度的根則輸出迭代失敗提示，如果出現(xiàn)主對角線元素aii=0給出Jacobi迭代法失效提示。程序中變量說明x0：存放初始向量和迭代過程中的向量x（k）x: 存放迭代過程中的向量x（k+1）nmax:存放迭代允許的最大次數(shù)err:存放誤差|x-x0|µt1:臨時變量注：迭代最大次數(shù)可以修改為其他數(shù)字。4)例題與實驗例1.用Jacobi 迭代法解如下線性方程組 5x1+2x2+x3= -12 -x1+4x2+2x3= 20 2x1-3x2+10x3= 3要求誤差|x(k+1)-x(k)|µ<10

20、-4，并用取不同初值的方法實驗觀察迭代收斂的情況。解：執(zhí)行Jacobi迭代法程序后在輸入的四個窗口中按提示分別輸入:3、5, 2, 1, -1, 4, 2, 2, -3, 10、-12, 20, 3、0, 0, 0、0.0001每次輸入后用鼠標(biāo)點擊窗口的“OK”按扭，得如下輸出結(jié)果： x=-2.4, 5., 0.3 i=1 err=5.x=-4.46, 4.25, 2.28 i=2 err=2.06x=-4.556, 2.745, 2.467 i=3 err=1.505x=-3.9914, 2.6275, 2.0347 i=4 err=0.5646x=-3.85794, 2.9848, 1.8

21、8653 i=5 err=0.3573x=-3.97123, 3.09225, 1.96703 i=6 err=0.113286x=-4.03031, 3.02368, 2.02192 i=7 err=0.0685705x=-4.01386, 2.98146, 2.01316 i=8 err=0.042216x=-3.99522, 2.98995, 1.99721 i=9 err=0.0186374x=-3.99542, 3.00259, 1.99603 i=10 err=0.0126367x=-4.00024, 3.00313, 1.99986 i=11 err=0.0048186x=-4.

22、00122, 3.00001, 2.00099 i=12 err=0.00312067x=-4.0002, 2.9992, 2.00025 i=13 err=0.00102319x=-3.99973, 2.99983, 1.9998 i=14 err=0.000625697x=-3.99989, 3.00017, 1.99989 i=15 err=0.00034136x=-4.00005, 3.00008, 2.00003 i=16 err=0.000155236x=-4.00004, 2.99997, 2.00003 i=17 err=0.000106099x=-4., 2.99997, 2

23、. i=18 err=0.0000414468此結(jié)果說明迭代18次，求得誤差為err=0.0000414468的近似解，最后顯示的近似解向量為x=-4., 2.99997, 2.，它表示所求解為 x1=-4，x2=2.99997，x3=2 。本題的準(zhǔn)確解為 x1=-4，x2=3，x3=2 。如果將如上輸入的初值改為21, -18, 30，執(zhí)行Jacobi迭代法程序后得如下輸出結(jié)果： x=-1.2, -4.75, -9.3 i=1 err=39.3x=1.36, 9.35, -0.885 i=2 err=14.1x=-5.963, 5.7825, 2.833 i=3 err=7.323x=-5.

24、2796, 2.09275, 3.22735 i=4 err=3.68975x=-3.88257, 2.06642, 1.98375 i=5 err=1.39703x=-3.62332, 3.03749, 1.69644 i=6 err=0.97106x=-3.95428, 3.24595, 1.93591 i=7 err=0.330963x=-4.08556, 3.04347, 2.06464 i=8 err=0.202475x=-4.03032, 2.94629, 2.03015 i=9 err=0.0971858x=-3.98455, 2.97734, 1.98995 i=10 err=

25、0.0457716x=-3.98893, 3.00889, 1.99011 i=11 err=0.0315451x=-4.00158, 3.00771, 2.00045 i=12 err=0.0126504x=-4.00318, 2.99938, 2.00263 i=13 err=0.00833246x=-4.00028, 2.99789, 2.00045 i=14 err=0.00289753x=-3.99925, 2.99971, 1.99942 i=15 err=0.0018145x=-3.99977, 3.00048, 1.99976 i=16 err=0.000770763x=-4.

26、00014, 3.00018, 2.0001 i=17 err=0.000375926x=-4.00009, 2.99992, 2.00008 i=18 err=0.000261658x=-3.99998, 2.99994, 1.99999 i=19 err=0.000107579x=-3.99997, 3.00001, 1.99998 i=20 err=0.0000714045從計算結(jié)果可以看到雖然所取的初值不同，但迭代總是收斂相同的結(jié)果?？疾毂绢}系數(shù)矩陣可以發(fā)現(xiàn)它是嚴(yán)格對角占優(yōu)矩陣，由收斂判別法，Jacobi迭代對任意初值都是收斂的，我們通過實際計算驗證了這個結(jié)果。例2.通過Jacobi

27、迭代序列觀察用Jacobi 迭代法解如下線性方程組 x1+2x2+2x3= -12 -x1+4x2+x3= 20 2x1-3x2+x3= 3的收斂性。解：執(zhí)行Jacobi迭代法程序后在輸入的四個窗口中按提示分別輸入:3、1, 2, 2, -1, 4, 1, 2, -3, 1、-12, 20, 3、0, 0, 0、0.0001每次輸入后用鼠標(biāo)點擊窗口的“OK”按扭，得如下輸出結(jié)果： x=-12., 5., 3. i=1 err=12.x=-28., 1.25, 42. i=2 err=39.x=-98.5, -12.5, 62.75 i=3 err=70.5x=-112.5, -35.3125,

28、 162.5 i=4 err=99.75.246.719, -120.176, 696.719 i=8 err=763.5x=-1165.09, -107.5, -850.965 i=9 err=1547.68x=1904.93, -73.5303, 2010.67 i=10 err=3070.02x=-3886.28, -21.4355, -4027.45 i=11 err=6038.12x=8085.77, 40.2917, 7711.26 i=12 err=11972.1x=-15515.1, 98.6279, -16047.7 i=13 err=23758.9x=31886.1, 13

29、8.141, 31329.1 i=14 err=47401.2x=-62946.5, 144.247, -63354.7 i=15 err=94832.5x=126409., 107.068, 126329. i=16 err=189683.x=-252883., 25.0775, -252494. i=17 err=379292.x=504925., -92.4305, 505845. i=18 err=758339.通過觀察迭代過程中的誤差是不斷變大的特點可以知道本題的Jacobi 迭代序列是不收斂的，因此，本題線性方程組不能用Jacobi 迭代法求解。這個實驗說明并不是每個線性方程組都能

30、用Jacobi 迭代法求解。2.     Seidel迭代1) Seidel迭代的構(gòu)造過程 為了加快收斂速度，同時節(jié)省計算機的內(nèi)存，對Jocobi迭代作如下的改進：每算出一個分量的近似值，立即用到下一個分量的計算中去，即用迭代格式： x1(k+1)= c12x2(k)+c13x3(k)+×××× +c1nxn(k)+g1 x2(k+1)=c21x1(k+1) +c23x3(k)+×××× +c2nxn(k)+g2 。 xn(k+1)=cn1x1(k+1) +cn2x2(k+1

31、)+××××+cn(n-1)xn-1(k+1) + gn如上迭代可以寫成如下緊湊格式：這樣所得的迭代法就稱為Seidel迭代法。 利用Ax=b 及A=L+D+U,其中D為對角矩陣,L,U分別為嚴(yán)格下，上三角矩陣則有Seidel迭代法的矩陣形式為 x(k+1)= Bs x(k)+gs 其中：Bs =-(L+D)-1U，g s=D-1bSeidel迭代矩陣為 Bs =-(L+D)-1U。 在給定初始迭代向量x(0)后就可以進行Seidel迭代求解了。 Jacobi迭代和Seidel迭代格式可表述為統(tǒng)一形式:2) Seidel迭代算法1.輸入變量個數(shù)n、初值向

32、量x（0）、迭代精度eps、系數(shù)矩陣A、常數(shù)項b 和迭代最大次數(shù)nmax2. For i=1,2,n2.1 如果|aii|<eps1,則輸出“迭代失敗”提示并終止3.For i=1,2,nmax 3.1 For i=1,2,n 3.2 如果|x(k+1)-x(k)|µ<eps ,輸出解向量x ，終止；否則x（k） Ü x(k+1)4. 如果|x-x0|µ<eps ,輸出迭代失敗，終止。3) Seidel 迭代法程序Cleara,b,x;nmax=500;n=Input“線性方程組階數(shù)n=”;a=Input"系數(shù)矩陣A="；b=

33、Input"常數(shù)項b="；x0=Input"輸入迭代初值向量x0"；eps1=0.000001;eps=Input"輸入精度控制eps="； x=x0; DoIfAbsai,i<eps1,t1=1;Return,t1=0,i,1,n;Ift1=1, Print"Seidel迭代法失效", Do Dou1=Sumai,j*xj,j,1,i-1; u2=Sumai,j*x0j,j,i+1,n; xi=(bi-u1-u2)/ai,i, i,1,n; err=MaxAbsx-x0; Print"x=&quo

34、t;,x/N," k=",k," err=",err/N; IfNerr<eps,Break,x0=x, k,1,50; Iferr>=eps,Print"迭代失敗 " 說明 本程序用于求線性方程組Ax=b的解。程序執(zhí)行后，先通過鍵盤輸入線性方程組階數(shù)n、系數(shù)矩陣A、常數(shù)項b、迭代初值向量x0和輸入精度控制eps，程序即可給出每次迭代的次數(shù)和對應(yīng)的迭代向量序列x(K),其中最后輸出的結(jié)果即為所求的根。如果迭代超出500次還沒有求出滿足精度的根則輸出迭代失敗提示，如果出現(xiàn)主對角線元素aii=0給出Jacobi迭代法失效提示。

35、程序中變量說明x0：存放初始向量和迭代過程中的向量x（k）x: 存放迭代過程中的向量x（k+1）nmax:存放迭代允許的最大次數(shù)err:存放誤差|x-x0|µt1,u1,u2:臨時變量注：迭代最大次數(shù)可以修改為其他數(shù)字4)例題與實驗例3.用Seidel 迭代法解如下線性方程組 5x1+2x2+x3= -12 -x1+4x2+2x3= 20 2x1-3x2+10x3= 3要求誤差|x(k+1)-x(k)|µ<10-4。解：執(zhí)行Seidel迭代法程序后在輸入的四個窗口中按提示分別輸入:3、5, 2, 1, -1, 4, 2, 2, -3, 10、-12, 20, 3、0,

36、 0, 0、0.0001每次輸入后用鼠標(biāo)點擊窗口的“OK”按扭，得如下輸出結(jié)果：x=-2.4, 4.4, 2.1 k=1 err=4.4x=-4.58, 2.805, 2.0575 k=2 err=2.18x=-3.9335, 2.98787, 1.98306 k=3 err=0.6465x=-3.99176, 3.01053, 2.00151 k=4 err=0.0582625x=-4.00451, 2.99812, 2.00034 k=5 err=0.0127509x=-3.99931, 3., 1.99986 k=6 err=0.00519946x=-3.99997, 3.00007,

37、2.00002 k=7 err=0.000659841x=-4.00003, 2.99998, 2. k=8 err=0.0000916628 此結(jié)果說明迭代8次，求得誤差為err=0.0000916628的近似解，最后顯示的近似解向量為x=-4.00003, 2.99998, 2.，與本題的準(zhǔn)確解為 x1=-4.，x2=3.，x3=2. 誤差很小。注意到本題在同樣迭代誤差和初值前提下，用Jacobi 迭代需要18次迭代才獲得具有同樣精度的解，這說明，在都收斂的前提下，Seidel 迭代比Jacobi 迭代收斂快。例4.設(shè)線性方程組為 x1+2x2-2x3= 1 x1+x2+x3= 1 2x1+2x2+x3= 1考察用Jacobi迭代法和Seidel 迭代法求解該線性方程組的收斂情況。如果收斂，給出誤差滿足|x(k+1)-x(k)|µ<10-4的解。解：執(zhí)行Seidel迭代法程序后在輸入的四個窗口中按提示分別輸入:3、1, 2, -2, 1,1,1, 2,2,1、1,1,1、0, 0, 0、0.0001每次輸入后用鼠標(biāo)點擊窗口的“OK”按扭，得如下輸出部分結(jié)果：.x=-5123., 5379., -511. k=9 err=3072