这个系列是学习《数值分析》时,将所学内容进行简单记录,以及自己用matlab的一些实现过程。
这次是这周学的三角分解,matlab有对应的函数,在这里主要记录一下大概的分解结果和方法,最后自己用代码实现了一下。具体的求解方法和直接三角分解的方法一致,都是写出LU的形式,然后展开LU并与A对应,求取L和U的元素。
直接三角分解法
如果A是n阶方阵,并且A的前n阶顺序主子式均不为0,则A存在唯一的分解A=LU
A
=
[
1
l
21
1
l
31
l
32
1
⋮
⋮
⋮
⋱
l
n
1
l
n
2
l
n
3
⋯
1
]
[
u
11
u
12
u
13
⋯
u
1
n
u
22
u
23
⋯
u
2
n
u
33
⋯
u
3
n
⋱
⋮
u
n
n
]
A= \begin{bmatrix} 1& & & & \\ l_{21}& 1& & & \\ l_{31}& l_{32} & 1& & \\ \vdots & \vdots& \vdots& \ddots& \\ l_{n1}& l_{n2}& l_{n3}& \cdots &1 \end{bmatrix} \begin{bmatrix} u_{11}& u_{12}& u_{13}& \cdots & u_{1n} \\ & u_{22}& u_{23}& \cdots &u_{2n} \\ & & u_{33}& \cdots & u_{3n}\\ & & & \ddots & \vdots \\ & & & &u_{nn} \end{bmatrix}
A=⎣⎢⎢⎢⎢⎢⎡1l21l31⋮ln11l32⋮ln21⋮ln3⋱⋯1⎦⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎡u11u12u22u13u23u33⋯⋯⋯⋱u1nu2nu3n⋮unn⎦⎥⎥⎥⎥⎥⎤
前边的下三角矩阵是L,后边的上三角矩阵是U。
求解时,先求U的第一行,用L的第一列与U做乘法,得到
A
(
1
,
:
)
=
L
(
1
,
:
)
U
=
U
(
1
,
:
)
A(1,:)=L(1,:)U=U(1,:)
A(1,:)=L(1,:)U=U(1,:),即U的第一行与A的第一行相等。
然后求L的第一列,用L与U的第一列做乘法,即
A
(
:
,
1
)
=
L
∗
U
(
:
,
1
)
=
L
(
:
,
1
)
∗
u
11
A(:,1)=L*U(:,1)=L(:,1)*u_{11}
A(:,1)=L∗U(:,1)=L(:,1)∗u11,所以,U的第一列就是
A
(
:
,
1
)
/
u
11
A(:,1)/u_{11}
A(:,1)/u11,第一个元素是1,也可以不用计算,
A
(
2
:
e
n
d
,
1
)
∗
U
(
:
,
1
)
A(2:end,1)*U(:,1)
A(2:end,1)∗U(:,1)。
对于L和U中间行列的计算,也是同理,假设前L的前r-1列和U的前r行已经得到,现在求取L的第r列和U的第r行。
求U的第r行第j列的元素(j>=r)
a
r
,
j
=
L
(
r
,
:
)
∗
U
(
:
,
j
)
=
[
l
r
,
1
⋯
l
r
,
r
−
1
1
0
⋯
0
]
∗
[
u
1
,
j
⋯
u
r
−
1
,
j
u
r
,
j
0
⋯
0
]
T
=
∑
m
=
1
r
−
1
(
l
r
,
m
u
m
,
j
)
+
u
(
r
,
j
)
a_{r,j}=L(r,:)*U(:,j)\\=\begin{bmatrix} l_{r,1}&\cdots &l_{r,r-1} & 1 &0&\cdots &0 \end{bmatrix}*\begin{bmatrix} u_{1,j}&\cdots &u_{r-1,j} & u_{r,j} &0&\cdots &0 \end{bmatrix}^T\\=\sum_{m=1}^{r-1}(l_{r,m}u_{m,j})+u(r,j)
ar,j=L(r,:)∗U(:,j)=[lr,1⋯lr,r−110⋯0]∗[u1,j⋯ur−1,jur,j0⋯0]T=m=1∑r−1(lr,mum,j)+u(r,j)
所以
u
r
,
j
=
a
r
,
j
−
∑
m
=
1
r
−
1
(
l
r
,
m
u
m
,
j
)
u_{r,j}=a_{r,j}-\sum_{m=1}^{r-1}(l_{r,m}u_{m,j})
ur,j=ar,j−m=1∑r−1(lr,mum,j)
对于L的第i行也是类似,求L的第r列的第i个元素,r<i<=n。
a
i
,
r
=
L
(
i
,
:
)
∗
U
(
:
,
r
)
=
[
l
i
,
1
⋯
l
i
,
r
⋯
l
i
,
i
0
⋯
0
]
∗
[
u
1
,
r
⋯
u
r
,
r
0
⋯
0
]
T
=
∑
m
=
1
r
−
1
(
l
i
,
m
∗
u
m
,
r
)
+
l
i
,
r
∗
u
r
,
r
a_{i,r}=L(i,:)*U(:,r)\\=\begin{bmatrix} l_{i,1}&\cdots &l_{i,r} & \cdots & l_{i,i}&0&\cdots &0 \end{bmatrix}*\\\begin{bmatrix} u_{1,r}&\cdots &u_{r,r} & 0&\cdots &0 \end{bmatrix}^T\\=\sum_{m=1}^{r-1}(l_{i,m}*u_{m,r})+l_{i,r}*u_{r,r}
ai,r=L(i,:)∗U(:,r)=[li,1⋯li,r⋯li,i0⋯0]∗[u1,r⋯ur,r0⋯0]T=m=1∑r−1(li,m∗um,r)+li,r∗ur,r
所以
l
i
,
r
=
a
i
,
r
−
∑
m
=
1
r
−
1
(
l
i
,
m
∗
u
m
,
r
)
u
r
,
r
l_{i,r}=\frac{a_{i,r}-\sum_{m=1}^{r-1}(l_{i,m}*u_{m,r})}{u_{r,r}}
li,r=ur,rai,r−∑m=1r−1(li,m∗um,r)
下边是一张图比较清晰地说明这个步骤:
比如求U的第r行第j列时,用这一行左边黄色部分(属于L的前r-1列)和第r行第j列上方的黄色部分(属于U的前r-1行)作乘积,然后用A(r,j)减去这个乘积就是
u
r
,
j
u_{r,j}
ur,j。
求L的第i行第r列时,用这一列上方的蓝色部分(属于U的前r-1行)与对应列的蓝色部分(属于L的前r-1列)做乘积,假设值是sum,然后
l
i
,
r
=
a
i
,
r
−
s
u
m
u
r
,
r
l_{i,r}=\frac{a_{i,r}-sum}{u_{r,r}}
li,r=ur,rai,r−sum
自己用matlab实现的代码:
%对应于高斯消元法的LU分解,结果都存在A中
function A = LU(A)
%U的第一行就是A的第一行
%一个循环处理L的第i列(i+1:n,i),U的第i+1行(i+1,i+1:n)
[rows,cols] = size(A);
if rows~=cols
error("rows != cols");
end
for i=1:rows-1
if i==1
A(i+1:end,i)=A(i+1:end,i)/A(i,i);%第一列
A(i+1,i+1:end) = A(i+1,i+1:end)- A(i+1,1:i)*A(1:i,i+1:end);
continue
end
%列
A(i+1:end,i) = (A(i+1:end,i) - A(i+1:end,1:i-1)*A(1:i-1,i))/A(i,i);
A(i+1,i+1:end) = A(i+1,i+1:end)- A(i+1,1:i)*A(1:i,i+1:end);
end
end
这里直接将LR合并写到了A中,方便和手动计算题目作比较。
列主元的三角分解
列主元的三角分解和之前列主元法一样,只不过是把消元过程的系数存储到A的下半部分,并且在做行变换时,对当前A的所有列都做行变换(即L也参与行变换)。列主元法的三角分解可以写为:
P
A
=
L
U
PA=LU
PA=LU
从这里也可以看出为什么要对L进行行变换,普通的三角分解假设是
A
=
l
U
A=lU
A=lU,则列主元法的三角分解为
P
A
=
P
l
U
=
L
U
PA=PlU=LU
PA=PlU=LU。
下边是实现的代码,在上一篇中列主元法的基础上稍作修改。
%对应于列主元高斯消元法的LU分解,结果都存在A中
%和普通的相比,就是多了一部初等行交换。
%PA=LU,P是多次初等行变换叠加的矩阵
%相对于LU分解,这里的L也要行交换,多了次找列主元
%L在A的下三角部分,U在A的对角线及以上
function [P,A] = lzy_LU(A)
%U的第一行就是A的第一行
%一个循环处理L的第i列(i+1:n,i),U的第i+1行(i+1,i+1:n)
[rows,cols] = size(A);
if rows~=cols
error("rows != cols");
end
P=eye(rows);
for i=1:rows-1
%第i行,第i列
[M,idx]=max(abs(A(i:end,i)));
idx=idx+i-1;
%注意这里是全部交换
A([i,idx],:) = A([idx,i],:);
P([i,idx],:)=P([idx,i],:);
A(i+1:end,i) = A(i+1:end,i)/A(i,i);%L的第i列
A(i+1:end,i+1:end) = A(i+1:end,i+1:end)-A(i+1:end,i)*A(i,i+1:end); %更新U及A
end
end
三对角矩阵的三角分解(追赶法)
如果矩阵A的形式为
[
b
1
c
1
a
2
b
2
c
2
⋱
⋱
⋱
a
n
−
1
b
n
−
1
c
n
−
1
a
n
b
n
]
\begin{bmatrix} b_1& c_1& & & \\ a_2& b_2& c_2& & \\ & \ddots& \ddots &\ddots & \\ & & a_{n-1}& b_{n-1}&c_{n-1} \\ & & & a_n &b_n \end{bmatrix}
⎣⎢⎢⎢⎢⎡b1a2c1b2⋱c2⋱an−1⋱bn−1ancn−1bn⎦⎥⎥⎥⎥⎤
并且满足
∣
b
1
∣
>
∣
c
1
∣
>
0
,
∣
b
i
∣
≥
∣
a
i
∣
+
∣
c
i
∣
,
∣
b
n
∣
>
∣
a
n
∣
>
0
|b_1|>|c_1|>0,|b_i|\ge|a_i|+|c_i|,|b_n|>|a_n|>0
∣b1∣>∣c1∣>0,∣bi∣≥∣ai∣+∣ci∣,∣bn∣>∣an∣>0,则A可分解为:
L
U
=
[
1
l
2
1
⋱
⋱
l
n
−
1
1
l
n
1
]
[
u
1
c
1
u
2
c
2
⋱
⋱
u
n
−
1
c
n
−
1
1
u
n
]
LU=\begin{bmatrix} 1& & & & \\ l_2& 1& & & \\ & \ddots& \ddots &\ & \\ & & l_{n-1}& 1&\\ & & & l_n &1 \end{bmatrix} \begin{bmatrix} u_1&c_1 & & & \\ & u_2& c_2 & & \\ & & \ddots &\ddots & \\ & & & u_{n-1}&c_{n-1}\\ & & & &1u_n \end{bmatrix}
LU=⎣⎢⎢⎢⎢⎡1l21⋱⋱ln−1 1ln1⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎡u1c1u2c2⋱⋱un−1cn−11un⎦⎥⎥⎥⎥⎤
根据这个形式可以仿照上边计算LU矩阵,并能简化。
u
1
=
b
1
,
l
2
=
a
2
/
u
1
,
u
2
=
b
2
−
c
1
∗
l
2
u_1=b_1,l_2=a_2/u_1,u_2=b_2-c_1*l_2
u1=b1,l2=a2/u1,u2=b2−c1∗l2
l
k
=
a
k
/
u
k
−
1
,
u
k
=
b
k
−
l
k
∗
c
k
−
1
l_k=a_k/u_{k-1},u_k=b_k-l_k*c_{k-1}
lk=ak/uk−1,uk=bk−lk∗ck−1
这里计算
l
k
l_k
lk时,
l
k
l_k
lk左边全为0,所以
a
k
a_k
ak不用减左边向量和当前列上边向量的乘积。计算
u
k
u_k
uk时,
u
k
u_k
uk左边为
l
k
l_k
lk,上边为
c
k
−
1
c_{k-1}
ck−1,内积也只有这两个值相乘。
对称正定矩阵的三角分解
后边几个晚上再写。