2021-02-27-数值分析-Day02-三角分解

续 2.1 - 2.4

矩阵三角分解法

初等行变换 \(\Leftrightarrow\) 矩阵左乘初等矩阵消元 \(\Leftrightarrow\) 矩阵左乘 \((n-1)\) 个初等矩阵

若 \(a_{11}^{(1)}\neq0\) ，令 \(l_{i1} = a_{i1}^{(1)}\div a_{11}^{11},\ i=2,3,\dots,n\) ，记： \(L_{1}=\begin{pmatrix} 1 \\ -l_{21} & 1 \\ -l_{31} && 1 \\ \vdots &&& \ddots \\ -l_{31} &&&& 1 \end{pmatrix}\)

则有 \(A^{(2)}=L_{1}A^{(1)}=\begin{pmatrix} a_11^{(1)} & a_{12}^{(1)} & \cdots & a_{1n}^{(1)} \\ 0 & a_{22}^{(2)} & \cdots & a_{2n}^{(2)} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & a_{n2}^{(2)} & \cdots & a_{nn}^{(2)} \end{pmatrix}\)

同理进行第二步消元 \(A^{(3)}=L_{2}A^{(2)},\ \dots\) 第 \((n-1)\) 步得到： \(A^{(n)}=L_{n-1}A^{(n-1)}=\begin{pmatrix} a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1n}^{(1)} \\ & a_{22}^{(2)} & \cdots & a_{2n}^{(2)} \\ && \ddots & \vdots \\ &&& a_{nn}^{(n)} \end{pmatrix}\)

其中： \(L_{n-1}=\begin{pmatrix}1\\&1\\&&\ddots\\&&&1\\&&&-l_{n,n-1}&1\end{pmatrix}\)

也就是： \(A^{(n)}=L_{n-1}A^{(n-1)}=L_{n-1}L_{n-2}A^{(n-2)}=\dots=L_{n-1}L_{n-2}\dots L_{2}L_{1}A^{(1)}\)

其中： \(L_{k}=\begin{pmatrix}1\\&\cdots\\&&1\\&&-l_{k+1k}&1\\&&\vdots&&\ddots\\&&-l_{k+(n-1)k}&&&1\end{pmatrix},\ k=1,2,\dots,n-1\)
\(L_{k}^{-1}=\begin{pmatrix}1\\&\cdots\\&&1\\&&l_{k+1k}&1\\&&\vdots&&\ddots\\&&l_{k+(n-1)k}&&&1\end{pmatrix}\)

所以有： \(A=A^{(1)}=L_{1}^{-1}L_{2}^{-1}\dots L_{n-1}^{-1}A^{(n)}=LU\) ，称为 \(A\) 的 \(LU\) 三角分解。

其中： \(L=L_{1}^{-1}L_{2}^{-1}\dots L_{n-1}^{-1}, U=A^{(n)}\) \(L=\begin{pmatrix}1\\l_{21}&1\\l_{31}&l_{32}&1\\\vdots&\vdots&\vdots&\ddots\\l_{n1}&l_{n2}&l_{n3}&\cdots&1\end{pmatrix}\) 即一个单位下三角矩阵

\(U=\begin{pmatrix} a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1n}^{(1)} \\ & a_{22}^{(2)} & \cdots & a_{2n}^{(2)} \\ && \ddots & \vdots \\ &&& a_{nn}^{(n)} \end{pmatrix}\) 为一个上三角矩阵

2.5 - 2.6 直接三角分解法

定理1：设 \(n\) 阶方阵 \(A\) 的各阶顺序主子式不为零，则存在唯一单位下三角矩阵 \(L\) 和上三角矩阵 \(U\) 使得 \(A=LU\) 。（条件与顺序Gauss消去法的执行条件一致。）

