科学计算(1):线性方程组解法

打印 上一主题 下一主题

主题 981|帖子 981|积分 2943

系列笔记为上科学计算课上随意记下的笔记,大概有一些缭乱,勉强拼集着看。以记载和梳理为主,证明较少。
一、直接法

1. 高斯消去法

计算复杂度分析


  • 乘除法次数
                                                              ∑                                               k                                     =                                     1                                              n                                                      [                                  (                                  n                                  −                                  k                                  +                                  2                                  )                                  (                                  n                                  −                                  k                                  )                                  +                                  (                                  n                                  −                                  k                                  +                                  1                                  )                                  ]                                          =                                                        n                                     3                                              3                                          +                                           n                                  2                                          −                                           n                                  3                                                 \sum_{k=1}^n \left[(n-k+2)(n-k) + (n-k+1)\right] = \frac{n^3}{3} + n^2 - \frac{n}{3}                        k=1∑n​[(n−k+2)(n−k)+(n−k+1)]=3n3​+n2−3n​
  • 加减法次数
                                                              ∑                                               k                                     =                                     1                                              n                                                      [                                  (                                  n                                  −                                  k                                  +                                  1                                  )                                  (                                  n                                  −                                  k                                  )                                  ]                                          +                               (                               n                               −                               k                               +                               1                               )                               =                                                        n                                     3                                              3                                          +                                           1                                  2                                                      n                                  2                                          +                                           1                                  6                                          n                                      \sum_{k=1}^n \left[(n-k+1)(n-k)\right] +(n-k+1)= \frac{n^3}{3} +\frac{1}{2}n^{2}+\frac{1}{6}n                        k=1∑n​[(n−k+1)(n−k)]+(n−k+1)=3n3​+21​n2+61​n
  • 总复杂度:                                        O                            (                                       n                               3                                      )                                  O(n^3)                     O(n3),以                                                                n                                  3                                          3                                            \frac{n^3}{3}                     3n3​ 为主项
2. 选主元消去法



  • 列主元法:在列中选择最大元素作为主元,避免小主元导致的数值不稳定
  • 全主元法:在整个子矩阵中选择最大主元(计算量更大,但稳定性最佳)
3. LU 分解

定理:若矩阵                                    A                         ∈                                   R                                       n                               ×                               n                                                 A \in \mathbb{R}^{n \times n}                  A∈Rn×n 的次序主子式均非零,则存在唯一的单元下三角矩阵                                    L                              L                  L 和上三角矩阵                                    U                              U                  U,使得                                    A                         =                         L                         U                              A = LU                  A=LU。
存储优化


  •                                         U                                  U                     U 存储于原矩阵的上三角部分
  •                                         L                                  L                     L 的非单元元素存储于下三角部分(对角元1无需存储)
Crout 分解:将                                    A                              A                  A 分解为下三角矩阵与单元上三角矩阵的乘积。
明白LU分解并推得递推表达式,可以从初等矩阵乘积进行消元的角度进行明白,也可以直接从LU的形式入手,直接进行计算。从初等矩阵乘积消元的角度看,利用列主元LU分解,也就是选取列主元消元,我们可以得到以下结论:
列主元 LU 分解
定理:对恣意非奇异矩阵                                    A                              A                  A,存在置换矩阵                                    P                              P                  P,使得:
                                         P                            A                            =                            L                            U                                  PA = LU                     PA=LU
