解線性方程組的直接方法2023/7/201第1頁，課件共155頁，創作于2023年2月這一章討論線性方程組引言與預備知識的數值解法.2023/7/202第2頁，課件共155頁，創作于2023年2月在自然科學和工程技術中,很多問題歸結為解線性方程組.例如電學中的網絡問題，船體數學放樣中建立三次樣條函數問題，用最小二乘法求實驗數據的曲線擬合問題，解非線性方程組問題，用差分法或者有限元方法解常微分方程、偏微分方程邊值問題等都導致求解線性方程組，而這些方程組的系數矩陣大致分為兩種，一種是低階稠密矩陣(例如，階數不超過150).另一種是大型稀疏矩陣(即矩陣階數高且零元素較多).一、引言2023/7/203第3頁，課件共155頁，創作于2023年2月有的問題的數學模型中雖不直接表現為含線性方程組,但它的數值解法中將問題“離散化”或“線性化”為線性方程組.因此線性方程組的求解是數值分析課程中最基本的內容之一.

關于線性方程組的解法一般有兩大類：2023/7/204第4頁，課件共155頁，創作于2023年2月但實際計算中由于舍入誤差的存在和影響，這種方法也只能求得線性方程組的近似解.本章將闡述這類算法中最基本的和具有代表性的算法就是高斯消去法,以及它的某些變形和應用.這類方法是解低階稠密矩陣方程組及某些大型稀疏矩陣方程組(例如，大型帶狀方程組)的有效方法.

1.直接法經過有限次的算術運算,可以求得方程組的精確解(假定計算過程沒有舍入誤差).如線性代數課程中提到的克萊姆算法就是一種直接法.但該法對高階方程組計算量太大,不是一種實用的算法.2023/7/205第5頁，課件共155頁，創作于2023年2月就是用某種極限過程去逐步逼近方程組精確解的方法.迭代法具有計算機的存儲單元較少、程序設計簡單、原始系數矩陣在計算過程中始終不變等優點，但存在收斂條件和收斂速度問題.迭代法是解大型稀疏矩陣方程組(尤其是由微分方程離散后得到的大型方程組)的重要方法(見第6章).為了討論線性方程組數值解法，需復習一些基本的矩陣代數知識.

2.迭代法2023/7/206第6頁，課件共155頁，創作于2023年2月二、向量和矩陣

基本概念：用Rm×n表示全部m×n實矩陣的向量空間；用Cm×n表示全部m×n復矩陣的向量空間.(由數排成的矩陣表，稱為m行n列矩陣).2023/7/207第7頁，課件共155頁，創作于2023年2月其中aj為A的第j列的m維列向量.同理(m維列向量).其中biT為A的第i行的n維行向量.2023/7/208第8頁，課件共155頁，創作于2023年2月

矩陣的基本運算：

(1)矩陣加法

(2)矩陣與標量的乘法

(3)矩陣與矩陣的乘法

(4)轉置矩陣2023/7/209第9頁，課件共155頁，創作于2023年2月(5)單位矩陣其中

(6)非奇異矩陣則稱B是A的逆矩陣，記為A-1，且(A-1)T=(AT)-1.如果A-1存在，則A稱為非奇異矩陣.如果A、B均為非奇異矩陣，則有(AB)-1=B-1A-1.2023/7/2010第10頁，課件共155頁，創作于2023年2月

(7)矩陣的行列式設A∈Rn×n，則A的行列式可按任一行(列)展開,其中Aij為aij的代數余子式，Aij=(-1)i+jMij，為Mij元素aij的余子式.

行列式性質：2023/7/2011第11頁，課件共155頁，創作于2023年2月三、矩陣的特征值與譜半徑

定義1設A=(aij)Rn×n，若存在數(實數或復數)和非零向量x=(x0,x1,…,xn)TRn，使Ax=x,(0.1)則稱為A的特征值,x為A對應的特征向量,A的全體特征值稱為A的譜,記作σ(A),即σ(A)={1,2,…,n}.記稱為A的譜半徑.由(0.1)式可知可使齊次線性方程組(I-A)x=0有非零解，故系數行列式det(I-A)=0，記2023/7/2012第12頁，課件共155頁，創作于2023年2月p()稱為矩陣A的特征多項式，方程(0.3)稱為A的特征方程.因為n次代數方程p()在復數域中有n個根1,2,…,n，故2023/7/2013第13頁，課件共155頁，創作于2023年2月由(1.3)式中的行列式展開可得故A=(aij)Rn×n的n個特征值1,2,…,n是它的特征方程(1.3)的n個根，并有稱trA為A的跡.及2023/7/2014第14頁，課件共155頁，創作于2023年2月

A的特征值和特征向量x還有以下性質：⑴AT與A有相同的特征值及特征向量x.⑵若A非奇異,則A-1的特征值為-1,特征向量為x.⑶相似矩陣B=S-1AS有相同的特征多項式.

例1求矩陣的特征值及譜半徑

解矩陣A的特征方程為得A的特征值為1=2=2,3=7，譜半徑(A)=7.2023/7/2015第15頁，課件共155頁，創作于2023年2月特殊矩陣設A=(aij)∈Rn×n.

(1)對角矩陣如果當i≠j時，aij=0.

(2)三對角矩陣如果當|i-j|>1時，aij=0.

(3)上三角矩陣如果當i>j時，aij=0.

(4)上海森伯格(Hessenberg)矩陣如果當i>j+1時，aij=0.

(5)對稱矩陣如果AT=A.

(6)埃爾米特(Hermit)矩陣設A∈Cn×n，如果當AH=A(AH是A的共軛轉置矩陣).2023/7/2016第16頁，課件共155頁，創作于2023年2月

