雅可比迭代法与矩阵的特征值

雅可比迭代法与矩阵的特征值 实验五 矩阵的lu分解法，雅可比迭代法 班级: 学号: 姓名: 实验五 矩阵的LU分解法，雅可比迭代 一、目的与要求: ? 熟悉求解线性方程组的有关理论和方法; ? 会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法 德程序; ? 通过实际计算，进一步了解各种方法的优缺点，选择合适的 数值方法。 二、实验 内容 财务内部控制制度的内容财务内部控制制度的内容人员招聘与配置的内容项目成本控制的内容消防安全演练内容 : ? 会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法 德程序，进一步了解 各种方法的优缺点。 三、程序与实例 ? 列主元高斯消去法 算法:将方程用增广矩阵[A?b]=(aij)n?(n?1)表示 1) 消元过程 对k=1,2,…,n-1 ?选主元，找ik??k,k?1,?,n?使得 aik,k=max k?i?naik ?如果aik,k?0，则矩阵A奇异，程序结束;否则执行?。 ?如果ik?k，则交换第k行与第ik行对应元素位置， akj?aikj j=k,?,n+1 ?消元，对i=k+1, ?,n计算 lik?aik/akk 对j=l+1, ?,n+1计算 aij?aij?likakj 2) 回代过程 ?若ann?0，则矩阵A奇异，程序结束;否则执行?。 ?xn?an,n?1/ann;对i=n-1, ?,2,1,计算 ?xi???ai,n?1? ? 程序与实例 程序 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 如下: ??ax?ijj?/aii j?i?1?n #include <iostream> #include <cmath> using namespace std; void disp(double** p,int row,int col){ for(int i=0;i<row;i++){ for(int j=0;j<col;j++) cout<<p[i][j]<<' '; cout<<endl; } } void disp(double* q,int n){ cout<<"====================================="<<endl; for(int i=0;i<n;i++) cout<<"X["<<i+1<<"]="<<q[i]<<endl; cout<<"====================================="<<endl; } void input(double** p,int row,int col){ for(int i=0;i<row;i++){ cout<<"输入第"<<i+1<<"行:"; for(int j=0;j<col;j++) cin>>p[i][j]; } } int findMax(double** p,int start,int end){ int max=start; for(int i=start;i<end;i++){ if(abs(p[i][start])>abs(p[max][start])) max=i; } return max; } void swapRow(double** p,int one,int other,int col){ double temp=0; for(int i=0;i<col;i++){ temp=p[one][i]; p[one][i]=p[other][i]; p[other][i]=temp; } } bool dispel(double** p,int row,int col){ for(int i=0;i<row;i++){ int flag=findMax(p,i,row); //找列主元行号 if(p[flag][i]==0) return false; swapRow(p,i,flag,col); //交换行 for(int j=i+1;j<row;j++){ double elem=p[j][i]/p[i][i]; //消元因子 p[j][i]=0; for(int k=i+1;k<col;k++){ p[j][k]-=(elem*p[i][k]); } } } return true; } double sumRow(double** p,double* q,int row,int col){ double sum=0; for(int i=0;i<col-1;i++){ if(i==row) continue; sum+=(q[i]*p[row][i]); } return sum; } void back(double** p,int row,int col,double* q){ for(int i=row-1;i>=0;i--){ q[i]=(p[i][col-1]-sumRow(p,q,i,col))/p[i][i]; } } int main() { cout<<"Input n:"; int n; //方阵的大小 cin>>n; double **p=new double* [n]; for(int i=0;i<n;i++){ p[i]=new double [n+1]; } input(p,n,n+1); if(!dispel(p,n,n+1)){ cout<<"奇异"<<endl; return 0; } double* q=new double[n]; for(int i=0;i<n;i++) q[i]=0; back(p,n,n+1,q); disp(q,n); delete[] q; for(int i=0;i<n;i++) delete[] p[i]; delete[] p; } 1. 用列主元消去法解方程 ?0.101x1 ?2.304x2?3.555x3?1.183???1.347x1?3712x2?4.623x3?2.137??2.835x?1.072x?5.643x?3.035123? 例2 解方程组 ?8.77B?2.40C?5.66D?1.55E?1.0F?-32.04?4.93B?1.21C?4.48D?1.10E?1.0F?-20.07 ???3.53B?1.46C?2.92D?1.21E?1.0F?-8.53?5.05B?4.04C?2.51D?2.01E?1.0F?-6.30???3.54B?1.04C?3.47D?1.02E?1.0F?-12.04计算结果如下 B=-1.461954 C= 1.458125 D=-6.004824 E=-2.209018 F= 14.719421 ? 矩阵直接三角分解法 算法:将方程组Ax=b 中的A分解为A=LU，其中L为单位下三角矩阵， U为上三角矩阵，则方程组Ax=b化为解2个方程组Ly=b，Ux=y，具体算 法如下: ?对j=1,2,3,…,n计算 u1j?a1j 对i=2,3,…,n计算 li1?ai1/a11 ?对k=1,2,3,…,n: a. 对j=k,k+1,…,n计算 ukj?akj??lkquqj q?1k?1 b. 对i=k+1,k+2,…,n计算 lik?(aik??liquqk)/ukk q?1k?1 ?y1?b1，对k=2,3,…,n计算 yk?bk??lkqyq q?1k?1 ?xn?yn/unn,对k=n-1,n-2,…,2,1计算 xk?(yk? q?k?1?unkqxq)/ukk 注:由于计算u的公式于计算y的公式形式上一样，故可直接对增广矩阵 ?a11a12?aa2221[A? b]=?????an1an2a1,n?1??a2na2,n?1?? ???????annan,n?1??a1n 施行算法?，?，此时U的第n+1列元素即为y。 程序与实例 例3 求解方程组Ax=b ?12?128??27??547?2??4??, b=?? A=???3795??11?????6?12?83???49? 结果为 X[0]= 3.000001 X[1]=-2.000001 X[2]= 1.000000 X[3]= 5.000000 #include <iostream> using namespace std; double** newMatrix(int row,int col){ double **p=new double* [row]; //行 for(int i=0;i<row;i++) //列 p[i]=new double [col]; return p; } void delMatrix(double** p,int row){ for(int i=0;i<row;i++) delete[] p[i]; delete[] p; } void inputMatrix(double** p,int row,int col){ for(int i=0;i<row;i++){ cout<<"输入第"<<i+1<<"行:"; for(int j=0;j<col;j++) cin>>p[i][j]; } } void dispMatrix(double** p,int row,int col){ for(int i=0;i<row;i++){ for(int j=0;j<col;j++) cout<<p[i][j]<<' '; cout<<endl; } } void dispVector(double* q,int n){ cout<<"====================================="<<endl; for(int i=0;i<n;i++) cout<<"X["<<i+1<<"]="<<q[i]<<endl; cout<<"====================================="<<endl; } void initMatrix(double** p,int row,int col){ for(int i=0;i<row;i++) for(int j=0;j<col;j++) p[i][j]=0; } double sum(double** L,double** U,int row,int col){ int min=(row>=col?col:row); double sum=0; for(int i=0;i<min;i++) sum+=(U[i][col]*L[row][i]); return sum; } void resolve(double** A,double** L,double** U,int row,int col){ for(int i=0;i<col;i++) U[0][i]=A[0][i]; //把A的第一行给U的第一行 L[0][0]=A[0][0]; for(int i=1;i<row;i++) //填充L的第一列 L[i][0]=A[i][0]/A[0][0]; for(int i=1;i<row;i++) for(int j=1;j<col;j++){ if(i<=j) U[i][j]=A[i][j]-sum(L,U,i,j); else L[i][j]=(A[i][j]-sum(L,U,i,j))/U[j][j]; } for(int i=0;i<row;i++) L[i][i]=1; } double sumRowY(double** L,double* y,int row){ double sum=0; for(int i=0;i<row;i++) sum+=(L[row][i]*y[i]); return sum; } void backY(double** L,double* b,double* y,int row){ for(int i=0;i<row;i++) y[i]=0; //初始化y向量全为零 for(int i=0;i<row;i++) y[i]=b[i]-sumRowY(L,y,i); } double sumRowX(double** U,double* x,int row,int col){ double sum=0; for(int i=row+1;i<col;i++) sum+=(U[row][i]*x[i]); return sum; } void backX(double** U,double* y,double* x,int row){ for(int i=0;i<row;i++) x[i]=0; //初始化x向量全为零 for(int i=row-1;i>=0;i--) x[i]=(y[i]-sumRowX(U,x,i,row))/U[i][i]; } int main() { cout<<"输入方阵大小 n :"; int n=0; cin>>n; cout<<"====================================="<<endl; cout<<"开始录入方阵数据....."<<endl; double** A=newMatrix(n,n); //开辟矩阵A inputMatrix(A,n,n); //录入数据到A中 cout<<"录 入方阵数据完毕......"<<endl; cout<<"====================================="<<endl<< endl; cout<<"开始录入列阵....."<<endl; cout<<"输入列阵:"; double* b=new double[n]; //开辟向量b for(int i=0;i<n;i++) cin>>b[i]; //录入数据到b中 cout<<"录入列阵数据完毕....."<<endl; double** L=newMatrix(n,n); //开辟方阵L initMatrix(L,n,n); // double** U=newMatrix(n,n); // initMatrix(U,n,n); // resolve(A,L,U,n,n); // double* y=new double[n]; // backY(L,b,y,n); // double* x=new double[n]; // backX(U,y,x,n); // dispVector(x,n); //释放空间 delMatrix(A,n); delMatrix(L,n); delMatrix(U,n); delete[] b; delete[] y; delete[] x; } 初始化L全为零 开辟方阵U 初始化U 分解A为L与U 开辟向量y 回代求出y 开辟向量X 回代求出X ? 迭代法 雅可比迭代法 算法:设方程组Ax=b系数矩阵的对角线元素aii?0(i?1,2,?,n),M为迭代次数容许的最大值,ε为容许误差。 ?取初始向量x=(x1,x2,?,xn)，令k=0。 ?对i=1,2,…,n 计算 n1?(bi??aijx(k) j) aiij?1 j?i(0)(0)(0)Tx(k?1)i ?如果?x i?1n(k?1)i?xi(k)??，则输出x(k?1)，结束;否则执行?。 ?如果k?M,则不收敛，终止程序;否则k?k?1,转?。 程序与实例 例4 用雅可比迭代法解方程组 ?5x1?2x2?x3?8??2x1?8x2?3x3?21 ?x?3x?6x?123?1 结果为 迭代次数为20 X[0]= 1.000000 X[1]= 2.000000 X[2]=-1.000000 #include <iostream> #include <cmath> using namespace std; double** newMatrix(int row,int col){ double **p=new double* [row]; //行 for(int i=0;i<row;i++) //列 p[i]=new double [col]; return p; } void delMatrix(double** p,int row){ for(int i=0;i<row;i++) delete[] p[i]; delete[] p; } void inputMatrix(double** p,int row,int col){ for(int i=0;i<row;i++){ cout<<"输入第"<<i+1<<"行:"; for(int j=0;j<col;j++) cin>>p[i][j]; } } void dispMatrix(double** p,int row,int col){ for(int i=0;i<row;i++){ for(int j=0;j<col;j++) cout<<p[i][j]<<' '; cout<<endl; } } void dispVector(double* q,int n){ for(int i=0;i<n;i++) cout<<"X["<<i+1<<"]="<<q[i]<<'\t'; cout<<endl; } double sumRow(double** A,double* x1,int row,int col){ double sum=0; for(int i=0;i<col;i++){ if(row==i) continue; sum+=(A[row][i]*x1[i]); } return sum; } void iterat(double** A,double* b,double* x1,double* x2,int row){ for(int i=0;i<row;i++) x2[i]=(b[i]-sumRow(A,x1,i,row))/A[i][i]; } bool blow_error(double* x1,double* x2,int row,double e){ double sum=0; for(int i=0;i<row;i++){ sum+=abs(x2[i]-x1[i]); } if(sum<e) return true; else return false; } int main() { cout<<"输入方阵大小 n :"; int n=0; cin>>n; cout<<"====================================="<<endl; cout<<"开始录入方阵数据....."<<endl; double** A=newMatrix(n,n); //开辟矩阵A inputMatrix(A,n,n); //录入数据到A中 cout<<"录入方阵数据完毕......"<<endl; cout<<"====================================="<<endl<< endl; cout<<"开始录入列阵....."<<endl; cout<<"输入列阵:"; double* b=new double[n]; //开辟向量b for(int i=0;i<n;i++) cin>>b[i]; //录入数据到b中 cout<<"录入列阵数据完毕....."<<endl; cout<<"====================================="<<endl; cout<<"输入迭代次数容许的最大值:"; int M=0; cin>>M; cout<<"输入容许误差:"; double e=0; cin>>e; cout<<"====================================="<<endl; double* x1=new double[n]; //开辟x1向量 double* x2=new double[n]; //开辟x2向量 cout<<"输入初始向量:"; for(int i=0;i<n;i++) cin>>x1[i]; cout<<"====================================="<<endl; int k=0; //迭代计数器 for(;;){ iterat(A,b,x1,x2,n); //迭代一次 if(blow_error(x1,x2,n,e)){ dispVector(x2,n); cout<<"迭代次数为:"<<k<<endl; return 0; }else{ if(k>=M){ cout<<"不收敛"<<endl; return 0; }else{ k++; for(int i=0;i<n;i++) x1[i]=x2[i]; } } } //释放空间 delMatrix(A,n); delete[] b; delete[] x1; delete[] x2; } ? 高斯-塞尔德迭代法 算法:设方程组Ax=b的系数矩阵的对角线元素,aii?0(i?1,2,?,n),M为迭代次数容许的最大值,ε为容许误差 ?取初始向量x?(x1,x2,?,xn),令k=0。 ?对i=1,2,…,n计算 (0)(0)(0)T x(k?1) ii?11?(bi??aijx( jk?1)?aiij?1 nj?i?1?anijx(jk)) ?如果?x i?1(k?1)i?xi(k)??，则输出x(k?1),结束;否则执行?。 ?如果k?M,则不收敛，终止程序;否则k?k?1,转?。 程序与实例 例5 用高斯-塞尔德迭代法解方程组 ?8x1?3x2?2x3?20??4x1?11x2?x3?33 ?6x?3x?12x?3623?1 结果为 X[0]=3.000000 X[1]=2.000000 X[2]=1.000000 #include <stdio.h> #include <conio.h> #include <malloc.h> #include <math.h> #define N 100 main() { int i; float *x; float c[12]={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]=%f\n",i,x[i]); 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++) x[i]=0.0; do { for (i=0;i<=n-1;i++) { d=0.0; for (j=0;j<=n-1;j++) d+=*(a+i*(n+1)+j)*x[j]; dx=(*(a+i*(n+1)+n)-d)/(*(a+i*(n+1)+i)); x[i]+=dx; } if (nu>=N) { printf("fold number\n"); } nu++; } while (fabs(dx)>1e-6); return x; } 例6 用雅可比迭代法解方程组 ?x1?2x2?2x3?7??x1?x2?x3?2 ?2x?2x?x?523?1 迭代4次得解(1,2,?1),若用高斯-塞尔德迭代法则发散。 #include<stdio.h> void main (void) { int k,n; double x[3]={7,2,5}; for(k=0;k<5;k++) { double a,b; a=x[0]; b=x[1]; x[0]=(7-2.0*x[1]+2*x[2])/1; x[1]=(2-a-x[2])/1; x[2]=(5-2*a-2*b)/1; } for(n=0;n<3;n++) { printf("x[%d]=%8.6f\n",n,x[n]); } } T 用高斯-塞尔德迭代法解方程组 ?x1?0.9x2?0.9x3?1.9??0.9x1?x2?0.9x3?2.0 ?0.9x?0.9x?x?1.7123? 迭代84次得解?1,2,3?，若用雅克比迭代法则发散。 T #include<stdio.h> #include<math.h> void LOOP(float a[10][10],float b[10],float x[10],int); void main(void) { float a[10][10],b[10],x[10],A[10]; 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",&a[i][j]); } } printf("请输入各方程值:"); for(i=0;i<=n-1;i++) { scanf("%f",&b[i]); } printf("请依次输入首次迭代x值:"); for(i=0;i<=n-1;i++) { scanf("%f",&x[i]); } do { S=0.0; for(i=0;i<=n-1;i++) { A[i]=x[i]; } LOOP(a,b,x,n); M--; for(i=0;i<=n-1;i++) { S=S+fabs(x[i]-A[i]); } }while(M>=0&&S>=0.000001); if(M>=0) { printf("迭代次数M=%d\n",M); for(i=0;i<=n-1;i++) { printf("x[%d]=%f\n",i,x[i]); } } else { printf("该迭代发散\n"); } } void LOOP(float a[10][10],float b[10],float x[10],int n) { float S1,S2,A[10]; int i,j; for(i=0;i<=n-1;i++) { A[i]=x[i]; } for(i=0;i<=n-1;i++) { S1=0.0; S2=0.0; for(j=0;j<=i-1;j++) { S1=S1+a[i][j]*x[j]; } for(j=i+1;j<=n-1;j++) { S2=S2+a[i][j]*A[j]; } x[i]=(b[i]-S1-S2)/a[i][i]; } } 实验六 矩阵的特征值与特征向量的计算 班级: 学号: 姓名: 实验六 矩阵的特征值与特征向量的计算 一、目的与要求: 领会求矩阵特征值及特征向量的幂法的理论及其方法; 会编制幂法的计算程序，并用来计算有关问题。 二、实验内容: 编制幂法的计算程序，并用来计算有关问题 三、算法概要 幂法是矩阵主特征值的一种迭代方法。设A?Rn?n有n个线性无关的特征向量x1，x2,?,xn,而相应的特征值满足?1>?2????n,则对任意非零初始向量 ?32?用幂法求矩阵A???的按模最大的特征值和相应的特征向量。精确至6位有效数字。45?? 取u0???。 ?1? ?1? 结果为A按模取最大的特征值为7.00000;相应的特征向量为? #include<stdio.h> #include<math.h> void LOOP(float a[5][5],float u[5],int); float MAX(float u[5],int); void main(void) { float a[5][5],u[5],x[5],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",&a[i][j]); } } printf("请输入初次迭代向量:"); for(i=0;i<=n-1;i++) { scanf("%f",&u[i]); } y=MAX(u,n); ?0.500000??。 1.000000?? do { z=y; LOOP(a,u,n); y=MAX(u,n); for(i=0;i<=n-1;i++) { x[i]=u[i]/y; u[i]=x[i]; } }while(fabs(z-y)>=0.0000005); printf("矩阵特征值λ=%f\n",y); printf("矩阵特征向量x:\n"); for(i=0;i<=n-1;i++) { printf("%10f\n",x[i]); } } void LOOP(float a[5][5],float u[5],int n) { float S,U[5]; int i,j; for(i=0;i<=n-1;i++) { U[i]=u[i]; } for(i=0;i<=n-1;i++) { S=0.0; for(j=0;j<=n-1;j++) { S=S+a[i][j]*U[j]; } u[i]=S; } } float MAX(float u[5],int n) { float max; int i; max=u[0]; for(i=0;i<=n-1;i++) { if(u[i]>max) } { max=u[i]; } } return(max);