此中                                    L                              L                  L 为单元下三角矩阵,                                   U                              U                  U 为上三角矩阵。
实现步骤
假设前r-1行都进行了换行消元的过程,即                                             l                                       i                               k                                                 l_{ik}                  lik​                                    (                         i                         =                         1                         ,                         ⋯                          ,                         n                         ;                         k                         =                              (i=1,\cdots,n;k=                  (i=1,⋯,n;k=                                    1                         ,                         ⋯                         r                         −                         1                         )                              1,\cdots r-1)                  1,⋯r−1) )及                                             u                                       k                               i                                                 u_{ki}                  uki​                                    (                         i                         =                         1                         ,                         ⋯                          ,                         n                         ;                         k                         =                         1                         ,                         ⋯                         r                         −                         1                         )                              (\boldsymbol{i}=1,\cdots,n;k=1,\cdots r-1)                  (i=1,⋯,n;k=1,⋯r−1) 均已得到。下面要得到必要互换的两行以及U的第r行。换行的尺度是希望形成的                                             u                                       r                               r                                                 u_{rr}                  urr​ 绝对值最大。考虑A的第r列的r至n元素                                             a                                       i                               r                                                 a_{ir}                  air​ ,i=r…n,有
                                                    a                                           i                                  r                                                 =                                       ∑                                           k                                  =                                  1                                                      r                                  −                                  1                                                            l                                           i                                  k                                                            u                                           k                                  r                                                 +                                                   u                                               r                                     r                                                      ^                                            a_{ir}=\sum_{k=1}^{r-1}l_{ik}u_{kr}+\widehat{u_{rr}}                     air​=k=1∑r−1​lik​ukr​+urr​             ​
此中                                             s                            i                                  =                              s_{i}=                  si​=                                                         u                                           r                                  r                                                 ^                                       \widehat{u_{rr}}                  urr​            ​ 为各行若称为主行后,大概得到的                                             u                                       r                               r                                                 u_{rr}                  urr​ .比较                                             ∣                                       s                                           i                                  r                                                 ∣                                  =                         max                         ⁡                         ∣                                   s                            i                                  ∣                              \left|s_{i_r}\right|=\max|s_i|                  ∣sir​​∣=max∣si​∣ ,互换第r行与                                             i                            r                                       i_r                  ir​ 行,
4. 对称矩阵的特殊分解

1. LDLᵀ 分解

定理:若对称矩阵                                    A                              A                  A 的次序主子式非零,则存在唯一分解:
                                         A                            =                            L                            D                                       L                               T                                            A = LDL^T                     A=LDLT
此中                                    L                              L                  L 为单元下三角阵,                                   D                              D                  D 为对角阵。
2. Cholesky 分解

定理:若                                    A                              A                  A 对称正定,则存在唯一的下三角矩阵                                    L                              L                  L(对角元全正),使得:
                                         A                            =                            L                                       L                               T                                            A = LL^T                     A=LLT
计算步骤
                                                                                           l                                                       k                                           k                                                                                                                                     =                                                                                        a                                                                   k                                                    k                                                                               −                                                               ∑                                                                   m                                                    =                                                    1                                                                                    k                                                    −                                                    1                                                                                                l                                                                   k                                                    m                                                                  2                                                                                                                                                                         l                                                       i                                           k                                                                                                                                     =                                                       1                                                           l                                                               k                                                 k                                                                                                    (                                                           a                                                               i                                                 k                                                                          −                                                           ∑                                                               m                                                 =                                                 1                                                                               k                                                 −                                                 1                                                                                          l                                                               i                                                 m                                                                                          l                                                               k                                                 m                                                                          )                                                                   (                                        i                                        >                                        k                                        )                                                                                \begin{aligned} l_{kk} &= \sqrt{a_{kk} - \sum_{m=1}^{k-1} l_{km}^2} \\ l_{ik} &= \frac{1}{l_{kk}} \left( a_{ik} - \sum_{m=1}^{k-1} l_{im}l_{km} \right) \quad (i > k) \end{aligned}                     lkk​lik​​=akk​−m=1∑k−1​lkm2​                     ​=lkk​1​(aik​−m=1∑k−1​lim​lkm​)(i>k)​

二、迭代法

1. 根本理论