(7)對稱正定矩陣如果AT=A且對任意非零向量，x∈Rn，(Ax,x)=xTAx>0.

(8)正交矩陣如果A-1=AT.(9)酉矩陣設A∈Cn×n，如果A-1=AH.(10)初等置換陣由單位矩陣I交換第i行與第j行(或第i列與第j列)，得到的矩陣記為Iij，且.

IijA=B(為交換A第i行與第j行得到的矩陣)；

AIij=C(為交換A第i列與第j列得到的矩陣).

(11)置換陣由初等置換陣的乘積得到的矩陣.2023/7/2017第17頁，課件共155頁，創作于2023年2月

定理1設A∈Rn×n，則下述命題等價：(1)對任何b∈Rn，方程組Ax=b有唯一解.(2)齊次方程組Ax=0只有唯一解零解x=0.(3)det(A)≠0.(4)A-1存在.(5)A的秩rank(A)=n.2023/7/2018第18頁，課件共155頁，創作于2023年2月

定理2設A∈Rn×n為對稱正定矩陣，則(1)A為非奇異矩陣，且A-1亦是對稱正定矩陣.(2)設Ak為A的順序主子陣，則Ak(k=1,2,,n)亦是對稱正定矩陣，其中(3)A的特征值i>0

(i=1,2,,n).(4)A的順序主子式都大于零，即Ak的行列式det(Ak)>0(k=1,2,,n).2023/7/2019第19頁，課件共155頁，創作于2023年2月

定理3設A∈Rn×n為對稱矩陣，如果det(Ak)>0

(k=1,2,,n)，或A的特征值i>0

(i=1,2,,n)，則A為對稱正定矩陣.有重特征值的矩陣不一定相似于對角矩陣，那么一般n階矩陣A在相似變換下能簡化到什么形狀.

定理4(Jordan標準型)設A為n階矩陣，則存在一個非奇異n階矩陣P使得2023/7/2020第20頁，課件共155頁，創作于2023年2月為若當(Jordan)塊.其中2023/7/2021第21頁，課件共155頁，創作于2023年2月(1)當A的若當標準型中所有若當塊Ji均為一階時，此標準型變為對角矩陣.(2)如果A的特征值各不相同，則其若當標準型必為對角矩陣diag(1,

2,,n).2023/7/2022第22頁，課件共155頁，創作于2023年2月6.1高斯消去法

本節介紹高斯消去法(逐次消去法)及消去法和矩陣三角分解之間的關系.雖然高斯消去法是一種古老的求解線性方程組的方法(早在公元前250年我國就掌握了解方程組的消去法)，但由它改進、變形得到的選主元素消去法、三角分解法仍然是目前計算機上常用的有效方法.我們在中學學過消去法,高斯消去法就是它的標準化的、適合在計算機上自動計算的一種方法.2023/7/2023第23頁，課件共155頁，創作于2023年2月6.1.1高斯消去法設有線性方程組或寫成矩陣形式簡記為Ax=b.2023/7/2024第24頁，課件共155頁，創作于2023年2月解為當m=n時，對三角形方程組的解非常容易的.如:上三角矩陣所對應的線性方程組2023/7/2025第25頁，課件共155頁，創作于2023年2月下三角矩陣所對應的線性方程組計算量（乘除法的主要部分）都為n2/2.解為因此，我們將一般的線性方程組化成等價的三角形方程組來求解.2023/7/2026第26頁，課件共155頁，創作于2023年2月首先舉一個簡單的例子來說明消去法的基本思想.

例2用消去法解線性方程組

解第1步，將方程(1.2)乘上-2加到方程(1.4)上去，消去(1.4)中的未知數x1，得到第2步，將方程(1.3)加到方程(1.5)上去，消去(1.5)中的未知數x2，得到與原方程組等價的三角形方程組2023/7/2027第27頁，課件共155頁，創作于2023年2月顯然，方程組是(1.6)是容易求解的，解為上述過程相當于對方程的增廣陣做初等行變換其中ri用表示矩陣的第i行.2023/7/2028第28頁，課件共155頁，創作于2023年2月由此看出，用消去法解方程組的基本思想是用逐次消去未知數的方法把原方程組Ax=b化為與其等價的三角形方程組，而求解三角形方程組可用回代的方法求解.換句話說，上述過程就是用初等行變換將原方程組系數矩陣化為簡單形式(上三角矩陣)，從而求解原方程組(1.1)的問題轉化為求解簡單方程組的問題.或者說，對系數矩陣A施行一些行變換(用一些簡單矩陣左乘A)將其約化為上三角矩陣.這就是高斯消去法.2023/7/2029第29頁，課件共155頁，創作于2023年2月下面討論求解一般線性方程組的高斯消去法.由2023/7/2030第30頁，課件共155頁，創作于2023年2月將(1.1)記為A(1)x=b(1)，其中2023/7/2031第31頁，課件共155頁，創作于2023年2月(1)第1步(k=1).

設a11(1)0，首先計算乘數

mi1=ai1(1)/a11(1)(i=2,3,…,m).用-mi1乘(1.1)的第一個方程，加到第i個(i=2,3,,m)方程上，消去(1.1)的從第二個方程到第m個方程中的未知數x1，得到與(1.1)等價的方程組(1.7)2023/7/2032第32頁，課件共155頁，創作于2023年2月簡記為A(2)x=b(2)，其中A(2),b(2)的元素計算公式為(2)第k次消元(k=1,2,,s,s=min(m-1,n)).

設上述第1步，

，第k-1步消元過程計算已經完成，即已計算好與(1.1)等價的方程組，簡記為A(k)x=b(k)，其中2023/7/2033第33頁，課件共155頁，創作于2023年2月

