数值分析：解线性方程组2

定理：

A是一个n*n阶的矩阵，如果A的顺序主子式不为0，则A可以唯一分解为两个三角矩阵相乘，即A=LU，其中L是单位下三角阵，U是上三角矩阵。

关于矩阵的LU分解，我以前是没有涉及过的。但其实也不难想明白，这和矩阵化为行阶梯矩阵的过程还是相通的。一个n阶矩阵A经过行变换，最终变成一个行阶梯矩阵，也就是一个上三角矩阵。那些行变换，可以记录下来（按照“左乘行变”的原理），所有变换的操作组成一个单位下三角矩阵。这个单位下三角矩阵左乘A，得到那个行阶梯矩阵。在对单位下三角求逆，A的LU分解就得到了。

LU分解的目的还是要解线性方程组。A = LU之后，那么A x = b可以写成LUx =b。接下来，令Ux=y,再代入得Ly = b,通过一次回代过程，就可以把y求出。然后，Ux=y，再用一次回代过程，就解出x了。

下面，来看A=LU，怎么计算L和U：

可见，选定主元（主对角线位置）为循环的参数[这样表述好像不太严谨]，那么该循环内，依次计算主元所在行的u，所在列的l。也就是rr位置的元素，计算第r行u,再计算第r列l。具体计算可代公式。要是再细致的看一下参数，可见不管是计算u还是计算l，计算所需的参数都在其行上，列前位置。所以，也为了节省存储空间，可以把计算得到的u 和 l存入到A中。不影响结果。

好了，下面就浅看一下代码：

//LU  直接三角分解法 

#include<stdio.h>
#define N 3

int main()
{
	double a[N][N] = {{1,2,3},{2,5,2},{3,1,5}};
	double b[N] = {14,18,20}; 
	int i,j,k;
	double sum;
	
	for(i = 0; i<N ; i++)  // 对于u的第i行，l的第i列 
	{
		//计算u (横着往后走) 
		for(j = i ; j<N ; j++){  //u的第i行,j代表列 
			sum = 0;
			for(k = 0 ; k<i ; k++){
				sum = sum + a[i][k]*a[k][j];
			}
			a[i][j] = a[i][j] -sum; 
		}
		
		//计算l (竖着往下走) 
		for(j = i+1 ; j<N ;j++){   //l的第i列,j代表行 
			sum = 0;
			for(k = 0; k<i ;k++)
			{
				sum = sum + a[j][k]*a[k][i]; 
			}
			a[j][i] = (a[j][i] - sum)/a[i][i]; 
		} 
	}
	
	printf("输出LV合起来的矩阵：\n");
	for(i=0;i<N;i++){
		for(j=0;j<N;j++)
		{
			printf(" %.2f ",a[i][j]);
		} 
		printf("\n");
	}
	
	//Ax=b LUx=b 令Ux=y
	//则 Ly=b ,回代求y
	double y[N];
	y[0] = b[0];
	for(i = 1; i<N ; i++)
	{
		double t=0;
		for(j=0 ; j<i ; j++)
		{
			t = t + a[i][j]*y[j];
		}
		y[i]=b[i] -t; //由上向下的计算解的过程 
	}

    printf("输出过程变量:\n");	
	for(i=0;i<N;i++)
	    printf("%.2f ",y[i]); 
	printf("\n");
	
	// Ux=y ,回代求x
	double x[N];
	x[N-1] = y[N-1]/a[N-1][N-1];
	for(i = N-2; i >= 0 ; i--)
	{
		double t=0;
		for(j=i+1 ; j<N ; j++)
		{
			t = t + a[i][j]*x[j];
		}
		x[i]=(y[i]-t) / a[i][i]; //由下向上的计算解的过程 
	}
	
	printf("输出解：\n"); 
	for(i=0;i<N;i++) //输出解 
	    printf("%.2f ",x[i]); 
	 
	return 0;
} 

输出结果为：

------------------------------------------------------------------------------------------------------------------

