yuzhifan.cpp

说什么说,你真上由于啊打算啊地方啊地方饿地方广泛

#include<iostream.h>
#include<iomanip.h>
#include<math.h>
#include<stdio.h>
#define  n  4
#define NULL 0

//输入方程组的控制精度
float inputeps( )
{
	float eps;
    cout<<"请输入控制精度(一个小的正数)eps:";
     cin>>eps;
	 if(eps<=0)
	 {
		cout<<"错误，请重新输入一个小的正数:";
		cin>>eps;
	 }
	 return eps;
}

//从磁盘文件写入系数矩阵A
void loadab(float (*a)[n],float b[n])
{FILE *fp;
 int i,j;
 float A[n][n+1];
 if((fp=fopen("数字.txt","rb"))==NULL)
 {
	 printf("cannot open infile\n");
	 return;
 }
 for(i=0;i<n;i++)
	 for(j=0;j<n+1;j++)
		 fscanf(fp,"%f",*(A+i)+j);   //调用fscanf函数用于从磁盘文件写入系数矩阵A
		 fclose(fp);
  for(i=0;i<n;i++)
	 {
	  for(j=0;j<n;j++)	 
		 a[i][j]=A[i][j];
		 b[i]=A[i][n];
	 }
		 return;
 }


//显示方程组
void show(float A[n][n],float b[n],float eps)
{
	int i,j;
	cout<<"要求解的方程组是:\n";
	for(i=0;i<n;i++)
	{
		cout<<A[i][0]<<"X0";
		for(j=1;j<n;j++)
		
		   if(A[i][j]>=0)
				cout<<"+"<<A[i][j]<<'X'<<j;
			 else
				 cout<<A[i][j]<<'X'<<j;
		cout<<"="<<b[i];
		cout<<endl;
	}
  cout<<"控制精度是eps="<<eps<<endl;
  return;
}

//输出方程组的解,以转置矩阵的形式输出
void output(float x[])
{
	int i;
	cout<<"方程组的解是:"<<endl;
	cout<<"x=(";
	for(i=0;i<n-1;i++)
		cout<<x[i]<<',';
	cout<<x[n-1]<<")T"<<endl;
	return;
}

//输出增广矩阵(A,b)
void outarray(float A[][n],float b[])
{
	int i,j;
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
			cout<<A[i][j]<<setw(7);
		cout<<b[i];
		cout<<endl;
	}
	return;
}

//输出矩阵，主要用于输出L矩阵和U矩阵
void outputa(float a[ ][n])
{
	int i,j;
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
			cout<<a[i][j]<<setw(7);
		cout<<endl;
	}
	return;
}


//将n行n+1列矩阵消元为上三角形式
void xiaoyuan(float A[][n],float b[],float eps)
{
	int i,j,k;
	float m;
	for(i=0;i<n-1;i++)
	{
		for(j=i+1;j<n;j++)
		{
			if(A[i][i]!=0)
			{
				m=A[j][i]/A[i][i];
			  for(k=i;k<n;k++) 
				  A[j][k]=A[j][k]-m*A[i][k];		  
		     	b[j]=b[j]-m*b[i];
			}
          else
		  {
		       cout<<"可选用列主元消去法\n";
			   return;
		  }
		}
	}
	return;
}

//回代函数
void  huidai(float A[][n],float b[],float x[])
{
	float sum;
	int i,j;
	if(A[n-1][n-1]!=0)
		x[n-1]=b[n-1]/A[n-1][n-1];
	else
		cout<<"方程组无解\n";
	for(i=n-2;i>=0;i--)
	{
		sum=0;
		for(j=i+1;j<n;j++)
			sum=sum+A[i][j]*x[j];
		if(A[i][i]!=0)
			x[i]=(b[i]-sum)/A[i][i];
		else
	       cout<<"x["<<i<<"]为任意实数\n";
	}
	return;
}

//Gauss消去法
void Gauss(float A[][n],float b[],float x[],float eps)
{
	cout<<"用高斯消去法解方程组\n";
	cout<<"增广阵(A,b)为:\n";
	outarray(A,b);             //输出增广矩阵(A,b)
	xiaoyuan(A,b,eps);
	cout<<"消元后的增广阵:\n";
	outarray(A,b);             //输出消元后的增广矩阵(A,b)
	huidai(A,b,x);
	output(x);
	return;
}

//选列主元并逐步消元求解
void xuanzhu(float A[][n],float b[],float eps)
{
	int i,j,k,row;
	float max,temp,x[n],m[n-1];
	for(i=0;i<n-1;i++)
	{	
		max=(float)fabs(A[i][i]);
		     row=i;	             //row用于标记列最大值的行数
		for(j=i+1;j<n;j++)
		{            
			if(fabs(A[j][i])>max)
			{
				max=(float)fabs(A[j][i]);
				row=j;
			}
			else
				continue;
		}
			if(max<=eps)
			{
				cout<<"方程组无解"<<endl;
				return;
			}
			else
			{
				for(k=0;k<n;k++)        //找出最大值后，互换两行的对应值
				{
					temp=A[i][k];          //temp为中间变量
					A[i][k]=A[row][k];
					A[row][k]=temp;
				}
				temp=b[i];
				b[i]=b[row];
				b[row]=temp;
				//消元
				for(j=i+1;j<n;j++)
				{
					m[i]=A[j][i]/A[i][i];
					for(k=i;k<n;k++)	
						A[j][k]=A[j][k]-m[i]*A[i][k];
					if(fabs(A[j][k])<=eps)
						A[j][k]=0;
					b[j]=b[j]-m[i]*b[i];
				}
			}
	}
	if(fabs(A[n-1][n-1])<=eps)
		cout<<"方程组无解"<<endl;
	else
	{
		cout<<"列主消元后的矩阵为:"<<endl;
		outarray(A,b);
		huidai(A,b,x);    //调用回代函数求解
		output(x);        //输出方程组的解
	}
	return;
}