設akk(k)0，計算乘數

mik=aik(k)/akk(k)

(i=k+1,,m).用-mik乘(1.8)的第k個方程加到第i個(i=k+1,,m)方程上，消去從第k+1個方程到第m個方程中的未知數xk，得到與(1.1)等價的方程組A(k+1)x=b(k+1)，2023/7/2034第34頁，課件共155頁，創作于2023年2月其中A(k+1),b(k+1)的元素計算公式為，顯然A(k+1)中從第1行到第k行與A(k)相同.

(3)繼續上述過程，且設aii(i)0(i=1,2,,s)，直到完成第s步消元計算.最后得到與原方程組等價的簡單方程組A(s+1)x=b(s+1)，其中A(s+1)為上階梯.2023/7/2035第35頁，課件共155頁，創作于2023年2月

特別當m=n時，與原方程組等價的簡單方程組為A(n)x=b(n)，即由(1.1)約化為(1.10)的過程稱為消元過程.

如果A是非奇異矩陣，且akk(k)0(k=1,2,,n-1)，求解三角形方程組(1.10)，得到求解公式(1.10)的求解過程(1.11)稱為回代過程.2023/7/2036第36頁，課件共155頁，創作于2023年2月

注意：設Ax=b，其中A∈Rn×n為非奇異矩陣，如果a11(1)=0，由于A為非奇異矩陣，所以A的第1列一定有元素不等于零，例如al10，于是可交換兩行元素(即r1?rl)，將al1

調到(1,1)位置，然后進行消元計算，這時A(2)右下角矩陣為n-1階非奇異矩陣，繼續這過程，高斯消去法照樣可進行計算.總結上述討論即有2023/7/2037第37頁，課件共155頁，創作于2023年2月

定理5設Ax=b，其中A∈Rn×n.

(1)如果akk(k)0(k=1,2,…,n)，則可通過高斯消去法將Ax=b約化為等價的三角形方程組(1.10)，且計算公式為

①消元計算(k=1,2,,n-1)2023/7/2038第38頁，課件共155頁，創作于2023年2月

②回代計算

(2)如果A為非奇異矩陣，則可通過高斯消去法(及交換兩行的初等變換)將方程組Ax=b約化為等價的三角形方程組(1.10).2023/7/2039第39頁，課件共155頁，創作于2023年2月

算法1(高斯算法)設A∈Rm×n(m>1),s=min(m-1,n)，如果akk(k)0(k=1,2,,s)，本算法用高斯方法將A約化為上三角形矩陣，且A(k)覆蓋A，乘數mik覆蓋aik.

對于k=1,2,,s

(1)如果akk=0，則計算停止；

(2)對于i=k+1,,m

(a)aik←mik=aik/akk

(b)對于j=k+1,,n

aij←aij-mik·akj

.2023/7/2040第40頁，課件共155頁，創作于2023年2月

顯然，算法1第k步需要m-k次除法，(m-k)(n-k)次乘法，因此，本算法(從第1步到第s步消元計算總的計算量)大約需要s3/3-(m-k)s2/2+mns次乘法運算(對相當大的s).當m=n時,總共大約需要n3/3次乘法運算.

數akk(k)在高斯消去法中有著突出的作用，稱為約化的主元素(簡稱主元).2023/7/2041第41頁，課件共155頁，創作于2023年2月

算法2(回代算法)設Ux=b，其中U=(uij)∈Rn×n為非奇異上三角矩陣，本算法計算Ux=b的解.

對于i=n,n-1,,1

(1)xi←bi

(2)對于j=i+1,,n

xi←xi-uij·xj(i>n不算)

(3)xi←xi/uii

這個算法需要n(n-1)/2次乘除法運算.2023/7/2042第42頁，課件共155頁，創作于2023年2月例子解方程組解：消元回代得2023/7/2043第43頁，課件共155頁，創作于2023年2月

高斯消去法對于某些簡單的矩陣可能會失敗，例如

由此，需要對算法1進行修改，首先研究原來矩陣A在什么條件下才能保證akk(k)0(k=1,2,,s).下面的定理給出了這個條件.2023/7/2044第44頁，課件共155頁，創作于2023年2月

定理6約化的主元素aii(i)0(i=1,2,,k)的充要條件是方陣A的順序主子式Di0(i=1,2,,k)，即

證明首先利用歸納法證明充分性.顯然,當k=1時，結論成立，假設結論對k-1是成立的，對k由條件設Di0(i=1,2,,k)，于是由歸納法假設我們有aii(i)0

(i=1,2,,k-1)，可用高斯消去法將A(1)約化到A(k)，即2023/7/2045第45頁，課件共155頁，創作于2023年2月利用行列式的性質，我們有2023/7/2046第46頁，課件共155頁，創作于2023年2月

由條件有Di0(i=1,2,,k)，利用(1.13)式，則有aii(i)0(i=1,2,,k)，定理6的對k成立.

必要性，由條件aii(i)0(i=1,2,,k)，利用(1.13)式亦可推出Di0(i=1,2,,k).

定理6得證.

推論如果A的順序主子式Di0(i=1,2,,n-1),則2023/7/2047第47頁，課件共155頁，創作于2023年2月6.1.2矩陣的三角分解

下面我們借助矩陣理論進一步對消去法做些分析，從而建立高斯消去法與矩陣因式分解的關系.

設(1.1)的系數矩陣A∈Rn×n的各順序主子式均不為零，對于A施行初等行變換相當于用初等矩陣左乘A，于是對(1.1)施行第1次消元后化為(1.7)，這時A(1)化為A(2)，b(1)化為b(2)，即L1A(1)=A(2),L1b(1)=b(2)，2023/7/2048第48頁，課件共155頁，創作于2023年2月其中

