1
3
n
3
+
n
2
−
1
3
n
2
\frac{1}{3}n^3+n^2-\frac{1}{3}n^2
31n3+n2−31n2
2.主元素高斯消元法
为了控制舍入误差的扩大和传播而提出的。
列主元素高斯消元法
将该列绝对值尽可能大的系数作为第k步消元的主元素
3.高斯-约当(Gauss-Jordan)消去法
每次消去对角线下方和上方的元素 每次消元过程中,选取列主元素。
总计算量
1
2
n
3
+
1
2
n
2
\frac{1}{2}n^3+\frac{1}{2}n^2
21n3+21n2 计算量比高斯消去法计算量大,但不需要回代。
二、矩阵三角分解法
1.直接三角分解法(LU分解、Doolittle分解)
条件:
矩阵A的各阶顺序主子式都不为0(可逆)。
计算顺序
先计算第一行,从左到右依次计算
u
u
u ,
y
y
y。
再计算第一列,从上到下计算
l
l
l。
按照前两部依次计算第二行到最后一行。
计算公式
设矩阵
A
A
A 的增广矩阵
A
A
A 为
A
[
n
+
1
]
[
n
+
2
]
;
A[n+1][n+2];
A[n+1][n+2]; 下标为0的元素位置为0.
A
=
(
a
11
a
12
⋯
a
1
n
∣
b
1
a
21
a
22
⋯
a
2
n
∣
b
2
⋮
⋮
⋮
⋮
a
n
1
a
n
2
⋯
a
n
n
∣
b
n
)
\boldsymbol{A}=\left(\begin{array}{cccc} a_{11} & a_{12} & \cdots & a_{1 n} & |& b_1 \\ a_{21} & a_{22} & \cdots & a_{2 n} & |& b_2 \\ \vdots & \vdots & & \vdots & & \vdots \\ a_{n 1} & a_{n 2} & \cdots & a_{n n} & |& b_n \end{array}\right)
A=⎝⎜⎜⎜⎛a11a21⋮an1a12a22⋮an2⋯⋯⋯a1na2n⋮ann∣∣∣b1b2⋮bn⎠⎟⎟⎟⎞
则有:
for(int i=1;i<=n;i++){//i表示计算的层数//对于第i层//先计算u,y值,用u值替换a值,用y值替代b值for(int k=i;k<=n+1;k++){for(int r =1; r < i; r++)
A[i][k]= A[i][k]- A[i][r]*A[r][k];}//再计算l值,替代下三角的a值for(int k = i +1; k <= n; k++){for(int r =1; r<i; r++){
A[k][i]= A[k][i]- A[k][r]*A[r][i];}
A[k][i]/= A[k][k];}}
替换后的值如下:
L
&
U
&
y
=
(
u
11
u
12
u
13
⋯
u
1
n
∣
y
1
l
21
u
22
u
23
⋯
u
2
n
∣
y
2
l
31
l
32
u
33
⋯
u
3
n
∣
y
3
⋮
⋮
⋮
⋮
⋮
l
n
1
l
n
2
⋯
⋯
u
n
n
∣
y
n
)
\boldsymbol{L\&U\&y}=\left(\begin{array}{ccccc} u_{11} & u_{12} & u_{13} & \cdots & u_{1n} & |&y_1\\ l_{21} & u_{22} & u_{23} & \cdots & u_{2n} & |&y_2\\ l_{31} & l_{32} & u_{33} & \cdots & u_{3n} & |&y_3\\ \vdots & \vdots & \vdots & & \vdots & &\vdots\\ l_{n 1} & l_{n 2} & \cdots & \cdots & u_{nn} & |&y_n \end{array}\right)
L&U&y=⎝⎜⎜⎜⎜⎜⎛u11l21l31⋮ln1u12u22l32⋮ln2u13u23u33⋮⋯⋯⋯⋯⋯u1nu2nu3n⋮unn∣∣∣∣y1y2y3⋮yn⎠⎟⎟⎟⎟⎟⎞
最后通过计算
U
x
=
y
Ux = y
Ux=y得出
x
x
x。用
x
x
x的值替代
y
y
y值。 有
f
o
r
(
k
=
i
+
1
;
k
<
n
;
k
+
+
)
for(k=i+1;k<n;k++)
for(k=i+1;k<n;k++)
A
[
i
]
[
n
+
1
]
=
A
[
i
]
[
n
+
1
]
−
A
[
i
]
[
k
]
∗
A
[
k
]
[
n
+
1
]
A[i][n+1] = A[i][n+1] - A[i][k]*A[k][n+1]
A[i][n+1]=A[i][n+1]−A[i][k]∗A[k][n+1]
A
[
i
]
[
n
+
1
]
/
=
A
[
i
]
[
i
]
A[i][n+1] /=A[i][i]
A[i][n+1]/=A[i][i] 即 :
黑
色
=
黑
色
−
红
色
的
积
和
黑色 = 黑色 - 红色的积和
黑色=黑色−红色的积和 再计算:
黑
色
=
黑
色
/
橘
色
黑色 = 黑色 / 橘色
黑色=黑色/橘色
计算量
也是
1
3
n
3
\frac{1}{3}n^3
31n3数量级,与高斯消元法计算量基本相同。
优点
若是有多个系数矩阵相同而右边项不同的一系列方程组,用直接三角分解法更简单。
三、Cholesky分解法(平方根法)
条件
系数矩阵
A
A
A对称,即
A
=
A
T
A = A^T
A=AT
系数矩阵
A
A
A正定,即
A
A
A的特征值均为正值。
分解结果
可以将
A
A
A唯一的分解为
A
=
L
L
T
A = LL^T
A=LLT
L
=
(
l
11
0
⋯
0
l
21
l
22
⋯
0
⋮
⋮
⋱
⋮
l
n
1
l
n
2
⋯
l
n
n
)
L=\left(\begin{array}{cccc} l_{11} & 0 & \cdots & 0 \\ l_{21} & l_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ l_{n 1} & l_{n 2} & \cdots & l_{n n} \end{array}\right)
L=⎝⎜⎜⎜⎛l11l21⋮ln10l22⋮ln2⋯⋯⋱⋯00⋮lnn⎠⎟⎟⎟⎞
计算L的步骤
从外层向内侧依次计算。同
L
U
LU
LU分解。
Step1,计算对角元素
l
k
k
=
a
k
k
−
∑
p
=
1
k
−
1
l
k
p
2
\begin{array}{c} l_{k k}=\sqrt{a_{k k}-\sum_{p=1}^{k-1} l_{k p}^{2}} \end{array}
lkk=akk−∑p=1k−1lkp2
l
i
k
=
(
a
i
k
−
∑
p
=
1
k
−
1
l
i
p
l
k
p
)
/
l
k
k
,
i
=
k
+
1
,
k
+
2
,
⋯
,
n
\begin{array}{c} l_{i k}=\left(a_{i k}-\sum_{p=1}^{k-1} l_{i p} l_{k p}\right) / l_{k k}, \\ i=k+1, k+2, \cdots, n \end{array}
lik=(aik−∑p=1k−1liplkp)/lkk,i=k+1,k+2,⋯,n
A
[
i
]
[
j
]
=
A
[
i
]
[
j
]
−
红
色
部
分
的
积
和
A[i][j]= A[i][j]- 红色部分的积和
A[i][j]=A[i][j]−红色部分的积和(注意:
因
为
A
13
并
未
更
新
为
l
13
,
所
以
要
将
A
[
1
]
[
3
]
替
换
为
A
[
3
]
[
1
]
这
是
与
U
L
分
解
不
同
的
地
方
因为A_{13}并未更新为l_{13},所以要将A[1][3]替换为A[3][1]这是与UL分解不同的地方
因为A13并未更新为l13,所以要将A[1][3]替换为A[3][1]这是与UL分解不同的地方)
A
[
i
]
[
j
]
=
A
[
i
]
[
j
]
/
对
角
元
素
(
橘
色
)
A[i][j] = A[i][j]/对角元素(橘色)
A[i][j]=A[i][j]/对角元素(橘色)
在完成 Cholesky 分解后, 我们可分别求解以下系数矩阵。即
U
=
L
T
U = L^T
U=LT:
L
Y
=
b
,
L
T
X
=
Y
\boldsymbol{L} \boldsymbol{Y}=\boldsymbol{b}, \quad \boldsymbol{L}^{\mathrm{T}} \boldsymbol{X}=\boldsymbol{Y}
LY=b,LTX=Y
{
y
1
=
b
1
y
i
=
b
i
−
∑
k
=
1
i
−
1
l
i
k
y
k
,
i
=
2
,
3
,
⋯
,
n
\left\{\begin{array}{l} y_{1}=b_{1} \\ y_{i}=b_{i}-\sum_{k=1}^{i-1} l_{i k} y_{k}, \quad i=2,3, \cdots, n \end{array}\right.
{y1=b1yi=bi−∑k=1i−1likyk,i=2,3,⋯,n
{
x
n
=
y
n
/
u
n
n
x
i
=
(
y
i
−
∑
k
=
i
+
1
n
u
i
k
x
k
)
/
u
i
i
,
i
=
n
−
1
,
n
−
2
,
⋯
,
1.
\left\{\begin{array}{l} x_{n}=y_{n} / u_{n n} \\ x_{i}=\left(y_{i}-\sum_{k=i+1}^{n} u_{i k} x_{k}\right) / u_{i i}, \quad i=n-1, n-2, \cdots, 1 . \end{array}\right.
{xn=yn/unnxi=(yi−∑k=i+1nuikxk)/uii,i=n−1,n−2,⋯,1.
由于上述线性方程组有大量开方运算,故改进Cholesky分解法。
改进的Cholesky分解法
将
A
A
A分解为:
A
=
L
D
L
T
A = LDL^T
A=LDLT
Step 1. 计算 D 值
d
j
=
a
j
j
−
∑
k
=
1
j
−
1
l
j
k
v
j
k
,
v
j
k
=
l
j
k
d
k
,
d_{j}=a_{j j}-\sum_{k=1}^{j-1} l_{j k} v_{j k}, \quad v_{j k}=l_{j k} d_{k},
dj=ajj−k=1∑j−1ljkvjk,vjk=ljkdk,
Step 2. 计算第j列的L值
注意:L为单位下三角矩阵。即:L矩阵的对角元素为1.
l
i
j
=
(
a
i
j
−
∑
k
=
1
j
−
1
l
i
k
v
j
k
)
/
d
j
,
i
=
j
+
1
,
j
+
2
,
⋯
,
n
.
\begin{array}{c} \\ l_{i j}=\left(a_{i j}-\sum_{k=1}^{j-1} l_{i k} v_{j k}\right) / d_{j}, \quad i=j+1, j+2, \cdots, n . \end{array}
lij=(aij−∑k=1j−1likvjk)/dj,i=j+1,j+2,⋯,n.
在完成分解后, 我们可分别求解下列系数矩。
L
Y
=
b
,
D
L
T
X
=
Y
\boldsymbol{L Y}=\boldsymbol{b}, \quad \boldsymbol{D} \boldsymbol{L}^{\mathrm{T}} \boldsymbol{X}=\boldsymbol{Y}
LY=b,DLTX=Y
追赶法
条件
三对角非线性方程组
对角占优(
∣
b
i
∣
>
=
∣
a
i
∣
+
∣
c
i
∣
|b_i|>=|a_i|+|c_i|
∣bi∣>=∣ai∣+∣ci∣)
(
b
1
c
1
a
2
b
2
c
2
⋱
⋱
⋱
a
i
b
i
c
i
⋱
⋱
⋱
a
n
−
1
b
n
−
1
c
n
−
1
a
n
b
n
)
(
x
1
x
2
⋮
x
i
⋮
x
n
−
1
x
n
)
=
(
d
1
d
2
⋮
d
i
⋮
d
n
−
1
d
n
)
\left(\begin{array}{ccccccc} b_{1} & c_{1} & & & & & \\ a_{2} & b_{2} & c_{2} & & & & \\ & \ddots & \ddots & \ddots & & & \\ & & a_{i} & b_{i} & c_{i} & & \\ & & & \ddots & \ddots & \ddots & \\ & & & & a_{n-1} & b_{n-1} & c_{n-1} \\ & & & & & a_{n} & b_{n} \end{array}\right)\left(\begin{array}{c} x_{1} \\ x_{2} \\ \vdots \\ x_{i} \\ \vdots \\ x_{n-1} \\ x_{n} \end{array}\right)=\left(\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{i} \\ \vdots \\ d_{n-1} \\ d_{n} \end{array}\right)
⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛b1a2c1b2⋱c2⋱ai⋱bi⋱ci⋱an−1⋱bn−1ancn−1bn⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛x1x2⋮xi⋮xn−1xn⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛d1d2⋮di⋮dn−1dn⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
回代求得方程组的解。
x
n
=
γ
n
,
x
i
=
γ
i
−
β
i
x
i
+
1
,
i
=
n
−
1
,
n
−
2
,
⋯
,
2
x_{n}=\gamma_{n}, \quad x_{i}=\gamma_{i}-\beta_{i} x_{i+1}, \quad i=n-1, n-2, \cdots, 2
xn=γn,xi=γi−βixi+1,i=n−1,n−2,⋯,2
误差分析
c
o
n
d
(
A
)
=
∣
∣
A
−
1
∣
∣
⋅
∣
∣
A
∣
∣
cond(A) = ||A^{-1}||\cdot||A||
cond(A)=∣∣A−1∣∣⋅∣∣A∣∣刻画了方程组
A
x
=
b
Ax=b
Ax=b的病态程度,及解对
A
,
b
A,b
A,b扰动的敏感程度。
性质
对可逆矩阵A,有:
c
o
n
d
(
A
)
=
c
o
n
d
(
A
−
1
)
cond(A) = cond(A^{-1})
cond(A)=cond(A−1)
c
o
n
d
(
k
A
)
=
c
o
n
d
(
A
)
cond(kA) = cond(A)
cond(kA)=cond(A)