矩陣的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