解線性方程組--數(shù)值分析

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 線性方程組的階數(shù)暫定為4階!"<<endl;cout<<"tt 源程序中可以更改方程階數(shù)!"<<endl;cout<<"tt 溫馨提示：更改MAXN!"<<

3、;endl;cout<<"t*"<<endl<<endl;cout<<"請先輸入有關(guān)該線性方程組的數(shù)據(jù):"<<endl<<endl; cout<<"線性方程組的系數(shù)矩陣:"<<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;/第一步：對線性方程組的系數(shù)矩陣進(jìn)行賦

6、值/右端項 以及初始向量的初始化double x0MAXN=0,0,0,0;double xMAXN;double eMAXN;/第二步：利用while循環(huán)進(jìn)行求解do for( i=0;i<MAXN;i+) /中間插入的 for循環(huán) 該循環(huán)用來 模仿求和公式for(j=0;j<MAXN && j!=i;j+) s+=aij*x0j; /外層循環(huán) 用來 求新的解 xi=(bi-s)/aii;/求誤差向量ei=fabs(xi-x0i); /調(diào)用 子函數(shù)取得 該誤差向量中的 最大分量M(e); /進(jìn)行判斷 將符合條件的 解向量輸出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+)/當(dāng)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 循環(huán)輸出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 /外層循環(huán)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算法 實現(xiàn)的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<