计算方法引论-第六章

Gauss消去法求行列式
• 行列式
– det(Ak)=det(Uk)=u11u22…ukk , k=1,2,…,n U是顺序消元过程结束时的上三角矩阵. Ak 和Uk分别是A和U的k阶主子阵
– 例3中系数矩阵的行列式等于2×3×6=36. – 两个等价的可顺序消元的充分条件
主元(消第k个元时左上角的元素)皆非零 顺序主子式非零,det(Ak)≠0,k=1,2,…,n-1
计算方法引论：数值代数
解线性方程组的直接法 解线性方程组最小二乘问题 解线性方程组的迭代法 矩阵特征值和特征向量的计算 非线性方程及非线性方程组解法
第六章 解线性方程组的直接法
• Gauss消去法 • 主元素法 • LU分解 • LLT分解和LDLT分解 • 误差分析
计算方法引论( 第三版)
• 例3(续) 2 2 3 1
2 2 3
4
7
7
2
1
3 1
2 4 5 1 2 1
6
• LU分解:顺序主子式非零,det(Ak)≠0,k=1,2,…,n-1则可
唯一分解A=LU,单位下三角阵与上三角阵之积
1
l21 1
(3) L l31
l32
1
u11 u12 u13 u1n
u22
找 p: apk max akk , ak1,k ,, ank
○p ○k
ik=p for j=k+1:n
ukj=akj-lk1u1j-lk2u2j -…- lk,k-1uk-1,j end（ukk=akk） for i=k+1:n
• 应用
– 解Ax=b:分解A=LU 解Ly=b求y 解Ux=y求x
– 计算det(A)= det(L) det(U)= u11u22…unn
计算方法引论( 第三版)
6.18
徐萃薇、孙绳武 高教2007
直接LU分解算例
• 例4(对照例3,右端项计算与U同)
223
3
7
477
1
18
-2 4
5
-7
7
223
end for i=1:j–1
aij = aij/aii end end
(计算 rij)
计算方法引论( 第三版)
6.20
徐萃薇、孙绳武 高教2007
小主元扩大误差
• 例5 顺序消去法,用精确运算:
0.0001 1 1 0.0001 1
1
1
1
2
0
9999 9998
得(10000/9999,9998/9999)≈(1.0001,0.9999)若
0
2
6
(1)
1 1 1 2
2 1 1 7 (1) 1 1 1 (1/ 2) 0 2 6 3
括号内是乘数,k=2时2,3行交换.
•LU分解
1
2 2 0 1
2 2 0
11
1
2
1
1
1 1
1 2 1 1 1/ 2 0 1
2
即列主元素法实现了LU分解:PA=LU ,P是行交换结
果的排列阵.
计算方法引论( 第三版)
分开进行 – 运算量: n(n-1)/2除法,加法乘法次数
对系数矩阵各为n(n-1/2)(n-1)/3 对右端项各为n(n-1)/2
•Gauss消去法
– 顺序消元.需加法乘法各为n3/3+O(n2) – 回代.需加法乘法各为n2/2 +O(n)
计算方法引论( 第三版)
6.11
徐萃薇、孙绳武 高教2007
6.2
徐萃薇、孙绳武 高教2007
引言
• 解线性代数方程组 a11x1 a12x2 a1n xn b1
(1)
a21x1
a22x2
a2n xn
b2
an1x1 an2x2 ann xn bn
• 矩阵方程
Ax=b
常用增广矩阵表示 (A∣b)
• 解法
– 直接法:用有限步计算得到准确解
– 回代加法和乘法运算各n(n-1)/2,除法n次 – 亦可解出一未知数即代入其它方程,消
去该未知数 – 其它形式三角方程组可类似计算
计算方法引论( 第三版)
6.7
徐萃薇、孙绳武 高教2007
Gauss消去法算例
• 例3
原方程组(二右边项) k=1: 2行减1行2倍
3行减1行(-1)倍 k=2:
3行减2行2倍 回代
6.24
徐萃薇、孙绳武 高教2007
列主元LU分解
• 列主元素法实现LU分
解
– 如上例,只要记住 交换历史
• 直接列主元LU分解
– 修改直接LU分解 加入选主元
– 算法如右(可就地 完成:lik,ukj置A中)
for k=1:n-1 for i=k:n aik=aik-li1u1k-li2u2k -…-li,k-1uk-1,k end
• 算法:当u11u22…unn≠0时,可解出
xn=gn/unn for k=n-1:1
xk=(gk- uk,k+1xk+1-…- uknxn)/ ukk end
计算方法引论( 第三版)
6.6
徐萃薇、孙绳武 高教2007
三角方程组(续)
•算法注记
– 程序实现时x,g可共用一组单元,即回代 就地完成
mik=aik/akk
{mik可置于aik,如(2)}
i行=i行-k行×mik {前k列元素不在内}
end
else stop end
end
计算方法引论( 第三版)
6.10
徐萃薇、孙绳武 高教2007
Gauss消去法运算量
• 顺序消元算法注记
– mik置于aik 则顺序消元就地完成 – 顺序消元算法也可将系数矩阵和右端项
l31=a31/u11 l32=(a32-l31u12)/ u22 u33=a33-l31u13- l32u23
• LDR
分解
for j = 1:n for i=2:j
aij aij ai1a1j ai2a2 j
end
ai,i1ai1, j
(计算 uij)
for i= j+1:n aij (aij ai1a1 j ai 2a2 j ai, j1a j1, j ) / ajj (计算 lij)
计算方法引论( 第三版)
6.17
徐萃薇、孙绳武 高教2007
直接LU分解 (续)
• 计算表格
u11=a11
u12=a12
u13=a13
l21=a21/u11
u22=a22-l21u12
u23=a23-l21u13
l31=a31/u11 l32=(a32-l31u12)/ u22 u33=a33-l31u13- l32u23 – 也可逐行算,或逐列算,或其它可行次序算
u23
u2
n
,U
u33
u3n
ln1 ln2 ln3 1
unn
计算方法引论( 第三版)
6.15
徐萃薇、孙绳武 高教2007
LU分解
• 方法
– 消去法实现LU分解 – 直接LU分解(紧凑Gauss消去法)
解方程A=LU确定lik,ukj 追踪顺序消元所得L,U元素的历史确定lik,ukj
ak,k1 akn
a a k1,k1
k 1,n
bk
bk 1
akk (mk1,k
)
ak ,k1
ak 1,k 1
kk
ak
n
ak1,n
bk bk 1
ank
an,k1 ann
bn (mn1)
an,k1 ann
bn
这里各行变换:i行-k行×mik ,结果置原地
其中 mik=aik/akk,i=k+1,…,n
对值最大的做主元(交换行及列). – 这对误差控制有利,但搜索太费时.通常列主 元素法误差控制就已可以了.
计算方法引论( 第三版)
6.22
徐萃薇、孙绳武 高教2007
列主元素法
for k=1:n-1
找 p: apk max akk , ak1,k ,, ank
○p ○k ik=p
if akk≠0
• 直接LU分解公式
ukj=akj-mk1u1j- mk2u2j -…- mk,k-1uk-1,j, j=k,k+1,…,n
mik=(aik-mi1u1k- mi2u2k -…-mi,k-1uk-1,k)/ukk i=k+1,k+2,…,n
计算方法引论( 第三版)
6.16
徐萃薇、孙绳武 高教2007
在十进三位尾数舍入的浮点计算机系统中运
算,第二行将是(0 -10000∣-10000)得到解
x2=1, x1=0.与真解相去甚远. – 把两个方程(两行)交换次序再消元,得解
x2=1, x1=1,与真解很近:
0.0001 1 1 1 1 2 1 1 2
1
1 2 0.0001 1 1 0 1 1
3
7
2
3
1
-5
4
-1 2 6
6
6
1
2
1
1
-2
1
1
1
1
det(A)= 2×3×6=36
计算方法引论( 第三版)
6.19
徐萃薇、孙绳武 高教2007
其它形式的分解
• A=LU(L下三角矩阵,U单位上三角矩阵 )
u11=a11 l21=a21/u11
u12=a12 u22=a22-l21u12
u13=a13 u23=a23-l21u13
计算方法引论( 第三版)
6.21
徐萃薇、孙绳武 高教2007
主元素法
• 列主元素法:在每次消元前先选该列中绝对值
最大的做主元(交换两行,每行包括右端项!). – 列 加主 误元 差素法乘数mik绝对值不大于1,不会增 – 列主元素法用来求行列式时要注意两行交 换行列式变号.
• 全主元素法:在整个右下(n-k)×(n-k)矩阵找绝
直接LU分解
• 算法:
for k=1:n-1 for j=k:n ukj=akj-lk1u1j-lk2u2j -…- lk,k-1uk-1,j end for i=k+1:n lik=(aik-li1u1k-li2u2k -…-li,k-1uk-1,k)/ukk end
end
ukj,li可置A中.
– 消元结果得上三角方程组
M(A∣b)=(U∣g),M=Mn-1Mn-2…M1
MA=U, Mb=g
1
A=LU,L=M-1
m21 1
=M1-1…Mn-2-1Mn-1-1= m31
m32
1
mn1 mn2 mn3 1
计算方法引论( 第三版)
6.14
徐萃薇、孙绳武 高教2007
消去法实现了LU分解
– 对k=2,对右下n-1阶方程组消第k个元
(akk≠0),结果置原地.
– 对k=3,…,n-1,对右下n-(k- 1)阶方程组消
第k个元(akk≠0),结果置原地.
计算方法引论( 第三版)
6.9
徐萃薇、孙绳武 高教2007
顺序消元(续)
• 算法
for k=1:n-1
if akk≠0 for i=k+1:n
0
3
1
2 1
0
0 0 6 5 2 1
1 0 0 7 / 36 1/18 7 / 36
0
1
0
Байду номын сангаас
17 /18
4/9
1 / 18
0 0 1 5 / 6 1/ 3 1/ 6
计算方法引论( 第三版)
6.13
徐萃薇、孙绳武 高教2007
Gauss消去法矩阵解释
•消去法矩阵解释
– 消第k个元,(2),相当于左乘矩阵Mk=I-mkekT mk=(0, … ,0, mk+1,k ,… , mnk) T, ek单位向量
徐萃薇、孙绳武 高教2007
三角方程组(例)
• 例1:向前代入消去法
1
3 1
3
2
1
1
1
5
1 2 1 7
1
6
• 例2:向后代入消去法(回代)
2 2 3 3 1
2
3
1
5
1
2
6
6
1 1
计算方法引论( 第三版)
6.5
徐萃薇、孙绳武 高教2007
三角方程组
• 一般情况
u11x1+u12x2+…+u1nxn=g1 u22x2+…+u2nxn=g2 …… unnxn=gn
– 迭代法:给出一个近似解序列
计算方法引论( 第三版)
6.3
徐萃薇、孙绳武 高教2007
直接法
• Crame法则
计算量太大,以(n+1)!计,不实用
(11!=39 916 800).
• 消去法
计算量以n3计 .
• 直接法
– 解线性代数方程组 – 求行列式 – 求逆矩阵
计算方法引论( 第三版)
6.4
for i=k+1:n mik=aik/akk ○i =○i -○k ×mik
end else stop end end
计算方法引论( 第三版)
6.23
徐萃薇、孙绳武 高教2007
列主元素法算例
•例6 列主元素法解方程组
2 2 0 6 2 2 0 6 2 2 0 6 1
1 1 2
9 (1/ 2)
计算方法引论( 第三版)
6.12
徐萃薇、孙绳武 高教2007
Gauss消去法求逆矩阵
• 逆矩阵
– 解n个方程组 (A∣I),其中I是单位矩阵.
– 需加法乘法各为n3+O(n2)
– 例3中系数矩阵的逆,写出增广矩阵,消元得
2 2 3 1 0 0
4
7
7
01
0
2
45 0
最后解出
0
1
2 2 3 1 0 0
223 3 7
4 7 7 1 18
-2 4 5 -7 7
223 3 7
2 3 1 -5 4
-1 6 8 -4 14
223 3 7
2 3 1 -5 4
-1 2 6 6 6
1
21
1
-2 1
11 1
计算方法引论( 第三版)
6.8
徐萃薇、孙绳武 高教2007
顺序消元
• 典型步
– 对k=1,对n阶方程组消第k个元(a ≠0): akk (2) ak1,k