工程数学 计算方法 第三章 线性方程组的数值解法
线性方程组的数值解法
迭代法
A x = b ⟹ x ‾ = B x ‾ + f ‾ ⟹ x ‾ ( k + 1 ) = B x ‾ ( k ) + f ‾ Ax=b\quad\Longrightarrow\quad\overline{x}=B\overline{x}+\overline{f}\quad\Longrightarrow\quad \overline{x}^{(k+1)}=B\overline{x}^{(k)}+\overline{f} Ax=b⟹x=Bx+f⟹x(k+1)=Bx(k)+f
问题:
B=?(迭代格式?)
能使吗?(收敛条件?)
好使吗?(收敛速度?误差程度?)
逐次逼近。从一个或多个初始量出发,按照一定的计算格式获得方程组数值解的方法。简单实用。
收敛性是迭代法的前提。
适用于高次稀疏矩阵(常为200阶以上)。
Jacobi迭代法
将AX=b改写为X=BX+f的形式,建立雅可比方法的迭代格式:
A
x
=
b
⟹
x
(
k
+
1
)
=
B
x
(
k
)
+
f
Ax=b\,\,\Longrightarrow\,\, x^{(k+1)}=Bx^{(k)}+f
Ax=b⟹x(k+1)=Bx(k)+f
其中,B称为迭代矩阵。
思想与不动点迭代类似,使用前需要判别收敛性。
问题:怎么用(迭代格式)?什么时候能用(收敛条件)?计算到什么程度(迭代停止条件)?
迭代格式
第n式留xn,其他移到等号右边,写为 xn=求和 的形式。
{
a
11
x
1
+
a
12
x
2
+
.
.
.
+
a
1
n
x
n
=
b
1
a
21
x
1
+
a
22
x
2
+
.
.
.
+
a
2
n
x
n
=
b
2
.
.
.
a
n
1
x
1
+
a
n
2
x
2
+
.
.
.
+
a
n
n
x
n
=
b
n
⟹
a
i
i
≠
0
{
x
1
=
1
a
11
(
−
a
12
x
2
−
.
.
.
−
a
1
n
x
n
+
b
1
)
x
2
=
1
a
22
(
−
a
21
x
1
−
.
.
.
−
a
2
n
x
n
+
b
2
)
.
.
.
x
n
=
1
a
n
n
(
−
a
n
1
x
1
−
.
.
.
−
a
n
(
n
−
1
)
x
n
−
1
+
b
n
)
A
x
=
b
⇔
(
D
+
L
+
U
)
x
=
b
⇔
D
x
=
−
(
L
+
U
)
x
+
b
⇔
x
=
−
D
−
1
(
L
+
U
)
x
+
D
−
1
b
⇔
x
=
B
x
+
f
\begin{cases} a_{11}x_1+a_{12}x_2+...+a_{1n}x_n=b_1\\ a_{21}x_1+a_{22}x_2+...+a_{2n}x_n=b_2\\ ...\\ a_{n1}x_1+a_{n2}x_2+...+a_{nn}x_n=b_n\\ \end{cases} \overset{a_{ii} \neq 0}{\Longrightarrow} \begin{cases} x_1=\frac{1}{a_{11}}(-a_{12}x_2-...-a_{1n}x_n+b_1)\\ x_2=\frac{1}{a_{22}}(-a_{21}x_1-...-a_{2n}x_n+b_2)\\ ...\\ x_n=\frac{1}{a_{nn}}(-a_{n1}x_1-...-a_{n(n-1)}x_{n-1}+b_n)\\ \end{cases} \\\,\\ \begin{aligned} Ax=b&\Leftrightarrow (D+L+U)x=b\\ &\Leftrightarrow Dx=-(L+U)x+b\\ &\Leftrightarrow x=-D^{-1}(L+U)x+D^{-1}b\\ &\Leftrightarrow x=Bx+f\\ \end{aligned}
⎩⎪⎪⎪⎨⎪⎪⎪⎧a11x1+a12x2+...+a1nxn=b1a21x1+a22x2+...+a2nxn=b2...an1x1+an2x2+...+annxn=bn⟹aii=0⎩⎪⎪⎪⎨⎪⎪⎪⎧x1=a111(−a12x2−...−a1nxn+b1)x2=a221(−a21x1−...−a2nxn+b2)...xn=ann1(−an1x1−...−an(n−1)xn−1+bn)Ax=b⇔(D+L+U)x=b⇔Dx=−(L+U)x+b⇔x=−D−1(L+U)x+D−1b⇔x=Bx+f
得Jacobi迭代公式:
X
(
k
+
1
)
=
−
D
−
1
(
L
+
U
)
X
(
k
)
+
D
−
1
b
,
k
=
0
,
1
,
.
.
.
X^{(k+1)}=-D^{-1}(L+U)X^{(k)}+D^{-1}b,\,k=0,1,...
X(k+1)=−D−1(L+U)X(k)+D−1b,k=0,1,...,其中
B
=
−
D
−
1
(
L
+
U
)
B=-D^{-1}(L+U)
B=−D−1(L+U)称为Jacobi迭代矩阵。
例:
{
a
11
x
1
+
a
12
x
2
+
a
13
x
3
=
b
1
a
21
x
1
+
a
22
x
2
+
a
23
x
3
=
b
2
a
31
x
1
+
a
32
x
2
+
a
33
x
3
=
b
3
⟹
a
i
i
≠
0
{
x
1
=
1
a
11
(
−
a
12
x
2
−
a
13
x
3
+
b
1
)
x
2
=
1
a
22
(
−
a
21
x
1
−
a
23
x
3
+
b
2
)
x
n
=
1
a
33
(
−
a
31
x
1
−
a
32
x
2
+
b
3
)
⟹
a
i
i
≠
0
{
x
1
(
k
+
1
)
=
1
a
11
(
−
a
12
x
2
(
k
)
−
a
13
x
3
(
k
)
+
b
1
)
x
2
(
k
+
1
)
=
1
a
22
(
−
a
21
x
1
(
k
)
−
a
23
x
3
(
k
)
+
b
2
)
x
n
(
k
+
1
)
=
1
a
33
(
−
a
31
x
1
(
k
)
−
a
32
x
2
(
k
)
+
b
3
)
\begin{cases} a_{11}x_1+a_{12}x_2+a_{13}x_3=b_1\\ a_{21}x_1+a_{22}x_2+a_{23}x_3=b_2\\ a_{31}x_1+a_{32}x_2+a_{33}x_3=b_3\\ \end{cases} \overset{a_{ii} \neq 0}{\Longrightarrow} \begin{cases} x_1=\frac{1}{a_{11}}(-a_{12}x_2-a_{13}x_3+b_1)\\ x_2=\frac{1}{a_{22}}(-a_{21}x_1-a_{23}x_3+b_2)\\ x_n=\frac{1}{a_{33}}(-a_{31}x_1-a_{32}x_{2}+b_3)\\ \end{cases}\\ \overset{a_{ii} \neq 0}{\Longrightarrow} \begin{cases} x_1^{(k+1)}=\frac{1}{a_{11}}(-a_{12}x_2^{(k)}-a_{13}x_3^{(k)}+b_1)\\ x_2^{(k+1)}=\frac{1}{a_{22}}(-a_{21}x_1^{(k)}-a_{23}x_3^{(k)}+b_2)\\ x_n^{(k+1)}=\frac{1}{a_{33}}(-a_{31}x_1^{(k)}-a_{32}x_2^{(k)}+b_3)\\ \end{cases}
⎩⎪⎨⎪⎧a11x1+a12x2+a13x3=b1a21x1+a22x2+a23x3=b2a31x1+a32x2+a33x3=b3⟹aii=0⎩⎪⎨⎪⎧x1=a111(−a12x2−a13x3+b1)x2=a221(−a21x1−a23x3+b2)xn=a331(−a31x1−a32x2+b3)⟹aii=0⎩⎪⎨⎪⎧x1(k+1)=a111(−a12x2(k)−a13x3(k)+b1)x2(k+1)=a221(−a21x1(k)−a23x3(k)+b2)xn(k+1)=a331(−a31x1(k)−a32x2(k)+b3)
收敛条件⭐
误差:
e
→
(
k
)
=
B
k
e
→
(
0
)
\overrightarrow e^{(k)}=B^k\overrightarrow e^{(0)}
e(k)=Bke(0)
得收敛的充分条件
∥
B
∥
<
1
\Vert B\Vert < 1
∥B∥<1。
向量: x \textbf{x} x是否收敛。向量收敛:分量各自收敛。
影响收敛性的因素:初值和方法本身。
Jacobi迭代法收敛三种判断方法:
-
A行对角占优;
-
充分条件: ∥ B ∥ < 1 \Vert B\Vert < 1 ∥B∥<1;
-
⭐充要条件: B k → 0 ⟺ ρ ( B ) < 1 B^k\rightarrow0\Longleftrightarrow \rho(B)<1 Bk→0⟺ρ(B)<1。含参判断收敛性必须用特征值来算。
对角占优:矩阵 A = ( a i j ) ∈ R n × n A=(a_{ij})\in \textbf{R}^{n\times n} A=(aij)∈Rn×n,若其满足 ∣ a i j ∣ ⩾ ∑ j = 1 , j ≠ i n ∣ a i j ∣ , i = 1 , 2 , . . . , n , |a_{ij}|\geqslant \sum_{j=1,j\neq i}^n |a_{ij}|,\,i=1,2,...,n, ∣aij∣⩾∑j=1,j=in∣aij∣,i=1,2,...,n,则称其对角占优。(行对角占优,即对角线上的数值比同行其他数绝对值之和大)
停止条件
一般精度要求为
∣
∣
X
−
X
(
k
)
∣
∣
i
⩽
ε
||\textbf{X}-\textbf{X}^{(k)}||_i\leqslant \varepsilon
∣∣X−X(k)∣∣i⩽ε。因为
R
n
R^n
Rn上范数等价,而且1-范数计算简单且更为严格,故实际计算中常取1-范数。且实际中无法得知精确值,故一般取两步之间的插值,即当
∣
∣
X
(
k
+
1
)
−
X
(
k
)
∣
∣
1
⩽
ε
||\textbf{X}^{(k+1)}-\textbf{X}^{(k)}||_1\leqslant \varepsilon
∣∣X(k+1)−X(k)∣∣1⩽ε
时停止迭代。
Gauss-Seidel迭代法
对于Jacobi迭代法,每一步计算的时候都已经计算出新值,但靠后的分量仍然是用旧值计算的。优化:将式子中的值都使用最新值代入迭代,第i式中将
x
1
x_1
x1至
x
i
−
1
x_i-1
xi−1都使用第k+1次的迭代值代入计算。即:
x
i
(
k
+
1
)
=
1
a
i
i
(
−
a
i
1
x
1
(
k
+
1
)
−
.
.
.
−
a
i
(
i
−
1
)
x
i
−
1
(
k
+
1
)
−
a
i
(
i
+
1
)
x
i
+
1
(
k
)
−
.
.
.
−
a
i
n
x
n
(
k
)
+
b
i
)
=
1
a
i
i
(
b
i
−
∑
j
=
1
i
−
1
a
i
j
x
j
(
k
+
1
)
−
∑
j
=
i
+
1
n
a
i
j
x
j
(
k
)
)
(
i
=
1
,
2
,
.
.
.
,
n
)
\begin{array}{c} x_i^{(k+1)}&=\frac{1}{a_{ii}}(-a_{i1}x_1^{(k+1)}-...-a_{i(i-1)}x_{i-1}^{(k+1)}-a_{i(i+1)}x_{i+1}^{(k)}-...-a_{in}x_{n}^{(k)}+b_i)\\ &=\frac{1}{a_{ii}}(b_i -\sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)} -\sum_{j=i+1}^na_{ij}x_j^{(k)})(i=1,2,...,n) \end{array}
xi(k+1)=aii1(−ai1x1(k+1)−...−ai(i−1)xi−1(k+1)−ai(i+1)xi+1(k)−...−ainxn(k)+bi)=aii1(bi−∑j=1i−1aijxj(k+1)−∑j=i+1naijxj(k))(i=1,2,...,n)
矩阵形式:
x
(
k
+
1
)
=
−
D
−
1
(
L
x
(
k
+
1
)
+
U
x
(
k
)
)
+
D
−
1
b
⇔
(
D
+
L
)
x
(
k
+
1
)
=
−
U
x
(
k
)
+
b
⇔
x
(
k
+
1
)
=
−
(
D
+
L
)
−
1
U
x
(
k
)
+
(
D
+
L
)
−
1
b
⇔
x
(
k
+
1
)
=
B
x
(
k
)
+
f
\begin{aligned} x^{(k+1)}=-D^{-1}(Lx^{(k+1)}+Ux^{(k)})+D^{-1}b &\Leftrightarrow (D+L)x^{(k+1)}=-Ux^{(k)}+b\\ &\Leftrightarrow x^{(k+1)} = -(D+L)^{-1}Ux^{(k)}+(D+L)^{-1}b\\ &\Leftrightarrow x^{(k+1)}=Bx^{(k)}+f\\ \end{aligned}
x(k+1)=−D−1(Lx(k+1)+Ux(k))+D−1b⇔(D+L)x(k+1)=−Ux(k)+b⇔x(k+1)=−(D+L)−1Ux(k)+(D+L)−1b⇔x(k+1)=Bx(k)+f
其中 B = − ( D + L ) − 1 U B=-(D+L)^{-1}U B=−(D+L)−1U称为Gauss-Seidel迭代矩阵。
与Jacobi迭代相比,这里每一次计算都是用的是最新计算得的数值。
可以原位运算。
收敛条件:
- A行对角占优;
- A为正定矩阵;
- 充分条件: ∥ B ∥ < 1 \Vert B\Vert < 1 ∥B∥<1;
- 充要条件: ρ ( B ) < 1 \rho(B)<1 ρ(B)<1。
判断正定矩阵常用条件:一切顺序主子式为正。i阶顺序主子式:取左上方i×i阶矩阵,这个矩阵对应的方程即为i阶顺序主子式。
精度判断与Jacobi迭代法相同。
SOR方法
逐次超松弛迭代法(Successive Over Relaxation Method),课看作带参数ω的Gauss-Seidel迭代法。是对G-S迭代法的修正或加速。
算法构造
x ‾ i ( k + 1 ) = 1 a i i ( b i − ∑ j = 1 i − 1 a i j x j ( k + 1 ) − ∑ j = i + 1 n a i j x j ( k ) ) ( i = 1 , 2 , . . . , n ) x i ( k + 1 ) = ( 1 − ω ) x i ( k ) + ω x ‾ i ( k + 1 ) = x i ( k ) + ω ( x ‾ i ( k + 1 ) − x i ( k ) ) \overline{x}_i^{(k+1)}=\frac{1}{a_{ii}}(b_i -\sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)} -\sum_{j=i+1}^na_{ij}x_j^{(k)})(i=1,2,...,n)\\\,\\ \begin{array}{c} x_i^{(k+1)} &=(1-\omega)x_i^{(k)}+\omega \overline{x}_i^{(k+1)}\\ &=x_i^{(k)}+\omega (\overline{x}_i^{(k+1)}-x_i^{(k)}) \end{array} xi(k+1)=aii1(bi−j=1∑i−1aijxj(k+1)−j=i+1∑naijxj(k))(i=1,2,...,n)xi(k+1)=(1−ω)xi(k)+ωxi(k+1)=xi(k)+ω(xi(k+1)−xi(k))
ω称为松弛因子。
由G-S迭代得修正前的值 x ‾ i ( k + 1 ) \overline{x}_i^{(k+1)} xi(k+1),并通过这个式子来修正,得到迭代值 x i ( k + 1 ) x_i^{(k+1)} xi(k+1)。
0<ω<1 :低松弛法
ω=1 :G-S迭代法
1<ω<2 :超松弛法
SOR迭代法的计算(不考这个计算)
(每个式子都是通式,代表一个方程组)
先移项:
x
i
=
1
a
i
i
(
−
a
i
1
x
1
−
.
.
.
−
a
i
(
i
−
1
)
x
i
−
1
−
a
i
(
i
+
1
)
x
i
+
1
−
.
.
.
−
a
i
n
x
n
+
b
i
)
x_i=\frac{1}{a_{ii}}(-a_{i1}x_1-...-a_{i(i-1)}x_{i-1}-a_{i(i+1)}x_{i+1}-...-a_{in}x_{n}+b_i)
xi=aii1(−ai1x1−...−ai(i−1)xi−1−ai(i+1)xi+1−...−ainxn+bi)
写G-S迭代公式:
x
i
(
k
+
1
)
=
1
a
i
i
(
−
a
i
1
x
1
(
k
+
1
)
−
.
.
.
−
a
i
(
i
−
1
)
x
i
−
1
(
k
+
1
)
−
a
i
(
i
+
1
)
x
i
+
1
(
k
)
−
.
.
.
−
a
i
n
x
n
(
k
)
+
b
i
)
x_i^{(k+1)}=\frac{1}{a_{ii}}(-a_{i1}x_1^{(k+1)}-...-a_{i(i-1)}x_{i-1}^{(k+1)}-a_{i(i+1)}x_{i+1}^{(k)}-...-a_{in}x_{n}^{(k)}+b_i)
xi(k+1)=aii1(−ai1x1(k+1)−...−ai(i−1)xi−1(k+1)−ai(i+1)xi+1(k)−...−ainxn(k)+bi)
改写为SOR公式:
x
i
(
k
+
1
)
=
(
1
−
ω
)
x
i
(
k
)
+
ω
1
a
i
i
(
−
a
i
1
x
1
(
k
+
1
)
−
.
.
.
−
a
i
(
i
−
1
)
x
i
−
1
(
k
+
1
)
−
a
i
(
i
+
1
)
x
i
+
1
(
k
)
−
.
.
.
−
a
i
n
x
n
(
k
)
+
b
i
)
即
x
i
(
k
+
1
)
=
(
1
−
ω
)
x
i
(
k
)
+
ω
x
‾
i
(
k
+
1
)
x_i^{(k+1)}=(1-\omega)x_i^{(k)}+\omega\frac{1}{a_{ii}}(-a_{i1}x_1^{(k+1)}-...-a_{i(i-1)}x_{i-1}^{(k+1)}-a_{i(i+1)}x_{i+1}^{(k)}-...-a_{in}x_{n}^{(k)}+b_i)\\ 即\,\,x_i^{(k+1)}=(1-\omega)x_i^{(k)}+\omega \overline{x}_i^{(k+1)}
xi(k+1)=(1−ω)xi(k)+ωaii1(−ai1x1(k+1)−...−ai(i−1)xi−1(k+1)−ai(i+1)xi+1(k)−...−ainxn(k)+bi)即xi(k+1)=(1−ω)xi(k)+ωxi(k+1)
选取ω并带入初值。
直接法
直接法是理论上可以求得精确解的方法。
Guess消元法
- 消元
- 回代
矩阵:
- 化上三角矩阵:前推过程
主对角线上有0:换行/换列, D i ≠ 0 D_i\neq 0 Di=0时可换行 - 自下而上解出自变量:回代过程
问题:化上三角矩阵时若有小主元会导致舍入误差的扩散,从而导致计算失败。
解决:全主元消去法和列主元消去法。全主元消去法更加稳定。
Gauss-Jordan消去法:主元变1,上下消为0。常用与求逆矩阵,不用于解方程。 [ A ∣ I ] ⇒ [ I ∣ A − 1 ] [A|I]\Rightarrow[I|A^{-1}] [A∣I]⇒[I∣A−1]。
Doolittle分解法⭐
有时系数矩阵是确定的,而对于不同的问题只是b不同,而Guess消元法需要重新计算。
此外前推过程计算繁琐且易造成误差,而回代过程较为简单:能不能只回代不前推?或者增加回代的过程减少前推的过程:用回代过程代替前推过程。(回代过程的条件:A化为上三角矩阵/下三角矩阵。)
A
x
→
=
b
→
A\overrightarrow{x}=\overrightarrow{b}
Ax=b,A表示
x
i
x_i
xi之间的线性组合,表现了
x
i
x_i
xi间的相关性。若A定,分解成
L
U
X
→
=
b
→
LU\overrightarrow{X}=\overrightarrow{b}
LUX=b,其中LU为三角阵。
L
U
X
→
=
b
→
⇒
{
L
y
→
=
b
→
U
x
→
=
y
→
LU\overrightarrow{X}=\overrightarrow{b} \,\,\Rightarrow\,\, \begin{cases} L \overrightarrow{y}=\overrightarrow{b}\\ U \overrightarrow{x}=\overrightarrow{y} \end{cases}
LUX=b⇒{Ly=bUx=y
分解然后二次迭代来避免前推过程。
可靠吗?
定理:A正定 ⇔ \Leftrightarrow ⇔A的顺序主子式均不为0 ⇒ \Rightarrow ⇒A的LU分解唯一
如何分解?
-
乘积递推:(麻烦,运算量浪费大,不这么写)
{ A ( n ) = L 1 L 2 . . . L n − 2 L n − 1 A ( 1 ) b ( n ) = L 1 L 2 . . . L n − 2 L n − 1 b ( 1 ) ⟹ { A = L 1 − 1 L 2 − 1 . . . L n − 2 − 1 L n − 1 − 1 A ( n ) b = L 1 − 1 L 2 − 1 . . . L n − 2 − 1 L n − 1 − 1 b ( n ) 则 L ≜ L 1 − 1 L 2 − 1 . . . L n − 2 − 1 L n − 1 − 1 为 下 三 角 矩 阵 , U = A ( n ) \begin{cases} A^{(n)}=L_1L_2...L_{n-2}L_{n-1}A^{(1)}\\ b^{(n)}=L_1L_2...L_{n-2}L_{n-1}b^{(1)} \end{cases} \Longrightarrow \begin{cases} A=L_1^{-1}L_2^{-1}...L_{n-2}^{-1}L_{n-1}^{-1}A^{(n)}\\ b=L_1^{-1}L_2^{-1}...L_{n-2}^{-1}L_{n-1}^{-1}b^{(n)} \end{cases}\\ 则L\triangleq L_1^{-1}L_2^{-1}...L_{n-2}^{-1}L_{n-1}^{-1}为下三角矩阵,U=A^{(n)} {A(n)=L1L2...Ln−2Ln−1A(1)b(n)=L1L2...Ln−2Ln−1b(1)⟹{A=L1−1L2−1...Ln−2−1Ln−1−1A(n)b=L1−1L2−1...Ln−2−1Ln−1−1b(n)则L≜L1−1L2−1...Ln−2−1Ln−1−1为下三角矩阵,U=A(n) -
比较法
根据A=LU导出L和U:
[ a 11 ⋯ ⋯ a 1 n ⋮ ⋮ ⋮ ⋮ a n 1 ⋯ ⋯ a n n ] = [ 1 l 21 1 ⋮ ⋯ l n 1 ⋯ ⋯ 1 ] [ u 11 ⋯ ⋯ u 1 n ⋯ ⋮ ⋮ u n n ] { u 1 i = a 1 i , i = 1 , 2 , . . . , n , l i 1 = a i 1 u 11 , i = 2 , 3 , . . . , n , u r i = a r i − ∑ k = 1 r − 1 l r k u k i i = r , r + 1 , . . . , n , r = 2 , 3 , . . . , n , l r i = ( a i r − ∑ k = 1 r − 1 l i k u k r ) / u r r i = r , r + 1 , . . . , n , r = 2 , 3 , . . . , n , \left[\begin{matrix} a_{11} &\cdots &\cdots &a_{1n}\\ \vdots & & &\vdots\\ \vdots & & &\vdots\\ a_{n1} &\cdots &\cdots &a_{nn}\\ \end{matrix}\right]= \left[\begin{matrix} 1 \\ l_{21} &1 \\ \vdots &\cdots \\ l_{n1} &\cdots &\cdots &1 \\ \end{matrix}\right] \left[\begin{matrix} u_{11} &\cdots &\cdots &u_{1n}\\ & &\cdots &\vdots\\ & & &\vdots\\ & & &u_{nn}\\ \end{matrix}\right]\\ \begin{cases} u_{1i}=a_{1i}, &i=1,2,...,n,\\ l_{i1}=\frac{a_{i1}}{u_{11}}, &i=2,3,...,n,\\ u_{ri}=a_{ri}-\sum_{k=1}^{r-1}l_{rk}u_{ki} &i=r,r+1,...,n,\,r=2,3,...,n,\\ l_{ri}=(a_{ir}-\sum_{k=1}^{r-1}l_{ik}u_{kr})/u_{rr} &i=r,r+1,...,n,\,r=2,3,...,n,\\ \end{cases} ⎣⎢⎢⎢⎢⎡a11⋮⋮an1⋯⋯⋯⋯a1n⋮⋮ann⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎡1l21⋮ln11⋯⋯⋯1⎦⎥⎥⎥⎤⎣⎢⎢⎢⎢⎡u11⋯⋯⋯u1n⋮⋮unn⎦⎥⎥⎥⎥⎤⎩⎪⎪⎪⎨⎪⎪⎪⎧u1i=a1i,li1=u11ai1,uri=ari−∑k=1r−1lrkukilri=(air−∑k=1r−1likukr)/urri=1,2,...,n,i=2,3,...,n,i=r,r+1,...,n,r=2,3,...,n,i=r,r+1,...,n,r=2,3,...,n,
u第一行,l第一列,u第二行,l第二列……(公式只是明确一点计算方法,实际计算之间按照这个顺序就能写出来了,写过程的时候也不需要写出分解的过程,直接得到UL就可以了)可以使用原位运算: a i r a_{ir} air只使用一次,用完后就不会在用到了,可把UL存在A的原地址中,节约空间。
扰动分析
实际求解中可能会有输入输出的微小扰动。不同方程组对扰动的敏感程度不同。分析对扰动的敏感程度:
定理:设线性方程组系数矩阵A及右端向量b分别带有微小扰动
δ
A
\boldsymbol{\delta_A}
δA和
δ
b
\boldsymbol{\delta_b}
δb,其导出的扰动解
X
+
δ
X
\boldsymbol{X}+\boldsymbol{\delta_X}
X+δX满足
(
A
+
δ
A
)
(
X
+
δ
X
)
=
b
+
δ
b
(\boldsymbol{A}+\boldsymbol{\delta_A}) (\boldsymbol{X}+\boldsymbol{\delta_X})=\boldsymbol{b}+\boldsymbol{\delta_b}
(A+δA)(X+δX)=b+δb
则其解具有如下局部相对误差估计:
∥
δ
X
∥
∥
X
∥
⩽
∥
A
−
1
∥
∥
A
∥
1
−
∥
A
−
1
∥
∥
δ
A
∥
(
∥
δ
b
∥
∥
b
∥
+
∥
δ
A
∥
∥
A
∥
)
\frac{\Vert\boldsymbol{\delta_X}\Vert}{\Vert\boldsymbol{X}\Vert} \leqslant \frac{\Vert\boldsymbol{A}^{-1}\Vert\Vert\boldsymbol{A}\Vert} {1-\Vert\boldsymbol{A}^{-1}\Vert\Vert\boldsymbol{\delta_A}\Vert} \left(\frac{\Vert\boldsymbol{\delta_b}\Vert}{\Vert\boldsymbol{b}\Vert} +\frac{\Vert\boldsymbol{\delta_A}\Vert}{\Vert\boldsymbol{A}\Vert}\right)
∥X∥∥δX∥⩽1−∥A−1∥∥δA∥∥A−1∥∥A∥(∥b∥∥δb∥+∥A∥∥δA∥)
当扰动
δ
A
\boldsymbol{\delta_A}
δA充分小时,有近似关系:
∥
A
−
1
∥
∥
A
∥
1
−
∥
A
−
1
∥
∥
δ
A
∥
≈
∥
A
−
1
∥
∥
A
∥
⟹
∥
δ
X
∥
∥
X
+
δ
X
∥
⩽
∥
A
−
1
∥
∥
δ
A
∥
=
∥
A
−
1
∥
∥
A
∥
∥
δ
A
∥
∥
A
∥
∴
∥
δ
A
∥
∥
A
∥
被
放
大
了
∥
A
−
1
∥
∥
A
∥
倍
。
\frac{\Vert\boldsymbol{A}^{-1}\Vert\Vert\boldsymbol{A}\Vert} {1-\Vert\boldsymbol{A}^{-1}\Vert\Vert\boldsymbol{\delta_A}\Vert} \approx \Vert\boldsymbol{A}^{-1}\Vert\Vert\boldsymbol{A}\Vert\\\,\\ \Longrightarrow\,\, \frac{\Vert\boldsymbol{\delta_X}\Vert}{\Vert\boldsymbol{X}+\boldsymbol{\delta_X}\Vert} \leqslant \Vert\boldsymbol{A}^{-1}\Vert\Vert\boldsymbol{\delta_A}\Vert = \Vert\boldsymbol{A}^{-1}\Vert\Vert\boldsymbol{A}\Vert \frac{\Vert\boldsymbol{\delta_A}\Vert}{\Vert\boldsymbol{A}\Vert}\\\,\\ \therefore \frac{\Vert\boldsymbol{\delta_A}\Vert}{\Vert\boldsymbol{A}\Vert} 被放大了\Vert\boldsymbol{A}^{-1}\Vert\Vert\boldsymbol{A}\Vert倍。
1−∥A−1∥∥δA∥∥A−1∥∥A∥≈∥A−1∥∥A∥⟹∥X+δX∥∥δX∥⩽∥A−1∥∥δA∥=∥A−1∥∥A∥∥A∥∥δA∥∴∥A∥∥δA∥被放大了∥A−1∥∥A∥倍。
定义:设A为非奇异矩阵(满秩 / |A|=0),称数
C
o
n
d
(
A
)
v
=
∥
A
−
1
∥
v
∥
A
∥
v
(
v
=
1
,
2
或
∞
)
Cond(A)_v=\Vert\boldsymbol{A}^{-1}\Vert_v\Vert\boldsymbol{A}\Vert_v(v=1,2或\infty)
Cond(A)v=∥A−1∥v∥A∥v(v=1,2或∞)为矩阵A的条件数。
C
o
n
d
1
(
A
)
≜
∥
A
−
1
∥
1
∥
A
∥
1
C
o
n
d
∞
(
A
)
≜
∥
A
−
1
∥
∞
∥
A
∥
∞
C
o
n
d
2
(
A
)
≜
∥
A
−
1
∥
2
∥
A
∥
2
=
max
1
⩽
i
⩽
n
{
λ
i
A
T
A
}
min
1
⩽
i
⩽
n
{
λ
i
A
T
A
}
若
A
为
对
称
矩
阵
,
C
o
n
d
2
(
A
)
=
max
∣
λ
A
∣
min
∣
λ
A
∣
Cond_1(A)\triangleq\Vert\boldsymbol{A}^{-1}\Vert_1\Vert\boldsymbol{A}\Vert_1\\\,\\ Cond_\infty(A)\triangleq\Vert\boldsymbol{A}^{-1}\Vert_\infty\Vert\boldsymbol{A}\Vert_\infty\\\,\\ Cond_2(A)\triangleq\Vert\boldsymbol{A}^{-1}\Vert_2\Vert\boldsymbol{A}\Vert_2= \sqrt{\frac{\underset{1\leqslant i\leqslant n}{\max}\{\lambda_i^{\boldsymbol{A}^T\boldsymbol{A}}\}} {\underset{1\leqslant i\leqslant n}{\min}\{\lambda_i^{\boldsymbol{A}^T\boldsymbol{A}}\}}}\\\,\\ 若A为对称矩阵, Cond_2(A)=\sqrt{\frac{\max\vert\lambda^A\vert}{\min\vert\lambda^A\vert}}
Cond1(A)≜∥A−1∥1∥A∥1Cond∞(A)≜∥A−1∥∞∥A∥∞Cond2(A)≜∥A−1∥2∥A∥2=1⩽i⩽nmin{λiATA}1⩽i⩽nmax{λiATA}若A为对称矩阵,Cond2(A)=min∣λA∣max∣λA∣
当
C
o
n
d
(
A
)
>
>
1
Cond(A)>>1
Cond(A)>>1时,视为方程组是病态的,即当输入有微小偏差时扰动解
X
^
\hat{\boldsymbol{X}}
X^存在较大误差。
病态与否为方程组的性质,病态矩阵没有任何补救方法,只能修改计算方法。(这一部分记这一句话就行了)
应试
迭代法的收敛性
收敛性完全取决于迭代矩阵的性质,与初值选取无关。
考试可能只考二阶,三阶不好求逆矩阵(线代忘光)。
迭代的求解计算量大,不考,只计算收敛性。记住迭代公式的形式,记矩阵形式即可,然后用
ρ
(
B
)
<
1
\rho(B)<1
ρ(B)<1 来判断收敛性。
ρ
(
B
)
\rho(B)
ρ(B) :
B
B
B 特征值绝对值的最大值。
迭代格式里要求逆矩阵:
二阶矩阵:主交换,副相反,除以行列式
三阶矩阵:
- 初等变换法
- 伴随式除以行列式
伴随式:伴随阵中 a i j a_{ij} aij为原矩阵第 j j j行和第 i i i列划掉,剩下的部分求行列式的值,并乘上 ( − 1 ) i + j (-1)^{i+j} (−1)i+j。例: a 22 = ( − 1 ) 4 ∣ a 11 a 13 a 31 a 33 ∣ a_{22}=(-1)^4\left\vert\begin{matrix}a_{11}&a_{13}\\a_{31}&a_{33}\end{matrix}\right\vert a22=(−1)4∣∣∣∣a11a31a13a33∣∣∣∣ 。
求
ρ
(
B
)
\rho(B)
ρ(B):求矩阵特征值:
解方程:
∣
λ
I
−
B
∣
=
0
\vert\lambda I-B\vert=0
∣λI−B∣=0 。若为三阶矩阵注意因式分解,先找一个因式才好解。
Doolittle分解法
考试考三阶的。
分解方法:u第一行,l第一列,u第二行,l第二列,u第三行,l第三列。
可能会考条件数的计算,和矩阵范数一起考。概念在上文扰动分析部分。