《雅可比迭代法与矩阵的特征值.doc》由会员分享，可在线阅读，更多相关《雅可比迭代法与矩阵的特征值.doc（21页珍藏版）》请在课桌文档上搜索。
1、实验五矩阵的LU分解法，雅可比迭代一、目的与要求： 熟悉求解线性方程组的有关理论和方法； 会编制列主元消去法、LU 分解法、雅可比及高斯塞德尔迭代法德程序； 通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法。二、实验容： 会编制列主元消去法、LU 分解法、雅可比及高斯塞德尔迭代法德程序，进一步了解各种方法的优缺点。三、程序与实例 列主元高斯消去法算法：将方程用增广矩阵Ab=(表示1) 消元过程对k=1,2,n-1选主元，找使得=如果，则矩阵A奇异，程序结束；否则执行。如果，则交换第k行与第行对应元素位置， j=k,n+1消元，对i=k+1, ,n计算对j=l+1, ,n+1计算2)
2、 回代过程若，则矩阵A奇异，程序结束；否则执行。；对i=n-1, ,2,1,计算程序与实例程序设计如下：#include #include using namespace std;void disp(double* p,int row,int col) for(int i=0;irow;i+) for(int j=0;jcol;j+) coutpij ; coutendl; void disp(double* q,int n) cout=endl; for(int i=0;in;i+) coutXi+1=qiendl; cout=endl;void input(double* p,int ro
3、w,int col) for(int i=0;irow;i+) cout输入第i+1行:; for(int j=0;jpij; int findMax(double* p,int start,int end) int max=start; for(int i=start;iabs(pmaxstart) max=i; return max;void swapRow(double* p,int one,int other,int col) double temp=0; for(int i=0;icol;i+) temp=ponei; ponei=potheri; potheri=temp; boo
4、l dispel(double* p,int row,int col) for(int i=0;irow;i+) int flag=findMax(p,i,row); /找列主元行号 if(pflagi=0) return false; swapRow(p,i,flag,col); /交换行 for(int j=i+1;jrow;j+) double elem=pji/pii; /消元因子 pji=0; for(int k=i+1;kcol;k+) pjk-=(elem*pik); return true;double sumRow(double* p,double* q,int row,in
5、t col) double sum=0; for(int i=0;i=0;i-) qi=(picol-1-sumRow(p,q,i,col)/pii; int main() coutn; double *p=new double* n; for(int i=0;in;i+) pi=new double n+1; input(p,n,n+1); if(!dispel(p,n,n+1) cout奇异endl; return 0;double* q=new doublen; for(int i=0;in;i+) qi=0; back(p,n,n+1,q); disp(q,n); delete q;
6、for(int i=0;in;i+) delete pi; delete p; 1. 用列主元消去法解方程例2 解方程组计算结果如下B=-1.461954C= 1.458125D=-6.004824E=-2.209018F= 14.719421 矩阵直接三角分解法算法：将方程组Ax=b 中的A分解为A=LU，其中L为单位下三角矩阵，U为上三角矩阵，则方程组Ax=b化为解2个方程组Ly=b，Ux=y，具体算法如下：对j=1,2,3,n计算对i=2,3,n计算对k=1,2,3,n:a. 对j=k,k+1,n计算b. 对i=k+1,k+2,n计算，对k=2,3,n计算,对k=n-1,n-2,2,1计
7、算注：由于计算u的公式于计算y的公式形式上一样，故可直接对增广矩阵Ab=施行算法，此时U的第n+1列元素即为y。程序与实例例3 求解方程组Ax=bA=, b=结果为X0= 3.000001X1=-2.000001X2= 1.000000X3= 5.000000#include using namespace std;double* newMatrix(int row,int col) double *p=new double* row; /行 for(int i=0;irow;i+) /列 pi=new double col; return p;void delMatrix(double* p
8、,int row) for(int i=0;irow;i+) delete pi; delete p; void inputMatrix(double* p,int row,int col) for(int i=0;irow;i+) cout输入第i+1行:; for(int j=0;jpij; void dispMatrix(double* p,int row,int col) for(int i=0;irow;i+) for(int j=0;jcol;j+) coutpij ; coutendl; void dispVector(double* q,int n) cout=endl; fo
9、r(int i=0;in;i+) coutXi+1=qiendl; cout=endl;void initMatrix(double* p,int row,int col) for(int i=0;irow;i+) for(int j=0;j=col?col:row); double sum=0; for(int i=0;imin;i+) sum+=(Uicol*Lrowi); return sum;void resolve(double* A,double* L,double* U,int row,int col) for(int i=0;icol;i+) U0i=A0i; /把A的第一行给
10、U的第一行 L00=A00; for(int i=1;irow;i+) /填充L的第一列 Li0=Ai0/A00; for(int i=1;irow;i+) for(int j=1;jcol;j+) if(i=j) Uij=Aij-sum(L,U,i,j); else Lij=(Aij-sum(L,U,i,j)/Ujj; for(int i=0;irow;i+) Lii=1;double sumRowY(double* L,double* y,int row) double sum=0; for(int i=0;irow;i+) sum+=(Lrowi*yi); return sum;void
11、 backY(double* L,double* b,double* y,int row) for(int i=0;irow;i+) yi=0; /初始化y向量全为零 for(int i=0;irow;i+) yi=bi-sumRowY(L,y,i);double sumRowX(double* U,double* x,int row,int col) double sum=0; for(int i=row+1;icol;i+) sum+=(Urowi*xi); return sum;void backX(double* U,double* y,double* x,int row) for(i
12、nt i=0;i=0;i-) xi=(yi-sumRowX(U,x,i,row)/Uii;int main() coutn; cout=endl; cout开始录入方阵数据.endl; double* A=newMatrix(n,n); /开辟矩阵A inputMatrix(A,n,n); /录入数据到A中 cout录入方阵数据完毕.endl; cout=endlendl; cout开始录入列阵.endl; cout输入列阵:; double* b=new doublen; /开辟向量b for(int i=0;ibi; /录入数据到b中 cout录入列阵数据完毕.endl; double*
13、L=newMatrix(n,n); /开辟方阵L initMatrix(L,n,n); /初始化L全为零 double* U=newMatrix(n,n); /开辟方阵U initMatrix(U,n,n); /初始化U resolve(A,L,U,n,n); /分解A为L与U double* y=new doublen; /开辟向量y backY(L,b,y,n); /回代求出y double* x=new doublen; /开辟向量X backX(U,y,x,n); /回代求出X dispVector(x,n); /释放空间 delMatrix(A,n); delMatrix(L,n);
14、 delMatrix(U,n); delete b; delete y; delete x; 迭代法雅可比迭代法算法：设方程组Ax=b系数矩阵的对角线元素,M为迭代次数容许的最大值,为容许误差。取初始向量x=，令k=0。对i=1,2,n 计算如果，则输出，结束；否则执行。如果kM,则不收敛，终止程序；否则,转。程序与实例例4 用雅可比迭代法解方程组结果为迭代次数为20X0= 1.000000X1= 2.000000X2=-1.000000#include #include using namespace std;double* newMatrix(int row,int col) double
15、 *p=new double* row; /行 for(int i=0;irow;i+) /列 pi=new double col; return p;void delMatrix(double* p,int row) for(int i=0;irow;i+) delete pi; delete p; void inputMatrix(double* p,int row,int col) for(int i=0;irow;i+) cout输入第i+1行:; for(int j=0;jpij; void dispMatrix(double* p,int row,int col) for(int
16、i=0;irow;i+) for(int j=0;jcol;j+) coutpij ; coutendl; void dispVector(double* q,int n) for(int i=0;in;i+) coutXi+1=qit; coutendl;double sumRow(double* A,double* x1,int row,int col) double sum=0; for(int i=0;icol;i+) if(row=i) continue; sum+=(Arowi*x1i); return sum;void iterat(double* A,double* b,dou
17、ble* x1,double* x2,int row) for(int i=0;irow;i+) x2i=(bi-sumRow(A,x1,i,row)/Aii;bool blow_error(double* x1,double* x2,int row,double e) double sum=0; for(int i=0;irow;i+) sum+=abs(x2i-x1i); if(sume) return true; else return false;int main() coutn; cout=endl; cout开始录入方阵数据.endl; double* A=newMatrix(n,
18、n); /开辟矩阵A inputMatrix(A,n,n); /录入数据到A中 cout录入方阵数据完毕.endl; cout=endlendl; cout开始录入列阵.endl; cout输入列阵:; double* b=new doublen; /开辟向量b for(int i=0;ibi; /录入数据到b中 cout录入列阵数据完毕.endl; cout=endl; coutM; coute; cout=endl; double* x1=new doublen; /开辟x1向量 double* x2=new doublen; /开辟x2向量 cout输入初始向量:; for(int i=
19、0;ix1i; cout=endl; int k=0; /迭代计数器 for(;) iterat(A,b,x1,x2,n); /迭代一次 if(blow_error(x1,x2,n,e) dispVector(x2,n); cout迭代次数为:k=M) cout不收敛endl; return 0; else k+; for(int i=0;in;i+) x1i=x2i; /释放空间 delMatrix(A,n); delete b; delete x1; delete x2; 高斯-塞尔德迭代法算法:设方程组Ax=b的系数矩阵的对角线元素,M为迭代次数容许的最大值,为容许误差取初始向量,令k=
20、0。对i=1,2,n计算如果，则输出,结束；否则执行。如果则不收敛，终止程序；否则,转。程序与实例例5 用高斯-塞尔德迭代法解方程组结果为X0=3.000000X1=2.000000X2=1.000000#include #include #include #include #define N 100main()int i;float *x;float c12=8.0,-3.0,2.0,20.0,4.0,11.0,-1.0,33.0,6.0,3.0,12.0,36.0;float *GauseSeidel(float *,int);x=GauseSeidel(c,3);for (i=0;i=2
21、;i+) printf(x%d=%fn,i,xi);getch();float *GauseSeidel(float *a,int n)int i,j,nu=0;float *x,dx,d;x=(float *)malloc(n*sizeof(float);for (i=0;i=n-1;i+) xi=0.0;dofor (i=0;i=n-1;i+)d=0.0;for (j=0;j=N)printf(fold numbern);nu+;while (fabs(dx)1e-6);return x;例6 用雅可比迭代法解方程组迭代4次得解,若用高斯-塞尔德迭代法则发散。#includevoid ma
22、in (void)int k,n;double x3=7,2,5;for(k=0;k5;k+)double a,b;a=x0;b=x1;x0=(7-2.0*x1+2*x2)/1;x1=(2-a-x2)/1;x2=(5-2*a-2*b)/1;for(n=0;n3;n+)printf(x%d=%8.6fn,n,xn); 用高斯-塞尔德迭代法解方程组迭代84次得解，若用雅克比迭代法则发散。#include#includevoid LOOP(float a1010,float b10,float x10,int);void main(void) float a1010,b10,x10,A10;doub
23、le S;int M,n,i,j;printf(请输入方阵阶数:);scanf(%d,&n);printf(请输入最大允许迭代次数:);scanf(%d,&M);printf(请按行输入各方程系数:);for(i=0;i=n-1;i+)for(j=0;j=n-1;j+)scanf(%f,&aij);printf(请输入各方程值:);for(i=0;i=n-1;i+)scanf(%f,&bi);printf(请依次输入首次迭代x值:);for(i=0;i=n-1;i+)scanf(%f,&xi);doS=0.0;for(i=0;i=n-1;i+)Ai=xi;LOOP(a,b,x,n);M-;fo
24、r(i=0;i=0&S=0.000001);if(M=0)printf(迭代次数M=%dn,M);for(i=0;i=n-1;i+)printf(x%d=%fn,i,xi);elseprintf(该迭代发散n);void LOOP(float a1010,float b10,float x10,int n)float S1,S2,A10;int i,j;for(i=0;i=n-1;i+)Ai=xi;for(i=0;i=n-1;i+)S1=0.0;S2=0.0;for(j=0;j=i-1;j+)S1=S1+aij*xj;for(j=i+1;j,则对任意非零初始向量用幂法求矩阵的按模最大的特征值和
25、相应的特征向量。精确至6位有效数字。取。结果为A按模取最大的特征值为7.00000；相应的特征向量为。#include#includevoid LOOP(float a55,float u5,int);float MAX(float u5,int);void main(void)float a55,u5,x5,y,z;int i,j,n; printf(请输入方阵阶数:);scanf(%d,&n);printf(请按行输入各矩阵元素值:);for(i=0;i=n-1;i+)for(j=0;j=n-1;j+)scanf(%f,&aij);printf(请输入初次迭代向量:);for(i=0;i=
26、n-1;i+)scanf(%f,&ui);y=MAX(u,n);doz=y;LOOP(a,u,n);y=MAX(u,n);for(i=0;i=0.0000005);printf(矩阵特征值=%fn,y);printf(矩阵特征向量x:n);for(i=0;i=n-1;i+)printf(%10fn,xi);void LOOP(float a55,float u5,int n)float S,U5;int i,j;for(i=0;i=n-1;i+)Ui=ui;for(i=0;i=n-1;i+)S=0.0;for(j=0;j=n-1;j+)S=S+aij*Uj;ui=S;float MAX(float u5,int n)float max;int i;max=u0;for(i=0;imax)max=ui;return(max);