由此可得： \(Ax=b \ \Rightarrow\ LUx=b\) ，令 \(Ux=y\) ，得 \(\left\{\begin{aligned}Ly=b\\Ux=y\end{aligned}\right.\) 从而将问题转换为两个三角型方程组的求解。

直接得到 \(L\ U\) 两个矩阵——Doolittle分解法 设： \[ \begin{pmatrix}a_{11}&a_{12}&\cdots&a_{1n}\\a_{21}&a_{22}&\cdots&a_{2n}\\\vdots&\vdots&\ddots&\vdots\\a_{n1}&a_{n2}&\cdots&a_{nn}\end{pmatrix} = \begin{pmatrix}1\\l_{21}&1\\l_{31}&l_{32}&1\\\vdots&\vdots&\vdots&\ddots\\l_{n1}&l_{n2}&l_{n3}&\cdots&1\end{pmatrix} \begin{pmatrix}u_{11}&u_{12}&\cdots&u_{1n}\\&u_{22}&\cdots&u_{2n}\\&&\ddots&\vdots\\&&&u_{nn}\end{pmatrix} \] 则得： \[ \left\{\begin{array}{l} u_{1j}=a_{1j}\quad j=1,2,\dots,n\\ l_{i1}=a_{i1}\div u_{11}\quad i=2,3,\dots,n\\\vdots\\ 对k=2,3,\dots,n,计算\\ u_{kj}=a_{kj}-\sum_{m=1}^{k-1}{l_{km}u_{mj}}\quad j=k,k+1,\dots,n\\ l_{ik}=(a_{ik}-\sum_{m=1}^{k-1}{l_{im}u_{mk}})\div u_{kk}\quad i=k+1,k+2,\dots,n \end{array}\right. \] 应用中为节省存储空间，将矩阵 \(U\) 覆盖矩阵 \(A\) 的上三角进行存放，矩阵 \(L\) 覆盖矩阵 \(A\) 除对角线的下三角进行存放，即在原位形成矩阵： \(\begin{pmatrix}u_{11}&u_{12}&\cdots&u_{1n}\\l_{21}&u_{22}&\cdots&u_{2n}\\\vdots&\vdots&\ddots&\vdots\\l_{n1}&l_{n2}&\cdots&u_{nn}\end{pmatrix}\)

回代得到： \[ \begin{array}{l} \left\{\begin{array}{l} \begin{pmatrix}1\\l_{21}&1\\\vdots&\vdots&\ddots\\l_{n1}&l_{n2}&\cdots&1\end{pmatrix} \begin{pmatrix}y_1\\y_2\\\vdots\\y_n\end{pmatrix} = \begin{pmatrix}b_1\\b_2\\\vdots\\b_n\end{pmatrix},\\ \begin{pmatrix}u_{11}&u_{12}&\cdots&u_{1n}\\&u_{22}&\cdots&u_{2n}\\&&\ddots&\vdots\\&&&u_{nn}\end{pmatrix} \begin{pmatrix}x_1\\x_2\\\vdots\\x_n\end{pmatrix} = \begin{pmatrix}y_1\\y_2\\\vdots\\y_n\end{pmatrix} \end{array}\right.\\ 解得\left\{\begin{array}{l} y_1=b_1\\ y_k=b_k-\sum_{i=1}^{k-1}{l_{ki}y_i},\quad k=2,3,\dots,n\\ x_n=y_n\div u_{nn}\\ x_i=(y_i-\sum_{j=i+1}^{n}{u_{ij}x_j})\div u_{ii},\quad i=n-1,n-2m\dots,1 \end{array}\right. \end{array} \]

2.7 平方根法

针对系数矩阵 \(A\) 为对称正定矩阵情况，有 \(A=LU,u_{kk}>0\) ，而：

\(\begin{pmatrix}u_{11}&u_{12}&\cdots&u_{1n}\\&u_{22}&\cdots&u_{2n}\\&&\ddots&\vdots\\&&&u_{nn}\end{pmatrix} = \begin{pmatrix}u_{11}\\&u_{22}\\&&\ddots\\&&&u_{nn}\end{pmatrix} \begin{pmatrix}1&u_{12}&\cdots&u_{1n}\\&1&\cdots&u_{2n}\\&&\ddots&\vdots\\&&&1\end{pmatrix}\)

令 \(D=\begin{pmatrix}u_{11}\\&u_{22}\\&&\ddots\\&&&u_{nn}\end{pmatrix},\ M=\begin{pmatrix}1&u_{12}&\cdots&u_{1n}\\&1&\cdots&u_{2n}\\&&\ddots&\vdots\\&&&1\end{pmatrix}\)

则有： \(A=LDM\) ，又因为 \((LDM)^T=M^TDL^T=LDM\) ，所以： \(M=L^T\) ，则有： \(A=LDM=LDL^T\) 。

令 \(D=\begin{pmatrix}\sqrt{u_{11}}\\&\sqrt{u_{22}}\\&&\ddots\\&&&\sqrt{u_{nn}}\end{pmatrix}^2=(D^{\frac{1}{2}})^2\)

则有： \(A=LD^{\frac{1}{2}}D^{\frac{1}{2}}L^T=(LD^{\frac{1}{2}})(LD^{\frac{1}{2}})^T=GG^T,\quad G=LD^{\frac{1}{2}}\)

分解 \(A=GG^T\) 称为对称正定矩阵的Cholesky分解， \(Ax=b\) 则转换为 \(\left\{\begin{aligned}Gy=b\\G^Tx=y\end{aligned}\right.\) ，称为平方根法。

求解方法：

若记 \(G=(g_{ij})\) ，则对 \(k=1,2,\dots,n\) 有： \(\left\{\begin{array}{l}g_{kk}=(a_{kk}-\sum_{m=1}^{k-1}{g_{km}^2})^\frac{1}{2}\\g_{ik}=(a_{ik}-\sum_{m=1}^{k-1}{g_{im}g_{km}})\div g_{kk},\ i=k+1,\dots,n\end{array}\right.\)

同样将 \(G\) 存放在 \(A\) 的位置来节省存储。

解三角方程 \(\left\{\begin{aligned}Gy=b\\G^Tx=y\end{aligned}\right.\) 得： \(\left\{\begin{array}{l}y_k=(b_k-\sum_{m=1}^{k-1}{g_{km}y_m})\div g_{kk},\quad k=1,2,\dots,n\\x_k=(y_k-\sum_{m=k+1}^{n}g_{mk})\div g_{kk},\quad k=n,n-1,\dots,1\end{array}\right.\)

2.8 追赶法

追赶法是求三对角线性方程组的三角分解法： \(\begin{pmatrix}a_1&c_1\\d_2&a_2&c_2\\&\ddots&\ddots&\ddots\\&&d_{n-1}&a_{n-1}&c_{n-1}\\&&&d_n&a_n\end{pmatrix} \begin{pmatrix}x_1\\x_2\\\vdots\\x_{n-1}\\x_n\end{pmatrix} = \begin{pmatrix}b_1\\b_2\\\vdots\\b_{n-1}\\b_n\end{pmatrix}\)

三对角矩阵 \(A\) 的各阶顺序主子式都不为零的一个充分条件是： \(|a_1|>|c_1|>0;\ |a_n|>|d_n|>0;\ |a_i|\ge|c_i|+|d_i|, c_id_i\neq0, i=2,3,\dots,n-1\) 在此条件下， \(A=LDM=TM\) ，称为矩阵 \(A\) 的Crout分解。

对三对角矩阵 \(A\) 进行Crout分解，有：\(A=\begin{pmatrix}\alpha_1\\\gamma_2&\alpha_2\\&\ddots&\ddots\\&&\gamma_{n-1}&\alpha_{n-1}\\&&&\gamma_n&\alpha_n\end{pmatrix} \begin{pmatrix}1&\beta_1\\&1&\beta_2\\&&\ddots&\ddots\\&&&1&\beta_{n-1}\\&&&&1\end{pmatrix}\)
\(\left\{\begin{array}{l}\alpha_1=a_1,\ \beta_1=c_1\div\alpha_1,\ \gamma_i=d_i,\ i=1,2,3,\dots,n\\ \alpha_i=a_i-d_i\beta_{i-1},\ i=2,3,\dots,n\\ \beta_i=c_i\div\alpha_i,\ i=2,3,\dots,n-1 \end{array}\right.\)

代码：LU分解法（Doolittle）

//
// Created by xa on 2021-02-27.
//
#include <iostream>
#include <vector>

using std::vector;

vector<vector<double>> A;
vector<double> b;
int size = 0;

vector<double> solve(vector<vector<double>> A, vector<double> b);
void getLU();

void printA();

int main(){
    vector<vector<double>> A_ {{8, 1, 7}, {3, 7, 9}, {9, 1, 5}};
    vector<double> b_ {3, 8, 6};
    vector<double> x = solve(A_, b_);

    for (int i = 0; i < size; ++i)
        std::cout << x[i] << " ";
}

vector<double> solve(vector<vector<double>> A_, vector<double> b_) {
    A.assign(A_.begin(), A_.end());
    b.assign(b_.begin(), b_.end());
    size = A.size();
    vector<double> x(size);
    vector<double> y(size);

    getLU();

    // solve y
    y[0] = b[0];
    for (int k = 1; k < size; ++k) {
        y[k] = b[k];
        for (int i = 0; i < k; ++i)
            y[k] -= A[k][i] * y[i];
    }

    // solve x
    x[size-1] = y[size-1] / A[size-1][size-1];
    for (int i = size-2; i >= 0; --i) {
        double sum = 0;
        for (int j = i+1; j < size; ++j)
            sum += A[i][j] * x[j];
        x[i] = (y[i] - sum) / A[i][i];
    }
    return x;
}

void getLU() {
    for (int k = 0; k < size; ++k) {
        if (k > 0) for (int j = k; j < size; ++j)
                for (int m = 0; m < k; ++m) {
                    A[k][j] -= A[k][m] * A[m][j];
                }
        if (k < size-1) for (int i = k+1; i < size; ++i) {
            double sum = 0;
            for (int m = 0; m < k; ++m)
                sum += A[i][m] * A[m][k];
            A[i][k] = (A[i][k] - sum) / A[k][k];
        }
    }
    printA();
}

void printA() {
    for (int i = 0; i < size; ++i) {
        for (int j = 0; j < size; ++j)
            std::cout << A[i][j] << ' ';
        std::cout << std::endl;
    }
    std::cout << std::endl;
}

续 2.1 - 2.4¶

矩阵三角分解法¶

2.5 - 2.6 直接三角分解法¶

2.7 平方根法¶

2.8 追赶法¶

代码：LU分解法（Doolittle）¶