cmpguassseuder.cpp

高斯赛德尔迭代法 通过矩阵的迭代计算 不断逼近 精确值

#include"iostream.h"
#include"math.h"
#include"stdio.h"
const double pi=3.141592653;
double **Alloc2D(int nx,int ny);
double magnet(double **a,double **b,double nx,double ny);
void CmpGuassSeuder(double **phi,int nx,int ny,double ratio,double eps);
void main()
{
	double **phi;
	int nx=50,ny=50;
	//double ratio=2/(1+sin(pi/(nx-1)));
	double ratio=1;
	double eps=1e-8;
	phi=Alloc2D(nx,ny);
	for(int i=0;i<nx;i++)
	{
		phi[ny-1][i]=10;
	}
	CmpGuassSeuder(phi,nx,ny,ratio,eps);	
}
void CmpGuassSeuder(double **phi,int nx,int ny,double ratio,double eps)
{
	FILE *p;
    p=fopen("data.txt","w+");
	double **b;
	int i,j=1;
	b=Alloc2D(nx,ny);
	while(1)
	{
		for(i=1;i<nx-1;i++)
		{
			phi[i-1][j]=b[i-1][j];//把phi矩阵更新，若注释了则成jacobi法了
			for(j=1;j<ny-1;j++)
			{
				phi[i][j-1]=b[i][j-1];//把phi矩阵更新，若注释了则成jacobi法了
				b[i][j]=.25*ratio*(phi[i+1][j]+phi[i][j+1]+phi[i-1][j]+phi[i][j-1]);
			}
		}
			if(magnet(phi,b,nx,ny)<eps) break;
		for( i=1;i<nx-1;i++)
		{
			for( j=1;j<ny-1;j++)
			{phi[i][j]=b[i][j];}
		}
	}
    for(i=1;i<nx-1;i++)
		{
			for(j=1;j<ny-1;j++)
			{
				cout<<b[i][j]<<"\t";
				fprintf(p,"%f\t",b[i][j]);
			}
			cout<<endl;
			fprintf(p,"\n");
		}

}
double magnet(double **a,double **b,double nx,double ny)//计算两数组的差的平方和在开方
{
	int i,j;
	double eps=0,**c;
	c=Alloc2D(nx,ny);
	for( i=1;i<nx-1;i++)
		{
			for(j=1;j<ny-1;j++)
			{c[i][j]=b[i][j]-a[i][j];}
		}
	for(i=1;i<nx-1;i++)
		{
			for(j=1;j<ny-1;j++)
			{eps=eps+c[i][j]*c[i][j];}
		}
	return sqrt(eps);
}
double **Alloc2D(int nx,int ny)
{
	double **D;
	D=new double*[nx];
	for(int i=0;i<nx;i++)
	{
		D[i]=new double[ny];
		for(int j=0;j<ny;j++)
		{
			D[i][j]=0;
		}		
	}	
	return D;
}