一般第k步消元，A(k)化為A(k+1)，b(k)化為b(k+1)，相當于LkA(k)=A(k+1),Lkb(k)=b(k+1)，2023/7/2049第49頁，課件共155頁，創作于2023年2月其中重復以上過程，最后得到(1.14)2023/7/2050第50頁，課件共155頁，創作于2023年2月將上三角矩陣A(n)記為U，由(1.14)得到其中為單位下三角矩陣.2023/7/2051第51頁，課件共155頁，創作于2023年2月

這就是說，高斯消去法實質上產生了一個將A分解為兩個三角形矩陣相乘的因式分解，于是我們得到如下重要定理，它在解方程組的直接法中起著重要作用.

定理7(矩陣的LU分解)設A為n階矩陣，如果A的順序主子式Di0(i=1,2,,n-1)，則A可分解為一個單位下三角矩陣L和一個上三角矩陣U的乘積，即A=LU,且這種分解是唯一的.2023/7/2052第52頁，課件共155頁，創作于2023年2月

證明根據以上高斯消去法的矩陣分析，A=LU的存在性已經得到證明，現僅在A為非奇異矩陣的假定下來證明唯一性，當A為奇異矩陣的情況留作練習，設A有兩種分解為A=LU=L1U1，其中L,L1為單位下三角矩陣，U,U1為上三角矩陣.由于L-1,U1-1存在，故有L-1L1=UU1-1.上式右邊為上三角矩陣，左邊為單位下三角矩陣，從而上式兩邊都必須等于單位矩陣I，故有

L=L1，U=U1.證畢.2023/7/2053第53頁，課件共155頁，創作于2023年2月

例3對于例2，系數矩陣由高斯消去法m21=0，m31=2，m32=-1，故從而得到求矩陣行列式的計算公式為2023/7/2054第54頁，課件共155頁，創作于2023年2月6.1.3列主元消去法

由高斯消去法知道,在消元過程中可能有akk(k)＝0的情況，這時消去法將無法進行；即使主元素akk(k)0但很小時，用其作除數，會導致其它元素數量級的嚴重增長和舍入誤差的擴散，最后也使得計算解不可靠.先看一個例子.

高斯消去法也稱主元素消去法(條件detA0)

即當akk(k)=0時，高斯消元法無法進行；或

|akk(k)|<<1時，帶來舍入誤差的擴散。（小主元）2023/7/2055第55頁，課件共155頁，創作于2023年2月解

(方法1)高斯消去法(用4位有效數字)例4(小主元產生了大誤差)2023/7/2056第56頁，課件共155頁，創作于2023年2月顯然計算得到的解x1=0是一個很壞的結果，不能作為方程組的近似解.其原因是我們在消元計算時用了小主元0.00001。使得約化后的方程組元素數量級大大增長，經再舍入使得在計算(0.3),(1.2),(1.3)元素時發生了嚴重的相消情況，因此經消元后得到的三角形方程組就相當不準確了，所以產生了很大的誤差，結果是不能用的.2023/7/2057第57頁，課件共155頁，創作于2023年2月

(方法2)列主元消元法，先交換行2023/7/2058第58頁，課件共155頁，創作于2023年2月

這個例子告訴我們，在采用高斯消去法解方程組時，小主元可能產生麻煩，故應避免絕對值小的主元素akk(k)，對一般矩陣來說，最好每一步選取系數矩陣(或消元后的低階矩陣)中絕對值最大的元素作為主元素，以使高斯消去法具有較好的數值穩定性，這就是全主元素消去法，在選主元時要花費較多機器時間，目前主要使用的是列主元素消去法，下面介紹列主元素消去法，并假定線性方程組(1.1)的系數矩陣A∈Rn×n為非奇異的.2023/7/2059第59頁，課件共155頁，創作于2023年2月設方程組(2.1)的增廣矩陣為首先在A的第1列中選取絕對值最大的元素作為主元素，例如2023/7/2060第60頁，課件共155頁，創作于2023年2月然后交換B的第1行與第i1行，經第1次消元計算得(A|b)→(A(2)|b(2))重復上述過程，設已完成第k-1步的選主元素，交換兩行及消元計算，(A|b)約化為其中A(k)的元素仍記為aij，b(k)的元素仍記為bi.2023/7/2061第61頁，課件共155頁，創作于2023年2月第k步的選主元素(在A(k)右下角方陣的第1列內選)，即確定ik，使交換(A(k)|b(k))第k行與ik行的元素，再進行消元計算，最后將原方程組化為

(k=1,2,,n-1)2023/7/2062第62頁，課件共155頁，創作于2023年2月

回代求解

算法1(列主元素消去法)設Ax=b.本算法用A的具有行交換的列主元素消去法，消元結果沖掉A，乘數mij沖掉aij，計算解x沖掉常數項b，行列式存放在det中.

1.det←1

2.對k=1,2,,n-1

(1)按列選主元2023/7/2063第63頁，課件共155頁，創作于2023年2月

(2)如果aik,k=0，則計算停止(det(A)=0)

(4)消元計算對于i=k+1,,n

①aik←mik=aik/akk

②對于j=k+1,,n

aij←aij-mik·akj

(3)如果ik=k，則轉(4)

換行akj?aik,j(j=k,k+1,,n)

bk?bik

det←-det

③bi←bi-mik·bk

(5)det←akk·

det2023/7/2064第64頁，課件共155頁，創作于2023年2月

3.如果ann=0，則計算停止(det(A)=0)

4.回代求解

(1)bn←bn/ann

