[数值计算c++]解线性方程的直接方法之改进平方根法

对于解线性方程的直接方法，本文一共采用五种方法：Gauss顺序消去法，Gauss列主元消去法、平方根法，改进平方根法，追赶法。
对于每种方法的具体原理请参考：https://zhuanlan.zhihu.com/p/195732053注：本文的部分方法来自书本，部分来自上述链接中matlab原理改编。

Gauss 顺序消元法

//高斯顺序消元法 ax=b n代表方阵维数
void GuassEliminination(double a[][1000], double b[],int n)
{
    int k=1,i,j;
    double s;
    while(k<=n)
    {
    	for(j=k+1;j<=n;j++)a[k][j] = a[k][j]*1.0/a[k][k];
    	b[k]=b[k]*1.0/a[k][k];
    	for(i=k+1;i<=n;i++)
    	{
    		for(j=k+1;j<=n;j++)a[i][j] = a[i][j] - a[i][k]*a[k][j];
    		b[i] = b[i] - a[i][k]*b[k];
		}
		k++;
	}
	for(i=n-1;i>=1;i--)
	{
		s=0;
		for(j=i+1;j<=n;j++)
			s+=a[i][j]*b[j];
		b[i] = b[i] -s;	
	} 
	for(i=1;i<=n;i++)cout<<b[i]<<endl; //三角分解式为a[][]上三角
}

Gauss列主元消去法

void GuassElimininationWithPartialPivoting(double a[][1000], double b[],int n)
{
    int k=1,i,j,t;//t记录每次最大的列 
    double ave=0,tmp=0,s=0;//记录最大的数值 
    while(k<=n)
    {
    	ave = a[k][k];
    	for(i=k+1;i<=n;i++)
    	{
    		if(fabs(a[i][k])>ave)
    		{
    			ave = fabs(a[i][k]);
    			t = i;
			}
		}
		if(a[k][k] == 0)cout<<"error"<<endl;
		if(k!=t)
		{
			for(i=k;i<=n;i++)
			{
				tmp = a[k][i];
				a[k][i] = a[t][i];
				a[t][i] = tmp;
			}
			tmp = b[k];
			b[k] = b[t];
			b[t] = tmp;
		}
    	for(j=k+1;j<=n;j++)a[k][j] = a[k][j]*1.0/a[k][k];
    	b[k]=b[k]*1.0/a[k][k];
    	for(i=k+1;i<=n;i++)
    	{
    		for(j=k+1;j<=n;j++)a[i][j] = a[i][j] - a[i][k]*a[k][j];
    		b[i] = b[i] - a[i][k]*b[k];
		}
		k++;
	}
	for(i=n-1;i>=1;i--)
	{
		s=0;
		for(j=i+1;j<=n;j++)
			s+=a[i][j]*b[j];
		b[i] = b[i] -s;	
	} 
	for(i=1;i<=n;i++)cout<<b[i]<<endl; //三角分解式为a[][]上三角
}

平方根法

void cholesky(double a[][1000],double b[],int n)
{
	int i,j,k;
	double s;
	double l[n+1][n+1];
	double x[n+1],y[n+1];
	for(i=0;i<=n;i++)
		for(j=0;j<=n;j++)
			l[i][j] =0;
	for(i=1;i<=n;i++)
	{
		for(j=1;j<=i-1;j++)
		{
			for(s=0,k=1;k<=j-1;k++) s+=l[i][k]*l[j][k];
			l[i][j] = (a[i][j]-s)*1.0/l[j][j];
		}
		s=0;
		for(k=1;k<=i-1;k++) s+=l[i][k]*l[i][k];
		l[i][i] =sqrt(a[i][i]-s);	
	}
	for(i=1;i<=n;i++)
	{
		s=0;
		for(j=1;j<=i-1;j++)
			s+=l[i][j]*y[j];
		y[i] = (b[i] -s)*1.0/l[i][i];
	} 
	for(i=n;i>=1;i--)
	{
		s=0;
		for(j=i+1;j<=n;j++)
			s+=l[j][i]*x[j];
		x[i] = (y[i] -s)*1.0/l[i][i];	
	} 
	for(i=1;i<=n;i++)cout<<x[i]<<endl; //三角分解是l[][]下三角
} 

改进平方根法

void cholesky2(double a[][1000],double b[],int n)
{
	int i,j,k;
	double s;
	double T[n+1][n+1],x[n+1],y[n+1];
	for(i=2;i<=n;i++)
	{
		
		for(j=1;j<=i-1;j++)
		{
			s=0;
			for(k=1;k<=j-1;k++) s+=T[i][k]*a[j][k];
			T[i][j] = (a[i][j]-s)*1.0;
		}
		for(j=1;j<=i-1;j++)a[i][j] = T[i][j]/a[j][j];
		s=0;
		for(k=1;k<=i-1;k++) s+=T[i][k]*a[i][k];
		a[i][i] =a[i][i]-s;	
	}
	y[1] = b[1];
	for(i=2;i<=n;i++)
	{
		s=0;
		for(k=1;k<=i-1;k++)
			s+=a[i][k]*y[k];
		y[i] = (b[i] -s);
	} 
	x[n] = y[n]*1.0/a[n][n];
	for(i=n-1;i>=1;i--)
	{
		s=0;
		for(k=i+1;k<=n;k++)
			s+=a[k][i]*x[k];
		x[i] = (y[i]/a[i][i] -s);	
	} 
	for(i=1;i<=n;i++)cout<<x[i]<<endl;三角分解是a[][]下三角
} 

追赶法

void myrun(double a[][1000],double b[],int n)
{
	int i;
	double aa[n+1],bb[n+1],cc[n+1];
	double beta[n];
	double x[n+1],y[n+1];
	aa[1] =0;
	for(int i=1;i<=n;i++)
	{
		bb[i] = a[i][i];
		if(i!=n)aa[i+1] = a[i+1][i];
		if(i!=1)cc[i-1] = a[i-1][i];
	}
	beta[1] = cc[1]*1.0/bb[1];
	for(i=2;i<=n-1;i++)
		beta[i] = cc[i]*1.0/(bb[i]-aa[i]*beta[i-1]);
	y[1] = b[1]/bb[1];
	for(i=2;i<=n;i++)
		y[i] = (b[i] - aa[i]*y[i-1])/(bb[i]-aa[i]*beta[i-1]);
	x[n] =y[n];
	for(i=n-1;i>=1;i--)
		x[i] = y[i] - beta[i]*x[i+1];
	for(i=1;i<=n;i++)cout<<x[i]<<endl;	
}