線性方程組的數(shù)值解法及其計(jì)算機(jī)實(shí)現(xiàn)

1、大連大學(xué)DALIAN UNIVERSITY2014屆畢業(yè)論文（設(shè)計(jì)）題目名稱：線性方程組的數(shù)值解法及其計(jì)算機(jī)實(shí)現(xiàn)所在學(xué)院： 信息工程學(xué)院 專業(yè)（班級(jí)）： 信息與計(jì)算科學(xué)10級(jí)1班 學(xué)生姓名： 李國慶 指導(dǎo)教師： 董學(xué)東 評(píng)閱人： 于萬波 院 長 ： 裴炳南 線性方程組的數(shù)值解法及其計(jì)算機(jī)實(shí)現(xiàn) 總計(jì)：畢業(yè)論文 43 頁 表 格 0 表 插 圖 11 幅指導(dǎo)老師 ：董學(xué)東（教授）評(píng) 閱 人 ：于萬波完成日期 ：2014年 5 月 21日線性方程組的數(shù)值解法及其計(jì)算機(jī)實(shí)現(xiàn)線性方程組的數(shù)值解法及其計(jì)算機(jī)實(shí)現(xiàn)摘 要現(xiàn)在在工程技術(shù)、自然和社會(huì)科學(xué)中常常會(huì)碰到許多問題，比如用最小二乘法來擬合實(shí)驗(yàn)數(shù)據(jù)的曲線、

2、電學(xué)中相關(guān)的網(wǎng)絡(luò)題目、計(jì)算建筑與機(jī)器中的策劃問題以及大地中的測量和經(jīng)濟(jì)中的投入產(chǎn)出問題，最后均可歸納為求解線性方程組。因此，在解決實(shí)際問題時(shí)，求解線性方程組起著不可替代的作用。由于在實(shí)際問題中，方程組的未知數(shù)個(gè)數(shù)多達(dá)上百個(gè)，用克萊姆法則解決起來過于麻煩，借助計(jì)算機(jī)來解決是一個(gè)非常可行而有效地方法。本文主要研究求解線性方程組的直接解法和迭代解法這兩種數(shù)值解法，并且利用Matlab編寫程序來實(shí)現(xiàn)這些方程組的求解。關(guān)鍵詞：線性方程組；數(shù)值解；直接法；迭代法；MatlabABSTRACTNow in engineering, natural sciences and social sciences,

3、many of the problems which frequently encountered can be attributed to the solution of linear equations final, such as the problems of curve fitting with experimental data for using the least squares method, network topics in the electrical, the design of building structures and the calculation of m

4、achine, measuring the earth and the economic operation of the input-output.Therefore, it is extremely significant to solve linear equations for practical problems. However, these equations are made of several hundred of unknown number, it will be viable and effective if solve linear equations using

5、computer than Cramer's rule. This article main researches the two numerical solutions of linear equations, including direct methods and iterative methods, and uses Matlab to achieve the solution of these equations.Keyword: Linear Equations; Numerical Solution; direct method; iterative methods; M

6、atlab目 錄摘要IABSTRACTII緒論11線性方程組數(shù)值解法的提出和定義21.1提出21.2定義22.線性方程組的直接解法及其計(jì)算機(jī)實(shí)現(xiàn)32.1 高斯消去法32.1.1 高斯消去法的定義及算法構(gòu)造32.1.2用Matlab 作高斯消去法52.1.3 列主元消去法的定義與算法62.1.4用Matlab 作列主元消去法72.2矩陣三角分解法82.2.1直接三角分解法的定義及算法構(gòu)造82.2.2 直接三角分解法的適用條件102.2.3 用Matlab 作直接三角分解法102.3平方根法與改進(jìn)的平方根法122.3.1平方根法的定義122.3.2改進(jìn)的平方根法的定義132.3.3 用Matlab

7、 作改進(jìn)的平方根法142.4 解三對(duì)角線方程組的追趕法152.4.1追趕法的定義152.4.2 用Matlab 實(shí)現(xiàn)追趕法163. 線性方程組的迭代解法及其計(jì)算機(jī)實(shí)現(xiàn)193.1 迭代法的定義193.2 雅克比迭代法193.2.1 雅克比迭代法定義193.2.2 用Matlab 實(shí)現(xiàn)雅克比迭代法203.3 高斯-賽德爾迭代法213.3.1高斯-賽德爾迭代法的定義213.3.2用Matlab 實(shí)現(xiàn)高斯-賽德爾迭代法223.4 超松弛迭代法233.4.1超松弛迭代法的定義233.5迭代法的收斂性253.5.1向量范數(shù)和矩陣范數(shù)253.5.2迭代法的收斂性判斷264.結(jié)論28參考文獻(xiàn)29附錄一30附錄

8、二38致謝43大連大學(xué)學(xué)位論文版權(quán)使用授權(quán)書44III緒論對(duì)于線性方程組的求解，在中國古代的九章算術(shù)里面，解線性方程組的“方程術(shù)”方法是世界上最完備、最先的求解方法，而劉徽則提出了相對(duì)系統(tǒng)的方程理論。線性方程組的研究在西方是在17世紀(jì)由萊布尼茨最先開始的。1-2隨著現(xiàn)在科學(xué)技術(shù)的不斷發(fā)展，有效地求解線性方程組在解決實(shí)際問題中越來越重要。由于在工程技術(shù)、科學(xué)研究中、實(shí)驗(yàn)設(shè)計(jì)等很多的方程組中，未知量往往過多，只能借助計(jì)算機(jī)才能快速解決，線性方程組的數(shù)值解法就是比較重要的方法。求解線性方程組的數(shù)值解的應(yīng)用非常廣泛，如：在建立物理現(xiàn)象的數(shù)學(xué)模型中，會(huì)間接出現(xiàn)一些數(shù)學(xué)模型的數(shù)值解。而在非線性方程組、常微

