平方根法求解线性方程组（C++实现）

背景：应用有限元法求解结构力学问题时，最后归结为求解线性方程组，该系数矩阵大多具有对称正定性质。平方根法就是利用对称正定的三角分解求解对称正定方程组的一种有效方法。
前提：系数矩阵A为对称矩阵，且A的所有顺序主子式均不为0。

算法1：将系数矩阵A分解为L和L的转置。则求解Ax=b的过程，就变为求解Ly=b，（L转置）x=y的过程。 但是该算法会出现开平方，在计算机中，开平方这个算法本身就对浮点数有误差，传递起来会很大，并且开平方里面的数为正数。限制较多，且误差大。本文第一个实现了此算法，但是求解过程不尽如人意。所以推荐改进的平方根法。

算法2：改进的平方根算法，将A分解为LD（L的转置）。同上然后求解就行。此算法比较稳定，且误差较小。
对于上述算法，由于都是对称矩阵，所以存储的时候可以存储在一维数组中，节省空间资源。对于A的n维矩阵，存储在一维数组的元素个数为n（n-1）/2。对于aij对应于一维数组的a[I(I+1)/2+j]。
本算法并没有用一维数组，由于时间关系等因素，等空闲在实现，思路就是将输入数组a[N][N]修改为a[N(N-1)/2]，关系就是刚才的。直接替换应该就可以。
下面贴代码：

#include<iostream>
#include<cmath>
using namespace std;

const int n = 3;


//平方根法求解线性方程组
//前提：系数矩阵必须是对称矩阵，且A的顺序主子式都不为0。
//可以分解为A=L和L的转置
//现在按这种方法来求解
//由于是对称矩阵，为了节省空间，就用一维数组进行存储

template<class T>
void Square(T a[n][n], T b[n])
{
	T l[n][n];             //存储矩阵
	T y[n];                //中间求解向量
	T x[n];                    //求解结果
	int i, j, k;

	//求解l[N][N]的过程
	for (j = 0; j < n; j++)
	{
		T sum1 = 0;
		for (k = 0; k < j; j++)
		{
			sum1 += l[j][k] * l[j][k];
		}
		l[j][j] = sqrt(a[j][j] - sum1);             //开平方根，会损失精度
		for (i = j + 1; j < n; j++)
		{
			T sum2 = 0;
			for (k = 0; k < j; k++)
			{
				sum2 += l[i][k] * l[j][k];
			}
			if (l[j][j] < 0)                                    //判断除数是否为0，为0抛出异常
			throw("division is zero,can't do!!");
			l[i][j] = (a[i][j] - sum2)/l[j][j];
		}
	}
	//求解y
	for (i = 0; i < n; i++)
	{
		T sum3 = 0;
		for ( k = 0; k < i; k++)
			sum3 += l[i][k] * y[k];
		y[i] = (b[i] - sum3) / l[i][i];
	}
	//输出y的值
	for (i = 0; i < n; i++)
		cout << y[i] << endl;
	
	
	//求解x
	for (i = n - 1; i >= 0; --i)
	{
		T sum4 = 0;
		for (k = i + 1; k < n; k++)
			sum4 += l[k][i] * x[k];
		x[i] = (y[i] - sum4) / l[i][i];		//原先是x[i] = (b[i] - sum4) / l[i][i];，经评论区提醒，已修改。
	}

	//输出x的值
	for (i = 0; i < n; i++)
		cout << "x[" << i + 1 << "]=" << x[i] << "	" << endl;
}

//改进的开平方根算法
template<class T=double>
void Improve_Square(T a[n][n], T b[n])
{
	int i, j, k;
	T t[n][n];
	T d[n];
	T l[n][n];
	//求解过程
	for (i = 0;i < n; i++)
	{  
		T sum1 = 0;
		for (j = 0; j <= i - 1; j++)
		{
			for (k = 0; k <= j - 1; k++)
				sum1 += t[i][k] * l[j][k];
			t[i][j] = a[i][j] - sum1;
			l[i][j] = t[i][j] / d[j];
		}
		T sum2 = 0;
		for (k = 0; k <= i - 1; k++)
			sum2 += t[i][k] * l[i][k];
		d[i] = a[i][i] - sum2;
	}

	//打印输出中间矩阵D
	for (i = 0; i < n; i++)
		cout << " d[n]" << d[i] << "    ";
	cout << endl;
	//Ly=b；
	T y[n];
	y[0] = b[0];
	for (i = 1; i < n; i++)
	{
		T sum3 = 0;
		for (k = 0; k <= i - 1; k++)
			sum3 += l[i][k] * y[k];
		y[i] = b[i] - sum3;
	}
	//打印输出中间矩阵y
	for (i = 0; i < n; i++)
		cout << " y[n]" << y[i] << "	";
	cout << endl;
	//求解D(L转置)x=y
	T x[n];
	x[n-1] = y[n - 1] / d[n - 1];
	for (i = n - 2; i >= 0; i--)
	{
		T sum4 = 0;
		for (k = i + 1; k < n; k++)
			sum4 += l[k][i] * x[k];
		x[i] = y[i] / d[i] - sum4;
	}

	//打印输出
	for(i=0;i<n;i++)
		cout << "x[" << i + 1 << "]=" << x[i] << "	" << endl;
}


int main()
{
	double a[3][3] = {2,-1,1,-1,-2,3,1,3,1};
	double b[3] = {4,5,6};
	Improve_Square(a,b);
	Square(a, b);
	return 0;
}

由于这个算法用的不多，也就少跑了点数据，如果有不正确的，请指出~