1.矩阵三角分解原理

设AAA为非奇异矩阵，且有分解式A=LU,A=LU,A=LU,
A=[a11a12…a1na21a21…a2n⋮⋮⋮an1an1…ann]A=\begin{bmatrix} a_{11}&a_{12}&\ldots&a_{1n}\\ a_{21}&a_{21}&\ldots&a_{2n}\\ \vdots&\vdots&&\vdots\\ a_{n1}&a_{n1}&\ldots&a_{nn} \end{bmatrix}A=⎣⎢⎢⎢⎡a11a21⋮an1a12a21⋮an1………a1na2n⋮ann⎦⎥⎥⎥⎤
其中LLL为单位下三角矩阵,UUU为上三角矩阵,即
A=[1l211⋮⋮⋱ln1ln2…1][u11u12…u1nu22…u2n⋱⋮unn]A=\begin{bmatrix} 1&&&\\ l_{21}&1&&\\ \vdots&\vdots&\ddots&\\ l_{n1}&l_{n2}&\ldots&1 \end{bmatrix} \begin{bmatrix} u_{11}&u_{12}&\ldots&u_{1n}\\ &u_{22}&\ldots&u_{2n}\\ &&\ddots&\vdots\\ &&&u_{nn} \end{bmatrix}A=⎣⎢⎢⎢⎡1l21⋮ln11⋮ln2⋱…1⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡u11u12u22……⋱u1nu2n⋮unn⎦⎥⎥⎥⎤
下面说明L,UL,UL,U的元素可以由nnn步直接定出,其中第rrr步定出UUU的第rrr行和LLL的第rrr列.由上述分解式有a1i=u1i,i=1,2,…,n,a_{1i}=u_{1i},i=1,2,\ldots,n,a1i=u1i,i=1,2,…,n,
得到UUU的第111行元素;
ai1=li1u11,li1=ai1/u11,i=2,3,…,n,a_{i1}=l_{i1}u_{11},l_{i1}=a_{i1}/u_{11},i=2,3,\ldots,n,ai1=li1u11,li1=ai1/u11,i=2,3,…,n,
得到LLL的第111列元素.
设已经定出UUU的第1行到第r−1r-1r−1列元素.有上述分解式，利用矩阵乘法(注意当r<k,lrk=0r <k,l_{rk}=0r<k,lrk=0),有ari=∑k=1nlrkuki=∑k=1r−1lrkuki+uri,a_{ri}=\sum_{k = 1}^{n}l_{rk}u_{ki}=\sum_{k = 1}^{r-1}l_{rk}u_{ki}+u_{ri},ari=k=1∑nlrkuki=k=1∑r−1lrkuki+uri,
故uri=ari−∑k=1r−1lrkuki,i=r,r+1,…,n.u_{ri}=a_{ri}-\sum_{k = 1}^{r-1}l_{rk}u_{ki},i=r,r+1,\ldots,n.uri=ari−k=1∑r−1lrkuki,i=r,r+1,…,n.
又由该分解式有air=∑k=1nlikukr=∑k=1r−1likukr+lirurr.a_{ir}=\sum_{k = 1}^{n}l_{ik}u_{kr}=\sum_{k = 1}^{r-1}l_{ik}u_{kr}+l_{ir}u_{rr}.air=k=1∑nlikukr=k=1∑r−1likukr+lirurr.
总结上述讨论，得到直接三角分解法解Ax=bAx=bAx=b(要求AAA的所有顺序主子式都不为零)的计算公式.

