病态方程组.cpp

病态方程组的运算

#include<iostream.h>
#include<math.h>
#include<process.h>
const int N=101;
double inner(double a[N],double b[N]){
	int i;
	double temp=0;
	for(i=1;i<N;i++)
		temp+=a[i]*b[i];
	if(temp==0)  {
		cout<<"程序中止！";
		exit(1);
	}
	return temp;                        //返回两向量的内积值
}
double maxe(double a[N]){               //求一个向量中的最大元素，以求两个相邻Ｘ值的最大元素差值
	int i;
	double temp=0;
	for(i=1;i<N;i++)
		if(temp<fabs(a[i]))
			temp=fabs(a[i]);
	return temp;
}
void main(){
	int i,j,k;
	double m,t,s,e=0.0001,maxd=1;       //e为要求的误差，maxd为两个相邻Ｘ值的最大元素差值
	double x[N],r[N],p[N],b[N],Ap[N],A[N][N];
    for(i=1;i<N;i++){                   //开始进行运算
		x[i]=0;
		b[i]=0;
	    for(j=1;j<N;j++){
			b[i]+=(double)1/(i+j-1);
		    A[i][j]=(double)1/(i+j-1);
		}
		r[i]=b[i];
		p[i]=r[i];
	}
	cout<<"当矩阵为"<<N-1<<"阶时："<<endl;
	for(k=1;maxd>=e;k++){             //当两个相邻Ｘ值的最大元素差值小于要求的误差值时，继续运算
		 for(i=1;i<N;i++){
			 Ap[i]=0;
			 for(j=1;j<N;j++)
			  Ap[i]+=A[i][j]*p[j];
		}
		m=inner(r,r);
		t=m/inner(p,Ap);
		for(i=1;i<N;i++){
		    x[i]=x[i]+t*p[i];
		    r[i]=r[i]-t*Ap[i];
		}
		s=inner(r,r)/m;
		for(i=1;i<N;i++)
		   p[i]=r[i]+s*p[i];
        maxd=t*maxe(p);             //maxd为两个相邻Ｘ值的最大元素差值
	 }
     cout<<endl<<"当要求精度为"<<e<<"时,"<<"需运算到第"<<k<<"次,"<<"此时误差maxe为"<<maxd<<endl;
	 cout<<"求出解为　X等于:"<<endl;
		cout<<"[  ";
		for(i=1;i<N;i++)
		   cout<<x[i]<<"  ";
        cout<<"]    的转置"<<endl<<endl;
}