-解線性方程組的直接解法

1、科學(xué)計(jì)算與科學(xué)計(jì)算與MATLABMATLAB第九講解線性方程組的直接解法第九講解線性方程組的直接解法內(nèi)容提要內(nèi)容提要n引言引言nGauss消元法消元法n列主元素消元法列主元素消元法n誤差分析誤差分析nMATLAB的線性方程組求解函數(shù)的線性方程組求解函數(shù)1n小結(jié)小結(jié)1、引言、引言例：小行星軌道問題：例：小行星軌道問題：天文學(xué)家要確定一小行星的軌道天文學(xué)家要確定一小行星的軌道,在軌道平面建立以太陽為在軌道平面建立以太陽為原點(diǎn)的直角坐標(biāo)系原點(diǎn)的直角坐標(biāo)系.在坐標(biāo)軸上取天文測量單位在坐標(biāo)軸上取天文測量單位(一天文單位為地一天文單位為地球到太陽的平均距離：球到太陽的平均距離：9300萬哩萬哩),對小行星

2、作對小行星作5次觀察次觀察, 測得軌道測得軌道上上5個(gè)點(diǎn)的坐標(biāo)數(shù)據(jù)如下：個(gè)點(diǎn)的坐標(biāo)數(shù)據(jù)如下： x 5.7640 6.2860 6.7590 7.1680 7.4800 y 0.6480 1.2020 1.8230 2.5260 3.3600 橢圓的一般方程橢圓的一般方程: a1x2 + a2xy + a3y2 + a4x + a5y + 1 = 0 將數(shù)據(jù)逐個(gè)代入將數(shù)據(jù)逐個(gè)代入,可得五個(gè)方程的方程組，求解該線性方程組即可得五個(gè)方程的方程組，求解該線性方程組即可得行星軌道方程?？傻眯行擒壍婪匠?。對一般線性方程組對一般線性方程組: A x = b, 其中其中nnnnnnaaaaaaaaaA2122

3、22111211nbbbb21nxxxX21當(dāng)且僅當(dāng)矩陣當(dāng)且僅當(dāng)矩陣A行列式不為行列式不為0時(shí)，即時(shí)，即A非奇異時(shí)，方程組非奇異時(shí)，方程組存在唯一解，可根據(jù)克萊姆法則求解，其算法設(shè)計(jì)如下：存在唯一解，可根據(jù)克萊姆法則求解，其算法設(shè)計(jì)如下：(1) 輸入系數(shù)矩陣輸入系數(shù)矩陣A和右端向量和右端向量b；(2)計(jì)算系數(shù)矩陣計(jì)算系數(shù)矩陣A的行列式值的行列式值D，如果，如果D=0，則輸出錯(cuò)誤信，則輸出錯(cuò)誤信息，結(jié)束，否則進(jìn)行第息，結(jié)束，否則進(jìn)行第(3)步；步；(3) 對對k=1,2,n，用用b替換替換A的第的第k列數(shù)據(jù)，并計(jì)算替換后矩列數(shù)據(jù)，并計(jì)算替換后矩陣的行列式值陣的行列式值Dk；(4) 計(jì)算并輸出計(jì)算

4、并輸出x1 = D1 / D，x2 = D2 / D，,xn=Dn/D， 結(jié)束。結(jié)束。 克萊姆法則只適用于低階方程組，高階方程組工作量太克萊姆法則只適用于低階方程組，高階方程組工作量太大，故一般用數(shù)值方法求解。數(shù)值方法分兩類：大，故一般用數(shù)值方法求解。數(shù)值方法分兩類： 直接解法：在計(jì)算過程中，如果所有運(yùn)算都是精確的，直接解法：在計(jì)算過程中，如果所有運(yùn)算都是精確的，在理論上，經(jīng)過有限次運(yùn)算就可以得到精確解，適用于在理論上，經(jīng)過有限次運(yùn)算就可以得到精確解，適用于變量較少的方程組。變量較少的方程組。 迭代解法：近似解法，運(yùn)算次數(shù)因要求的計(jì)算精度而變迭代解法：近似解法，運(yùn)算次數(shù)因要求的計(jì)算精度而變化。

