gs.cpp

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

#include "iostream.h"
#include "math.h"
#include  "stdio.h"
#define N 3 //方程的阶数
#define MaxK   100 //最大迭代次数
#define EPS   0.5e-4 //精度控制

static double aa[N][N]={{10,-1,-2},{-1,10,-2},{-1,-1,5}};
static double bb[N]={7.2,8.3,4.2};

void main()
{
	int i,j;
    double x[N+1];   
	double a[N+1][N+1],b[N+1];
	int GaussSeidel(double a[][N+1],double b[],double x[]);

	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];
		 x[i]=0;
	 }
	 if(GaussSeidel(a,b,x)==1)
	 {
		 printf("the roots is:");
		 for(i=1;i<=N;i++)
			 printf("  x[%d]=%f   ",i,x[i]);printf("\n");
	 }
	 else printf("\nthe G-S method failed!\n");  
}


int GaussSeidel(double a[][N+1],double b[],double x[])
//a传入系数矩阵(有效元素为a[1,1]...a[n,n]),b为方程组右边常数列，x传入迭代初值并返回解向量
//x**(k+#x)=Gx**(k)+f
{int k=1,i,j;
 double m,max,t[N+1];
 while(true)
 {printf("k=%d:",k);
  max=0.0;
  for(i=1;i<=N;i++){if(a[i][i]==0)
                       return -1;//存在元素akk为0
                    m=0.0;
					t[i]=x[i];
					for(j=1;j<=N;j++)
						if(j!=i)
						  m+=a[i][j]*x[j];
					x[i]=(b[i]-m)/a[i][i];
					if(max<fabs(x[i]-t[i]))
						max=fabs(x[i]-t[i]);
                    printf("   x[%d]=%f   ",i,x[i]);}
 printf("\n");
 if(max<EPS)
	 return 1;//正常结束
 else
	 k++;
 if(k>MaxK)
	 return -2;//迭代次数越界
 }
}