定理1:下列命题等价

  •                                                                lim                                  ⁡                                                      k                                  →                                  ∞                                                            B                               k                                      =                            0                                  \lim_{k\to\infty}B^k=0                     limk→∞​Bk=0
  •                                         ρ                            (                            B                            )                            <                            1                                  \rho(B)<1                     ρ(B)<1
  • 存在                                        ∥                            ⋅                                       ∥                               ε                                            \|{\cdot}\|_{\varepsilon}                     ∥⋅∥ε​,使                                                               ∥                                  B                                  ∥                                          ε                                      <                            1                                  \left\|\boldsymbol{B}\right\|_{\varepsilon}<\mathbf{1}                     ∥B∥ε​<1
定理2:对恣意一种范数,有
                                                                lim                                  ⁡                                                      k                                  →                                  ∞                                                                        ∣                                               ∣                                                   B                                        k                                                  ∣                                              ∣                                                      1                                  /                                  k                                                 =                            ρ                            (                            B                            )                                  \lim_{k\to\infty}\left|\left|B^{k}\right|\right|^{1/k}=\rho(B)                     k→∞lim​                ​                 ​Bk                 ​                ​1/k=ρ(B)
这个定理告诉我们,无论我们选取任何一种范数,其收敛举动都是由谱半径决定的。
下面先容迭代法的根本定理。
迭代格式:                                             x                                       (                               k                               +                               1                               )                                            =                         B                                   x                                       (                               k                               )                                            +                         g                              x^{(k+1)} = Bx^{(k)} + g                  x(k+1)=Bx(k)+g
收敛定理:迭代法收敛当且仅当                                    ρ                         (                         B                         )                         <                         1                              \rho(B) < 1                  ρ(B)<1(谱半径小于1)
毛病估计
                                                                                           ∥                                                       x                                           ∗                                                      −                                                       x                                                           (                                              k                                              )                                                                     ∥                                                                                                                       ≤                                        ∥                                        B                                                       ∥                                           k                                                      ⋅                                        ∥                                                       x                                                           (                                              ∗                                              )                                                                     −                                                       x                                                           (                                              0                                              )                                                                     ∥                                                                                                                            ∥                                                       x                                           ∗                                                      −                                                       x                                                           (                                              k                                              )                                                                     ∥                                                                                                                       ≤                                                                       ∥                                              B                                                               ∥                                                 k                                                                                          1                                              −                                              ∥                                              B                                              ∥                                                                     ∥                                                       x                                                           (                                              1                                              )                                                                     −                                                       x                                                           (                                              0                                              )                                                                     ∥                                                                                                                            ∥                                                       x                                           ∗                                                      −                                                       x                                                           (                                              k                                              )                                                                     ∥                                                                                                                       ≤                                                                       ∥                                              B                                              ∥                                                                          1                                              −                                              ∥                                              B                                              ∥                                                                     ∥                                                       x                                                           (                                              k                                              )                                                                     −                                                       x                                                           (                                              k                                              −                                              1                                              )                                                                     ∥                                                                                \begin{aligned} \|x^* - x^{(k)}\| &\leq \|B\|^k\cdot \|x^{(*)} - x^{(0)}\| \\ \|x^* - x^{(k)}\| &\leq \frac{\|B\|^k}{1 - \|B\|} \|x^{(1)} - x^{(0)}\| \\ \|x^* - x^{(k)}\| &\leq \frac{\|B\|}{1 - \|B\|} \|x^{(k)} - x^{(k-1)}\| \end{aligned}                     ∥x∗−x(k)∥∥x∗−x(k)∥∥x∗−x(k)∥​≤∥B∥k⋅∥x(∗)−x(0)∥≤1−∥B∥∥B∥k​∥x(1)−x(0)∥≤1−∥B∥∥B∥​∥x(k)−x(k−1)∥​