9、分方程邊值間題、最優(yōu)化問題數(shù)值解等問題都涉及到求解線性方程組，故在科技、醫(yī)藥和經(jīng)濟(jì)等各個(gè)領(lǐng)域，線性方程組的數(shù)值解法都有廣泛的應(yīng)用。3-5從參考文獻(xiàn)中可以看出，很多先者已對(duì)線性方程組的數(shù)值解法有所研究，比如關(guān)于病態(tài)線性方程組的迭代解法6、線性方程組解結(jié)構(gòu)的歷史研究、線性方程組的迭代解法及其Matlab實(shí)現(xiàn)程序7等。但是，幾乎都沒有對(duì)直接方法和迭代方法給出詳細(xì)的解法程序，又由于大學(xué)生在長期的學(xué)習(xí)中對(duì)于計(jì)算機(jī)的實(shí)踐環(huán)節(jié)往往忽視，不重視如何將編程語言與教材中給出的算法結(jié)合起來，這樣就很難對(duì)一些問題有較為深刻的理解，所以，使用計(jì)算機(jī)編程完成求解線性方程組，無論是對(duì)提高大學(xué)生實(shí)踐能力還是對(duì)解決實(shí)際問題都是

10、非常必要的。1線性方程組數(shù)值解法的提出和定義1.1提出求解線性方程組是代數(shù)里的非常重要的一個(gè)部分，廣泛應(yīng)用于數(shù)學(xué)以及其它科學(xué)領(lǐng)域。現(xiàn)代工程技術(shù)、科學(xué)研究中、實(shí)驗(yàn)設(shè)計(jì)等很多計(jì)算，最后都化為求解線性方程組，而這些方程組中的未知量往往多達(dá)幾百個(gè)，故不能用克萊姆法則來求解，只有借助計(jì)算機(jī)。因此，我們介紹兩種線性方程組的數(shù)值解法：直接法和迭代法，并進(jìn)行一些探討。1.2定義線性方程組的數(shù)值解法一般有兩類：直接法和迭代法直接法：便是通過有限步算術(shù)的運(yùn)算，可求得方程組的精確解的方法（如果計(jì)算過程當(dāng)中沒有舍入誤差），比如克萊姆法則就是一種直接法，高斯(Gauss)消去法是直接法最中具有代表性的算法。迭代法: 便

11、是用某種極限過程去逐步的逼近線性方程組的精確解的方法。也就是從解的某一個(gè)近似值出發(fā)，然后通過構(gòu)造一個(gè)無窮序列去逼近精確解的方法。8(一般有限步內(nèi)是得不到精確解)。2.線性方程組的直接解法及其計(jì)算機(jī)實(shí)現(xiàn)常見的線性方程組的形式是方程個(gè)數(shù)以及未知量個(gè)數(shù)相同的都為n階的線性方程組，一般形式為： 簡記為Ax=b，其中下面介紹幾種線性方程組的直接方法：高斯消去法、列主元消去法、直接三角分解法、改進(jìn)的平方根法和追趕法。2.1 高斯消去法2.1.1 高斯消去法的定義及算法構(gòu)造我們知道，線性方程組的矩陣表示為：第一步：設(shè)，把方程組的實(shí)數(shù)矩陣的第一列元素消為零，令用乘以第一個(gè)方程后加到第i個(gè)方程上去，消去第2n個(gè)

12、方程的未知數(shù)，得到，即其中：即為其中： 只有,消元過程才能夠繼續(xù)進(jìn)行下去,直到經(jīng)過n-1次消元后，消元過程結(jié)束后，會(huì)得到與原方程組等價(jià)的上三角形方程組,記為，或?qū)懗杉吹诙剑夯卮^程，便是對(duì)上三角方程組從下到上逐步回代求解方程組，即： 高斯消去法的算法： 若，覆蓋，1. 消去：對(duì)于，(1)若，則計(jì)算停止 (2)對(duì)于 (i) (ii)對(duì)于 2.回代：對(duì)于， (1) (2)對(duì)于，(3) 2.1.2用Matlab 作高斯消去法這是順序高斯消去法的M文件：function x=gauss(A,b)n=length(b);A=A,b;for k=1:(n-1) A(k+1):n,(k+1):(n+1)=

13、A(k+1):n,(k+1):(n+1). -A(k+1):n,k)/A(k,k)*A(k,(k+1):(n+1); A(k+1):n,k)=zeros(n-k,1);endx=zeros(n,1);x(n)=A(n,n+1)/A(n,n);for k=n-1:-1:1 x(k,:)=(A(k,n+1)-A(k,(k+1):n)*x(k+1):n)/A(k,k);end取A=1 -1 1;5 -4 3;2 1 1; b=-4 -12 11' 運(yùn)行結(jié)果為：圖2.12.1.3 列主元消去法的定義與算法一般線性方程組在用高斯消去法求解時(shí)，消元過程當(dāng)中有可能存在的情況，這時(shí)候消去法將沒法進(jìn)行；