列选主元/ 全选主元是考虑了计算因子计算的除法过程中，除数不能太小的问题。这里，在计算l的过程中，同样涉及到除法。也同样，考虑选最大主元的必要性。在主元位置每次计算第r行的u和第r列的l之前，选列选主元，当前主元所在列最大的元素作为主元，进行行交换。其实，只是加了这一步骤，其余不变。

代码如下：

//LU  选主元三角分解法 

#include<stdio.h>
#include<math.h>
#define N 3
void Columnselect(double a[][N],double b[N],int i);
void Change2row(double a[][N],double b[N],int r1,int r2);

int main()
{
	
	double a[N][N] = {{0.001,2,3},{-1,3.712,4.623},{-2,1.072,5.643}};  
	double b[N] = {1,2,3}; 
	/*
	double a[N][N] = {{1,2,3},{2,5,2},{3,1,5}};
	double b[N] = {14,18,20}; 
	*/
	int i,j,k;
	double sum;
	 
	for(i = 0; i<N ; i++)  // 对于u的第i行，l的第i列 
	{
	    Columnselect(a,b,i);//列选主元 
	     
		//计算u (横着往后走) 
		for(j = i ; j<N ; j++){  //u的第i行,j代表列 
			sum = 0;
			for(k = 0 ; k<i ; k++){
				sum = sum + a[i][k]*a[k][j];
			}
			a[i][j] = a[i][j] -sum; 
		}
			
		
		//计算l (竖着往下走) 
		for(j = i+1 ; j<N ;j++){   //l的第i列,j代表行 
			sum = 0;
			for(k = 0; k<i ;k++)
			{
				sum = sum + a[j][k]*a[k][i]; 
			}
			a[j][i] = (a[j][i] - sum)/a[i][i]; 
		} 
	}
	
	//Ax=b LUx=b 令Ux=y
	//则 Ly=b ,回代求y
	double y[N];
	y[0] = b[0];
	for(i = 1; i<N ; i++)
	{
		double t=0;
		for(j=0 ; j<i ; j++)
		{
			t = t + a[i][j]*y[j];
		}
		y[i]=b[i] -t; //由上向下的计算解的过程 
	}
	
	// Ux=y ,回代求x
	double x[N];
	x[N-1] = y[N-1]/a[N-1][N-1];
	for(i = N-2; i >= 0 ; i--)
	{
		double t=0;
		for(j=i+1 ; j<N ; j++)
		{
			t = t + a[i][j]*x[j];
		}
		x[i]=(y[i]-t) / a[i][i]; //由下向上的计算解的过程 
	}
	
	printf("输出解：\n"); 
	for(i=0;i<N;i++)   //输出解 
	    printf("%.2f ",x[i]); 
	
	return 0;
} 


void Columnselect(double a[][N],double b[N],int i) //列选主元 
{
	double sum; int j;
	double s[N] = {0}; //选主元用的
	for(j = i ; j<N ;j++){   //i代表列,j代表行 
		sum = 0;
		for(int k = 0; k<i ;k++) {
			sum = sum + a[j][k]*a[k][i]; 
		}
		s[j] = a[j][i] - sum ; 
	}	
		
	//最大主元 ,主元位置i i 
	int max_index = i;
	double max = fabs(s[i]);
	for(int t = i+1 ; t < N ; t++)
	{
		if (fabs(s[t]) > max)  {  max = fabs(s[t]);  max_index = t; }  //记录最大s[i]所在的序号 
	} 
	if( fabs(s[i]) != max) Change2row(a,b,max_index,i); //交换两行 
}


void Change2row(double a[][N],double b[N],int r1,int r2)//交换两行(r1和r2) 
{
	double tmp;
	for(int i=0;i<N;i++)
	{
		tmp=a[r1][i];
		a[r1][i]=a[r2][i];
		a[r2][i]=tmp;
	}
	tmp = b[r1];
	b[r1] = b[r2];
	b[r2] = tmp; 
}

输出结果：

此外，因为我是将A和b分开来写的，换行的时候同时都要换，一开始我就是遗忘这一点。

