矩阵的lu分解算法

Ⅰ 矩阵A 的LU 分解

将矩阵A分解为单位下三角阵L和上三角阵U，就是找出L的元素l_ij以及U的元素u_ij与A的元素a_ij的关系式。下面说明矩阵L，U的元素可以通过n步由矩阵A的元素计算出来，其中第k步求出U的第k行和L的第k列元素。

根据式（4－1）中，由矩阵乘法可知：

（1）由L的第1行和U的第j列元素对应相乘后与A的对应元素相等，即

a_1j＝u_1j

同理可知U的第1行元素：

u_1j＝a_1j （j＝1，2，…，n）

由L的第i行和U的第1列元素对应相乘后与A的对应元素相等，即

a_i1＝l_i1·u₁₁

从而得到L的第一列元素

l_i1＝a_i1/u₁₁ （i＝2，3，…，n）

这样就确定了U的第1行元素及L的第1列元素。类似地，用矩阵乘法再确定U的第2行及L的第2列，如此继续。

假设已经得出了U的前k－1行及L的前k－1列（1≤k≤n）的全部元素，现在来确定U的第k行元素和L的第k列元素。

（2）在式（4－1）中，由矩阵乘法知

地球物理数据处理基础

注意，由于L是单位下三角阵，当r>k时，l_kr＝0，而l_kk＝1，所以

地球物理数据处理基础

从而有

地球物理数据处理基础

如此即可得到上三角矩阵U的第k行元素。

同理，可确定L的第k列元素。

由 并且当r>k时，u_rk＝0，则有

地球物理数据处理基础

从而有

地球物理数据处理基础

即可计算出矩阵L的第k列元素。

上述步骤进行，就可以定出L及U的全部元素，完成矩阵A的LU分解，即对k＝1，2，…，n，计算A＝LU分解的公式为

地球物理数据处理基础

这里约定

上述这种矩阵A的LU分解计算顺序也可按图4－1所示逐步进行。

由于以上计算公式（4－7）中不含消去法的中间结果a^（k）_ij，可直接逐框计算，所以称为紧凑格式。

从计算过程中可以看出，对A进行LU分解时，在求出u_1j（j＝1，2，…，n）后不再需要a_1j（j＝1，2，…，n），求出l_i1（i＝2，3，…，n）后不再需要a_i1（i＝2，3，…，n）。同样求出u_kj和l_ik后不再需要a_kj和a_ik（j＝k，k＋1，…，n；i＝k＋1，k＋2，…，n；k＝1，2，…，n），因此，在计算过程中可以不另设存放L和U的数组，而是将u_kj和l_ik算出后直接存入矩阵A对应元素a_kj和a_ik的存贮单元中。这也称为动态算法，它可节省大量内存。

图4－1 矩阵LU分解顺序图

对A进行LU分解后，可按式（4－5）和式（4－6）求解线性方程组。

同样，求解线性方程组时，由式（4－5）和式（4－6）式可知，在求yi（i＝1，2，…，n）时，只用到对应的b_i（i＝1，2，…，n），故可不另设存放y的数组，而将y的元素存入列向量b对应元素b_i的存贮单元中。但应注意，这时的b已变成y。同理求x时，也不用另设存放单元，计算结果直接存入b中，因而在求解线性方程组的最后结果中，b存储的是x。

另外，由于求yi的顺代公式（4－5）与求u_kj的公式（4－7）形式完全相同，故可合并到k循环之内进行，这样可以简化程序结构。所以，在进行LU分解的同时，Ly＝b的求解也完成了。

［例1］用直接三角分解法求解线性方程组

解：（1）分解A＝LU，

地球物理数据处理基础

即

地球物理数据处理基础

（2）求解Ly＝b，

地球物理数据处理基础

得到y＝（3，－5，6）^T

（3）求解

得到x＝（2，－2，1）^T