g_s_iter.cpp

计算方法程序常微分方程的数值解法课堂讲义

//////////////////////////////////////////////////////////////////////
//	程序6.2   Gauss-Seidel 迭代

#include	<stdio.h>
#include	<stdlib.h>
#include	<conio.h>
#include	<math.h>

#define	MAX_n	100                //本程序能处理的方程最大阶数
#define PRECISION	0.0000001      //迭代精度,相邻二步解之差的最大范数是不是小于这一值，小于则迭代结束
#define MAX_Number	1000           //最大迭代次数,超过则算发散

//输入n元向量到数组x中
void VectorInput(float x[],int n)
{
	int i;

	for(i=1;i<=n;++i)
	{
		printf("x[%d]=",i);
		scanf("%f",&x[i]);
	}
}

//输入m X n 阶阵 到二元数组 A[][] 中
void MatrixInput(float A[][MAX_n],int m,int n)
{
	int i,j;
	printf("\n===Begin input Matrix elements===\n");
	for(i=1;i<=m;++i)
	{
		printf("Input_Line %d : ",i);
		for(j=1;j<=n;++j)
			scanf("%f",&A[i][j]);
	}
}

//输出n元向量数组x
void VectorOutput(float x[],int n)
{
	int i;
	for(i=1;i<=n;++i)
		printf("\nx[%d]=%f",i,x[i]);
}

//前后二步迭代之差的最大范数是不是小于精度要求
int IsSatisfyPricision(float x1[],float x2[],int n)
{
	int i;

	for(i=1;i<=n;++i)
		if(fabs(x1[i]-x2[i])>PRECISION) return 1;
	return 0;
}

//G-S迭代,n为阶数，x为初始向量，A为方程组系数的扩展矩阵
//请与Jacobi迭代的程序比较，只有一点差别
int Jacobi_(float A[][MAX_n],float x[],int n)
{
	float x_former[MAX_n]; //放前一步结果向量
	int i,j,k;

	printf("\nInput vector x0:\n");
	VectorInput(x,n);     

	k=0;//当前迭代步
	do{
		//显示上一步结果
		for(i=1;i<=n;++i)
		{
			printf("\nx[%d]=%f",i,x[i]);
			x_former[i]=x[i];
		}
		printf("\n");

		//迭代计算
		for(i=1;i<=n;++i)
		{
			x[i]=A[i][n+1];
			for(j=1;j<=n;++j)
				if(j!=i)x[i]-=A[i][j]*x[j];//Only defference:here x[j] replace x_former[j]
			if(fabs(A[i][i])>PRECISION)
				x[i]/=A[i][i];
			else
				return 1; //满足精度，迭代结束
		}
		++k;
	}while(IsSatisfyPricision(x,x_former,n) && k<MAX_Number);

	if(k>=MAX_Number)
		return 1; //发散
	else
	{
		printf("\nG-S %d times!",k); //收敛，显示迭代了多少步
		return 0;
	}
}

void main()
{
	int n;
	float A[MAX_n][MAX_n],x[MAX_n];

	//输入阶数
	printf("\nInput n=");
	scanf("%d",&n);
	if(n>=MAX_n-1) //不能超过最大充许阶数，否则数组A放不下
	{
		printf("\n\007n must <%d!",MAX_n);
		exit(0);
	}

	//输入方程组增广（扩展）矩阵
	MatrixInput(A,n,n+1);

	//调用G-S迭代函数
	if(Jacobi_(A,x,n))
		printf("\nG-S Failed!");
	else
	{
		printf("\nOutput Solution:");
		VectorOutput(x,n); //输出结果
	}
	printf("\n\n\007Press any key to quit!\n");
	getch();
}

/*
	运行实例：(注意：输入的是方程组的增广系数矩阵)

Input n=3

===Begin input Matrix elements===
Input_Line 1 : 2 -1 -1 4
Input_Line 2 : 3 4 -2 11
Input_Line 3 : 3 -2 4 11

Input vector x0:
x[1]=0
x[2]=0
x[3]=0
  .
  .
  .
x[1]=3.000000
x[2]=1.000000
x[3]=1.000000

G-S 21 times!
Output Solution:
x[1]=3.000000
x[2]=1.000000
x[3]=1.000000

Press any key to quit!
*/