14、即便，如果它的絕對(duì)值很小時(shí)，用其作除數(shù)，將會(huì)致使其余元素的數(shù)量級(jí)的嚴(yán)重增加以及舍入誤差的擴(kuò)散，這將嚴(yán)重影響計(jì)算結(jié)果的精度。在實(shí)際計(jì)算時(shí)必須避免這種情形的產(chǎn)生。而主元素消去法就可填補(bǔ)這個(gè)弊端。主元素：經(jīng)過方程或交換變量的次序，使其在對(duì)角線位置上獲取絕對(duì)值盡可能大的系數(shù)作為，稱這樣的 即為主元素。9-10主元素法：使用主元素的消元法為主元素法。按照主元素選取范圍分為：列主元素法、行主元素法、全主元素法。下面主要講解列主元素法。列主元素法：在待消元的所在列中選取主元，經(jīng)過方程的行交換，置主元素于對(duì)角線位置后再進(jìn)行消元的方法。列主元素的計(jì)算方法與高斯消去法完全的一樣,差別是在每步消元之前都要按列選出主

15、元。1. det1 2. 對(duì)于k=1,2，n-1（1） 按列選主元 （2） 假如=0，則計(jì)算停止(det(A)=0）（3） 如果=k則轉(zhuǎn)(4) 換行： （4） 消元計(jì)算 對(duì)于， 對(duì)于， （5）2.1.4用Matlab 作列主元消去法下面是列主元消去法的M文件：function x=gausslzy(A,b)n=length(b);A=A,b;for k=1:(n-1) Ap,p=max(abs(A(k:n,k);p=p+k-1; if p>k, t=A(k,:);A(k,:)=A(p,:);A(p,:)=t; end A(k+1):n,(k+1):(n+1)=A(k+1):n,(k+1)

16、:(n+1). -A(k+1):n,k)/A(k,k)*A(k,(k+1):(n+1); A(k+1):n,k)=zeros(n-k,1);endx=zeros(n,1);x(n)=A(n,n+1)/A(n,n);for k=n-1:-1:1 x(k,:)=(A(k,n+1)-A(k,(k+1):n)*x(k+1):n)/A(k,k);end取A=1 -1 1;5 -4 3;2 1 1; b=-4 -12 11' 運(yùn)行結(jié)果為：圖2.22.2矩陣三角分解法2.2.1直接三角分解法的定義及算法構(gòu)造矩陣三角分解法是以高斯消去法解線性方程組為基礎(chǔ)的一種變形解法。應(yīng)用高斯消去法去解n階線性方程組

17、Ax=b,再通過n步消元之后, 將會(huì)得出一個(gè)等價(jià)的上三角型方程組A(n) x=b(n),然后對(duì)上三角形方程組用逐步回代就可以求出解來。上述過程可通過矩陣分解來實(shí)現(xiàn)。將非奇異陣A分解成一個(gè)下三角陣L以及一個(gè)上三角陣U的乘積A=LU稱為對(duì)矩陣A的三角分解，又稱LU分解。11-13 其中 將A分解成一個(gè)單位上三角陣L和一個(gè)下三角陣U的乘積，稱為杜利特爾（Doolittle）分解。其中：若把A分解成一個(gè)下三角陣L和一個(gè)單位上三角陣U的乘積稱為克洛特分解Crout）。 其中 ：在求解線性方程組Ax=b時(shí),要先對(duì)非奇異矩陣A進(jìn)行LU分解使A=LU，則方程組就化為LUx=b。那么原方程組轉(zhuǎn)化為求解兩個(gè)簡單的

18、三角方程組：Ly=b求解y，Ux=y求解x。這就是求解線性方程組的三角分解法的基本思想。14下面只介紹杜利特爾（Doolittle）分解法。 設(shè)A=LU為： 由矩陣乘法規(guī)則：由此可得U的第1行元素以及L的第1列元素：再確定U的第k行元素和L的第k列元素,對(duì)于k=2,3, ,n計(jì)算： 計(jì)算U的第k行元素 （j=k,k+1,n） 計(jì)算L的第k列元素 （i=k,k+1,n） 利用上述計(jì)算公式即可逐步求出U與L的各元素求解 Ly=b , 即計(jì)算: 求解 Ux=y , 即計(jì)算：2.2.2 直接三角分解法的適用條件很明顯, 當(dāng)時(shí),解Ax=b直接三角分解法計(jì)算才有可能實(shí)現(xiàn)。設(shè)A為非奇異矩陣, 當(dāng)時(shí)計(jì)算將中斷

19、或當(dāng)絕對(duì)值很小時(shí)，按分解公式計(jì)算可能會(huì)引起舍入誤差的積累，于是可采用與列主元消去法相似的方法，對(duì)矩陣進(jìn)行行交換，則再實(shí)現(xiàn)矩陣的三角分解。2.2.3 用Matlab 作直接三角分解法下面是直接三角分解法的M文件：function x=LU(A,b)A=1 2 3;2 5 2;3 1 5;b=14 18 20'n=length(b);x=zeros(n,1);A(2:n,1)=A(2:n,1)./A(1,1);for i=2:n-1 A(i,i)=A(i,i)-sum(A(i,1:i-1)'.*A(1:i-1,i); for j=i+1:n A(i,j)=A(i,j)-sum(A(

20、i,1:i-1)'.*A(1:i-1,j); A(j,i)=(A(j,i)-sum(A(j,1:i-1)'.*A(1:i-1,i)/A(i,i); endendA(n,n)=A(n,n)-sum(A(n,1:n-1)'.*A(1:n-1,n);AU=A;L=A;for i=1:n L(i,i)=1;endfor i=1:n-1 for j=i+1:n L(i,j)=0; endendLfor i=2:n for j=1:i-1 U(i,j)=0; endendUy=zeros(n,1);y(1)=b(1);for i=2:n y(i)=b(i)-sum(L(i,1:i-

21、1)'.*y(1:i-1);endyx(n)=y(n)/U(n,n);for i=n-1:-1:1 x(i)=(y(i)-sum(U(i,i+1:n)'.*x(i+1)/U(i,i);endx運(yùn)行結(jié)果為：圖2.3圖2.42.3平方根法與改進(jìn)的平方根法在工程現(xiàn)實(shí)情況中,線性方程組的系數(shù)矩陣經(jīng)常具有對(duì)稱正定性，其各階順序主子式及全部特征值均大于0。矩陣的這一特性使它的三角分解具備更簡單的形式，從而導(dǎo)出一些特殊的解法，比如平方根法與改進(jìn)的平方根法。2.3.1平方根法的定義由于A是正定矩陣, A的順序主子式i>0, i=1,2,n ，因此存在惟一的分解A=LU，L是單位下三角陣,

22、 U是上三角陣, 將U再分解： 其中D為對(duì)角陣, U0為單位上三角陣，因此A = L U = L D U0，又A = AT = U0TD LT ，由分解惟一性, 即得U0T=L，A=L D LT 。記：又由于det(Ak)0,(k=1,2,n),因此，于是對(duì)角陣D還可分解：，其中為下三角陣。將A=LLT展開，寫成按矩陣乘法展開，可逐行求出分解矩陣L的元素，計(jì)算公式是對(duì)于i=1,2,n ， ， j=i+1, i+2,n ，這一方法稱為平方根法,又稱喬累斯基(Cholesky)分解。2.3.2改進(jìn)的平方根法的定義因?yàn)槠椒礁ń庹ǚ匠探M的弊端是要進(jìn)行開方運(yùn)算。為避免開方運(yùn)算，我們改用單位三角陣作為

23、分解陣，把對(duì)稱正定矩陣A分解成的形式，其中為對(duì)角陣，而是單位下三角陣,這里分解公式為利用這類矩陣分解法,方程組Ax=b即可歸結(jié)為求解兩個(gè)上三角方程組和。其計(jì)算公式分別為： 求解方程組的上述算法稱為改進(jìn)的平方根法。2.3.3 用Matlab 作改進(jìn)的平方根法下面是改進(jìn)的平方根法的M文件：function x=gaijin(A,b,n) A=2 -1 1;-1 -2 3;1 3 1;b=4 5 6'n=length(b);L=zeros(n,n); D=diag(n,0); S=L*D; for i=1:n L(i,i)=1;end for i=1:n for j=1:n if (eig(

24、A)<=0)|(A(i,j)=A(j,i) disp('wrong');break;endendend D(1,1)=A(1,1); for i=2:n for j=1:i-1 S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)'); L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1); end D(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');end y=zeros(n,1); x=zeros(n,1);for i=1:n y(i)=(b(i)-sum(L(i,1:i-1

25、)*D(1:i-1,1:i-1)*y(1:i-1)/D(i,i); endfor i=n:-1:1 x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n); end運(yùn)行結(jié)果為：圖2.52.4 解三對(duì)角線方程組的追趕法2.4.1追趕法的定義在數(shù)值計(jì)算中,有一種系數(shù)矩陣是三對(duì)角方程組： 將上述方程組簡記為Ax=f,其中A滿足條件(對(duì)角占優(yōu)） （1）（2）（3）按乘法展開則可計(jì)算可依次計(jì)算：方程組Ax=f等價(jià)解兩個(gè)方程組Ly=f和Ux=y。從而得到解三角形方程組的追趕公式：（1） 計(jì)算的遞推公式 =/， （2） 解Ly=f ， （3） 解Ux=y ， 我們將計(jì)算系數(shù)以及的過程稱

26、為追的過程，將計(jì)算方程組的解得過程稱為趕的過程。152.4.2 用Matlab 實(shí)現(xiàn)追趕法下面是追趕法的M文件：function x=zhuigan(A,b)A=2,-1,0,0,0;-1,2,-1,0,0;0,-1,2,-1,0;0,0,-1,2,-1;0,0,0,-1,2;d=1 0 0 0 0'n=length(d);U=zeros(n);L=eye(n);y=zeros(n,1);x=y;U(1,1)=A(1,1);y(1)=d(1);for i=2:n L(i,i-1)=A(i,i-1)/U(i-1,i-1); U(i,i)=A(i,i)-L(i,i-1)*A(i-1,i);

27、 y(i)=d(i)-L(i,i-1)*y(i-1); U(i-1,i)=A(i-1,i);endLUyx(n)=y(n);for i=n-1:-1:1 x(i)=(y(i)-A(i,i+1)*x(i+1)/U(i,i);endx運(yùn)行結(jié)果為：圖2.6圖2.73. 線性方程組的迭代解法及其計(jì)算機(jī)實(shí)現(xiàn)3.1 迭代法的定義迭代法基本思想是將線性方程組轉(zhuǎn)化成便于迭代的等價(jià)方程組，對(duì)任選的一組初始值，按照某種計(jì)算規(guī)則，不斷對(duì)所得到的值進(jìn)行修正，最后將得到滿足精度要求的方程組的近似解。設(shè)非奇異，線性方程組有惟一解，通過變換構(gòu)造出一個(gè)等價(jià)同解方程組，將上式改寫成迭代式: 選定初始向量,重復(fù)地使用迭代式去逐步

28、逼近方程組的精確解,直到滿足精度要求為止。這種方法稱作迭代法。16假如存在極限，則稱迭代法是收斂的，不然便是發(fā)散的。收斂時(shí)，在迭代公式中，當(dāng)時(shí)，則，故是方程組的解。3.2 雅克比迭代法3.2.1 雅克比迭代法定義考察一般的方程組，將n元線性方程組，據(jù)此建立迭代公式：即：設(shè)方程組的系數(shù)矩陣A是非奇異的，并有主對(duì)角元素，那么可將A分裂成：記作：A = D - L - U 則Ax=b等價(jià)于，即，因?yàn)?，則，這樣便得到一個(gè)迭代公式，令，則有（k = 0,1,2），稱為雅可比迭代公式, B稱為雅可比迭代矩陣。3.2.2 用Matlab 實(shí)現(xiàn)雅克比迭代法下面是雅克比迭代法的M文件：e=100;disp(&#

29、39;輸入方程組維數(shù)n');n=input('n');for(i=1:n) for(j=1:n) fprintf('輸入方程組第%d行第%d列的系數(shù)（從左到右,從上到下） n',i,j); a(i,j)=input('n'); endend for(i=1:n)fprintf('輸入等號(hào)右邊第%d個(gè)系數(shù)n',i);y(i,1)=input('n');end for(i=1:n)fprintf('輸入第%d個(gè)變量的初值n',i);x(i,1)=input('n');end fp

30、rintf('輸入精確度n');p=input('n'); for(i=1:n) for(j=1:n) if i=j b(i,j)=-a(i,j)/a(i,i); else b(i,j)=0; end g(i,1)=y(i,1)/a(i,i); end g(i,1)=y(i,1)/a(i,i);end while e>p|e<-p x1=b*x+g; e1=x1(1,1)-x(1,1); e2=x1(2,1)-x(2,1); e3=x1(3,1)-x(3,1); e4=min(e1,e2); e5=min(e1,e3); e=min(e4,e5);

31、x=x1;end for(i=1:n) fprintf('x%d=%fn',i,x(i,1);end運(yùn)行結(jié)果是：圖3.13.3 高斯-賽德爾迭代法3.3.1高斯-賽德爾迭代法的定義在雅克比的迭代法中，每次迭代只用到前一次的迭代值，若每次迭代充分利用當(dāng)前的最新的迭代值，即在求時(shí)用新分量取代舊分量就會(huì)得到高斯-賽德爾迭代法。其迭代法格式為： (i=1,2,n k=0,1,2,)將A分裂成A =D-L-U，則等價(jià)于(D-L-U )x = b，則高斯塞德爾迭代過程，故，令，則則高斯-塞德爾迭代形式為：。3.3.2用Matlab 實(shí)現(xiàn)高斯-賽德爾迭代法下面是高斯賽德爾迭代法的M文件：fu

32、nction x=GSdiedai(A,b)A= 8 -3 2;4 11 -1;6 3 12;b=20 33 36'N=length(b); fprintf('庫函數(shù)計(jì)算結(jié)果：');x=inv(A)*bx=zeros(N,1);D=diag(diag(A);L=-tril(A,-1); U=-triu(A,1); B=inv(D-L)*U;g=inv(D-L)*b;eps=0.001;for k=1:100 fprintf('第%d次迭代',k); y=B*x+g; if norm(x-y)<eps break; end x=yendx運(yùn)行結(jié)果為：

33、圖3.23.4 超松弛迭代法3.4.1超松弛迭代法的定義超松弛迭代法的目的是為了提升迭代法的收斂速度，在高斯塞德爾迭代公式的基礎(chǔ)上做一些修改。這種方法是將前一步結(jié)果和高斯-塞德爾迭代方法的迭代值適當(dāng)?shù)募訖?quán)平均，希望能得到更好的近似值。超松弛迭代法是解大型稀疏矩陣方程組的有效方法之一，有著廣泛的應(yīng)用。17-18用高斯塞德爾迭代法定義輔助量：把取為與的加權(quán)平均，即：合并表示為：其中稱為松弛因子，當(dāng)=1時(shí)，便是高斯-塞德爾迭代法。為了保證迭代過程收斂性，必須要有0<< 2。 當(dāng)0<< 1時(shí)，低松弛法；當(dāng)1<< 2時(shí)稱為超松弛法。但通常統(tǒng)稱為超松弛法(SOR)。19

34、設(shè)線性方程組Ax=b的系數(shù)矩陣A為非奇異，且主對(duì)角元素, 則將A分裂成A=d-L-U, 那么超松弛迭代公式用矩陣表示為： 或故 明顯，對(duì)任何一個(gè)值,(D+L)非奇異, (因?yàn)榧俣?于是超松弛迭代公式為：。令 則超松弛迭代可寫成：。3.4.2用Matlab 實(shí)現(xiàn)超松弛迭代法下面是超松弛迭代法的M文件：function x=chaosongchi(A,b)A=8,-3,2;4,11,-1;6,3,12;b=20,33,36'N=length(b); fprintf('庫函數(shù)計(jì)算結(jié)果：');x=inv(A)*b x=zeros(N,1);D=diag(diag(A);L=-t

35、ril(A,-1);U=-triu(A,1);w=1.5;B=inv(D-w*L)*(1-w)*D+w*U;g=w*inv(D-w*L)*b;eps=0.00001;for k=1:100 fprintf('第%d次迭代： ',k); y=B*x+g; if abs(x-y)<eps break; end x=yendx當(dāng)w=1時(shí)，運(yùn)行結(jié)果和高斯-賽德爾迭代法一樣為：圖3.3當(dāng)w=1.5時(shí)，運(yùn)行結(jié)果為：圖3.43.5迭代法的收斂性在給出迭代法的收斂性判斷之前，我們首先給出一些基礎(chǔ)性知識(shí)：向量范數(shù)和矩陣范數(shù)。3.5.1向量范數(shù)和矩陣范數(shù)向量范數(shù)：假如向量xÎRn

36、（或）的某個(gè)實(shí)值函數(shù)N(x)= |x|，滿足下列條件條件：(1) |x|0（|x|=0當(dāng)且僅當(dāng)x=0）（正定條件）(2) |x|= |x|，（或）(3) （三角不等式）則稱N(x)是Rn（或）上的一個(gè)向量范數(shù)（或模）向量范數(shù)有以下幾種常用的類型：（1）向量的范數(shù)：（2）向量的1-范數(shù)：（3）向量的2-范數(shù)：（4）向量的p-范數(shù)：矩陣范數(shù)：假如矩陣的某個(gè)非負(fù)的實(shí)值函數(shù)，滿足條件：（1）且僅當(dāng)A=0時(shí)，（正定條件）（2），為實(shí)數(shù)（齊次條件）（3）對(duì)于任意的矩陣A,B有（三角不等式）（4）則稱是上的一個(gè)矩陣范數(shù)(或模)。矩陣范數(shù)有以下幾種類型：設(shè)，（1）（稱為A的行范數(shù)） （2）（稱為A的列范數(shù)）（

37、3）（稱為A的2-范數(shù)），其中表示的最大特征值，即。3.5.2迭代法的收斂性判斷我們知道, 對(duì)于給定的方程組能夠構(gòu)造成簡單的迭代公式如：雅可比迭代公式、高斯-塞德爾迭代公式和超松弛迭代公式，但并不是一定會(huì)收斂。現(xiàn)在分析它們的收斂性。對(duì)于方程組經(jīng)過等價(jià)變換構(gòu)造出的等價(jià)方程組。定理1. 迭代公式收斂的充分必要條件是迭代矩陣G的譜半徑。由此定理可知，不論是雅可比迭代法、高斯 塞德爾迭代法還是超松弛迭代法，它們收斂的充要條件是其迭代矩陣的譜半徑。定理2. 若迭代矩陣G的一種范數(shù),則迭代公式收斂。嚴(yán)格對(duì)角占優(yōu)矩陣：設(shè)A=，假如A的元素滿足 >,稱A為嚴(yán)格對(duì)角占優(yōu)矩陣。弱對(duì)角占優(yōu)矩陣：設(shè)A=，假如A

38、的元素滿足 ,，且上式至少有一個(gè)不等式嚴(yán)格成立，則稱A為弱對(duì)角占優(yōu)矩陣?？杉s與不可約矩陣：設(shè)A= ，如果存在置換矩陣P使 ，其中為r階矩陣，為n-r階矩陣（），則A為可約矩陣，否則如果不存在這樣的置換矩陣P是其成立，則稱A為不可約矩陣。定理3. 設(shè)Ax=b,如果A為嚴(yán)格對(duì)角占優(yōu)矩陣，則解Ax=b的雅克比迭代法，高斯賽德爾迭代法均收斂。如果A為弱對(duì)角占優(yōu)矩陣，且A為不可約矩陣，則解Ax=b的雅克比迭代法，高斯賽德爾迭代法均收斂。定理4. 設(shè)矩陣A對(duì)稱，且對(duì)角元，則解線性方程組Ax=b的雅克比迭代法收斂的充分必要條件是A及2D-A均為正定矩陣，其中。解線性方程組Ax=b的高斯-賽德爾法收斂的充分條

39、件是A正定。定理5. 對(duì)于線性方程組Ax=b，若A為對(duì)稱正定矩陣，則當(dāng)0<<2時(shí)，超松弛迭代法收斂。定理6. 對(duì)于線性代數(shù)方程組Ax=b， 若A按行(或列)嚴(yán)格對(duì)角占優(yōu)，或按行(或列)弱對(duì)角占優(yōu)不可約；則當(dāng)0<w1時(shí)，超松弛迭代法收斂。4.結(jié)論求解線性方程組是代數(shù)里的一個(gè)非常重要的部分?，F(xiàn)代工程技術(shù)、科學(xué)研究中、實(shí)驗(yàn)設(shè)計(jì)等很多計(jì)算，最后都化為求解線性方程組。本論文總結(jié)了直接方法中的高斯消去法、列主元消去法、直接三角分解法、改進(jìn)的平方根發(fā)和追趕法；迭代方法中的雅克比迭代法、高斯塞德爾迭代法、和超松弛迭代法，并討論迭代法的收斂性和收斂速度的問題。在此基礎(chǔ)上應(yīng)用Matlab程序?qū)崿F(xiàn)

40、這些方法，給出了直接方法和迭代方法詳細(xì)的解法程序，這對(duì)提高大學(xué)生實(shí)踐能力和解決實(shí)際問題有所幫助。參考文獻(xiàn)1李慶揚(yáng).數(shù)值分析第五版M.清華大學(xué)出版社，2008.2丙辛. 線性方程組的解法J.數(shù)學(xué)通報(bào)，1961，（1）:32-36.3曹學(xué)濱. 線性方程組解結(jié)構(gòu)的歷史研究D.山東:山東大學(xué)數(shù)學(xué)系基礎(chǔ)數(shù)學(xué)，2009.4王永寶.關(guān)于病態(tài)線性方程組的迭代解法J.航空計(jì)算技術(shù)，1995，(3): 20-25.5周建興，豈興明，矯津毅，張延偉.MATLAB從入門到精通（第2版）M.人民郵電出版社，2012.6 Wilkinson J H. Rounding errors in algebraic prenti

41、ces. London: H. M. Stationery office, 1963.7Golub G H, Van Loan C F. Matrix computations. 3rd ed. Johns Hopkins University Press, Baltimore MD, 1996.8Young D M. Iterative solution of large linear systems. New York: Academic, 1971.9Hackbusch W. Iterative solution of large sparse systems of equation.

42、New York: Springer-Verlag, 1994.10姜啟源.數(shù)學(xué)建模方法及其應(yīng)用M.北京:高等教育出版社，2003.11杜廷松.數(shù)值分析及實(shí)驗(yàn)M.科學(xué)出版社，2006.12陳延梅.計(jì)算方法指導(dǎo)M.科學(xué)出版社，2003.13龐遼軍,姜正濤,王育民. 基于一般訪問結(jié)構(gòu)的多重秘密共享方案J.計(jì)算機(jī)研究與發(fā)展， 2006, 43(1): 33-38. 14曹戈,MATLAB教程及實(shí)訓(xùn)M,機(jī)械工業(yè)出版社,2014.15周義倉、赫孝良.數(shù)學(xué)建模實(shí)驗(yàn)M.西安:西安交通大學(xué)出版社，1999.16黃進(jìn)陽. 論線性方程組的數(shù)值解法J. 邵陽高專學(xué)報(bào)，1995，(4): 292-297.17楊鳳霞.

43、 線性方程組的數(shù)值解法J. 滄州師范?？茖W(xué)校學(xué)報(bào)，2000，(3): 38-40.18 Xuedong Dong，Liqiang Hou，Yan Zhang，Public Key Cryptosystem and Signature Schemes over the Ring of Eisenstein Integers，2013 Fourth International Conference on Intelligent Control and Information Processing (ICICIP) June 9 11, 2013, Beijing, China.19童小紅,秦新強(qiáng).

44、線性方程組的解法探討及MATLAB實(shí)現(xiàn)J.數(shù)學(xué)學(xué)習(xí)與研究，2013，(13): 23-25.20花威.線性方程組的迭代解法及其Matlab實(shí)現(xiàn)程序J. 長江工程職業(yè)技術(shù)學(xué)院學(xué)報(bào),2009,(4): 12-14.附錄一 在愛森斯坦整數(shù)環(huán)下的公鑰密碼體制和簽名方案董學(xué)東摘要在本文中，使用愛森斯坦整環(huán)模愛森斯坦素?cái)?shù)的有限域的乘法群，給出了恢復(fù)信息的擴(kuò)展的ElGamal公鑰密碼體制和數(shù)字簽名方案。這種方法比經(jīng)典的ElGamal密碼體制和數(shù)字簽名方案的優(yōu)點(diǎn)多。文中給出許多例子說明所給出的ElGamal公鑰密碼體制和數(shù)字簽名方案的有效性。一引言基于整數(shù)環(huán)模素?cái)?shù)的單位群的離散對(duì)數(shù)問題，ElGamal1提出公