5、化。迭代法也稱輾轉(zhuǎn)法，是一種不斷用變量的舊值遞推新值的迭代法也稱輾轉(zhuǎn)法，是一種不斷用變量的舊值遞推新值的過程。過程。迭代算法是用計(jì)算機(jī)解決問題的一種基本方法。它利用計(jì)算迭代算法是用計(jì)算機(jī)解決問題的一種基本方法。它利用計(jì)算機(jī)運(yùn)算速度快、適合做重復(fù)性操作的特點(diǎn)，讓計(jì)算機(jī)對一組機(jī)運(yùn)算速度快、適合做重復(fù)性操作的特點(diǎn)，讓計(jì)算機(jī)對一組指令（或一定步驟）進(jìn)行重復(fù)執(zhí)行，在每次執(zhí)行這組指令指令（或一定步驟）進(jìn)行重復(fù)執(zhí)行，在每次執(zhí)行這組指令（或這些步驟）時(shí)，都從變量的原值推出它的一個(gè)新值。（或這些步驟）時(shí)，都從變量的原值推出它的一個(gè)新值。迭代方法的使用：迭代方法的使用：確定迭代變量。確定迭代變量。建立迭代關(guān)系式。

6、建立迭代關(guān)系式。給定初值給定初值對迭代過程進(jìn)行控制。對迭代過程進(jìn)行控制。對迭代結(jié)果判定對迭代結(jié)果判定2、Gauss消元法消元法n2.1 基本思想：基本思想：逐步消去未知元，將方程組化為與逐步消去未知元，將方程組化為與其等價(jià)的上三角方程組求解。其等價(jià)的上三角方程組求解。n2.2 分兩步：分兩步： 第一步第一步: 消元過程消元過程，將方程組消元化為等價(jià)的上三角形，將方程組消元化為等價(jià)的上三角形方程組方程組; 第二步第二步: 回代過程回代過程， 解上三角形方程組解上三角形方程組,得原方程組的得原方程組的解。解。 )1()1()1(2)1(22)1(2211212111nnnnnnnnnnbxabxa

7、xabxaxaxa2.3 Gauss消元的目的消元的目的a11x1+ a12x2+ a1nxn = b1a21x1+ a22x2+ a2nxn = b2an1x1+ an2x2+ annxn = bn原始方程組原始方程組約化方程組約化方程組2.4.1 消元過程消元過程(化一般方程組為上三角方程組化一般方程組為上三角方程組)以四階為例：以四階為例：其系數(shù)增廣矩陣為：其系數(shù)增廣矩陣為：4444343242141343433323213124243232221211414313212111bxaxaxaxabxaxaxaxabxaxaxaxabxaxaxaxa44443424133433323122

8、4232221114131211baaaabaaaabaaaabaaaaA第一輪消元：第一輪消元：計(jì)算計(jì)算3個(gè)數(shù)個(gè)數(shù): m21 m31 m41T = a21 a31 a41T / a11 用用-m21乘矩陣第一行后加到矩陣第二行乘矩陣第一行后加到矩陣第二行; 用用-m31乘矩陣第一行后加到矩陣第三行乘矩陣第一行后加到矩陣第三行; 用用-m41乘矩陣第一行后加到矩陣第四行乘矩陣第一行后加到矩陣第四行;其系數(shù)增廣矩陣變?yōu)椋浩湎禂?shù)增廣矩陣變?yōu)椋?()()()()()()()()()()()()(1414414314213134133132121241231221141312111baaabaaaba

9、aabaaaaA)0(11a第二輪消元：第二輪消元：計(jì)算計(jì)算2個(gè)數(shù)：個(gè)數(shù)：m32 m42T = a32(1) a42(1)T / a22(1) 用用-m32乘矩陣第二行后加到矩陣第三行乘矩陣第二行后加到矩陣第三行; 用用-m42乘矩陣第二行后加到矩陣第四行乘矩陣第二行后加到矩陣第四行;其系數(shù)增廣矩陣變?yōu)椋浩湎禂?shù)增廣矩陣變?yōu)椋?()()()()()()()()()()(2424424323234233121241231221141312112baabaabaaabaaaaA)0()1(22a第三輪消元：第三輪消元：計(jì)算計(jì)算: m43=a43(2)/a33(2)用用-m43乘矩陣第三行后加到矩陣第

10、四行乘矩陣第三行后加到矩陣第四行;其系數(shù)增廣矩陣變?yōu)椋浩湎禂?shù)增廣矩陣變?yōu)椋?()()()()()()()()()(3434423234233121241231221141312113babaabaaabaaaaA)0()2(33a其對應(yīng)的上三角方程組為其對應(yīng)的上三角方程組為)()()()()()()()()(3443442342343233124124312321221414313212111bxabxaxabxaxaxabxaxaxaxa 若對于一般的線性方程組若對于一般的線性方程組Ax=b，其消元過程的計(jì)算，其消元過程的計(jì)算公式為公式為： （k=1,2,n-1）), 1(), 1,(),

