數值分析解線性方程組的直接方法演示文稿現在是1頁\一共有91頁\編輯于星期三優選數值分析解線性方程組的直接方法Ppt現在是2頁\一共有91頁\編輯于星期三

在自然科學和工程技術中有很多問題的解決常常歸結為解線性代數方程組．如三次樣條函數問題，用最小二乘法求實驗數據的曲線擬合問題，解非線性方程組問題，用差分法或者有限元方法解常微分方程、偏微分方程的邊值問題等都導致求解線性代數方程組，而這些方程組的系數矩陣大致分為兩種，一種是低階稠密矩陣，另一種是大型稀疏矩陣。

關于線性方程組的數值解法一般有兩類：1.直接法2.迭代法§5.1引言與預備知識5.1.1引言現在是3頁\一共有91頁\編輯于星期三

本章討論n元線性方程組

（5.1）的直接解法。方程組(5.1)的矩陣形式為

Ax=b其中

若矩陣A非奇異,即det(A)≠0,則方程組(2.1)有唯一解。現在是4頁\一共有91頁\編輯于星期三

所謂直接解法是指，若不考慮計算過程中的舍入誤差，經過有限次算術運算就能求出線性方程組的精確解的方法。但由于實際計算中舍入誤差的存在，用直接解法一般也只能求出方程組的近似解。Cramer法則是一種不實用的直接法，本章將介紹幾種實用的直接法。現在是5頁\一共有91頁\編輯于星期三5.1.2預備知識M行n列矩陣.n維列向量.現在是6頁\一共有91頁\編輯于星期三矩陣的基本運算:(1)矩陣的加法(7)矩陣的行列式行列式性質:(a)det(AB)=det(A)det(B)(6)非奇異矩陣(5)單位矩陣(4)轉置矩陣(3)矩陣與矩陣的乘法(2)矩陣與標量的乘法現在是7頁\一共有91頁\編輯于星期三矩陣特征值與譜半徑定義1設若存在一個數λ(實數或復數)和非零向量使(1.1)則稱λ為A的特征值,x為A對應λ的特征向量,A的全體特征值稱為A的譜，記作稱為A的譜半徑.(1.2)現在是8頁\一共有91頁\編輯于星期三由式(1.1)知,λ

可使齊次方程組有非零解,故系數行列式記稱為矩陣Ａ的特征多項式，方程(1.3)稱為Ａ的特征方程．(1.3)現在是9頁\一共有91頁\編輯于星期三在復數域中有n個根故由行列式(1.3)展開可知：的n個特征值故是它的特征方程(1.3)的幾個根,并有(1.4)(1.5)A的跡.現在是10頁\一共有91頁\編輯于星期三A的特征值λ和特征向量x還有以下性質:(1)AT與A有相同的特征值λ及相同的特征向量x.(2)若A非奇異，則A-1的特征值為λ-1,特征向量為x.

(3)相似矩陣B=S-1AS有相同的特征多項式.現在是11頁\一共有91頁\編輯于星期三例1求的特征值及譜半徑.解:A的特征方程為故A的特征值為A的譜半徑為現在是12頁\一共有91頁\編輯于星期三5.1.4特殊矩陣現在是13頁\一共有91頁\編輯于星期三現在是14頁\一共有91頁\編輯于星期三定理1.現在是15頁\一共有91頁\編輯于星期三定理2.現在是16頁\一共有91頁\編輯于星期三定理3.定理4(Jordan標準型)設A為n階矩陣,則存在一個非奇異矩陣P使得現在是17頁\一共有91頁\編輯于星期三其中:(1)當A的若當標準型中所有若當塊Ji均為一階時,此標準型變成對角矩陣.返回主頁現在是18頁\一共有91頁\編輯于星期三求解高斯消去法(逐次消去法)及消去法和矩陣三角分解之間的關系：§5.2高斯消去法現在是19頁\一共有91頁\編輯于星期三例2

用消去法解方程組解第1步.將方程(2.2)乘上-2加到方程(2.4)上,消去未知數x1,得到(2.2)(2.3)(2.4)第2步.將方程(2.3)加到方程(2.5)上去,消去方程(2.5)中的x2,得到與原方程組等價的三角形方程組解為:首先舉一個簡單的例子來說明消去法的基本思想.現在是20頁\一共有91頁\編輯于星期三上述過程相當于思路首先將A化為上三角陣/*upper-triangularmatrix*/，再回代求解/*backwardsubstitution*/。=現在是21頁\一共有91頁\編輯于星期三消元記Step1：設，計算因子將增廣矩陣/*augmentedmatrix*/第i行mi1

