「プログラムソース」
４−１、-------------------------------------------------------------
#include<stdio.h>
#include<math.h>

#define N 2	/* ここの数字をかえる事で行列数を変更できる、今は2×2 */

int main(){
	int t, i, k;
	int gyou, retu;
	double r, r1, r2, kaihi_r;
	double A[N][N];
	double x[N], kaihi_x[N];

	for(gyou=0;gyou<N;gyou++){
		for(retu=0;retu<N;retu++){
			printf("%d行%d列を入力して下さい。\n",gyou+1, retu+1);
			scanf("%lf",&A[gyou][retu]);
		}
	}
	
	printf("入力した行列は以下です。\n");
	for(gyou=0;gyou<N;gyou++){
		for(retu=0;retu<N;retu++){
			printf("%f　",A[gyou][retu]);
		}
		printf("\n");
	}
	printf("\n");
	printf("累乗法により固有値を求めます。\n");
	
	t = 0;
	for(k=0; k<N; k++)
		x[k] = 1.0;

	while(1){
		do{
			t = t + 1;

			for(i=0; i<N; i++){
				kaihi_x[i] = 0;
				for( k=0; k<N; k++ ){
					kaihi_x[i] += A[i][k] * x[k];
				}
			}

			r1 = r2 = 0.0;
	
			for( k=0; k<N; k++ ){
				r1 += kaihi_x[k] * x[k];
			}
			for( k=0; k<N; k++ ){
				r2 += x[k] * x[k];
			}
			for( k=0; k<N; k++ ){
				x[k] = kaihi_x[k];
			}
		
			r = r1 / r2;
			printf( "r=%g\n", r );

		}while(t == 1);

		if(fabs( r - kaihi_r ) < 10e-6){
			printf( "固有値＝%g\n", r );
			break;
		}
		kaihi_r = r;
	}
	return 0;	
}

--------------------------------------------------------------------------------
４−２、

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define N 2		/* ここの数字をかえる事で行列数を変更できる、今は2×2 */

double kannsuu1(double []);
double kannsuu2(double [],double []);

int main()
{
	
	double a[N][N],b[N][N],q[N][N],r[N][N];
	double koyuu[N][N],koyuu_kaihi[N][N];
	double cr;

	int gyou,retu;
	int h=0;
	int k,m,p;
    
	
	printf("ＱＲ法により、固有値を求めます。\n");
    	printf("2×2の対称行列を入力して下さい。\n");
	for(gyou=0;gyou<N;gyou++){
		for(retu=0;retu<N;retu++){
			printf("\n%d行%d列を入力して下さい。\n",gyou+1, retu+1);
			scanf("%lf",&a[gyou][retu]);
		}
	}
	
	printf("入力した行列は以下です。\n");
	for(gyou=0;gyou<N;gyou++){
		for(retu=0;retu<N;retu++){
			printf("%f　",a[gyou][retu]);
		}
		printf("\n");
	}
	printf("\n");
	 

	for(gyou = 0; gyou < N;gyou++) {
		for(retu = 0; retu < N;retu++) {
			koyuu[gyou][retu] = 0;
		}
	}
     
	for(gyou = 0; gyou < N;gyou++){
		koyuu[gyou][gyou] = 1;
	}

	while(h < 100){
		for(gyou = 0; gyou < N;gyou++) {
			for(retu = 0; retu < N;retu++) {
				b[gyou][retu] = 0;
				r[gyou][retu] = 0;
				q[gyou][retu] = 0;
			}
		}

		for(gyou = 0;gyou < N;gyou++) {
			b[0][gyou] = a[0][gyou];
		}
		
		cr = r[0][0] = kannsuu1(b[0]);
		
		for(gyou = 0;gyou < N;gyou++) {
			q[0][gyou] = b[0][gyou]/cr;
		}

		for(k = 1;k < N;k++) {
			for(retu=0;retu < k;retu++) {
				r[k][retu] = kannsuu2(a[k],q[retu]);
			}
			
			for(m = 0;m < N;m++) {
				b[k][m] = a[k][m];
			}
			
			for(m = 0;m < N;m++) {
				for(p = 0;p < k;p++) {
					b[k][m] -= r[k][p]*q[p][m];
				}
			}

			cr = r[k][k] = kannsuu1(b[k]);
			
			for(gyou = 0;gyou < N;gyou++) {
				q[k][gyou] = b[k][gyou]/cr;
			}
		}

		for(gyou = 0;gyou < N;gyou++) {
			for(retu = 0;retu < N;retu++) {
				a[gyou][retu] = 0.0;
				for(k = 0;k < N;k++) {
					a[gyou][retu] += r[k][gyou]*q[retu][k];
				}
			}
		}

		for(gyou = 0;gyou < N;gyou++) {
			for(retu = 0;retu < N;retu++) {
				koyuu_kaihi[gyou][retu] = 0;
				for(k = 0;k < N;k++) {
					koyuu_kaihi[gyou][retu] += koyuu[gyou][k]*q[retu][k];
				}
			}
		}

		for(gyou = 0;gyou < N;gyou++) {
			for(retu = 0;retu < N;retu++) {
				koyuu[gyou][retu] = koyuu_kaihi[gyou][retu];
			}
		}

		printf("\nベクトルＡ\n");

		for(gyou = 0; gyou < N;gyou++) {
			for(retu = 0; retu < N;retu++) {
				printf("%f　",a[retu][gyou]);
			}
			printf("\n");
		}
		printf("\n\n");
		printf("\n固有ベクトル\n");

		for(gyou = 0; gyou < N;gyou++) {
			for(retu = 0; retu < N;retu++) {
				printf("%f　",koyuu[gyou][retu]);
			}
			printf("\n");
		}         
		printf("\n\n");
		h++;
	}
	printf("ＱＲ法最終結果（ループ回数１００回）\n");
	printf("\nベクトルＡ\n\n");

		for(gyou = 0; gyou < N;gyou++) {
			for(retu = 0; retu < N;retu++) {
				printf("%f　",a[retu][gyou]);
			}
			printf("\n");
		}
		printf("\n\n");
		printf("\n固有ベクトル\n\n");

		for(gyou = 0; gyou < N;gyou++) {
			for(retu = 0; retu < N;retu++) {
				printf("%f　",koyuu[gyou][retu]);
			}
			printf("\n");
		}
	return 0;
}


