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

<?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;
}

?>