收敛速度分析
对恣意给定的范数,对恣意                                   δ                         >                         0                              \delta>0                  δ>0 ,存在着正整数K使得                                   k                         ≥                         K                              k\geq K                  k≥K 时,
                                                    ∣                                           ∣                                               ε                                                   (                                        k                                        )                                                           ∣                                          ∣                                      =                                       ∣                                           ∣                                               B                                     k                                                           ε                                                   (                                        0                                        )                                                           ∣                                          ∣                                      ≤                                       ∣                                           ∣                                               B                                     k                                              ∣                                          ∣                                                 ∣                                           ∣                                               ε                                                   (                                        0                                        )                                                           ∣                                          ∣                                      ≤                                                   (                                  ρ                                  (                                  B                                  )                                  +                                  δ                                  )                                          k                                                 ∣                                           ∣                                               ε                                                   (                                        0                                        )                                                           ∣                                          ∣                                            \left|\left|\varepsilon^{(k)}\right|\right|=\left|\left|B^{k}\varepsilon^{(0)}\right|\right|\leq\left|\left|B^{k}\right|\right|\left|\left|\varepsilon^{(0)}\right|\right|\leq\left(\rho(B)+\delta\right)^{k}\left|\left|\varepsilon^{(0)}\right|\right|                                    ​                ​ε(k)                ​               ​=               ​                ​Bkε(0)                ​               ​≤               ​                ​Bk                ​               ​               ​                ​ε(0)                ​               ​≤(ρ(B)+δ)k               ​                ​ε(0)                ​               ​
于是有:
                                                                ∣                                  ∣                                               ε                                                   (                                        k                                        )                                                           ∣                                  ∣                                                      ∣                                  ∣                                               ε                                                   (                                        0                                        )                                                           ∣                                  ∣                                                 ≤                            (                            ρ                            (                            B                            )                            +                            δ                                       )                               k                                      ≤                            τ                                  \frac{||\varepsilon^{(k)}||}{||\varepsilon^{(0)}||}\leq(\rho(B)+\delta)^k\leq\tau                     ∣∣ε(0)∣∣∣∣ε(k)∣∣​≤(ρ(B)+δ)k≤τ
于是为达到精度                                    τ                              \tau                  τ,需迭代步数:
                                         k                            ≥                                                   −                                  ln                                  ⁡                                  τ                                                      −                                  ln                                  ⁡                                  (                                  ρ                                  (                                  B                                  )                                  +                                  δ                                  )                                                       k \geq \frac{-\ln \tau}{-\ln (\rho(B) + \delta)}                     k≥−ln(ρ(B)+δ)−lnτ​
2. 矩阵分裂法

将                                    A                              A                  A 分裂为                                    A                         =                         M                         −                         N                              A = M - N                  A=M−N,迭代公式:
                                                    x                                           (                                  k                                  +                                  1                                  )                                                 =                                       M                                           −                                  1                                                 N                                       x                                           (                                  k                                  )                                                 +                                       M                                           −                                  1                                                 b                                  x^{(k+1)} = M^{-1}Nx^{(k)} + M^{-1}b                     x(k+1)=M−1Nx(k)+M−1b
Richardson迭代

将A分裂为:
                                         A                            =                            I                            −                            (                            I                            −                            A                            )                                  A = \mathcal{I} - (\mathcal{I} - A)                     A=I−(I−A)
于是:
                                                    x                                           (                                  k                                  )                                                 =                            (                            I                            −                            A                            )                                       x                                           (                                  k                                  −                                  1                                  )                                                 +                            b                                  x^{(k)} = (\mathcal{I}-A)x^{(k-1)}+b                     x(k)=(I−A)x(k−1)+b
可以写成更有启发性的形式:
                                                    x                                           (                                  k                                  )                                                 =                                       x                                           (                                  k                                  −                                  1                                  )                                                 +                                       (                                           b                                  −                                  A                                               x                                                   (                                        k                                        −                                        1                                        )                                                                   )                                            x^{(k)} = x^{(k-1)} + \left({b-Ax^{(k-1)}}\right)                     x(k)=x(k−1)+(b−Ax(k−1))
