jacobi.cpp

一般问题的线性插值和JACOBI算法

/*
  Name: Jacobi
  Copyright: 
  Author: 陈国栋  学号 200131000001  01级数学基地班 
  Date: 29-09-04 23:55
  Description: 
*/
# include <iostream.h>
# include <math.h>
# define  max(a,b) (a)>(b)?(a):(b)
# define  MAX  100
double a[MAX][MAX],b[MAX],x0[MAX],x1[MAX];
void Jacobi(int n,double a[][MAX],double b[],double x0[],double delta,double eps,int k,int M)
{
	double s,temp;
	int i,j;
lap1:if(k<=M){
		delta=0;
		for(i=1;i<=n;i++){
			s=0;
			for(j=1;j<=n;j++) s+=a[i][j]*x0[j];
			x1[i]=(b[i]-s)/a[i][i]+x0[i];
			if((temp=fabs(x0[i]-x1[i]))>delta) delta=temp;
		}
		k++;
		if(delta+0.00005>eps){
			for(i=1;i<=n;i++) x0[i]=x1[i];
			goto lap1;
		}
		else{
			cout<<"ans: ";
			for(i=1;i<n;i++) cout<<x1[i]<<" ";
			cout<<x1[i]<<endl;
			cout<<"k: "<<k<<endl;
			goto lap2;
		}
	}
	else cout<<"too much times!"<<endl;
lap2: ;
}
int main()
{
	int i,j,n,M;
	double eps;
	n=3;  eps=1e-1; M=100;
	for(i=1;i<=n;i++){
		b[i]=0;
		for(j=1;j<=n;j++){
			a[i][j]=max(i,j);
			b[i]+=a[i][j];
		}
	}
	for(i=1;i<=n;i++) x0[i]=0;
	Jacobi(n,a,b,x0,0,eps,0,M);
	return 0;
}