解析

• ヤコビ法による対称行列の固有値問題

　特に物理現象において、対称行列の固有値、固有ベクトルを求める場面が多く、対称行列の場合のみ、ヤコビ法により簡単に固有値問題を解くことができます。非対称行列の場合は「べき乗法」という解法がありますが、絶対値最大の固有値1個しか求められません。

行列、対称行列の実装についてはこちらを参照してください。
行列クラス
対称行列クラス

		
		
SquareMatrix a = new SquareMatrix( 3 );
a.SetRow( 0, 1, 0.446, -0.56 );
a.SetRow( 1, 0.446, 1, -0.239 );
a.SetRow( 2, -0.56, -0.239, 1 );

Vector eigenValue;
List<Vector> eigenVectors;
eigenVectors = GetEigenVectorsByJacobi( a, 100000, 0.00000001, out eigenValue );



/// <summary>正方行列の固有値、固有ベクトルを算出（ヤコビ法・対称行列のみ有効）</summary>
/// <param name="matrix">対称行列</param>
/// <param name="nMaxStep">最大評価ステップ数</param>
/// <param name="dSettleValue">収束値</param>
/// <param name="eigenValue">固有値リスト</param>
/// <returns>固有ベクトルのリスト</returns>
/// <remarks>
/// 非対称行列の場合は「べき乗法」という解法があるが、絶対値最大の固有値しか求まらない。
/// http://www.asahi-net.or.jp/~UC3K-YMD/Lesson/Section03/section03_09.html
/// </remarks>
public static List<Vector> GetEigenVectorsByJacobi( SquareMatrix matrix, int nMaxStep, double dSettleValue, out Vector eigenValue )
{
　　eigenValue = null;
　　SquareMatrix A1 = ( SquareMatrix )matrix.Clone();
　　SquareMatrix X1 = new SquareMatrix( matrix.Size );
　　X1.SetUnitMatrix();

　　SquareMatrix A2 = new SquareMatrix( matrix.Size );
　　SquareMatrix X2 = new SquareMatrix( matrix.Size );
　　int nStep = 0, p = 0, q = 0;

　　while( nStep < nMaxStep )
　　{
　　　　// 最大要素の収束判定
　　　　double dMaxElement = A1.GetMaxElement( out p, out q );
　　　　if( dMaxElement < dSettleValue )
　　　　{
　　　　　　// 収束した
　　　　　　break;
　　　　}

　　　　double t = 0.5 * ( A1[ p, p ] - A1[ q, q ] );
　　　　double v = Math.Abs( t ) / Math.Sqrt( A1[ p, q ] * A1[ p, q ] + t * t );
　　　　double sn = Math.Sqrt( 0.5 * ( 1.0 - v ) );
　　　　if( -A1[ p, q ] * t < 0.0 )
　　　　{
　　　　　　sn = -sn;
　　　　}
　　　　double cs = Math.Sqrt( 1.0 - sn * sn );

　　　　// Akの計算
　　　　for( int i = 0; i < matrix.Size; i++ )
　　　　{
　　　　　　if( i == p )
　　　　　　{
　　　　　　　　for( int j = 0; j < matrix.Size; j++ )
　　　　　　　　{
　　　　　　　　　　if( j == p )
　　　　　　　　　　　　A2[ p, p ] = A1[ p, p ] * cs * cs + A1[ q, q ] * sn * sn -
　　　　　　　　　　　　　　　　   2.0 * A1[ p, q ] * sn * cs;
　　　　　　　　　　else if( j == q )
　　　　　　　　　　　　A2[ p, q ] = 0.0;
　　　　　　　　　　else
　　　　　　　　　　　　A2[ p, j ] = A1[ p, j ] * cs - A1[ q, j ] * sn;
　　　　　　　　}
　　　　　　}
　　　　　　else if( i == q )
　　　　　　{
　　　　　　　　for( int j = 0; j < matrix.Size; j++ )
　　　　　　　　{
　　　　　　　　　　if( j == q )
　　　　　　　　　　　　A2[ q, q ] = A1[ p, p ] * sn * sn + A1[ q, q ] * cs * cs +
　　　　　　　　　　　　　　　　   2.0 * A1[ p, q ] * sn * cs;
　　　　　　　　　　else if( j == p )
　　　　　　　　　　　　A2[ q, p ] = 0.0;
　　　　　　　　　　else
　　　　　　　　　　　　A2[ q, j ] = A1[ q, j ] * cs + A1[ p, j ] * sn;
　　　　　　　　}
　　　　　　}
　　　　　　else
　　　　　　{
　　　　　　　　for( int j = 0; j < matrix.Size; j++ )
　　　　　　　　{
　　　　　　　　　　if( j == p )
　　　　　　　　　　　　A2[ i, p ] = A1[ i, p ] * cs - A1[ i, q ] * sn;
　　　　　　　　　　else if( j == q )
　　　　　　　　　　　　A2[ i, q ] = A1[ i, q ] * cs + A1[ i, p ] * sn;
　　　　　　　　　　else
　　　　　　　　　　　　A2[ i, j ] = A1[ i, j ];
　　　　　　　　}
　　　　　　}
　　　　}
　　　　// Xkの計算
　　　　for( int i = 0; i < matrix.Size; i++ )
　　　　{
　　　　　　for( int j = 0; j < matrix.Size; j++ )
　　　　　　{
　　　　　　　　if( j == p )
　　　　　　　　　　X2[ i, p ] = X1[ i, p ] * cs - X1[ i, q ] * sn;
　　　　　　　　else if( j == q )
　　　　　　　　　　X2[ i, q ] = X1[ i, q ] * cs + X1[ i, p ] * sn;
　　　　　　　　else
　　　　　　　　　　X2[ i, j ] = X1[ i, j ];
　　　　　　}
　　　　}

　　　　A1 = ( SquareMatrix )A2.Clone();
　　　　X1 = ( SquareMatrix )X2.Clone();

　　　　// 次のステップへ
　　　　nStep++;
　　}

　　if( nStep >= nMaxStep ) return null;

　　// 行列A1,X1に固有値と固有ベクトルが設定された。

　　// 固有ベクトルリスト
　　List<Vector> eigenVectors = new List<Vector>();

　　// 固有値リスト
　　eigenValue = new Vector( matrix.Size );
　　for( int i = 0; i < matrix.Size; i++ )
　　{
　　　　// 固有値
　　　　eigenValue[ i ] = A1[ i, i ];

　　　　// 固有ベクトル
　　　　Vector eigenVector = new Vector( matrix.Size );
　　　　for( int j = 0; j < matrix.Size; j++ )
　　　　{
　　　　　　eigenVector[ j ] = X1[ j, i ];
　　　　}
　　　　eigenVectors.Add( eigenVector );
　　}

　　return eigenVectors;
}