情報学部 菅沼ホーム 目次 索引

最大固有値と固有ベクトル（べき乗法）

    1. A. C++
    2. B. Java
    3. C. JavaScript
    4. D. PHP
    5. E. Ruby
    6. F. Python
    7. G. C#
    8. H. VB

　　プログラムは，行列の固有値と固有ベクトルを，固有値の絶対値が最大のものから順に求めていく方法（べき乗法）です．

  1. C++ 

    /****************************************/
/* 最大固有値と固有ベクトル（べき乗法） */
/*      coded by Y.Suganuma             */
/****************************************/
#include <stdio.h>

int power(int, int, int, int, double, double **, double **, double **,
          double *, double **, double *, double *, double *);

int main()
{
	double **A, **B, **C, *r, **v, *v0, *v1, *v2, eps;
	int i1, i2, ind, ct, m, n;
					// データの設定
	ct  = 200;
	eps = 1.0e-10;
	n   = 3;
	m   = 15;

	A   = new double * [n];
	B   = new double * [n];
	C   = new double * [n];
	v   = new double * [n];
	for (i1 = 0; i1 < n; i1++) {
		A[i1] = new double [n];
		B[i1] = new double [n];
		C[i1] = new double [n];
		v[i1] = new double [n];
	}

	r  = new double [n];
	v0 = new double [n];
	v1 = new double [n];
	v2 = new double [n];

	A[0][0] = 11.0;
	A[0][1] = 6.0;
	A[0][2] = -2.0;

	A[1][0] = -2.0;
	A[1][1] = 18.0;
	A[1][2] = 1.0;

	A[2][0] = -12.0;
	A[2][1] = 24.0;
	A[2][2] = 13.0;

	for (i1 = 0; i1 < n; i1++) {
		for (i2 = 0; i2 < n; i2++)
			B[i1][i2] = 0.0;
		B[i1][i1] = 1.0;
		v0[i1]    = 1.0;
	}
					// 計算
	ind = power(0, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
					// 出力
	if (ind == 0)
		printf("固有値が求まりませんでした！\n");
	else {
		for (i1 = 0; i1 < ind; i1++) {
			printf("固有値 %f", r[i1]);
			printf(" 固有ベクトル ");
			for (i2 = 0; i2 < n; i2++)
				printf(" %f", v[i1][i2]);
			printf("\n");
		}
	}

	for (i1 = 0; i1 < n; i1++) {
		delete [] A[i1];
		delete [] B[i1];
		delete [] C[i1];
		delete [] v[i1];
	}
	delete [] A;
	delete [] B;
	delete [] C;
	delete [] v;
	delete [] r;
	delete [] v0;
	delete [] v1;
	delete [] v2;

	return 0;
}

/****************************************/
/* 最大固有値と固有ベクトル（べき乗法） */
/*      i : 何番目の固有値かを示す      */
/*      n : 次数                        */
/*      m : 丸め誤差の修正回数          */
/*      ct : 最大繰り返し回数           */
/*      eps : 収束判定条件              */
/*      A : 対象とする行列              */
/*      B : nxnの行列，最初は，単位行列 */
/*      C : 作業域，nxnの行列           */
/*      r : 固有値                      */
/*      v : 各行が固有ベクトル（nxn）   */
/*      v0 : 固有ベクトルの初期値       */
/*      v1,v2 : 作業域（n次元ベクトル） */
/*      return : 求まった固有値の数     */
/*      coded by Y.Suganuma             */
/****************************************/
#include <math.h>

int power(int i, int n, int m, int ct, double eps, double **A, double **B,
          double **C, double *r, double **v, double *v0, double *v1, double *v2)
{
	double l1, l2, x;
	int i1, i2, i3, ind, k, k1;
					// 初期設定
	ind = i;
	k   = 0;
	l1  = 0.0;
	for (i1 = 0; i1 < n; i1++) {
		l1     += v0[i1] * v0[i1];
		v1[i1]  = v0[i1];
	}
	l1 = sqrt(l1);
					// 繰り返し計算
	while (k < ct) {
						// 丸め誤差の修正
		if (k%m == 0) {
			l2 = 0.0;
			for (i1 = 0; i1 < n; i1++) {
				v2[i1] = 0.0;
				for (i2 = 0; i2 < n; i2++)
					v2[i1] += B[i1][i2] * v1[i2];
				l2 += v2[i1] * v2[i1];
			}
			l2 = sqrt(l2);
			for (i1 = 0; i1 < n; i1++)
				v1[i1] = v2[i1] / l2;
		}
						// 次の近似
		l2 = 0.0;
		for (i1 = 0; i1 < n; i1++) {
			v2[i1] = 0.0;
			for (i2 = 0; i2 < n; i2++)
				v2[i1] += A[i1][i2] * v1[i2];
			l2 += v2[i1] * v2[i1];
		}
		l2 = sqrt(l2);
		for (i1 = 0; i1 < n; i1++)
			v2[i1] /= l2;
						// 収束判定
							// 収束した場合
		if (fabs((l2-l1)/l1) < eps) {
			k1 = -1;
			for (i1 = 0; i1 < n && k1 < 0; i1++) {
				if (fabs(v2[i1]) > 0.001) {
					k1 = i1;
					if (v2[k1]*v1[k1] < 0.0)
						l2 = -l2;
				}
			}
			k    = ct;
			r[i] = l2;
			for (i1 = 0; i1 < n; i1++)
				v[i][i1] = v2[i1];
			if (i == n-1)
				ind = i + 1;
			else {
				for (i1 = 0; i1 < n; i1++) {
					for (i2 = 0; i2 < n; i2++) {
						C[i1][i2] = 0.0;
						for (i3 = 0; i3 < n; i3++) {
							x          = (i1 == i3) ? A[i1][i3] - l2 : A[i1][i3];
							C[i1][i2] += x * B[i3][i2];
						}
					}
				}
				for (i1 = 0; i1 < n; i1++) {
					for (i2 = 0; i2 < n; i2++)
						B[i1][i2] = C[i1][i2];
				}
				for (i1 = 0; i1 < n; i1++) {
					v1[i1] = 0.0;
					for (i2 = 0; i2 < n; i2++)
						v1[i1] += B[i1][i2] * v0[i2];
				}
				for (i1 = 0; i1 < n; i1++)
					v0[i1] = v1[i1];
				ind = power(i+1, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
			}
		}
							// 収束しない場合
		else {
			for (i1 = 0; i1 < n; i1++)
				v1[i1] = v2[i1];
			l1 = l2;
			k++;
		}
	}

	return ind;
}

  2. Java 

    /****************************************/
/* 最大固有値と固有ベクトル（べき乗法） */
/*      coded by Y.Suganuma             */
/****************************************/
import java.io.*;

public class Test {
	public static void main(String args[]) throws IOException
	{
		double A[][], B[][], C[][], r[], v[][], v0[], v1[], v2[], eps;
		int i1, i2, ind, ct, m, n;
					// データの設定
		ct  = 200;
		eps = 1.0e-10;
		n   = 3;
		m   = 15;

		A   = new double [n][n];
		B   = new double [n][n];
		C   = new double [n][n];
		v   = new double [n][n];

		r  = new double [n];
		v0 = new double [n];
		v1 = new double [n];
		v2 = new double [n];

		A[0][0] = 11.0;
		A[0][1] = 6.0;
		A[0][2] = -2.0;

		A[1][0] = -2.0;
		A[1][1] = 18.0;
		A[1][2] = 1.0;

		A[2][0] = -12.0;
		A[2][1] = 24.0;
		A[2][2] = 13.0;

		for (i1 = 0; i1 < n; i1++) {
			for (i2 = 0; i2 < n; i2++)
				B[i1][i2] = 0.0;
			B[i1][i1] = 1.0;
			v0[i1]    = 1.0;
		}
					// 計算
		ind = power(0, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
					// 出力
		if (ind == 0)
			System.out.println("固有値が求まりませんでした！");
		else {
			for (i1 = 0; i1 < ind; i1++) {
				System.out.print("固有値 " + r[i1]);
				System.out.print(" 固有ベクトル ");
				for (i2 = 0; i2 < n; i2++)
					System.out.print(" " + v[i1][i2]);
				System.out.println();
			}
		}
	}

	/****************************************/
	/* 最大固有値と固有ベクトル（べき乗法） */
	/*      i : 何番目の固有値かを示す      */
	/*      n : 次数                        */
	/*      m : 丸め誤差の修正回数          */
	/*      ct : 最大繰り返し回数           */
	/*      eps : 収束判定条件              */
	/*      A : 対象とする行列              */
	/*      B : nxnの行列，最初は，単位行列 */
	/*      C : 作業域，nxnの行列           */
	/*      r : 固有値                      */
	/*      v : 各行が固有ベクトル（nxn）   */
	/*      v0 : 固有ベクトルの初期値       */
	/*      v1,v2 : 作業域（n次元ベクトル） */
	/*      return : 求まった固有値の数     */
	/*      coded by Y.Suganuma             */
	/****************************************/
	static int power(int i, int n, int m, int ct, double eps, double A[][], double B[][],
	                 double C[][], double r[], double v[][], double v0[], double v1[], double v2[])
	{
		double l1, l2, x;
		int i1, i2, i3, ind, k, k1;
					// 初期設定
		ind = i;
		k   = 0;
		l1  = 0.0;
		for (i1 = 0; i1 < n; i1++) {
			l1     += v0[i1] * v0[i1];
			v1[i1]  = v0[i1];
		}
		l1 = Math.sqrt(l1);
					// 繰り返し計算
		while (k < ct) {
						// 丸め誤差の修正
			if (k%m == 0) {
				l2 = 0.0;
				for (i1 = 0; i1 < n; i1++) {
					v2[i1] = 0.0;
					for (i2 = 0; i2 < n; i2++)
						v2[i1] += B[i1][i2] * v1[i2];
					l2 += v2[i1] * v2[i1];
				}
				l2 = Math.sqrt(l2);
				for (i1 = 0; i1 < n; i1++)
					v1[i1] = v2[i1] / l2;
			}
						// 次の近似
			l2 = 0.0;
			for (i1 = 0; i1 < n; i1++) {
				v2[i1] = 0.0;
				for (i2 = 0; i2 < n; i2++)
					v2[i1] += A[i1][i2] * v1[i2];
				l2 += v2[i1] * v2[i1];
			}
			l2 = Math.sqrt(l2);
			for (i1 = 0; i1 < n; i1++)
				v2[i1] /= l2;
						// 収束判定
							// 収束した場合
			if (Math.abs((l2-l1)/l1) < eps) {
				k1 = -1;
				for (i1 = 0; i1 < n && k1 < 0; i1++) {
					if (Math.abs(v2[i1]) > 0.001) {
						k1 = i1;
						if (v2[k1]*v1[k1] < 0.0)
							l2 = -l2;
					}
				}
				k    = ct;
				r[i] = l2;
				for (i1 = 0; i1 < n; i1++)
					v[i][i1] = v2[i1];
				if (i == n-1)
					ind = i + 1;
				else {
					for (i1 = 0; i1 < n; i1++) {
						for (i2 = 0; i2 < n; i2++) {
							C[i1][i2] = 0.0;
							for (i3 = 0; i3 < n; i3++) {
								x          = (i1 == i3) ? A[i1][i3] - l2 : A[i1][i3];
								C[i1][i2] += x * B[i3][i2];
							}
						}
					}
					for (i1 = 0; i1 < n; i1++) {
						for (i2 = 0; i2 < n; i2++)
							B[i1][i2] = C[i1][i2];
					}
					for (i1 = 0; i1 < n; i1++) {
						v1[i1] = 0.0;
						for (i2 = 0; i2 < n; i2++)
							v1[i1] += B[i1][i2] * v0[i2];
					}
					for (i1 = 0; i1 < n; i1++)
						v0[i1] = v1[i1];
					ind = power(i+1, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
				}
			}
							// 収束しない場合
			else {
				for (i1 = 0; i1 < n; i1++)
					v1[i1] = v2[i1];
				l1 = l2;
				k++;
			}
		}

		return ind;
	}
}

  3. JavaScript

    　　ここをクリックすると，任意のデータに対して，画面上で解を得ることができます． 
    <!DOCTYPE HTML>

<HTML>

<HEAD>

	<TITLE>最大固有値と固有ベクトル（べき乗法）</TITLE>
	<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=utf-8">
	<SCRIPT TYPE="text/javascript">
		function main()
		{
					// データの設定
			let eps = 1.0e-10;
			let ct  = parseInt(document.getElementById("trial").value);
			let n   = parseInt(document.getElementById("order").value);
			let m   = parseInt(document.getElementById("mod").value);
			let A  = new Array();
			for (let i1 = 0; i1 < n; i1++)
				A[i1] = new Array();
			let s = (document.getElementById("ar").value).split(/ {1,}|\n{1,}/);
			let k = 0;
			for (let i1 = 0; i1 < n; i1++) {
				for (let i2 = 0; i2 < n; i2++) {
					A[i1][i2] = parseFloat(s[k]);
					k++;
				}
			}
			let B = new Array();
			for (let i1 = 0; i1 < n; i1++) {
				B[i1] = new Array();
				for (let i2 = 0; i2 < n; i2++)
					B[i1][i2] = 0.0;
				B[i1][i1] = 1.0;
			}
			let C = new Array();
			for (let i1 = 0; i1 < n; i1++)
				C[i1] = new Array();
			let v = new Array();
			for (let i1 = 0; i1 < n; i1++)
				v[i1] = new Array();
			let r = new Array();
			let v0 = (document.getElementById("init").value).split(/ {1,}/);
			let v1 = new Array();
			let v2 = new Array();

			ind = power(0, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
					// 出力
			if (ind == 0)
				document.getElementById("ans").value = "解を求めることができません！";
			else {
				let str = "";
				for (let i1 = 0; i1 < ind; i1++) {
					str = str + "固有値： " + r[i1] + "\n";
					str = str + "  固有ベクトル：";
					for (let i2 = 0; i2 < n; i2++) {
						if (i2 == n-1)
							str = str + v[i1][i2] + "\n";
						else
							str = str + v[i1][i2] + " ";
					}
				}
				document.getElementById("ans").value = str;
			}
		}

		/****************************************/
		/* 最大固有値と固有ベクトル（べき乗法） */
		/*      i : 何番目の固有値かを示す      */
		/*      n : 次数                        */
		/*      m : 丸め誤差の修正回数          */
		/*      ct : 最大繰り返し回数           */
		/*      eps : 収束判定条件              */
		/*      A : 対象とする行列              */
		/*      B : nxnの行列，最初は，単位行列 */
		/*      C : 作業域，nxnの行列           */
		/*      r : 固有値                      */
		/*      v : 各行が固有ベクトル（nxn）   */
		/*      v0 : 固有ベクトルの初期値       */
		/*      v1,v2 : 作業域（n次元ベクトル） */
		/*      return : 求まった固有値の数     */
		/*      coded by Y.Suganuma             */
		/****************************************/
		function power(i, n, m, ct, eps, A, B, C, r, v, v0, v1, v2)
		{
			let l1 = 0.0;
			let l2;
			let x;
			let i1;
			let i2;
			let i3;
			let ind = i;
			let k = 0;
			let k1;
						// 初期設定
			for (i1 = 0; i1 < n; i1++) {
				l1     += v0[i1] * v0[i1];
				v1[i1]  = v0[i1];
			}
			l1 = Math.sqrt(l1);
						// 繰り返し計算
			while (k < ct) {
							// 丸め誤差の修正
				if (k%m == 0) {
					l2 = 0.0;
					for (i1 = 0; i1 < n; i1++) {
						v2[i1] = 0.0;
						for (i2 = 0; i2 < n; i2++)
							v2[i1] += B[i1][i2] * v1[i2];
						l2 += v2[i1] * v2[i1];
					}
					l2 = Math.sqrt(l2);
					for (i1 = 0; i1 < n; i1++)
						v1[i1] = v2[i1] / l2;
				}
							// 次の近似
				l2 = 0.0;
				for (i1 = 0; i1 < n; i1++) {
					v2[i1] = 0.0;
					for (i2 = 0; i2 < n; i2++)
						v2[i1] += A[i1][i2] * v1[i2];
					l2 += v2[i1] * v2[i1];
				}
				l2 = Math.sqrt(l2);
				for (i1 = 0; i1 < n; i1++)
					v2[i1] /= l2;
							// 収束判定
								// 収束した場合
				if (Math.abs((l2-l1)/l1) < eps) {
					k1 = -1;
					for (i1 = 0; i1 < n && k1 < 0; i1++) {
						if (Math.abs(v2[i1]) > 0.001) {
							k1 = i1;
							if (v2[k1]*v1[k1] < 0.0)
								l2 = -l2;
						}
					}
					k    = ct;
					r[i] = l2;
					for (i1 = 0; i1 < n; i1++)
						v[i][i1] = v2[i1];
					if (i == n-1)
						ind = i + 1;
					else {
						for (i1 = 0; i1 < n; i1++) {
							for (i2 = 0; i2 < n; i2++) {
								C[i1][i2] = 0.0;
								for (i3 = 0; i3 < n; i3++) {
									x          = (i1 == i3) ? A[i1][i3] - l2 : A[i1][i3];
									C[i1][i2] += x * B[i3][i2];
								}
							}
						}
						for (i1 = 0; i1 < n; i1++) {
							for (i2 = 0; i2 < n; i2++)
								B[i1][i2] = C[i1][i2];
						}
						for (i1 = 0; i1 < n; i1++) {
							v1[i1] = 0.0;
							for (i2 = 0; i2 < n; i2++)
								v1[i1] += B[i1][i2] * v0[i2];
						}
						for (i1 = 0; i1 < n; i1++)
							v0[i1] = v1[i1];
						ind = power(i+1, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
					}
				}
								// 収束しない場合
				else {
					for (i1 = 0; i1 < n; i1++)
						v1[i1] = v2[i1];
					l1 = l2;
					k++;
				}
			}

			return ind;
		}
	</SCRIPT>

</HEAD>

<BODY STYLE="font-size: 130%; background-color: #eeffee;">

	<H2 STYLE="text-align:center"><B>最大固有値と固有ベクトル（べき乗法）</B></H2>

	<DL>
		<DT>　　テキストフィールドおよびテキストエリアには，例として，以下に示す行列の固有値及び固有ベクトル（固有ベクトルは，正規化されて表示）を求める場合に対する値が設定されています．他の問題を実行する場合は，それらを適切に修正してください．
		<P STYLE="text-align:center"><IMG SRC="power.gif"></P>
	</DL>

	<DIV STYLE="text-align:center">
		次数：<INPUT ID="order" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="3"> 
		修正：<INPUT ID="mod" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="15"> 
		最大繰り返し：<INPUT ID="trial" STYLE="font-size: 100%;" TYPE="text" SIZE="4" VALUE="200"> 
		初期値：<INPUT ID="init" STYLE="font-size: 100%;" TYPE="text" SIZE="20" VALUE="1 1 1"> 
		<BUTTON STYLE="font-size: 100%; background-color: pink" onClick="main()">OK</BUTTON><BR><BR>
		<TEXTAREA ID="ar" COLS="50" ROWS="15" STYLE="font-size: 100%">11 6 -2
-2 18 1
-12 24 13</TEXTAREA><BR><BR>
		<TEXTAREA ID="ans" COLS="80" ROWS="15" STYLE="font-size: 100%"></TEXTAREA>
	</DIV>

</BODY>

</HTML>

  4. PHP 

    <?php

/****************************************/
/* 最大固有値と固有ベクトル（べき乗法） */
/*      coded by Y.Suganuma             */
/****************************************/

					// データの設定
	$ct  = 200;
	$eps = 1.0e-10;
	$n   = 3;
	$m   = 15;

	$A   = array($n);
	$B   = array($n);
	$C   = array($n);
	$v   = array($n);
	for ($i1 = 0; $i1 < $n; $i1++) {
		$A[$i1] = array($n);
		$B[$i1] = array($n);
		$C[$i1] = array($n);
		$v[$i1] = array($n);
	}

	$r  = array($n);
	$v0 = array($n);
	$v1 = array($n);
	$v2 = array($n);

	$A[0][0] = 11.0;
	$A[0][1] = 6.0;
	$A[0][2] = -2.0;

	$A[1][0] = -2.0;
	$A[1][1] = 18.0;
	$A[1][2] = 1.0;

	$A[2][0] = -12.0;
	$A[2][1] = 24.0;
	$A[2][2] = 13.0;

	for ($i1 = 0; $i1 < $n; $i1++) {
		for ($i2 = 0; $i2 < $n; $i2++)
			$B[$i1][$i2] = 0.0;
		$B[$i1][$i1] = 1.0;
		$v0[$i1]     = 1.0;
	}
					// 計算
	$ind = power(0, $n, $m, $ct, $eps, $A, $B, $C, $r, $v, $v0, $v1, $v2);
					// 出力
	if ($ind == 0)
		printf("固有値が求まりませんでした！\n");
	else {
		for ($i1 = 0; $i1 < $ind; $i1++) {
			printf("固有値 %f", $r[$i1]);
			printf(" 固有ベクトル ");
			for ($i2 = 0; $i2 < $n; $i2++)
				printf(" %f", $v[$i1][$i2]);
			printf("\n");
		}
	}

/****************************************/
/* 最大固有値と固有ベクトル（べき乗法） */
/*      i : 何番目の固有値かを示す      */
/*      n : 次数                        */
/*      m : 丸め誤差の修正回数          */
/*      ct : 最大繰り返し回数           */
/*      eps : 収束判定条件              */
/*      A : 対象とする行列              */
/*      B : nxnの行列，最初は，単位行列 */
/*      C : 作業域，nxnの行列           */
/*      r : 固有値                      */
/*      v : 各行が固有ベクトル（nxn）   */
/*      v0 : 固有ベクトルの初期値       */
/*      v1,v2 : 作業域（n次元ベクトル） */
/*      return : 求まった固有値の数     */
/*      coded by Y.Suganuma             */
/****************************************/
function power($i, $n, $m, $ct, $eps, $A, $B, $C, &$r, &$v, $v0, $v1, $v2)
{
					// 初期設定
	$ind = $i;
	$k   = 0;
	$l1  = 0.0;
	for ($i1 = 0; $i1 < $n; $i1++) {
		$l1      += $v0[$i1] * $v0[$i1];
		$v1[$i1]  = $v0[$i1];
	}
	$l1 = sqrt($l1);
					// 繰り返し計算
	while ($k < $ct) {
						// 丸め誤差の修正
		if ($k%$m == 0) {
			$l2 = 0.0;
			for ($i1 = 0; $i1 < $n; $i1++) {
				$v2[$i1] = 0.0;
				for ($i2 = 0; $i2 < $n; $i2++)
					$v2[$i1] += $B[$i1][$i2] * $v1[$i2];
				$l2 += $v2[$i1] * $v2[$i1];
			}
			$l2 = sqrt($l2);
			for ($i1 = 0; $i1 < $n; $i1++)
				$v1[$i1] = $v2[$i1] / $l2;
		}
						// 次の近似
		$l2 = 0.0;
		for ($i1 = 0; $i1 < $n; $i1++) {
			$v2[$i1] = 0.0;
			for ($i2 = 0; $i2 < $n; $i2++)
				$v2[$i1] += $A[$i1][$i2] * $v1[$i2];
			$l2 += $v2[$i1] * $v2[$i1];
		}
		$l2 = sqrt($l2);
		for ($i1 = 0; $i1 < $n; $i1++)
			$v2[$i1] /= $l2;
						// 収束判定
							// 収束した場合
		if (abs(($l2-$l1)/$l1) < $eps) {
			$k1 = -1;
			for ($i1 = 0; $i1 < $n && $k1 < 0; $i1++) {
				if (abs($v2[$i1]) > 0.001) {
					$k1 = $i1;
					if ($v2[$k1]*$v1[$k1] < 0.0)
						$l2 = -$l2;
				}
			}
			$k     = $ct;
			$r[$i] = $l2;
			for ($i1 = 0; $i1 < $n; $i1++)
				$v[$i][$i1] = $v2[$i1];
			if ($i == $n-1)
				$ind = $i + 1;
			else {
				for ($i1 = 0; $i1 < $n; $i1++) {
					for ($i2 = 0; $i2 < $n; $i2++) {
						$C[$i1][$i2] = 0.0;
						for ($i3 = 0; $i3 < $n; $i3++) {
							$x            = ($i1 == $i3) ? $A[$i1][$i3] - $l2 : $A[$i1][$i3];
							$C[$i1][$i2] += $x * $B[$i3][$i2];
						}
					}
				}
				for ($i1 = 0; $i1 < $n; $i1++) {
					for ($i2 = 0; $i2 < $n; $i2++)
						$B[$i1][$i2] = $C[$i1][$i2];
				}
				for ($i1 = 0; $i1 < $n; $i1++) {
					$v1[$i1] = 0.0;
					for ($i2 = 0; $i2 < $n; $i2++)
						$v1[$i1] += $B[$i1][$i2] * $v0[$i2];
				}
				for ($i1 = 0; $i1 < $n; $i1++)
					$v0[$i1] = $v1[$i1];
				$ind = power($i+1, $n, $m, $ct, $eps, $A, $B, $C, $r, $v, $v0, $v1, $v2);
			}
		}
							// 収束しない場合
		else {
			for ($i1 = 0; $i1 < $n; $i1++)
				$v1[$i1] = $v2[$i1];
			$l1 = $l2;
			$k++;
		}
	}

	return $ind;
}

?>

  5. Ruby 

    #***************************************/
# 最大固有値と固有ベクトル（べき乗法） */
#      coded by Y.Suganuma             */
#***************************************/

#***************************************/
# 最大固有値と固有ベクトル（べき乗法） */
#      i : 何番目の固有値かを示す      */
#      n : 次数                        */
#      m : 丸め誤差の修正回数          */
#      ct : 最大繰り返し回数           */
#      eps : 収束判定条件              */
#      a : 対象とする行列              */
#      b : nxnの行列，最初は，単位行列 */
#      c : 作業域，nxnの行列           */
#      r : 固有値                      */
#      v : 各行が固有ベクトル（nxn）   */
#      v0 : 固有ベクトルの初期値       */
#      v1,v2 : 作業域（n次元ベクトル） */
#      return : 求まった固有値の数     */
#      coded by Y.Suganuma             */
#***************************************/
def power(i, n, m, ct, eps, a, b, c, r, v, v0, v1, v2)
					# 初期設定
	ind = i
	k   = 0
	l1  = 0.0
	for i1 in 0 ... n
		l1     += v0[i1] * v0[i1]
		v1[i1]  = v0[i1]
	end
	l1 = Math.sqrt(l1)
					# 繰り返し計算
	while k < ct
						# 丸め誤差の修正
		if k%m == 0
			l2 = 0.0
			for i1 in 0 ... n
				v2[i1] = 0.0
				for i2 in 0 ... n
					v2[i1] += b[i1][i2] * v1[i2]
				end
				l2 += v2[i1] * v2[i1]
			end
			l2 = Math.sqrt(l2)
			for i1 in 0 ... n
				v1[i1] = v2[i1] / l2
			end
		end
						# 次の近似
		l2 = 0.0
		for i1 in 0 ... n
			v2[i1] = 0.0
			for i2 in 0 ... n
				v2[i1] += a[i1][i2] * v1[i2]
			end
			l2 += v2[i1] * v2[i1]
		end
		l2 = Math.sqrt(l2)
		for i1 in 0 ... n
			v2[i1] /= l2
		end
						# 収束判定
							# 収束した場合
		if ((l2-l1)/l1).abs() < eps
			k1 = -1
			for i1 in 0 ... n
				if v2[i1].abs() > 0.001
					k1 = i1
					if v2[k1]*v1[k1] < 0.0
						l2 = -l2
					end
				end
				if k1 >= 0
					break
				end
			end
			k    = ct
			r[i] = l2
			for i1 in 0 ... n
				v[i][i1] = v2[i1]
			end
			if i == n-1
				ind = i + 1
			else
				for i1 in 0 ... n
					for i2 in 0 ... n
						c[i1][i2] = 0.0
						for i3 in 0 ... n
							x          = (i1 == i3) ? a[i1][i3] - l2 : a[i1][i3]
							c[i1][i2] += x * b[i3][i2]
						end
					end
				end
				for i1 in 0 ... n
					for i2 in 0 ... n
						b[i1][i2] = c[i1][i2]
					end
				end
				for i1 in 0 ... n
					v1[i1] = 0.0
					for i2 in 0 ... n
						v1[i1] += b[i1][i2] * v0[i2]
					end
				end
				for i1 in 0 ... n
					v0[i1] = v1[i1]
				end
				ind = power(i+1, n, m, ct, eps, a, b, c, r, v, v0, v1, v2)
			end
							# 収束しない場合
		else
			for i1 in 0 ... n
				v1[i1] = v2[i1]
			end
			l1 = l2
			k += 1
		end
	end

	return ind
end

				# データの設定
ct  = 200
eps = 1.0e-10
n   = 3
m   = 15

a   = Array.new(n)
b   = Array.new(n)
c   = Array.new(n)
v   = Array.new(n)
for i1 in 0 ... n
	a[i1] = Array.new(n)
	b[i1] = Array.new(n)
	c[i1] = Array.new(n)
	v[i1] = Array.new(n)
end

r  = Array.new(n)
v0 = Array.new(n)
v1 = Array.new(n)
v2 = Array.new(n)

a[0][0] = 11.0
a[0][1] = 6.0
a[0][2] = -2.0

a[1][0] = -2.0
a[1][1] = 18.0
a[1][2] = 1.0

a[2][0] = -12.0
a[2][1] = 24.0
a[2][2] = 13.0

for i1 in 0 ... n
	for i2 in 0 ... n
		b[i1][i2] = 0.0
	end
	b[i1][i1] = 1.0
	v0[i1]    = 1.0
end
				# 計算
ind = power(0, n, m, ct, eps, a, b, c, r, v, v0, v1, v2)
				# 出力
if ind == 0
	printf("固有値が求まりませんでした！\n")
else
	for i1 in 0 ... ind
		printf("固有値 %f", r[i1])
		printf(" 固有ベクトル ")
		for i2 in 0 ... n
			printf(" %f", v[i1][i2])
		end
		printf("\n")
	end
end

  6. Python 

    # -*- coding: UTF-8 -*-
import numpy as np
from math import *

############################################
# 最大固有値と固有ベクトル（べき乗法）
#      i : 何番目の固有値かを示す
#      n : 次数
#      m : 丸め誤差の修正回数
#      ct : 最大繰り返し回数
#      eps : 収束判定条件
#      A : 対象とする行列
#      B : nxnの行列，最初は，単位行列
#      C : 作業域，nxnの行列
#      r : 固有値
#      v : 各行が固有ベクトル（nxn）
#      v0 : 固有ベクトルの初期値
#      v1,v2 : 作業域（n次元ベクトル）
#      return : 求まった固有値の数
#      coded by Y.Suganuma
############################################

def power(i, n, m, ct, eps, A, B, C, r, v, v0, v1, v2) :
			# 初期設定
	ind = i
	k   = 0
	l1  = 0.0
	for i1 in range(0, n) :
		l1     += v0[i1] * v0[i1]
		v1[i1]  = v0[i1]
	l1 = sqrt(l1)
			# 繰り返し計算
	while k < ct :
				# 丸め誤差の修正
		if k%m == 0 :
			l2 = 0.0
			for i1 in range(0, n) :
				v2[i1] = 0.0
				for i2 in range(0, n) :
					v2[i1] += B[i1][i2] * v1[i2]
				l2 += v2[i1] * v2[i1]
			l2 = sqrt(l2)
			for i1 in range(0, n) :
				v1[i1] = v2[i1] / l2
				# 次の近似
		l2 = 0.0
		for i1 in range(0, n) :
			v2[i1] = 0.0
			for i2 in range(0, n) :
				v2[i1] += A[i1][i2] * v1[i2]
			l2 += v2[i1] * v2[i1]
		l2 = sqrt(l2)
		for i1 in range(0, n) :
			v2[i1] /= l2
				# 収束判定
					# 収束した場合
		if abs((l2-l1)/l1) < eps :
			k1 = -1
			for i1 in range(0, n) :
				if abs(v2[i1]) > 0.001 :
					k1 = i1
					if v2[k1]*v1[k1] < 0.0 :
						l2 = -l2
					break
			k    = ct
			r[i] = l2
			for i1 in range(0, n) :
				v[i][i1] = v2[i1]
			if i == n-1 :
				ind = i + 1
			else :
				for i1 in range(0, n) :
					for i2 in range(0, n) :
						C[i1][i2] = 0.0
						for i3 in range(0, n) :
							if i1 == i3 :
								x = A[i1][i3] - l2
							else :
								x = A[i1][i3]
							C[i1][i2] += x * B[i3][i2]
				for i1 in range(0, n) :
					for i2 in range(0, n) :
						B[i1][i2] = C[i1][i2]
				for i1 in range(0, n) :
					v1[i1] = 0.0
					for i2 in range(0, n) :
						v1[i1] += B[i1][i2] * v0[i2]
				for i1 in range(0, n) :
					v0[i1] = v1[i1]
				ind = power(i+1, n, m, ct, eps, A, B, C, r, v, v0, v1, v2)
					# 収束しない場合
		else :
			for i1 in range(0, n) :
				v1[i1] = v2[i1]
			l1 = l2
			k += 1

	return ind

############################################
# 最大固有値と固有ベクトル（べき乗法）
#      coded by Y.Suganuma
############################################
			# データの設定
ct  = 200
eps = 1.0e-10
n   = 3
m   = 15
A   = np.array([[11, 6, -2],[-2, 18, 1],[-12, 24, 13]], np.float)
B   = np.array([[1, 0, 0],[0, 1, 0],[0, 0, 1]], np.float)
C   = np.empty((n, n), np.float)
v   = np.empty((n, n), np.float)
r   = np.empty(n, np.float)
v0  = np.array([1, 1, 1], np.float)
v1  = np.empty(n, np.float)
v2  = np.empty(n, np.float)
			# 計算
ind = power(0, n, m, ct, eps, A, B, C, r, v, v0, v1, v2)
			# 出力
if ind == 0 :
	print("固有値が求まりませんでした！")
else :
	for i1 in range(0, ind) :
		print("固有値 " + str(r[i1]), end="")
		print(" 固有ベクトル ", end="")
		for i2 in range(0, n) :
			print(" " + str(v[i1][i2]), end="")
		print()

  7. C# 

    /****************************************/
/* 最大固有値と固有ベクトル（べき乗法） */
/*      coded by Y.Suganuma             */
/****************************************/
using System;

class Program
{
	static void Main()
	{
		Test1 ts = new Test1();
	}
}

class Test1
{
	public Test1()
	{
					// データの設定
		int ct       = 200;
		int n        = 3;
		int m        = 15;
		double eps   = 1.0e-10;
		double[][] A = new double [][] {
			new double[] {11.0, 6.0, -2.0},
			new double[] {-2.0, 18.0, 1.0},
			new double[] {-12.0, 24.0, 13.0}
		};
		double[][] B = new double [n][];
		B[0] = new double[] {1.0, 0.0, 0.0};
		B[1] = new double[] {0.0, 1.0, 0.0};
		B[2] = new double[] {0.0, 0.0, 1.0};
		double[][] C = new double [n][];
		double[][] v = new double [n][];
		for (int i1 = 0; i1 < n; i1++) {
			C[i1] = new double [n];
			v[i1] = new double [n];
		}
		double[] v0 = {1.0, 1.0, 1.0};
		double[] v1 = new double [n];
		double[] v2 = new double [n];
		double[] r  = new double [n];
					// 計算
		int ind = power(0, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
					// 出力
		if (ind == 0)
			Console.WriteLine("固有値が求まりませんでした！");
		else {
			for (int i1 = 0; i1 < ind; i1++) {
				Console.WriteLine("固有値 " + r[i1]);
				Console.Write("   固有ベクトル ");
				for (int i2 = 0; i2 < n; i2++)
					Console.Write(" " + v[i1][i2]);
				Console.WriteLine();
			}
		}
	}

	/****************************************/
	/* 最大固有値と固有ベクトル（べき乗法） */
	/*      i : 何番目の固有値かを示す      */
	/*      n : 次数                        */
	/*      m : 丸め誤差の修正回数          */
	/*      ct : 最大繰り返し回数           */
	/*      eps : 収束判定条件              */
	/*      A : 対象とする行列              */
	/*      B : nxnの行列，最初は，単位行列 */
	/*      C : 作業域，nxnの行列           */
	/*      r : 固有値                      */
	/*      v : 各行が固有ベクトル（nxn）   */
	/*      v0 : 固有ベクトルの初期値       */
	/*      v1,v2 : 作業域（n次元ベクトル） */
	/*      return : 求まった固有値の数     */
	/*      coded by Y.Suganuma             */
	/****************************************/
	static int power(int i, int n, int m, int ct, double eps, double[][] A,
	                 double[][] B, double[][] C, double[] r, double[][] v,
	                 double[] v0, double[] v1, double[] v2)
	{
					// 初期設定
		int ind   = i;
		int k     = 0;
		double l1 = 0.0;
		for (int i1 = 0; i1 < n; i1++) {
			l1     += v0[i1] * v0[i1];
			v1[i1]  = v0[i1];
		}
		l1 = Math.Sqrt(l1);
					// 繰り返し計算
		while (k < ct) {
						// 丸め誤差の修正
			double l2;
			if (k%m == 0) {
				l2 = 0.0;
				for (int i1 = 0; i1 < n; i1++) {
					v2[i1] = 0.0;
					for (int i2 = 0; i2 < n; i2++)
						v2[i1] += B[i1][i2] * v1[i2];
					l2 += v2[i1] * v2[i1];
				}
				l2 = Math.Sqrt(l2);
				for (int i1 = 0; i1 < n; i1++)
					v1[i1] = v2[i1] / l2;
			}
						// 次の近似
			l2 = 0.0;
			for (int i1 = 0; i1 < n; i1++) {
				v2[i1] = 0.0;
				for (int i2 = 0; i2 < n; i2++)
					v2[i1] += A[i1][i2] * v1[i2];
				l2 += v2[i1] * v2[i1];
			}
			l2 = Math.Sqrt(l2);
			for (int i1 = 0; i1 < n; i1++)
				v2[i1] /= l2;
						// 収束判定
							// 収束した場合
			if (Math.Abs((l2-l1)/l1) < eps) {
				int k1 = -1;
				for (int i1 = 0; i1 < n && k1 < 0; i1++) {
					if (Math.Abs(v2[i1]) > 0.001) {
						k1 = i1;
						if (v2[k1]*v1[k1] < 0.0)
							l2 = -l2;
					}
				}
				k    = ct;
				r[i] = l2;
				for (int i1 = 0; i1 < n; i1++)
					v[i][i1] = v2[i1];
				if (i == n-1)
					ind = i + 1;
				else {
					for (int i1 = 0; i1 < n; i1++) {
						for (int i2 = 0; i2 < n; i2++) {
							C[i1][i2] = 0.0;
							for (int i3 = 0; i3 < n; i3++) {
								double x   = (i1 == i3) ? A[i1][i3] - l2 : A[i1][i3];
								C[i1][i2] += x * B[i3][i2];
							}
						}
					}
					for (int i1 = 0; i1 < n; i1++) {
						for (int i2 = 0; i2 < n; i2++)
							B[i1][i2] = C[i1][i2];
					}
					for (int i1 = 0; i1 < n; i1++) {
						v1[i1] = 0.0;
						for (int i2 = 0; i2 < n; i2++)
							v1[i1] += B[i1][i2] * v0[i2];
					}
					for (int i1 = 0; i1 < n; i1++)
						v0[i1] = v1[i1];
					ind = power(i+1, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
				}
			}
							// 収束しない場合
			else {
				for (int i1 = 0; i1 < n; i1++)
					v1[i1] = v2[i1];
				l1 = l2;
				k++;
			}
		}

		return ind;
	}
}

  8. VB 

    ''''''''''''''''''''''''''''''''''''''''
' 最大固有値と固有ベクトル（べき乗法） '
'      coded by Y.Suganuma             '
''''''''''''''''''''''''''''''''''''''''
Module Test

	Sub Main()
					' データの設定
		Dim ct As Integer  = 200
		Dim n As Integer   = 3
		Dim m As Integer   = 15
		Dim eps As Double  = 1.0e-10
		Dim A(,) As Double = {{11.0, 6.0, -2.0}, {-2.0, 18.0, 1.0}, {-12.0, 24.0, 13.0}}
		Dim B(,) As Double = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}
		Dim C(n,n) As Double
		Dim v(n,n) As Double
		Dim v0() As Double = {1.0, 1.0, 1.0}
		Dim v1(n) As Double
		Dim v2(n) As Double
		Dim r(n) As Double
					' 計算
		Dim ind As Integer = power(0, n, m, ct, eps, A, B, C, r, v, v0, v1, v2)
					' 出力
		If ind = 0
			Console.WriteLine("固有値が求まりませんでした！")
		Else
			For i1 As Integer = 0 To ind-1
				Console.WriteLine("固有値 " & r(i1))
				Console.Write("   固有ベクトル ")
				For i2 As Integer = 0 To n-1
					Console.Write(" " & v(i1,i2))
				Next
				Console.WriteLine()
			Next
		End If

	End Sub

	''''''''''''''''''''''''''''''''''''''''
	' 最大固有値と固有ベクトル（べき乗法） '
	'      i : 何番目の固有値かを示す      '
	'      n : 次数                        '
	'      m : 丸め誤差の修正回数          '
	'      ct : 最大繰り返し回数           '
	'      eps : 収束判定条件              '
	'      A : 対象とする行列              '
	'      B : nxnの行列，最初は，単位行列 '
	'      C : 作業域，nxnの行列           '
	'      r : 固有値                      '
	'      v : 各行が固有ベクトル（nxn）   '
	'      v0 : 固有ベクトルの初期値       '
	'      v1,v2 : 作業域（n次元ベクトル） '
	'      return : 求まった固有値の数     '
	'      coded by Y.Suganuma             '
	''''''''''''''''''''''''''''''''''''''''
	Function power(i As Integer, n As Integer, m As Integer, ct As Integer,
	               eps As Double, A(,) As Double, B(,) As Double, C(,) As Double,
	               r() As Double, v(,) As Double, v0() As Double, v1() As Double,
	               v2() As Double)
					' 初期設定
		Dim ind As Integer = i
		Dim k As Integer   = 0
		Dim l1 As Double   = 0.0
		For i1 As Integer = 0 To n-1
			l1     += v0(i1) * v0(i1)
			v1(i1)  = v0(i1)
		Next
		l1 = Math.Sqrt(l1)
					' 繰り返し計算
		Do While k < ct
						' 丸め誤差の修正
			Dim l2 As Double
			If (k Mod m) = 0
				l2 = 0.0
				For i1 As Integer = 0 To n-1
					v2(i1) = 0.0
					For i2 As Integer = 0 To n-1
						v2(i1) += B(i1,i2) * v1(i2)
					Next
					l2 += v2(i1) * v2(i1)
				Next
				l2 = Math.Sqrt(l2)
				For i1 As Integer = 0 To n-1
					v1(i1) = v2(i1) / l2
				Next
			End If
						' 次の近似
			l2 = 0.0
			For i1 As Integer = 0 To n-1
				v2(i1) = 0.0
				For i2 As Integer = 0 To n-1
					v2(i1) += A(i1,i2) * v1(i2)
				Next
				l2 += v2(i1) * v2(i1)
			Next
			l2 = Math.Sqrt(l2)
			For i1 As Integer = 0 To n-1
				v2(i1) /= l2
			Next
						' 収束判定
							' 収束した場合
			If Math.Abs((l2-l1)/l1) < eps
				Dim k1 As Integer = -1
				Dim ii As Integer = 0
				Do While ii < n and k1 < 0
					If Math.Abs(v2(ii)) > 0.001
						k1 = ii
						If v2(k1)*v1(k1) < 0.0
							l2 = -l2
						End If
					End If
					ii += 1
				Loop
				k    = ct
				r(i) = l2
				For i1 As Integer = 0 To n-1
					v(i,i1) = v2(i1)
				Next
				If i = n-1
					ind = i + 1
				Else
					For i1 As Integer = 0 To n-1
						For i2 As Integer = 0 To n-1
							C(i1,i2) = 0.0
							For i3 As Integer = 0 To n-1
								Dim x As Double = A(i1,i3) - l2
								If i1 <> i3
									x = A(i1,i3)
								End If
								C(i1,i2) += x * B(i3,i2)
							Next
						Next
					Next
					For i1 As Integer = 0 To n-1
						For i2 As Integer = 0 To n-1
							B(i1,i2) = C(i1,i2)
						Next
					Next
					For i1 As Integer = 0 To n-1
						v1(i1) = 0.0
						For i2 As Integer = 0 To n-1
							v1(i1) += B(i1,i2) * v0(i2)
						Next
					Next
					For i1 As Integer = 0 To n-1
						v0(i1) = v1(i1)
					Next
					ind = power(i+1, n, m, ct, eps, A, B, C, r, v, v0, v1, v2)
				End If
							' 収束しない場合
			Else
				For i1 As Integer = 0 To n-1
					v1(i1) = v2(i1)
				Next
				l1 = l2
				k += 1
			End If
		Loop

		Return ind

	End Function

End Module

情報学部 菅沼ホーム 目次 索引