(2)對于i=n-1,,2,1

5.det←ann·

det

例4的方法2用的就是列主元素消去法.2023/7/2065第65頁，課件共155頁，創作于2023年2月

下面用矩陣運算來描述解(1.1)的列主元素消去法.列主元素消去法為其中Lk的元素滿足|mik|≤1

(k=1,2,,n-1)，是初等置換矩陣.

利用(1.15)得到(1.15)2023/7/2066第66頁，課件共155頁，創作于2023年2月簡記為其中下面就n=4來考察一下矩陣其中可知亦為單位下三角矩陣，其元素的絕對值不超過1.記(1.16)2023/7/2067第67頁，課件共155頁，創作于2023年2月由(1.16)得到PA=LU.其中P為排列矩陣，L為單位下三角矩陣，U為上三角矩陣.這說明對(1.1)應用列主元消去法相當于對(A|b)先進行一系列行交換后對PAx=Pb再應用高斯消去法.在實際計算中我們只能在計算過程中做行的交換.總結以上的討論我們有2023/7/2068第68頁，課件共155頁，創作于2023年2月

定理8(列主元素的三角分解定理)如果A為非奇異矩陣，則存在排列矩陣P使PA=LU，其中L為單位下三角矩陣，U為上三角矩陣.在編程過程中，L元素存放在數組A的下三角部分，U元素存放在數組A的上三角部分，由紀錄主行的整型數組Ip(n)可知P的情況.2023/7/2069第69頁，課件共155頁，創作于2023年2月消元法的應用1.求行列式高斯消元法2.求逆矩陣(用高斯-若當消去法)2023/7/2070第70頁，課件共155頁，創作于2023年2月例子分別用高斯消元法和列主元消元法求矩陣的行列式.解：高斯消元法大數“吃掉”了小數！2023/7/2071第71頁，課件共155頁，創作于2023年2月列主元消元法2023/7/2072第72頁，課件共155頁，創作于2023年2月

全主元素消去法介紹在消元計算過程中，若每次選主元素不局限在方框列中，而在整個主子矩陣中選取,便稱為全主元高斯消去法,此時增加的步驟為:

①，確定ik

,jk,若=0，給出|Ａ|=0的信息，退出計算.2023/7/2073第73頁，課件共155頁，創作于2023年2月

行交換:(kjn)

列交換:(kin)②作如下行交換和列交換值得注意的是,在全主元的消去法中，由于進行了列交換，x各分量的順序已被打亂.因此必須在每次列交換的同時，讓機器“記住”作了一次怎樣的交換，在回代解出后將x各分量換回原來相應的位置，這樣增加了程序設計的復雜性.此外作①步比較大小時，全主元消去法將耗用更多的機時.2023/7/2074第74頁，課件共155頁，創作于2023年2月但全主元消去法比列主元消去法數值穩定性更好一些.實際應用中，這兩種選主元技術都在使用.選主元素的高斯消元法是一種實用的算法，它可以應用于任意的方程組Ax=b，只要|A|≠0.2023/7/2075第75頁，課件共155頁，創作于2023年2月6.2.追趕法在一些實際問題中,例如解常微分方程邊值問題,解熱傳導方程以及船體數學放樣中建立三次樣條插值函數，都會要求解系數矩陣為對角占優的三對角線方程組Ax=f，即:(2.1)2023/7/2076第76頁，課件共155頁，創作于2023年2月其中,當|i-j|>1時,aij=0,且滿足如下的對角占優條件:①|b1|>|c1|>0;②|bi|≥|ai|+|ci|,ai,ci≠0,(i=2,3,…,n-1).③|bn|>|an|>0.我們利用矩陣的直接三角分解法來推導解三對角線方程組(2.1)的計算公式.由系數陣A的特點，可以將A分解為兩個三角矩陣的乘積，即A=LU,其中取L下三角陣,取U為單位上三角陣,這樣求解方程組Ax=f的方法稱為追趕法.設2023/7/2077第77頁，課件共155頁，創作于2023年2月其中為待定系數，比較(2.2)兩邊即得(2.2)2023/7/2078第78頁，課件共155頁，創作于2023年2月下面我們用歸納法證明(2.4)對i=1是成立的.現設(2.4)對i-1成立，求證對i亦成立.(2.3)(2.3)(2.4)2023/7/2079第79頁，課件共155頁，創作于2023年2月由歸納法假設0<|i-1|<1，又由(2.4)及A的假設條件有(2.4)2023/7/2080第80頁，課件共155頁，創作于2023年2月簡記為Ax=f,A滿足條件（1）（2）（3）小結：2023/7/2081第81頁，課件共155頁，創作于2023年2月用歸納法可以證明，滿足上述條件的三對角線性方程組的系數矩陣A非奇異，所以可以利用矩陣的直接三角分解法來推導解三對角線性方程組的計算公式，用克洛特分解法，將A分解成兩個三角陣的乘積，設A=LU：

按乘法展開

則可計算

可依次計算當，由上式可惟一確定L和U2023/7/2082第82頁，課件共155頁，創作于2023年2月例6.9用追趕法求解三對角方程組

解由Ly=f

解出y又由Ux=y解出x2023/7/2083第83頁，課件共155頁，創作于2023年2月求解Ax=f等價于求解兩個三角形方程組.①Ly=f,求y;②Ux=y,求x.從而得到解三對角線方程組的追趕法公式.(1).計算{i},{i}的遞推公式小結：2023/7/2084第84頁，課件共155頁，創作于2023年2月我們將計算系數及的過程稱為追的過程,將計算方程組的解的過程稱為趕的過程.合起來就是追趕法.追趕法公式實際上就是把高斯消去法用到求解三對角線方程組上去的結果.這時由于A特別簡單,因此使得求解的計算公式非常簡單,而且計算量僅為5n-4次乘除法,而另外增加一個方程組Ax=f2僅增加3n-2次乘除法運算,易見追趕法的計算量是比較小的.2023/7/2085第85頁，課件共155頁，創作于2023年2月總結上述討論有下面定理

