lu.cpp

这是计算方法上的所有算法

#include "math.h"
#include "stdio.h"
#define N 3
static aa[N][N]={{3,1,6},{2,1,3},{1,1,1}};
static bb[N]={2,7,4};

void main()
{
	int i,j;
	double a[N+1][N+1],b[N+1];
	void solve(double a[][N+1],double b[N+1]);int LU(double a[][N+1]);
	for(i=1;i<=N;i++)
	{
		for(j=1;j<=N;j++)a[i][j]=aa[i-1][j-1];
		b[i]=bb[i-1];
	}
	if(LU(a)==1)
	{
		printf("矩阵L如下\n");
		for(i=1;i<=N;i++)
		{
			for(j=1;j<=i;j++)
				if(i==j)printf("%12d  ",1);
				else printf("%12f",a[i][j]);
			printf("\n");
		}
		
		printf("\n矩阵U如下\n");
		for(i=1;i<=N;i++)
		{
			for(j=1;j<=N;j++) 
				if(i<=j)printf("%12f",a[i][j]);
				else printf("%12s","");
			printf("\n");
		}	
		solve(a,b);
		for(i=1;i<=N;i++)printf("x[%d]=%f  ",i,b[i]);printf("\n");
	}
}

int LU(double a[][N+1])//对N阶矩A阵进行LU分解，结果L、U存放在A的相应位置
{int i,j,k,s;
 double m,n;
 for(i=2;i<=N;i++)
	 a[i][1]/=a[1][1];
 for(k=2;k<=N;k++)
 {for(j=k;j<=N;j++){m=0;
                    for(s=1;s<k;s++)
						m+=a[k][s]*a[s][j];
					a[k][j]-=m;}
  if(a[k][k]==0)
	  return -1;//存在元素akk为0
  for(i=k+1;i<=N;i++)
  {n=0;
   for(s=1;s<k;s++)
	   n+=a[i][s]*a[s][k];
   a[i][k]=(a[i][k]-n)/a[k][k];
  }
 }
 return 1;//正常结束
}

void solve(double a[][N+1],double b[N+1])//利用分解的LU求x
//回代求解，L和U元素均在A矩阵相应位置；b存放常数列，返回时存放方程组的解
{double y[N+1],F;
 int i,j;
 y[1]=b[1];
 for(i=2;i<=N;i++)
 {F=0.0;
  for(j=1;j<i;j++)
	  F+=a[i][j]*y[j];
  y[i]=b[i]-F;
 }
 b[N]=y[N]/a[N][N];
 for(i=N-1;i>0;i--)
 {F=0.0;
  for(j=N;j>i;j--)
	  F+=a[i][j]*b[j];
  b[i]=(y[i]-F)/a[i][i];
 }
}