平方根法求解线性方程组- -函数模块化编程极易调用（C++实现）

• 算法思想

由于，A矩阵可以表示为,

此时，线性方程组Ax=b变为下面两个三角形方程组:.

(1)Gy=b，求y；
(2)G^Tx=y，求x。                                                

• 算法理论：

如果A为n阶对称正定矩阵，则存在唯一的主对角线元素都是正数的下三角阵G，使得A=GG^T

乔列斯基分解矩阵G计算公式：

 

• 算法实现流程图

 

• 代码：

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

const int n = 3;
void solve_G(const double a[][n], double G[][n], const int n);//使用乔列斯基分解求解矩阵G
void zhuanzhi(const double G[][n], double G2[][n], const int n);//矩阵的转置
void fillarr2(double a[][n], const int n);//填充二维数组系数矩阵元素
void fillarr1(double* b, const int n);//填充一维数组b矩阵
void showarr2(const double a[][n], const int n);//显示二维数组元素内容
void showarr1(const double* arr, const int n);//显示一维数组元素内容
void solve_y(const double G[][n], const int n, double* y, double* b);//求解中间矩阵y，其中，G是下三角方阵
void solve_x(const double G2[][n], const int n, const double* y, double* x);//求解矩阵x，其中，矩阵G2（矩阵G的转置）是上三角方阵

int main()
{

	double a[n][n];
	/* 测试数据 ：double a[n][n] = {{4,2,-2},{2,2,-3},{-2,-3,14}};
					double b[n] = { 10,5,4 };
		结果：			x=[2,2,1];	*/
	
	double b[n];// = { 10,5,4 };		//

	//初始化系数矩阵A 和参数矩阵b
	cout << "请输入系数矩阵a： " << endl;
	fillarr2(a, n);
	cout << "系数矩阵a：" << endl;
	showarr2(a, n);
	cout << "请输入矩阵b： " << endl;
	fillarr1(b, n);
	cout << "参数矩阵b：" << endl;
	showarr1(b, n);
	double G[n][n] = {};
	//求矩阵G
	solve_G(a, G, n);
	cout << "矩阵G= :" << endl;
	showarr2(G, n);

	//求矩阵G的转置
	cout << "矩阵G的转置G'= :" << endl;
	double G2[n][n];
	zhuanzhi(G, G2, n);
	showarr2(G2, n);
	//求解y
	double y[n] = { };
	solve_y(G, n, y, b);
	cout << "矩阵y的值为：" << endl;
	showarr1(y, n);
	//求解x
	double x[n]={ };
	solve_x(G2, n, y, x);
	cout << "矩阵x的值为：" << endl;
	showarr1(x, n);




	return 0;
}


void solve_G(const double a[][n], double G[][n], const int n)
{
	int i, j, r;
	double sum = 0;
	G[0][0] = sqrt(a[0][0]);//第一行第一列元素a[0][0]等于自身开根号
	for (i = 1; i < n; i++)
	{
		G[i][0] = a[i][0] / G[0][0];//从第二行开始第一列元素等于自身除以a[0][0]
	}
	for (i = 1; i < n; i++)
	{

		for (j = 1; j <= i; j++)
		{
			if (i == j)
			{
				sum = 0;
				for (r = 0; r <= j - 1; r++)
					sum = sum + G[i][r] * G[i][r];
				G[i][j] = sqrt(a[i][j] - sum);
			}
			else if (i != j)
			{
				sum = 0;
				for (r = 0; r <= j - 1; r++)
					sum = sum + G[i][r] * G[j][r];
				G[i][j] = (a[i][j] - sum) / G[j][j];
			}
		}
	}
}

void zhuanzhi(const double G[][n], double G2[][n], const int n)
{
	int i, j;
	for (i = 0; i < n; i++)
	{
		for (j = 0; j < n; j++)
			G2[i][j] = G[j][i];
	}

}
//填充二维数组系数矩阵元素
void fillarr2(double a[][n], const int n)
{
	int i, j;
	cout << "请输入系数矩阵A：" << endl;
	for (i = 0; i < n; i++)
	{
		cout << "请输入第" << i + 1 << "行系数:";
		for (j = 0; j < n; j++)
			cin >> a[i][j];
	}
}
//填充一位数组b矩阵
void fillarr1(double* b, const int n)
{
	int i;
	for (i = 0; i < n; i++)
		cin >> b[i];
}
//显示二维数组元素内容
void showarr2(const double a[][n], const int n)
{

	int i, j;
	for (i = 0; i < n; i++)
	{
		for (j = 0; j < n; j++)
			cout << a[i][j] << "  ";
		cout << endl;
	}

}
//显示一维数组元素内容
void showarr1(const double* arr, const int n)
{
	int i;
	for (i = 0; i < n; i++)
	{
		cout << arr[i] << "  ";
	}
	cout << endl;
}

• 实例验证：

代码运行界面：