jacobi.cpp

从出版社求得的经典数值算法

//////////////////////////////////////////////////////
//有关注释请参考Gauss-Saidel迭代的程序 G_s_iter.cpp
//

/*	§6  线性方程组的迭代法

	程序6.1  Jacobi 迭代
*/
#include	<stdio.h>
#include	<stdlib.h>
#include	<conio.h>
#include	<math.h>

#define	MAX_n	100
#define PRECISION	0.0000001
#define MAX_Number	1000

void VectorInput(float x[],int n)
{
	int i;

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

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]);
	}
}

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;
}

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_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("\nJacobi %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)
	{
		printf("\n\007n must <%d!",MAX_n);
		exit(0);
	}

	MatrixInput(A,n,n+1);

	if(Jacobi_(A,x,n))
		printf("\nJacobi 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

Jacobi 112 times!
Output Solution:
x[1]=3.000000
x[2]=1.000000
x[3]=1.000000

Press any key to quit!
*/