学习笔记十六——线性方程组的数值解法

本章主要讲的是求解方程组
AX=b() AX=b\qquad\qquad\qquad\qquad (*)
其中 ARn×nA\in R^{n\times n} 为非奇异矩阵

Gauss消元法

前提条件

消元过程的所有主元素 akk(k)0a_{kk}^{(k)}\neq0\Leftarrow \Rightarrow 系数矩阵 AAkk 阶顺序主子阵 det(Ak)(k=1,2, ,m)det(A_k)(k=1,2,\cdots,m) 均非奇异

列选主元

我们从子块(如果是构造上三角矩阵,它的左边全是零)
(ak+1,k+1(k+1)ak+2,k+1(k+1)an,k+1(k+1)) \left(\begin{array}{ccccc} a^{(k+1)}_{k+1,k+1}\\ a^{(k+1)}_{k+2,k+1}\\ \vdots\\ a^{(k+1)}_{n,k+1} \end{array}\right)
中找到绝对值最大的元素 ap,k+1(k+1)a^{(k+1)}_{p,k+1} ,将整个矩阵的第 k+1k+1 行与第 pp 行互换,从而使每次做消元时,主元素最大。

前推过程

构造形式如下:
A(n)=(a11(1)a12(1)a1n(1)a22(2)a2n(2)ann(n)),b(n)=(b1(1)b1(2)b1(n)) A^{(n)}=\left(\begin{array}{ccccc} a^{(1)}_{11}&a^{(1)}_{12}&\cdots&a^{(1)}_{1n}\\ &a^{(2)}_{22}&\cdots&a^{(2)}_{2n}\\ &&&\vdots\\ &&&a^{(n)}_{nn}\\ \end{array}\right) ,\quad b^{(n)}= \left(\begin{array}{ccccc} b_1^{(1)}\\ b_1^{(2)}\\ \vdots\\ b_1^{(n)} \end{array}\right)

回代过程

我们从第 nn 个方程开始,自下而上依次解出 xn,xn1, ,x1x_n,x_{n-1},\cdots,x_{1}

Doolittle分解法

我们记
A=LU A=LU
定理: 若矩阵 ARn×nA\in R^{n\times n} 的顺序主子式 det(Ai)0(i=1,2, ,n),det(A_i)\neq0(i=1,2,\cdots,n), 则存在唯一的下三角矩阵 LL 及上三角矩阵 UU 使得上式成立。
求解过程可以分为下列子过程:
LY=bY=(y1,y2, ,yn)TUX=YX. LY=b\Rightarrow Y=(y_1,y_2,\cdots,y_n)^T\Rightarrow UX=Y\Rightarrow X.
步骤:

  1. LL 的第一列与 AA 的第一列相同;
  2. UU 的第一行;
  3. LL 的第二列;
  4. UU 的第二行;
  5. \cdots\cdots

最后可得到 LLUU ,在得到解 XX

改进的Cholesky分解法

没看懂,建议直接看《计算方法(第二版)》的P60 。

追赶法

也就是Gauss消元法的特殊应用,没什么难,《计算方法(第二版)》的P62。

扰动分析

条件数 Cond(A):A1ACond(A):||A^{-1}||||A||。当 Cond(A)>>1Cond(A)>>1 时,方程组 ()(*) 视为病态的。常用的条件数有:
Cond1(A)=A11A1,Cond(A)=A1A. Cond_1(A)=||A^{-1}||_1||A||_1,\\ Cond_\infty(A)=||A^{-1}||_\infty||A||_\infty.
上述方式就是一般的直接法,而迭代法比直接法更适合于现代大规模科学工程计算。

一般单步迭代法

设线性方程 ()(*) 有如下迭代格式:
X(k+1)=BK(k)+F,k=0,1,2, ,() X^{(k+1)}=BK^{(k)}+F,\quad k=0,1,2,\cdots,\qquad(**)
定理(重要): 当给定初始向量 X(0)X^{(0)} 时,迭代格式 ()(**) 收敛的充要条件是其迭代矩阵 BB 的谱半径 ρ(B)<1\rho(B)<1

Jacobi迭代法

