平方根法计求解线性方程组

—Cholesky方法

解n阶线性方程组Ax=b的choleskly方法也叫做平方根法，这里对系数矩阵A是有要求的，需要A是对称正定矩阵，根据数值分析的相关理论，如果A对称正定，那么系数矩阵就可以被分解为的A=LLT形式，其中L是下三角矩阵，将其代入Ax=b中，可得：LLTx=b 进行如下分解：

yLTx Lyb那么就可先计算y,再计算x，由于L是下三角矩阵，是LT上三角矩阵，这样的计算比直接使用A计算简便，同时你应该也发现了工作量就转移到了矩阵的分解上面，

那么对于对称正定矩阵A进行Cholesky分解，我再描述一下过程吧： 如果你对原理很清楚那么这一段可以直接跳过的。

设A=LLT，即

a11a12a21a22an1an2...a1nl11l11l21...ln1l...a2nll...l212222n2 .....................annln1ln2...lnnlnn其中aijaji,i,j1,2,...,n

2,ai1li1l11故求得 第1步，由矩阵乘法，a11l11l11a11,li1ai1,i2,3,...n a11一般的，设矩阵L的前k-1列元素已经求出 第k步，由矩阵乘法得

akkll, 　aiklimlkmliklkk2km2kkm1m1k1k1

于是

k12lkkakklkmm1(k2,3,...,n) k1l1(all)　　,　ik1,k2,...n ikimkmiklm1kk2注意到akklkm，于是有

m1k

2lkmakkmaxaii

1in2这充分说明分解过程中元素lkm的平方不会超过系数矩阵A的最大对角元，

因而分解过程中舍入误差的放大收到了控制，用平方根法解对称正定方程组时可以不考虑选主元问题。

再说明一点，我的程序代码中采取了紧促格式存贮，就是将L和L转置的结果再放入矩阵A中存贮，由于Ｌ和Ｌ转置的主对角元元素是相同的这种存储格式可以节省存贮空间，在以后的很多代码中我都将采用这种存贮格式。

好了，终于讲完概念这一部分，下来再说下我的代码中要注意的一个问题，我用到了求解矩阵行列式求解的.cpp和.h文件来计算所给的行列式是否正定，这一部分代码可以参见我对于矩阵行列式求所提供的Determinant.h个Determinant.cpp文件，为了减小Word页数这里不再附上，而是直接调用其.h文件

好了下来直接上代码，代码我个人感觉还是编写的比较清晰，又都含有注释，如果你还有什么不明白的地方，或者我哪里有说不清楚的地方，你都可以给我留言哦。

Cholesky.h

#ifndef _CHOLRSKY_H #define _CHOLESKY_H

//---------------------------------------------------- // 使用平方法求解Ax=b

// 对于n阶 对称正定 矩阵A，称A=L(L转置)为矩阵A的Cholesky分解 // Cholesky分解不用进行选主元操作

//----------------------------------------------------

//---------------------------------------------------- // *a是A[n][n]的首地址&A[0][0]即A[0]

// 平方法计算线性非齐次方程组的结果存入数组x中 //----------------------------------------------------

bool Cholesky_Run(double *a, double *b, double *x, int n); #endif

Cholesky.cpp

#include #include #include \"Cholesky.h\" #include \"Determinant.h\"

//------------------------------------------------------

// 判断a[n][n]行列式的第m个顺序主子式是否大于0 // 返回值： true 大于0， false小于0 //------------------------------------------------------ bool CheckDeterminant(double *a, int m, int n) {

double result = CalDeterminant(a, m, n); if(result>=0) return true; else return false; }

//------------------------------------------------------ // 判断a[n][n]是否是对称正定矩阵

// 返回值： true 对称正定， false 不对称正定 //------------------------------------------------------ bool CheckArrayProperity(double *a, int n) {

int i, j;

for(i=0; i//检测是否对称 for(j=0; jfor(i=0; i//检测是否正定，即n个顺序主子式是否大于零 if(!CheckDeterminant(a, i+1, n)) { printf(\"系数矩阵不正定！\\n\"); return false;

} }

return true; }

//------------------------------------------------------ // 使用Cholesky分解法计算Ax=b // A是n*n矩阵

// 返回值 为真 可以进行Cholesky分解 // 为假 不可以进行Cholesky分解 //------------------------------------------------------

bool Cholesky_Run(double *a, double *b, double *x, int n) {

int i, j, k, m;

double ProcData;

//创建动态arr[n][n+1]作为(A|b) double **arr = new double *[n]; for(i=0; iarr[i] = new double [n+1]();

for(i=0; i//判断a[n][n]是否是对称正定矩阵 if(!CheckArrayProperity(a, n)) {

return false; } else {

//进行Cholesky分解 for(k=0; k { for(m=0; m//计算Ux =y,讲结果存入y(即存入arr[k][n]); for(k=n-1; k>=0; k--) { if(k==n-1) arr[k][n] /= arr[k][k]; else { ProcData = 0; for(m=k+1; m<=n-1; m++) ProcData += arr[k][m]*arr[m][n]; arr[k][n] -= ProcData; arr[k][n] /= arr[k][k]; } }

//将计算结果存入x数组 for(j=0; j//释放arr[n][n+1] for(i=0; ireturn true; } }

Main.cpp

#include #include \"Cholesky.h\" #include \"Determinant.h\"

void main() {

int i, j;

/*double A[3][3] = {

3, 2, 3, 2, 2, 0, 3, 0, 12 };

double b[3] = {5, 3, 7}; double x[3] = {0}; */

double A[5][5] = {

6, 1, 2, 4, 1, 1, 7, 2, 3, 5, 2, 2, 9, 6, 2, 4, 3, 6, 8, 5, 1, 5, 2, 5, 10 };

double b[5] = {5, 7, 3, 3, 4}; double x[5] = {0};

if(Cholesky_Run(A[0], b, x, 5))

}

{

double result[5] ; printf(\"\\n\");

printf(\"--------------------------\\n\"); for(i=0; i<5; i++) { result[i]=0; for(j=0; j<5; j++) { result[i] += A[i][j]*x[j]; } printf(\"x = %.7f result = %7.f\\n\ }

printf(\"\\n\"); } else {

printf(\"请修改您的系数矩阵为对称正定矩阵！\\n\"); }

——柴门野雪