double kannsuu1(double a[]){
	int i;
	double s = 0.0;

	for(i = 0;i < N;i++) {
		s += pow(a[i],2);
	}
	
	return sqrt(s);
}


double kannsuu2(double a[],double b[]){
	int i;
	double s = 0.0;

	for(i = 0;i < N;i++) {
		s += a[i]*b[i];
	}
	return s;
}


-------------------------------------------------------------------------

「実行結果」
２−１--------------------------------------------------------------------------
C:\Documents and Settings\i865gpen4\デスクトップ\色んな資料\数値計算プログラミン
グ\数値計算課題４>ruijyouhou_t
1行1列を入力して下さい。
1
1行2列を入力して下さい。
2
2行1列を入力して下さい。
3
2行2列を入力して下さい。
1
入力した行列は以下です。
1.000000　2.000000　
3.000000　1.000000　

累乗法により固有値を求めます。
r=3.5
r=3.4
r=3.46552
r=3.44189
r=3.45253
r=3.44819
r=3.45003
r=3.44926
r=3.44959
r=3.44945
r=3.44951
r=3.44948
r=3.44949
r=3.44949
固有値＝3.44949
---------------------------------------------------------------------------
２−２（ループは実際のソースは１００回していますが、それだとメールに実行結果を記述できないので、対角化が完了する２０回程度で実行結果を出力しておきました。）

C:\Documents and Settings\i865gpen4\My Documents\プログラムソース>qrhou_t
ＱＲ法により、固有値を求めます。
2×2の対称行列を入力して下さい。

1行1列を入力して下さい。
1

1行2列を入力して下さい。
2

2行1列を入力して下さい。
2

2行2列を入力して下さい。
1
入力した行列は以下です。
1.000000　2.000000　
2.000000　1.000000　


ベクトルＡ
2.600000　1.200000　
1.200000　-0.600000　



固有ベクトル
0.447214　0.894427　
0.894427　-0.447214　



ベクトルＡ
2.951220　0.439024　
0.439024　-0.951220　



固有ベクトル
0.780869　-0.624695　
0.624695　0.780869　



ベクトルＡ
2.994521　0.147945　
0.147945　-0.994521　



固有ベクトル
0.680451　0.732793　
0.732793　-0.680451　