11、1()1()1()()1()1()()1(nkibmbbnkjiamaankiaamkkikkikikkjikkijkijkkkikik）（記ijijaa)0(2.4.2 回代過程回代過程(解上三角方程組解上三角方程組)上三角方程組的一般形式為：上三角方程組的一般形式為： 其中其中a11ann0nnnnnnnnbxabxaxabxaxaxa2222211212111) 1 , 2 , 1()(1niaxabxabxiinijjijiinnnn回代過程的計(jì)算公式：回代過程的計(jì)算公式：2.5 工作量計(jì)算：工作量計(jì)算： 消去過程：消去過程：n“”：第：第k步，步，n-k次，共次，共(n-1)+(n-

12、2)+1=n(n-1)/2n“”：第：第k步，步，(n-k)(n-k+1)次，共次，共 (n-1)n+(n-2)(n-1)+12= (n3-n)/3n總工作量：總工作量：s1=n(n-1)/2+ (n3-n)/3 回代過程：回代過程：n“”：nn“”：1+2+(n-1)=n(n-1)/2n總工作量：總工作量：s2=n+ n(n-1)/2= n(n+1)/2故故Gauss消元法的總工作量為：消元法的總工作量為：s=s1+s2= n2+(n3-n)/3克萊姆法則求解的工作量為：克萊姆法則求解的工作量為：“”：（n+1個(gè)個(gè)n階行列式的值）階行列式的值）(n+1)(n-1)n!“”：n故總工作量為：故

13、總工作量為： (n+1)(n-1) n!+n 當(dāng)當(dāng)n=6時(shí)時(shí), Gauss消元法工作量為消元法工作量為106 ;而克萊姆而克萊姆法則求解工作量為法則求解工作量為25206。function U,x=Gauss(A,b)% 順序Gauss消去法求解線性方程組% 輸入?yún)?shù)：% -A：線性方程組的系數(shù)矩陣% -b：線性方程組的右端項(xiàng)% 輸出參數(shù)：% -U：消元后的上三角方程組的增廣矩陣% -x：線性方程組的解n=length(b);% 消元過程for k=1:n-1 m=A(k+1:n,k)/A(k,k); A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n); b(

14、k+1:n)=b(k+1:n)-m*b(k); A(k+1:n,k)=zeros(n-k,1);endU=A,b;% 回代過程x=zeros(n,1);x(n)=b(n)/A(n,n); % 求x_nfor k=n-1:-1:1 x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)/A(k,k); % 求x_k，k=n-1,n-2,1end定理：定理： 約化的主元素約化的主元素ak+1,k+1(k) 0 (k=0,1,n-1)的充分必要條件是的充分必要條件是 矩陣矩陣A的各階順序主子式的各階順序主子式不為零。即不為零。即01111 kkkkkaaaaD注：注：對角線上的元素對角線上的元

