https://zhuanlan.zhihu.com/p/30485749
设 n n n阶线性方程组:
{ 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 m 1 x 1 + a m 2 x 2 + . . . + a m n x n = b n ( 2 − 1 ) \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\\ \dots \dots \dots\dots\dots\dots\\ a_{m1}x_1+a_{m2}x_2+...+a_{mn}x_n=b_n\\ \end {cases}\quad\quad\quad\quad\quad(2-1) ⎩⎪⎪⎪⎨⎪⎪⎪⎧a11x1+a12x2+...+a1nxn=b1a21x1+a22x2+...+a2nxn=b2………………am1x1+am2x2+...+amnxn=bn(2−1)
其矩阵形式为:
A x = b ( 2 − 2 ) Ax=b\quad\quad\quad\quad\quad\quad(2-2) Ax=b(2−2)
如果线性方程组(2-1)的系数矩阵 A A A的行列式不为0,即 det A ≠ 0 \det A\ne0 detA=0,则该方程组有唯一解。由克拉默法则,其解为:
x i = det A i det A ( i = 1 , 2 , … , n ) x_i=\frac{\det A_i}{\det A}\quad(i=1,2,\dots,n) xi=detAdetAi(i=1,2,…,n)
其中 A i A_i Ai为用方程组(2-1)的右端向量 b b b代替 A A A中第 i i i列向量得到的矩阵。这种方法需要计算n+1个n阶行列式并作n次除法,而每个n阶行列式计算需做 ( n − 1 ) × n ! (n-1)\times n! (n−1)×n!次乘法,计算量十分惊人。因此当n较大时,在实际计算中是不可行的。
解线性方程组的方法大致可分为两类:直接方法和迭代法。直接方法是指假设计算过程中不产生舍入误差,经过有限次运算可求得方程组的精确解的方法。迭代法是从解的某个近似值出发,通过构造一个无穷序列去逼近精确解的方法。一般地,有限步内得不到精确解。
在本章的介绍中,均假设方程组(2-1)的系数矩阵非奇异(?),方程组存在唯一解。
一、消去法(什么时候用小数什么时候用分数?)
1.1 Gauss消去法
解线性方程组的直接方法的基本思想是将方程组(2-1)变形为等价的三角形方程组,然后求解。Gauss消去法的基本思想和过程是:将方程组逐列逐行消去变量,转化为等价的上三角形方程组,此过程称为消元过程。然后按着方程组相反顺序求解上三角形方程组,得到原方程组的解,此过程称为回代过程。
为统一起见,将方程组(2-1)改写成:
{
a
11
(
1
)
x
1
+
a
12
(
1
)
x
2
+
.
.
.
+
a
1
n
(
1
)
x
n
=
b
1
(
1
)
a
21
(
1
)
x
1
+
a
22
(
1
)
x
2
+
.
.
.
+
a
2
n
(
1
)
x
n
=
b
2
(
1
)
…
…
…
…
…
…
a
m
1
(
1
)
x
1
+
a
m
2
(
1
)
x
2
+
.
.
.
+
a
m
n
(
1
)
x
n
=
b
n
(
1
)
(
2
−
3
)
\begin {cases} a_{11}^{(1)}x_1+a_{12}^{(1)}x_2+...+a_{1n}^{(1)}x_n=b_1^{(1)}\\ a_{21}^{(1)}x_1+a_{22}^{(1)}x_2+...+a_{2n}^{(1)}x_n=b_2^{(1)}\\ \dots \dots \dots\dots\dots\dots\\ a_{m1}^{(1)}x_1+a_{m2}^{(1)}x_2+...+a_{mn}^{(1)}x_n=b_n^{(1)}\\ \end {cases}\quad\quad\quad\quad\quad(2-3)
⎩⎪⎪⎪⎨⎪⎪⎪⎧a11(1)x1+a12(1)x2+...+a1n(1)xn=b1(1)a21(1)x1+a22(1)x2+...+a2n(1)xn=b2(1)………………am1(1)x1+am2(1)x2+...+amn(1)xn=bn(1)(2−3)
简记为 A ( 1 ) x = b ( 1 ) A^{(1)}x=b^{(1)} A(1)x=b(1),其中 A ( 1 ) = A , b ( 1 ) = b A^{(1)}=A,b^{(1)}=b A(1)=A,b(1)=b。
一般地,求解 n n n阶方程组(2-1)的Gauss消去法的步骤如下:
消元过程
第一步:设 a 11 ( 1 ) ≠ 0 a_{11}^{(1)}\ne0 a11(1)=0,记 l i 1 = a i 1 ( 1 ) / a 11 ( 1 ) l_{i1}=a_{i1}^{(1)}/a_{11}^{(1)} li1=ai1(1)/a11(1),将方程组(2-3)中第 i i i个方程减去第1个方程乘以 l i 1 ( i = 2 , 3 , … , n ) l_{i1}\ (i=2,3,\dots,n) li1 (i=2,3,…,n),完成第一次消元,得(2-3)的同解方程组:
待补充 16
1.2 主元素消去法
Gauss消去法虽然简单易行,但计算过程要求 a k k ( k ) a_{kk}^{(k)} akk(k)(称为主元素)均不为零,否则Gauss消去法不能进行。另一方面,有时主元素 a k k ( k ) ≠ 0 a_{kk}^{(k)}\ne0 akk(k)=0,但 ∣ a k k ( k ) ∣ |a_{kk}^{(k)}| ∣akk(k)∣很小,用其作除数,会使:
∣ l i k ( k ) ∣ = ∣ a i k ( k ) / a k k ( k ) ∣ |l_{ik}^{(k)}|=|a_{ik}^{(k)}/a_{kk}^{(k)}| ∣lik(k)∣=∣aik(k)/akk(k)∣很大,导致舍入误差的严重扩大以及计算结果的不正确。
待补充 19
1.2.1 列主元素法
为简便起见,用方程组(2-1)的增广矩阵:
[ A , b ] = [ a 11 … a 1 n b 1 ⋮ ⋮ ⋮ a n 1 … a n n b n ] ( 2 − 9 ) [A,b]=\begin{bmatrix}a_{11}&&\dots&&a_{1n}&&b_1\\\vdots&&&&\vdots&&\vdots\\a_{n1}&&\dots&&a_{nn}&&b_n\end{bmatrix}\quad\quad\quad\quad\quad(2-9) [A,b]=⎣⎢⎡a11⋮an1……a1n⋮annb1⋮bn⎦⎥⎤(2−9)
表示方程组,并直接在增广矩阵上运算。
列主元素法的基本思想是在每次消元前,在要消去未知数的系数中(即一列中)找到绝对值最大的系数作为主元,通过方程对换将其换到对角线上,然后进行消元,其具体步骤如下:
待补充 20
1.2.2 全主元素法
如果不是按列选主元,而是在全体待选系数 a i j ( k ) ( i , j = k , k + 1 , … , n ) a_{ij}^{(k)}(i,j=k,k+1,\dots,n) aij(k)(i,j=k,k+1,…,n)中选取主元,则得到全主元素法,其计算过程如下:
第一步:在全体系数 a i j ( i m j = 1 , … , n ) a_{ij}(imj=1,\dots,n) aij(imj=1,…,n)中选取绝对值最大的元作为主元,并通过行与列的互换把它换到 a 11 ( 1 ) a_{11}^{(1)} a11(1)的位置,然后进行第一次消元,得到矩阵 [ A ( 2 ) , b ( 2 ) ] [A^{(2)},b^{(2)}] [A(2),b(2)]。
第k步:在矩阵 A ( k ) A^{(k)} A(k)的右下方 n − k + 1 n-k+1 n−k+1阶子矩阵的所有元素 a i j ( k ) ( i , j = k , k + 1 , … , n ) a_{ij}^{(k)}(i,j=k,k+1,\dots,n) aij(k)(i,j=k,k+1,…,n)中,选取绝对值最大的元作为主元,并通过行与列的互换将它换到 a k k ( k ) a_{kk}^{(k)} akk(k)的位置,然后进行第k次消元。
经过n-1次消元后,得到与方程组(2-1)同解的上三角形方程组,再由回代过程求解。
待补充 21
二、直接三角分解法
2.1 Gauss消去法的矩阵形式
如果用增广矩阵 [ A , b ] [A,b] [A,b]表示线性方程组(2-1),Gauss消去法的消元过程可以通过用一串初等矩阵左乘增广矩阵表示。例如,第一次消元等价于用初等矩阵:
L 1 = [ 1 0 0 … 0 − l 21 1 0 … 0 − l 31 0 1 … 0 ⋮ ⋮ ⋮ ⋮ − l n 1 0 0 … 1 ] ( 2 − 10 ) L_1=\begin{bmatrix}1&&0&&0&&\dots&&0\\-l_{21}&&1&&0&&\dots&&0\\-l_{31}&&0&&1&&\dots&&0\\\vdots&&\vdots&&\vdots&&&&\vdots\\-l_{n1}&&0&&0&&\dots&&1\end{bmatrix}\quad\quad\quad\quad\quad(2-10) L1=⎣⎢⎢⎢⎢⎢⎡1−l21−l31⋮−ln1010⋮0001⋮0…………000⋮1⎦⎥⎥⎥⎥⎥⎤(2−10)
左乘矩阵 [ A ( 1 ) , b ( 1 ) ] = [ A , b ] [A^{(1)},b^{(1)}]=[A,b] [A(1),b(1)]=[A,b],其中 l i 1 = a i 1 ( 1 ) / a 11 ( 1 ) ( i = 2 , 3 , … , n ) l_{i1}=a_{i1}^{(1)}/a_{11}^{(1)}\ (i=2,3,\dots,n) li1=ai1(1)/a11(1) (i=2,3,…,n),即:
[ A ( 2 ) , b ( 2 ) ] = L 1 [ A ( 1 ) , b ( 1 ) ] [A^{(2)},b^{(2)}]=L_1[A^{(1)},b^{(1)}] [A(2),b(2)]=L1[A(1),b(1)]
一般地,第 k k k次消元等价于用初等矩阵:
L k = [ 1 ⋱ 1 − l k + 1 k 1 ⋮ ⋱ − l n k 1 ] ( 2 − 11 ) L_k=\begin{bmatrix}1&&&&&\\&\ddots&&&&\\&&1&&&\\&&-l_{k+1\ k}&1&&\\&&\vdots&&\ddots&\\&&-l_{nk}&&&1\end{bmatrix}\quad\quad\quad\quad(2-11) Lk=⎣⎢⎢⎢⎢⎢⎢⎢⎡1⋱1−lk+1 k⋮−lnk1⋱1⎦⎥⎥⎥⎥⎥⎥⎥⎤(2−11)
左乘矩阵 [ A ( k ) , b ( k ) ] [A^{(k)},b^{(k)}] [A(k),b(k)],其中 l i k = a i k ( k ) / a k k ( k ) ( i = k + 1 , … , n ) l_{ik}=a_{ik}^{(k)}/a_{kk}^{(k)}\ (i=k+1,\dots,n) lik=aik(k)/akk(k) (i=k+1,…,n),经过 n − 1 n-1 n−1次消元后得到:
待补充 23
2.2 矩阵的直接三角分解
Gauss消去法要求主元素 a k k ( k ) ≠ 0 ( k = 1 , … , n − 1 ) a_{kk}^{(k)}\ne0\ (k=1,\dots,n-1) akk(k)=0 (k=1,…,n−1),将这一条件与矩阵 A A A自身性质联系起来,有:
定理2.1
设 A A A为 n n n阶方阵,若 A A A的顺序主子式 A i ( i = 1 , 2 , … , n − 1 ) A_i\ (i=1,2,\dots,n-1) Ai (i=1,2,…,n−1)均不为零,则矩阵 A A A存在唯一的Doolittle分解。
下面讨论如何对 A A A进行LU分解。
两个矩阵相等就是它们的对应元素都相等,比较 A A A与LU的对应元素,即可得出L,U中元素的计算公式 A = L U A=LU A=LU,即:
[ a 11 a 12 a 13 … a 1 n a 21 a 22 a 23 … a 2 n a 31 a 32 a 33 … a 3 n ⋮ ⋮ ⋮ ⋮ a n 1 a n 2 a n 3 … a n n ] = [ 1 l 21 1 l 31 l 32 1 ⋮ ⋮ ⋮ ⋱ l n 1 l n 2 l n 3 … 1 ] [ u 11 u 12 … u 1 n u 22 … u 2 n ⋱ u n n ] \begin{bmatrix}a_{11}&&a_{12}&&a_{13}&&\dots&&a_{1n}\\a_{21}&&a_{22}&&a_{23}&&\dots&&a_{2n}\\a_{31}&&a_{32}&&a_{33}&&\dots&&a_{3n}\\\vdots&&\vdots&&\vdots&&&&\vdots\\a_{n1}&&a_{n2}&&a_{n3}&&\dots&&a_{nn}\end{bmatrix}=\begin{bmatrix}1&&&&\\l_{21}&1&&&\\l_{31}&l_{32}&1&&\\\vdots&\vdots&\vdots&\ddots&\\l_{n1}&l_{n2}&l_{n3}&\dots&1\end{bmatrix}\begin{bmatrix}u_{11}&u_{12}&\dots&u_{1n}\\&u_{22}&\dots&u_{2n}\\&&\ddots&\\&&&u_{nn}\end{bmatrix} ⎣⎢⎢⎢⎢⎢⎡a11a21a31⋮an1a12a22a32⋮an2a13a23a33⋮an3…………a1na2na3n⋮ann⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡1l21l31⋮ln11l32⋮ln21⋮ln3⋱…1⎦⎥⎥⎥⎥⎥⎤⎣⎢⎢⎡u11u12u22……⋱u1nu2nunn⎦⎥⎥⎤
由矩阵的乘法规则,得:
a
1
j
=
u
1
j
(
j
=
1
,
2
,
…
,
n
)
a_{1j}=u_{1j}\quad(j=1,2,\dots,n)
a1j=u1j(j=1,2,…,n)
a
i
j
=
{
∑
k
=
1
i
=
1
l
i
k
u
k
j
+
u
i
j
,
j
≥
i
∑
k
=
1
j
l
i
k
u
k
j
,
j
<
i
(
i
=
2
,
…
,
n
)
a_{ij}=\begin{cases}\sum\limits_{k=1}^{i=1}l_{ik}u_{kj}+u_{ij},\quad j\ge i\\\\\sum\limits_{k=1}^{j}l_{ik}u_{kj},\quad\quad\quad\ \ j<i\end{cases}\quad\quad\quad(i=2,\dots,n)
aij=⎩⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎧k=1∑i=1likukj+uij,j≥ik=1∑jlikukj, j<i(i=2,…,n)
由此可得计算 l i j l_{ij} lij和 u i j u_{ij} uij的公式:
{ u 1 j = a 1 j ( j = 1 , 2 , … , n ) u i j = a i j − ∑ k = 1 i − 1 l i k u k j ( i = 2 , … , n ; j = i , … , n ) l i j = ( a i j − ∑ k = 1 j − 1 l i k u k j ) / u j j ( j = 1 , 2 , … , n − 1 ; i = j + 1 , … , n ) ( 2 − 16 ) \begin{cases}u_{1j}=a_{1j}\quad\quad\quad\quad\quad\quad\quad\quad(j=1,2,\dots,n)\\\\u_{ij}=a_{ij}-\sum\limits_{k=1}^{i-1}l_{ik}u_{kj}\quad\quad\quad\ (i=2,\dots,n;\ j=i,\dots,n)\\\\l_{ij}=(a_{ij}-\sum\limits_{k=1}^{j-1}l_{ik}u_{kj})/u_{jj}\quad\ (j=1,2,\dots,n-1;\ i=j+1,\dots,n)\end{cases}\quad\quad\quad(2-16) ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧u1j=a1j(j=1,2,…,n)uij=aij−k=1∑i−1likukj (i=2,…,n; j=i,…,n)lij=(aij−k=1∑j−1likukj)/ujj (j=1,2,…,n−1; i=j+1,…,n)(2−16)
计算过程应按第一行,第一列,第二行,第二列, … … \dots\dots …… 的顺序进行,具体步骤如下:
- 计算
U
U
U的第一行,
L
L
L的第一列:
u 1 j = a 1 j ( j = 1 , 2 , … , n ) l i 1 = a i 1 / u 11 ( i = 2 , 3 , … , n ) u_{1j}=a_{1j}\quad\quad(j=1,2,\dots,n)\\l_{i1}=a_{i1}/u_{11}\quad\quad(i=2,3,\dots,n) u1j=a1j(j=1,2,…,n)li1=ai1/u11(i=2,3,…,n) - 计算
U
U
U的第
r
r
r行,
L
L
L的第
r
r
r列(
r
=
2
,
…
,
n
)
r=2,\dots,n)
r=2,…,n)
u r j = a r j − ∑ k = 1 r − 1 l r k u k j ( j = r , r + 1 , … , n ) l i r = ( a i r − ∑ k = 1 r − 1 l i k u k r ) / u r r ( i = r + 1 , … , n ) u_{rj}=a_{rj}-\sum_{k=1}^{r-1}l_{rk}u_{kj}\quad\quad\quad\quad(j=r,r+1,\dots,n)\\l_{ir}=(a_{ir}-\sum\limits_{k=1}^{r-1}l_{ik}u_{kr})/u_{rr}\quad\quad(i=r+1,\dots,n) urj=arj−k=1∑r−1lrkukj(j=r,r+1,…,n)lir=(air−k=1∑r−1likukr)/urr(i=r+1,…,n)
实际计算时,常将 u r j u_{rj} urj存放在 a r j a_{rj} arj的位置上,将 l i r l_{ir} lir存放在 a i r a_{ir} air的位置上,按表2-1所示的紧凑格式及顺序进行。
2.3 直接三角分解法
如果线性方程组 A x = b Ax=b Ax=b的系数矩阵已进行三角分解 A = L U A=LU A=LU,则解方程组 A x = b Ax=b Ax=b等价于求解两个三角形方程组 L y = b , U x = y Ly=b,Ux=y Ly=b,Ux=y,即由:
L y = [ 1 l 21 1 l 31 l 32 1 ⋮ ⋮ ⋮ ⋱ l n 1 l n 2 l n 3 … 1 ] [ y 1 y 2 ⋮ y n ] = [ b 1 b 2 ⋮ b n ] ( 2 − 17 ) Ly=\begin{bmatrix}1&&&&\\l_{21}&1&&&\\l_{31}&l_{32}&1&&\\\vdots&\vdots&\vdots&\ddots&\\l_{n1}&l_{n2}&l_{n3}&\dots&1\end{bmatrix}\begin{bmatrix}y_1\\y_2\\\vdots\\y_n\end{bmatrix}=\begin{bmatrix}b_1\\b_2\\\vdots\\b_n\end{bmatrix}\quad\quad\quad\quad(2-17) Ly=⎣⎢⎢⎢⎢⎢⎡1l21l31⋮ln11l32⋮ln21⋮ln3⋱…1⎦⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎡y1y2⋮yn⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡b1b2⋮bn⎦⎥⎥⎥⎤(2−17)
可求出
y k = { b 1 ( k = 1 ) b k − ∑ j = 1 k − 1 l k j y j ( k = 2 , … , n ) ( 2 − 18 ) y_k=\begin{cases}b_1\quad\quad\quad\quad\quad\quad\quad(k=1)\\b_k-\sum\limits_{j=1}^{k-1}l_{kj}y_j\quad\quad\ (k=2,\dots,n)\end{cases}\quad\quad\quad\quad\quad\quad(2-18) yk=⎩⎪⎨⎪⎧b1(k=1)bk−j=1∑k−1lkjyj (k=2,…,n)(2−18)
再由:
U x = [ u 11 u 12 … u 1 n u 22 … u 2 n ⋱ u n n ] [ x 1 x 2 ⋮ x n ] = [ y 1 y 2 ⋮ y n ] ( 2 − 19 ) Ux=\begin{bmatrix}u_{11}&u_{12}&\dots&u_{1n}\\&u_{22}&\dots&u_{2n}\\&&\ddots&\\&&&u_{nn}\end{bmatrix}\begin{bmatrix}x_1\\x_2\\\vdots\\x_n\end{bmatrix}=\begin{bmatrix}y_1\\y_2\\\vdots\\y_n\end{bmatrix}\quad\quad\quad\quad\quad(2-19) Ux=⎣⎢⎢⎡u11u12u22……⋱u1nu2nunn⎦⎥⎥⎤⎣⎢⎢⎢⎡x1x2⋮xn⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡y1y2⋮yn⎦⎥⎥⎥⎤(2−19)
解得:
待补充 28
三、特殊矩阵的三角分解法
如果矩阵 A A A事系数的三对角形矩阵或 A A A是对称正定阵,这使得 A A A的三角分解也具有更简洁的形式,从而导出一些特殊的算法。
3.1 解三对角方程组的追赶法
方程组:
[ b 1 c 1 a 2 b 2 c 2 ⋱ ⋱ ⋱ a n − 1 b n − 1 c n − 1 a n b n ] [ x 1 x 2 ⋮ x n − 1 x n ] = [ d 1 d 2 ⋮ d n − 1 d n ] ( 2 − 21 ) \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}\begin{bmatrix}x_1\\x_2\\\vdots\\x_{n-1}\\x_n\end{bmatrix}=\begin{bmatrix}d_1\\d_2\\\vdots\\d_{n-1}\\d_n\end{bmatrix}\quad\quad\quad\quad\quad(2-21) ⎣⎢⎢⎢⎢⎡b1a2c1b2⋱c2⋱an−1⋱bn−1ancn−1bn⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎡x1x2⋮xn−1xn⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡d1d2⋮dn−1dn⎦⎥⎥⎥⎥⎥⎤(2−21)
称为三对角方程组,其中系数矩阵 A A A是一种特殊的系数矩阵,它的非零元素集中分布在主对角线及其相邻的两条次对角线上,称为三对角矩阵。
如果 A A A满足三角分解的条件,可将 A A A分解成 L U LU LU,按Doolittle分解法可以证明如下定理。
定理2.3
设方程组(2-21)中的系数矩阵 A A A满足条件:
{ ∣ b 1 ∣ > ∣ c 1 ∣ > 0 ∣ b i ∣ ≥ ∣ a i ∣ + ∣ c i ∣ ( a i c i ≠ 0 , i = 2 , 3 , … , n − 1 ) ∣ b n ∣ > ∣ a n ∣ > 0 \begin{cases}|b_1|>|c_1|>0\\|b_i|\ge|a_i|+|c_i|\quad(a_ic_i\ne0,\ i=2,3,\dots,n-1)\\|b_n|>|a_n|>0\end{cases} ⎩⎪⎨⎪⎧∣b1∣>∣c1∣>0∣bi∣≥∣ai∣+∣ci∣(aici=0, i=2,3,…,n−1)∣bn∣>∣an∣>0
则 A A A可唯一分解为:
A = L U = [ 1 l 2 1 l 3 1 ⋱ ⋱ l n 1 ] [ u 1 c 1 u 2 c 2 u 3 c 3 ⋱ ⋱ ⋱ c n − 1 u n ] ( 2 − 22 ) A=LU=\begin{bmatrix}1&&&&\\l_2&1&&&\\&l_3&1&&\\&&\ddots&\ddots&\\&&& l_n&1\end{bmatrix}\begin{bmatrix}u_1&c_1&&&&\\&u_2&c_2&&&\\&&u_3&c_3&&\\&&&\ddots&\ddots&\\&&&&\ddots&c_{n-1}\\&&&&&u_n\end{bmatrix}\quad\quad\quad\quad\quad(2-22) A=LU=⎣⎢⎢⎢⎢⎡1l21l31⋱⋱ln1⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎡u1c1u2c2u3c3⋱⋱⋱cn−1un⎦⎥⎥⎥⎥⎥⎥⎤(2−22)
其中 c i ( i = 1 , 2 , … , n − 1 ) c_i(i=1,2,\dots,n-1) ci(i=1,2,…,n−1)由原矩阵给出,且:
{ u 1 = b 1 l i = a i / u i − 1 ( i = 2 , 3 , … , n ) u i = b i − c i − 1 l i ( 2 − 23 ) \begin{cases}u_1=b_1\\l_i=a_i/u_{i-1}\quad\quad(i=2,3,\dots,n)\\u_i=b_i-c_{i-1}l_i\end{cases}\quad\quad\quad\quad(2-23) ⎩⎪⎨⎪⎧u1=b1li=ai/ui−1(i=2,3,…,n)ui=bi−ci−1li(2−23)
当矩阵 A A A三角分解后,解线性方程组(2-21)化为求解方程组 L y = d Ly=d Ly=d和 U x = y Ux=y Ux=y。解 L y = d Ly=d Ly=d,得:
{ y 1 = d 1 y k = d k − l k y k − 1 ( 2 − 24 ) \begin{cases}y_1=d_1\\y_k=d_k-l_k y_{k-1}\end{cases}\quad\quad\quad\quad(2-24) {y1=d1yk=dk−lkyk−1(2−24)
再解 U x = y Ux=y Ux=y,得三对角方程组(2-21)的解:
{ x n = y n / u n x k = ( y k − c k x k + 1 ) / u k ( 2 − 25 ) \begin{cases}x_n=y_n/u_n\\x_k=(y_k-c_kx_{k+1})/u_k\end{cases}\quad\quad\quad\quad(2-25) {xn=yn/unxk=(yk−ckxk+1)/uk(2−25)
称上述求解三对角方程组的方法为追赶法。
算法
3.2 平方根法(Cholesky分解法)
设 A A A是对称正定矩阵,则存在唯一的非奇异下三角形矩阵 L L L,使得:
A = L L T ( 2 − 26 ) A=LL^T\quad\quad\quad\quad\quad\quad(2-26) A=LLT(2−26)
且 L L L的对角元素皆为正数。
四、误差分析
4.1 向量和矩阵的范数
向量的范数
向量的范数是衡量向量大小的度量概念,它可以用来刻画 R n R^n Rn中向量序列的收敛性。
一般地,定义范数如下:
定义2.1
设对任意向量 x ∈ R n x\in R^n x∈Rn,按一定的规则有一实数与之对应,记为 ∥ x ∥ \Vert x\Vert ∥x∥,若 ∥ x ∥ \Vert x\Vert ∥x∥满足:
- ∥ x ∥ ≥ 0 \Vert x\Vert\ge0 ∥x∥≥0,且 ∥ x ∥ = 0 \Vert x\Vert=0 ∥x∥=0当且仅当 x = 0 x=0 x=0
- 对任意实数 α \alpha α,都有 ∥ α x ∥ = ∣ α ∣ ∥ x ∥ \Vert\alpha x\Vert=|\alpha|\Vert x\Vert ∥αx∥=∣α∣∥x∥
- 对任意
x
,
y
∈
R
n
x,y\in R^n
x,y∈Rn,都有:
∥ x + y ∥ ≤ ∥ x ∥ + ∥ y ∥ \Vert x+y\Vert\le\Vert x\Vert+\Vert y\Vert ∥x+y∥≤∥x∥+∥y∥
则称 ∥ x ∥ \Vert x\Vert ∥x∥为向量 x x x的范数
按照定义2.1,向量 x x x的范数的具体形式可以有多种,常用的有以下三种:
- 2-范数(向量的模):
∥ x ∥ 2 = x 1 2 + x 2 2 + ⋯ + x n 2 ( 2 − 33 ) \Vert x\Vert_2=\sqrt{x_1^2+x_2^2+\dots+x_n^2}\quad\quad\quad\quad(2-33) ∥x∥2=x12+x22+⋯+xn2(2−33) - 1-范数:
∥ x ∥ 1 = ∑ i = 1 n ∣ x i ∣ ( 2 − 34 ) \Vert x\Vert_1=\sum\limits_{i=1}^n|x_i|\quad\quad\quad(2-34) ∥x∥1=i=1∑n∣xi∣(2−34) -
∞
\infty
∞-范数:
∥ x ∥ ∞ = max 1 ≤ i ≤ n { ∣ x i ∣ } ( 2 − 35 ) \Vert x\Vert_{\infty}=\max\limits_{1\le i\le n}\{|x_i|\}\quad\quad\quad(2-35) ∥x∥∞=1≤i≤nmax{∣xi∣}(2−35)
如果 R n R^n Rn中有两个范数 ∥ ⋅ ∥ \Vert \cdot\Vert ∥⋅∥和 ∥ ⋅ ∥ ′ \Vert\cdot\Vert' ∥⋅∥′,存在实数 n , M > 0 n,M>0 n,M>0,使得对任意 n n n维向量 x x x,都有:
m ∥ x ∥ ≤ ∥ x ∥ ′ ≤ M ∥ x ∥ m\Vert x\Vert\le\Vert x\Vert'\le M\Vert x\Vert m∥x∥≤∥x∥′≤M∥x∥
待补充 42
矩阵的范数
向量是特殊的矩阵,矩阵是推广了的向量。矩阵满足线性空间的8条性质(?),所以我们可以说矩阵是线性空间。同样的我们也可以验证向量也满足线性空间的要求(?),这是矩阵和向量的共性。在Kronecker积那节中(?)介绍了Vector转化的概念。 m × n m\times n m×n维矩阵可以转化成 m × n m\times n m×n维空间的向量。因此在学习过程中可以将矩阵和向量结合起来学习。
不要忘记,引入范数的目的是为了进行度量。就如同内积的概念一样。所有的向量空间都可以定义内积空间(?),引入内积不是目的,引入内积之后,就可以引入夹角长度等概念。这个空间就变得可以度量了。
额外说一点,并不是只有线性空间才有范数的定义,任意空间都可以引进范数(这样的空间我们称为赋范空间),使得这个空间可以被度量。如希尔伯特空间等。
定义2.2
对任意 n n n阶方阵 A A A,按一定的规则有一实数与之对应,记为 ∥ A ∥ \Vert A\Vert ∥A∥,若 ∥ A ∥ \Vert A\Vert ∥A∥满足:
- 正定性: ∥ A ∥ ≥ 0 \Vert A\Vert\ge0 ∥A∥≥0,且 ∥ A ∥ = 0 \Vert A\Vert=0 ∥A∥=0当且仅当 A = 0 A=0 A=0
- 齐次性:对任意实数 α \alpha α,都有 ∥ α A ∥ = ∣ α ∣ ∥ A ∥ \Vert\alpha A\Vert=|\alpha|\Vert A\Vert ∥αA∥=∣α∣∥A∥
- 三角不等式:对任意两个
n
n
n阶方阵
A
,
B
A,B
A,B,都有
∥ A + B ∥ ≤ ∥ A ∥ + ∥ B ∥ \Vert A+B\Vert\le\Vert A\Vert+\Vert B\Vert ∥A+B∥≤∥A∥+∥B∥ - ∥ A B ∥ ≤ ∥ A ∥ ⋅ ∥ B ∥ \Vert AB\Vert\le\Vert A\Vert\cdot\Vert B\Vert ∥AB∥≤∥A∥⋅∥B∥(相容性条件)
则称 ∥ A ∥ \Vert A\Vert ∥A∥为矩阵 A A A的范数
这里 ( 1 ) ∼ ( 3 ) (1)\sim(3) (1)∼(3)是向量范数定义的直接推广,条件(4)则使矩阵范数在数值计算中使用更为方便
定理2.5
设 A A A为 n n n阶方阵, ∥ ⋅ ∥ \Vert\cdot\Vert ∥⋅∥是 R n R^n Rn中的向量范数,则:
∥ A ∥ = max x ≠ 0 ∥ A x ∥ ∥ x ∥ \Vert A\Vert=\max\limits_{x\ne0}\frac{\Vert Ax\Vert}{\Vert x\Vert} ∥A∥=x=0max∥x∥∥Ax∥
是一种矩阵范数,称其为由向量范数 ∥ ⋅ ∥ \Vert\cdot\Vert ∥⋅∥诱导出的矩阵范数,也称算子范数。
由定理2.5的证明知,由向量范数诱导出的矩阵范数除满足定义2.2中的 ( 1 ) ∼ ( 4 ) (1)\sim(4) (1)∼(4)外,还满足性质
- 对任意 n n n维向量 x x x,都有 ∥ A x ∥ ≤ ∥ A ∥ ∥ x ∥ \Vert Ax\Vert\le\Vert A\Vert\Vert x\Vert ∥Ax∥≤∥A∥∥x∥
这一性质称为矩阵范数与向量范数的相容性,在误差分析中经常会用到。
常用的矩阵范数有四种,前三种是由向量的1-范数、2-范数、 ∞ \infty ∞-范数诱导出的矩阵范数。设 A = ( a i j ) A=(a_{ij}) A=(aij)为 n n n阶方阵:
待补充 44
4.2 方程组的状态与条件数
一个实际问题转化为数学问题,初始数据往往会有误差,即有扰动,从而使计算结果产生误差,因此需要研究扰动对解的影响。
待补充 45
当一个方程组,由于系数矩阵或右端项的微小扰动,而引起解发生巨大变化时,称该方程组是“病态”的。为了定量刻画方程组“病态”的程度,下面对方程组 A x = b Ax=b Ax=b就系数矩阵或右端项分别有扰动的两种情形进行讨论。
首先考虑右端项 b b b的扰动对解的影响。设 b b b有扰动 δ b \delta b δb,相应的解 x x x的扰动记为 δ x \delta x δx,即:
A ( x + δ x ) = b + δ b A(x+\delta x)=b+\delta b A(x+δx)=b+δb
由 A x = b Ax=b Ax=b,得 A δ x = δ b A\delta x=\delta b Aδx=δb,从而有:
δ x = A − 1 δ b \delta x=A^{-1}\delta b δx=A−1δb
两边取范数,得:
∥ δ x ∥ = ∥ A − 1 δ b ∥ ≤ ∥ A − 1 ∥ ∥ δ b ∥ \Vert\delta x\Vert=\Vert A^{-1}\delta b\Vert\le\Vert A^{-1}\Vert\Vert\delta b\Vert ∥δx∥=∥A−1δb∥≤∥A−1∥∥δb∥
又因为
∥ x ∥ > \Vert x\Vert> ∥x∥>
五、超定线性方程组的最小二乘解
如果方程(2-1)中方程个数多于变量的个数,称之为超定线性方程组。这种方程组一般无精确解,只能在各种规定的意义下求它的广义解。
设有超定线性方程组:
{ 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 m 1 x 1 + a m 2 x 2 + . . . + a m n x n = b m ( 2 − 43 ) \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\\ \dots \dots \dots\dots\dots\dots\\ a_{m1}x_1+a_{m2}x_2+...+a_{mn}x_n=b_m\\ \end {cases}\quad\quad\quad\quad\quad(2-43) ⎩⎪⎪⎪⎨⎪⎪⎪⎧a11x1+a12x2+...+a1nxn=b1a21x1+a22x2+...+a2nxn=b2………………am1x1+am2x2+...+amnxn=bm(2−43)
其矩阵形式仍记为 A x = b Ax=b Ax=b,其中 A A A 为 m × n m\times n m×n 阶矩阵, x ∈ R n , b ∈ R m ( m > n ) x\in R^n,b \in R^m\ (m>n) x∈Rn,b∈Rm (m>n)。
考虑方程组(2-43)在下述意义下求解:求 x x x,使得由它产生的每个方程的偏差 ( b i − ∑ j = 1 n a i j x j ) (b_i-\sum\limits_{j=1}^na_{ij}x_j) (bi−j=1∑naijxj) 的平方和最小,即求问题:
min F ( x ) = ∑ i = 1 m ( b i − ∑ j = 1 n a i j x j ) 2 = ∥ b − A x ∥ 2 2 ( 2 − 44 ) \min F(x)=\sum\limits_{i=1}^m(b_i-\sum\limits_{j=1}^na_{ij}x_j)^2=\Vert b-Ax\Vert_2^2\quad\quad\quad\quad(2-44) minF(x)=i=1∑m(bi−j=1∑naijxj)2=∥b−Ax∥22(2−44)
的解。问题(2-44)称为线性最小二乘问题,其解称为问题(2-43)的最小二乘解。
由多元微分学,问题(2-44)的解必定是多元函数 F ( x ) F(x) F(x)的驻点(?)。因为:
∂ F ∂ x k = 2 ∑ i = 1 m ( b i − ∑ j = 1 n a i j x j ) ( − a i k ) = − 2 ∑ i = 1 m a i k ( b i − ∑ j = 1 n a i j x j ) ( k = 1 , 2 , … , n ) \frac{\partial F}{\partial x_k}=2\sum\limits_{i=1}^m(b_i-\sum\limits_{j=1}^na_{ij}x_j)(-a_{ik})\\=-2\sum\limits_{i=1}^ma_{ik}(b_i-\sum\limits_{j=1}^na_{ij}x_j)\ (k=1,2,\dots,n) ∂xk∂F=2i=1∑m(bi−j=1∑naijxj)(−aik)=−2i=1∑maik(bi−j=1∑naijxj) (k=1,2,…,n)
因此,求解问题(2-44)可转化为求解方程组: