LU分解法求解线性方程组

算法名称：LU分解法

算法描述：

假定能够把矩阵A写成两个矩阵相乘的形式

                                                     （1）

其中L为下三角矩阵，U为上三角矩阵。

例如，4×4矩阵A的情况，（1）式如下：

 （2）

可以用如（1）式分解来求解线性方程组

                                  （3）

首先求解向量y使得

                                                                               （4）

然后再来求解

                                                      （5）

此拆分方法的优点在于求解一个三角形方程组相当容易，这样，（4）式可用向前替代过程求解，如下：

                                                （6）

（5）式可用回代过程求解，这与（2）式~（3）式一样，

                        （7）

（6）式和（7）式共需执行N²次内层循环（对每个右端项b），每个内层循环包括一次乘法和一次加法。如果有N个右端项，它们是单位列向量（在求矩阵逆时就是这种情况），考虑这些零元素可把（6）式的总执行次数从N³/2减少到N³/6，而（7）式的执行次数不变，仍为N³/2。

       注意：一点对A进行了LU分解，就可以一次求解所有要解的右端项。

 

算法实现：

首先，写出（1）式或（2）式的第i，j分量。它总是一个和式，开始部分形式如下：

和式中的项数依赖于i和j中较小的数。事实上有三种形式：

                         （8，9，10）

显然，（8）~（10）式共有N²个方程，而要求N²+N个未知的α和β（因对角线的未知元素有两套），既然未知数的个数比方程个数多，就人为指定N各位指数，然后再来求解其他的未知数。事实上，总是令

                                                                          （11）

有一个算法称为Crout算法，它仅按某种次序排列方程，就能容易的求出（8）式~（11）式的N²+N各方程中的所有α和β。步骤如下：

设
    ，即（11）式

对每个j=0,1,2,...,N-1进行以下两步：

第一步，对每个i=0,1,...,j用（8）式、（9）式和（11）式来解β_ij，即

                                                                          （12）

第二步，对每个i=j+1,j+2,...,N-1用（10）式来求解α_ij，即

                              （13）

在求解下一个j之前要保证进行了以上两步。

 

如果按上述过程进行几次迭代后，就会发现（12）式和（13）式右端的α和β在需要时已经得到，还会发现，每一个a_ij仅被使用一次就不再使用了。这意味着分解是“同址”进行的。简言之Crout算法得到的矩阵是混合矩阵，对本例排列如下：

注：不是把矩阵A分解成LU形式，而是将其按行置换的方式分解。

 

Crout算法的精妙之处：

l         （12）式，在i=j（最后一次应用）时，与（13）式（除后者还要做一次除法外）是完全一样的，这两种情况要求和的上线都是k=j-1（=i-1）。这意味着，不必费心去考虑对角线元素βjj是否会正落在对角线上，也不必考虑该列中，它下面的某个元素（未做除法的）αij，i=j+1，j+2,...,N-1是否会提升成为对角线元素β。

l         它首先找到每行的最大元素，而后（在找最大主元时）乘以一个比例系数，这就实现了“隐式主元法”。

 

运行示例：

Origin coefficient matrix:

 | 0.0 2.0 0.0 1.0 |

 | 2.0 2.0 3.0 2.0 |

 | 4.0 -3.0 0.0 1.0 |

 | 6.0 1.0 -6.0 -5.0 |

-----------------------------------------------

LU mixed matrix:

 | 6.0 1.0 -6.0 -5.0 |

 | 0.0 2.0 0.0 1.0 |

 | 0.3333333333333333 0.8333333333333334 5.0 2.833333333333333 |