北航数值分析第二次大作业--QR分解

一、算法设计方案：按题目要求，本程序运用带双步位移的QR方法求解给定矩阵的特征值，并对每一实特征值，求解其相应的特征向量。

总体思路：1）初始化矩阵首先需要将需要求解的矩阵输入程序。

为了防止矩阵在后面的计算中被破坏保存A[][]。

2）对给定的矩阵进行拟上三角化为了尽量减少计算量，提高程序的运行效率，在对矩阵进行QR分解之前，先进行拟上三角化。

由于矩阵的QR 分解不改变矩阵的结构，所以具有拟上三角形状的矩阵的QR分解可以减少大量的计算量。

这里用函数void QuasiTriangularization()来实现，函数形参为double型N维方阵double a[][N]。

3）对拟上三角化后的矩阵进行QR分解对拟上三角化的矩阵进行QR分解会大大减小计算量。

用子程序void QR_decomposition()来实现，将Q、R设为形参，然后将计算出来的结果传入Q和R，然后求出RQ乘积。

4）对拟上三角化后的矩阵进行带双步位移的QR分解为了加速收敛，对QR分解引入双步位移，适当选取位移量，可以避免进行复数运算。

为了进一步减少计算量，在每次进行QR分解之前，先判断是否可以直接得到矩阵的一个特征值或者通过简单的运算得到矩阵的一对特征值。

若可以，则得到特征值，同时对矩阵进行降阶处理；若不可以，则进行QR分解。

这里用函数intTwoStepDisplacement_QR()来实现。

这是用来存储计算得到的特征值的二维数组。

考虑到特征值可能为复数，因此将所有特征值均当成复数处理。

此函数中，QR分解部分用子函数void QR_decompositionMk()实现。

这里形参有三个，分别用来传递引入双步位移后的Mk阵，A矩阵，以及当前目标矩阵的维数m。

5）计算特征向量得到特征值后，计算实特征值相应的特征向量。

这里判断所得特征值的虚数部分是否为零。

求实特征值的特征向量采用求解相应的方程组((A-λI)x=0)的方法。

因此先初始化矩阵Array，计算(A-λI)，再求解方程组。

应用数理统计作业二学号：姓名：电话：二〇一四年十二月对NBA球队的聚类分析和判别分析摘要：NBA联盟作为篮球的最高殿堂深受广大球迷的喜爱，联盟的30支球队大家也耳熟能详，本文选取NBA联盟30支球队2013-2014常规赛赛季场均数据。

利用spss软件通过聚类分析对27个地区进行实力类型分类，并利用判断分析对其余3支球队对分类结果进行验证。

可以看出各球队实力类型与赛季实际结果相吻合。

关键词：聚类分析，判别分析，NBA目录1. 引言 (4)2、相关统计基础理论 (5)2.1、聚类分析 (5)2.2，判别分析 (6)3.聚类分析 (7)3.1数据文件 (7)3.2聚类分析过程 (9)3.3 聚类结果分析 (11)4、判别分析 (12)4.1 判别分析过程 (12)4.2判别检验 (17)5、结论 (20)参考文献 (21)致谢 (22)1. 引言1896年，美国第一个篮球组织"全国篮球联盟（简称NBL）"成立，但当时篮球规则还不完善，组织机构也不健全，经过几个赛季后，该组织就名存实亡了。

1946年4月6日，由美国波士顿花园老板沃尔特.阿.布朗发起成立了“美国篮球协会”（简称BAA）。

1949年在布朗的努力下，美国两大篮球组织BAA和NBL合并为“全国篮球协会”（简称NBA）。

NBA季前赛是 NBA各支队伍的热身赛，因为在每个赛季结束后，每支球队在阵容上都有相当大的变化，为了让各队磨合阵容，熟悉各自球队的打法，确定各队新赛季的比赛阵容、同时也能增进队员、教练员之间的沟通，所以在每个赛季开始之前，NBA就举办若干场季前赛，使他们能以比较好的状态投入到漫长的常规赛的比赛当中。

为了扩大NBA在全球的影响，季前赛有约三分之一的球队在美国以外的国家举办。

从总体上看，NBA的赛程安排分为常规赛、季后赛和总决赛。

常规赛采用主客场制，季后赛和总决赛采用七场四胜制的淘汰制。

[31]NBA常规赛从每年的11月的第一个星期二开罗，到次年的4月20日左右结束。

北航数值分析全部三次大作业第一次大作业是关于解线性方程组的数值方法。

我们被要求实现各种常用的线性方程组求解算法，例如高斯消元法、LU分解法和迭代法等。

我首先学习了这些算法的原理和实现方法，并借助Python编程语言编写了这些算法的代码。

在实验中，我们使用了不同规模和条件的线性方程组进行测试，并比较了不同算法的性能和精度。

通过这个作业，我深入了解了线性方程组求解的原理和方法，提高了我的编程和数值计算能力。

第二次大作业是关于数值积分的方法。

数值积分是数值分析中的重要内容，它可以用于计算曲线的长度、函数的面积以及求解微分方程等问题。

在这个作业中，我们需要实现不同的数值积分算法，例如矩形法、梯形法和辛普森法等。

我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。

在实验中，我们计算了不同函数的积分值，并对比了不同算法的精度和效率。

通过这个作业，我深入了解了数值积分的原理和方法，提高了我的编程和数学建模能力。

第三次大作业是关于常微分方程的数值解法。

常微分方程是数值分析中的核心内容之一，它可以用于描述众多物理、化学和生物现象。

在这个作业中，我们需要实现不同的常微分方程求解算法，例如欧拉法、龙格-库塔法和Adams法等。

我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。

在实验中，我们解决了一些具体的常微分方程问题，并比较了不同算法的精度和效率。

通过这个作业，我深入了解了常微分方程的原理和方法，提高了我的编程和问题求解能力。

总的来说，北航数值分析课程的三次大作业非常有挑战性，但也非常有意义。

通过这些作业，我在数值计算和编程方面得到了很大的提升，也更加深入地了解了数值分析的理论和方法。

虽然这些作业需要大量的时间和精力，但我相信这些努力将会对我未来的学习和工作产生积极的影响。

北京航空航天大学数值分析课程知识点总结1.2 误差知识与算法知识1.2.2 绝对误差、相对误差与有效数字设a 是准确值x 的一个近似值，记e x a =-，称e 为近似值a 的绝对误差，简称误差。

如果||e 的一个上界已知，记为ε，即||e ε≤，则称ε为近似值a 的绝对误差限或绝对误差界，简称误差限或误差界。

记r e x ae x x-==，称r e 为近似值a 的相对误差。

由于x 未知，实际上总把e a 作为a 的相对误差，并且也记为r e x a e a a -==，相对误差一般用百分比表示。

r e 的上界，即||r a εε=称为近似值a 的相对误差限或相对误差界。

定义设数a 是数x 的近似值。

如果a 的绝对误差限时它的某一位的半个单位，并且从该位到它的第一位非零数字共有n 位，则称用a 近似x 时具有n 位有效数字。

1.2.3 函数求值的误差估计设()u f x =存在足够高阶的导数，a 是x 的近似值，则~()u f a =是()u f x =的近似值。

若'()0f a ≠且|''()|/|'()|f a f a 不很大，则有误差估计~~()'()()()'()()e uf a e a u f a a εε≈≈。

若(1)()'()''()...()0,()0k k f a f a fa f a -====≠，且比值(1)()()/()k k f a fa +不很大，则有误差估计[][]()~()~()()()!()()()!k kk k f a e u e a k f a u a k εε≈≈。

