用LU分解法解非齐次线性方程组Ax=b

LU.h

#include <iostream>

double **mat_L, **mat_U; //矩阵L U
int N; //方阵的阶
class LU {
public :
	LU();

public :
	void or_mat(int N); //初始矩阵(系数矩阵mat_A, 右值矩阵mat_b)初始化,N为方阵阶
	void creatLU(void); //将矩阵A 分解为L U
	void cal_mat(void); //用LU分解法求Ax=b的解

private :
	double **mat_A, *mat_b;
	double *result; //结果数组
};

 LU_File.h

#include <iostream>
#include <malloc.h>

#include "LU.h"
////////////////////////////
LU::LU() {
	mat_A = NULL;
	mat_b = NULL;
}

/////////////////////////
void LU::or_mat(int n) {
	N = n;
	std::cout << "方阵的维数为 " << N << '\n';
	std::cout << "input the mat_A" << '\n';
	mat_A = (double **)malloc(N*sizeof(double));
	for(int i=0; i<N; i++)
		*(mat_A+i) = (double *)malloc(N*sizeof(double));
	for(i=0; i<N; i++)
		for(int j=0; j<N; j++)
			std::cin >> *(*(mat_A+i)+j);
	std::cout << "input the mat_b" << '\n';
	mat_b = (double *)malloc(N*sizeof(double));
	for(i=0; i<N; i++)
		std::cin >> *(mat_b+i);
}

////////////////////////////////////////
void LU::creatLU() {
	mat_L = (double **)malloc(N*sizeof(double));
	for(int i=0; i<N; i++)
		*(mat_L+i) = (double *)malloc(N*sizeof(double));
	mat_U = (double **)malloc(N*sizeof(double));
	for(i=0; i<N; i++)
		*(mat_U+i) = (double *)malloc(N*sizeof(double));

	for(int k=0; k<N; k++) { //初始化L U的已知元素
		*(*(mat_U+0)+k) = *(*(mat_A+0)+k); //U的0行=A的0行
		*(*(mat_L+k)+k) = 1.0; //L的对角线元素为0
		if(k!=0)
			*(*(mat_L+k)+0) = *(*(mat_A+k)+0)/(*(*(mat_U+0)+0)); //L的0列
	}
	for(i=0; i<N; i++) {
		for(int j=i+1; j<N; j++)
			*(*(mat_L+i)+j) = 0; //L上三角为0
		for(j=0; j<i; j++)
			*(*(mat_U+i)+j) = 0; //U下三角为0	
	}

	for(i=1; i<N; i++) {
		//求解U的i行
		for(int j=i; j<N; j++) {
			double sum =0.0;
			for(k=0; k<i; k++)
				sum += *(*(mat_L+i)+k) * (*(*(mat_U+k)+j));
			*(*(mat_U+i)+j) = *(*(mat_A+i)+j)-sum; //j>=i时U元素值
		}
		//求解L的i列
		for(j=i+1; j<N; j++) {
			double sum =0.0;
			for(k=0; k<i; k++)
				sum += *(*(mat_L+j)+k) * (*(*(mat_U+k)+i));
			*(*(mat_L+j)+i) = (*(*(mat_A+j)+i)-sum)/(*(*(mat_U+i)+i)); //j<i时L元素值
		}
	
	}
	//L U创建结束
	std::cout << "output mat_L \n";
	for(i=0; i<N; i++) {
		for(int j=0; j<N; j++)
			std::cout << *(*(mat_L+i)+j) << "  ";
		std::cout << '\n';
	}
	std::cout << "output mat_U \n";
	for(i=0; i<N; i++) {
		for(int j=0; j<N; j++)
			std::cout << *(*(mat_U+i)+j) << "  ";
		std::cout << '\n';
	}
}

///////////////////////////////////////////////////////////////////////////////////
void LU::cal_mat() { //Ax=b 转化成LUx=b, 先求Ly=b,再求Ux=y
	/*----求解Ly=b----*/
	double *mat_y = (double *)malloc(N*sizeof(double));
	*(mat_y+0) = *(mat_b+0);
	for(int i=1; i<N; i++) {
		double sum=0.0;
		for(int j=0; j<i; j++)
			sum += *(*(mat_L+i)+j)* (*(mat_y+j));
		*(mat_y+i) = (*(mat_b+i)-sum)/1.0;
	}
	std::cout << "output mat_y \n";
	for(i=0; i<N; i++)
		std::cout << *(mat_y+i) << "  ";
	std::cout << '\n';
	/*----求解Ly=b----*/

	/*---求解Ux=y---*/
	result = (double *)malloc(N*sizeof(double));
	*(result+N-1) = *(mat_y+N-1)/(*(*(mat_U+N-1)+N-1));
	for(i=N-2; i>=0; i--) {
		double sum=0.0;
		for(int k=N-1; k>i; k--)
			sum += *(*(mat_U+i)+k) * (*(result+k));
		*(result+i) = (*(mat_y+i)-sum)/(*(*(mat_U+i)+i));
	}
	//输出结果
	std::cout << "output result \n";
	for(i=0; i<N; i++)
		std::cout << *(result+i) << "  ";
	std::cout << '\n';
	/*---求解Ux=y---*/
}

 LU.cpp

#include "LU_File.h"

int main() {
	LU lu;
	int n;
	std::cout << "input the dimetion of mat_A ";
	std::cin >> n;
	lu.or_mat(n);
	lu.creatLU();
	lu.cal_mat();

	return 0;
}

 结果验证