2002年2月14日 htmlファイル作成
関谷トップページへ
数値計算2001s トップページへ
数値計算2001s No.18 実対称行列の固有値･固有ベクトル解法ページへ
作者へのメッセージ
数値計算2001s No.18 実対称行列の固有値･固有ベクトル解法 ヤコビ法プログラム
/* jacobipq.c program for jacobiの実対称行列の固有値･固有ベクトル解法
応用数値計算法入門 pp.80-89 1996.9.26 z.sekiya( k134084 関谷順太)
1997.10.27 正方行列に限定して関数名の変更、p'ap の検算の追加
1997.11.4 戸川の式(6.19,6.20)により、対角化の計算を行う。 */
/* プログラム仕様
1. インプット
１） データのタイトル 半角で49文字まで、(例："１．９．３項"）
２） n:実対称行列の次数
３） 実対称行列のａの要素（ｎ行数 × ｎ列数個）
行列の添え字は、１からとし、実数データとする
２．処理
１） パラメータデータを読み込む
２） 実対称行列ａを読み込み、その表示（チェックﾌﾟﾘﾝﾄ）を行う
３） ヤコビ法で、固有値･固有ベクトルを解く。
４） 収束すれば、解，検算結果１）ｐａｐ、２）逆行列ｐｉを出力する。
５） 収束しないときは収束しないと出力する。
３．アウトプット
１） ヤコビ法とデータのタイトル（学籍番号と氏名もつけること）
２） 入力した実対称行列ａ（行列の名前と、行番号をつける）
３） 計算過程(繰り返し数､ 非対角要素(k,l)の絶対値最大値,sin,cos,tan,
直交変換後の要素,繰り返し途中の行列a,p）
４） 収束すれば、解（固有値･固有ベクトル）、検算1)ｐａｐ, 2)ｐi
５） 収束しないときは「収束しない」
*/
#define NM 20 /* 最大の行・列 */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
// /* 関数のプロトタイプ宣言 */
void symInput(int, double [NM][NM]); /* 実対称行列の入力関数 */
void sOutput(char [], int, double [NM][NM]);/* 正方行列の出力関数 */
void sOutput2(char [], int, double [][NM], double [][NM]);/* 2 行列/行 */
void jacobif(int, double [NM][NM], double [NM][NM], int*);
void sMatrixCopy(int, double [][NM], double [][NM]);
void sTranspose(char [], int, double [][NM], double [][NM]);
void sProduct(char [], int, double [][NM], double [][NM], double[][NM]);
void inverse(int, double [][NM], double [][NM], double[][NM], double*, int*);
void main(void)
{
double a[NM][NM], p[NM][NM]; /* ａ:行列（戻りは対角に固有値）p:ﾓｰﾄﾞ行列 */
double ac[NM][NM], pta[NM][NM], ptap[NM][NM]; /* ａc:ac=a, pa=pt*a, pap=pa*p */
double pt[NM][NM], pi[NM][NM], work[NM][NM], det; /* pt=p' */
int n, i, j, cond; /* cond： ０なら、収束した */
char title[50];
char *pa; /* *pa:日付用文字列(ポインタ変数）*/
long nowtime;
time(&nowtime);
pa = ctime(&nowtime);
*(pa+24)='\0';
printf("＜%s＞\n", pa); /* 実行年月日などをコマンドプロンプトに表示する。 */
scanf("%s", title); /* データタイトルの入力 */
printf("<<< ヤコビ法（実対称行列の固有値解法）戸川の式 %s by k134084 関谷順太 >>>\n",title); /* プログラム、データタイトル表示 */
scanf("%d", &n); /* 次数ｎの入力 */
if(n>19) {
printf("プログラムの制限を超えています。");
exit(1);
}
symInput(n, a); /* 実対称行列Ａの入力 */
sMatrixCopy(n, a, ac);
jacobif(n, a, p, &cond); /* ヤコビ法関数を呼ぶ */
if(cond == 0){ /* 結果の出力 */
printf("解が求まりました\n 固有値 固有ベクトル\n");
for(j=1; j<=n; j++){
printf("%8.3f ", a[j][j]);
/* 固有ベクトルの表示 */
for(i=1; i<=n; i++)
printf(" %8.3f",p[i][j]);
printf("\n");
}
sTranspose("転置行列ｐt", n, p, pt);
sProduct("積ｐ'＊ａ", n, pt, ac, pta);
sProduct("積ｐ'ａｐ", n, pta, p, ptap);
inverse(n, p, pi, work, &det, &cond);
if (cond == 0)
sOutput(" 逆行列ｐi ", n, pi);
}else
printf("解が求まりません（収束しない）\n");
}
void jacobif(int n, double a[][NM],double p[][NM], int* cond)
{ /* ヤコビ法（実対称行列の固有値･固有ベクトル計算）戸川のp.117-124) */
/* n: 実対称行列の次数 ａ:行列（結果として対角に固有値）p:モード行列
cond： ０なら、収束した */
int i, j, k, km, l, m, nm1, ip1;
double eps, amax, akk, all,akl, wk, sd, t,c,s,pik, pil,aik, ail, s2,ss,cc;
double akj, alj;
eps = 1.0e-6;
for(i=1;i<=n; i++) /* 行列ｐの初期値をセット */
for(j=1;j<=n; j++)
if(i!=j) p[i][j]=0.0; else p[i][j]=1.0;
sOutput2("ヤコビ法、初期値ａ,ｐ", n, a, p);
nm1=n-1; m=100; /* m=(n*n-n)/2; m:繰り返し数 */
for(km=1; km<=m; km++) { /* 絶対値最大の非対角要素を探す */
amax=0.0;
for(i=1;i<=nm1; i++){
ip1=i+1;
for(j=ip1;j<=n; j++){
wk=fabs(a[i][j]);
if(wk<=amax) continue;
else { /* 最大値amax、行ｋ、列ｌのセット */
amax=wk; k=i; l=j;
}
}
}
if(amax<eps) { /* 収束した */
printf("(yacobi法で収束までの直交変換の繰り返し数: %d)\n", km-1);
*cond=0; return;
}
printf("(繰り返し数: %d､ 非対角要素の絶対値最大値:%f )", km-1,amax);
printf(" その要素のｋ：%d ｌ：%d\n", k, l);
akk=a[k][k]; all=a[l][l]; akl=a[k][l]; wk=akk-all;
sd=sqrt(wk*wk+4.0*akl*akl);
t=2.0*akl/(wk-sd); /* (1.9.59) */
c=1.0/sqrt(1.0+t*t); /* (1.9.60) */
s=c*t;
/* 直交変換 (戸川、p.119) */
for (j = 1; j <= n; j++) { /* 左からのＲtＡ乗算(6．19) */
akj = a[k][j]; alj =a[l][j];
???
}
for (i = 1; i <= n; i++) { /* 右からの（ＲtＡ）Ｒ乗算(6.20)*/
aik = a[i][k]; ail =a[i][l];
???
}
for(i=1; i<=n; i++){ /* 固有ベクトル乗算（右から）*/
pik=p[i][k]; pil=p[i][l];
???
}
sOutput2(" 繰り返し途中のａ,ｐ", n, a, p);
}
*cond=1;
}
void symInput(int n, double x[][NM])
{ /* 実対称行列の入力 関数
n:次数、x:実対称行列 */
int i,j;
for(i=1;i<=n; i++)
for(j=1;j<=n; j++){
scanf("%lf",&x[i][j]); /* 行iの成分ijの入力 */
if(i>j) /* 実対称か否かのチェック */
if(x[i][j] != x[j][i]){
printf("行列が対称でない x[%d][%d]\n", i,j);
exit(1);
}
}
sOutput("実対称行列Ａ", n, x); /* チェックプリント */
}
void sOutput(char matrixName[], int n, double x[][NM])
{ /* 正方行列の出力 関数 */
int i,j;
printf("\nMatrix: %s \n", matrixName); /*行列名の表示*/
printf("行/列 ");
for (j = 1; j <= n; j++) /* 列番号のループ */
printf(" %3d ", j); /* 見出し 列jの表示 */
printf("\n");
for(i=1;i<=n; i++){
printf("%3d ", i); /* 行番号の表示 */
for(j=1;j<=n; j++)
printf(" %8.3f", x[i][j]); /* 行iの成分ijの表示 */
printf("\n");
}
}
void sMatrixCopy(int n, double a[][NM], double ac[][NM])
{ /* 行列のコピー関数 */
int i,j; /* n:元数 a[n][n], ac[n][n]:行列 */
for(i=1;i<=n; i++) /* 行のループ */
for(j=1;j<=n; j++) /* 列のループ */
ac[i][j] = a[i][j]; /* 成分ijの代入 */
}
/* 正方行列(n,n)の積の計算 関数, 結果の行列の表示も行う。cc = aa * bb */
void sProduct(char title[], int n, double aa[][NM],
double bb[][NM], double cc[][NM])
{
int i,j,k; double x;
for(i=1;i<=n;i++)
for(k=1;k<=n;k++){
cc[i][k]=x=0.0;
for(j=1; j<=n; j++)
x += aa[i][j] * bb[j][k]; /* 成分の積和 */
cc[i][k]=x;
}
sOutput(title, n, cc); /* 積の出力 */
}
/* 正方行列の転置関数, 結果の転置行列の表示も行う。
入力行列 x[n][n]、ｘの転置行列xt[n][n] */
void sTranspose(char matrixName[], int n, double x[][NM], double xt[][NM])
{
int i,j;
for (i = 1; i<=n; i++)
for (j = 1; j <= n; j++) {
xt[j][i] = x[i][j]; /* 転置行列の成分j,i */
}
sOutput(matrixName, n, xt);
}
/* 行列の出力 関数2 正方行列の名前 matrixName[], x[n][n],x2[n][n]*/
void sOutput2(char matrixName[], int n, double x[][NM], double x2[][NM])
{
int i,j;
printf("Matrix: %s \n", matrixName); /*行列名の表示*/
printf("行/列 ");
for (j = 1;j <= n; j++) /* 列番号のループ for x */
printf(" %3d ", j); /* 列jの表示 */
printf(" ");
for (j = 1; j <= n; j++) /* 列番号のループ for x2 */
printf(" %3d ", j); /* 列jの表示 */
printf("\n");
for (i = 1; i <= n; i++){
printf("%3d ", i); /* 行番号の表示 */
for (j = 1; j <= n; j++) /* 列番号のループ for x */
printf(" %8.5lf", x[i][j]); /* xの成分ijの表示 */
printf(" ");
for (j = 1; j <= n; j++) /* 列番号のループ for x2 */
printf(" %8.5lf", x2[i][j]); /* x2の成分ijの表示 */
printf("\n");
}
}
/* 逆行列, 行列式の値計算 関数(消去法による） */
void inverse(int n, double a[][NM], double b[][NM], double cc[][NM],
double* det, int* cond)
{
/* n: 正方行列の次数､
ａ(input matrix)，ｂ(inversed matrix, a-1)，cc(work)：正方行列
det: 行列式|a|の値、 cond： ０なら逆行列がある。 */
int i,j,k,imax, kp1, m;
double eps, amax, absaik, work, akk, aik;
eps = 1.0e-5;
*det = 1.0;
m=0; /* 行の交換回数 */
for(i=1;i<=n;i++)
for(j=1;j<=n;j++){
cc[i][j]=a[i][j]; /* 行列ａのｃｃへの保存 */
if(i==j) b[i][j]=1.0; else b[i][j]=0.0;/* ｂを単位行列にセット */
}
for(k=1; k<=n; k++) { /* 行ｋの１からｎまでの繰り返し */
imax=k; /* 枢軸（ pivot )imaxの初期値k */
amax= fabs(a[k][k]); /* amaxに a[k][k]の絶対値をセット */
if(k != n){
kp1=k+1;
for(i=kp1; i<=n; i++){ /* 最大値の探索ループ */
absaik=fabs(a[i][k]); /* 列の要素をabsaikとする */
if(absaik <= amax) continue;/* absaik <= amaxなら、ループ端へ */
imax = i; amax= absaik; /* else、最大行imax=i,amax= absaik*/
}
}
if(amax <= eps) {
*cond=k; return; /* 逆行列が存在しない */
}
if(imax != k){ /* imaxとkが違うなら、行の入れ替え */
m++;
for(j=1; j<=n; j++){
work = a[imax][j]; a[imax][j] = a[k][j]; a[k][j] = work;
work = b[imax][j]; b[imax][j] = b[k][j]; b[k][j] = work;
}
}
akk = a[k][k]; /* 主対角要素をakkとする */
*det *= akk; /* 行列式の計算 */
for(j=1; j<=n; j++){ /* 消去ステップ１ ｋ行 */
a[k][j] = a[k][j]/akk; b[k][j] = b[k][j]/akk;
}
for(i=1; i<=n; i++){ /* 消去ステップ２ ｋでない行*/
if(i==k) continue;
aik = a[i][k];
for(j=1; j<=n; j++){
a[i][j] -= aik*a[k][j]; b[i][j] -= aik*b[k][j];
}
}
// sOutput2("枢軸(pivot)行入替えによる消去途中の行列ａ,ｂ ", n, n, a, b);
}
*det *= pow(-1.0,(double)m); /* (1.7.8)行列式の符号 */
*cond = 0; /* 逆行列があった */
for(i=1;i<=n;i++) /* 行列aをccを使い、元に戻す */
for(j=1;j<=n;j++)
a[i][j]=cc[i][j];
}
/*
(1) 収束が早い例
＜Wed Mar 13 15:28:16 2002＞
<<< ヤコビ法（実対称行列の固有値解法）戸川の式 戸川p.121-124 by k134084 関谷順太 >>>
Matrix: 実対称行列Ａ
行/列 1 2 3 4
1 5.000 4.000 1.000 1.000
2 4.000 5.000 1.000 1.000
3 1.000 1.000 4.000 2.000
4 1.000 1.000 2.000 4.000
Matrix: ヤコビ法、初期値ａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 5.00000 4.00000 1.00000 1.00000 1.00000 0.00000 0.00000 0.00000
2 4.00000 5.00000 1.00000 1.00000 0.00000 1.00000 0.00000 0.00000
3 1.00000 1.00000 4.00000 2.00000 0.00000 0.00000 1.00000 0.00000
4 1.00000 1.00000 2.00000 4.00000 0.00000 0.00000 0.00000 1.00000
(繰り返し数: 0､ 非対角要素の絶対値最大値:4.000000 ) その要素のｋ：1 ｌ：2
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 1.00000 0.00000 0.00000 0.00000 0.70711 0.70711 0.00000 0.00000
2 0.00000 9.00000 1.41421 1.41421 -0.70711 0.70711 0.00000 0.00000
3 0.00000 1.41421 4.00000 2.00000 0.00000 0.00000 1.00000 0.00000
4 0.00000 1.41421 2.00000 4.00000 0.00000 0.00000 0.00000 1.00000
(繰り返し数: 1､ 非対角要素の絶対値最大値:2.000000 ) その要素のｋ：3 ｌ：4
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 1.00000 0.00000 0.00000 0.00000 0.70711 0.70711 0.00000 0.00000
2 0.00000 9.00000 0.00000 2.00000 -0.70711 0.70711 0.00000 0.00000
3 0.00000 0.00000 2.00000 0.00000 0.00000 0.00000 0.70711 0.70711
4 0.00000 2.00000 0.00000 6.00000 0.00000 0.00000 -0.70711 0.70711
(繰り返し数: 2､ 非対角要素の絶対値最大値:2.000000 ) その要素のｋ：2 ｌ：4
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 1.00000 0.00000 0.00000 0.00000 0.70711 0.31623 0.00000 0.63246
2 0.00000 5.00000 0.00000 -0.00000 -0.70711 0.31623 0.00000 0.63246
3 0.00000 0.00000 2.00000 0.00000 0.00000 -0.63246 0.70711 0.31623
4 0.00000 -0.00000 0.00000 10.00000 0.00000 -0.63246 -0.70711 0.31623
(yacobi法で収束までの直交変換の繰り返し数: 3)
解が求まりました
固有値 固有ベクトル
1.000 0.707 -0.707 0.000 0.000
5.000 0.316 0.316 -0.632 -0.632
2.000 0.000 0.000 0.707 -0.707
10.000 0.632 0.632 0.316 0.316
Matrix: 転置行列ｐt
行/列 1 2 3 4
1 0.707 -0.707 0.000 0.000
2 0.316 0.316 -0.632 -0.632
3 0.000 0.000 0.707 -0.707
4 0.632 0.632 0.316 0.316
Matrix: 積ｐ'＊ａ
行/列 1 2 3 4
1 0.707 -0.707 0.000 0.000
2 1.581 1.581 -3.162 -3.162
3 0.000 0.000 1.414 -1.414
4 6.325 6.325 3.162 3.162
Matrix: 積ｐ'ａｐ
行/列 1 2 3 4
1 1.000 0.000 0.000 0.000
2 0.000 5.000 0.000 -0.000
3 0.000 0.000 2.000 0.000
4 0.000 0.000 0.000 10.000
Matrix: 逆行列ｐi
行/列 1 2 3 4
1 0.707 -0.707 -0.000 0.000
2 0.316 0.316 -0.632 -0.632
3 0.000 0.000 0.707 -0.707
4 0.632 0.632 0.316 0.316
(2) これは収束が遅い例である。
<<< ヤコビ法（実対称行列の固有値解法）戸川の式 多変量p.85 by k134084 関谷順太 >>>
Matrix: 実対称行列Ａ
行/列 1 2 3 4
1 1.000 0.863 0.732 0.920
2 0.863 1.000 0.896 0.883
3 0.732 0.896 1.000 0.783
4 0.920 0.883 0.783 1.000
Matrix: ヤコビ法、初期値ａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 1.00000 0.86320 0.73210 0.92050 1.00000 0.00000 0.00000 0.00000
2 0.86320 1.00000 0.89650 0.88270 0.00000 1.00000 0.00000 0.00000
3 0.73210 0.89650 1.00000 0.78290 0.00000 0.00000 1.00000 0.00000
4 0.92050 0.88270 0.78290 1.00000 0.00000 0.00000 0.00000 1.00000
(繰り返し数: 0､ 非対角要素の絶対値最大値:0.920500 ) その要素のｋ：1 ｌ：4
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 0.07950 -0.01379 -0.03592 0.00000 0.70711 0.00000 0.00000 0.70711
2 -0.01379 1.00000 0.89650 1.23454 0.00000 1.00000 0.00000 0.00000
3 -0.03592 0.89650 1.00000 1.07127 0.00000 0.00000 1.00000 0.00000
4 0.00000 1.23454 1.07127 1.92050 -0.70711 0.00000 0.00000 0.70711
(繰り返し数: 1､ 非対角要素の絶対値最大値:1.234538 ) その要素のｋ：2 ｌ：4
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 0.07950 -0.01133 -0.03592 -0.00786 0.70711 -0.40332 0.00000 0.58080
2 -0.01133 0.14271 0.12533 0.00000 0.00000 0.82138 0.00000 0.57038
3 -0.03592 0.12533 1.00000 1.39126 0.00000 0.00000 1.00000 0.00000
4 -0.00786 0.00000 1.39126 2.77779 -0.70711 -0.40332 0.00000 0.58080
(繰り返し数: 2､ 非対角要素の絶対値最大値:1.391264 ) その要素のｋ：3 ｌ：4
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 0.07950 -0.01133 -0.02773 -0.02415 0.70711 -0.40332 -0.27903 0.50939
2 -0.01133 0.14271 0.10992 0.06021 0.00000 0.82138 -0.27402 0.50025
3 -0.02773 0.10992 0.23791 -0.00000 0.00000 0.00000 0.87704 0.48041
4 -0.02415 0.06021 -0.00000 3.53988 -0.70711 -0.40332 -0.27903 0.50939
(繰り返し数: 3､ 非対角要素の絶対値最大値:0.109922 ) その要素のｋ：2 ｌ：3
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 0.07950 0.00575 -0.02939 -0.02415 0.70711 -0.18397 -0.45462 0.50939
2 0.00575 0.07052 0.00000 0.05033 0.00000 0.83699 0.22182 0.50025
3 -0.02939 -0.00000 0.31010 0.03305 0.00000 -0.48142 0.73310 0.48041
4 -0.02415 0.05033 0.03305 3.53988 -0.70711 -0.18397 -0.45462 0.50939
(繰り返し数: 4､ 非対角要素の絶対値最大値:0.050330 ) その要素のｋ：2 ｌ：4
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 0.07950 0.00610 -0.02939 -0.02407 0.70711 -0.19133 -0.45462 0.50667
2 0.00610 0.06979 -0.00048 -0.00000 0.00000 0.82964 0.22182 0.51234
3 -0.02939 -0.00048 0.31010 0.03305 0.00000 -0.48834 0.73310 0.47338
4 -0.02407 0.00000 0.03305 3.54061 -0.70711 -0.19133 -0.45462 0.50667
(繰り返し数: 5､ 非対角要素の絶対値最大値:0.033048 ) その要素のｋ：3 ｌ：4
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 0.07950 0.00610 -0.02914 -0.02437 0.70711 -0.19133 -0.45978 0.50199
2 0.00610 0.06979 -0.00048 -0.00000 0.00000 0.82964 0.21657 0.51458
3 -0.02914 -0.00048 0.30976 0.00000 0.00000 -0.48834 0.72822 0.48086
4 -0.02437 -0.00000 0.00000 3.54095 -0.70711 -0.19133 -0.45978 0.50199
(繰り返し数: 6､ 非対角要素の絶対値最大値:0.029145 ) その要素のｋ：1 ｌ：3
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 0.07587 0.00600 0.00000 -0.02418 0.64483 -0.19133 -0.54369 0.50199
2 0.00600 0.06979 -0.00123 -0.00000 0.02678 0.82964 0.21491 0.51458
3 0.00000 -0.00123 0.31339 0.00301 0.09005 -0.48834 0.72263 0.48086
4 -0.02418 -0.00000 0.00301 3.54095 -0.75853 -0.19133 -0.36882 0.50199
(繰り返し数: 7､ 非対角要素の絶対値最大値:0.024181 ) その要素のｋ：1 ｌ：4
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 0.07570 0.00600 0.00002 -0.00000 0.64831 -0.19133 -0.54369 0.49748
2 0.00600 0.06979 -0.00123 -0.00005 0.03037 0.82964 0.21491 0.51438
3 0.00002 -0.00123 0.31339 0.00301 0.09340 -0.48834 0.72263 0.48022
4 0.00000 -0.00005 0.00301 3.54112 -0.75501 -0.19133 -0.36882 0.50727
(繰り返し数: 8､ 非対角要素の絶対値最大値:0.005996 ) その要素のｋ：1 ｌ：2
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 0.06606 -0.00000 0.00106 0.00004 0.50495 0.44938 -0.54369 0.49748
2 0.00000 0.07943 -0.00063 -0.00002 -0.68837 0.46408 0.21491 0.51438
3 0.00106 -0.00063 0.31339 0.00301 0.46397 -0.17869 0.72263 0.48022
4 0.00004 -0.00002 0.00301 3.54112 -0.23641 -0.74213 -0.36882 0.50727
(繰り返し数: 9､ 非対角要素の絶対値最大値:0.003013 ) その要素のｋ：3 ｌ：4
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 0.06606 -0.00000 0.00106 0.00004 0.50495 0.44938 -0.54415 0.49697
2 0.00000 0.07943 -0.00063 -0.00003 -0.68837 0.46408 0.21443 0.51458
3 0.00106 -0.00063 0.31339 -0.00000 0.46397 -0.17869 0.72218 0.48089
4 0.00004 -0.00003 0.00000 3.54112 -0.23641 -0.74213 -0.36929 0.50693
(繰り返し数: 10､ 非対角要素の絶対値最大値:0.001056 ) その要素のｋ：1 ｌ：3
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 0.06606 0.00000 -0.00000 0.00004 0.50727 0.44938 -0.54199 0.49697
2 0.00000 0.07943 -0.00063 -0.00003 -0.68928 0.46408 0.21149 0.51458
3 0.00000 -0.00063 0.31339 0.00000 0.46089 -0.17869 0.72415 0.48089
4 0.00004 -0.00003 0.00000 3.54112 -0.23484 -0.74213 -0.37030 0.50693
(繰り返し数: 11､ 非対角要素の絶対値最大値:0.000632 ) その要素のｋ：2 ｌ：3
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 0.06606 0.00000 -0.00000 0.00004 0.50727 0.44791 -0.54320 0.49697
2 0.00000 0.07943 0.00000 -0.00003 -0.68928 0.46465 0.21023 0.51458
3 -0.00000 -0.00000 0.31339 0.00000 0.46089 -0.17673 0.72463 0.48089
4 0.00004 -0.00003 0.00000 3.54112 -0.23484 -0.74313 -0.36829 0.50693
(繰り返し数: 12､ 非対角要素の絶対値最大値:0.000041 ) その要素のｋ：1 ｌ：4
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 0.06606 0.00000 -0.00000 -0.00000 0.50727 0.44791 -0.54320 0.49698
2 0.00000 0.07943 0.00000 -0.00003 -0.68929 0.46465 0.21023 0.51457
3 -0.00000 -0.00000 0.31339 0.00000 0.46088 -0.17673 0.72463 0.48090
4 -0.00000 -0.00003 0.00000 3.54112 -0.23484 -0.74313 -0.36829 0.50692
(繰り返し数: 13､ 非対角要素の絶対値最大値:0.000025 ) その要素のｋ：2 ｌ：4
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 0.06606 0.00000 -0.00000 -0.00000 0.50727 0.44792 -0.54320 0.49697
2 0.00000 0.07943 0.00000 -0.00000 -0.68929 0.46465 0.21023 0.51457
3 -0.00000 0.00000 0.31339 0.00000 0.46088 -0.17673 0.72463 0.48090
4 -0.00000 0.00000 0.00000 3.54112 -0.23484 -0.74313 -0.36829 0.50693
(繰り返し数: 14､ 非対角要素の絶対値最大値:0.000003 ) その要素のｋ：1 ｌ：2
Matrix: 繰り返し途中のａ,ｐ
行/列 1 2 3 4 1 2 3 4
1 0.06606 -0.00000 -0.00000 -0.00000 0.50718 0.44802 -0.54320 0.49697
2 0.00000 0.07943 0.00000 -0.00000 -0.68938 0.46452 0.21023 0.51457
3 -0.00000 0.00000 0.31339 0.00000 0.46092 -0.17663 0.72463 0.48090
4 -0.00000 -0.00000 0.00000 3.54112 -0.23469 -0.74317 -0.36829 0.50693
(yacobi法で収束までの直交変換の繰り返し数: 15)
解が求まりました
固有値 固有ベクトル
0.066 0.507 -0.689 0.461 -0.235
0.079 0.448 0.465 -0.177 -0.743
0.313 -0.543 0.210 0.725 -0.368
3.541 0.497 0.515 0.481 0.507
Matrix: 転置行列ｐt
行/列 1 2 3 4
1 0.507 -0.689 0.461 -0.235
2 0.448 0.465 -0.177 -0.743
3 -0.543 0.210 0.725 -0.368
4 0.497 0.515 0.481 0.507
Matrix: 積ｐ'＊ａ
行/列 1 2 3 4
1 0.034 -0.046 0.030 -0.016
2 0.036 0.037 -0.014 -0.059
3 -0.170 0.066 0.227 -0.115
4 1.760 1.822 1.703 1.795
Matrix: 積ｐ'ａｐ
行/列 1 2 3 4
1 0.066 0.000 -0.000 -0.000
2 0.000 0.079 0.000 -0.000
3 -0.000 0.000 0.313 0.000
4 -0.000 -0.000 0.000 3.541
Matrix: 逆行列ｐi
行/列 1 2 3 4
1 0.507 -0.689 0.461 -0.235
2 0.448 0.465 -0.177 -0.743
3 -0.543 0.210 0.725 -0.368
4 0.497 0.515 0.481 0.507
*/