这里还提供使用类的写法，在project中：

Select_LU.h

#include <iostream>

class LU {

private :
	double** mat_A, *mat_b; //mat_A是系数矩阵，mat_b是结果矩阵
	double* result; //解数组
	int N; //方阵的阶数

public :
	LU(); //构造函数

public : //一些成员函数
	void or_mat(void); //初始矩阵(系数矩阵mat_A, 右值矩阵mat_b)初始化,N为方阵阶
	void creatLU(void); //将矩阵A 分解为L U
	void cal_mat(void); //用LU分解法求Ax=b的解
	void Columnselect(int i); //列选主元
	void Change2row(int r1,int r2); //交换两行

};

 Select_LU.cpp

#include <iostream>
#include <string.h>
#include <malloc.h>
#include <stdlib.h>
#include <math.h>
#include "Select_LU.h"

using namespace std;

//构造函数
LU::LU() {
	mat_A = NULL;
	mat_b = NULL;
}

//初始矩阵(系数矩阵mat_A, 右值矩阵mat_b)初始化,N为方阵阶
void LU::or_mat() {
	char fileName[20];
	cout << "input filrname of Matrix" << endl;
	cin >> fileName;
	FILE *fp;
	char ch;
	if((fp=fopen(fileName,"rt"))==NULL) { //covariance文件中每一行代表一个样本向量,列数为维数
		printf("Cannot open file strike any key exit!");
		exit(0);
	}
	ch=fgetc(fp);
	string str="";

	//============read the rowNum and colNum=============//
	while(ch==' ') ch=fgetc(fp); //往后读，直到非空格
	while(ch!=' ') {
		str+=ch;
		ch=fgetc(fp);
	}
	double aa=atof(str.c_str()); //将字符串转换成double型
	N = aa; //阶数
	str="";
	while(ch==' ') ch=fgetc(fp);
	while(ch!=' ' && ch!='\n') {
		str+=ch;
		ch=fgetc(fp);
	}
	aa=atof(str.c_str()); //将字符串转换成double型
	N = aa;
	str="";
	//============read the rowNum and colNum=============//

	std::cout << "方阵的阶数为 " << N << endl;
	std::cout << "read the Matix file" << endl;

	//开辟空间
	mat_A = (double **)malloc(N*sizeof(double));
	for(int i=0; i<N; i++)
		*(mat_A+i) = (double *)malloc(N*sizeof(double));
	mat_b = (double *)malloc(N*sizeof(double));

	//============read matrix of index and result=============//
	ch=fgetc(fp);
	while(ch != EOF) //没有读取到最后
    {
		for(int i=0; i<N; i++) {
			while(ch==' ' || ch=='\n')
				ch=fgetc(fp);

			//读取index
			for(int j=0; j<N; j++) {
				while(ch==' ')
					ch=fgetc(fp);
				while(ch!=' ') {
					str+=ch;
					ch=fgetc(fp);
				}
				aa=atof(str.c_str()); //将字符串转换成double型
				*(*(mat_A+i)+j) = aa;
				str="";
			}
			//读取index结束

			while(ch==' ')
				ch=fgetc(fp);
			while(ch!=' ' &&  ch!='\n') {
				str+=ch;
				ch=fgetc(fp);
			}
			aa=atof(str.c_str()); //将字符串转换成double型
			*(mat_b+i) = aa; //读取mat_b+i
			str="";
		}
		ch=fgetc(fp);
	}
	//============read matrix of index and result=============//
	std::cout << "read the matix over" << '\n';
	fclose(fp);
}