以及更一般的形式:
                                                    x                                           (                                  k                                  )                                                 =                                       x                                           (                                  k                                  −                                  1                                  )                                                 +                            α                                       (                                           b                                  −                                  A                                               x                                                   (                                        k                                        −                                        1                                        )                                                                   )                                            x^{(k)} = x^{(k-1)} + \alpha\left({b-Ax^{(k-1)}}\right)                     x(k)=x(k−1)+α(b−Ax(k−1))
一般来说,在Richardson迭代中,我们要求矩阵                                   A                              A                  A是一个对称正定矩阵,此中正定矩阵是因为要求谱半径绝对值要求严酷小于1。于是迭代矩阵的特征值:
                                         λ                            [                            I                            −                            α                            A                            ]                            =                            1                            −                            α                            A                                  \lambda[\mathcal{I}-\alpha A] = 1-\alpha A                     λ[I−αA]=1−αA
收敛条件为:
                                         α                                       λ                                           max                                  ⁡                                                  (                            A                            )                            <                            2                                  \alpha \lambda_{\max\: }(A) <2                     αλmax​(A)<2
外推法:
                                                                                                                                                                             x                                                           (                                              k                                              )                                                                      =                                         ν                                        [                                        (                                        I                                        −                                        A                                        )                                                       x                                                           (                                              k                                              −                                              1                                              )                                                                      +                                         b                                        ]                                        +                                        (                                        1                                        −                                        ν                                        )                                                       x                                                           (                                              k                                              −                                              1                                              )                                                                                                                                                                                                                            充要条件:                                        0                                        <                                        ν                                        <                                                       2                                                           λ                                                               m                                                 a                                                 x                                                                                     .                                                                                \begin{aligned}&x^{(k)}\:=\:\nu\Big[(I-A)x^{(k-1)}\:+\:b\Big]+(1-\nu)x^{(k-1)}\\&\text{充要条件:}0<\nu<\frac{2}{\lambda_{max}}.\end{aligned}                     ​x(k)=ν[(I−A)x(k−1)+b]+(1−ν)x(k−1)充要条件:0<ν<λmax​2​.​
Jacobi迭代与Gauss-Seidel迭代



  • Jacobi 迭代:                                             M                               =                               D                                      M = D                        M=D(对角阵)
                                                                    x                                     i                                                   (                                        k                                        )                                                           =                                               1                                                   a                                                       i                                           i                                                                                      (                                                   b                                        i                                                  −                                                   ∑                                                       j                                           ≠                                           i                                                                              a                                                       i                                           j                                                                              x                                        j                                                       (                                           k                                           −                                           1                                           )                                                                )                                                      x_i^{(k)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j \neq i} a_{ij}x_j^{(k-1)} \right)                           xi(k)​=aii​1​                 ​bi​−j=i∑​aij​xj(k−1)​                 ​

    • 可并行计算,需存储两代迭代值
    • 存储需求:需保存两个向量                                                                x                                                   (                                        k                                        )                                                           ,                                               x                                                   (                                        k                                        +                                        1                                        )                                                                   x^{(k)},x^{(k+1)}                           x(k),x(k+1)。
    • 收敛速度:通常较慢,但对角占优或对称正定矩阵包管收敛。

  • Gauss-Seidel 迭代:                                             M                               =                               D                               −                               L                                      M = D - L                        M=D−L(下三角部分)
                                                                    x                                     i                                                   (                                        k                                        )                                                           =                                               1                                                   a                                                       i                                           i                                                                                      (                                                   b                                        i                                                  −                                                   ∑                                                       j                                           <                                           i                                                                              a                                                       i                                           j                                                                              x                                        j                                                       (                                           k                                           )                                                                −                                                   ∑                                                       j                                           >                                           i                                                                              a                                                       i                                           j                                                                              x                                        j                                                       (                                           k                                           −                                           1                                           )                                                                )                                                      x_i^{(k)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j < i} a_{ij}x_j^{(k)} - \sum_{j > i} a_{ij}x_j^{(k-1)} \right)                           xi(k)​=aii​1​(bi​−j<i∑​aij​xj(k)​−j>i∑​aij​xj(k−1)​)

    • 串行计算,节省存储空间