定理11設有三對角線方程組Ax=f,其中A滿足條件(1),(2),(3),則A為非奇異矩陣且追趕法計算公式中的{i},{i}滿足由定理11的(1),(2)說明追趕法計算公式中不會出現中間結果數量級的巨大增長和舍入誤差的嚴重累積.2023/7/2086第86頁，課件共155頁，創作于2023年2月補充：矩陣三角分解法高斯消去法有很多變形，有的是高斯消去法的改進、改寫，有的是用于某一類特殊矩陣的高斯消去法的簡化.下面我們將介紹矩陣的直接三角分解法，解特殊方程組用的平方根法及追趕法.

定義

如果L為單位下三角陣，U為上三角陣，則稱A=LU為杜里特爾(Doolittle)分解；如果L為下三角陣，U為單位上三角陣，則稱A=LU為克勞特(Crout)分解.2023/7/2087第87頁，課件共155頁，創作于2023年2月1、直接三角分解法(LU分解)

在6.1已經通過高斯消去法得到一個將A分解為一個單位下三角矩陣L和一個上三角矩陣U的乘積，A=LU，其中并由定理7得到這種分解是唯一的.2023/7/2088第88頁，課件共155頁，創作于2023年2月

將高斯消去法改寫為緊湊形式，可以直接從矩陣A的元素得到計算L、U元素的遞推公式，而不需要任何中間步驟，這就是所謂直接三角分解法.一旦實現了矩陣A的LU分解，那么求解Ax=b的問題就等價于求兩個三角形方程組.①Ly=b，求y；②Ux=y，求x.

1.不選主元的三角分解法設A為非奇異矩陣，且有分解式A=LU，其中L為單位下三角矩陣，U為上三角矩陣，即2023/7/2089第89頁，課件共155頁，創作于2023年2月其中2023/7/2090第90頁，課件共155頁，創作于2023年2月a11a12a1k

a1n

u11u12u1k

u1n

第1框a21a22a2k

a2n

l21u22u2k

u2n第2框

ak1ak2

akk

akn

lk1lk2

ukk

ukn第k框

an1an2ankann

ln1ln2

lnk

unn第n框比較式A=LU

兩端的元素,按下圖所示順序逐框進行，先求ukj，后求lik.由第一框可得2023/7/2091第91頁，課件共155頁，創作于2023年2月得假設前k-1框元素已求出，則由2023/7/2092第92頁，課件共155頁，創作于2023年2月有了矩陣A

的LU分解計算公式,解線性方程組Ax=b

就轉化為依次解下三角方程組

Ly=b

與上三角方程組Ux=y

其計算公式如下：

Ly=b

Ux=y2023/7/2093第93頁，課件共155頁，創作于2023年2月例6.10求矩陣解用緊湊形式計算的LU(Doolittle)分解.

u11=2

u12=1u13=4

2023/7/2094第94頁，課件共155頁，創作于2023年2月所以u11=2u12=1u13=4

得行列式2023/7/2095第95頁，課件共155頁，創作于2023年2月由于在計算機實現時當uri計算好后ari就不用了，因此計算好L,U的元素后就存放在A的相應位置.例如最后在存放A的數組中得到L,U的元素.由直接三角分解計算公式，需要計算形如∑aibi的式子，可采用“雙精度累加”，以提高精度.2023/7/2096第96頁，課件共155頁，創作于2023年2月如果已經實現了A=LU的分解計算，且L,U保存在A的相應位置，則用直接三角分解法解具有相同系數的方程組Ax=(b1,b2,…,bm)是相當方便的，每解一個方程組Ax=bj僅需要增加次乘除法運算.矩陣A的分解公式(3.2),(3.3)又稱為杜里特爾(Doolittle)分解.直接分解法大約需要n3/3次乘除法，和高斯消去法計算量基本相同.2023/7/2097第97頁，課件共155頁，創作于2023年2月

2.選主元的三角分解法從直接三角分解公式了看出當urr=0時計算將中斷，或者當urr絕對值很小時，按分解公式直接計算可能引起舍入誤差的積累.但如果A非奇異，我們可通過交換A的行實現矩陣PA的LU分解，因此可采用與列主元消去法類似的方法(可以證明下述方法與列主元消去法等價)，將直接三角分解法修改為(部分)選主元的三角分解法.2023/7/2098第98頁，課件共155頁，創作于2023年2月設第r-1步分解已完成，這時有設r步分解需用到(3.2)式及(3.3)式，為了避免用小的數urr作除數，引進量2023/7/2099第99頁，課件共155頁，創作于2023年2月于是有取交換A的r行與ir行元素，將調到(r,r)位置(將(i,j)位置的新元素仍記為lij及aij

),于是有|lir|1(i=r+1,…,n).由此再進行第r步分解計算.2023/7/20100第100頁，課件共155頁，創作于2023年2月算法2(選主元的三角分解法)設Ax=b，其中A為非奇異矩陣.本算法采用選主元的三角分解法，用PA的三角分解沖掉A，用整型數組Ip(n)記錄主行，解x存放在b內.

1.對于r=1,2,…,n

(1)計算si

(2)選主元2023/7/20101第101頁，課件共155頁，創作于2023年2月上述計算過程完成后就實現了PA的LU分解，且U保存在A的上三角部分，L保存在A的下三角部分，排列矩陣P由Ip(n)最后記錄可知.(這時有|lir|1)

