lu.cpp

用LU分解法求解线性方程组的c++程序

//LU.cpp
#define N 5
#include "math.h"
#include "iostream.h"
#include "time.h"
void main()
{
	int i,j,r,k;
	double a[N][N],b[N],x[N],y[N],m,n,t,duration;
	clock_t start,finish;
	static double l[N][N]={0,0};
	static double u[N][N]={0,0};
	
	start=clock();
	b[0]=-27;
	for(i=1;i<N;i++)
		b[i]=-15;
	for(i=0;i<N;i++)
	{
		for(j=0;j<N;j++)
		{
			if(i==j)
				a[i][j]=-4;
			else if(j==i+1)
				a[i][j]=1;
			else if(j==i-1)
				a[i][j]=1;
			else
				a[i][j]=0;
		}
	}

	for(j=0;j<N;j++)
		u[0][j]=a[0][j];
	for(i=1;i<N;i++)
	    l[i][0]=a[i][0]/u[0][0];
	for(r=1;r<N;r++)
		{
		 for(j=r;j<N;j++)
			{
			 m=0;
			 for(k=0;k<=r-1;k++)
				{
				 m+=l[r][k]*u[k][j];
				 u[r][j]=a[r][j]-m;
				}
			}

		 for(i=r+1;i<N;i++)
			 {
			 n=0;
			 for(k=0;k<=r-1;k++)
				{
				 n+=l[i][k]*u[k][r];
				 l[i][r]=(a[i][r]-n)/u[r][r];
				}
			}			 
		}
		y[0]=b[0];
		for(k=1;k<N;k++)
			{
			 t=0;
			 for(j=0;j<=k-1;j++)
				{
				 t+=l[k][j]*y[j];
				 y[k]=b[k]-t;
				}
			}
		x[N-1]=y[N-1]/u[N-1][N-1];
		for(k=N-2;k>=0;k--)
			{
			 t=0;
			 for(j=k+1;j<N;j++)
				 t+=u[k][j]*x[j];
			 x[k]=(y[k]-t)/u[k][k];				 
			}
		finish=clock();
		duration=(double)(finish-start)/CLOCKS_PER_SEC;
		cout<<"利用LU分解法解方程组得："<<endl;
		for(i=0;i<N;i++)
		{
			cout<<"x["<<i+1<<"]="<<x[i]<<" ";
			if(((i+1)%3)==0)
				cout<<endl;
		}
		cout<<endl;
		cout<<"程序运行时间为"<<duration<<"秒"<<endl;
		
}