45、鑰密碼體制和數(shù)字簽名方案。ElGamal加密方案和數(shù)字簽名方案可以很容易地推廣到任何有限循環(huán)群中。群G應(yīng)慎重選擇，使在群G中的運(yùn)算相對(duì)容易以便提高效率。在計(jì)算G中離散對(duì)數(shù)時(shí)，樸素的算法是增加G的生成元的冪，直到所需的元素被找到。該算法需要G群的線性空間，因此是群G的指數(shù)冪。因此，若G的階是足夠大，G中的離散對(duì)數(shù)問題應(yīng)該是計(jì)算上不可行的，以確保使用的ElGamal公鑰密碼體制和數(shù)字簽名方案的協(xié)議的安全性。在密碼學(xué)中最感興趣的群是有限域的乘法群。把ElGamal公鑰密碼體制推廣到高斯整數(shù)環(huán)已在文2中給出。ELKamchouchi等在文34中提出了高斯整數(shù)環(huán)上的擴(kuò)展RSA密碼體制和數(shù)字簽名方案。注意

46、到高斯整數(shù)環(huán)是分圓域的代數(shù)整數(shù)環(huán)，艾森斯坦整數(shù)環(huán)是分圓域的代數(shù)整數(shù)環(huán)，我們想知道ElGamal公鑰密碼體制和數(shù)字簽名方案是否可以擴(kuò)展到艾森斯坦整數(shù)環(huán)。愛森斯坦的環(huán)像整數(shù)環(huán)和高斯整數(shù)環(huán)有許多良好的性質(zhì)（5,6,7）。在艾森斯坦整數(shù)環(huán)求最大公約數(shù)有高效的算法8。本文給出了艾森斯坦整數(shù)環(huán)上的ElGamal公鑰密碼體制和數(shù)字簽名方案。這種方法具有很多優(yōu)點(diǎn)。首先，如果在擴(kuò)展的ElGamal密碼系統(tǒng)中使用的循環(huán)群是艾森斯坦整數(shù)環(huán)模具有形式的素?cái)?shù)，則其階是比在古典ElGamal密碼所使用的平方大。其次，由于艾森斯坦整數(shù)環(huán)的運(yùn)算，一個(gè)素?cái)?shù)的擴(kuò)展歐拉函數(shù)是，而在整數(shù)環(huán)中相應(yīng)的歐拉函數(shù)是。由于加冪指數(shù)是被選擇的與