收敛性分析

设                                   A                         =                                              (                                           a                                               i                                     j                                                      )                                                 n                               ×                               n                                                 A=\left(a_{ij}\right)_{n\times n}                  A=(aij​)n×n​
如果满足
                                                                                           ∣                                                       a                                                           i                                              i                                                                     ∣                                                                                                                       >                                                       ∑                                                           j                                              =                                              1                                              ,                                              j                                              ≠                                              i                                                          n                                                                     a                                                           i                                              j                                                                                  (                                        i                                        =                                        1                                        ,                                        2                                        ,                                        ⋯                                         ,                                        n                                        )                                                                                \begin{aligned}|a_{ii}|&>\sum_{j=1,j\neq i}^{n}a_{ij}\quad(i=1,2,\cdots,n)\end{aligned}                     ∣aii​∣​>j=1,j=i∑n​aij​(i=1,2,⋯,n)​
即每一行的对角线元素绝对值都严酷大于同行其他元素的绝对值之和, 则称A为严酷对角占优阵。如果满足
                                         ∣                                       a                                           i                                  i                                                 ∣                            ≥                                       ∑                                                        j                                     =                                     1                                                           j                                     ≠                                     i                                                      n                                                 a                                           i                                  j                                                          (                            i                            =                            1                            ,                            2                            ,                            ⋯                             ,                            n                            )                                  |a_{ii}|\geq\sum_{j=1\atop j\neq i}^{n}a_{ij}\quad(i=1,2,\cdots,n)                     ∣aii​∣≥j=ij=1​∑n​aij​(i=1,2,⋯,n)
且至少有一个不等式严酷建立,则称A为弱对角占优阵。如果存在                                   n                              n                  n 阶置换矩阵                                   P                              P                  P 使
                                                    P                               T                                      A                            P                            =                                       [                                                                                                     A                                              11                                                                                                                            A                                              12                                                                                                                                  0                                                                                                             A                                              22                                                                                                ]                                            P^TAP=\begin{bmatrix}A_{11}&A_{12}\\0&A_{22}\end{bmatrix}                     PTAP=[A11​0​A12​A22​​]
建立,此中                                             A                            11                                       A_{11}                  A11​ 为r阶子矩阵,                                              A                            22                                       A_{22}                  A22​ 为                                   n                         −                         r                              n-r                  n−r 阶子矩阵,则称A 是可约矩阵。否则,称A为不可约矩阵。
定理1:如果                                   A                         =                         (                                   a                                       i                               j                                            )                              A=(a_{ij} )                  A=(aij​)是严酷对角占优矩阵或不可约弱对角占优矩阵,则Jacobi和Gauss-Seidel方法均收敛
起首分析Jacobi方法。等价于证明矩阵:
                                                    B                               J                                      =                            I                            −                                       D                                           −                                  1                                                 A                                  B_{J} = \mathcal{I}-D^{-1}A                     BJ​=I−D−1A
的所有特征值的模都小于1。等价于证明对于恣意模大于1的特征值                                   λ                              \lambda                  λ:
                                         ∣                            λ                            I                            −                                       B                               J                                      ∣                            ≠                            0                                  \lvert \lambda \mathcal{I} - B_J\rvert\neq 0                     ∣λI−BJ​∣=0
