第三章-数值分析(08)用矩阵分解法解线性代数方程组

f (1) b(P(1)) b(3), f (2) b(P(2)) b(4)
f (3) b(P(3)) b(1), f (4) b(P(4)) b(2)
数值分析
数值分析
1
1
解LY
1 2
0
12
1
2 7
17
1 0
y1
y2
1
y3 y4
1
6
0
2
Pb
x2
d
2
OOO
M M
O
O
cn1
M
M
an bn xn dn
用矩阵表示 Ax d
应用追赶法求解三对角线性方程组。追赶法仍然 保持LU分解特性,它是一种特殊的LU分解。充分利用 了系数矩阵的特点，而且使之分解更简单,得到对三对 角线性方程组的快速解法。
数值分析
数值分析
数值分析
数值分析
求解 Ly=b, 得 y1=1, y2=1.5, y3=1, y4=0.5 求解 Ux=y , 得 x4=0.3333, x3=-0.3333, x2=-1, x1=-1
数值分析
数值分析
周期三对角方程组的一般形式
b1 c1 a2 b2 c2
x1 d1
x2
d2
数值分析
数值分析
2 x1 x2
1
例：求解方程组：
x1
2 x2
1
x2 2 x3 x4 0
x3
2
x4
1
2 1
1
u1 r1
解：A
1
2
LU
l2
1
u2 r2
1 2 1
1 2
l3
1 l4
1
u3
r3 u4
1
2 1
1
2
1 23
1 12
1
32 0 2
1
3
2
z(i)=(y(i)-LU(i,i+1:n)*z(i+1:n)')/LU(i,i); x(q(i))=z(i); end
数值分析
数值分析
四、利用Cholesky分解A LLT求解Ax b
A Rnn是对称正定的矩阵，有Cholesky分解式
A LLT代入原方程 LLT x b,可分成两步
(1)求解下三角方程组
LY b
i 1
yi (bi lik yk ) / lii , (i 1, 2,L , n) k 1
(2)求解上三角方程组
LT x Y
n
xi ( yi lki xk ) / lii , (i n, n 1,L ,1) k i 1 数值分析
数值分析
五、利用正交分解A QR求解 Ax b Ax b QRx b Rx QTb，于是可分成三步
OO O
M M
an1
bn1
cn1
xn1
d
n1
an bn xn dn
矩阵表示 Ax d
ck rk ,
k 2, 3,L , n 1
u1 b1, rk ck (k 1, 2,L , n 1) 对k 2, 3,L , n
lk ak uk1 , uk bk lkck1
追赶法的计算过程分为三步
(1)A LU (2)LY d (3)Ux Y
数值分析
数值分析
求解三对角方程组的追赶法 ①计算分解因子阵
较常见带状矩阵为带宽为3（p=q=2,w=3)的矩阵。
a11
a12
0
0
a21
a22
a23
A 0
0
A称为三对
角矩阵。 an1,n2
an1,n1
an1,n
0
0
an,n1
an,n
系数矩阵为三对角矩阵的线性方程组称为三对角方程组。
数值分析
数值分析
三对角线性方程组
b1 c1
a2
b2
c2
x1 d1
U
0 0 11 7
数值分析
数值分析
L2 P2 L1P1 A U , A P11L11P21L21U P1 L11P2 L21U 令P P2 P1 PA P2 P1P1 L11P2 L21U P2 L11P2 L21U ， L P2 L11P2 L21
PA LU
1
L
P2 L11P2 L21
x(i)=(y(i)-LU(i,i+1:n)*x(i+1:n)')/LU(i,i); end
数值分析
数值分析
三、用全主元的三角分解PAQT LU求解Ax b Ax b PAQT (Qx) Pb LU(Qx) Pb
lupqdsv.m
%功能：调用全主元三角分解函数[LU,p,q]=lupqd(A)
数值分析
数值分析
function x=lupqdsv(A,b) n=length(b); [LU,p,q]=lupqd(A); y(1)=b(p(1)); for i=2:n
y(i)=b(p(i))-LU(i,1:i-1)*y(1:i-1)'; end z(n)=y(n)/LU(n,n);x(q(n))=z(n); for i=(n-1):-1:1
数值分析
第三节 用矩阵分解法求解线性方程组
一、利用三角分解A LU求解Ax b 二、用列主元的三角分解PA LU求解Ax b
三、用全主元的三角分解PAQT LU求解Ax b
四、利用Cholesky分解A LLT求解Ax b 五、利用正交分解A QR求解 Ax b 六、利用矩阵奇异值分解A UV T求解Ax b 七、 三对角方程组的解法
1,n
O
O
O M
0
ln,n p1 L 1 0
unn
数值分析
数值分析
当A为三对角阵,且 b1 c1 , bi ci ai ，(i 1, 2,L , n 1),
bn an 时，A有LU分解展开式
b1 c1
a2
b2
c2
1
l2
1
u1 r1
u2 r2
OOO
O
O
cn1
l3 O OO
OO
O
rn1
an bn
ln 1
un
b1 u1 , c1 r1, a2 l2u1
b2 l2r1 u2 l2c1 u2, c2 r2 ,
ak lk uk1 ,
k 2, 3,L , n
bk lk rk1 uk lkck1 uk , k 2, 3,L , n
数值分析
数值分析
function x=lupdsv(A,b) n=length(b); [LU,p]=lupd(A); y(1)=b(p(1)); for i=2:n
y(i)=b(p(i))-LU(i,1:i-1)*y(1:i-1)'; end x(n)=y(n)/LU(n,n); for i=(n-1):-1:1
ck rk ,
k 2, 3,L , n 1
数值分析
数值分析
b1 u1 , c1 r1, a2 l2u1
b2 l2r1 u2 l2c1 u2, c2 r2 ,
ak lk uk1 ,
k 2, 3,L , n
bk lk rk1 uk lkck1 uk , k 2, 3,L , n
u1 b1 , rk ck (k 1, 2,L , n 1) 对k 2, 3,L , n lk ak uk 1 , uk bk lk ck 1 ②求解LY d , y1 d1 yk dk lk yk1(k 2, 3,L , n) ③求解Ux Y , xn yn un xk ( yn ck xk1 ) uk (k n 1,L , 2,1)
（1）正交分解 A QR （2）矩阵乘法 b QT b （3）回代求解上三角方程组 Rx b
这种方法没有选主元问题，计算稳定，但计算量略
大，消元法的计算量大约为 大约为 2 n3。
3
1 n3，正交分解法的计算量 3
数值分析
数值分析
(1) A QR
q11 q12 L
Q
q21
q22
L
L L L
于是，求解 Ax b x A1b
可化为 x V diag( 1 , 1 ,..., 1 )UTb
1 2 n
数值分析
数值分析
七、 三对角方程组的解法
定义1 若n 阶矩阵A=(aij)的元素满足:对于1<p,q<n的
正整数p、q,有j≥i+p及i≥j+q时，aij=0，则A称为带
状矩阵. 带宽为w=p+q-1。
y1
f
,
得
y2
y3 y4
113327
117
2 解Ux
1 7
2 0
0
0 0 1 0
1
1 2
13 7
11 7
x1 x2 x3 x4
=
-1
13Байду номын сангаас2
13 7 117
,得
x1 x2 x3 x4
=
u11 u12 u1n x1 y1
u22
u2n
x2
y2
unn
xn
yn
数值分析
数值分析
二、用列主元的三角分解PA LU求解Ax b
Ax
b
PAx
Pb
LUx
Pb
LY
Ux
Pb Y
例：用列主元三角分解求解Ax=b
0 1 1 2
其中A 1 0 0 1 2 1 0 1
qn1
qn 2
L
(2)计算b QT b
q1n
r11 r12 L
q2n
,
R
r22
L
L
qnn
r1n
L
L
rnn
n
bk qikbi , k 1, 2,L , n i 1
(3)回代求解Rx b
xn
bn
rnn
xk
(bn
n
rkj x j )
, rkk
(k n 1,L , 2,1)
1
2
0
1
数值分析
数值分析
lupdsv.m %功能：调用列主元三角分解函数 [LU,p]=lupd(A) % 求解线性方程组Ax=b。 %解法：PA=LU, Ax=b←→PAx=Pb % LUx=Pb, y=Ux % Ly=f=Pb, f(i)=b(p(i)) %输入：方阵A，右端项b（行或列向量均可） %输出：解x（行向量）
l31
l32
1
j1
1
ln1 ln2 ln,n1 1 yn bn
数值分析
数值分析
第 二 步: 求 解 上 三 角 方 程 组Ux Y ,向 后 回 代 求 出x
xn yn unn
n
xk ( yk ukj x j ) ukk j k 1
(k n 1, n 2, ,1)
% 求解线性方程组Ax=b。
%解法：PAQ-1=LU, Ax=b←→(PAQ-1)(Qx)=Pb
% LU(Qx)=Pb, z=Qx, y=Uz
% Ly=f=Pb, f(i)=b(p(i))
% Uz=y, z=Qx , x(q(i))=z(i).
%输入：方阵A，右端项b（行或列向量均可）
%输出：解x（行向量）
jk 1
数值分析
数值分析
六、利用矩阵奇异值分解A UV T求解Ax b
对于A Rnn非奇异情况，此时有奇异值分解
A UV T
1
其中U、V为正交阵，
O
,
i 0, i 1, 2, ..., n
n
显然
A1 V 1U T V diag( 1 , 1 , ..., 1 )U T
1 2 n
数值分析
数值分析
一、利用三角分解A LU求解Ax b 设已有 A LU 代入原方程 Ax b 得
LUx
b
LY Ux
b Y
第一步 : 求解下三角方程组LY b,向前回代求出
Y ( y1, y2 ,L , yn )T
1
y1 b1
y1 b1
l21
1
y2
b2
k 1
yk bk lkj y j (k 2, 3,L , n)
定理 如果带宽为 w=p+q-1 的n阶带状矩阵A有LU 分解:A=LU，则L是带宽为p的下三角矩阵，U是带宽 为q的上三角矩阵。
例：
a11 L M a22
a
p1
O
0
a1q O
O O
an,n p1 L
0
anq
1,n
M
ann
1
M
1
0 u11 L u1q
u22
O
0
l
p1
O
O
unq
1
3
0 1
解： i1 3, i1 1
0
b 2 1
6
0 0 1 0
P1
0 1
1 0
0 0
0 0
0 0 0 1
2 1 0 1
P1
A
1
0
0 1
0 1 1 2
1
3
0 1
数值分析
数值分析
1
L1
1 2 0
1 2
i2 4,
1
1
1
i2 2
2 1 0 1
L1P1 A
12 1
0 1
3 2 2
7 2
0
1
2
1 0 0 0
P2
0 0
0 0
0 1
1 0
0 1 0 0
2 1 0 1
P2 L1 P1 A
7 2 1
0 1
1 2 2
12
0
3
2
1
1
L2
27 1
1 7
1
2 1 0 1
L2 P2 L1 P1 A
7 2 0
0 1
1 2
13 7
1 2 0
1
2 7
1
1 2 17 0 1
数值分析
数值分析
P为排列阵，在计算机中用向量表示
例 P (1 2 3 4)T , P1 (3 2 1 4)T ,
P2 (3 4 1 2)T ,
P (3 4 1 2)T
Ax b, PA LU ,
PAx Pb,
LUx Pb f
f (i) b(P(i))