u1i=a1i(i=1,2,…,n),li1=ai1/u11,i=2,3,…,n.u_{1i}=a_{1i}(i=1,2,\ldots,n),l_{i1}=a_{i1}/u_{11},i=2,3,\ldots,n.u1i=a1i(i=1,2,…,n),li1=ai1/u11,i=2,3,…,n.
计算UUU的第rrr行,LLL的第rrr列元素(r=2,3,…,n):(r=2,3,\ldots,n):(r=2,3,…,n):
uri=ari−∑k=1r−1lrkuki,i=r,r+1,…,n;u_{ri}=a_{ri}-\sum\limits_{k = 1}^{r-1}l_{rk}u_{ki},i=r,r+1,\ldots,n;uri=ari−k=1∑r−1lrkuki,i=r,r+1,…,n;
lir=(air−∑k=1r−1likukr)/urr,i=r+1,…,n,且r≠n.l_{ir}=(a_{ir}-\sum\limits_{k = 1}^{r-1}l_{ik}u_{kr})/u_{rr},i=r+1,\ldots,n,\text{且}r \ne n.lir=(air−k=1∑r−1likukr)/urr,i=r+1,…,n,且r=n.
求解Ly=b,Ux=yLy=b,Ux=yLy=b,Ux=y的计算公式:
{y1=b1,yi=bi−∑k=1i−1likyk,i=2,3,…,n;\begin{cases} y_1=b_1,\\ y_i=b_i-\sum\limits_{k = 1}^{i-1}l_{ik}y_k,i=2,3,\ldots,n; \end{cases}⎩⎨⎧y1=b1,yi=bi−k=1∑i−1likyk,i=2,3,…,n;
{xn=yn/unn,xi=(yi−∑k=i+1nuikxk)/uii,i=n−1,n−2,…,1.\begin{cases} x_n=y_n/u_{nn},\\ x_i=(y_i-\sum\limits_{k=i+1}^{n}u_{ik}x_k)/u_{ii},i=n-1,n-2,\ldots,1. \end{cases}⎩⎨⎧xn=yn/unn,xi=(yi−k=i+1∑nuikxk)/uii,i=n−1,n−2,…,1.

2.Python实现不选主元矩阵三角分解

# 自己原创
def triangle_decomposition(coefficient_matrix: np.ndarray):"""实现方程Ax=b系数矩阵A的Doolittle decomposition:param coefficient_matrix: 初始系数矩阵A:return: 单位下三角矩阵L,上三角矩阵U"""# first step: evaluate the decomposition conditionrows, columns = coefficient_matrix.shapeif rows == columns:  # judge if it is a square matrixfor k in range(rows):  # 判断各阶顺序主子式是否为0if det(coefficient_matrix[:k + 1, :k + 1]) == 0:raise Exception("cannot decompose ")else:# directive triangle  decomposition# 初始化square_ones_matrix = np.ones((rows, columns))lower_triangle_matrix = np.tril(square_ones_matrix)upper_triangle_matrix = np.triu(square_ones_matrix)# 计算第1行的U(1,i),第1列的L(i,1)upper_triangle_matrix[0, :] = coefficient_matrix[0, :]lower_triangle_matrix[1:, 0] = coefficient_matrix[1:, 0] / upper_triangle_matrix[0, 0]# 计算第row行的U(row,i),第row列的L(i,row)for row in range(1, rows - 1):for i in range(row, rows):# 定义一个累积和变量row_i_prod = 0for k in range(row):row_i_prod += lower_triangle_matrix[row, k] * upper_triangle_matrix[k, i]upper_triangle_matrix[row, i] = coefficient_matrix[row, i] - row_i_prodfor i in range(row + 1, rows):# 定义一个累积和变量i_row_prod = 0for k in range(row):i_row_prod += lower_triangle_matrix[i, k] * upper_triangle_matrix[k, row]lower_triangle_matrix[i, row] = (coefficient_matrix[i, row] - i_row_prod) / upper_triangle_matrix[row, row]# 定义一个累积和变量,单独计算n行n列的U(n,n)row_i_prod = 0for k in range(rows - 1):row_i_prod += lower_triangle_matrix[rows - 1, k] * upper_triangle_matrix[k, rows - 1]upper_triangle_matrix[rows - 1, rows - 1] = coefficient_matrix[rows - 1, rows - 1] - row_i_prodreturn lower_triangle_matrix, upper_triangle_matrixelse:raise Exception("ERROR:please pass a square matrix.")
# 自己原创
def directive_triangle_decomposition(coefficient_matrix: np.ndarray, right_hand_side_vector: np.ndarray,decomposition_function=triangle_decomposition):"""实现方程Ax=b系数矩阵A directive triangle decomposition:param decomposition_function: 三角分解函数:param coefficient_matrix: 初始系数矩阵A:param right_hand_side_vector: 初始常数列向量b:return: 解向量x"""l, u = decomposition_function(coefficient_matrix)rows, columns = coefficient_matrix.shape# Ly=b，此处不指定双精度浮点数据类型，则可能会损失精度，导致结果偏差太大y = np.ones_like(right_hand_side_vector, dtype=np.float64)y[0, 0] = right_hand_side_vector[0, 0]for row in range(1, rows):# 定义一个累积和变量i_k_prod = 0for k in range(row):i_k_prod += l[row, k] * y[k, 0]y[row, 0] = right_hand_side_vector[row, 0] - i_k_prod# Ux=y，此处不指定双精度浮点数据类型，则可能会损失精度，导致结果偏差太大x = np.ones_like(right_hand_side_vector, dtype=np.float64)x[rows - 1] = y[rows - 1] / u[rows - 1, rows - 1]for row in range(rows - 2, 0, -1):i_k_prod = 0for k in range(row + 1, rows):i_k_prod += u[row, k] * x[k, 0]x[row, 0] = (y[row, 0] - i_k_prod) / u[row, row]return x

3.解方程

用直接三角分解法解[123252315][x1x2x3]=[141820]\begin{bmatrix} 1&2&3\\ 2&5&2\\ 3&1&5 \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ x_3 \end{bmatrix}= \begin{bmatrix} 14\\ 18\\ 20 \end{bmatrix}⎣⎡123251325⎦⎤⎣⎡x1x2x3⎦⎤=⎣⎡141820⎦⎤

4.测试

if __name__ == '__main__':
# 直接三角分解法-不选主元三角分解法测试成功，来源详见李庆扬数值分析第5版P153，e.g.5square_matrix = np.array([[1, 2, 3], [2, 5, 2], [3, 1, 5]])column_vector = np.array([14, 18, 20]).reshape((3, 1))print(directive_triangle_decomposition(square_matrix, column_vector))

5.运行结果

不选主元的矩阵三角分解法相关推荐

解线性方程组的python实现(2)——矩阵三角分解法
解线性方程组的python实现2--矩阵三角分解法 1. 矩阵三角分解法实现代码 2. LU分解 2.1 基本步骤 2.2 LU分解的计算公式 2.3 LU分解的结果表示实现代码 3. 选主元的L ...
矩阵三角分解matlab,4矩阵三角分解法.ppt
* 一.直接法概述直接法是将原方程组化为一个或若干个三角形方程组的方法,共有若干种．对于线性方程组其中系数矩阵未知量向量常数项根据Cramer(克莱姆)法则,若 determinant ...
Guass列选主元消去法和三角分解法
最近数值计算学了Guass列主消元法和三角分解法解线性方程组,具体原理如下: 1.Guass列选主元消去法对于AX =B 1).消元过程:将(A|B)进行变换为,其中是上三角矩阵.即: k从1到n-1 ...
2.3 数值分析: 矩阵三角分解法
本文内容为东北大学数值分析国家精品慕课课程的课程讲义, 将其整理为OneNote笔记同时添加了本人上课时的课堂笔记, 且主页中的思维导图就是根据课件内容整理而来, 为了方便大家和自己查看,特将此上传到 ...
Python02 雅克比迭代法 Gauss-Seidel迭代法列选主元法 LU分解法(附代码)
1. 实验结果 (1)在定义的矩阵类中设置需要求解的方程为: (2)在 test.py 中选择雅克比迭代法求解: 输入:最大容许迭代次数和精度要求: 输出:根据谱半径判断方法是否收敛,收敛时得到满足精 ...
doolittle分解matlab,如何理解选主元的Doolittle分解法
书中讲解不是很详细,理解之后总结一下首先说一下,之所以要理解选主元的Doolittle分解是因为书中对于该分解过程的讲解比较违和.本文的目的是为了说明:选主元的Doolittle分解法分解得到的LU ...
Matlab实现 LU分解法解线性方程组（全选主元列选主元）
选主元LU分解实验内容:列选主元LU分解和全选主元LU分解求解线性方程组计算方法: 全选主元消元法 1.1 初始化根据参数A.b,记录下矩阵.右端项的尺寸n: 以得到的尺寸n初始化解向量x: 同 ...
lu分解法matlab_MIT 18.065—机器学习中的矩阵方法02 矩阵乘法与矩阵分解
数据分析.信号处理和机器学习中的矩阵方法第02讲矩阵乘法与矩阵分解新MIT 线性代数|机器学习(中英机翻字幕)18.065 by Gilbert Strang_哔哩哔哩 (゜-゜)つロ干杯~- ...
产品需求分析与市场分析方法汇总（SWOT+PDCA+波士顿矩阵BCG+5W2H分析法+STAR关键事件分析法+目标管理SMART+时间管理紧急重要矩阵+WBS任务分解法）
产品需求分析与市场分析方法汇总(SWOT+PDCA+波士顿矩阵BCG+5W2H分析法+STAR关键事件分析法+目标管理SMART+时间管理紧急重要矩阵+WBS任务分解法) 产品需求分析与市场分析方法汇 ...

不选主元的矩阵三角分解法