而:
                                                                                           ∣                                        λ                                        I                                        −                                                       B                                           J                                                      ∣                                                                                                                       =                                        ∣                                        λ                                        I                                        −                                                       D                                                           −                                              1                                                                                    (                                                           L                                              +                                              U                                                          )                                                      ∣                                                                                                                                                                                               =                                        ∣                                                       D                                                           −                                              1                                                                     ∣                                        ⋅                                        ∣                                        λ                                        D                                        −                                                       (                                                           L                                              +                                              U                                                          )                                                      ∣                                                                                \begin{align*} \lvert \lambda \mathcal{I}-B_{J}\rvert &= \lvert \lambda \mathcal{I} - D^{-1}\left({L+U}\right)\rvert\\ & = \lvert D^{-1}\rvert\cdot\lvert \lambda D - \left({L+U}\right)\rvert \end{align*}                     ∣λI−BJ​∣​=∣λI−D−1(L+U)∣=∣D−1∣⋅∣λD−(L+U)∣​
于是这是显然的。
定理2:如果                                   A                              A                  A 对称,                                              a                                       i                               i                                            >                         0                              a_{ii}>0                  aii​>0则:
1.Gauss-Seidel方法收敛的充分须要条件是A正定;
2.Jacobi方法收敛的充分须要条件是A与2D-A均正定。
这里我们只证明Jacobi方法。
                                                    B                               J                                      =                            I                            −                                       D                                           −                                  1                                                 A                              ⟹                              λ                                       (                                           B                                  J                                          )                                      =                            λ                                       (                                           I                                  −                                               D                                                   −                                        1                                                           A                                          )                                      =                            1                            −                            λ                                       (                                                        D                                                   −                                        1                                                           A                                          )                                            B_{J} = \mathcal{I} - D^{-1}A\implies\lambda \left({B_{J}}\right) = \lambda \left({\mathcal{I} - D^{-1}A}\right) = 1 - \lambda \left({D ^{-1}A}\right)                     BJ​=I−D−1A⟹λ(BJ​)=λ(I−D−1A)=1−λ(D−1A)
而:
                                         λ                            (                                       D                                           −                                  1                                                 A                            )                            =                            λ                                       (                                                        D                                                   −                                        1                                        /                                        2                                                           A                                               D                                                   −                                        1                                        /                                        2                                                                   )                                            \lambda(D^{-1}A) = \lambda \left({D^{-1/2}AD^{-1/2}}\right)                     λ(D−1A)=λ(D−1/2AD−1/2)
于是:
                                         λ                            (                                       D                                           −                                  1                                                 A                            )                            >                            0                              ⟺                              A                            正定                                  \lambda(D^{-1}A)>0 \iff A正定                     λ(D−1A)>0⟺A正定
                                         λ                            (                                       D                                           −                                  1                                                 A                            )                            =                            λ                            (                                       D                                           −                                  1                                  /                                  2                                                 A                                       D                                           −                                  1                                  /                                  2                                                 )                            <                            2                              ⟺                              2                            D                            −                            A                            正定                                  \lambda(D^{-1}A) = \lambda(D^{-1/2}AD^{-1/2})<2 \iff 2\mathcal{D} - A正定                     λ(D−1A)=λ(D−1/2AD−1/2)<2⟺2D−A正定
SOR 迭代

在Gauss-Seidel迭代的底子上引入松懈因子                                   ω                              \omega                  ω,通过调整                                   ω                              \omega                  ω加速收敛。 迭代格式为Gauss-Seidel效果的加权平均:
                                                    x                               i                                           (                                  k                                  )                                                 =                            (                            1                            −                            ω                            )                                       x                               i                                           (                                  k                                  −                                  1                                  )                                                 +                                       ω                                           a                                               i                                     i                                                                        (                                           b                                  i                                          −                                           ∑                                               j                                     =                                     1                                                           i                                     −                                     1                                                                  a                                               i                                     j                                                                  x                                  j                                               (                                     k                                     )                                                      −                                           ∑                                               j                                     =                                     i                                     +                                     1                                              n                                                      a                                               i                                     j                                                                  x                                  j                                               (                                     k                                     −                                     1                                     )                                                      )                                            x_i^{(k)} = (1-\omega)x_i^{(k-1)} + \frac{\omega}{a_{ii}} \left( b_i - \sum_{j=1}^{i-1} a_{ij}x_j^{(k)} - \sum_{j=i+1}^n a_{ij}x_j^{(k-1)} \right)                     xi(k)​=(1−ω)xi(k−1)​+aii​ω​(bi​−j=1∑i−1​aij​xj(k)​−j=i+1∑n​aij​xj(k−1)​)