ベクトルＡ
2.999390　0.049375　
0.049375　-0.999390　



固有ベクトル
0.715782　-0.698324　
0.698324　0.715782　



ベクトルＡ
2.999932　0.016461　
0.016461　-0.999932　



固有ベクトル
0.704191　0.710011　
0.710011　-0.704191　



ベクトルＡ
2.999992　0.005487　
0.005487　-0.999992　



固有ベクトル
0.708076　-0.706136　
0.706136　0.708076　



ベクトルＡ
2.999999　0.001829　
0.001829　-0.999999　



固有ベクトル
0.706783　0.707430　
0.707430　-0.706783　



ベクトルＡ
3.000000　0.000610　
0.000610　-1.000000　



固有ベクトル
0.707215　-0.706999　
0.706999　0.707215　



ベクトルＡ
3.000000　0.000203　
0.000203　-1.000000　



固有ベクトル
0.707071　0.707143　
0.707143　-0.707071　



ベクトルＡ
3.000000　0.000068　
0.000068　-1.000000　



固有ベクトル
0.707119　-0.707095　
0.707095　0.707119　



ベクトルＡ
3.000000　0.000023　
0.000023　-1.000000　



固有ベクトル
0.707103　0.707111　
0.707111　-0.707103　



ベクトルＡ
3.000000　0.000008　
0.000008　-1.000000　



固有ベクトル
0.707108　-0.707105　
0.707105　0.707108　



ベクトルＡ
3.000000　0.000003　
0.000003　-1.000000　



固有ベクトル
0.707106　0.707107　
0.707107　-0.707106　



ベクトルＡ
3.000000　0.000001　
0.000001　-1.000000　



固有ベクトル
0.707107　-0.707107　
0.707107　0.707107　



ベクトルＡ
3.000000　0.000000　
0.000000　-1.000000　



固有ベクトル
0.707107　0.707107　
0.707107　-0.707107　



ベクトルＡ
3.000000　0.000000　
0.000000　-1.000000　



固有ベクトル
0.707107　-0.707107　
0.707107　0.707107　



ベクトルＡ
3.000000　0.000000　
0.000000　-1.000000　



固有ベクトル
0.707107　0.707107　
0.707107　-0.707107　



ベクトルＡ
3.000000　0.000000　
0.000000　-1.000000　



固有ベクトル
0.707107　-0.707107　
0.707107　0.707107　



ベクトルＡ
3.000000　0.000000　
0.000000　-1.000000　



固有ベクトル
0.707107　0.707107　
0.707107　-0.707107　



ベクトルＡ
3.000000　0.000000　
0.000000　-1.000000　



固有ベクトル
0.707107　-0.707107　
0.707107　0.707107　


ＱＲ法最終結果（ループ回数２０回）

ベクトルＡ

3.000000　0.000000　
0.000000　-1.000000　



固有ベクトル

0.707107　-0.707107　
0.707107　0.707107　

--------------------------------------------------------------------------------

「考察」
問題２の検討
累乗法は絶対最大固有値と、その次に絶対値の大きな固有値との関係で、大きく収束状況が影響される。

また、ＱＲ法に関しては今回きっちりと固有値を出力するところまでは至らなかったが、実際にもっと完成度の高いプログラムだと、ループ回数は良好な結果が得られるが、累乗のオーダーの際に実行時間が延びてしまうようである。
実行時間は累乗法と比べて大幅に時間がかかるようだ。
ＱＲ法は行列の大きさよりも、固有値の大きさに影響を受けるようだ。

ＱＲ法のプログラムはとても煩雑で、とてもむずかしかった。
ウェブ上の様々なプログラムソースを参考にして、自分なりに色々工夫してやってみた。
対称行列Ａの対角化が完了した際の、終了条件をうまく表すことができず、しょうがないのでループ回数を強引に指定し、それを終了条件とした。固有値を求めることが目的であったが、固有ベクトルを求める所までしかわからなかった。あとちょっと工夫すれば、固有値を求めることができそうなんだが。


累乗法とＱＲ法の理論的なことも、恐らくほとんど理解できていないので、これを提出した後にもう一度考え直そうと思う。
後、行列数を任意に指定できるようにしようと思ったが、ソースの一番上のdefineのNを変えることで、それと同等の意味を持ちえるプログラムを作ればいいと思ったので、あえて自分で入力するようにするのをやめた。