C++学习日记#7——方程组的直接分解法求解

直接分解得出的LU矩阵与追赶法分解得出的矩阵有较大的区别，直接分解法得出的L矩阵的主对角线元素都为1，而且下三角形元素(除了特殊情况)都非0，U矩阵上三角形没有固定的规律。当分解矩阵时候只能逐一通过LU交替求解。毕竟直接法不像追赶法没有有规律性。

分解公式：

LU=A

程序思路

虽然直接法需要逐一求解每一个元素，但是由于C++中含有+=的累加计算符号，使得通过循环语句就可以逐一求出LU矩阵的每个元素成为了可能。不过在求解U矩阵的第一行元素时，与求解其他行元素有区别，主要是第一行元素求解时是没有之前元素的累加，也就是说，可以通过A矩阵的第一行元素依次除以L矩阵的第一行元素，也就是除以1,换句话来说A矩阵的第一行与U矩阵的第一行元素是相同的。为了将第一行求解公式与后面行的元素求解公式归纳在一起，程序将赋值U矩阵一第0行，以及L矩阵的第0列。这些元素均为0。

	for (int i = 1; i <= n; i++)
	{
		L[i][i] = 1;
		L[i][0] = 0;
		U[0][i] = 0;
	}//定义LU矩阵的初始值：L对角线为1；L矩阵添加第0列；U矩阵添加第0行；

附上实现程序（Visual Studio 2019）

#include <iostream>
#include <iomanip>//保留小数

using namespace std;
double A[100][100];//定义系数矩阵A
double L[100][100];//定义分解矩阵L
double U[100][100];//定义分解矩阵U
double Y[100]; //定义矩阵Y
double b[100]; //定义矩阵b
double x[100]; //定义求解结果x


int main()
{
	int n;
	cout << "请输入A矩阵的阶数" << endl;
	cin >> n;
	cout << "***********请输入A矩阵***********" << endl;
	for (int i = 1; i <= n; i++)
	{
		cout << "请输入第" << i << "行的元素" << endl;
		for (int j = 1; j <= n; j++)
		{
			cin >> A[i][j];
		}
	}
	cout << "***********请输入b矩阵***********" << endl;
	for (int i = 1; i <= n; i++)
	{
		cin >> b[i];
	}//输入模块
	for (int i = 1; i <= n; i++)
	{
		L[i][i] = 1;
		L[i][0] = 0;
		U[0][i] = 0;
	}//定义LU矩阵的初始值：L对角线为1；L矩阵添加第0列；U矩阵添加第0行；主要是因为L和U矩阵分别按照列元素和行元素递推计算；
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			L[i][j + i] = 0;
			U[j + i][i] = 0;
		}
	}
	//矩阵上三角元素为0以及主对角线为1、赋值了U矩阵下三角元素为0；
	for (int i = 1; i <= n; i++)
	{
		for (int k = i; k <= n; k++)
		{
			double S_u = 0;//重置求和S_u，用于下一次循环
			for (int p = 0; p <= i - 1; p++)
			{
				S_u += L[i][p] * U[p][k];//求和S_u
			}
			U[i][k] = (A[i][k] - S_u) / L[i][i];//求解U行元素
			for (int j = i + 1; j <= n; j++)
			{
				double S_l = 0;//重置求和S_l，用于下一次循环
				for (int p = 0; p <= i - 1; p++)
				{
					S_l += L[j][p] * U[p][i];///求和S_l
				}
				L[j][i] = (A[j][i] - S_l) / U[i][i];//求解L列元素
			}
		}
	}//L和U矩阵分别按照列元素和行元素递推计算；
	for (int i = 1; i <= n; i++)
	{
		double S_y = 0; //重置求和S_y，用于下一次循环
		Y[0] = 0;
		for (int p = 0; p <= i - 1; p++)
		{
			S_y += L[i][p] * Y[p];///求和S_y
		}
		Y[i] = (b[i] - S_y) / L[i][i];
	}
	for (int i = 1; i <= n; i++)
	{
		double S_x = 0;//重置求和S_x，用于下一次循环
		x[n + 1] = 0;
		U[n][n + 1] = 0;
		for (int p = 0; p <= i - 1; p++)
		{
			S_x += U[n - i + 1][n + 1 - p] * x[n + 1 - p];///求和S_y
		}
		x[n - i + 1] = (Y[n - i + 1] - S_x) / U[n - i + 1][n - i + 1];
	}
	cout << "*************输出结果************" << endl;
	cout << setiosflags(ios::fixed) << setprecision(5);//结果保留5位小数
	cout << "输出L矩阵" << endl;
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			cout << L[i][j] << " ";
		}
		cout << endl;
	}
	cout << "输出U矩阵" << endl;
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			cout << U[i][j] << " ";
		}
		cout << endl;
	}
	cout << "输出方程组的解x" << endl;
	for (int i = 1; i <= n; i++)
	{
		cout << "x" << i << ": " << x[i] << endl;
	}
	return 0;
}

例题

与追赶法的例题对比

https://blog.csdn.net/ALEX_2814/article/details/121050311https://blog.csdn.net/ALEX_2814/article/details/121050311https://blog.csdn.net/ALEX_2814/article/details/121050311

 

 后话

其实#6和#7这两篇日记在周一（11.8)就已经把程序弄好了，直到周四(11.11)才开始总结成日记，导致没有很好介绍程序里面的一些主要求解函数。而且将直接法编成程序的过程还是比较需要功夫的，因为LU矩阵的互为前置条件，导致了求解LU时必须将函数汇集在同一组循环中才能得出结果。当理解了这一个注释后，对于求解LU矩阵就比较简单了。

//定义LU矩阵的初始值：L对角线为1；L矩阵添加第0列；U矩阵添加第0行；主要是因为L和U矩阵分别按照列元素和行元素递推计算；