北航数值分析大作业第二题2014秋.docx

数值分析大作业2 参考前人程序基础上自己编写的，大家参考(copy)的时候又多了一个版本了，不要问我是谁，请叫我红领巾SY1405*1. 算法设计方案1.1. 定义矩阵、向量、复数的类为了能将矩阵、复数、向量作为函数的返回值，也方便计算，定义矩阵类Matrix、向量类Vector、复数类Complex；1.1.1. 矩阵类变量：double dataNN;用于存放数据int m;阶数，默认为N阶函数：Matrix(int x=N)构造x阶零矩阵，不输入x时，x默认为N阶，下同Matrix(double aNN,int x=N)用数组aNN的前x阶构造一个x阶矩阵Matrix(string statement,int x=N)当表述（statement）为I时，构造一个单位矩阵Matrix(string statement,double offset)当表述为InitializeA时，构造题中给出的A，注意此时函数的第二个参数写为0.0，不然会执行上一个函数；当表述为InitializeAWithOffset时，构造A-Ivoid show(ofstream &out)以精确到小数点后3位、固定小数点位置的形式在屏幕输出矩阵的结果；并以精确到小数点后11位、科学计数法的形式将结果输出到文件result.txt中Matrix operator-(Matrix a)同阶矩阵的减法运算Matrix operator*(Matrix B)同阶矩阵的乘法运算Matrix operator*(double a)矩阵与数a相乘Matrix operator/(double a)矩阵与数a相除Matrix& operator=(Matrix a)在程序中可直接写成this=a实现矩阵a向矩阵this赋值1.1.2. 向量类与建立矩阵类的思路基本相同，仅列出需要说明的：Vector(Matrix A,Vector u)利用矩阵A与向量u构造当前向量this；即this=Audouble operator*(Vector a)Vector operator*(double b)当*号右边是个向量a时，返回值为数，即返回thisTa；当*号右边是个数时，返回值为向量b thisMatrix operator/(Vector b)Vector operator/(double h)当/号右边是个向量b时，返回值为矩阵，即返回this bT；当/号右边是个数h时，返回值为向量this/h1.1.3. 复数类复数类用于存储特征值。1.2. 除主要函数外方便编程的小函数Matrix T(Matrix A)计算矩阵A的转置。double _sgn(double a,double dr)计算-sgnadr的值。void Root(Matrix A,Complex &s1,Complex &s2)计算m阶矩阵A最右下角的二阶子式的特征值。并将特征值写入复数s1, s2。1.3. 拟上三角化Matrix Hessenberg(Matrix A0)按照教材P61-62所述步骤将A0拟上三角化后返回Hessenberg矩阵A0(n-1)，重载了运算符之后在此函数中可直接调用运算符，代码形式上更加简洁明了。1.4. QR分解void QR(Matrix A0,Matrix &Q,Matrix &R,ofstream &out)按照教材P59所述步骤将A0做QR分解并储存Q与R。将第一步做QR分解的RQ矩阵输出以便分析。1.5. 带双步位移的QR方法Complex* sbQR(Matrix An_1)该函数返回类型为指向复数类的指针，这样做可以返回一组复数。上面流程图按照教材P63-64结合老师优化过的步骤整理而成，转换为代码形式，得：1.6. 高斯消去法Vector* MainRowGauss(double realN,int count)利用列主元素高斯消去法求解count个实特征值（储存在数组real中，real为占用count个实数的空间，使用动态申请空间new double实现，详见程序）对应的特征向量并返回特征向量组的指针。2. 源代码#include#include#include#include#define N 10#define L 10000using namespace std;class Matrixpublic:double dataNN;int m;Matrix(int x=N)m=x;for(int i=0;iN;i+)for(int j=0;jN;j+)dataij=0;Matrix(double aNN,int x=N)m=x;for(int i=0;iN;i+)for(int j=0;jN;j+)dataij=aij;Matrix(string statement,int x=N)/用于定义x阶单位阵Im=x;if(statement=I)for(int i=0;iN;i+)for(int j=0;jN;j+)if(i=j) dataij=1;else dataij=0;else cout请检查statementn;Matrix(string statement,double offset)/用于定义原始A与带偏移量offset的原始Am=N;double offsettemp=0;if(statement=InitializeA)offsettemp=0;else if(statement=InitializeAWithOffset)offsettemp=offset;else cout请检查statementn;for(int i=1;i=N;i+)for(int j=1;j=N;j+)if(i!=j) datai-1j-1=sin(0.5*i+0.2*j);else datai-1j-1=1.52*cos(i+1.2*j)-offsettemp;void show(ofstream &out)for(int i=0;im;i+)for(int j=0;jm;j+)/setiosflags(ios:fixed)表示不用科学计数法表示，即0.00001记为0.000而不是1.000e-5coutsetprecision(3)setiosflags(ios:fixed)dataijt;coutendl;cout-endl;for(int i=0;im;i+)for(int j=0;jm;j+)outsetprecision(11)setw(20)setiosflags(ios:scientific)dataijt;outendl;out-endl;Matrix operator-(Matrix a)Matrix Result(m);if(m!=a.m)cout阶数不同，两矩阵不能相减！n;return Result;for(int i=0;im;i+)for(int j=0;jm;j+)Result.dataij=dataij-a.dataij;return Result;Matrix operator*(Matrix B)Matrix Result(m);if(m!=B.m)cout阶数不同，两矩阵不能相乘！n;return Result;for(int i=0;im;i+)for(int j=0;jm;j+)double temp=0;for(int t=0;tm;t+)temp+=datait*B.datatj;Result.dataij=temp;return Result;Matrix operator*(double a)Matrix Result(m);for(int i=0;im;i+)for(int j=0;jm;j+)Result.dataij=dataij*a;return Result;Matrix operator/(double a)Matrix Result(m);for(int i=0;im;i+)for(int j=0;jm;j+)Result.dataij=dataij/a;return Result;Matrix& operator=(Matrix a)m=a.m;for(int i=0;im;i+)for(int j=0;jm;j+)dataij=a.dataij;return *this;class Vectorpublic:double dataN;int m;Vector(int x=N)m=x;for(int i=0;iN;i+)datai=0;Vector(double aN,int x=N)m=x;for(int i=0;iN;i+)datai=ai;Vector(Matrix A,Vector u)m=A.m;if(A.m!=u.m)cout阶数不同，矩阵与向量不能相乘！n;for(int i=0;iN;i+)datai=0;for(int i=0;im;i+)double sum=0;for(int j=0;jm;j+)sum+=A.dataij*u.dataj;datai=sum;Vector& operator=(Vector a)m=a.m;for(int i=0;im;i+)datai=a.datai;return *this;double operator*(Vector a)/this*a返回值为行向量thisT 乘以 列向量a，为一个数if(m!=a.m)cout阶数不同，矩阵与矩阵不能相乘返回一个数！n;return 0;double sum=0;for(int i=0;im;i+)sum+=datai*a.datai;return sum;Vector operator*(double b)Vector Result(m);for(int i=0;im;i+)Result.datai=datai*b;return Result;Matrix operator/(Vector b)/this/b返回值为列向量this 乘以 行向量b，为一个矩阵Matrix Result(m);if(m!=b.m)cout阶数不同，向量与向量不能相乘返回一个矩阵！n;return Result;for(int i=0;im;i+)for(int j=0;jm;j+)Result.dataij=datai*b.dataj;return Result;Vector operator/(double h)Vector Result(m);for(int i=0;im;i+)Result.datai=datai/h;return Result;Vector operator-(Vector b)Vector Result(m);if(m!=b.m)cout阶数不同，向量与向量不能相减！n;return Result;for(int i=0;im;i+)Result.datai=datai-b.datai;return Result;void show(ofstream &out)for(int i=0;im;i+)coutsetprecision(3)setiosflags(ios:fixed)datait;coutendl-endl;for(int i=0;im;i+)outsetprecision(11)setw(20)setiosflags(ios:scientific)datait;outendl-endl;class Complexpublic:double Re,Im;Complex()Re=0;Im=0;Complex(double re,double im)Re=re;Im=im;void show(ofstream &out)if(Im=0)coutReendl;else if(Im0)coutReImin;else coutRe+Imin;if(Im=0)outsetprecision(11)setw(20)setiosflags(ios:scientific)Reendl;else if(Im0)outReImin;else outRe+Imin;Matrix T(Matrix A)Matrix Result(A.m);for(int i=0;iA.m;i+)for(int j=0;j0) return -dr;else return dr;Matrix Hessenberg(Matrix A0)/*矩阵的阶数功能是在编程至双步位移时添加的，在此之前矩阵的阶数一直是N，本函数中阶数仍用N，未改为m，注意：此函数仅适用于A0为N阶矩阵*/Matrix A=A0;Vector Ur,Wr;double dr,hr,tr;bool flag;for(int r=0;r=N-3;r+)flag=0;/ar+1,r以下元素全零时flag=0for(int i=r+2;iN;i+)if(A.datair!=0) flag=1;break;if(flag=0) continue;/结束本次循环，执行下一次循环dr=0;for(int i=r+1;iN;i+)dr+=A.datair*A.datair;dr=sqrt(dr);dr=_sgn(A.datar+1r,dr);/cr与dr共用内存hr=dr*(dr-A.datar+1r);for(int i=0;ir+1;i+)Ur.datai=0;/*第r+1以前的元素置0，此步不能省略,因为Ur是在for循环外定义的，上次的使用的Ur残留有值存在。要不就在for循环里定义Ur。！做循环一定要注意是否此次循环所有变量都已重写再使用而不使用之前的残留值！*/Ur.datar+1=A.datar+1r-dr;for(int i=r+2;iN;i+)Ur.datai=A.datair;Vector Pr(T(A),Ur);Pr=Pr/hr;Vector Qr(A,Ur);Qr=Qr/hr;tr=Pr*Ur/hr;Wr=Qr-Ur*tr;A=A-Wr/Ur-Ur/Pr;return A;void QR(Matrix A0,Matrix &Q,Matrix &R,ofstream &out)/A0=QRR=A0;Matrix Qr(I);Vector Ur;bool flag;double dr,hr;for(int r=0;rN-1;r+)flag=0;/ar,r以下元素全零时flag=0for(int i=r+1;iN;i+)if(R.datair!=0) flag=1;break;if(flag=0) continue;/结束本次循环，执行下一次循环dr=0;for(int i=r;iN;i+)dr+=R.datair*R.datair;dr=sqrt(dr);dr=_sgn(R.datarr,dr);/cr与dr共用内存hr=dr*(dr-R.datarr);for(int i=0;ir;i+)Ur.datai=0;Ur.datar=R.datarr-dr;for(int i=r+1;iN;i+)Ur.datai=R.datair;Vector Wr(Qr,Ur);Qr=Qr-Wr/Ur/hr;Vector Pr(T(R),Ur);Pr=Pr/hr;R=R-Ur/Pr;if(r=0)cout对An_1做第一步QR时的RQ矩阵为：n;out=0)s1.Re=(-b+sqrt(d)/2;s2.Re=(-b-sqrt(d)/2;s1.Im=0;s2.Im=0;elses1.Re=(-b)/2;s2.Re=(-b)/2;s1.Im=(sqrt(-d)/2;s2.Im=(-sqrt(-d)/2;Complex* sbQR(Matrix An_1)Complex *s=new ComplexN;int k=1;Matrix A=An_1;while(k1e-12)if(fabs(A.datam-1m-2)1e-12)Root(A,sm,sm-1);A.m=A.m-2;m=A.m-1;if(m=0)sm.Re=A.datamm;sm.Im=0;return s;else if(m=1)Root(A,sm,sm-1);return s;else break;elsedouble ss=A.datam-1m-1+A.datamm;double tt=A.datam-1m-1*A.datamm-A.datamm-1*A.datam-1m;Matrix C=A;Matrix I(I,A.m);Matrix B=A*A-A*ss-I*(-tt);bool flag;double dr,hr,tr;Vector Ur(A.m),Wr(A.m);for(int r=0;rm;r+)flag=0;/ar,r以下元素全零时flag=0for(int i=r+1;i=m;i+)if(B.datair!=0) flag=1;break;if(flag=0) continue;/结束本次循环，执行下一次循环dr=0;for(int i=r;i=m;i+)dr+=B.datair*B.datair;dr=sqrt(dr);dr=_sgn(B.datarr,dr);/cr与dr共用内存hr=dr*(dr-B.datarr);for(int i=0;ir;i+)Ur.datai=0;Ur.datar=B.datarr-dr;for(int i=r+1;i=m;i+)Ur.datai=B.datair;Vector Vr(T(B),Ur);Vr=Vr/hr;B=B-Ur/Vr;Vector Pr(T(C),Ur);Pr=Pr/hr;Vector Qr(C,Ur);Qr=Qr/hr;tr=Pr*Ur/hr;Wr=Qr-Ur*tr;C=C-Wr/Ur-Ur/Pr;A=C;k+;sm.Re=A.datamm;sm.Im=0;A.m-;m-;if(m=0)sm.Re=A.datamm;sm.Im=0;return s;if(m=1)Root(A,sm,sm-1);return s;Vector* MainRowGauss(double real,int count)Vector *x=new Vectorcount;for(int num=0;numcount;num+)/消元过程Matrix A(InitializeAWithOffset,realnum);for(int k=0;kN-1;k+)int wheremax=k;for(int i=k+1;iA.datawheremaxk)wheremax=i;for(int j=k;jN;j+)double temp=A.datawheremaxj;A.datawheremaxj=A.datakj;A.datakj=temp;for(int i=k+1;iN;i+)double mik=A.dataik/A.datakk;for(int j=k+1;j=0;k-)double sum=0;for(int j=k+1;jN;j+)sum+=A.datakj*xnum.dataj;xnum.datak=-sum/A.datakk;return x;void main()ofstream out(result.txt);Matrix A(InitializeA,0.0);/*一定是0.0而不是0，因为重载了两个函数Matrix(string statement,double offset)/用于定义原始A与带偏移量offset的原始A与Matrix(string statement,int x=N)/用于定义x阶单位阵I*/cout初始矩阵A：n;out初始矩阵A：n;A.show(out);Matrix An_1=Hessenberg(A);cout拟上三角化后得到的矩阵A(n-1)：n;out拟上三角化后得到的矩阵A(n-1)：n;An_1.show(out);Matrix Q,R;QR(An_1,Q,R,out);cout对拟上三角阵A(n-1)进行QR分解后得到的Q：n;out对拟上三角阵A(n-1)进行QR分解后得到的Q：n;Q.show(out);cout对拟上三角阵A(n-1)进行QR分解后得到的R：n;out对拟上三角阵A(n-1)进行QR分解后得到的R：n;R.show(out);Complex *eigenvalue=sbQR(An_1);cout所有特征值：n;out所有特征值：n;for(int i=0;iN;i+)eigenvaluei.show(out);int count=0;for(int i=0;iN;i+)if(eigenvaluei.Im=0)count+;/先数出实特征值的个数double *real=new doublecount;/有几个实特征值real就申请几个空间count=0;for(int i=0;iN;i+)if(eigenvaluei.Im=0)realcount=eigenvaluei.Re;count+;Vector *Eigenvector=MainRowGauss(real,count);cout实特征值对应的特征向量：n;out实特征值对应的特征向量：n;for(int i=0;i=3的元素没变。如果执行for(int r=0;r=x+3的元素没有变，如程序注释所述。对拟上三角矩阵进行QR分解得到的Q矩阵：对拟上三角阵A(n-1)进行QR分解后得到的Q：-0.356 0.444 -0.694 0.066 0.370 0.187 -0.016 0.114 0.048 -0.054-0.934 -0.169 0.264 -0.025 -0.141 -0.071 0.006 -0.044 -0.018 0.021-0.000 -0.880 -0.401 0.038 0.214 0.108 -0.009 0.066 0.028 -0.031-0.000 0.000 -0.537 -0.126 -0.707 -0.358 0.031 -0.218 -0.093 0.1040.000 -0.000 0.000 0.989 -0.127 -0.064 0.006 -0.039 -0.017 0.0190.000 -0.000 0.000 -0.000 0.531 -0.685 0.059 -0.418 -0.177 0.1990.000 -0.000 0.000 -0.000 -0.000 -0.589 -0.096 0.677 0.287 -0.322-0.000 0.000 -0.000 0.000 0.000 0.000 -0.993 -0.100 -0.042 0.0480.000 -0.000 0.000 -0.000 -0.000 -0.000 0.000 0.538 -0.561 0.629-0.000 0.000 -0.000 0.000 -0.000 -0.000 -0.000 0.000 0.746 0.665-R矩阵：对拟上三角阵A(n-1)进行QR分解后得到的R：2.513 -2.181 -1.317 -0.032 -0.255 -1.263 -0.143 0.855 -0.406 0.3670.000 -1.973 0.228 0.701 0.595 0.534 -0.330 0.061 -0.284 -0.1020.000 0.000 2.407 1.723 -0.451 3.392 -0.112 -1.449 0.731 -0.4850.000 -0.000 0.000 1.596 0.640 0.350 -0.065 -0.399 0.239 -0.0910.000 -0.000 0.000 0.000 -1.456 -1.416 -0.274 0.282 0.031 0.2180.000 -0.000 0.000 0.000 0.000 1.240 0.144 -0.197 -0.550 -0.1560.000 0.000 -0.000 0.000 0.000 0.000 -0.800 0.324 -0.435 0.1300.000 -0.000 0.000 0.000 0.000 0.000 0.000 1.310 -0.452 -0.2540.000 -0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.659 0.4900.000 0.000 -0.000 0.000 0.000 0.000 0.000 -0.000 0.000 0.064-特征值：所有特征值：3.390-2.337-0.893i-2.337+0.893i1.590-1.493-0.989-0.109i-0.989+0.109i0.9430.0500.649实特征值对应的特征向量（按上面输出时的顺序）：实特征值对应的特征向量：-0.438 -0.909 -1.982 -1.083 -1.272 -1.083 0.363 1.692 2.128 1.000-0.622 -0.112 -2.521 -1.306 -3.809 8.133 -1.230 -0.675 2.712 1.000-15.996 22.177 0.410 -7.915 0.103 -0.072 -0.628 -0.335 -0.376 1.000-2.758 1.574 -0.632 -1.662 -12.110 7.179 -5.277 28.422 -12.308 1.000-4.889 -4.729 8.847 -0.712 -8.713 -2.863 14.745 -7.049 -6.769 1.000-4.456 2.932 15.720 -1.936 -29.499 7.460 -9.288 15.961 11.906 1.000-请按任意键继续. . .4. 编程时的问题与思考动手编程前详细阅读了教材相关知识以及网上搜索的前人所编的程序。发现大多数都程序都没有使用类，这样就不能使函数返回值为向量或矩阵。而且大多数前人的程序在计算矩阵与向量的乘法时也是直接拆成元素之间的运算。本程序用运算符重载将元素之间的运算封装起来，形式上简洁，计算量上并没有增加，只是在执行时多出一小部分调用重载运算符的时间。运算符重载也是类带来的便利之处。1) 遇到的第一个问题，就是按教材P239的表达式初始化矩阵A后得到的矩阵跟别人的结果不一样。于是推测其原因是没有#includ