15、素ak+1,k+1(k)在在Gauss消去法中作消去法中作用突出，稱約化的主元素。用突出，稱約化的主元素。推論推論: 如果如果A的順序主子式的順序主子式Dk 0 (k=1,n-1)，則則Gauss消元法中的約化主元可以表示為消元法中的約化主元可以表示為) 1,.,3 , 2(/1)(1, 1111nkDDaDakkkkk例例 用高斯消元法求解方程組用高斯消元法求解方程組 32563034253432321xxx 32303452536432202230445 . 0064324200445 . 006432 446245 . 0432321xxxx1= -13, x2 = 8, x3 = 2m

16、21=3/2m31=4/2m32=-3/0.5矩陣的三角分解矩陣的三角分解：1111211nmmL對線性方程組對線性方程組Ax=b的系數(shù)矩陣的系數(shù)矩陣A施行初等行變換相施行初等行變換相當(dāng)于用初等矩陣左乘當(dāng)于用初等矩陣左乘A，故第一次消元后方程組化，故第一次消元后方程組化為為A(1)x=b(1)，即，即L1A=A(1)x，L1b=b(1)，其中，其中同理同理LkA(k-1)=A(k) Lkb(k-1)=b (k)11111,1,knkkkmmL其中其中將將A 分解為單位下三角矩陣分解為單位下三角矩陣 L 和上三角矩陣和上三角矩陣 U的乘積的乘積的算法稱為矩陣的算法稱為矩陣A的三角分解算法。的三角

17、分解算法。重復(fù)該過程，最后得重復(fù)該過程，最后得)1(1221nnnAALLLL)1(1221nnnbbLLLL記記U=A(n-1)，則，則LUULLLAn111211其中其中1112121nnmmmL clear all;A=2 3 4;3 5 2;4 3 30;b=6 5 32;L,U=lu(A);x=U(Lb)x = -13.0000 8.0000 2.0000MATLAB 處理：處理：例例 對矩陣對矩陣A 作作LU分解。分解。解解 由由Gauss消去法可得，消去法可得， m21=0， m31=2， m32=1。故故A LU122140111112010001200140111如果已經(jīng)有如

18、果已經(jīng)有 A = L U 則則 AX = b = L U X = b，（1）求解方程組）求解方程組 LY = b 得向量得向量 Y 的值；的值； L 是下三角矩陣是下三角矩陣,用順代算法用順代算法（2）求解方程組）求解方程組 UX = Y 得向量得向量 X 的值。的值。U 是上三角矩陣是上三角矩陣, 用回代算法用回代算法記記UX = Y ， LY = b ，則求解方程組分兩步進(jìn)行：，則求解方程組分兩步進(jìn)行：3、Gauss列主元素消元法列主元素消元法000. 3000. 2000. 1643. 5072. 1000. 2623. 4712. 3000. 1000. 3000. 2001. 032

19、1xxx例：求解方程組例：求解方程組（用四位浮點(diǎn)數(shù)計(jì)算，精確解舍入到（用四位浮點(diǎn)數(shù)計(jì)算，精確解舍入到4位有效數(shù)字為：位有效數(shù)字為：x1*=-0.4904,x2*=-0.05104,x3*=0.3675）000. 3643. 5072. 1000. 2000. 2623. 4712. 3000. 1000. 1000. 3000. 2001. 020036006400101002300520040000. 1000. 3000. 2001. 0解：方法一解：方法一Gauss消去法消去法（A/b）=其中，其中， m21=-1.000/0.001=-1000 m31=-2.000/0.001=-20

20、00 m32=4001/2004=1.997 解為解為x1=-0.4000,x2=-0.09980,x3=0.4000(x1*=-0.4904,x2*=-0.05104,x3*=0.3675)顯然，此解并不準(zhǔn)確。顯然，此解并不準(zhǔn)確。000. 2000. 5001002300520040000. 1000. 3000. 2001. 0方法二：方法二：交換行，避免絕對值小的主元作除數(shù)。交換行，避免絕對值小的主元作除數(shù)。000. 1000. 3000. 2001. 0000. 2623. 4712. 3000. 1000. 3643. 5072. 1000. 2002. 1003. 3001. 20

21、5000. 0801. 1176. 30000. 3643. 5072. 1000. 2（A/b）=其中，其中， m21=0.5000 m31=-0.0005 m32=0.6300 解為解為x1=-0.4900,x2=-0.05113,x3=0.3678(x1*=-0.4904,x2*=-0.05104,x3*=0.3675)與方法一相比，此解顯然要精確得多。與方法一相比，此解顯然要精確得多。6870. 0868. 1005000. 0801. 1176. 30000. 3643. 5072. 1000. 23.1 基本思想：設(shè)基本思想：設(shè)Axb的增廣矩陣為的增廣矩陣為nnnnnnnbaaab

22、aaabaaa21222221111211在在A的第一列中選絕對值最大的元素作主元，設(shè)該元素的第一列中選絕對值最大的元素作主元，設(shè)該元素所在行為第所在行為第i1行，交換第一行與第行，交換第一行與第i1行，進(jìn)行第一次消行，進(jìn)行第一次消元；再在第元；再在第2n行的第二列中選絕對值最大的元素作行的第二列中選絕對值最大的元素作主元，設(shè)該元素所在行為第主元，設(shè)該元素所在行為第i2行，交換第二行與第行，交換第二行與第i2行，行，進(jìn)行第二次消元，進(jìn)行第二次消元，直到消元過程完成為止。直到消元過程完成為止。例：例：用列主元法解用列主元法解 6745150710623321xxx 6515707104623 6

23、515462370710解解：第一列中絕對值最大是：第一列中絕對值最大是10，取，取10為主元為主元第二列的后兩個(gè)數(shù)中選出主元第二列的后兩個(gè)數(shù)中選出主元 2.5 5 . 255 . 21 . 661 . 070710 2 . 62 . 65 . 255 . 270710 x3=6.2/6.2=1x2=(2.5-5x3)/2.5=-1x1=(7+7x2-0 x3)/10=0 x1=0 x2=-1x3=1function U,x=Gauss_column(A,b)% 列主元Gauss消去法求解線性方程組% 輸入?yún)?shù)：% -A：線性方程組的系數(shù)矩陣% -b：線性方程組的右端項(xiàng)% 輸出參數(shù)：% -U：