//列主元消去法
void liezhuyuan(float A[][n],float b[],float eps)
{
   cout<<"用列主元消去法解方程组\n";
  cout<<"增广阵(A,b)为:\n";
  outarray(A,b);
  xuanzhu(A,b,eps);
  return;
}

//将A分解为LU矩阵
void fenjie(float A[][n],float l[][n],float u[][n])
{
	int i,k,r;
	float sum=0;
	for(i=0;i<n;i++)              //求解LU矩阵
	{
		u[0][i]=A[0][i];
		l[i][0]=A[i][0]/u[0][0];
	}
	for(r=1;r<n-1;r++)
	{
		for(i=r;i<n;i++)
		{
			for(k=0;k<r;k++)
				sum=sum+l[r][k]*u[k][i];
			u[r][i]=A[r][i]-sum;
			sum=0;
		}
		for(i=r+1;i<n;i++)
		{
			for(k=0;k<r;k++)
				sum=sum+l[i][k]*u[k][r];
			l[i][r]=(A[i][r]-sum)/u[r][r];
			sum=0;
		}
	}
	if(r==n-1)
	{
		for(k=0;k<r;k++)
			sum=sum+l[r][k]*u[k][r];
		u[r][r]=A[r][r]-sum;
		sum=0;
	}
	return;
}

//根据LU矩阵求解
void qiujie(float l[][n],float u[][n],float b[],float x[])
{
	float y[n],sum;
	int i,r;
	y[0]=b[0];
	for(r=1;r<n;r++)          //求LY=b的解Y
	{
		sum=0;
		for(i=0;i<r;i++)
			sum=sum+l[r][i]*y[i];
		y[r]=b[r]-sum;
	}
	//输出Y的解，以转置矩阵的形式输出
	cout<<"LY=b,UX=Y,其中Y的解为:"<<endl;
    cout<<"y=(";
	for(i=0;i<n-1;i++)
		cout<<y[i]<<',';
	cout<<y[n-1]<<")T"<<endl;
    //求UX=Y的解X
	x[n-1]=y[n-1]/u[n-1][n-1];
	for(r=n-2;r>=0;r--)
	{
		sum=0;
		for(i=r+1;i<n;i++)
			sum=sum+u[r][i]*x[i];
		x[r]=(y[r]-sum)/u[r][r];
	}
	return;
}

//Doolittle分解法
void Doolittle(float A[][n],float b[],float x[])
{
	float l[n][n],u[n][n];
	int i,j;
	for(i=0;i<n;i++)     //给LU矩阵赋初值
	{
		for(j=0;j<n;j++)
		{
			u[i][j]=0;
			if(j==i)
				l[i][j]=1;
			l[i][j]=0;
		}
	}
	cout<<"用分解的方法解方程组\n";
	fenjie(A,l,u);
	cout<<"分解为LU矩阵积，其中L矩阵为:\n";
	outputa(l);           //输出L矩阵
	cout<<"U矩阵为:\n";
	outputa(u);           //输出U矩阵
    qiujie(l,u,b,x);
	output(x);             //输出方程组的解
	return;
}

//主函数
void main( )
{
	int choic,flag=1;
	float A[n][n],b[n];
	float x[n]={0},eps;
	char ch;
	loadab(A,b);
	eps=inputeps( );
	cout<<"1:Gauss消去法"<<endl; 
	cout<<"2:列主元高斯消去法"<<endl;
	cout<<"3:三角分解法(Doolittle分解)"<<endl;
	cout<<"4:退出"<<endl;
    show(A,b,eps);    //显示方程组
	while(flag)
	{
		cout<<"请输入你的选择choic(1~4):";
			cin>>choic;
		if(choic>4||choic<1)
		{
			cout<<"错误!请重新输入choic(1~4):"<<endl;
			cin>>choic;
		}
		switch(choic)
		{
		    case 1: loadab(A,b);                                                 
			        Gauss(A,b,x,eps);        
				    break;			                
			case 2: loadab(A,b);
                    liezhuyuan(A,b,eps); 
				    break;
			case 3: loadab(A,b);                  
				    Doolittle(A,b,x);  
				    break;
			case 4: return;
		}
		cout<<"想再用其他方法吗(Y/N):";
		cin>>ch;
		if(ch=='Y'||ch=='y')
			flag=1;
		else 
			flag=0;
	}
	return;
}