Doolittle分解C++实现

• 1、三角分解法基本思想
• 2、Doolittle分解步骤
• 3、代码

当线性方程组AX=b时：

1. 先将非奇异矩阵A进行LU（这里指Doolittle分解）分解，使A=LU。则方程化为LUX=b。
2. 从而将问题转化为两个简单地三角方程组：

LY=b 求解Y
UX=Y 求解X

#include
#include
#include
using namespace std;

//Doolittle分解法
void Doolittle(double** a, double* b, int n)
{
	double** u = new double* [n];
	for (int i = 0; i < n; i++)
		u[i] = new double[n];
	double** l = new double* [n];
	for (int i = 0; i < n; i++)
		l[i] = new double[n];
	double* jie_x = new double[n];
	double* jie_y = new double[n];
	for (int k = 1; k < n + 1; k++)
	{
		for (int j = k, i = k + 1; j < n + 1, i < n + 1; j++, i++)
		{
			double sum_1 = 0;
			for (int t = 1; t < k; t++)
			{
				sum_1 += l[k - 1][t - 1] * u[t - 1][j - 1];
			}
			u[k - 1][j - 1] = a[k - 1][j - 1] - sum_1;
			double sum_2 = 0;
			for (int t = 1; t < k; t++)
			{
				sum_2 += l[i - 1][t - 1] * u[t - 1][k - 1];
			}
			l[i - 1][k - 1] = (a[i - 1][k - 1] - sum_2) / u[k - 1][k - 1];
		}
		double sum_1 = 0;
		for (int t = 1; t < k; t++)
		{
			sum_1 += l[k - 1][t - 1] * u[t - 1][n - 1];
		}
		u[k - 1][n - 1] = a[k - 1][n - 1] - sum_1;
	}
	jie_y[0] = b[0];
	for (int i = 2; i < n + 1; i++)
	{
		double sum = 0;
		for (int t = 1; t < i; t++)
			sum += l[i - 1][t - 1] * jie_y[t - 1];
		jie_y[i - 1] = b[i - 1] - sum;
	}
	jie_x[n - 1] = jie_y[n - 1] / u[n - 1][n - 1];
	for (int i = n - 1; i > 0; i--)
	{
		double sum = 0;
		for (int t = i + 1; t < n + 1; t++)
			sum += u[i - 1][t - 1] * jie_x[t - 1];
		jie_x[i - 1] = (jie_y[i - 1] - sum) / u[i - 1][i - 1];
	}
	cout << "使用Doolittle分解法的解为n";
	for (int i = 0; i < n - 1; i++)
		cout << "x" << i + 1 << " = " << jie_x[i] << endl;
	cout << "x" << n << " = " << jie_x[n - 1] << endl;
	cout << endl;

	delete[]jie_x;
	delete[]jie_y;
	for (int i = 0; i < n; i++)
		delete[]u[i];
	delete[]u;
	for (int i = 0; i < n; i++)
		delete[]l[i];
	delete[]l;
}


void main()
{

	//输入矩阵A和b
    int n;
	cout << "请输入系数矩阵的行数(或列数，列数等于行数):n";
	cin >> n;
	cin.clear();//这3行代码重置输入
	while (cin.get() != 'n')
		continue;
	cout << "请输入系数矩阵A:n";

	//要用到二维动态数组
	double** a;
	a = new double* [n];  //申请行的空间
	for (int i = 0; i < n; i++)
		a[i] = new double[n];
	for (int i = 0; i < n; i++)
	{
		cout << "第 " << i + 1 << " 行:";
		for (int j = 0; j < n; j++)
			while (!(cin >> a[i][j]))
			{
				cin.clear();
				while (cin.get() != 'n')
					continue;
				cout << "请输入数字";
			}
		cin.clear();
		while (cin.get() != 'n')
			continue;
	}
	cout << "请输入矩阵b:n";
	double* b = new  double[n];
	for (int j = 0; j < n; j++)
		cin >> b[j];
	cin.clear();
	while (cin.get() != 'n')
		continue;

	//将系数矩阵A与矩阵b存储在临时矩阵中，临时矩阵中的元素参与行列变换
	double** temp_a = new double* [n];
	double* temp_b = new double[n];
	for (int i = 0; i < n; i++)
		temp_a[i] = new double[n];
	for (int i = 0; i < n; i++)
		for (int j = 0; j < n; j++)
			temp_a[i][j] = a[i][j];
	for (int i = 0; i < n; i++)
		temp_b[i] = b[i];

	Doolittle((double**)temp_a, temp_b, n);

	for (int i = 0; i < n; i++)
		delete[]temp_a[i];
	delete[]temp_a;
	delete[]temp_b;
	for (int i = 0; i < n; i++)
		delete[]a[i];
	delete[]a;
	delete[]b;
}

待博主有空后，再将上述基础知识补齐。