(3)交換A的r行與ir行元素

(4)計算U的第r行元素，L的第r列元素2023/7/20102第102頁，課件共155頁，創作于2023年2月

求解方程組Ly=Pb及Ux=y.

2.對于i=1,2,,n-1

(1)t←Ip(i)

(2)如果i=t則轉(3)

bi

←→bt

(3)(繼續循環)

3.

4.2023/7/20103第103頁，課件共155頁，創作于2023年2月

利用算法2的結果(實現PA=LU三角分解)，則可以計算A的逆矩陣A-1=U-1L-1P.

利用PA的三角分解計算A-1步驟：

(1)計算上三角矩陣的逆陣U-1；

(2)計算U-1L-1；

(3)交換U-1L-1列(利用Ip(n)最后記錄).

上述方法求A-1大約需要n3次乘除運算.2023/7/20104第104頁，課件共155頁，創作于2023年2月6.3平方根法實際問題中Ax=b，系數矩陣A大多是對稱正定矩陣，所謂平方根法，就是利用對稱正定矩陣的三角分解而得到求解對稱正定方程組的一種有效方法，目前在計算機上廣泛應用平方根法解此類方程組.

對稱矩陣的LDLT

分解法當A為對稱矩陣,且其順序主子式均不為0時,由定理7可知A有唯一的LU分解為(3.1)式.為了利用A的對稱性，將U再分解為2023/7/20105第105頁，課件共155頁，創作于2023年2月其中D為對角矩陣，U0為上三角矩陣，于是

A=LU=LDU0.(3.6)又A=AT=U0TDLT,由分解的唯一性即得U0T=L.代入(3.6)得到對稱矩陣A的分解式A=LDLT.總結上述討論有下面定理.2023/7/20106第106頁，課件共155頁，創作于2023年2月

定理9(對稱矩陣的三角分解定理)設A為n階對稱矩陣，且A的各階順序主子式均不為零，則A可以唯一分解為A=LDLT,其中L是單位下三角矩陣，D是對角矩陣.如果A為對稱正定矩陣，則A的分解式A=LDLT中D的對角元素di均為正數.事實上,由A的對稱正定性,定理6的推論成立,即2023/7/20107第107頁，課件共155頁，創作于2023年2月于是有由定理6得到其中L1=LD1/2為下三角矩陣.2023/7/20108第108頁，課件共155頁，創作于2023年2月

定理10(對稱正定矩陣的三角分解或喬雷斯基(Cholesky)分解)如果A為n階對稱正定矩陣，則存在一個實的非奇異下三角矩陣L使A=LLT，當限定L的對角元素為正時，這種分解是唯一的.下面我們用直接分解方法來確定計算L元素的遞推公式，因為2023/7/20109第109頁，課件共155頁，創作于2023年2月其中lii>0(i=1,2,,n).由矩陣乘法及ljk=0(當j<k時),得于是得到解對稱正定方程組Ax=b的平方根法計算公式對于j=1,2,,n2023/7/20110第110頁，課件共155頁，創作于2023年2月求解Ax=b，即求解兩個三角形方程組①Ly=b，求y；②LTx=y，求x.由計算公式(3.7)知所以有2023/7/20111第111頁，課件共155頁，創作于2023年2月于是上面分析說明，分析過程中元素ljk的數量級不會增長且對角元素ljj恒為正數.于是有結論不選主元素的平方根法是一個數值穩定的方法.當求出L的第j列元素時，LT的第j行元素亦算出.所以平方根約需n3/6次乘除法，大約為一般直接LU分解法計算量的一半.2023/7/20112第112頁，課件共155頁，創作于2023年2月例6.11用平方根法求解對稱正定方程組解首先對A進行Cholesky分解求解Ly=b，得y1=2,y2=3.5,y3=1.求解LTx=y，得x1=1,x2=1,x3=1.2023/7/20113第113頁，課件共155頁，創作于2023年2月例6.12平方根法求解方程組

解:因方程組系數矩陣對稱正定,設A=,即：由Ly=b解得由解得

由此例可以看出，平方根法解正定方程組的缺點是需要進行開方運算。為避免開方運算，我們改用單位三角陣作為分解陣，即把對稱正定矩陣A分解成

的形式，其中

2023/7/20114第114頁，課件共155頁，創作于2023年2月為對角陣，而

是單位下三角陣,這里分解公式為2023/7/20115第115頁，課件共155頁，創作于2023年2月據此可逐行計算

運用這種矩陣分解方法,方程組Ax=b即可歸結為求解兩個上三角方程組

和其計算公式分別為

和求解方程組的上述算法稱為改進的平方根法。這種方法總的計算量約為，即僅為高斯消去法計算量的一半。2023/7/20116第116頁，課件共155頁，創作于2023年2月總之，用平方根法解對稱正定方程組時，計算L的元素ljj需要用到開方運算.為了避免開方，采用定理10的分解式A=LDLT.即由矩陣乘法，并注意ljj=1,ljk=0(j<k)，得2023/7/20117第117頁，課件共155頁，創作于2023年2月于是得到L的元素及D的對角元素公式：為了避免重復計算，我們引進tij=lijdj

.由(3.9)得到按行計算L,T(T=LD)元素的公式：d1=a11