其矩阵分裂形式如下


  • 分裂矩阵
                                                  A                               =                                           1                                  ω                                          D                               −                               L                               −                                           (                                                             1                                        −                                        ω                                                  ω                                              D                                  +                                  U                                  )                                                 A = \frac{1}{\omega}D - L - \left( \frac{1-\omega}{\omega}D + U \right)                        A=ω1​D−L−(ω1−ω​D+U)
    于是对应迭代矩阵:
                                                              B                                  SOR                                          =                               (                               D                               −                               ω                               L                                           )                                               −                                     1                                                                  [                                  (                                  1                                  −                                  ω                                  )                                  D                                  +                                  ω                                  U                                  ]                                                 B_{\text{SOR}} = (D - \omega L)^{-1} \left[ (1-\omega)D + \omega U \right]                        BSOR​=(D−ωL)−1[(1−ω)D+ωU]
收敛性定理


  • 须要条件
                                                       0                                  <                                  ω                                  <                                  2                                          0 < \omega < 2                           0<ω<2
    证明可以通过计算:
                                                       det                                  ⁡                                               B                                                   S                                        O                                        R                                                           =                                               1                                                   det                                        ⁡                                        D                                                           ⋅                                  det                                  ⁡                                  D                                                             (                                                       1                                           −                                           w                                                      )                                                  n                                              =                                                             (                                                       1                                           −                                           w                                                      )                                                  n                                                      \det B_{SOR} = \frac{1}{\det D}\cdot\det D\left({1-w}\right)^{n} = \left({1-w}\right)^{n}                           detBSOR​=detD1​⋅detD(1−w)n=(1−w)n
    如果                                             w                               ∉                               (                               0                               ,                               2                               )                                      w \not \in (0,2)                        w∈(0,2),则行列式必然大于1,于是必然存在大于1的特征值,从而不收敛。但是也应该明白如果行列式小于1,也不能包管收敛。
  • 对称正定矩阵的收敛性(Ostrowski定理)
    若                                             A                                      A                        A对称正定且                                             0                               <                               ω                               <                               2                                      0 < \omega < 2                        0<ω<2,则SOR方法收敛。
    这个定理的证明比较复杂,这里略去。
  • 最优松懈因子
    当                                             A                                      A                        A为对称正定的三对角矩阵时,最优松懈因子为:
                                                                    ω                                     opt                                              =                                               2                                                   1                                        +                                                                       1                                              −                                              ρ                                              (                                                               B                                                 J                                                                               )                                                 2                                                                                                                \omega_{\text{opt}} = \frac{2}{1 + \sqrt{1 - \rho(B_J)^2}}                           ωopt​=1+1−ρ(BJ​)2                      ​2​
    此中                                             ρ                               (                                           B                                  J                                          )                                      \rho(B_J)                        ρ(BJ​)为Jacobi迭代矩阵的谱半径。

免责声明:如果侵犯了您的权益,请联系站长,我们会及时删除侵权内容,谢谢合作!更多信息从访问主页:qidao123.com:ToB企服之家,中国第一个企服评测及商务社交产业平台。
回复

使用道具 举报

0 个回复

倒序浏览

快速回复

您需要登录后才可以回帖 登录 or 立即注册

本版积分规则

大号在练葵花宝典

金牌会员
这个人很懒什么都没写!
快速回复 返回顶部 返回列表