第1行，得到其中Stepk：設，計算因子且計算現在是22頁\一共有91頁\編輯于星期三回代注意1：只要A

非奇異，即A1

存在，則可通過逐次消元及行交換，將方程組化為三角形方程組，求出唯一解。共進行？步n

1現在是23頁\一共有91頁\編輯于星期三注意2:設Ax=b,其中A為非奇異矩陣,如果由于A為非奇異矩陣,所以A的第一列一定有元素不等于零.例如現在是24頁\一共有91頁\編輯于星期三定理5

設Ax=b,其中(1)如果則可通過高斯消去法將Ax=b約化為等價的三角形線性方程組(2.10),且計算公式為:①消元計算(k=1,2,…,n-1)②回代計算現在是25頁\一共有91頁\編輯于星期三(2)如果A為非奇異矩陣,則可通過高斯消去法(及交換兩行的初等變換)將方程組Ax=b約化為方程組(2.10).定理6

約化的主元素aii(i)

≠0(i=1,2,…,k)的充要條件是矩陣A的所有順序主子式

/*determinantofleadingprincipalsubmatrices*/Di≠0(i=1,2,…,k).即(2.12)推論如果A的順序主子式Dk≠0(k=1,2,…,n-1),則現在是26頁\一共有91頁\編輯于星期三§5.2.2三角分解法

/*MatrixFactorization*/高斯消元法的矩陣形式

/*MatrixFormofG.E.*/：Step1:記L1=，則Stepn

1:其中

Lk=現在是27頁\一共有91頁\編輯于星期三記為L單位下三角陣/*unitarylower-triangularmatrix*/記

U=現在是28頁\一共有91頁\編輯于星期三定理7

若A的所有順序主子式/*determinantofleadingprincipalsubmatrices*/

均不為0，則A

的

LU

分解唯一（其中L

為單位下三角陣）。證明：由§1中定理可知，LU分解存在。下面證明唯一性。若不唯一，則可設A=L1U1=L2U2

，推出Upper-triangularLower-triangularWithdiagonalentries1注：L

為一般下三角陣而U

為單位上三角陣的分解稱為Crout分解。實際上只要考慮A*的LU

分解，即

，則即是A的Crout分解。現在是29頁\一共有91頁\編輯于星期三例3

對于例2,系數矩陣由高斯消去法,返回主頁現在是30頁\一共有91頁\編輯于星期三5.3.1列主元消去法:在順序消元過程中,只要即可進行計算,但如果很小,則將導致舍入誤差增長,使解的誤差很大.例4用Gauss消去法求解方程組§5.3高斯主元素消去法現在是31頁\一共有91頁\編輯于星期三解:因故方程有唯一解,且精確解為若用Gauss消去法取四位有效數字計算,可得解比較,誤差很大,若將兩個方程互換為用Gauss消去法取四位有效數字計算,可得解現在是32頁\一共有91頁\編輯于星期三

本例表明通過行交換可避免舍入誤差增長,這就是列主元消去法的基本思想.其計算步驟如下:第1步，在中的第1列選主元，即行為主元,若將

的第i1行與第1行互換，再按消元公式計算得到假定上述過程已進行(k-1)步，得到第k步，在中的第k列選主元，即若則在中將ik行與第k行互換，再按消元公式計算得到對k=1,2,…,n-1,重復以上過程則得如果某個k出現主元方程沒有唯一解或嚴重病態，否則可由(3.2.4)求得解.則表明detA=0，現在是33頁\一共有91頁\編輯于星期三它也表明當Ａ非奇異時，存在排列矩陣Ｐ(若干初等排列矩陣的乘積),使PA=LU,其中L為單位下三角矩陣,其元素|lij|<=1,U為上三角矩陣.上述每步行交換后再消元相當于其中是指標為k的初等下三角矩陣,為初等排列矩陣時,表示不換行,經過(n-1)步換行與消元,A化為上三角矩陣.即:現在是34頁\一共有91頁\編輯于星期三解:例5用列主元消去法解Ａx=b,其中消元消元現在是35頁\一共有91頁\編輯于星期三消元結束．由回代公式求得解此例的精確解為可見結果精度較高．若不選列主元Ｇauss消去法，求得解，誤差較大．除列主元消去法外，還有一種消去法，是在Ａ的所有元素aij中選主元，稱為全主元消去法．因計算量較大且應用列主元已能滿足實際要求，故不再討論．目前很多數學軟件庫都有列主元消去法，可直接調用．現在是36頁\一共有91頁\編輯于星期三注:為了減少計算的舍入誤差,使用消去法通常都要選主元.目前最常用的是列主元消去法,也就是每步消元之前選主元,當A=(aij)第一步選A中第1列的主元,即max|ai1|=ai1.然后將i1行與第１行互換，再進行消元，以后每步消元做法類似，先選主元，再消元．現在是37頁\一共有91頁\編輯于星期三5.3.2高斯若當消去法消去對角線上方和下方的元素.假設已經完成k-1步,得到與方程Ax=b等價的方程組返回主頁現在是38頁\一共有91頁\編輯于星期三