//进行LU分解
void LU::creatLU()
{
    /*
    改：从（0，0）就开始选主元
	//U的第一行就是原矩阵的第一行, 不用处理
	//L的第一列为源矩阵的第一列除以原矩阵的第一个元素

	for(int i=0; i<N; i++)
		*(*(mat_A+i)+0) = *(*(mat_A+i)+0) / *(*(mat_A+0)+0); //处理第一列
    */

	//从i(i=0~N-1)开始依次处理U的i行和L的i列
	double sumU;
	double sumL;

	for(int i=0; i<N; i++) {

        Columnselect(i);  //进行选主元

		//U的i行
		for(int j=i; j<N; j++) {
			sumU=0.0;
			for(int k=0; k<i; k++)
				sumU += *(*(mat_A+i)+k)* (*(*(mat_A+k)+j));
			*(*(mat_A+i)+j) = *(*(mat_A+i)+j)-sumU;
		}

		//L的i列
		for(int k=i+1; k<N; k++) {
			sumL=0.0;
			for(int j=0; j<i; j++)
				sumL += *(*(mat_A+k)+j)* (*(*(mat_A+j)+i));
			*(*(mat_A+k)+i) = (*(*(mat_A+k)+i)-sumL) / (*(*(mat_A+i)+i));

		}
	}//L U创建结束
	//注意：由于在原矩阵上进行l u的存储，故到此求出的结果对角线上的元素为U的对角线的元素;
	//L 的对角线元素都为1，L对角线下的元素已经保存在此时的对角线下的相应位置了

	std::cout << "output LU \n";
	for(int i=0; i<N; i++) {
		for(int j=0; j<N; j++)
			std::cout << *(*(mat_A+i)+j) << "  ";
		std::cout << '\n';
	}
}

//做两次迭代然后求解
//Ax=b 转化成LUx=b, 先求Ly=b,再求Ux=y
void LU::cal_mat() {

	/*----求解Ly=b----*/

	//开辟y空间
	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_A+i)+j)* (*(mat_y+j));
		*(mat_y+i) = (*(mat_b+i)-sum)/1.0;
	}
	std::cout << "output mat_y \n";
	for(int 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_A+N-1)+N-1));

	for(int i=N-2; i>=0; i--) {
		double sum=0.0;
		for(int k=N-1; k>i; k--)
			sum += *(*(mat_A+i)+k) * (*(result+k));
		*(result+i) = (*(mat_y+i)-sum)/(*(*(mat_A+i)+i));
	}

	//输出结果
	std::cout << "output result \n";
	for(int i=0; i<N; i++)
		std::cout << *(result+i) << "  ";
	std::cout << '\n';
	/*---求解Ux=y---*/
}

void LU::Columnselect(int i) //列选主元
{
	double sum; int j;
	double s[N] = {0}; //选主元用的
	for(j = i ; j<N ;j++){   //i代表列,j代表行
		sum = 0;
		for(int k = 0; k<i ;k++) {
			sum = sum + *(*(mat_A + j)+k) * (*(*(mat_A + k)+i));
		}
		s[j] = *(*(mat_A+j)+i) - sum ;
	}

	//最大主元 ,主元位置i i
	int max_index = i;
	double max = fabs(s[i]);
	for(int t = i+1 ; t < N ; t++)
	{
		if (fabs(s[t]) > max)  {  max = fabs(s[t]);  max_index = t; }  //记录最大s[i]所在的序号
	}
	if( fabs(s[i]) != max) Change2row(max_index,i); //交换两行
	printf("%d---\n",max_index);
}

void LU::Change2row(int r1,int r2)//交换两行(r1和r2)
{
	double tmp;
	for(int i=0;i<N;i++)
	{
		tmp= *(*(mat_A + r1) + i);
		*(*(mat_A + r1)+ i)= *(*(mat_A + r2)+i);
		*(*(mat_A + r2)+i)=tmp;
	}
	tmp = *(mat_b+r1);
	*(mat_b + r1) = *(mat_b+r2);
	*(mat_b + r2) = tmp;
}

main.cpp

#include "Select_LU.h"

int main() {

	LU lu;
	lu.or_mat();  //读入系数矩阵和结果矩阵
	lu.creatLU();  //进行选主元LU分解
	lu.cal_mat();  //进行两次回代
	return 0;

}

结果也是正确的。

.txt文件 的举例如下（最后也要有换行，否则读入那里就出问题，和代码中条件有关）