47、擴(kuò)展的歐拉函數(shù)互質(zhì)，所以它比經(jīng)典的ElGamal密碼系統(tǒng)更安全。第三，如果模數(shù)是兩個(gè)非實(shí)的艾森斯坦素?cái)?shù)的乘積，分解模的難度提高，因?yàn)樵诎固拐麛?shù)環(huán)中，元素的分解比在整數(shù)環(huán)元素的分解更復(fù)雜。最后，所使用的運(yùn)行模式與使用整數(shù)環(huán)是相似的，可以使用負(fù)整數(shù)。文中給出許多例子說明所給出的ElGamal公鑰密碼體制和數(shù)字簽名方案的有效性。本文的其余部分安排如下：在第2節(jié)，我們給出關(guān)于艾森斯坦整數(shù)環(huán)的一些預(yù)備知識(shí)。在第3節(jié)，我們?cè)诎固拐麛?shù)環(huán)構(gòu)建ElGamal公鑰密碼體制。在第4節(jié)，在艾森斯坦整數(shù)環(huán)中給出了擴(kuò)展ElGamal數(shù)字簽名方案。最后，第5節(jié)給出了結(jié)束語。二預(yù)備知識(shí)定義域是一個(gè)包括兩個(gè)操作：，叫加

48、，和，所謂乘法的集合，滿足下列公理。如果在+下有零元，并記為零，集合是一個(gè)阿貝爾群。的所有非零元素的集合*，如果有幺元，記為1，那么也是阿貝爾群，在除法之上的乘法分布。 滿足相同公理的乘法分布交換環(huán)作為一個(gè)定義域，除了非零元素不一定有乘法的逆。一個(gè)交換環(huán)的元素被稱為劃分如果中的元素（符號(hào)：）存在使。中的元素被稱為有聯(lián)系的，如果并且，交換環(huán)R的元素是素?cái)?shù)，條件是：（1）是非零并且沒有任何乘法逆；(2) 或。艾森斯坦整數(shù)環(huán)是,是它的 三重復(fù)根。一個(gè)艾森斯坦整數(shù)的范數(shù)被定義為，表示它的共軛復(fù)數(shù)，有6個(gè)元逆素，。元素與有聯(lián)系。愛森斯坦素?cái)?shù)和它的相關(guān)素?cái)?shù)，理性素?cái)?shù)和它的相關(guān)素?cái)?shù)，考慮的理性素?cái)?shù)

