§2 方程组

C1 高斯消去法

1）高斯消去法：利用行变换将系数矩阵化为下三角矩阵之后，再化为对角矩阵

• 共需要进行$\frac{n(n-1)}{2}+n$次行变换，$O(n^3)$，回代过程$O(n^2)$
• $a_{ii}$称为主元，若算法过程中主元为0，算法终止

2）LU分解法：写出高斯消去法的行变换逆矩阵$L$，得$A=LU$，先求解$Ly=b$得到$y$。再求解$Ux=y$得到$x$

• 对于多个$b$且不是同时已知的计算，节省了重复计算$L$的时间
• MATLAB：x = A \ b

3）部分主元法：选择$\underset{k\le i\le n}{\max }|a_{ik}|$为第$k$个主元

• 目的是使乘子小于1，避免淹没和0主元问题(若仍有0主元则A奇异)
• 交换过程对应的逆矩阵为$P$，在运算中，将乘子保存在对应的0上，随行交换得到$L$
• $PA=LU;LUx=Pb$，先求解$Lc=Pb$，再求解$Ux=c$

C2 误差分析

1）余项：$r=b-A\hat{x}$；后向误差：$\Vert r{\Vert }_{\infty }$；前向误差：$\Vert x-\hat{x}{\Vert }_{\infty }$；

2）相对后向误差：$\frac{\Vert r{\Vert }_{\infty }}{\Vert b{\Vert }_{\infty }}$；相对前向误差：$\frac{\Vert x-\hat{x}{\Vert }_{\infty }}{\Vert x{\Vert }_{\infty }}$

3）误差放大因子：相对前向误差÷相对后向误差

4）条件数：求解$Ax=b$时可能出现的最大误差放大因子$\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}(A)=\Vert A{\Vert }_{\infty }\cdot \Vert A^{-1}{\Vert }_{\infty }$

证明：利用无穷范数是算子范数$\frac{\Vert x-\hat{x}{\Vert }_{\infty }}{\Vert x{\Vert }_{\infty }}\le \Vert A^{-1}{\Vert }_{\infty }\Vert r{\Vert }_{\infty },\frac{1}{\Vert b\Vert }\ge \frac{1}{\Vert A{\Vert }_{\infty }\Vert x{\Vert }_{\infty }}$

• 若$\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}(A)\approx 10^k$，将丢失$k$位精度

5）淹没(swamp)：由于乘子过大和舍入误差导致方程丢失

C3 迭代方法

1）Jacobi方法：

• 步骤：

  • 将矩阵分解为对角线，上三角，下三角三部分：$A=L+D+U$
  • $x_{k+1}=D^{-1}(b-(L+U)x_k)$

• 收敛性：

  • 严格对角占优：$\forall i,|a_{ii}|>\sum _{j\ne i}|a_{ij}|$

  • 定理：严格对角占优矩阵非奇异且使用雅可比方法收敛

    $x_{k+1}=D^{-1}(b-(L+U)x_k)=D^{-1}b-D^{-1}(L+U)x_k)$

    谱半径$\rho (D^{-1}(L+U))<1$收敛

2）Gauss-Seidel方法：

• $x_{k+1}=D^{-1}(b-U_k-L{X}_{k+1})$
• 收敛性：若Jacobi方法收敛，则Gauss-Seidel方法收敛得更快

3）连续过松弛(SOR)方法：

• $x_{k+1}=(\omega L+D{)}^{-1}[(1-\omega )D{x}_k-\omega U{x}_k]+\omega (D+\omega L^{-1})b$
• $\omega$称为松弛参数，$\omega >1$称过松弛，$\omega <1$称欠松弛。$\omega$的选择影响最终误差

4）应用：

• 修饰：当$b$发生微小变换时，利用原解作为初始估计进行迭代
• 稀疏矩阵在使用LU分解时常发生填充，使用迭代法则可避免

C4 对称正定矩阵解法

定理：对称矩阵正定$⟺$只有正特征值

定理：$A$对称正定且$X$满秩，则$X^{T}AX$对称正定

定理：若$A$对称正定，则$A$的任意主子矩阵对称正定

1）Cholesky分解：

• 任意对称正定矩阵$A$，存在上三角矩阵$R$使得$A=R^{T}R$
• 步骤：$R_A=\left[\begin{array}{cc}\sqrt{a_{11}}& u\\ 0& R_{A_{2:n}-u^{T}u}\end{array}\right],u=b/\sqrt{(}a),A_{2:n}$指$A$去除第一行第一列的主子矩阵
• $R^{T}Rx=b$求解，即先求解$R^{T}c=b$，再求解$Rx=c$

2）共轭梯度算法：

• $x_0$为初始估计，$d_0=b-A{x}_0$为初始搜索方向,$r$表示余项,$\alpha$为步长

• ${\alpha }_k=\frac{r_k^T{r}_k}{d_k^{T}A{d}_k}$

  $x_{k+a}=x_k+{\alpha }_k{d}_k$

  $r_{k+1}=r_k-{\alpha }_{k}A{d}_k$

  ${\beta }_k=\frac{r_{k+1}^T{r}_{k+1}}{r_k^T{r}_k}$

  $d_{k+1}=r_{k+1}+{\beta }_k{d}_k$

• 在n步之内一定能找到解使得$r_k=0$

• 对于非稀疏矩阵，可能意味着更大的计算代价

• 在得到足够小的$r_k$后，也可以停止计算

3）预条件子：通过$M^{-1}Ax=M^{-1}b$减少条件数，$M^{-1}$称为预条件子

• Jacobi预条件子：$D$

• $M=A$，则条件数为1；$M=E$，并不减少条件数

• 预条件共轭梯度算法

  • ${\alpha }_k=\frac{r_k^T{z}_k}{d_k^{T}A{d}_k}$

    $x_{k+1}=x_k+{\alpha }_k{d}_k$

    $z_{k+1}=M^{-1}{r}_{k+1},z_0=d_0=M^{-1}{r}_0$

    ${\beta }_k=\frac{r_{k+1}^T{z}_{k+1}}{r_k^T{z}_k}$

    $d_{k+1}=z_{k+1}+{\beta }_k{d}_k$

• Gauss-Seidel预条件子：$M=(D+\omega L)D^{-1}(D+\omega U)$

C5 非线性方程组

1）多元牛顿方法：$x_{k+1}=x_k+s,DF(x_k)s=-F(x_k)$

2）Broyden方法：难以求出雅可比矩阵时的次优方法

• 法Ⅰ：使用近似雅可比矩阵$A_i,A_{i+1}=A_i+\frac{({\Delta }_{i+1}-A_i{\delta }_i){\delta }_{i+1}^T}{{\delta }_{i+1}^T{\delta }_{i+1}}$

• 法Ⅱ：使用近似雅可比矩阵逆$B_i,B_{i+1}=B_i+\frac{({\delta }_{i+1}-B_i{\Delta }_{i+1})\delta {i+1}^T{B}_i}{\delta {i+1}^T{B}_i{\Delta }_{i+1}},{\delta }_i=x_i-x_{i-1},{\Delta }_i=F(x_i)-F(x_{i-1})$

  $x_{k+1}=x_k-B_{k}F(x_k)$。$B_0$使用初始估计的近似雅可比矩阵，最差是单位阵