高斯消去法有很多變形,有的是高斯消去法的改進、改寫，有的是用于某一類特殊性質矩陣的高斯消去法的簡化。5.3.1直接三角分解法.

將高斯消去法改寫為緊湊形式,可以直接從A的元素得到計算L,U元素的遞推公式,而不需要任何中間步驟,這就是所謂直接三角分解法.

一旦實現A的LU分解,那么求解Ax=b的問題就等價于求解兩個三角形方程組

(1)Ly=b,求y;(2)Ux=y,求x.§５.4矩陣三角分解法

/*MatrixFactorization*/現在是39頁\一共有91頁\編輯于星期三１.不選主元的三角分解法設A為非奇異矩陣,且有分解式A=LU,其中L為單位下三角,U為上三角即L,U元素可以由n步直接計算定出,其中第r步定出U的第r行和L的第r列元素.由上式有:現在是40頁\一共有91頁\編輯于星期三故同樣有:

設已經定出U的第1行到第r-1行元素與L的第1列到第r-1列元素.利用矩陣乘法(注意當r<k時,lrk=0),有得現在是41頁\一共有91頁\編輯于星期三通過比較法直接導出L和

U的計算公式。思路現在是42頁\一共有91頁\編輯于星期三固定r:對i=r,r+1,…,n

有lii=1a將r

，i

