数值分析 解线性方程组的直接法(二)

直接三角分解法及列主元三角分解法

    一个方阵可分解成两个以上矩阵的乘积，当给定一些约束后，这些分解将是唯一的.对矩阵进行三角分解实际上只是将矩阵分解成两个三角形矩阵的乘积，这样对任意方程组的求解问题就转化为三角形方程组的求解问题,方便计算,
定义1对于n阶矩阵$A$
$$A=\left(\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ \cdots & \cdots & \cdots & \cdots \\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right)$$
若存在n阶下三角方阵L和上三角方阵U,使得A=LU﹐则称方阵A有三角分解或LU分解.特别地，若L为单位下三角矩阵，则称它为 Doolittle分解.若U为单位上三角矩阵，则称它为Grout 分解.本文只介绍Doolittle三角分解.
定理2（存在唯一性定理)设$A$为$n$阶方阵,若$A$的所有顺序主子式${\Delta }_k,\ne 0(k=1,2,.\dots ,n)$,则$A$可以分解为一个单位下三角矩阵L和一个上三角矩阵U的乘积，且这种分解是唯一的.

LU分解法原理

如果把非奇异矩阵$A$分解成一个单位下三角矩阵$L$和一个上三角矩阵$U$的乘积
$$A=LU$$
则写成矩阵元素如下

$$\left(\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ \cdots & \cdots & \cdots & \cdots \\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right)=\left(\begin{array}{cccc}1& & & \\ l_{21}& 1& & \\ ⋮& ⋮& \ddots & \\ l_{n1}& l_{n2}& \cdots & 1\end{array}\right)\left(\begin{array}{cccc}{u}_{11}& u_{12}& \cdots & u_{1n}\\ & u_{22}& \cdots & u_{2n}\\ & & \ddots & ⋮\\ & & & u_{nn}\end{array}\right)$$
式中矩阵$L$与$U$的元素待定.根据矩阵乘法法则,对$LU$三角分解式分析如下
    首先比较等式两边第1行元素,可得矩阵$U$的第1行元素
$$u_{1j}=a_{1j},k=1,2,\cdots ,n$$
然后考察等式两边第1列元素，有
$$a_{i1}=l_{i1}{u}_{11},i=1,2,\cdots ,n$$
可得$L$的第1列元素
$$l_{i1}=\frac{a_{i1}}{u_{11}},i=1,2,\cdots ,n$$
接着比较等式两边，可得矩阵$A$的元素$a$等于L的第$k$行与$U$的第j列对应元素乘积之和,因此有
$$a_{kj}=\sum _{q=1}^n{l}_{kq}{u}_{qj}=\sum _{q=1}^{k-1}{l}_{kq}{u}_{qj}+u_{kj},j=k,k+1,\cdots ,n$$
转换后得$U$的第$k$行元素
$$u_{kj}=a_{kj}-\sum _{q=1}^{k-1}{l}_{kq}{u}_{qj},j=k,k+1,\cdots ,n\text{\hspace{1em}(3.10)}$$
同理，由
$$a_{ik}=\sum _{q=1}^n{l}_{iq}{u}_{qk}=\sum _{q=1}^{k-1}{l}_{iq}{u}_{qk}+l_{ik}{u}_{kk},i=k+1,k+2,\cdots ,n$$
可得与$L$的第$k$列元素
$$l_{ik}=\frac{a_{ik}-\sum _{q=1}^{k-1}{l}_{iq}{u}_{qk}}{u_{kk}},i=k+1,k+2,\cdots ,n\text{\hspace{1em}(3.11)}$$

最后交替使用式$(3.10)$和式$(3.11)$，可以逐步求出$U$与$L$的元素.
    直接三角分解法就是将方程组
$$Ax=b$$
转换成
$$LUx=b$$
令$Ux=y$，则$Ly=b$,即有
$$\left(\begin{array}{cccc}{u}_{11}& u_{12}& \cdots & u_{1n}\\ & u_{22}& \cdots & u_{2n}\\ & & \ddots & ⋮\\ & & & u_{nn}\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right)=\left(\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_n\end{array}\right)$$
$$\left(\begin{array}{cccc}1& & & \\ l_{21}& 1& & \\ ⋮& ⋮& \ddots & \\ l_{n1}& l_{n2}& \cdots & 1\end{array}\right)\left(\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_n\end{array}\right)=\left(\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right)$$
故有
$$\{\begin{array}{l}{y}_1=b_1\\ y_i=b_i-\sum _{j=1}^{i-1}{l}_{ij}{y}_j,i=2,3,\cdots ,n\end{array}\text{\hspace{1em}(3.12)}$$
$$\{\begin{array}{l}{x}_n=y_n/u_{nn}\\ x_i=(y_i-\sum _{j=i+1}^n{u}_{ij}{y}_j)/u_{ii},i=n-1,n-2,\cdots ,1\end{array}\text{\hspace{1em}(3.13)}$$
说明 计算u与的过程是先行后列交替进行的，即依次先计算出$u_{1j},(j=1,2,\dots ,n)$，再算出$L_{i1}(i=2,3,\dots ,n)$;先算出$u_{2j}(j=2,3,\dots ,n)$，再算出$l_{i_2},(i=3,4,\dots ,n);$其余类推，即可计算出$L、U$的各元素.

三角分解法的代码实例

def Trig_decomposition(coef_mtrix, b_mtix):
    n = coef_mtrix.shape[0]  # 获取方程组的行数和列数
    L = np.eye(n)  # 初始化矩阵L
    U = np.zeros([n, n])  # 初始化矩阵U
    U[0] = coef_mtrix[0]
    L[1:, 0] = coef_mtrix[1:, 0] / U[0, 0]
    for k in range(1, n):
        # 计算L、U矩阵
        for j in range(k, n):
            U[k, j] = coef_mtrix[k, j] - np.dot(L[k], U[:, j])
        for i in range(k + 1, n):
            L[i, k] = (coef_mtrix[i, k] - np.dot(L[i], U[:, k])) / U[k, k]

    y = np.zeros([n, 1])  # 初始化矩阵y
    for i in range(n):  # 计算矩阵y
        if i == 0:
            y[i] = b_mtix[i]
        else:
            y[i] = b_mtix[i] - np.dot(L[i, :i], y[:i])

    x = np.zeros([n, 1])  # 初始化矩阵X
    for i in range(n - 1, -1, -1): # 求解矩阵X
        x[i, 0] = (y[i, 0] - np.dot(U[i, i + 1:], x[i + 1:, 0])) / U[i, i]
    return x

需要注意的是公式3.11和3.10中的求和刚好对应的是L和U中的某一行或某一列的元素或部分元素，所以我们可以利用numpy数组切片的方式将对应矩阵部分取出并点积刚好对应公式中的元素求和。
同样的，求解矩阵y和矩阵x我们也利用了这个特点