发布时间 : 星期四 文章平方根法计求解线性方程组更新完毕开始阅读2bbfcf1a31126edb6f1a10da
解线性n阶方程组直接法
—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
//------------------------------------------------------
// 判断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; j for(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; i arr[i] = new double [n+1](); for(i=0; i //判断a[n][n]是否是对称正定矩阵 if(!CheckArrayProperity(a, n)) { return false; } else {