北航数值分析大作业-第二题-QR分解.docx

数值分析第二题梁进明SY0906529算法设计方案。一 矩阵的QR分解。把矩阵A分解为一个正交矩阵Q与一个上三角矩阵R的乘积，称为矩阵A的正三角分解，简称QR分解。QR分解的算法如下：记，并记，令（n阶单位矩阵）对于r=1,2,n-1执行(1) 若全为零，则令，转(5);否则转(2)(2) 计算(若,则取)(3) 令(4) 计算(5) 继续 当此算法执行完后就得到正交矩阵和上三角矩阵且有。二 矩阵的拟上三角化。对实矩阵A的拟上三角化具体算法如下：记，并记的第r列到第n列的元素为。对于执行（1） 若全为零，则令,转(5);否则转(2)。（2） 计算（3） 令（4） 计算（5） 继续算法执行完后，就得到与原矩阵A相似的拟上三角矩阵。三 带双步位移的QR方法求实矩阵全部特征值的具体算法如下：（1） 使用矩阵的拟上三角化的算法把矩阵化为拟上三角矩阵；给定精度水平和迭化最大次数L。（2） 记，令。（3） 如果，则得到A的一个特征值，置（降阶），转（4）；否则转（5）。（4） 如果,则得到A的一个特征值，转(11)；如果,则直接转(11);如果转（3）。（5） 求二阶子阵的两个特征值和，即计算二次方程的两个根和（6） 如果,则得到A的两个特征值和,转(11)；否则转(7)。（7） 如果,则得到A的两个特征值和，置(降阶),转(4);否则转(8)。（8） 如果k=L,则计算终止，未得到A的全部特征值；否则转9。（9） 记,计算（10） 置，转（3）。（11） A的全部特征值已计算完毕，停止计算。四 算法的主过程（1） 根据公式生成原始矩阵（2） 对矩阵A进行拟上三角化（3） 用带双步位移的QR方法求矩阵的特征值（4） 对每个特征值 ，解出线性方程组的所有非零解即得A 的属于的全部特征向量。全部源代码#include#includeusing namespace std;const int n = 10;const double precision = 1e-12;double a n n ;FILE *file;struct Eigenvalue double r;double i;Eigenvalue eigenvalue n ;/生成矩阵Avoid getMatrix() for( int i = 0; i n; i+ ) for( int j = 0; j n; j+ ) if( i != j ) a i j = sin( 0.5 * ( i + 1 ) + 0.2 * ( j + 1 ) ); else a i j = 1.5 * cos( ( i + 1 ) + 1.2 * ( j + 1 ) );/设单位矩阵void setUnitMatrix( double u n n ) for( int i = 0; i n; i+ ) for( int j = 0; j n; j+ ) if( i = j ) u i j = 1.0; else u i j = 0.0;/矩阵乘法void matrixMultiply( double ma n n , double mb n n , double mc n n ) double res n n ;int i, j, k;for( i = 0; i n; i+ ) for( j = 0; j n; j+ ) res i j = 0.0;for( k = 0; k n; k+ ) res i j += ma i k * mb k j ;for( i = 0; i n; i+ ) for( j = 0; j n; j+ ) mc i j = res i j ;/得到转置矩阵void getTransposedMatrix( double sMatrix n n , double dMatrix n n ) for( int i = 0; i n; i+ ) for( int j = 0; j n; j+ ) dMatrix i j = sMatrix j i ;/拟上三角化void triangulation() int i, j , k;double s n , u n , e n , h n n ;double sum, c;for( i = 0; i n - 2; i+ ) for( j = 0; j i ) s j = a j i ; else s j = 0;sum = 0.0;for( j = 0; j 0 ? -sqrt( sum ) : sqrt( sum );memset( e, 0, sizeof( e ) );e i + 1 = 1;sum = 0.0;for( j = 0; j precision ) for( j = 0; j n; j+ ) for( k = 0; k n; k+ ) if( j = k ) h j k = 1.0 - 2.0 *u j * u k / sum ; else h j k = - 2.0 * u j * u k / sum; else setUnitMatrix( h );matrixMultiply( h, a, a );matrixMultiply( a, h, a );/QR 分解void qrDecomposition( double aMatrix n n , double qMatrix n n , double rMatrix n n ) int i, j , k;double s n , u n , e n , h n n ;double sum, c;setUnitMatrix( qMatrix );for( i = 0; i n - 1; i+ ) for( j = 0; j = i ) s j = aMatrix j i ; else s j = 0;sum = 0.0;for( j = 0; j 0 ? -sqrt( sum ) : sqrt( sum );for( j = 0; j n; j+ ) e j = 0.0;e i = 1;sum = 0.0;for( j = 0; j precision ) for( j = 0; j n; j+ ) for( k = 0; k n; k+ ) if( j = k ) h j k = 1 - 2.0 *u j * u k / sum ; else h j k = - 2.0 * u j * u k / sum; else setUnitMatrix( h );matrixMultiply( qMatrix, h, qMatrix );matrixMultiply( h, aMatrix, aMatrix );for( i = 0; i n; i+ ) for( j = 0; j n; j+ ) rMatrix i j = aMatrix i j ;/解二元一次方程组void getEigenvalue( int m, Eigenvalue &ea , Eigenvalue &eb ) double A = 1.0, B = -( a m - 1 m - 1 + a m m ), C = a m - 1 m - 1 * a m m - a m m - 1 * a m - 1 m ;double delta = B * B - 4 * A * C;if( fabs( delta ) 0 ) ea.r = ( -B + sqrt( delta ) ) / ( 2.0 * A );ea.i = 0.0;eb.r = ( -B - sqrt( delta ) ) / ( 2.0 * A );eb.i = 0.0; else ea.r = -B / ( 2.0 * A );ea.i = sqrt( -delta ) / ( 2.0 * A );eb.r = -B / ( 2.0 * A );eb.i = -sqrt( -delta ) / ( 2.0 * A );/带双步位移的QR方法void qrAlgorithm() double mMatrix n n , qMatrix n n , rMatrix n n , tMatrix n n ;double dMatrix 2 2 ;double s, t;int i, j, k;int m = n - 1;while( true ) if( m 0 ) break;if( m = 0 ) eigenvalue m .r = a m m ;eigenvalue m .i = 0.0;break;if( fabs( a m m - 1 ) precision ) eigenvalue m .r = a m m ;eigenvalue m .i = 0.0;m -= 1;continue;if( m = 1 ) getEigenvalue( m, eigenvalue m-1 , eigenvalue m );break;if( fabs( a m - 1 m - 2 ) precision ) getEigenvalue( m, eigenvalue m-1 , eigenvalue m );m -= 2;continue;s = 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 ;for( i = 0; i n; i+ ) for( j = 0; j n; j+ ) mMatrix i j = 0.0;for( k = 0; k n; k+ ) mMatrix i j += a i k * a k j ;for( i = 0; i n; i+ ) for( j = 0; j n; j+ ) mMatrix i j -= s * a i j ;for( i = 0; i m; i+ ) mMatrix i i += t;qrDecomposition( mMatrix, qMatrix, rMatrix );getTransposedMatrix( qMatrix, tMatrix );matrixMultiply( tMatrix, a, a );matrixMultiply( a, qMatrix, a );/高斯消元法templatebool gaussianElimination( double (&a) n n , double (&x) n , double (&b) n ) int i, j, k, maxline;double maxUnit, m;for( i = 0; i n; i+ ) maxUnit = fabs( a i i ), maxline = i;for( j = i + 1; j maxUnit ) maxUnit = fabs( a j i );maxline = j;for( j = 0; j n; j+ ) swap( a maxline j , a i j );swap( b maxline , b i );if( fabs( a i i ) precision ) continue;for( j = i + 1; j n; j+ ) m = a j i / a i i ;for( k = i; k n; k+ ) a j k -= m * a i k ;b j -= m * b i ;bool flag = false;double solution 10 10 ;for( i = 0; i n; i+ ) for( j = 0; j n; j+ ) solution i j = 0.0;for( i = 0; i n; i+ ) if( fabs( a i i ) = 0; i- ) if( fabs( a i i ) precision ) for( j = 0; j n; j+ ) for( k = i + 1; k n; k+ ) solution i j -= a i j * solution k j ;for( j = 0; j n; j+ ) solution i j /= a i i ;for( i = 0; i n; i+ ) flag = false;for( j = 0; j precision ) flag = true;if( flag = true ) fprintf( file, );for( j = 0; j n; j+ ) fprintf( file, %0.12e , , solution j i );fprintf( file, * x%dn, i );return true;int main() double c 10 10 , b 10 , x 10 , tmp 10 10 ;file = fopen( second.txt, w );getMatrix();for( int i = 0; i n; i+ ) for( int j = 0; j n; j+ ) tmp i j = a i j ;triangulation(); for( int i = 0; i n; i+ ) for( int j = 0; j n; j+ ) fprintf( file, %0.12e , a i j );fprintf( file, n );qrAlgorithm();for( int i = 0; i n; i+ ) fprintf( file, %d %0.12e %0.12en, i, eigenvalue i .r, eigenvalue i .i );if( fabs( eigenvalue i .i ) precision ) for( int j = 0; j n; j+ ) b j = 0;for( int k = 0; k n; k+ ) if( j = k ) c j k = tmp j k - eigenvalue i .r; else c j k = tmp j k ;gaussianElimination( c, x, b );fflush( file );fclose( file );return 0;拟上三角化后得到的矩阵Aijaijijaij11-8.827516758830e-00112-9.933136491826e-00213-1.103349285994e+00014-7.600443585637e-001151.549101079914e-00116-1.946591862872e+00017-8.782436382927e-00218-9.255889387184e-001196.032599440534e-0011101.518860956469e-00121-2.347878362416e+000222.372370104937e+000231.819290822208e+000243.237804101546e-001252.205798440320e-001262.102692662546e+000271.816138086098e-001281.278839089990e+00029-6.380578124405e-001210-4.154075603804e-00131-2.409448106601e-016321.728274599968e+00033-1.171467642785e+00034-1.243839262699e+00035-6.399758341743e-00136-2.002833079037e+000372.924947206124e-00138-6.412830068395e-001399.783997621285e-0023102.557763574160e-00141-1.137050343952e-016426.849619556674e-01843-1.291669534130e+00044-1.111603513396e+000451.171346824096e+00046-1.307356030021e+000471.803699177750e-00148-4.246385358369e-001497.988955239304e-0024101.608819928069e-00151-3.218079046866e-017527.361962688705e-01753-9.111415398526e-017541.560126298527e+000558.125049397515e-001564.421756832923e-00157-3.588616128137e-002584.691742313671e-00159-2.736595050092e-001510-7.359334657750e-002611.645514213628e-016623.986588458133e-017631.795101023347e-016641.162374785919e-01765-7.707773755194e-00166-1.583051425742e+00067-3.042843176799e-001682.528712446035e-00169-6.709925401449e-0016102.544619929082e-001711.771418473192e-01672-2.213674804097e-01673-3.644564058372e-017746.835111609609e-017751.748248360870e-01776-7.463453456938e-00177-2.708365157019e-00278-9.486521893682e-001791.195871081495e-0017101.929265617952e-002813.004355963132e-016821.263510892920e-01683-2.122437351083e-01684-1.393244144293e-017851.103277909913e-017867.956540121416e-01787-7.701801374364e-00188-4.697623990618e-001894.988259468008e-0018101.137691603776e-001911.733070882891e-01792-6.859192396325e-017931.431139517732e-01694-2.959918849953e-01795-1.679690935074e-017966.936176716053e-01897-7.498270272400e-017987.013167092107e-001991.582180688475e-0019103.862594614233e-0011011.166891323731e-016102-1.182725878649e-016103-2.943370580897e-0171041.609893155796e-0161055.090843430892e-017106-9.377222276666e-019107-4.613183027682e-0171080.000000000000e+0001094.843807602783e-00110103.992777995177e-001QR分解方法后得到的矩阵ijaijijaij113.383039617436e+00012-6.549939826982e-00113-1.083524130978e+00014-8.466149434806e-002152.612724951410e-00116-1.037221495828e+000171.601955828153e+000181.204278823958e-001191.014804983176e+0001104.960407023480e-00121-2.910244719154e-02022-2.607230325699e+000232.342162457137e+000246.088185449434e-002254.982796671050e-002261.990753512493e+00027-1.099309977681e+000281.667218318714e-00129-1.893430300111e+000210-9.618234274457e-00131-7.239072923059e-02132-3.748785281332e-00133-2.039762094724e+000346.355514251060e-001351.295757273877e-00236-1.228105975101e+000372.939582953174e-00138-2.703389028865e-001398.759420577424e-0013106.887037680994e-00241-4.956319424587e-022421.697858876186e-01243-5.868046856648e-013441.577548559723e+000451.396154902502e-00246-5.683586669215e-001474.324633432180e-00148-2.037617784376e-002495.122305756965e-0014102.736674654195e-001519.296461122248e-02952-3.188315445277e-019531.107706655675e-01954-5.724228551970e-00755-1.484039824869e+000567.208773358910e-00257-7.942613126537e-00258-7.985487831331e-00359-3.384477335678e-002510-4.414179642373e-002612.773836628398e-03062-1.280643245291e-020639.519575558055e-021647.313838064850e-018651.175778294417e-01966-9.272351359411e-00167-6.846334354260e-00168-2.175327429496e-00169-2.433999342340e-001610-4.296654261792e-001712.421176050913e-03072-1.070413382049e-020737.043491179201e-021743.274954414022e-018751.344270057474e-019762.311426581385e-00277-1.033826776702e+000781.642245514305e-00179-3.654673010054e-001710-8.923018276994e-00281-9.423424275349e-030824.144194657899e-02083-2.687848528422e-02084-1.245736700019e-01785-5.085013679545e-01986-1.825709393012e-010874.426935739629e-010889.355889077604e-00189-1.876929916381e-001810-1.360314863118e-00191-1.504471484290e-029926.615591745028e-02093-4.290541157588e-02094-1.989442023600e-01795-8.118899381244e-01996-2.927314845537e-010977.096068621591e-01098-9.364060635126e-011996.360627876959e-0019102.736788651917e-0011019.853610797099e-035102-5.218295742846e-0251034.557627389873e-0251041.618363931390e-0241054.754782347687e-0261068.662762125027e-022107-2.181749650711e-0211082.704236690829e-022109-5.021431940651e-01710105.650488993501e-002矩阵A的全部特征值iRiIi13.383039617436e+0000.000000000000e+0002-2.323496210212e+0008.930405177196e-0013-2.323496210212e+000-8.930405177196e-00141.577548557113e+0000.000000000000e+0005-1.484039822259e+0000.000000000000e+0006-9.805309558452e-0011.139489139816e-0017-9.805309558452e-001-1.139489139816e-00189.355889078188e-0010.000000000000e+00096.360627875745e-0010.000000000000e+000105.650488993501e-0020.000000000000e+000A对应于实特征值的特征向量特征值