利用直接發(fā)求解線性方程組及其應(yīng)用.doc

利用直接法解線性方程組及其應(yīng)用摘要: 在科學(xué)研究和工程技術(shù)中有許多問題可歸結(jié)為求解線性代數(shù)方程組的問題和常微分方程的定解問題。本文主要討論了解線性方程組的各種直接解法，常微分方程的邊值問題及其應(yīng)用，最后介紹了病態(tài)方程和及其求解方法。關(guān)鍵字: 線性方程組 直接法 應(yīng)用舉例一、引言隨著科學(xué)技術(shù)的發(fā)展，提出了大量復(fù)雜的數(shù)值計(jì)算問題，在建立電子計(jì)算機(jī)成為數(shù)值計(jì)算的主要工具以后，它以數(shù)字計(jì)算機(jī)求解數(shù)學(xué)問題的理論和方法為研究對象。有可靠的理論分析，要有數(shù)值實(shí)驗(yàn)，并對計(jì)算的結(jié)果進(jìn)行誤差分析。數(shù)值分析的主要內(nèi)容包括插值法，函數(shù)逼近，曲線擬和，數(shù)值積分，數(shù)值微分，解線性方程組的直接方法，解線性方程組的迭代法，非線性方程求根，常微分方程的數(shù)值解法。運(yùn)用數(shù)值分析解決問題的過程包括：實(shí)際問題數(shù)學(xué)建模數(shù)值計(jì)算方法程序設(shè)計(jì)上機(jī)計(jì)算求出結(jié)果。在自然科學(xué)研究和工程技術(shù)中有許多問題可歸結(jié)為求解線性代數(shù)方程組的問題,線性方程組求解是科學(xué)計(jì)算中最常遇到的問題。如在應(yīng)力分析、電路分析、分子結(jié)構(gòu)、測量學(xué)中都會遇到解線性方程組問題。在很多廣泛應(yīng)用的數(shù)學(xué)問題的數(shù)值方法中，如三次樣條、最小二乘法、微分方程邊值問題的差分法與有限元法也都涉及到求解線性方程組。直接法是在沒有舍入誤差的情況下，通過有限步四則運(yùn)算來求的方程組精確解的方法。直接法包括高斯列主元消去法，三角分解法，追趕法等的基本思想和原理，學(xué)習(xí)它們各自的優(yōu)缺點(diǎn)及適用范圍。能夠利用直接法解線性方程組解決一些如道路的交通流量，電網(wǎng)的電流流量，解析幾何等的實(shí)際問題。二、解線性方程組的方法關(guān)于線性代數(shù)方程組的數(shù)值解法一般分為兩大類：直接法和迭代法。1、 直接法 直接法就是經(jīng)過有限步算術(shù)運(yùn)算，可求的線性方程組精確解的方法（若計(jì)算過程沒有舍入誤差），但實(shí)際猶如舍入誤差的存在和影響，這種方法也只能求得近似解，這類方法是解低階稠密矩陣方程組級某些大型稀疏矩陣方程組的有效方法。直接法包括高斯消元法，矩陣三角分解法、追趕法、平方根法。2、 迭代法迭代法就是利用某種極限過程去逐步逼近線性方程組精確解的方法。將方程組的解看作是某極限過程的極限值，且計(jì)算這一極限值的每一步是利用前一步所得結(jié)果施行相同的演算步驟而進(jìn)行。迭代法具有需要計(jì)算機(jī)的存儲單元少，程序設(shè)計(jì)簡單，原始系數(shù)矩陣在計(jì)算過程始終不變等優(yōu)點(diǎn)，但存在收斂性級收斂速度問題。迭代法是解大型稀疏矩陣方程組（尤其是微分方程離散后得到的大型方程組）的重要方法。迭代法包括Jacobi法SOR法、SSOR法等多種方法。本文主要介紹用高斯消元法，矩陣三角分解法求解線性方程組及其應(yīng)用。三、具體解法及其實(shí)例1、Gauss消元法通過一系列的加減消元運(yùn)算，也就是代數(shù)中的加減消去法，以使A對角線以下的元素化為零，將方程組化為上三角矩陣；然后，再逐一回代求解出x向量。現(xiàn)舉例說明如下：1.1消元過程第一步：將(1)/3使x1的系數(shù)化為1,再將(2)、(3)式中x1的系數(shù)都化為零，即由(2)-2(1)(1)得 由(3)-4(1)(1)得 第二步：將(2)(1)除以2/3，使x2系數(shù)化為1，得再將(3)(1)式中x2系數(shù)化為零，由(3)(1)-(-14/3)*(2)(2) ，得第三步：將(3)(2)除以18/3，使x3系數(shù)化為1，得經(jīng)消元后，得到如下三角代數(shù)方程組：1.2 回代過程由(3)(3)得 x3=1,將x3代入(2)(2)得x2=-2,將x2 、x3代入(1)(1)得x2=1,所以，本題解為x=1,2,-1T 1.3 用矩陣演示進(jìn)行消元過程第一步： 先將方程寫成增廣矩陣的形式第二步：然后對矩陣進(jìn)行初等行變換第三步：將增廣矩陣變換成上三角矩陣，即主對角線全為1，左下三角矩陣全為0，形式如下：即原方程組被等價(jià)轉(zhuǎn)化成為上三角方程組，然后，逐步回代得原方程組的解即可。1.4高斯消元的公式綜合以上討論，不難看出，高斯消元法解方程組的公式為第一步，消元（1） 令aij(1) = aij , （i,j=1,2,3,n）bi(1) =bi ,（i=1,2,3,n）（2） 對k=1到n-1，若akk(k)0，進(jìn)行l(wèi)ik = aik (k) / akk(k) ,（i=k+1,k+2,n）aij(k+1) = aij(k) - lik * akj(k), （i,j= k+1,k+2,n）bi(k+1) = bi(k) - lik * bk(k), （i= k+1,k+2,n）第二步，回代若ann(n) 0xn = bn(n) / ann(n)xi = (bi(i) sgm(aij(i) * xj )/- aii(i) ,（i = n-1,n-2,1）,( j = i+1,i+2,n )1.5 解的判斷 設(shè)方程組的增廣矩陣記為，則經(jīng)過初等行變換可化為如下的階梯形矩陣(必要是可重新排列未知量的順序)：其中cii0(i=1，2，r)于是可知：(1). 當(dāng)dr+1=0，且r=n時(shí)，原方程組有唯一解(2). 當(dāng)dr+1=0，且rk時(shí)，lkr=0, 且lkk=1，因?yàn)樗?，同理可推出?jì)算L的第k列的公式：因此得到如下算法杜利特（Doolittle）算法:(1) 將矩陣分解為A=LU,對k=1,2,n；j=k,k+1,n； i=k,k+1,n；公式1 (2) 解Ly=b (3) 解Ux=y對大規(guī)模稀疏問題，如果能夠通過調(diào)整方程及未知量的順序使得方程組的系數(shù)矩陣成帶狀結(jié)構(gòu)，則對系數(shù)矩陣使用通常的LU分解，可以保障單位下三角矩陣L及上三角矩陣U仍為帶狀結(jié)構(gòu).3. 直接法解線性方程組的應(yīng)用例1：確定下圖電網(wǎng)中的回路電流。解 在回路1中,電流I1流過三個(gè)電阻,且電壓降為；在回路2中的電流也流經(jīng)回路1的一部分,即從D到A的分支,對應(yīng)的電壓降為6I2伏特然而,回路1中電流在DA段的方向與回路2中選定的方向相反,因此,回路1中所有電壓降的代數(shù)和為由于回路1中的電壓為+60伏特,由基爾霍夫電壓定律,可得回路1的方程為,同理,可得回路2的方程為 , 其中, -6I1是回路1中流經(jīng)DA分支的電流(因?yàn)殡娏髋c回路2中的電流方向相反,所以電壓為負(fù))；12I2是回路2中所有的電阻乘上回路電流的和； -2I3是回路3中流經(jīng)CB分支上2歐姆電阻的電流,方向與回路2中該段的電流方向相反。回路3的方程為 注意,在CB分支上10伏特的電池被當(dāng)作是回路2和回路3中的一部分,但是由于回路3中電流方向,電池在回路3中為-10伏特出于同樣的道理,40伏特的電池也應(yīng)取負(fù)值.綜合上述討論,上述電網(wǎng)的回路電流滿足下列線性方程組 寫成矩陣形式為 (*)對增廣矩陣進(jìn)行行變換,得從而解得I1=3安培, I2=1安培, I3=-8安培I3取負(fù)值說明回路3中的實(shí)際電流與圖中顯示的電流方向相反.在方程組(*)中,如果將其系數(shù)矩陣記為R,右端列向量記為u,i =(I1,I2,I3)T,則可得到以矩陣形式表示的歐姆定律:.4 由于計(jì)算機(jī)的發(fā)展，人們開始考慮在計(jì)算機(jī)上求解線性方程組Ax=b的近似解，并用某種極限過程去逐漸逼近精確解。如初等函數(shù)及其他特殊函數(shù)，當(dāng)函數(shù)只在有限的點(diǎn)集上給定函數(shù)值，要在包含該點(diǎn)集的區(qū)間上用公式寫出函數(shù)簡單表達(dá)式，這些都涉及在去見a,b上用簡單的函數(shù)逼近一直復(fù)雜函數(shù)的問題，這也就是函數(shù)逼近問題。下面通過實(shí)例說明問題。在提拉單晶硅的過程中，假定晶體導(dǎo)電，熔化物和晶體的平衡方程式是：該模型還包括坩堝的傳導(dǎo),熔化物、晶體、坩堝、加熱元件表面和環(huán)境的熱輻射,以及融化物和周圍的環(huán)境的對流。利用格林第二定理，將面積分變成線積分Xi是邊界節(jié)點(diǎn)i出的方向矢量。是邊界線，A是橫截面積，對軸對稱的幾何解決方法，i擴(kuò)散方程式：Km是橢圓部分。（ri,zi）是節(jié)點(diǎn)i的坐標(biāo)。在數(shù)值實(shí)現(xiàn)方法、邊界被分為溫度變化和接近溫度梯度的每個(gè)值，通過插值逼近中間節(jié)點(diǎn)值:是k的函數(shù)，是坐標(biāo)，Ne是每個(gè)元素的節(jié)點(diǎn)數(shù)。方程(10) 、(12)代入(9)，形成N每個(gè)節(jié)點(diǎn)的微分/代數(shù)方程：T是節(jié)點(diǎn)處的溫度矢量，Q是節(jié)點(diǎn)的變化矢量。為了獲得一個(gè)邊界積分相當(dāng)于瞬時(shí)的域積分,一般利用插值(23)。在這種方法中,對時(shí)間的導(dǎo)數(shù)形式是:X是位置矢量，x= r z T,N是節(jié)點(diǎn)數(shù)，fkx; xk是幾何的已知函數(shù)，kt是時(shí)間的未知系數(shù)，方程（13）變?yōu)椋壕仃嘙的整數(shù)數(shù)域，注意到這矩陣的計(jì)算在一個(gè)固定的領(lǐng)域只計(jì)算一次,方程(16)應(yīng)用到N結(jié)點(diǎn),結(jié)果是:矩陣F是通過fkx; xk在節(jié)點(diǎn)處形成的，方程（18）代入方程（17），整數(shù)域里方程（9）變成：另外,為了避免整數(shù)域M在較復(fù)雜的幾何域,利用雙重定理，在這種方法中,節(jié)點(diǎn)fk x;Xk 是取代拉普拉斯節(jié)點(diǎn)的一個(gè)新的節(jié)點(diǎn), ：再將方程（21）代入（17），得到一個(gè)積分方程：如果節(jié)點(diǎn)，估測有相同的插值多項(xiàng)式k，利用溫度，矩陣M可寫成：, （23）,矩陣 ,是通過節(jié)點(diǎn)形成的。該模型的非線性介紹是通過自由表面?zhèn)鲗?dǎo)結(jié)合熱對流及熱輻射效應(yīng):例2 計(jì)算模型本文運(yùn)用有限元的方法對生長的硅單晶和晶體生長的熱場作為整體進(jìn)行研究,得到了不同條件下單晶爐內(nèi)的熱場分布。對晶體實(shí)際生產(chǎn)有一定的指導(dǎo)作用。 根據(jù)熱動力學(xué)第一定律即熱能守恒定律和熱擴(kuò)散定律,可以得到(1)式:將(1)式進(jìn)行展開得到其常見形式:(2)式是有限元計(jì)算中用到的偏微分方程。式中為物質(zhì)密度,C為熱容量,T為溫度(T=T(x, y,z, t), t為時(shí)間, L是矢量運(yùn)算符, V是熱量物質(zhì)輸運(yùn)矢量, q為熱流矢量,q為單位體積內(nèi)熱量產(chǎn)生率, D=是熱量傳導(dǎo)矩陣,Kxx, Kyy, Kzz是熱量在各坐標(biāo)軸方向的導(dǎo)熱系數(shù)。在計(jì)算中采取了近似方法5,假定: (1)固-液界面處的溫度為硅的熔點(diǎn),且界面形狀為平面;(2)物性參數(shù)如熱傳導(dǎo)系數(shù)、密度、熱容量等不隨溫度變化,都為常數(shù); (3)單晶爐內(nèi)熱場主要以熱傳導(dǎo)6和輻射7,8為主,不考慮熔體的熱對流問題。計(jì)算過程中用到的材料屬性和處理參數(shù)如表1所示。4、方程的性態(tài)與誤差分析實(shí)際上由于誤差的存在，直接法只能求得近似解。誤差分析是計(jì)算方法中一個(gè)既重要復(fù)雜的問題。 因?yàn)閹缀趺恳徊竭\(yùn)算都有誤差，而實(shí)際問題往往需要進(jìn)行千百萬次計(jì)算，所以每步運(yùn)算即分析誤差幾乎是不可能的，也是不必要的。設(shè)Axb（矩陣A非奇異）中，系數(shù)矩陣有小擾動和b，擾動后的方程組為A(A+A)(x +x) = b +b若A為n階非奇異矩陣，稱Cond(A) = A1 A為矩陣A的條件數(shù)。當(dāng)條件數(shù)Cond(A)比較大時(shí)，系數(shù)矩陣和右端向量的小擾動會引起解的相對誤差較大。如果矩陣A的條件數(shù)比較大，就說方程組是“病態(tài)”的；如果矩陣A的條件數(shù)比較小，就說方程組是“良態(tài)”的；方程組的性態(tài)是方程組的固有性質(zhì)，它與求解方法無關(guān)。在計(jì)算機(jī)上求解的方程組都是所給方程組的擾動方程（字長限制將產(chǎn)生舍入誤差）。對于良態(tài)方程組，只要求解方法穩(wěn)定，即可得到比較滿意的計(jì)算結(jié)果。但對于病態(tài)方程組，即使使用穩(wěn)定性好的算法求解也未必理想。求解病態(tài)方程組：采用高精度算術(shù)運(yùn)算 在求解方程組之前，對原方程