24、消元后的上三角方程組的增廣矩陣% -x：線性方程組的解n=length(b);for k=1:(n-1) % 選主元 Y,I=max(abs(A(k:n,k); I=I+k-1; if Ik t=A(k,:); A(k,:)=A(I,:); A(I,:)=t; t=b(k); b(k)=b(I); b(I)=t; end % 消元 m=A(k+1:n,k)/A(k,k); A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n); b(k+1:n)=b(k+1:n)-m*b(k); A(k+1:n,k)=zeros(n-k,1);endU=A,b;% 回代x=zer

25、os(n,1);x(n)=b(n)/A(n,n);for k=n-1:-1:1 x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)/A(k,k);end3.2 列主元矩陣的三角分解列主元矩陣的三角分解2) 11032(EEE 3) 11053(EEE 21EE 5150710623A 51562307101APA解：交換行變換解：交換行變換 55 . 261 . 00710111APFAP例：對矩陣?yán)簩仃嘇做列主做列主元三角分解：元三角分解： 1000010101P 100100131211mmF 0101000012P 10010001322mF3)25 . 21 . 03(E

26、EE 61 . 055 . 2071011211APFPAPF32EE 2 . 655 . 207101122APFPF則列主元的則列主元的Gauss變換可記為變換可記為A(2)=F2P2F1P1A記記 U=A(2)=F2(P2F1P2)P2P1A （因（因P2P2=I）P = P2P12121PFPF 則有則有1211FFL若記APFFU 12LUPA可得對于一般的對于一般的n階矩陣的列主元三角分解，通常令階矩陣的列主元三角分解，通常令12PPPPn)(1122111nnFPFPFPPL)(nAU LUPA 定理：定理：（列主元素的三角分解定理）若（列主元素的三角分解定理）若A非奇異，非奇異

27、，則存在排列陣則存在排列陣P使使PALU，其中，其中L為下三角陣，為下三角陣，U為上三角陣。為上三角陣。矩陣分解關(guān)系為矩陣分解關(guān)系為3.3 全主元素消元法全主元素消元法)1()1(maxkijnjknikkjiaakk定義：定義：則稱則稱) 1( kjikka為全主元素。為全主元素。經(jīng)過行列互換，使得經(jīng)過行列互換，使得 位于經(jīng)交換行和列位于經(jīng)交換行和列后的等價(jià)方程組中的后的等價(jià)方程組中的 位置，然后再實(shí)施消元。位置，然后再實(shí)施消元。)1( kjikka)1( kkka 全主元素消元法的基本思想：全主元素消元法的基本思想：若若注：全主元素消元法有可能改變未知數(shù)的順序。注：全主元素消元法有可能改變

28、未知數(shù)的順序。4、誤差分析、誤差分析例例 記方程組（記方程組（1） 為為Axb，其精確解為：，其精確解為：x1*=2，x2*=0現(xiàn)考察方程組（現(xiàn)考察方程組（2）可將其表示為：可將其表示為：A(x+ x)=b+ b，其中，其中 b= (0,0.0001)T ，x為（為（1）的解。）的解。 顯然（顯然（2）的）的解為：解為：x+ x= (1, 1)T 結(jié)論：結(jié)論：（1）的常數(shù)項(xiàng)）的常數(shù)項(xiàng)b的第二個(gè)分量只有的第二個(gè)分量只有1/1000的微小變化，方程組的解變化卻很大。的微小變化，方程組的解變化卻很大。220001. 111121xx0001. 220001. 111121xx定義定義 若矩陣若矩陣A或常數(shù)項(xiàng)或常數(shù)項(xiàng)b的微小變化引起方程的微小變化引起方程組組Ax=b的解的巨大變化，則稱此方程組為病的解的巨大變化，則稱此方程組為病態(tài)方程組，態(tài)方程組，A為病態(tài)矩陣（相對方程組而為病態(tài)矩陣（相對方程組而言）；否則稱方程組為良態(tài)方程組，言）；否則稱方程組為良態(tài)方程組，A為良為良態(tài)矩陣。態(tài)矩陣。研究方程組中研究方程組中A或或b的微小誤差對解的影響的分的微小誤差對解的影響的分析稱析稱“擾動分析擾動分析”。設(shè)設(shè)Axb的擾動方程組為的擾動方程組為