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

解线性n 阶方程组直接法
—Cholesky 方法
解n 阶线性方程组Ax=b 的choleskly 方法也叫做平方根法，这里对系数矩阵A 是有要求的，需要A 是对称正定矩阵，根据数值分析的相关理论，如果A 对称正定，那么系数矩阵就可以被分解为的T A=L L ∙形式，其中L 是下三角矩阵，将其代入Ax=b 中，可得：T LL x=b
进行如下分解：
T L x L b y y ⎧=⎨=⎩
那么就可先计算y,再计算x ，由于L 是下三角矩阵，是T L 上三角矩阵，这样的计算比直接使用A 计算简便，同时你应该也发现了工作量就转移到了矩阵的分解上面，
那么对于对称正定矩阵A 进行Cholesky 分解，我再描述一下过程吧： 如果你对原理很清楚那么这一段可以直接跳过的。

设T A=L L ∙，即
1112111112112122
221222221212....................................n n n n n n nn n n nn nn a a a l l l l a a a l l l l a a a l l l l ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦
其中,,1,2,...,ij ji a a i j n ==
第1步，由矩阵乘法，21111
1111,i i a l a l l == 故求得
111111
,2,3,...i i a l l i n a === 一般的，设矩阵L 的前k-1列元素已经求出
第k 步，由矩阵乘法得
112211k k kk km kk
ik im km ik kk m m a l l a l l l l --===+=+∑∑,
于是
11(2,3,...,n)1(),1,2,...kk k ik ik im km m kk l k l a l l i k k n l -=⎧=⎪⎪=⎨⎪=-=++⎪⎩
∑ 注意到21k
kk km m a l ==∑，于是有
21max km kk ii i n
l a a ≤≤≤≤ 这充分说明分解过程中元素2km l 的平方不会超过系数矩阵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<stdio.h>
#include<math.h>
#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<n; i++)
{
//检测是否对称
for(j=0; j<n; j++)
{
double x1 = *(a + i*n + j);
double x2 = *(a + j*n + i);
if(x1!=x2)
{
printf("系数矩阵不对称！\n");
return false;
}
}
}
for(i=0; i<n;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<n; i++)
arr[i] = new double [n+1]();
for(i=0; i<n;i++)
for(j=0; j<n+1; j++)
{
if(j==n) arr[i][j] = b[i];
else arr[i][j] = *(a + n*i + j);
}
//判断a[n][n]是否是对称正定矩阵
if(!CheckArrayProperity(a, n))
{
return false;
}
else
{
//进行Cholesky分解
for(k=0; k<n; k++)
{
for(j=k; j<n; j++)
{
ProcData =0;
if(j==k)
{
if(k!=0)
{
for(m=0; m<k; m++)
ProcData += arr[k][m]*arr[m][k];
}
arr[k][j] = sqrt((double)(arr[k][k] - ProcData));
}
else
{
if(k!=0)
{
for(m=0; m<k; m++)
ProcData += arr[k][m]*arr[m][j];
}
arr[k][j] -= ProcData;
arr[k][j] /= arr[k][k];
}
}
for(i=k; i<n; i++)
arr[i][k] = arr[k][i];
//计算Ly = b;将y的结果存入b(即存入arr[k][n]);
if(k==0) arr[k][n] /= arr[k][k];
else
{
ProcData = 0;
for(m=0; m<k; m++)
ProcData += arr[k][m]*arr[m][n];
arr[k][n] -= ProcData;
arr[k][n] /= arr[k][k];
}
}
//计算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<n; j++)
x[j] = arr[j][n];
//释放arr[n][n+1]
for(i=0; i<n; i++)
delete []arr[i];
delete []arr;
return true;
}
}
Main.cpp
#include <stdio.h>
#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", x[i], result[i]);
}
printf("\n");
}
else
{
printf("请修改您的系数矩阵为对称正定矩阵！\n");
}
}
——柴门野雪。