パソコンで数値計算

実対称行列の対角化

実対称行列の性質

実対称行列とは、行列A の転置行列をとすると
   (1.1)
を満たす行列のことです。実対称行列は、ある直行行列で対角化可能で、固有値は必ず実数となる性質を持っています。

固有値問題

物理学の問題において、固有値問題は数多く出てきますが、そのとき現れる方程式は
   (2.1)
でありこれは固有方程式と呼ばれます。ここで、Aは任意の数nのn×n次行列で、λは固有値、Xは列ベクトルで固有ベクトルといいます。

もしnが大きい値となると人の手による計算でこの方程式を解くのは大変であり、また、扱う物理学の問題によってはnの値が無限となる場合もあります。そこでこのような問題はコンピュータによって解く方法がとられます。

解く方法において行列Aを対角化して解く方法があります。その方法には様々な形式がありますが、ここでは最も簡単な方法で解く方法を紹介します。

実対称行列を対角化する手段

ある直行行列をUとすると、n次実対称行列Aは必ず次の計算
   (3.1)
により対角化されます。しかし、実対称行列Aからすぐに対角化する直行行列Uを見つけるのは難しいです。そこでまず実対称行列の非対角成分のうち最もその絶対値が大きい成分に注目します。その成分がであったとします。そして次に第i行i列と第j行j列の値がcosθであとの対角成分は全て1であり、第i行j列の値がsinθで第j行i列の値が-sinθで、残りの成分は全て0である直行行列を考えます。

そして、式(3.1)の計算を行い計算結果の行列をCとすると各成分は
   (3.2)
となります。次にとなるようにθの値を決めます。場合は
   (3.3)
となります。の場合、θの値は何通りにもとれますが、ここではπ/4とします。

これで行列Cの全ての成分の値が決まります。そして今度は行列Cに同様な計算を行います。以降これらの計算を繰り返し最終的に実対称行列の全ての非対角成分がある値より小さくなったら計算を終わりにします。これは、コンピュータでは全ての非対角行列を0にはできない場合があるからです。よって完全に0になるまで計算を行うとすると、無限ループになってしまう危険性があります。

プログラムソース

このプログラムではある適当に決めた3次実対称行列の対角化を行います。あらかじめ実対称行列の値を配列aに格納するようにしていますが、プログラムを実行するたびに任意の実対称行列を対角化したい場合は、実対称行列が書かれたファイルから読み取るようにする等とプログラムを改変すると良いでしょう。



#include <stdio.h>
#include <math.h>
 
#define N 3                     //実対象行列の次数
#define PI 3.1415926535897932385                //円周率
 
/*実対称行列とそれを対角化する直行行列を
格納する配列を受け取り対角化を行う関数。*/
void taikakuka(double a[N][N],double U[N][N]);  
 
int main()
{
        int i,j;
//配列aには実対称行列を入れる。また結果の直行行列を入れる配列Uを作っておく。
        double U[N][N],a[N][N]={{0.0,1.0,-1.0},{1.0,0.0,1.0},{-1.0,1.0,0.0}};
 
        taikakuka(a,U);         //対角化関数の呼び出し
 
//計算結果の対角行列の表示
        printf("対角行列は\n");
        for(i=0;i<N;i++) {
                for(j=0;j<N;j++) {
                        printf("%f ",a[i][j]);
                }
                printf("\n");
        }
//計算結果の直交行列の表示
        printf("直行行列は\n");
        for(i=0;i<N;i++) {
                for(j=0;j<N;j++) {
                        printf("%f ",U[i][j]);
                }
                printf("\n");
        }
 
        return 0;
}
 
