计算实习报告一 实习目的（1）了解矩阵特征值与相应特征向量求解的意义，理解幂法和反幂法的原理，能编制此算法的程序，并能求解实际问题。

（2）通过对比非线性方程的迭代法，理解线性方程组迭代解法的原理，学会编写Jacobi 迭代法程序，并能求解中小型非线性方程组。

初始点对收敛性质及收敛速度的影响。

（3）理解 QR 法计算矩阵特征值与特征向量的原理，能编制此算法的程序，并用于实际问题的求解。

二 问题定义及题目分析1. 分别用幂法和幂法加速技术求矩阵2.5 2.53.00.50.0 5.0 2.0 2.00.50.5 4.0 2.52.52.55.03.5-⎛⎫⎪-⎪= ⎪--⎪--⎝⎭A 的主特征值和特征向量．2. 对于实对称矩阵n n ⨯∈A R ，用Jacobi 方法编写其程序，并用所编程序求下列矩阵的全部特征值．151541141144114114⨯-⎛⎫ ⎪-- ⎪ ⎪- ⎪=⎪ ⎪- ⎪-- ⎪ ⎪-⎝⎭A3. 对于实矩阵n n ⨯∈A R ，用QR 方法编写其程序，并用所编程序求下列矩阵的全部特征值:11121113,4,5,62311111n n n nn n ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥==+⎢⎥⎢⎥⎢⎥⎢⎥+⎣⎦A三 概要设计（1） 幂法用于求按模最大的特征值及其对应的特征向量的一种数值算法，它要求矩阵 A 的特征值有如下关系： 12n ...λλλ>≥≥ ，对于相应的特征向量。

其算法如下：Step 0：初始化数据0,, 1.A z k = Step 1：计算1k k y A z +=。

Step 2：令 kkm y ∞=。

Step 3：令 k k k z y m = ；如果1k k m m +≈或1k k z z +≈，则 gotoStep 4；否则 ， k = k + 1 ，goto Step 1。

Step 4：输出结果算法说明与要求输入参数为实数矩阵、初始向量、误差限与最大迭代次数。

输出参数为特征值及相对应的特征向量。

注意初始向量不能为“0”向量。

（2） 迭代法的原理如果能将方程 Ax =b 改写成等价形式：x=Bx+f 。

如果B 满足：ρ(B )<1，则对于任意初始向量 x (0) ，由迭代 x ( k + 1) = Bx (k ) + f 产生的序列均收敛到方程组的精确解。

迭代法中两种最有名的迭代法就是Jacobi 迭代法，它的迭代矩阵 B 为：1()J DL U -=-+,1f D b -=其中，D 为系数矩阵 A 的对角元所组成对角矩阵，L 为系数矩阵 A 的对角元下方所有元素所组成的下三角矩阵，U 为系数矩阵 A 的对角元上方所有元素所组成的上三角矩阵。

算法如下：Step 0：初始化数据 00,,,,k A b x δ=和ε。

Step 1：计算D,L,U,J 或G, 得到迭代矩阵B. Step 2：:1kk =+0x B x f *=+0x x =如果0x x δ-<或()fx ε≤，goto Step 3;否则 gotoStep ２。

Step 3：输出结果。

程序说明与要求程序的输入参数为系数矩阵与常量、初始向量及误差控制，输出参数为方程组的近似解。

要求输入系数 A 中对角元不能存在 0，如果对角元出现 0 元素，则可以通过交换方程组中方程的次序解决。

（3） QR 法原理QR 算法是求实对称矩阵全部特征的最有效方法。

其基本过程就是利用矩阵的QR 分解，即将任一实数矩阵分解为一个正交矩阵Q 与一个上三角矩阵 R 的乘积。

QR 算法的计算过程如下Step 0：初始化数据0000,,,,1A Q R A Q R k == Step 1：令 11k k k A Q R ++=。

Step 2：计算 k k k A Q R =Step 3：令k m 为k R 所有对角元下方元素绝对值之和；如果k m ε< ，则 goto Step 4；否则， k = k + 1 ，goto Step 1。

Step 4：输出结果。

算法的要求与说明先利用函数 house 将对称矩阵 A 通过Householder 相似变换为对称三对角阵 T ，然后通过程序 QR_method 求出矩阵 T 的特征值。