对于n 元函数，有误差估计~121~121(,,...,)()()(,,...,)()()nn i i i nn i i if a a a e u e a x f a a a u a x εε==?≈??≈?∑∑；若一阶偏导全为零或很小，则要使用高阶项。

《数值分析B》课计算实习第一题设计文档与源程序姓名：杨彦杰学号：SY10171341 算法的设计方案（1）运行平台操作系统：Windows XP；开发平台：VC6.0++；工程类型：文档视图类；工程名：Numanalysis；（2）开发描述首先新建类CMetrix，该类完成矩阵之间的相关运算，包括相乘、加减等，以主程序方便调用；题目的解算过程在视图类CNumanalysisView中实现，解算结果在视图界面中显示；（3）运行流程（4）运行界面2、全部源代码（1）类CMetrixMetrix.h文件：class CMetrix{public:double** MetrixMultiplyConst(double**A,int nRow,int nCol,double nConst);//矩阵乘常数double** MetrixMultiplyMetrix(double**A,double**mA,int nRow,int nCol);//矩阵相乘double** MetrixSubtractMetrix(double **A, double **subA, int nRow,int nCol);//矩阵减矩阵double VectorMultiplyVector(double*V,double*mulV,int nV);//向量点积double** VectorMultiplyVectortoMetrix(double*V,double*VT,int nV);//向量相乘为矩阵double* VectorSubtractVector(double*V,double*subV,int nV);//向量相减double* VectorMultiplyConst(double *V, int nV, double nConst);//向量乘常数double LengthofVector(double *V,int nV);//求向量的长度double* MetrixMultiplyVector(double**A,int nRow,int nCol,double*V,int nV);//矩阵与向量相乘double** AtoAT(double **A,int Row,int Col);//矩阵转置运算void FreeMem();CMetrix(int nRow,int nCol);uCMetrix();virtual ~CMetrix();double* vector; //过渡向量double** B; //过渡矩阵};Metrix.cpp文件：CMetrix::CMetrix(int nRow, int nCol){B = new double*[nRow];for (int i = 0;i < nCol;i++){B[i] = new double[nCol];}vector = new double[nRow];}CMetrix::~CMetrix(){delete vector;B = NULL;delete B;}double** CMetrix::AtoAT(double **A, int nRow, int nCol){for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){B[col][row] = A[row][col];}}return B;}double* CMetrix::MetrixMultiplyVector(double **A, int nRow, int nCol, double *V, int nV) {if (nCol != nV){AfxMessageBox("矩阵列数和向量维数不等,不能相乘!");return 0;}double sum = 0.0;for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){sum += A[row][col]*V[col];}vector[row] = sum;sum = 0.0;}return vector;}double CMetrix::LengthofVector(double *V, int nV){double length = 0.0;for (int col = 0;col < nV;col++){length += V[col]*V[col];}return length;}double* CMetrix::VectorMultiplyConst(double *V, int nV, double nConst){for (int col = 0;col < nV;col++){vector[col] = V[col]*nConst;}return vector;}double* CMetrix::VectorSubtractVector(double *V, double *subV, int nV){for (int col = 0;col < nV;col++){vector[col] = V[col]-subV[col];}return vector;}double** CMetrix::VectorMultiplyVectortoMetrix(double*V, double *VT, int nV){for (int row = 0;row < nV;row++){for (int col = 0;col < nV;col++){B[row][col] = V[row]*VT[col];}}return B;}double CMetrix::VectorMultiplyVector(double *V, double *mulV, int nV){double length = 0.0;for (int col = 0;col < nV;col++){length += V[col]*mulV[col];}return length;}double** CMetrix::MetrixSubtractMetrix(double **A, double **subA, int nRow, int nCol) {for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){B[row][col] = A[row][col]-subA[row][col];}}return B;}double** CMetrix::MetrixMultiplyMetrix(double **A, double **mA, int nRow, int nCol) {double sum = 0.0;for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){for(int n = 0;n < nCol;n++){sum += A[row][n]*mA[n][col];}B[row][col] = sum;sum = 0.0;}}return B;}double** CMetrix::MetrixMultiplyConst(double **A, int nRow, int nCol, double nConst) {for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){B[row][col] = A[row][col]*nConst;}}return B;}（2）类CNumanalysisViewNumanalysisview.hclass CNumanalysisView : public CEditView{…………public:double Sign(double x);void DisplayVector(double*V,int nV); // 显示向量数据void DisplayMetrix(double **A,int Row,int Col); //显示矩阵void DisplayText(CString str); //显示文本protected://{{AFX_MSG(CNumanalysisView)afx_msg void OnQRanalyze(); //运行主函数…………};Numanalysisview.cppvoid CNumanalysisView::OnQRanalyze(){//开辟空间int nRow = 10;int nCol = 10;CString str;CMetrix Metrix(nRow,nCol);double tempa = 0.0;double *V = new double[nCol]; //分配10*10矩阵空间double *ur = new double[nCol];double *pr = new double[nCol];double *qr = new double[nCol];double *wr = new double[nCol];double *tempV = new double[nCol];double **Ar = new double*[nRow];double **C = new double*[nRow];double **Cr = new double*[nRow];double **tempA = new double*[nRow];double **A = new double*[nRow];double **R = new double*[nRow];for (int col = 0;col < nRow;col++){A[col] = new double[nCol];Ar[col] = new double[nCol];C[col] = new double[nCol];Cr[col] = new double[nCol];tempA[col] = new double[nCol];R[col] = new double[nCol];}//矩阵A求解for (int i = 0;i < nRow;i++){for (int j = 0;j < nCol;j++){if(i == j)A[i][j] = 1.5*cos((i+1.0)+1.2*(j+1.0));elseA[i][j] = sin(0.5*(i+1.0)+0.2*(j+1.0));}}//--------------------拟上三角化-------------------------// double dr = 0.0,cr = 0.0,hr = 0.0,tr = 0.0;for (int r = 0;r < nCol - 2;r++){dr = 0.0;for (i = r+1;i < nCol;i++) //dr{dr += A[i][r]*A[i][r];}dr = sqrt(dr);for (i = r+2;i < nCol;i++) //判断air是否全为零tempa += fabs(A[i][r]);if (tempa <= IPSLEN)continue;if (A[r+1][r] == 0.0) //crcr = dr;elsecr = -1*Sign(A[r+1][r])*dr;hr = cr*cr - cr*A[r+1][r]; //hrstr.Format("dr = %.6e, cr = %.6e, hr = %.6e",dr,cr,hr);for (int row = 0;row < nRow;row++) //ur{if (row < r+1)ur[row] = 0.0;else if (row == r+1)ur[row] = A[row][r]-cr;elseur[row] = A[row][r];}tempA = Metrix.AtoAT(A,nRow,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)Ar[row][col] = tempA[row][col];}tempV = Metrix.MetrixMultiplyVector(Ar,nRow,nCol,ur,nCol); //pr memcpy(pr,tempV,nCol*8);tempV = Metrix.VectorMultiplyConst(pr,nCol,1.0/hr);memcpy(pr,tempV,nCol*8);tempV = Metrix.MetrixMultiplyVector(A,nRow,nCol,ur,nCol); //qr memcpy(qr,tempV,nCol*8);tempV = Metrix.VectorMultiplyConst(qr,nCol,1.0/hr);memcpy(qr,tempV,nCol*8);tr = Metrix.VectorMultiplyVector(pr,ur,nCol)/hr; //trtempV = Metrix.VectorMultiplyConst(ur,nCol,tr); //wr memcpy(wr,tempV,nCol*8);tempV = Metrix.VectorSubtractVector(qr,wr,nCol);memcpy(wr,tempV,nCol*8);tempA = Metrix.VectorMultiplyVectortoMetrix(wr,ur,nCol); //Arfor (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(A,Ar,nRow,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)A[row][col] = tempA[row][col];}tempA = Metrix.VectorMultiplyVectortoMetrix(ur,pr,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(A,Ar,nRow,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++){A[row][col] = tempA[row][col];if (fabs(A[row][col]) < IPSLEN){A[row][col] = 0.0;}}}}DisplayText("矩阵A拟上三角化后所得的矩阵为：");DisplayMetrix(A,nRow,nCol);for (int row = 0;row < nRow;row++) //用于计算特征向量{for (col = 0;col < nCol;col++)R[row][col] = A[row][col];}// -------------------------------------------------////--------------------带双步位移的QR分解-------------------------// int m = nCol;struct EigenVal //定义特征值结构，实数和虚数{double Realnum;double Imagnum;};EigenVal *eigenvalue = new EigenVal[m];EigenVal tmpEigen1,tmpEigen2;double b = 0.0,c = 0.0,delta = 0.0,s = 0.0,t = 0.0;double *vr = new double[m];for (int k = 1;k < 100; k++){//m代表矩阵阶数，判断式中直接用，运算中需要-1while (m > 1 && fabs(A[m-1][m-2]) <= IPSLEN)//第三步和第四步{eigenvalue[m-1].Realnum = A[m-1][m-1];eigenvalue[m-1].Imagnum = 0.0;m = m - 1;}if (m == 1){eigenvalue[m-1].Realnum = A[m-1][m-1];eigenvalue[m-1].Imagnum = 0.0;DisplayText("已求出A的全部特征值：");break;}b = -(A[m-2][m-2]+A[m-1][m-1]); //第五步求一元二次方程式的根s1,s2c = A[m-2][m-2]*A[m-1][m-1]-A[m-2][m-1]*A[m-1][m-2];delta =b*b - 4*c;if (delta >= 0.0){tmpEigen1.Realnum = (-b-sqrt(delta))/2;tmpEigen1.Imagnum = 0.0;tmpEigen2.Realnum = (-b+sqrt(delta))/2;tmpEigen2.Imagnum = 0.0;}else{tmpEigen1.Realnum = -b/2;tmpEigen1.Imagnum = -sqrt(fabs(delta))/2 ;tmpEigen2.Realnum = -b/2;tmpEigen2.Imagnum = sqrt(fabs(delta))/2;}if (m == 2) //第六步 m=2时结束运算{eigenvalue[m-1] = tmpEigen1;eigenvalue[m-2] = tmpEigen2;DisplayText("已求出A的全部特征值：");break;}else //第七步 m > 1{if (fabs(A[m-2][m-3]) <= IPSLEN){eigenvalue[m-1] = tmpEigen1;eigenvalue[m-2] = tmpEigen2;m = m - 2;continue;}}for (int row = 0;row < m;row++) //Mk求之前需要把A付给C{for (int col = 0;col < m;col++)C[row][col] = A[row][col];}double **I = new double*[m]; //第九步求Mk和Mk的QR分解for (int i = 0;i < m;i++) //求单位矩阵I,分配m*m矩阵空间{I[i] = new double[m];}for (i = 0;i < m;i++){for (int j = 0;j < m;j++){if(i == j)I[i][j] = 1;else I[i][j] = 0;}}s = A[m-2][m-2]+A[m-1][m-1];t = A[m-2][m-2]*A[m-1][m-1] - A[m-2][m-1]*A[m-2][m-1];tempA = Metrix.MetrixMultiplyMetrix(A,A,m,m);//A*Afor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixMultiplyConst(A,m,m,s);//s*Afor (row = 0;row < m;row++){for (col = 0;col < m;col++)A[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(Ar,A,m,m);//A*A-s*Afor (row = 0;row < m;row++){for (col = 0;col < m;col++)A[row][col] = tempA[row][col]; }tempA = Metrix.MetrixMultiplyConst(I,m,m,-1*t);//-t*Ifor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col]; }tempA = Metrix.MetrixSubtractMetrix(A,Ar,m,m);//A*A - s*A + r*I for (row = 0;row < m;row++){for (col = 0;col < m;col++){A[row][col] = tempA[row][col];if (fabs(A[row][col]) < IPSLEN){A[row][col] = 0.0;}}}delete I;//Mk的QR分解for (int r = 0;r < m - 1;r++){dr = 0.0;for (i = r;i < m;i++) //dr{dr += A[i][r]*A[i][r];}dr = sqrt(dr);for (i = r+1;i < m;i++) //判断air是否全为零tempa += fabs(A[i][r]);if (tempa <= IPSLEN)continue;if (A[r][r] == 0.0) //crcr = dr;elsecr = -1*Sign(A[r][r])*dr;hr = cr*cr - cr*A[r][r]; //hrfor (int row = 0;row < m;row++) //ur{if (row < r)ur[row] = 0.0;else if (row == r)ur[row] = A[row][r]-cr;elseur[row] = A[row][r];}tempA = Metrix.AtoAT(A,m,m); //Btfor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col];}tempV = Metrix.MetrixMultiplyVector(Ar,m,m,ur,m); //Bt*ur memcpy(vr,tempV,m*8);tempV = Metrix.VectorMultiplyConst(vr,m,1.0/hr); //vr = Bt*ur/hr memcpy(vr,tempV,m*8);tempA = Metrix.VectorMultiplyVectortoMetrix(ur,vr,m);//Ur*vrfor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(A,Ar,m,m); //Br-ur*vrfor (row = 0;row < m;row++){for (col = 0;col < m;col++){A[row][col] = tempA[row][col];if (fabs(A[row][col]) < IPSLEN){A[row][col] = 0.0;}}}tempA = Metrix.AtoAT(C,m,m); //Ctfor (row = 0;row < m;row++){for (col = 0;col < m;col++)Cr[row][col] = tempA[row][col]; }tempV = Metrix.MetrixMultiplyVector(Cr,m,m,ur,m); //pr memcpy(pr,tempV,m*8);tempV = Metrix.VectorMultiplyConst(pr,m,1.0/hr);memcpy(pr,tempV,m*8);tempV = Metrix.MetrixMultiplyVector(C,m,m,ur,m); //qr memcpy(qr,tempV,m*8);tempV = Metrix.VectorMultiplyConst(qr,m,1.0/hr);memcpy(qr,tempV,m*8);tr = Metrix.VectorMultiplyVector(pr,ur,m)/hr; //trtempV = Metrix.VectorMultiplyConst(ur,m,tr); //wr memcpy(wr,tempV,m*8);tempV = Metrix.VectorSubtractVector(qr,wr,m);memcpy(wr,tempV,m*8);tempA = Metrix.VectorMultiplyVectortoMetrix(wr,ur,m);//Cr+1for (row = 0;row < m;row++){for (col = 0;col < m;col++)Cr[row][col] = tempA[row][col]; }tempA = Metrix.MetrixSubtractMetrix(C,Cr,m,m);for (row = 0;row < m;row++){for (col = 0;col < m;col++)C[row][col] = tempA[row][col]; }tempA = Metrix.VectorMultiplyVectortoMetrix(ur,pr,m);for (row = 0;row < m;row++){for (col = 0;col < m;col++)Cr[row][col] = tempA[row][col]; }tempA = Metrix.MetrixSubtractMetrix(C,Cr,m,m);for (row = 0;row < m;row++){for (col = 0;col < m;col++){C[row][col] = tempA[row][col];if (fabs(C[row][col]) < IPSLEN){C[row][col] = 0.0;}}}}str.Format("矩阵A%d QR分解结束后所得到的矩阵为：",m);//计算结果输出DisplayText(str);DisplayMetrix(A,m,m);for (row = 0;row < m;row++) //Mk的QR分解后需要把C付给A{for (col = 0;col < m;col++)A[row][col] = C[row][col];}str.Format("迭代完成后的矩阵A%d = ",k);DisplayText(str);DisplayMetrix(A,m,m);}DisplayText("矩阵A的全体特征值如下: ");for (i = 0;i<nCol;i++){str.Format("%.6e + j%.6e",eigenvalue[i].Realnum,eigenvalue[i].Imagnum);DisplayText(str);}// -------------------------------------------------//求实特征值的特征向量，在拟上三角矩阵基础上直接求解即可////(A-egiI)X = 0.0;m = nRow;for (row = 0;row < nRow;row++) //用于计算特征向量{for (col = 0;col < nCol;col++)A[row][col] = R[row][col];}double **I = new double*[m]; //求单位矩阵I,分配m*m矩阵空间double sum = 0.0;for (i = 0;i < m;i++){I[i] = new double[m];}for (i = 0;i < m;i++){for (int j = 0;j < m;j++){if(i == j)I[i][j] = 1;else I[i][j] = 0;}}for (i = 0;i < nRow;i++){if (eigenvalue[i].Imagnum != 0.0){str.Format("特征值%.6e+j%.6e为虚数，不需要求特征向量。

数值分析第二题 梁进明SY0906529算法设计方案。

一．矩阵的QR 分解。

把矩阵A 分解为一个正交矩阵Q 与一个上三角矩阵R 的乘积，称为矩阵A 的正三角分解，简称QR 分解。

QR 分解的算法如下：记1A A =，并记[]rij n n Ar a ⨯=，令1Q I =（n 阶单位矩阵） 对于r=1,2,…,n-1执行(1) 若(1,2,...,)rir a i r r n =++全为零，则令1r r Q Q +=，1r r A A +=转(5);否则转(2)(2) 计算2r d =()sgn()r r rr r c a d =-(若()0r rr a =,则取r r c d =) 2()r r r r rrh c c a =- (3) 令()()()1,(0,...,0,,,...,)r r r T nr rr r r r nr u a c a a R +=-∈ (4) 计算11//r r rTr r r r rT r rr rT r r r rQ u Q Q u h p A u h A A u p ωω++==-==-(5) 继续当此算法执行完后就得到正交矩阵n Q Q =和上三角矩阵n R A =且有A QR =。

二．矩阵的 拟上三角化。

对实矩阵A 的拟上三角化具体算法如下：记(1)AA =，并记()r A 的第r 列到第n 列的元素为(1,2,...,;,1,...,)rij a i n j r r n ==+。

对于1,2,...,2r n =-执行（1） 若()(2,3,...,)r ir a i r r n =++全为零，则令(1)()r r AA +=,转(5);否则转(2)。

（2） 计算2()()1,1,2()1,sgn()(0,)r r r r r r r r r r r r r r r r rd c a d a c d h c c a +++==-===-若则取 （3） 令()()()1,2,(0,...,0,,,...,)r r r T nr r r r r r nr u a c a a R ++=-∈ （4） 计算()()(1)()///r T r r r r r r rTr r r rr r r rr r T T r r r rp A u h q A u h t p u h q t u A A u u p ωω+====-=--（5） 继续算法执行完后，就得到与原矩阵A 相似的拟上三角矩阵(1)n A -。

数值分析大作业QR分解题目:设计求取n n ?实数矩阵A 的所有特征值及其特征向量的数值算法，并以矩阵20010-1-24A=0-2131431?? ? ? ? ???为例进行具体的求解计算。

一、算法分析：一般而言，求取实数矩阵所有特征值的方法有雅克比法和QR 分解法，两者都是变换法。

其中雅克比法只能求解对称矩阵的全部特征值和特征向量，而QR 则可用于更一般的矩阵特征值的求解，结合反幂法可进而求出各个特征向量。

二、算法设计：1、原始实矩阵A Rn n∈拟上三角化为了减少求特征值和特征向量过程中的计算量，对生成的矩阵A 进行拟上三角化，得到拟上三角化矩阵A ’记A (1)=A ，并记A (r)的第r 列到第n 列的元素(1,2,...,;,1,...,)rij a i n j r r n ==+。

对于r=1,2,…,n -2执行(1) 若()(2,3,...,)r ir a i r r n =++全为零，则令A (r+1)= A (r)，转(5)；否则转(2)。

(2) 计算令()2()()1,1,1,sgn(0,sgn()=1)r r r r r rr r r r r c aa a ρ+++=-=，（若则取(3) 令-0=r n r u u ?? ?，()()()-1,2,1(,,...,)r r r Tn r r r r r r nr u a c a a ρ++=- (4) 计算(r)(r)(r)Tn-r r+1,r r r+2,r nr r Tn-r T n-r n-r n-r n-r r+1r 1u =(a -c ,a ,...,a )ρI H =I -2uu =H H =I -2u u A =HA H(5) 继续算法执行完后，就得到与原矩阵A 相似的拟上三角矩阵A (n-1)。

2、拟上三角矩阵QR 分解的求原矩阵的全部特征值记A k 是对拟上三角矩阵进行QR 算法，产生的矩阵序列，A 0是起始拟上三角矩阵，对于k =1,2,…,n -1执行： (1) 分解k-1k-1k-1A =R Q选取旋转矩阵1P =R(2,1,θ)，使(1)1k-1A =PA ，从而使第一列次对角元(1)2,1=0a ；再选取旋转矩阵2P =R(3,2,θ)，使(2)(1)1A =PA ，从而使第一列次对角元(2)3,2=0a ……如此进行下去，最多经过n-1步，(n-1)A必然变为上三角矩阵k-1R ，即k-1n-121k-1-1=-1-1-1k-1n-12112n-1R =P P P A =P P P P P PQ （）k-1Q 必为正交矩阵，且满足k-1k-1k-1A =R Q(2) 计算k k-1k-1A =R Q(3) 上述两过程反复进行直到k A 变为近似舒尔矩阵，对角线元素即为A 的近似特征值。

数值分析实习二院（系）名称航空科学与工程学院专业名称动力工程及工程热物理学号SY0905303学生姓名解立垚1. 题目试用带双步位移QR 的分解法求矩阵A=[a ij ]10*10的全部特征值，并对其中的每一个实特征值求相应的特征向量。

已知()sin 0.50.2,1.5cos 1.2,ij i j i j a i j i j ⎧⎫+≠⎪⎪=⎨⎬+=⎪⎪⎩⎭(),1,2,...,10i j =。

说明：1、求矩阵特征值时，要求迭代的精度水平为1210ε-=。

2、打印以下内容：算法的设计方案；全部源程序(要求注明主程序和每个子程序的功能)； 矩阵A 经过拟上三角话之后所得的矩阵()1n A -；对矩阵()1n A-进行QR 分解方法结束后所得的矩阵；矩阵A 的全部特征值()(),1,2,......10i i iR I i λ=，和A 的相应于实特征值的特征向量；其中()(),.i e i m i R R I I λλ==如果i λ是实数，则令0.i I =3、采用e 型输出数据，并且至少显示12位有效数字。

2. 算法设计方案本题采用带双步位移的QR 分解方法。

为了使程序简洁，自定义类Xmatrix ，其中封装了所需要的函数方法。

在Xmatrix 类中封装了运算符重载的函数，即定义了矩阵的加、减、乘、除、数乘运算及转置运算(T())。

同时为了避免传递数组带来的额外内存开销，使用引用(&)代替值传递，以节省内存空间，避免溢出.（1）此程序的主要部分为Xmatrix 中的doubleQR()方法，具体如下：Step1：使用矩阵拟上三角化的算法将A 化为拟上三角阵A （n-1）(此处调用Xmatrix 中的preQR()方法)Step2：令121,,10k m n ε-===， 其中k 为迭代次数。

Step3：如果,1m m a ε-≤，则得到A 的一个特征值,m m a ，令1m m =-，goto Step4;否则goto Step5.Step4： 如果1m =，则得到A 的一个特征值11a ，goto Step11;如果0m =，则goto Step11;如果1m >，则goto Step3;Step5(Step6)：如果2m =，则得到A 的两个特征值12s s 和（12s s 和为右下角两阶子阵对应的特征方程21,1,()det 0m m m m a a D λλ---++=的两个根。

数理统计大作业（二）全国各省发展程度的聚类分析及判别分析指导教师院系名称材料科学与工程院学号学生姓名2015 年 12 月21 日目录全国各省发展程度的聚类分析及判别分析 (1)摘要： (1)引言 (1)1实验方案 (2)1.1数据统计 (2)1.2聚类分析 (3)1.3判别分析 (4)2结果分析与讨论 (5)2.1聚类分析结果 (5)2.2聚类分析结果分析： (8)2.3判别分析结果 (9)2.4 Fisher判别结果分析： (11)参考文献： (16)全国各省发展程度的聚类分析及判别分析摘要：利用SPSS软件对全国31个省、直辖市、自治区(浙江、安徽、甘肃除外)的主要经济指标进行多种聚类分析，分析选择最佳聚类类数，并对浙江、湖南、甘肃进行类型判别分析。

通过这两个方法对全国各省进行发展分类。

本文选取了7项社会发展指标作为决定发展程度的影响因素，其中经济因素为主要因素，同时评估城镇化率和人口素质因素。

各项数据均来自2014年国家统计年鉴。

分析结果表明：北京市和上海市和天津市为同一类；江苏省和山东省和广东省为同一类型；河北、湖北、河南、湖南、四川、辽宁为同一类；其余的为另一类。

关键词：聚类分析、判别分析、发展引言聚类分析是根据研究对象的特征对研究对象进行分类的多元统计分析技术的总称。

它直接比较各事物之间的性质，将性质相近的归为一类，将性质差别较大的归入不同的类。

系统聚类分析又称集群分析，是聚类分析中应用最广的一种方法，它根据样本的多指标（变量）、多个观察数据，定量地确定样品、指标之间存在的相似性或亲疏关系，并据此连结这些样品或指标，归成大小类群，构成分类树状图或冰柱图。

判别分析是根据多种因素(指标)对事物的影响来实现对事物的分类，从而对事物进行判别分类的统计方法。

判别分析适用于已经掌握了历史上分类的每一个类别的若干样品，希望根据这些历史的经验（样品），总结出分类的规律性（判别函数）来指导未来的分类。

2020级数值分析第二次作业(非线性方程求根)参考答案和评分标准班级学号姓名一(20分)用二分法求方程3()10f x x x =--=在区间[1.0,1.5]内的一个实根，且要求有3位有效数字。

试完成：(1)估计需要二分的次数；(8分)解：容易知道方程在[1.0,1.5]有且仅有一个实根。

记此实根为*x ，根据二分法误差估计公式有12()1*22k k k b a x x ++--≤=要使得近似解有3位有效数字，只需要有22111022k -+≤⨯从而可得6k ≥，即满足精度要求的二分次数为6次。

(2)将计算过程中数据填入表1.(中间过程填写到小数点后面3位)(12分，每个k x 得2分，其它空不计分)表1题1计算过程kk a k b kx 0 1.0 1.5 1.251 1.25 1.5 1.3752 1.25 1.375 1.3163 1.313 1.375 1.3494 1.313 1.344 1.3285 1.313 1.328 1.32061.3201.3281.324二.(10分)为了计算方程()3sin 2120f x x x =--=的根，某同学将()0f x =改写为14sin 23x x =+，并建立迭代公式114sin 23k k x x +=+。

请问此迭代公式在R 上是否全局收敛的吗？说明理由。

证明：(1)对任意的x R ∈，有11113()4sin 2,333x x R ϕ⎡⎤=+∈⊆⎢⎣⎦；(2)对任意的x R ∈，有22'()cos 2133x x ϕ=≤<;从而可知，迭代格式在R 上全局收敛。

三.(20分)设有方程3()10f x x x =--=，试回答下列问题：(1)确定方程3()10f x x x =--=实根的数目；(4分)解：由2'()31f x x =-可知函数3()1f x x x =--的单调递增区间是,33⎛⎛⎫-∞-⋃+∞ ⎪ ⎪⎝⎭⎝⎭，单调递减区间是,33⎛- ⎝⎭。

§10 矩阵的正交分解(QR 分解)设nm RA ⨯∈，则存在初等反射阵s H H 1使得)1(2+=s s A A H H (上梯形)[]nmn m m n n a a a a a a a a a a a a A ,,,21212222111211=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=(按列分块) (1)(1)第1步：当01=a 时，取I H =1这一步不需约化，不妨设01≠a ，于是有初等反射阵1H 使1111e a H σ-=，其中Tu u I H 11111--=β。

于是],,,[21)1(1n Ha Ha Ha A H =⎥⎦⎤⎢⎣⎡-=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡-=222)2(121)2()2(2)2(2)2(22)2(2)2(121000D c B a a a a a a a mn m n n σσ)2(A =其中)2()1(21)2(2)2(222,),,(-⨯--∈∈=n m m T m R D Ra a c (2)第k 步：设已完成对A 上述第1步~第k-1步约化，即存在初等反射阵11,,-k H H 使)(121k k A A H H H =-其中 ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---=-)()()()(12)2(1)2(1)2(121)(k mn k mkk knk kkk n k k a a a a a a a Aσσσ ⎥⎦⎤⎢⎣⎡=k k k k k D c B r R 0其中)()1(1)()(,],,[k n k m k k n T k m k k kk k R D Ra c c -⨯+-=+-∈∈= ，为 EMBED Equation.3阶上三角阵。

如果0=k c ，这一步不需约化，取I H k =。

不妨设0≠k c ，于是存在初等反射阵kH '使 1e c H k k kσ-='计算kH '的公式： T k k k ku u I H ''-='-1β ⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎨⎧+='=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡='=+=∑)(21)()()(22)(,)(,1)(2)()(k kk k k k k k k m k k k k kk k m k i k ik k kk k a u a a a u a a sign σσβσ ………………(2) 令mm k k m k k k R H I H ⨯-+--∈⎥⎦⎤⎢⎣⎡'=111第k 步约化：)1(1)(+==k k k k A A H H A H⎥⎦⎤⎢⎣⎡''⎥⎦⎤⎢⎣⎡'=-k k k k k k kk k D H c H B r R H I 01 )1(121+-=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡'----=k k k k k k k AD H B r σσσσ方框内为第k 步约化需要计算的部分，其中)1(+k A 左上角子阵，1+k R 为k阶上三角阵，这样就使A 三角化过程前进了一步。

北航数值分析大作业第二题本页仅作为文档封面，使用时可以删除This document is for reference only-rar21year.March数值分析第二次大作业史立峰SY1505327一、 方案（1）利用循环结构将sin(0.50.2)()1.5cos( 1.2)(){i j i j ij i j i j a +≠+==(i,j=1,2,……,10)进行赋值，得到需要变换的矩阵A ；（2）然后，对矩阵A 利用Householder 矩阵进行相似变换，把A 化为上三角矩阵A (n-1)。

对A 拟上三角化，得到拟上三角矩阵A (n-1)，具体算法如下： 记A(1)=A ，并记A(r)的第r 列至第n 列的元素为()n r r j n i a r ij,,1,;,,2,1)( +==。

对于2,,2,1-=n r 执行 1. 若()n r r i a r ir,,3,2)( ++=全为零，则令A(r+1) =A(r),转5；否则转2。

2. 计算()∑+==nr i r irr a d 12)(()()r r r r r r r r r r d c a d a c ==-=++则取,0sgn )(,1)(,1若 )(,12r rr r r r a c c h +-=3. 令()nTr nrr r r r r r r r R a a c a u ∈-=++)()(,2)(,1,,,,0,,0 。

4. 计算r r T r r h u A p /)(= r r r r h u A q /)(=r r Tr r h u p t /=r r r r u t q -=ωT rr T r r r r p u u A A --=+ω)()1(5. 继续。

（3）使用带双步位移的QR 方法计算矩阵A (n-1)的全部特征值，也是A 的全部特征值，具体算法如下：1. 给定精度水平0>ε和迭代最大次数L 。

数值分析大作业（二）学院名称宇航学院专业名称航空宇航推进理论与工程学生姓名段毓学号SY16153062016年11月5日1 算法设计方案首先将矩阵A 进行拟上三角化，把矩阵A 进行QR 分解，计算出RQ 。

要得出矩阵A 的全部特征值，首先对A 进行QR 的双步位移得出特征值。

最后，采用列主元的高斯消元法求解特征向量。

1.1 A 的拟上三角化因为对矩阵进行QR 分解并不改变矩阵的结构，因此在进行QR 分解前对矩阵A 进行拟上三角化可以大大减少计算机的计算量，提高程序的运行效率。

具体算法如下所示，记A A =)1(，并记)(r A 的第r 列至第n 列的元素为()n r r j n i a r ij,,1,;,,2,1)(ΛΛ+==。

对于2,,2,1-=n r Λ执行 若()n r r i a r ir,,3,2)(Λ++=全为零，则令)()1(r r A A =+,转5；否则转2。

计算()∑+==nri r ir r a d 12)(()()r r r r r r r r r r d c a d a c ==-=++则取,0sgn )(,1)(,1若)(,12r rr r r r a c c h +-=令()nTr nrr r r r r r r r R a a c a u ∈-=++)()(,2)(,1,,,,0,,0ΛΛ。

计算r r T r r h u A p /)(=r r rr r Tr r h u p t /=r r r r u t q -=ωT rr T r r r r p u u A A --=+ω)()1(继续。

1.2 A 的QR 分解具体算法如下所示，记)1(1-=n A A ，并记[]nn r ij r a A ⨯=)(，令I Q =1 对于1,,2,1-=n r Λ执行 1.若()n r r i a r ir ,,3,1)(Λ++=全为零，则令r r Q Q =+1r r A A =+1,转5；否则转2。

一、算法设计方案：按题目要求，本程序运用带双步位移的QR方法求解给定矩阵的特征值，并对每一实特征值，求解其相应的特征向量。

总体思路：1）初始化矩阵首先需要将需要求解的矩阵输入程序。

为了防止矩阵在后面的计算中被破坏保存A[][]。

2）对给定的矩阵进行拟上三角化为了尽量减少计算量，提高程序的运行效率，在对矩阵进行QR分解之前，先进行拟上三角化。

由于矩阵的QR 分解不改变矩阵的结构，所以具有拟上三角形状的矩阵的QR分解可以减少大量的计算量。

这里用函数void QuasiTriangularization()来实现，函数形参为double型N维方阵double a[][N]。

3）对拟上三角化后的矩阵进行QR分解对拟上三角化的矩阵进行QR分解会大大减小计算量。

用子程序void QR_decomposition()来实现，将Q、R设为形参，然后将计算出来的结果传入Q和R，然后求出RQ乘积。

4）对拟上三角化后的矩阵进行带双步位移的QR分解为了加速收敛，对QR分解引入双步位移，适当选取位移量，可以避免进行复数运算。

为了进一步减少计算量，在每次进行QR分解之前，先判断是否可以直接得到矩阵的一个特征值或者通过简单的运算得到矩阵的一对特征值。

若可以，则得到特征值，同时对矩阵进行降阶处理；若不可以，则进行QR分解。

这里用函数intTwoStepDisplacement_QR()来实现。

这是用来存储计算得到的特征值的二维数组。

考虑到特征值可能为复数，因此将所有特征值均当成复数处理。

此函数中，QR分解部分用子函数void QR_decompositionMk()实现。

这里形参有三个，分别用来传递引入双步位移后的Mk阵，A矩阵，以及当前目标矩阵的维数m。

5）计算特征向量得到特征值后，计算实特征值相应的特征向量。

这里判断所得特征值的虚数部分是否为零。

求实特征值的特征向量采用求解相应的方程组((A-λI)x=0)的方法。

因此先初始化矩阵Array，计算(A-λI)，再求解方程组。

“数值分析“计算实习大作业第二题——SY1415215孔维鹏一、计算说明本程序采用带双步位移的QR方法求解矩阵A的所有特征值，然后采用反幂法求解矩阵A的实特征值对应的特征向量。

在采用带双步位移的QR方法求解特征值时，对教材上所提供的具体算法作稍微的改动，以简化程序，具体算法如下所示：1、计算出A拟上三角化后的矩阵，给定精度水平和最大迭代次数L；2、记，令k=1，m=n；3、如果，则可直接计算出最后1或2个特征值，转8，否则转4；4、如果，则可得一个特征值，置m=m-1；转3，否则转5；5、如果，则可得两个特征值，置m=m-2；转3，否则转6；6、记，计算7、k=k+1,转38、A的全部特征值已经求出，停止计算。

二、计算源程序（FORTRAN）PROGRAM SY1415215_2PARAMETER (N=10)DIMENSION A(N,N),A1(N,N),A2(N,N),C(2,N),Q(N,N),R(N,N),CR(N),CM(N)!C为存储特征值的数组,1为实部，为虚部REAL(8) A,A1,A2,C,Q,R,CME=1E-12 !精度水平L=1000 !迭代最大次数OPEN(1,FILE='数值分析大作业第二题计算结果.TXT')DO I=1,NDO J=1,NIF(I==J) THENA(I,J)=*COS(I+*J)ELSEA(I,J)=SIN*I+*J)ENDIFENDDOENDDOA1=AWRITE(*,"('矩阵A为：')")WRITE(1,"('矩阵A为：')")DO I=1,NDO J=1,NWRITE(*,"(2X,,2X,\)") A(I,J)WRITE(1,"(2X,,2X,\)") A(I,J)ENDDOWRITE(*,"(' ')")WRITE(1,"(' ')")ENDDO!使用矩阵的拟上三角化的算法将矩阵A化为拟上三角矩阵A（n-1）CALL HESSENBERG(A,N)WRITE(*,"('拟上三角化后矩阵A（n-1）为：')")WRITE(1,"('拟上三角化后矩阵A（n-1）为：')")DO I=1,NDO J=1,NWRITE(*,"(2X,,2X,\)") A(I,J)WRITE(1,"(2X,,2X,\)") A(I,J)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDO!计算对矩阵A（n-1）实行QR方法迭代结束后所得矩阵A2=ACALL QRD(A2,N,Q,R)WRITE(*,"('对矩阵A（n-1）实行QR方法迭代结束后所得Q为：')") WRITE(1,"('对矩阵A（n-1）实行QR方法迭代结束后所得Q为：')") DO I=1,NDO J=1,NWRITE(*,"(2X,,2X,\)") Q(I,J)WRITE(1,"(2X,,2X,\)") Q(I,J)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDOWRITE(*,"('对矩阵A（n-1）实行QR方法迭代结束后所得R为：')") WRITE(1,"('对矩阵A（n-1）实行QR方法迭代结束后所得R为：')") DO I=1,NDO J=1,NWRITE(*,"(2X,,2X,\)") R(I,J)WRITE(1,"(2X,,2X,\)") R(I,J)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDO!使用带双步位移的QR方法求解矩阵A（n-1）的特征值K=1M=NDO WHILE(K<=L)IF(M<=2) THENIF(M==1) THENC(1,M)=A(M,M)ELSE IF(M==2) THENCALL CALCUS(A,N,M,C)ENDIFEXITELSE IF(ABS(A(M,M-1))<E) THENC(1,M)=A(M,M)M=M-1ELSE IF(ABS(A(M-1,M-2))<E) THENCALL CALCUS(A,N,M,C)M=M-2ELSECALL CALM(A,M,N)ENDIFK=K+1ENDDOWRITE(*,"('矩阵A的全部特征值为：')")WRITE(1,"('矩阵A的全部特征值为：')")DO J=1,NWRITE(*,",'+',,'i')") C(1,J),C(2,J)WRITE(1,",'+',,'i')") C(1,J),C(2,J)ENDDO!使用反幂法求解A的相应于实特征值的特征向量J=1DO I=1,NIF(C(2,I)==0)THENCR(J)=C(1,I)J=J+1ENDIFENDDOJC=J-1WRITE(*,"('矩阵A的实特征值为：')")WRITE(1,"('矩阵A的实特征值为：')")DO I=1,JCWRITE(*,"") CR(I)WRITE(1,"") CR(I)ENDDODO II=1,JCDO I=1,NDO J=1,NIF(I==J) THENA(I,J)=A1(I,J)-CR(II)ELSEA(I,J)=A1(I,J)ENDIFENDDOENDDOCALL INPOVERMETHOD(A,N,CM)WRITE(*,"('与实特征值',,'对应的特征向量为：')") CR(II) WRITE(1,"('与实特征值',,'对应的特征向量为：')") CR(II) DO I=1,NWRITE(*,"(2X,,2X,\)") CM(I)WRITE(1,"(2X,,2X,\)") CM(I)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDOCLOSE(1)END!***************拟上三角化子函数*************************！SUBROUTINE HESSENBERG(A,N)DIMENSION A(N,N),P(N),Q(N),W(N),U(N),AT(N,N)REAL(8) A,P,Q,W,U,ATREAL(8) S0,S1,S2,S3,S4,TDO L=1,N-2JUDGE=0DO I=L+2,NIF(A(I,L)/=0) THENJUDGE=1EXITENDIFENDDOIF(JUDGE==0) THENA=ACYCLEELSE IF(JUDGE/=0) THEN!计算DRS0=0DO I=L+1,NS0=S0+A(I,L)**2ENDDODR=SQRT(S0)!计算CRIF(A(L+1,L)==0)THENCR=DRELSECR=-SGN(A(L+1,L))*DRENDIF!计算HRHR=CR**2-CR*A(L+1,L)!给u赋值IF(I<L+1) THENU(I)=0ELSE IF(I==L+1) THEN U(I)=A(I,L)-CRELSE IF(I>L+1) THEN U(I)=A(I,L)ENDIFENDDO!计算PDO I=1,NDO J=1,NAT(I,J)=A(J,I)ENDDOENDDODO I=1,NS1=0DO J=1,NS1=S1+AT(I,J)*U(J)ENDDOP(I)=S1/HRENDDO!计算QDO I=1,NS2=0DO J=1,NS2=S2+A(I,J)*U(J)ENDDOQ(I)=S2/HRENDDO!计算TS3=0DO I=1,NS3=S3+P(I)*U(I)ENDDOT=S3/HR!计算WDO I=1,NW(I)=Q(I)-T*U(I)!计算A（r+1）DO I=1,NDO J=1,NA(I,J)=A(I,J)-W(I)*U(J)-U(I)*P(J)ENDDOENDDOENDIFENDDORETURNEND!***************符号函数子程序*****************！FUNCTION SGN(X)REAL(8) XIF(X>0) THENSGN=1ELSE IF(X<0) THENSGN=-1ELSE IF(X==0) THENSGN=0ENDIFEND!*********计算二阶子阵特征值s1，s2子函数*******！SUBROUTINE CALCUS(X,N,M,Y)DIMENSION X(N,N),Y(2,N)REAL(8) A,B,C,D,X,YA=1C=X(M-1,M-1)*X(M,M)-X(M-1,M)*X(M,M-1)B=-(X(M-1,M-1)+X(M,M))D=B**2-4*CIF(D>=0) THENY(1,M)=(-B-SQRT(D))/2Y(1,M-1)=(-B+SQRT(D))/2ELSEIF(D<0) THENY(1,M)=-B/2Y(1,M-1)=-B/2Y(2,M)=-SQRT(-D)/2Y(2,M-1)=-SQRT(-D)/2ENDIFRETURNEND!*********计算Mk，Ak+1子函数************！SUBROUTINE CALM(A,M,N)DIMENSION A(N,N),MK(M,M),X(M,M),QK(M,M),RK(M,M),S1(M,M),S2(M,M),QKT(M,M) REAL(8) A,MK,X,QK,RK,QKTREAL(8) S0,S1,S2DO I=1,MDO J=1,MIF(I==J) THENX(I,J)=1ELSEX(I,J)=0ENDIFENDDOENDDOS=A(M-1,M-1)+A(M,M)T=A(M-1,M-1)*A(M,M)-A(M,M-1)*A(M-1,M)DO I=1,MDO J=1,MS0=0DO K=1,MS0=S0+A(I,K)*A(K,J)ENDDOMK(I,J)=S0-S*A(I,J)+T*X(I,J)ENDDOENDDO!对Mk做QR分解CALL QRD(MK,M,QK,RK)DO I=1,MDO J=1,MQKT(I,J)=QK(J,I)ENDDOENDDODO I=1,MDO J=1,MS1(I,J)=0DO K=1,MS1(I,J)=S1(I,J)+QKT(I,K)*A(K,J)ENDDOENDDOENDDOA=S1DO I=1,MDO J=1,MS2(I,J)=0DO K=1,MS2(I,J)=S2(I,J)+A(I,K)*QK(K,J)ENDDOENDDOENDDOA=S2RETURNEND!************QR分解子程序***************！SUBROUTINE QRD(A,N,Q,R)DIMENSION A(N,N),AT(N,N),Q(N,N),U(N),W(N),P(N),R(N,N) REAL(8) A,AT,Q,U,W,P,RREAL(8) DR,S0,CR,HR,S1,S2DO I=1,NDO J=1,NIF(I==J) THENQ(I,J)=1ELSEQ(I,J)=0ENDIFENDDOENDDODO L=1,N-1JUDGE=0DO I=L+1,NIF(A(I,L)/=0) THENJUDGE=1EXITENDIF!A(I,L)中有一个不为零,判断条件为真，跳出循环转ENDDOIF(JUDGE==0) THENQ=QA=ACYCLE!A(I,L)全为零，结束本循环，进入下一个ELSE IF(JUDGE/=0) THEN!计算DRS0=0DO I=L,NS0=S0+A(I,L)**2ENDDODR=SQRT(S0)!计算CRIF(A(L,L)==0)THENCR=DRELSECR=-SGN(A(L,L))*DRENDIF!计算HRHR=CR**2-CR*A(L,L)!给u赋值DO I=1,NIF(I<L) THENU(I)=0ELSE IF(I==L) THENU(I)=A(I,L)-CRELSE IF(I>L) THENU(I)=A(I,L)ENDIFENDDO!计算WDO I=1,NS1=0DO J=1,NS1=S1+Q(I,J)*U(J)ENDDOW(I)=S1ENDDO!计算Q（r+1）DO I=1,NDO J=1,NQ(I,J)=Q(I,J)-W(I)*U(J)/HRENDDOENDDO!计算PDO I=1,NDO J=1,NAT(I,J)=A(J,I)ENDDOENDDODO I=1,NS2=0DO J=1,NS2=S2+AT(I,J)*U(J)ENDDOP(I)=S2/HRENDDO!计算A（r+1）DO I=1,NDO J=1,NA(I,J)=A(I,J)-U(I)*P(J)ENDDOENDDOENDIFENDDOQ=QR=ARETURNEND!*************运用反幂法求解矩阵A实特征值的特征向量***********！SUBROUTINE INPOVERMETHOD(A,N,Y)DIMENSION A(N,N),U(N),Y(N),U1(N,N),L1(N,N)REAL(8) E,Z,Z1,Z2,S1,S2,BREAL(8) A,U,Y,U1,L1U(I)=1ENDDO!任取非零向量UCALL DETA(A,N,U1,L1)Z2=EIZ1=E=1K=1DO WHILE (E>1E-12)S1=0DO I=1,NS1=S1+U(I)**2ENDDOB=SQRT(S1) !1DO I=1,NY(I)=U(I)/BENDDO!2CALL DOOLITTLE(U1,L1,Y,N,U) !3利用DOOLITTLE分解法法求解Au=yS2=0DO I=1,NS2=S2+Y(I)*U(I)ENDDOZ1=Z2Z2=S2 !4E=ABS(1/Z2-1/Z1)/ABS(1/Z2) !判断是否满足精度K=K+1ENDDORETURNENDSUBROUTINE DOOLITTLE(U,L,B1,N,X)DIMENSION B(N),U(N,N),X(N),Y(N),B1(N)REAL(8) L(N,N)REAL(8) B,U,X,Y,B1REAL(8) S1,S2,S3,S4B=B1Y(1)=B(1)S3=0DO M=1,I-1S3=S3+L(I,M)*Y(M)ENDDOY(I)=B(I)-S3ENDDOX(N)=Y(N)/U(N,N)DO I=N-1,1,-1S4=0DO M=I+1,NS4=S4+U(I,M)*X(M)ENDDOX(I)=(Y(I)-S4)/U(I,I)ENDDORETURNENDSUBROUTINE DETA(A1,N,U,L) DIMENSION A(N,N),U(N,N),A1(N,N) REAL(8) L(N,N)REAL(8) X,S1,S2REAL(8) A,U,A1X=1A=A1!对矩阵A进行Doolittle分解DO K=1,NDO J=K,NS1=0DO M=1,K-1S1=S1+L(K,M)*U(M,J)ENDDOU(K,J)=A(K,J)-S1A(K,J)=U(K,J)ENDDOIF (K==N) THENEXITELSEDO I=K+1,NS2=0DO M=1,K-1S2=S2+L(I,M)*U(M,K)ENDDOL(I,K)=(A(I,K)-S2)/U(K,K)A(I,K)=L(I,K)ENDDOENDIFENDDORETURNEND三、计算结果矩阵A为：+00 +00 +00 +00 +00 +00 +00 +00 +00 +00+00 +00 +00 +00 +00 +00 +00 +00 +00 +00+00 +00 +01 +00 +00 +00 +00 +00 +00+00 +00 +00 +01 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00+00 +00 +00 +00 +01 +00 +00 +00 +00+00 +00 +00 +00 +00 +00 +01 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00+00 +00 +00 +00 +00 +00 +00 +00 +01拟上三角化后矩阵A（n-1）为：+00 +01 +00 +00 +01 +00 +00 +00+01 +01 +01 +00 +00 +01 +00 +01 +00 +00+01 +01 +01 +00 +01 +00 +00 +00 +00+01 +01 +01 +01 +00 +00 +00 +00+01 +00 +00 +00 +00 +00+00 +01 +00 +00 +00 +00+00 +00 +00+00 +00 +00 +00+00 +00 +00+00 +00对矩阵A（n-1）实行QR方法迭代结束后所得Q为：+00 +00 +00 +00 +00 +00+00 +00 +00 +00+00 +00 +00 +00+00 +00 +00 +00 +00 +00+00 +00+00 +00 +00 +00 +00+00 +00 +00 +00+00+00 +00 +00+00 +00对矩阵A（n-1）实行QR方法迭代结束后所得R为：+01 +01 +01 +00 +01 +00 +00 +00 +00 +01 +00 +00 +00 +00 +00 +00 +00 +00 +01 +01 +00 +01 +00 +01 +00 +00 +00 +00 +01 +00 +00 +00 +00+00 +00 +01 +01 +00 +00 +00+00 +00 +01 +00 +00 +00 +00+00 +00 +00 +00 +00 +00 +00+00 +00 +00 +00 +01 +00 +00+00 +00 +00 +00 +00 +00+00 +00 +00 +00矩阵A的全部特征值为：+01+ +00i+01++00i+01++00i+01+ +00i+01+ +00i+00++00i+00++00i+00+ +00i+00+ +00i+ +00i矩阵A的实特征值为：+01+01+01+00+00与实特征值 +01对应的特征向量为：+00 +00 +00 +00 +00 +00 +00 +00 +00与实特征值 +01对应的特征向量为：+00 +00 +00 +00 +00 +00 +00与实特征值+01对应的特征向量为：+00 +00 +00与实特征值 +00对应的特征向量为：+00 +00 +00 +00 +00与实特征值 +00对应的特征向量为：+00 +00 +00 +00 +00 +00 +00与实特征值对应的特征向量为：+00 +00 +00 +00 +00 +00 +00 +00。

QR分解、SVD分解、最小二乘问题数值上机报告练习6.13 随机构造一个可逆方阵，利用不同的方法给出它的QR分解，观察所得列向量的正交性、CPU时间和所谓的向后稳定性。

解：本题中所得列向量的正交性由||Q T Q−I||2来进行刻画。

向后稳定性由||A−QR||2来进行刻画。

实验步骤：1、首先分别编写出CGS、MGS、Householder、Givens 等不同方法进行矩阵的QR分解的程序。

2、在主程序中调用这些子程序。

其中在编写利用Householder矩阵进行矩阵的QR分解时，可以利用H矩阵作为秩一修正矩阵在计算复杂度方面的优势，即 H∗g=g−b−1(u T g)u,只需2n+1次乘除法运算。

实验数据：（1）取实验矩阵A为满秩的200*200的随机矩阵时得到的数据如下：（2）取实验矩阵A为满秩的400*400的随机矩阵时得到的数据如下：实验数据分析：1.从向后稳定性上来看，CGS和MGS给出的QR乘积均能很好地近似原始矩阵A，因此都是向后稳定的算法。

而Householder方法和Givens方法的向后稳定性略差。

2.从所求得列向量的正交性来看，CGS方法最差，MGS比CGS稍好，Householder方法和Givens方法的表现最好。

Householder方法给出的Q具有更好的列正交性。

3.从CPU时间上来看，Givens方法所耗时间最长，这是合理的，因为对应于一个Householder矩阵，Givens方法中要构造一系列Givens矩阵。

理论分析，对于稠密阵，Householder方法所需的乘除法次数是Givens方法的一半。

且由讲义所指出，有如下结果：N opt House≈∑2(n+1−k)(m+1−k)≈mn2−13n3nk=1N opt GS ≈∑2(k −1)m ≈mn 2n k=1.而由Householder 镜像变换法的计算复杂度略低于MGS 方法，Householder 方法所耗的CPU 时间略低于MGS 方法，故是合理的。