void taikakuka(double a[N][N],double U[N][N])
{
        int i,j,l,m,p,q,count;
        double max,theta;
        double u[N][N],oldU[N][N],newa[N][N];
 
        FILE *output;
        output=fopen("output.dat","w");
 
/*計算結果としてだされた直行行列を
格納するための配列を単位行列に初期化しておく。*/
        for(p=0;p<N;p++) {
                for(q=0;q<N;q++) {
                        if(p==q) U[p][q]=1.0;
                        else U[p][q]=0.0;
                }
        }
 
        for(count=0;count<=10000;count++) {
 
/*配列olduは新たな対角化計算を行う前に
かけてきた直行行列を保持する。*/
                for(p=0;p<N;p++) {
                        for(q=0;q<N;q++) {
                                oldU[p][q]=U[p][q];
/*対角化するときの個々の直交行列を入れる配列uの初期化。(単位行列にする。)*/
                                if(p==q) u[p][q]=1.0;
                                else u[p][q]=0.0;
                        }
                }
 
/*非対角要素の中から絶対値の最大のものを見つける*/
                max=0.0;
                for(p=0;p<N;p++) {
                        for(q=0;q<N;q++) {
                                if(max<fabs(a[p][q]) && p!=q) {
                                        max=fabs(a[p][q]);
/*その最大のものの成分の行と列にあたる数を記憶しておく。*/
                                        i=p;
                                        j=q;
                                }
                        }
                }
/*先ほど選んだ最大のものが指定の値より小さければ対角化終了*/
                if(fabs(a[i][j]) < 1.0e-10) {
                        break;
                }
/*条件によってシータの値を決める*/
                if(a[i][i]==a[j][j]) theta=PI/4.0;
                else theta=atan(-2*a[i][j]/(a[i][i]-a[j][j]))/2.0;
 
/*ここでこのときに実対称行列にかける個々の直行行列uが決まるが
特にここでの計算の意味はない。(する必要はない。)*/
                u[i][j]=sin(theta);
                u[j][i]=-sin(theta);
                u[i][i]=u[j][j]=cos(theta);
 
/*ここでいままで実対称行列にかけてきた直行行列を配列Uに入れる。*/
                for(p=0;p<N;p++) {
                        U[p][i]=oldU[p][i]*cos(theta)-oldU[p][j]*sin(theta);
                        U[p][j]=oldU[p][i]*sin(theta)+oldU[p][j]*cos(theta);
                }
 
//対角化計算によってでた新たな実対称行列の成分を配列newaに入れる。
                newa[i][i]=a[i][i]*cos(theta)*cos(theta)
                        +a[j][j]*sin(theta)*sin(theta)
                        -2.0*a[i][j]*sin(theta)*cos(theta);
                newa[j][j]=a[i][i]*sin(theta)*sin(theta)
                        +a[j][j]*cos(theta)*cos(theta)
                        +2.0*a[i][j]*sin(theta)*cos(theta);
                newa[i][j]=newa[j][i]=0.0;
                for(l=0;l<N;l++) {
                        if(l!=i && l!=j) {
                                newa[i][l]=a[i][l]*cos(theta)
                                    -a[j][l]*sin(theta);
                                newa[l][i]=newa[i][l];
                                newa[j][l]=a[i][l]*sin(theta)
                                    +a[j][l]*cos(theta);
                                newa[l][j]=newa[j][l];
                        }
                }
                for(l=0;l<N;l++) {
                        for(m=0;m<N;m++) {
                            if(l!=i && l!=j && m!=i && m!=j) {
                                    newa[l][m]=a[l][m];
                            }
                        }
                }
 
//次の対角化計算を行う行列の成分を配列aへ上書きする。
                for(p=0;p<N;p++) {
                        for(q=0;q<N;q++) {
                                a[p][q]=newa[p][q];
                        }
                }
 
                if(count==10000) {
                        printf("対角化するために作業を繰り返す必要がまだあり.\n");
                }
 
//対角化の途中計算をファイルに記録する。
                fprintf(output,"選択した非対角要素は %d 行 %d 列です。\n",i+1,j+1);
                fprintf(output,"このときのシータの値は %f です。\n",theta);
                fprintf(output,"また計算後の実対称行列は\n");
                for(i=0;i<N;i++) {
                        for(j=0;j<N;j++) {
                                fprintf(output,"%f ",newa[i][j]);
                        }
                        fprintf(output,"\n");
                }
                fprintf(output,"\nこのときまでに計算してきた直行行列は\n");
                for(i=0;i<N;i++) {
                        for(j=0;j<N;j++) {
                                fprintf(output,"%f ",U[i][j]);
                        }
                        fprintf(output,"\n");
                }
                fprintf(output,"\n\n");
        }
        
        fclose(output);
}
 


数値計算実行結果

上記のプログラムを実行すると以下のように表示されます。これから固有値は縮退ありの-2と1である事がわかります。


対角行列は
-2.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000
直行行列は
0.577350 0.707107 -0.408248
-0.577350 0.707107 0.408248
0.577350 0.000000 0.816497

2009/04/19 更新