（二）计算机数值方法之Cholesky分解法

数学问题：

利用 Cholesky分解法解决线性方程问题Ax=b,其中(A,b)分别为:

解决代码：

​
#include "windows.h"
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;
void SolverEqCholesky(double** A, double* b, int n, double* x);
int main()
{
	
	int i, n;
	double** A, * b, * x;
	//输入维数
	cin >> n;
	A = new double* [n];
	b = new double[n];
	x = new double[n];
	for (i = 0; i < n; i++)
	{
		A[i] = new double[n];
	}
	//输入系数矩阵和向量
	for (i = 0; i < n; i++)
	{
		for (int j = 0; j < n; j++)
		{
			cin >> A[i][j];
		}
	}
	for (int i = 0; i < n; i++)
	{
		cin >> b[i];
	}
	//计算
	SolverEqCholesky(A, b, n, x);
	//输出
	for (i = 0; i < n; i++)
	{
		cout << fixed << setprecision(6) << x[i] << " ";
	}
	delete[]b;//释放内存
	delete[]x;
	for (i = 0; i < n; i++)
	{
		delete[]A[i];
	}
	
	return 0;
}
double** two_array_malloc(int m, int n)   //二维数组分配
{
	double** a;
	int i, j;
	a = (double**)malloc(m * sizeof(double*));
	for (i = 0; i < m; i++)
	{
		a[i] = (double*)malloc(n * sizeof(double));
		for (j = 0; j < n; j++)
		{
			a[i][j] = 0;
		}
	}
	return a;
}
void two_array_clear(double** a, int m)   //二维数组释放
{
	int i;
	for (i = 0; i < m; i++)
	{
		delete[] a[i];
	}
	delete[] a;
}
bool blabel = true;
void SolverEqCholesky(double** A, double* b, int n, double* x)
{
	int  i, j, k;
	double** a;
	double temp;
	a = two_array_malloc(n, n + 1);
	for (i = 0; i < n; i++)
	{
		for (j = 0; j < n; j++)
		{
			a[i][j] = A[i][j];
		}
		a[i][n] = b[i];
	}
	double** L;
	L = two_array_malloc(n, n);
	//分解过程
	//==============begin1===============
	for (j = 0; j < n; j++)
	{
		double qiuhe = 0;
		for (k = 0; k < j ; k++)
		{
			qiuhe = qiuhe + L[j][k] * L[j][k];
		}
		L[j][j] = sqrt(abs(a[j][j] - qiuhe));
		
		for (i = j+1; i < n; i++)
		{
			double qiuh = 0;
			for (k = 0; k <j; k++)
			{
				qiuh=qiuh+ L[i][k] * L[j][k];
			}
			L[i][j] = (a[i][j] - qiuh) / L[j][j];
		}
	}
	//==============end1=================
	double** tem;
	tem = two_array_malloc(2, n);
	//回代过程
	//==============begin2===============
	for (i = 0; i < n; i++)
	{
		double qiu = 0;
		for (int j = 0; j < i; j++)
		{
			qiu = qiu + L[i][j] * tem[0][j];
		}
		tem[0][i] = (a[i][n] - qiu) / L[i][i];
	}
	for (i = n - 1; i >= 0; i--)
	{
		double qi = 0;
		for (int j = i + 1; j < n; j++)
		{
			qi = qi + L[j][i] * tem[1][j];
		}
		tem[1][i] = (tem[0][i] - qi) / L[i][i];
	}
	//==============end2=================
	//输出解(在编辑区填写代码)
	//==============begin3===============
	for (i = 0; i < n; i++)
	{
		x[i] = tem[1][i];
	}
	two_array_clear(a, n);
	two_array_clear(L, n);
	two_array_clear(tem, 2);
	//==============end3=================
}

​

使用方法：

第一行是矩阵维度：n ；

后面的n行，输入n维矩阵A；

最后一行输入未知向量的值b

例如：

矩阵为：A=[[10,7,8,7],[7,5,6,5],[8,6,10,9],[7,5,9,10]]

向量为：b=[32,23,33,31]

测试输入：

4

10 7 8 7

7 5 6 5

8 6 10 9

7 5 9 10

32 23 33 31

预期输出： 1.000000 1.000000 1.000000 1.000000

问题解决：

注意：

Cholesky分解法的使用对象与（一）中的Gauss列主元素消去法略有不同。它主要用于对称正定的线性方程组，其矩阵必须是对称且正定的。该方法将对称正定矩阵分解成下三角矩阵的转置和该下三角矩阵之间的乘积。其计算效率较高，特别在处理大型对称正定矩阵时效果更为明显。由于其对称正定矩阵的性质，Cholesky分解通常更为稳定，能够减少数值误差的传播。