49、，其中。前幾個(gè)艾森斯坦素?cái)?shù)等于一個(gè)天然素3k+2：2, 5, 11, 17, 23, 29, 41, 47, 53, 59, 71, 83, 89, 101, 107, 113, 131, 137, 149, 167, 173, 179, 191, 197, 227, 233, 239, 251, 257, 263, 269, 281, 293, 311, 317, 347, 353, 359, 383, 389, 401, 419, 431, 443, 449, 461, 467, 479, 491, 503, 509, 521, 557, 563, 569, 587。一些非實(shí)數(shù)的艾森斯坦素

50、數(shù)是。下圖1顯示了小艾森斯坦素?cái)?shù)9。那些綠色軸是與天然素?cái)?shù)形式是相關(guān)聯(lián)的。其它的是絕對(duì)值平方等于天然素?cái)?shù)。一個(gè)艾森斯坦整數(shù)表達(dá)成艾森斯坦素?cái)?shù)的乘積是獨(dú)一無二的，除了素?cái)?shù)的順序，相關(guān)素?cái)?shù)之間的逆元素存在，存在著歧義。對(duì)于，其中是有理數(shù)，令，符號(hào)表示最大整數(shù)小于或等于，。接下來我們用來表示，令，表示，其中，是的共軛復(fù)數(shù)。如果是艾森斯坦素?cái)?shù)，這樣，是有理元素，然后這很容易證明是和的一個(gè)同構(gòu)。特別的，如果是的一個(gè)原始的元素，那么是的一個(gè)原始的元素。費(fèi)馬大定理對(duì)的模擬如下。定理1：讓是艾森斯坦素?cái)?shù)，其中和沒有直接關(guān)系，。然后對(duì)于任何艾森斯坦整數(shù)成立，其中，對(duì)于任何艾森斯坦整數(shù)成立。證明，應(yīng)為是艾森斯坦素

51、數(shù)，商環(huán)分別是有限域與元素和10,p.14。假設(shè)是艾森斯坦整數(shù)，注意，艾森斯坦整數(shù)環(huán)是歐幾里德域。因此，任何兩個(gè)整數(shù)愛森斯坦有最大公約數(shù)。如果，那么，因此，因?yàn)楹蜎]有直接關(guān)系，它遵從,因此。如果或2，然后，。因此我們有，如果，那么。這就是定理的證明。三在艾森斯坦整數(shù)環(huán)下的ElGamal公鑰密碼體制密碼學(xué)的基本目標(biāo)是為了讓兩個(gè)人，通常被稱為Alice和Bob，在不安全的信道進(jìn)行交流，這種方法使對(duì)手Eve不明白你在說什么。ElGamal加密系統(tǒng)中基于離散對(duì)數(shù)問題。離散對(duì)數(shù)問題在有限域內(nèi)一直是許多研究的對(duì)象。如果質(zhì)數(shù)是經(jīng)過精心選擇的，該問題通常被認(rèn)為是困難的。特別是，沒有已知的多項(xiàng)式的時(shí)間算法的離散對(duì)數(shù)問題。為了阻止已知的攻擊，應(yīng)該至少有一個(gè)