對換，對r=i,i+1,…,n有b現在是43頁\一共有91頁\編輯于星期三結論:用直接三角分解法解Ax=b(要求A的所有順序主子式都不等于0)的計算公式如下.Step1:u1i=a1i;li1=ai1/u11;(i=2,…,n)計算U的第r行,L的第r列元素(r=2,3,…，n）Step2:求解Ly=b,Ux=y

的計算公式：Step3:Step4:Step5:現在是44頁\一共有91頁\編輯于星期三例5

用直接三角分解法解解用分解公式計算得求解現在是45頁\一共有91頁\編輯于星期三2.選主元的三角分解法當urr=0時計算中斷，或者當urr絕對值很小時，按分解公式計算可能引起舍入誤差的積累。但如果A非奇異，可以通過交換A的行實現矩陣A的LU分解，因此可采用與列主元消去法類似的方法，將直接三角分解法修改為（部分）選主元的三角分解法。設第r-1步分解已完成，這時有現在是46頁\一共有91頁\編輯于星期三第r步分解需用到（3.2）及（3.3）式，為避免用小的數urr做除數，引進量于是有：取交換A的r行與ir行元素，將調到(r,r)位置(將（i，j）位置的新元素仍記為lij

及

aii),于是有|lir|<=1(i=r+1,…，n）.由此再進行第r步分解計算。現在是47頁\一共有91頁\編輯于星期三5.3.2平方根法當A對稱正定時，A的順序主子式故由定理知，A=LU的分解存在且唯一，其中L為單位下三角為了A利用對稱性其中D為對角陣,U0為單位上三角陣,于是又代入到上式,就得到對稱矩陣A的分解式矩陣，U為上三角矩陣，且現在是48頁\一共有91頁\編輯于星期三定理9(對稱陣的三角分解定理)設A為n階對稱陣,且A的順序主子式則A可唯一分解為，其中L為單位下三角矩陣，.

D為對角矩陣定理10(對稱正定矩陣的三角分解或Cholesky分解)

如果A為n階對稱正定矩陣,則存在一個實的非奇異下三角陣L使A=LLT,當限定L的對角元素為正時,這種分解是唯一的.將求方程組的解轉化為求方程的解LLTx=b.令

，求得方程的解由根據矩陣乘法，由現在是49頁\一共有91頁\編輯于星期三得

i=j有

當i＞j,得現在是50頁\一共有91頁\編輯于星期三例6用平方根法求以下方程組的解.

解先驗證系數矩陣A對稱正定，對稱顯然，

故A對稱正定，可用Cholesky分解計算，求得

求解

得再求得現在是51頁\一共有91頁\編輯于星期三5.3.3追趕法解三對角方程組

/*CroutReductionforTridiagonalLinearSystem*/在一些實際問題中,例如解常微分方程邊值問題,解熱傳導方程以及船體數學放樣中建立三次樣條函數等,都會要求解系數矩陣為對角占優的三對角線方程組.簡記為Ax=f.其中,當|i-j|>1時,aij=0,且：現在是52頁\一共有91頁\編輯于星期三Step1:對A作Crout分解直接比較等式兩邊的元素，可得到計算公式。Ｌ為下三角，Ｕ為單位上三角現在是53頁\一共有91頁\編輯于星期三注意當j=1時有對j=2，3，…，n求得L的元素，這就是A的Cholesky分解，然后再解兩個三角方程組，得這就是對稱正定方程組的平方根法

另外，由于

故有這表明分解過程中矩陣L中元素因此平方根法計算是數值穩定的.

的數量級不增長，現在是54頁\一共有91頁\編輯于星期三Step2:追——即解Step3:趕——即解：求解Ａx＝f等價于Step１:計算的遞推公式將計算系數為追的過程將計算方程組的解為趕的過程現在是55頁\一共有91頁\編輯于星期三定理：設有三對角線方程組Ａx＝f，其中Ａ滿足對角占優的條件，則Ａ為非奇異矩陣且追趕法計算公式中的滿足：現在是56頁\一共有91頁\編輯于星期三定理

若A

為對角占優

/*diagonallydominant*/的三對角陣，且滿足，則追趕法可解以A

為系數矩陣的方程組。注：

如果A是嚴格對角占優陣，則不要求三對角線上的所有元素非零。

根據不等式可知：分解過程中，矩陣元素不會過分增大，算法保證穩定。

運算量為O(6n)。返回主頁現在是57頁\一共有91頁\編輯于星期三5.5.1內積與向量范數為了研究方程組Ax=b解的誤差和迭代法收斂性，需對向量及矩陣的"大小"引進一種度量，就要定義范數，它是向量"長度"概念的直接推廣，通常用表示n維實向量空間，表示n維復向量空間.定義2

設將實數稱為向量x,y的數量積.非負實數稱為向量x的歐氏范數或2-范數.§5.5向量和矩陣范數現在是58頁\一共有91頁\編輯于星期三定理12設則內積有以下性質：(1)(2)(3)(4)(5)(柯西-施瓦茨不等式)等式當且僅當x與y線性相關時成立;(6)三角不等式現在是59頁\一共有91頁\編輯于星期三定義3(向量范數)如果向量的某個實值函數滿足條件:則稱現在是60頁\一共有91頁\編輯于星期三對于由內積性質可知它滿足定義2的三個條件，故它是一種向量范數.此外還有以下幾種常用的向量范數.容易驗證均滿足定義2的三個條件.更一般的還可定義但只有p=1,2,∞時的三種范數是常用的向量范數.例如給定則可求出現在是61頁\一共有91頁\編輯于星期三定理14設是則N(x)是向量x的分量上任一種向量范數，的連續函數.定理15(向量范數的等價性)設是上任意兩種向量范數，則存在常數使現在是62頁\一共有91頁\編輯于星期三5.5.2矩陣范數

矩陣可看成n×n維向量，如果直接將向量的2-范數用于矩陣A，則可定義稱為矩陣A的Frobenius范數，簡稱F-范數.它顯然滿足向量范數的三條性質，但由于矩陣還有乘法運算，因此矩陣范數的定義中應增加新條件.現在是63頁\一共有91頁\編輯于星期三定義4

如果的某個非負實函數N(A)，記作‖A‖，滿足條件：則稱現在是64頁\一共有91頁\編輯于星期三顯然滿足定義中的四個條件,(3)，(4)兩條均可由Cauchy-Schwarz不等式證明，故是一種矩陣范數.

除矩陣自身的運算外，在解方程中矩陣乘向量的運算即Ax，也是必不可少的.因此要求所引進的范數應滿足條件：上式稱為相容性條件.為使引進的矩陣范數滿足條件(4.5)，我們給出以下定義.(4.5)現在是65頁\一共有91頁\編輯于星期三定義6(矩陣的算子范數)設當給定向量范數時可定義稱為矩陣的算子范數或從屬范數.(4.6)定理17

設上的一種向量范數，則由(4.6)定義的是一種矩陣范數，且滿足相容性條件現在是66頁\一共有91頁\編輯于星期三證明因中有界閉集上的連續函數，故在D上有最大值，即使而對故所以從而當成立，而x=0時顯然也成立.現在是67頁\一共有91頁\編輯于星期三定理17

設則這里為矩陣的譜半徑.現在是68頁\一共有91頁\編輯于星期三例7

已知解從定理可以看出，計算較容易，而計算

時因為要求的特征值，所以較為困難.但當A對稱時，有現在是69頁\一共有91頁\編輯于星期三定理19定理18

對任何為任一種從屬范數則反之，對任意ε＞0，至少存在一種從屬范數使證明:設為A的特征值，則由得現在是70頁\一共有91頁\編輯于星期三非奇異，且證明用反證法.假定(I+B)奇異，則齊次方程有非零解而與‖B‖＜1的假設矛盾，故(I+B)非奇異.

又得取范數得定理20設

返回主頁現在是71頁\一共有91頁\編輯于星期三矩陣條件數與擾動方程組誤差界

在解方程組Ax=b時，由于各種原因，A或b往往有誤差，從而使得解也產生誤差.例8

方程組

的準確解為

,當A與b有微小變化時，如變為方程則準確解為

它表明A,b的微小擾動引起方程解x的很大變化，這就是病態方程.§5.6誤差分析與病態方程組現在是72頁\一共有91頁\編輯于星期三定義7

求解線性方程組Ax=b時，若A或b有微小擾動

解x的誤差

很大,，則稱此方程組為病態方程組，相應的系數矩陣A稱為病態矩陣.

反之，若此時

很小,，則稱此方程組為良態方程組，相應的系數矩陣A稱為良態矩陣.

注意方程組是否病態與用什么數值方法無關，它是由方程自身性質決定的.

在例8中因為行列式

因此出現病態.但有時A從表面上看性質很好，也可能是病態的.現在是73頁\一共有91頁\編輯于星期三那么如何判斷A是否病態?先給出如下定義.例9

方程組Ax=b表示為它的準確解

A對稱正定且

表面看性質"較好"，但若對右端b作微小變化，如方程改為

則解變為

這里b的相對誤差大約只有

但解的相對誤差卻很大，故A也是病態矩陣.現在是74頁\一共有91頁\編輯于星期三定義8

設

非奇異，‖·‖v為矩陣的任一種從屬范數，則

稱為矩陣A的條件數.

從定義看到矩陣條件數依賴于范數的選取，如范數為2-范數，

則記為

同理有

等等.現在是75頁\一共有91頁\編輯于星期三條件數有以下性質：

(1)(2)(3)U為正交矩陣，則

(4)若

與為A的按模最大與最小特征值，則若A對稱，則

現在是76頁\一共有91頁\編輯于星期三下面給出擾動方程組解的誤差分析.先考察b有擾動

則擾動方程為由于Ax=b,故得

于是再由Ax=b,有

即故得現在是77頁\一共有91頁\編輯于星期三下面再研究方程Ax=b，當A有擾動

時，其解

的誤差分析.

此時擾動方程為

因Ax=b,故有

因存在,若假定

則由定理20可知

非奇異,并有:(5.6)由(5.6)可得

因此(5.7)現在是78頁\一共有91頁\編輯于星期三定理22

設A為非奇異矩陣,Ax=b≠0,且如果則(5.7)式成立.從(5.7)看到，當A的條件數Cond(A)很大時，解的相對誤差

也很大，故方程組為病態.在例9中

而于是條件數很大，故方程是嚴重病態的.

現在是79頁\一共有91頁\編輯于星期三例10

Hibert矩陣是一個著名的病態矩陣，記作

它是一個對稱正定矩陣，當n≥3時它是病態矩陣.例如

故另外還有

等等.因此Hn是嚴重的病態矩陣，且n越大

Cond(Hn)越大.現在是80頁\一共有91頁\編輯于星期三例11

在例10的方程組中可算出A的特征值

故例中實際相對誤差是而根據(5.6)的誤差估計為這與實際相差不大，即相對誤差放大了將近3000倍.故方程為病態方程組.現在是81頁\一共有91頁\編輯于星期三定理23(事后誤差估計)設方程組

，則若實際求得解為證明記剩余則它表明如果方程組病態，即使剩余‖r‖很小，解的相對誤差仍可能很大.現在是82頁\一共有91頁\編輯于星期三5.6.2病態方程組的解法

如果A的條件數Cond(A)>>1，則Ax=b為病態方程，但計算Cond(A)時需要求A-1,計算量很大，相當于解方程組，在實際中常可通過求解過程直觀地判斷方程組的病態性質，如果解方程時出現下述情況之一，則可能是"病態"方程組.(1)在列主元消去法中出現小主元；(2)在計算過程中行或列幾乎線性相關或三角分解中對角元出現近似零的元素；(3)矩陣A的元素數量級相差很大且無規律；(4)剩余很小，而解很大，又達不到精度要求.現在是83頁\一共有91頁\編輯于星期三

(1)采用高精度運算，減輕病態影響，例如用雙倍字長運算.對病態方程組求解可采用以下措施：(2)用預處理方法改善A的條件數，即選擇非奇矩陣與Ax=b等價,而(3)平衡方法，當A中元素的數量級相差很大，可采用行均衡或列均衡的方法改善A的條件數.設非奇異，計算于是求Ax=b等價于求的條件數可得到改善，這就是行均衡法.現在是84頁\一共有91頁\編輯于星期三例12

給定方程組Ax=b為A的條件數若用行均衡法可取則平衡后的方程用三位有效數字的列主元消去法求解得現在是85頁\一共有91頁\編輯于星期三functionx=threedia(a,b,c,f)N=length(f);x=zeros(1,N);y=zeros(1,N);beta=zeros(1,N);gramma=zeros(1,N);beta(1)=b(1);fori=1:N-1gramma(i)=c(i)/beta(i);beta(i+1)=b(i+1)-a(i+1)*gramma(i);end%追的過程y(1)=f(1)/beta(1);fori=2:Ny(i)=(f(i)-a(i)*y(i-1))/beta(i);end%趕的過程x(N)=y(N);fori=N-1:-1:1x(i)=y(i)-gramma(i)*x(i+1);end

a=[0,-1,-1,-3];>>b=[2,3,2,5];>>c=[-1,-2,-1,0];>>f=[6,1,0,1]';>>x=threedia(a,b,c,f)追趕法求解三對角方程組現在是86頁\一共有91頁\編輯于星期三Cholesky方法：

A=[4,-2,4;-2,17,10;4,10,9];

b=[8.7,13.7,-0.7]';

[x,L,D]=Chol_decompose(A,b)L=1.000000-0.50001.000001.00000.75001.0000D=416-4x=-5.1457-3.17275.7344%用Cholesky分解求解%A是對稱矩陣%L是單位下三角陣%D是對角陣%對稱陣A進行三角分解：A=LDL'現在是87頁\一共有91頁\編輯于星期三function[x,L,D]=Chol_decompose(A,b)N=length(A);L=zeros(N,N);D=zeros(1,N);fori=1:NL(i,i)=1;endD(1)=A(1,1);fori=2:Nforj=1:i-1ifj==1L(i,j)=A(i,j)/D(j);elsesum1=0;fork=1:j-1sum1=sum1+L(i,k)*D(k)*L(j,k);endL(i,j)=(A(i,j)-sum1)/D(j);endendsum2=0;fork=1:i-1sum2=sum2+L(i,k)^2*D(k);endD(i)=A(i,i)-sum2;end%分別求解線性方程組Ly=b;L'x=y/Dy=zeros(1,N);y(1)=b(1);fori=2:Nsumi=0;fork=1:i-1sumi=sumi+L(i,k)*y(k);endy(i)=b(i)-sumi;endx=zeros(1,N);x(N)=y(N)/D(N);fori=N-1:-1:1sumi=0;fork=i+1:Nsumi=sumi+L(k,i)*x(k);endx(i)=y(i)/D(i)-sumi;end現在是88頁\一共有91頁\編輯于星期三用Dollittle三角分解法求解方程組：>>A=[0.001,2,3;-1,3.712,4.623;-2,1.072,5.643];>>b=[1,2,3]';>>[x,L,U]=lu_decompose(A,b)x=-0.4904-0.05100.3675L=1.0e+003*0.001000-1.00000.00100-2.00000.00200.0010U=1.0e+003*0.00000.00200.003002.00373.0046000.0059現在是89頁\一共有91頁\編輯于星期三function[x,L,U]=lu_decompose(A,b)%用Dollittle三角分解