《雅可比迭代法与矩阵地特征值.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,
3、int row,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=tem
4、p; bool 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
5、row,int 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); dele
6、te q; for(int i=0;in;i+) delete pi; delete p; 1. 用列主元消去法解方程例2 解方程组计算结果如下 矩阵直接三角分解法算法：将方程组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计算注：由于计算u的公式于计算y的公式形式上一样，故可直接对增广矩阵Ab=施行算法，此时U的第n+1
7、列元素即为y。程序与实例例3 求解方程组Ax=bA=, b=结果为#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,int row) for(int i=0;irow;i+) delete pi; delete p; void inputMatrix(double* p,int row,int col) f
8、or(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; for(int i=0;in;i+) coutXi+1=qiendl; cout=endl;void initMatrix(double* p,int row,int col) for(int i=
9、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的第一行给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
10、+) 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 backY(double* L,double* b,double* y,int row) for(int i=0;irow;i+) yi=0; /初始化y向量全为零 for(int i=0;i
11、row;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(int i=0;i=0;i-) xi=(yi-sumRowX(U,x,i,row)/Uii;int main() coutn; cout=endl; cout开始录入方阵数据.endl; doub
12、le* 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* L=newMatrix(n,n); /开辟方阵L initMatrix(L,n,n); /初始化L全为零 double* U=newMatrix(n,n); /开辟方阵U initMatrix(
13、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); delMatrix(U,n); delete b; delete y; delete x; 迭代法雅可比迭代法算法：设方程组Ax=b系数矩阵的对角线元素,M为迭代次数容许的最大值,为容许误差。
14、取初始向量x=，令k=0。对i=1,2,n 计算如果，如此输出，完毕；否如此执行。如果kM,如此不收敛，终止程序；否如此,转。程序与实例例4 用雅可比迭代法解方程组结果为迭代次数为20#include #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,int row) for(int i=0;irow;i+
15、) 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) for(int i=0;in;i+) coutXi+1=qit; coutendl;
16、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,double* 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,dou
17、ble 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,n); /开辟矩阵A inputMatrix(A,n,n); /录入数据到A中 cout录入方阵数据完毕.endl; cout=endlendl; cout开始录入列阵.endl; cout输入列阵:; double* b=new doublen; /开辟向
18、量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=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) co
19、ut不收敛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=0。对i=1,2,n计算如果，如此输出,完毕；否如此执行。如果如此不收敛，终止程序；否如此,转。程序与实例例5 用高斯-塞尔德迭代法解方程组结果为#include #include #include #include #define N 100main()i
20、nt 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;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;dof
21、or (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 main (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
22、); 用高斯-塞尔德迭代法解方程组迭代84次得解，假设用雅克比迭代法如此发散。#include#includevoid LOOP(float a1010,float b10,float x10,int);void main(void) float a1010,b10,x10,A10;double 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);prin
23、tf(请输入各方程值:);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-;for(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)
24、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,如此对任意非零初始向量用幂法求矩阵的按模最大的特征值和相应的特征向量。准确至6位有效数字。取。结果为A按模取最大的特征值为7.00000；相应的特征向量为。#include#includevoid LOOP(float a55,float u5,int);float MAX(float u5,int);void main(void)float a55,u5,
25、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=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);