對于i=1,2,…,n.2023/7/20118第118頁，課件共155頁，創作于2023年2月求解Ly=b及DLTx=y的計算公式為計算公式(3.10),(3.11)稱為改進的平方根法.對于i=1,2,…,n.2023/7/20119第119頁，課件共155頁，創作于2023年2月6.4誤差分析6.4.1矩陣的條件數由于A(或b)元素是測量得到的,或者是計算的結果,在第一種情況A(或b)常帶有某些觀測誤差,在后一種情況A(或b)又包含有舍入誤差.因此我們處理的實際矩陣是A+A(或b+b),下面我們來研究數據A(或b)的微小誤差對解的影響.即考慮估計x-y,其中y是(A+A)y=b的解.考慮線性方程組Ax=b,其中設A為非奇異矩陣,x為方程組的精確解.2023/7/20120第120頁，課件共155頁，創作于2023年2月首先考察一個例子.

例8設有方程組記為Ax=b,它的精確解為x=(2,0)T

.現在考慮常數項的微小變化對方程組解的影響,即考察方程組2023/7/20121第121頁，課件共155頁，創作于2023年2月也可表示為A(x+x)=b+b,其中b=(0,0.0001)T,y=x+x,x為(6.1)的解.顯然方程組(5.2)的解為x+x=(1,1)T.我們看到(5.1)的常數項b的第2個分量只有1/1000的微小變化,方程組的解卻變化很大.這樣的方程組稱為病態方程組.

定義7如果矩陣A或常數項b

的微小變化(小擾動),引起方程組Ax=b解的巨大變化,則稱此方程組為“病態”方程組,其系數矩陣A稱為“病態”矩陣(相對于方程組而言),否則稱方程組為“良態”方程組,A稱為“良態”矩陣.2023/7/20122第122頁，課件共155頁，創作于2023年2月應該注意,矩陣的“病態”性質是矩陣本身的特性,下面我們希望找出刻畫矩陣“病態”性質的量.設有方程組

Ax=b,(5.3)其中A為非奇異矩陣,x為(5.3)的精確解.以下我們研究方程組的系數矩陣A(或b)的微小誤差(小擾動)時對解的影響.2023/7/20123第123頁，課件共155頁，創作于2023年2月

(1)現設A是精確的,x為Ax=b的精確解，當方程組右端有誤差b,受擾解為

x+x,則

A(x+x)=b+b,Ax=b,x=A-1b,

||x||≤||A-1||||b||.(5.4)由(5.3)有||b||≤||A||||x||.于是由(5.4)及(5.5),得2023/7/20124第124頁，課件共155頁，創作于2023年2月定理21設A是非奇異矩陣,Ax=b≠0

,且

A(x+x)=b+b,則上式給出了解x的相對誤差的上界,常數項b的相對誤差在解中放大||A-1||||A||倍.

(2)現設b是精確的,x為Ax=b的精確解，當A有微小誤差(小擾動)A,受擾解為

x+x,則(A+A)(x+x)=b,(A+A)x=-(A)x.(5.6)2023/7/20125第125頁，課件共155頁，創作于2023年2月如果A不受限制的話,可能A+A奇異,而(A+A)=A(I+A-1A).由定理20知,||A-1A||<1時,(I+A-1A)-1存在.由(5.6)式x=-(I+A-1A)-1A-1(A)x.因此設||A-1||||A||<1,即得2023/7/20126第126頁，課件共155頁，創作于2023年2月定理22設A是非奇異矩陣,Ax=b≠0

,且

(A+A)(x+x)=b.如果||A-1||||A||<1,則(5.7)式成立.如果A充分小,且在條件||A-1||||A||<1下,那么(5.7)式說明矩陣A的相對誤差||A||/||A||在解中可能放大||A-1||||A||倍.總之,量||A-1||||A||越小,由A(或b)的相對誤差引起的解的相對誤差就越小;量||A-1||||A||越大,解的相對誤差就可能越大.所以量||A-1||||A||事實上刻畫了解對原始數據變化的靈敏程度,即刻畫了方程組的“病態”程度,于是引進下述定義:2023/7/20127第127頁，課件共155頁，創作于2023年2月

(3)現設x為Ax=b的精確解，當A有微小誤差(小擾動)A,而b同時也有微小誤差b(小擾動)時,受擾解為

x+x,則還可以推出相對誤差估計式為

補充2023/7/20128第128頁，課件共155頁，創作于2023年2月定義8設A是非奇異矩陣,稱數cond(A)v=||A-1||v||A||v(v=1,2或)為矩陣A的條件數

.由此看出矩陣的條件數與范數有關.矩陣的條件數是一個十分重要的概念.由上面討論知,當A的條件數相對的大,即cond(A)>>1時,則(5.3)是“病態”的(即A是“病態”矩陣,或者說A是壞條件的,相對于方程組),當A的條件數相對的小,則(5.3)是“良態”的(或者說A是好條件的).注意,方程組病態性質是方程組本身的特性.A的條件數越大,方程組的病態程度越嚴重,也就越難用一般的計算方法求得比較準確的解.2023/7/20129第129頁，課件共155頁，創作于2023年2月例如對前面例8的矩陣作分析由于條件數cond(A)1很大,可見矩陣A的病態程度十分嚴重,故由此方程組的解誤差非常大.計算其條件數2023/7/20130第130頁，課件共155頁，創作于2023年2月通常使用的條件數,有(1)cond(A)∞

=||A-1||||A||;(2)A的譜條件數;當A為對稱矩陣時其中1,n為A的絕對值最大和絕對值最小的特征值.2023/7/20131第131頁，課件共155頁，創作于2023年2月

條件數的性質:(1).對任何非奇異矩陣,都有cond(A)v≥1.事實上

cond(cA)v=cond(A)v.(2).設A為非奇異矩陣且c≠0(常數),則(3).如果A為正交矩陣,則cond(A)2=1;如果A為非奇異矩陣,R為正交矩陣,則

cond(RA)2=cond(AR)2=cond(A)2.2023