gauss-seidel.cpp

数值分析或计算方法中的Gauss－Seidel法求解线性方程组的c++程序

//用Gauss-Seidel求解线性方程组
//G_S.cpp
#include "iostream.h"
#include "math.h"

#define N 6
#define errorrange 1e-5

void main()
{
	double y[N],z[N],temp[N],t,m,p,q,max;
	int i,j,k;
	double a[N][N]={{4,-1,0,-1,0,0},{-1,4,-1,0,-1,0},{0,-1,4,-1,0,-1},
		{-1,0,-1,4,-1,0},{0,-1,0,-1,4,-1},{0,0,-1,0,-1,4}};
	double b[N]={0,5,-2,5,-2,6};
	double x[N]={0,0,0,0,0,0};
	k=0;
	max=1.0;         //为2范数赋大于errorrange的初值
	while(max>errorrange)     //当结果小于误差范围时退出循环
	{
		k+=1;
		for(i=0;i<N;i++)
			  temp[i]=x[i];
		t=0;
		for(j=1;j<N;j++)
			t+=a[0][j]*x[j];
		y[0]=(b[0]-t)/a[0][0];
		x[0]=y[0];
		for(i=1;i<N-1;i++)
			{
				m=0;
				for(j=0;j<=i-1;j++)
					 m+=a[i][j]*x[j];
				p=0;
				for(j=i+1;j<N;j++)
					  p+=a[i][j]*x[j];
				y[i]=(b[i]-m-p)/a[i][i];
				x[i]=y[i];
			}
		q=0;
		for(j=0;j<N-1;j++)
			q+=a[N-1][j]*x[j];
		y[N-1]=(b[N-1]-q)/a[N-1][N-1];
		x[N-1]=y[N-1];
		max=0;
		for(i=0;i<N;i++)
		{
			 z[i]=fabs(y[i]-temp[i]);  
			max+=z[i]*z[i];
			max=sqrt(max);
		}
	}

	cout<<"用Gauss-Seidel法求解本方程组得:"<<endl;
	for(i=0;i<6;i++)
		cout<<"y["<<i<<"]="<<y[i]<<endl;
	cout<<"迭代次数k="<<k<<endl;
}