LDLT分解法求解对称线性方程组.docVIP

LDLT 分解法求解线性方程组 南京师范大学 - PAGE 4 - LDLT 分解法求解对称线性方程组 朱松盛 041002045 南京师范大学 LDLT 分解法原理 利用矩阵 A 的分解来解Ax = b，的方法称为分解法。当求解方程组的系数矩阵是对称矩阵时，则用分解法可以简化程序设计并减少计算量。 当矩阵 A 的各阶顺序主子式不为零时，A 有唯一的 Doolittle 分解 A = LU，其中 ，。 此时，当然有，所以矩阵U 的对角线元素，（），将矩阵U 每行依次提出，则有，其中 ， 因而，显然这种A 的分解也是唯一的，当A 为对称矩阵时，由，得，由分解的唯一性，有，即。 由此可得，若对称矩阵A 的各阶顺序主子式不为零时，A 可唯一分解为，其中 ，，为L 的转置矩阵。 当A 有分解时，利用矩阵运算法则及相等原理易得计算及的公式为 ， 为减少乘法次数，引入辅助量，则上面公式可写成 ， （1） 当设计程序时，为了减少存储空间，又由于矩阵A 在经过计算后就不需要再使用了，因此可以采用原位存储的方法。也就是可以将存于矩阵A 的对角线元素即中，存于矩阵A 的上三角部分即中，而则可以存于矩阵A 的下三角部分即中。 当需要求解对称线性方程组时，由分解，有，这可以分解成三个简单的方程组，和由此易得求解公式为 （2） 同样，向量b 在经过计算后也不再需要，因此也可以采用相同的方法处理。可以在计算时将存于中，然后计算时也将存于中，最后计算时还是存于中，由输出计算结果即可。 公式（1）和（2）就是分解的求解公式。 用分解法求解对称方程组的计算量约为，这比用Gauss 消元法求解的计算量要少。对一般的对称方程组，会出现因某个而不能进行分解的情况，不过当矩阵A 的各阶顺序主子式都不为零时，则分解总能进行到底，从而可用分解法来求解。 LDLT 分解法程序主要部分（C 语言） 程序中所有矩阵的下标都是从0～n－1，与公式中的1～n不同，这只是习惯上的不同，其余都是一致的。 程序主要部分有两个模块LDLTDCMP 和LDLTBKSB。第一个是分解，第二个是解方程组，A 经LDLTDCMP计算后输出给LDLTBKSB 使用。两个模块分开，可用于解不同b，相同A 的多个对称线性方程组。 void LDLTDCMP (unsigned int n, double * *a) { int k; int m; int i; for (k = 0; k n; k++) { /* Calculate d[k], and saved in a[k][k] */ for (m = 0; m k; m++) a[k][k] = a[k][k] - a[m][k] * a[k][m]; if (a[k][k] == 0) { printf (\n\nERROR: LDL\ decompose failed !!\n); printf ( Main submatrix of every rank of A cannot be zero.\n\n); return; } for (i = k + 1; i n; i++) { /* temporary varible u[i][k] is saved in a[k][i]*/ for (m = 0; m k; m++) a[k][i] = a[k][i] - a[m][i] * a[k][m]; /* Calculate l[i][k], and saved in a[i][k]*/ a[i][k] = a[k][i] / a[k][k]; } } } void LDLTBKSB (unsigned int n, double * *a, double *b) { int i; int k; for (i = 0; i n; i++) { /* Calculate y[i], and saved in b[i] */ for (k = 0; k i; k++) b[i] =