C语言实现方程组LU分解法和列主元LU分解法

根据老师所讲LU分解法和列主元LU分解法本质就是高斯消去法和列主元高斯消去法，

将L部分变为都变为0，则就化为高斯消去法，我们在方程求回代的过程中只用到了U和y，不需要用到L，那个L部分实际保存了每一次运算的倍数关系，详细可自己举数学例子验证便可知。

以例子解释规律，我们以（1，1）表示第一行第一列数,则（3，1）表示第三行第一列，那么第一行所有数保持不变，第一列数全部除上第一行第一列。 然后原来（2，2）-(2,1)*(1,2)即：5-0.7*7。 6-0.7*8，5-0.7*7而第二列为（6-0.8*7）/0.1,    (5-0.7*7)/0.1.先进性行运算，再进行列运算，第（3，3）则为2=（10-0.8*8-4*0.4），行则:3=(9-0.8*7-4*0.1),-6=(6-0.8*10-4*1）列则：1.5=（9-0.7*8-1*0.4）/2

大家应该能总结出规律了：所有数字更新完的上三角为U，下三角为L，L的对角线都为1.

代码实现：

#include<stdio.h>
# define M 100
int main()
{
	int i,j,n,k,s;
	double a[M][M],w,x[M],y[M],c,f;
	printf("输入方程未知数的个数：");
	scanf("%d",&n);
	printf("先输入一行，再输入每行。");
	/*录入数据*/ 
	for(i=0;i<n;i++)
	{
		for(j=0;j<n+1;j++)
		{
			scanf("%lf",&a[i][j]);
		}
	}
	/*处理数据*/
	for(i=1;i<n;i++)
	{
		a[i][0]=a[i][0]/a[0][0];//对第一列数据进行处理
	} 
	for(k=1;k<n;k++)
	{
		for(j=k;j<n+1;j++)
		{   w=0;
			for(s=0;s<k;s++)
			{
				w=w+a[k][s]*a[s][j];
			}
			a[k][j]=a[k][j]-w;//对每一行数据进行更新
		}
		for(i=k+1;i<n;i++)
		{
			w=0;
			for(s=0;s<k;s++)
			{
				w=w+a[i][s]*a[s][k];
			}
			a[i][k]=(a[i][k]-w)/a[k][k];//对每一列数据进行更新
		}
	}
	/*检验结果*/ 
	for(i=0;i<n;i++)
	{
		for(j=0;j<n+1;j++)
		{
			printf("%lf  ",a[i][j]);
		}
		printf("\n");
	}
	for(i=0;i<n;i++)
	{
		y[i]=a[i][n];
	}
	/*回代*/ 
	x[n-1]=1.0*y[n-1]/a[n-1][n-1];
    for(i=n-2;i>=0;i--)
    {   f=0;
    	for(j=i+1;j<n;j++)
    	{
    		c=1.0*a[i][j]*x[j];
    		f=c+f;
		}
		x[i]=1.0*(y[i]-f)/a[i][i];	
	}
	for(i=1;i<n+1;i++)
	{
	printf("x[%d]=%lf\n",i,x[i-1]); 
	}
	return 0;
}

 列主元高斯消元法，只是增加了一个判断每一列数的大小情况，先将列变化后的最大值所在那一行进行全部数据换行，然后在进行行和列的变化。

代码实现：

#include<stdio.h>
#include<math.h>
# define M 100
int main()
{
	int i,j,n,k,s,e;
	double a[M][M],w,x[M],y[M],c,f,g[M],h,r[M],max,d;
	printf("输入方程未知数的个数：");
	scanf("%d",&n);
	printf("先输入一行，再输入每行。");
	/*录入数据*/ 
	for(i=0;i<n;i++)
	{
		for(j=0;j<n+1;j++)
		{
			scanf("%lf",&a[i][j]);
		}
	}
	//处理数据
	for(i=1;i<n;i++)
	{
		a[i][0]=a[i][0]/a[0][0];
	}

	for(k=1;k<n;k++)
	{
        for(i=k;i<n;i++)
        {    c=0;
        	for(s=0;s<k;s++)
        	{
			c=c+a[i][s]*a[s][k];//先处理每一列的数
	    	}
	    	a[i][k]=a[i][k]-c;	
		}
		max=a[k][k];
		e=k;
		for(i=k;i<n-1;i++)
		{
			if(fabs(a[i][k])<fabs(a[i+1][k]))//找出处理后的每一列最大值
			   {
			   max=a[i+1][k];	
			   e=i+1;
			   }
		}
		if(e>k)//进行换行
		{
			for(j=0;j<n+1;j++)
			{
				h=a[e][j];
				a[e][j]=a[k][j];
				a[k][j]=h;
			}
	   } 
		for(i=k+1;i<n;i++)
		{
			a[i][k]=a[i][k]/a[k][k];//列变化
		}
		for(j=k+1;j<n+1;j++)
		{   c=0;
			for(s=0;s<k;s++)
			{
				c=c+a[k][s]*a[s][j];
			}
			a[k][j]=a[k][j]-c;//行变化
		}
    } 
	 for(i=0;i<n;i++)
		{
			for(j=0;j<n+1;j++)
			{
				printf("%lf ",a[i][j]);
			}
			printf("\n");
		}
	for(i=0;i<n;i++)
	{
		y[i]=a[i][n];
	}
	//回代 
	x[n-1]=1.0*y[n-1]/a[n-1][n-1];
    for(i=n-2;i>=0;i--)
    {   f=0;
    	for(j=i+1;j<n;j++)
    	{
    		c=1.0*a[i][j]*x[j];
    		f=c+f;
		}
		x[i]=1.0*(y[i]-f)/a[i][i];	
	}
	for(i=1;i<n+1;i++)
	{
	printf("x[%d]=%lf\n",i,x[i-1]); 
	}
	return 0;
	 
}