北航數(shù)值分析大作業(yè)-第二題-QR分解.docx

數(shù)值分析第二題梁進明SY0906529算法設(shè)計方案。一 矩陣的QR分解。把矩陣A分解為一個正交矩陣Q與一個上三角矩陣R的乘積，稱為矩陣A的正三角分解，簡稱QR分解。QR分解的算法如下：記，并記，令（n階單位矩陣）對于r=1,2,n-1執(zhí)行(1) 若全為零，則令，轉(zhuǎn)(5);否則轉(zhuǎn)(2)(2) 計算(若,則取)(3) 令(4) 計算(5) 繼續(xù) 當此算法執(zhí)行完后就得到正交矩陣和上三角矩陣且有。二 矩陣的擬上三角化。對實矩陣A的擬上三角化具體算法如下：記，并記的第r列到第n列的元素為。對于執(zhí)行（1） 若全為零，則令,轉(zhuǎn)(5);否則轉(zhuǎn)(2)。（2） 計算（3） 令（4） 計算（5） 繼續(xù)算法執(zhí)行完后，就得到與原矩陣A相似的擬上三角矩陣。三 帶雙步位移的QR方法求實矩陣全部特征值的具體算法如下：（1） 使用矩陣的擬上三角化的算法把矩陣化為擬上三角矩陣；給定精度水平和迭化最大次數(shù)L。（2） 記，令。（3） 如果，則得到A的一個特征值，置（降階），轉(zhuǎn)（4）；否則轉(zhuǎn)（5）。（4） 如果,則得到A的一個特征值，轉(zhuǎn)(11)；如果,則直接轉(zhuǎn)(11);如果轉(zhuǎn)（3）。（5） 求二階子陣的兩個特征值和，即計算二次方程的兩個根和（6） 如果,則得到A的兩個特征值和,轉(zhuǎn)(11)；否則轉(zhuǎn)(7)。（7） 如果,則得到A的兩個特征值和，置(降階),轉(zhuǎn)(4);否則轉(zhuǎn)(8)。（8） 如果k=L,則計算終止，未得到A的全部特征值；否則轉(zhuǎn)9。（9） 記,計算（10） 置，轉(zhuǎn)（3）。（11） A的全部特征值已計算完畢，停止計算。四 算法的主過程（1） 根據(jù)公式生成原始矩陣（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 ) );/設(shè)單位矩陣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 ;/得到轉(zhuǎn)置矩陣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對應(yīng)于實特征值的特征向量特征值