解线性方程组--数值分析

1、# include<iostream># include<cmath># define MAXN 4using namespace std;void M(double eMAXN);void Jacobi(double aMAXNMAXN, double bMAXN, double TOL);void Gauss( double aMAXNMAXN, double bMAXN, double TOL);void SOR( double aMAXNMAXN, double bMAXN, double w, double TOL);int main()int i,j,cho

2、ice;double aMAXNMAXN, bMAXN,TOL,w;TOL=pow(10,-6);cout<<"t*"<<endl; cout<<"tt本程序运用三种方法解简单的线性方程组!"<<endl;cout<<"tt 线性方程组的阶数暂定为4阶!"<<endl;cout<<"tt 源程序中可以更改方程阶数!"<<endl;cout<<"tt 温馨提示：更改MAXN!"<<

3、;endl;cout<<"t*"<<endl<<endl;cout<<"请先输入有关该线性方程组的数据:"<<endl<<endl; cout<<"线性方程组的系数矩阵:"<<endl; for(i=0;i<MAXN;i+)for(j=0;j<MAXN;j+)cin>>aij; cout<<"线性方程的右端项:"<<endl;for(i=0;i<MAXN;i+)cin

4、>>bi;cout<<endl; cout<<"以下有三种方法可供选择:"<<endl;cout<<"1-Jacobi 2-Gauss seidl 3-SOR "<<endl;cout<<"请选择要运用的方法:"<<endl;cin>>choice; switch(choice)case 1:Jacobi(a,b,TOL);break;case 2:Gauss(a,b,TOL);break;case 3:cout<<

5、"请输入w的值："cin>>w;SOR(a,b,w,TOL);break;default:cout<<"该线性方程组无法运用以上的三种方法求的解!"<<endl;break;return 0;void M(double eMAXN)for(int i=1;i<MAXN;i+)if(ei>e0)e0=ei;void Jacobi(double aMAXNMAXN,double bMAXN, double TOL)int i,j; /声明一个中间变量 double s=0;/第一步：对线性方程组的系数矩阵进行赋

6、值/右端项 以及初始向量的初始化double x0MAXN=0,0,0,0;double xMAXN;double eMAXN;/第二步：利用while循环进行求解do for( i=0;i<MAXN;i+) /中间插入的 for循环 该循环用来 模仿求和公式for(j=0;j<MAXN && j!=i;j+) s+=aij*x0j; /外层循环 用来 求新的解 xi=(bi-s)/aii;/求误差向量ei=fabs(xi-x0i); /调用 子函数取得 该误差向量中的 最大分量M(e); /进行判断 将符合条件的 解向量输出if(e0<TOL)for( i=

7、0;i<MAXN;i+)cout<<"x"<<i<<"="<<xi<<endl; /如果不满足条件 则更新 初始向量for(i=0;i<MAXN;i+)x0i=xi; /清空 中介值 ss=0;while(e0>TOL);void Gauss(double aMAXNMAXN, double bMAXN, double TOL) int i,j;double s=0,s1=0,s2=0;double x0MAXN=0,0,0,0;double xMAXN;double eMAX

8、N;dofor( i=0;i<MAXN;i+)/当i等于0 的时候 值是利用Jacobi 公式求出来的 if(i=0)for( j=1; j<MAXN;j+)s+=aij*x0j; xi=(bi-s)/aii; else /利用Gauss seidel 方法计算 求和公式for(j=0; j<MAXN&& j!=i; j+) if(j<i) s1+=aij*xj; else s2+=aij*x0j; s=s1+s2; xi=(bi-s)/aii; /误差向量的 计算 ei=fabs(xi-x0i); /求向量的 范式M(e);/判断是否满足条件 将满足条

9、件的 解向量输出if(e0<TOL) /利用for 循环输出for(i=0;i<MAXN;i+)cout<<"x"<<i<<"="<<xi<<endl; /否则 就更换 初始解向量的 值for(i=0;i<MAXN;i+)x0i=xi; /更新解向量的 值 s1=s2=s=0;while(e0>TOL);void SOR(double aMAXNMAXN, double bMAXN, double w,double TOL)int i,j;/中介值double s=0,s

10、1=0,s2=0;double x0MAXN=0,0,0,0;double xMAXN;double eMAXN; do /外层循环for( i=0;i<4;i+) /解向量 第一个分量的 值if(i=0)for(j=1; j<MAXN;j+)s+=aij*x0j;xi=(1-w)*x0i+w*(bi-s)/aii;else/其余分量是利用 SOR算法 实现的for(j=0;j<MAXN && j!=i; j+) /用的 新值if(j<i)s1+=aij*xj;/用旧值elses2+=aij*x0j; /中介值的改变s=s1+s2;xi=(1-w)*x0i+w*(bi-s)/aii; /误差向量 ei=fabs(xi-x0i);/求向量的 范式M(e); /判断是否满足条件 if(e0<TOL) /若满足条件 则就爱那个解向量求出for(i=0;i<