将线性方程组 ()(*) 的系数矩阵 AA 分解为
A=L+D+U, A=L+D+U,
其中 D=diag(a11,a22, ,ann),D=diag(a_{11},a_{22},\cdots,a_{nn}),
L=(0000a21000a31a3200an1an2an,n10), L=\left(\begin{array}{ccccc} 0&0&\cdots&0&0\\ a_{21}&0&\cdots&0&0\\ a_{31}&a_{32}&\cdots&0&0\\ \vdots&\vdots&&\vdots&\vdots\\ a_{n1}&a_{n2}&\cdots&a_{n,n-1}&0\\ \end{array}\right) ,\\
U=(0a12a13a1n00a23a2n000an1,n0000). U=\left(\begin{array}{ccccc} 0&a_{12}&a_{13}&\cdots&a_{1n}\\ 0&0&a_{23}&\cdots&a_{2n}\\ \vdots&\vdots&\vdots&&\vdots\\ 0&0&0&\cdots&a_{n-1,n}\\ 0&0&0&\cdots&0\\ \end{array}\right) .\\
于是有
(L+D+U)X=bDX=(L+U)X+bX=D1(L+D)X+D1b (L+D+U)X=b\\ \Rightarrow DX=-(L+U)X+b\\ \Rightarrow X=-D^{-1}(L+D)X+D^{-1}b
Jacobi迭代公式:
X(k+1)=D1(L+D)X(k)+D1b,k=0,1, , X^{(k+1)}=-D^{-1}(L+D)X^{(k)}+D^{-1}b,\quad k=0,1,\cdots,
定理(重要): 若线性方程组 ()(*) 的系数矩阵 AA 严格对角占优,则Jacobi迭代法是收敛的。

Gauss-Seidel迭代法

方程组 ()(*) 也可以等价地写为
(D+L)X=UX+b (D+L)X=-UX+b
类似Jacobi迭代法可以得到Gauss-Seidel迭代法:
X(k+1)=(D+L)1UX(k)+(D+L)1b X^{(k+1)}=-(D+L)^{-1}UX^{(k)}+(D+L)^{-1}b
定理(重要): 若线性方程组 ()(*) 的系数矩阵 AA 严格对角占优,则Gauss-Seidel迭代法是收敛的。

JOR迭代法

JOR迭代法是由Jacobi迭代法加入松弛因子 ww 得到。
由:
X(k+1)=X(k)+w X^{(k+1)}=X^{(k)}+w步长
可以得到JOR迭代法:
X(k+1)=X(k)wD1(AX(k)b). X^{(k+1)}=X^{(k)}-wD^{-1}(AX^{(k)}-b).
JOR迭代法有最佳松弛因子
wopt=22λmaxBJλminBJ, w_{opt}=\frac{2}{2-\lambda^{B_J}_{max}-\lambda^{B_J}_{min}},
其中 λmaxBJ,λminBJ\lambda^{B_J}_{max},\lambda^{B_J}_{min} 分别表示Jacobi迭代矩阵 BJ=D1(L+U)B_J=-D^{-1}(L+U) 的最大和最小特征值。此外,当 λmaxBJλminBJ\lambda^{B_J}_{max}\neq\lambda^{B_J}_{min} 时,JOR迭代法的收敛速度相较于对应的Jacobi迭代法的收敛速度快。
定理(重要): 若线性方程组 ()(*) 的系数矩阵 AA 严格对角占优,则松弛因子 w(0,1]w\in (0,1] 的JOR迭代法是收敛的。

SOR迭代法

SOR迭代法是由Gauss-Seidel迭代法加入松弛因子 ww 得到。
由:
DX(k+1)=DX(k)+w DX^{(k+1)}=DX^{(k)}+w步长
得到SOR迭代法:
X(k+1)=(D+wL)1{[(1w)DwU]X(k)+wb}. X^{(k+1)}=(D+wL)^{-1}\{[(1-w)D-wU]X^{(k)}+wb\}.
SOR迭代法的最佳松弛因子
wopt=21+1ρ2(BJ) w_{opt}=\frac{2}{1+\sqrt{1-\rho^2(B_J)}}
定理(重要): 若线性方程组 ()(*) 的系数矩阵 AA 严格对角占优,则松弛因子 w(0,1]w\in (0,1] 的SOR迭代法是收敛的。

下面是自己推导Jacobi,Gauss-Seidel,JOR,SOR的过程:
学习笔记十六——线性方程组的数值解法