C/C++ 语言实现使用LU分解求解线性方程组

应用计算方法C语言程序：02

接应用计算方法C语言程序：01：C/C++ 语言实现矩阵LU分解：Doolittle（A）

以计算方法课本例题4-6为例。

设矩阵A,b分别为：A[3][3] = { {1,3,3},{2,1,1},{2,3,4} }; b[3] = { 1,2,1 };
利用C语言程序求解方程组如下：（其中LU分解程序见应用计算方法C语言程序：01）
程序中矩阵求逆部分程序使用了博友（重中之重小星星）的程序（感谢博主的分享）：原文链接

void Matrix_inverse(double arc[3][3], int n, double ans[3][3])//计算矩阵的逆
{
	int i, j, k;//列
	double max, tempA, tempB, P;
	int max_num;
	double arcs[3][3];
	memcpy(arcs, arc, n*n*sizeof(double));
	for (i = 0; i < n; i++)
	{
		ans[i][i] = 1;
	}
	for (i = 0; i < n; i++)//第i列
	{
		max = fabs(arcs[i][i]);
		max_num = i;
		for (j = i + 1; j < n; j++)//选出主元
		{
			if (fabs(arcs[j][i]) > max)
			{
				max = fabs(arcs[j][i]);
				max_num = j;
			}
		}
		for (k = 0; k < n; k++)//交换行
		{
			tempA = arcs[i][k];
			arcs[i][k] = arcs[max_num][k];
			arcs[max_num][k] = tempA;
			tempB = ans[i][k];
			ans[i][k] = ans[max_num][k];
			ans[max_num][k] = tempB;
		}
		for (k = i + 1; k < n; k++)
		{
			P = arcs[k][i] / arcs[i][i];
			for (j = 0; j < n; j++)
			{
				arcs[k][j] = arcs[k][j] - arcs[i][j] * P;
				ans[k][j] = ans[k][j] - ans[i][j] * P;
			}
		}
	}
	for (i = 0; i < n; i++)//行
	{
		P = arcs[i][i];
		for (j = i; j < n; j++)
		{
			arcs[i][j] = arcs[i][j] / P;
		}
		for (j = 0; j < n; j++)
		{
			ans[i][j] = ans[i][j] / P;
		}
	}
	for (i = n - 1; i > 0; i--)
	{
		for (j = i - 1; j >= 0; j--)
		{
			for (k = 0; k < n; k++)
			{
				ans[j][k] = ans[j][k] - ans[i][k] * arcs[j][i];
			}
		}
	}
}
int main()
{
	double A[3][3] = { {1,3,3},{2,1,1},{2,3,4} };
	double b[3] = { 1,2,1 };

	Doolittle(A);

	double b1[3][3] = { 0 };
	double b2[3][3] = { 0 };
	Matrix_inverse(U, 3, b1);//求矩阵U的逆
	Matrix_inverse(L, 3, b2);//求矩阵L的逆
	
	
	//x=U^(-1) * L^(-1) * b
	double b3[3][3] = { 0 };
	//U^(-1) * L^(-1)
	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			for (int k = 0; k < 3; k++)
			{
				b3[i][j] += b1[i][k] * b2[k][j];
			}
		}
	}
	//x=U^(-1) * L^(-1) * b
	double b4[3] = { 0 };
	for (int i = 0; i < 3; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			b4[i]+= b3[i][j] * b[j];
		}
	}

	//输出方程的解
	cout << endl;
	for (int i = 0; i < 3; i++)
	{
		cout << "x" << i + 1 << " = " << b4[i] << " ";
	}

	cout << endl << endl << endl;

	return 0;
}

运行结果如下：