程序输入参数为矩阵 A ，输出参数为矩阵 A 的所有特征值。

四 详细设计源程序 （1）#include"math.h" #include"stdio.h"void main() /*主函数*/ {int i,j; int n=4,k=0; floata[4][4]={{2.5,-2.5,3.0,0.5},{0.0,5.0,-2.0,2.0},{-0.5,-0.5,4.0,2.5},{-2.5,-2.5,5.0,3.5}},x[4]={1,1,1,1};double max,m,y[4],z[4],ej=0.0001,e;printf("k max(xk) y(k)=x(k)/max(x(k))x(k+1)=a*y(k)\n"); /*输出格式*/ cyclo:{for(i=0;i<n;i++) /*循环次数统计并从0次开始输出，并缓存临时向量组用于后面的循环条件比较*/z[i]=x[i];k=k+1;printf("%d ",k-1);{max=0.0; /*比较求出向量组x[]的元素最大值*/for(i=0;i<n;i++)if(fabs(x[i])>max)max=fabs(x[i]);}for(i=0;i<n;i++) /*求出向量组y[]的各元素值*/y[i]=x[i]/max;printf("%f ",max);printf("("); /*输出y[i]*/for(i=0;i<n;i++)printf("%f ",y[i]);printf(")");for(i=0;i<n;i++) /*计算求出x(k+1)的元素值*/{m=0.0;for(j=0;j<n;j++)m+=a[i][j]*y[j];x[i]=m;}printf("("); /*输出x(k+1)的元素值*/for(i=0;i<n;i++)printf("%f ",x[i]);printf(")\n");e=0.0; /*以下是判定循环条件如果前后两次结果差值大于给定的误差限，则继续循环，此处也可以采用真值来进行比较，采用真值就需要再最开始给定真值*/for(i=0;i<n;i++)if(e<fabs(x[i]-z[i]))e=fabs(x[i]-z[i]);if(e>ej) goto cyclo;else{printf("at last %f ",max); /*输出循环中止时的最终结果，即矩阵的主特征值和对应的特征向量y[]*/printf("(");for(i=0;i<n;i++)printf("%f ",y[i]);printf(")\n");}}}（2）#include <stdio.h>#include <math.h>//a 存放n阶实对称阵，返回时对角线存放n个特征值//n 矩阵阶数//v n阶矩阵，第j列为对应第j个特征值的特征向量//eps 控制精度//jt 最大迭代次数//返回值：大于0则计算成功，小于0则失败int jacobi(double *a,int n,double *v,double eps,int jt) {int i,j,p,q,u,w,t,s,l;double fm,cn,sn,omega,x,y,d;l=1;for (i=0; i<=n-1; i++){ v[i*n+i]=1.0;for (j=0; j<=n-1; j++)if (i!=j) v[i*n+j]=0.0;}while (1){ fm=0.0;for (i=1; i<=n-1; i++)for (j=0; j<=i-1; j++){ d=fabs(a[i*n+j]);if ((i!=j)&&(d>fm)){ fm=d; p=i; q=j;}}if (fm<eps) return(1);if (l>jt) return(-1);l=l+1;u=p*n+q; w=p*n+p; t=q*n+p; s=q*n+q;x=-a[u]; y=(a[s]-a[w])/2.0;omega=x/sqrt(x*x+y*y);if (y<0.0) omega=-omega;sn=1.0+sqrt(1.0-omega*omega);sn=omega/sqrt(2.0*sn);cn=sqrt(1.0-sn*sn);fm=a[w];a[w]=fm*cn*cn+a[s]*sn*sn+a[u]*omega;a[s]=fm*sn*sn+a[s]*cn*cn-a[u]*omega;a[u]=0.0; a[t]=0.0;for (j=0; j<=n-1; j++)if ((j!=p)&&(j!=q)){ u=p*n+j; w=q*n+j;fm=a[u];a[u]=fm*cn+a[w]*sn;a[w]=-fm*sn+a[w]*cn;}for (i=0; i<=n-1; i++)if ((i!=p)&&(i!=q)){ u=i*n+p; w=i*n+q;fm=a[u];a[u]=fm*cn+a[w]*sn;a[w]=-fm*sn+a[w]*cn;}for (i=0; i<=n-1; i++){ u=i*n+p; w=i*n+q;fm=v[u];v[u]=fm*cn+v[w]*sn;v[w]=-fm*sn+v[w]*cn;}}return(1);}int main(){int i,j;double a[15][15];double v[15][15];printf("示例矩阵：\n");for(i=0;i<=14;i++){for(j=0;j<=14;j++){if(i==j)a[i][j]=4;else if(fabs(i-j)==1)a[i][j]=-1;elsea[i][j]=0;printf("%lf\t",a[i][j]);}printf("\n");}jacobi((double*)a,15,(double*)v,0.0001,100);for(i=0;i<=14;i++){printf("特征值为：%lf\t",a[i][i]);}return 0;}（3）#include<iostream>#include<fstream>#include<iomanip>using namespace std;#include<math.h>#define N 3 //矩阵的维数#define NUM 15 //QR分解次数#define SNUM 13 //用来控制输出格式void Print(long double A[N][N],long double Q[N][N],long double R[N][N]); //矩阵输出void QR(long double A[N][N],long double Q[N][N],long double R[N][N]); //QR分解void Multiplicate(long double A[N][N],long double R[N][N],long doubleQ[N][N]); //迭代，获得下一次矩阵A=QRint main(){int i,j;long double A[N][N];long double Q[N][N]={0};long double R[N][N]={0};cout<<"*************************************************************"<<endl;cout<<"输入初始矩阵:\n";cout<<"-------------------------------------------------------------"<<endl;for(i=0;i<N;i++)for(j=0;j<N;j++)cin>>A[i][j];cout<<"*************************************************************"<<endl;cout<<"输出每一步迭代过程: \n";cout<<"*************************************************************"<< endl;for(i=1;i<=NUM;i++){QR(A,Q,R);cout<<"第"<<i<<"步QR分解:\n";cout<<"-------------------------------------------------------------"<< endl;Print(A,Q,R);Multiplicate(A,R,Q);}cout<<"矩阵特征值为:\n";cout<<"-------------------------------------------------------------"<< endl;for(i=0;i<N;i++) //输出特征值cout<<"r["<<(i+1)<<"]="<<R[i][i]<<endl;cout<<"*************************************************************"<< endl;return 0;}void QR(long double A[N][N],long double Q[N][N],long double R[N][N]) {int i,j,k,m;long double temp;long double a[N],b[N];for(j=0;j<N;j++){for(i=0;i<N;i++){a[i]=A[i][j];b[i]=A[i][j];}for(k=0;k<j;k++){R[k][j]=0;for(m=0;m<N;m++)R[k][j]+=a[m]*Q[m][k];for(m=0;m<N;m++)b[m]-=R[k][j]*Q[m][k];}temp=0;for(i=0;i<N;i++)temp+=b[i]*b[i];R[j][j]=sqrt(temp);for(i=0;i<N;i++)Q[i][j]=b[i]/sqrt(temp);}}void Multiplicate(long double A[N][N],long double R[N][N],long doubleQ[N][N]){int i,j,k;long double temp;for(i=0;i<N;i++)for(j=0;j<N;j++){temp=0;for(k=0;k<N;k++)temp+=R[i][k]*Q[k][j];A[i][j]=temp;}}void Print(long double A[N][N],long double Q[N][N],long double R[N][N]) {int i,j;cout<<left;cout<<"矩阵A:\n";for(i=0;i<N;i++){for(j=0;j<N;j++)cout<<setw(SNUM)<<A[i][j];cout<<endl;}cout<<"-------------------------------------------------------------"<< endl;cout<<"矩阵Q:\n";for(i=0;i<N;i++){for(j=0;j<N;j++)cout<<setw(SNUM)<<Q[i][j];cout<<endl;}cout<<"-------------------------------------------------------------"<< endl;cout<<"矩阵R:\n";for(i=0;i<N;i++){for(j=0;j<N;j++)cout<<setw(SNUM)<<R[i][j];cout<<endl;}cout<<"*************************************************************"<< endl;}五调试与测试结果分析及总结（1）（2）（3）六分析及总结（1）参考幂法程序，写出求矩阵按模最小特征值的反幂法程序，并用计算上例的按模最小特征值与其对应的特征向量。