计算方法——C语言实现——LU分解法求解非线性方程

最近在上计算方法这门课，要求是用MATLAB做练习题，但是我觉得C语言也很棒棒啊~

题目：

高斯消元法实际上可以看成是将系数矩阵A分解为一个单位上三角矩阵L和一个下三角矩阵U的乘积，当高斯消元无法使用的时候当然LU分解也不能用，只要A各阶顺序主子式不为0就行了。

使用VS2017，代码如下：

//使用LU分解法求解线性方程组
#include "stdafx.h"
#include<stdlib.h>
#include "math.h"

double **A, *b, *x, *y, **L, **U;
unsigned int RANK = 4;
unsigned int makematrix()
{
	unsigned int r, c;

	printf("请输入矩阵行列数，用空格隔开：");
	scanf_s("%d %d", &r, &c);

	A = (double**)malloc(sizeof(double*)*r);//创建一个指针数组，把指针数组的地址赋值给a ,*r是乘以r的意思
	for (int i = 0; i < r; i++)
		A[i] = (double*)malloc(sizeof(double)*c);//给第二维分配空间
	for (int i = 0; i < r; i++) {
		for (int j = 0; j < c; j++)
			A[i][j] = 0.0;
	}

	b = (double*)malloc(sizeof(double)*r);
	for (int i = 0; i < r; i++)
	{
		b[i] = 0.0;
	}
	x = (double*)malloc(sizeof(double)*c);
	for (int i = 0; i < c; i++)
	{
		x[i] = 0.0;
	}

	L = (double**)malloc(sizeof(double*)*r);//创建一个指针数组，把指针数组的地址赋值给a ,*r是乘以r的意思
	for (int i = 0; i < r; i++)
		L[i] = (double*)malloc(sizeof(double)*c);//给第二维分配空间
	for (int i = 0; i < r; i++) {
		for (int j = 0; j < c; j++)
			L[i][j] = 0.0;
	}
	U = (double**)malloc(sizeof(double*)*r);//创建一个指针数组，把指针数组的地址赋值给a ,*r是乘以r的意思
	for (int i = 0; i < r; i++)
		U[i] = (double*)malloc(sizeof(double)*c);//给第二维分配空间
	for (int i = 0; i < r; i++) {
		for (int j = 0; j < c; j++)
			U[i][j] = 0.0;
	}
	y = (double*)malloc(sizeof(double)*c);
	for (int i = 0; i < c; i++)
	{
		y[i] = 0.0;
	}
	return r;
}

void getmatrix(void)//输入矩阵并呈现
{
	printf("请按行从左到右依次输入系数矩阵A，不同元素用空格隔开\n");
	for (int i = 0; i < RANK; i++)
	{
		for (int j = 0; j<RANK; j++)
		{
			scanf_s("%lf", &A[i][j]);
		}
	}
	printf("系数矩阵如下\n");
	for (int i = 0; i < RANK; i++)
	{
		for (int j = 0; j<RANK; j++)
		{
			printf("%g\t", A[i][j]);
		}
		printf("\n");
	}
	printf("请按从上到下依次输入常数列b，不同元素用空格隔开\n");
	for (int i = 0; i<RANK; i++)
	{
		scanf_s("%lf", &b[i]);
	}
	printf("常数列如下\n");
	for (int i = 0; i<RANK; i++)
	{
		printf("%g\t", b[i]);
	}printf("\n");
}

void LUrush_calculation(void)//LU法解线性方程组
{
	double get_add=0.0;
	printf("利用以上A与b组成的增广阵进行LU分解法计算方程组\n");

	for (int i = 0; i < RANK; i++)
	{
		U[0][i] = A[0][i];
		L[i][i] = 1;
	}
	for (int i = 1; i < RANK; i++)
	{
		L[i][0] = A[i][0]/U[0][0];
	}
	for (int i = 1; i < RANK; i++)
	{
		for (int j = 0; j<RANK; j++)
		{
			if (i <= j)
			{
				get_add = 0.0;
				for (int k = 0; k <= i - 1; k++) get_add = get_add + L[i][k] * U[k][j];
				U[i][j] = A[i][j] - get_add;
			}
			if (i < j)
			{
				get_add = 0.0;
				for (int k = 0; k <= i - 1; k++) get_add = get_add + L[j][k] * U[k][i];
				L[j][i] = (A[j][i] - get_add) / U[i][i];
			}
		}
	}
	printf("由系数矩阵A得到的L U矩阵如下\n");
	for (int i = 0; i < RANK; i++)
	{
		for (int j = 0; j<RANK; j++)
		{
			printf("%g\t", L[i][j]);
		}		printf("\t\t");
		for (int j = 0; j<RANK; j++)
		{
			printf("%g\t", U[i][j]);
		}		printf("\n");
	}
	printf("利用 Ly = d求解y，解得：\n");
	y[0] = b[0];
	for (int i = 1; i<RANK; i++)
	{
		get_add = 0.0;
		for (int k = 0; k <= i - 1; k++) get_add = get_add + L[i][k] *y[k];
		y[i] = b[i] - get_add;
	}
	for (int i = 0; i<RANK; i++)
	{
		printf("y%d = %g\n", i + 1, y[i]);
	}
	printf("利用 Ux = y求解x，解得：\n");
	x[RANK - 1] = y[RANK - 1] / U[RANK - 1][RANK - 1];
	for (int i = RANK - 2; i >= 0; i--)
	{
		get_add = 0.0;
		for (int k = i+1; k <= RANK - 1; k++) get_add = get_add + U[i][k] * x[k];
		x[i] = (y[i] - get_add) / U[i][i];
	}
	for (int i = 0; i<RANK; i++)
	{
		printf("x%d = %g\n", i + 1, x[i]);
	}
	printf("计算完成，按回车退出程序或按1重新输入矩阵\n");
}

int main()
{
_again:
	RANK = makematrix();
	getmatrix();
	LUrush_calculation();
	getchar();

	if ('1' == getchar())
		goto _again;
	return 0;
}

           

按设计的提示老老实实 输入题目的系数矩阵和常数向量后，得到运行结果：

可以看出来，X4其实就是0，但是显示不是0，这并不是程序有问题，而是这就是直接法求线性方程组会失真的体现。