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

実対称行列の固有値・固有ベクトル（ヤコビ法）

    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 Jacobi(int, int, double, double **, double **, double **, double **, double **);

int main()
{
	double **A, **A1, **A2, **X1, **X2, eps;
	int i1, i2, ind, ct, n;
					// データの設定
	ct   = 1000;
	eps  = 1.0e-10;
	n    = 3;
	A    = new double * [n];
	A1   = new double * [n];
	A2   = new double * [n];
	X1   = new double * [n];
	X2   = new double * [n];
	for (i1 = 0; i1 < n; i1++) {
		A[i1]  = new double [n];
		A1[i1] = new double [n];
		A2[i1] = new double [n];
		X1[i1] = new double [n];
		X2[i1] = new double [n];
	}

	A[0][0] = 1.0;
	A[0][1] = 0.0;
	A[0][2] = -1.0;

	A[1][0] = 0.0;
	A[1][1] = 1.0;
	A[1][2] = -1.0;

	A[2][0] = -1.0;
	A[2][1] = -1.0;
	A[2][2] = 2.0;
					// 計算
	ind = Jacobi(n, ct, eps, A, A1, A2, X1, X2);
					// 出力
	if (ind > 0)
		printf("収束しませんでした！\n");
	else {
		printf("固有値\n");
		for (i1 = 0; i1 < n; i1++)
			printf(" %f", A1[i1][i1]);
		printf("\n固有ベクトル（各列が固有ベクトル）\n");
		for (i1 = 0; i1 < n; i1++) {
			for (i2 = 0; i2 < n; i2++)
				printf(" %f", X1[i1][i2]);
			printf("\n");
		}
	}

	for (i1 = 0; i1 < n; i1++) {
		delete [] A[i1];
		delete [] A1[i1];
		delete [] A2[i1];
		delete [] X1[i1];
		delete [] X2[i1];
	}
	delete [] A;
	delete [] A1;
	delete [] A2;
	delete [] X1;
	delete [] X2;

	return 0;
}

/*************************************************************/
/* 実対称行列の固有値・固有ベクトル（ヤコビ法）              */
/*      n : 次数                                             */
/*      ct : 最大繰り返し回数                                */
/*      eps : 収束判定条件                                   */
/*      A : 対象とする行列                                   */
/*      A1, A2 : 作業域（nxnの行列），A1の対角要素が固有値   */
/*      X1, X2 : 作業域（nxnの行列），X1の各列が固有ベクトル */
/*      return : =0 : 正常                                   */
/*               =1 : 収束せず                               */
/*      coded by Y.Suganuma                                  */
/*************************************************************/
#include <math.h>

int Jacobi(int n, int ct, double eps, double **A, double **A1, double **A2,
           double **X1, double **X2)
{
	double max, s, t, v, sn, cs;
	int i1, i2, k = 0, ind = 1, p = 0, q = 0;
					// 初期設定
	for (i1 = 0; i1 < n; i1++) {
		for (i2 = 0; i2 < n; i2++) {
			A1[i1][i2] = A[i1][i2];
			X1[i1][i2] = 0.0;
		}
		X1[i1][i1] = 1.0;
	}
					// 計算
	while (ind > 0 && k < ct) {
						// 最大要素の探索
		max = 0.0;
		for (i1 = 0; i1 < n; i1++) {
			for (i2 = 0; i2 < n; i2++) {
				if (i2 != i1) {
					if (fabs(A1[i1][i2]) > max) {
						max = fabs(A1[i1][i2]);
						p   = i1;
						q   = i2;
					}
				}
			}
		}
						// 収束判定
							// 収束した
		if (max < eps)
			ind = 0;
							// 収束しない
		else {
								// 準備
			s  = -A1[p][q];
			t  = 0.5 * (A1[p][p] - A1[q][q]);
			v  = fabs(t) / sqrt(s * s + t * t);
			sn = sqrt(0.5 * (1.0 - v));
			if (s*t < 0.0)
				sn = -sn;
			cs = sqrt(1.0 - sn * sn);
								// Akの計算
			for (i1 = 0; i1 < n; i1++) {
				if (i1 == p) {
					for (i2 = 0; i2 < n; i2++) {
						if (i2 == p)
							A2[p][p] = A1[p][p] * cs * cs + A1[q][q] * sn * sn -
                                       2.0 * A1[p][q] * sn * cs;
						else if (i2 == q)
							A2[p][q] = 0.0;
						else
							A2[p][i2] = A1[p][i2] * cs - A1[q][i2] * sn;
					}
				}
				else if (i1 == q) {
					for (i2 = 0; i2 < n; i2++) {
						if (i2 == q)
							A2[q][q] = A1[p][p] * sn * sn + A1[q][q] * cs * cs +
                                       2.0 * A1[p][q] * sn * cs;
						else if (i2 == p)
							A2[q][p] = 0.0;
						else
							A2[q][i2] = A1[q][i2] * cs + A1[p][i2] * sn;
					}
				}
				else {
					for (i2 = 0; i2 < n; i2++) {
						if (i2 == p)
							A2[i1][p] = A1[i1][p] * cs - A1[i1][q] * sn;
						else if (i2 == q)
							A2[i1][q] = A1[i1][q] * cs + A1[i1][p] * sn;
						else
							A2[i1][i2] = A1[i1][i2];
					}
				}
			}
								// Xkの計算
			for (i1 = 0; i1 < n; i1++) {
				for (i2 = 0; i2 < n; i2++) {
					if (i2 == p)
						X2[i1][p] = X1[i1][p] * cs - X1[i1][q] * sn;
					else if (i2 == q)
						X2[i1][q] = X1[i1][q] * cs + X1[i1][p] * sn;
					else
						X2[i1][i2] = X1[i1][i2];
				}
			}
								// 次のステップへ
			k++;
			for (i1 = 0; i1 < n; i1++) {
				for (i2 = 0; i2 < n; i2++) {
					A1[i1][i2] = A2[i1][i2];
					X1[i1][i2] = X2[i1][i2];
				}
			}
		}
	}

	return ind;
}

  2. Java 

    /************************************************/
/* 実対称行列の固有値・固有ベクトル（ヤコビ法） */
/*      coded by Y.Suganuma                     */
/************************************************/
import java.io.*;

public class Test {
	public static void main(String args[]) throws IOException
	{
		double A[][], A1[][], A2[][], X1[][], X2[][], eps;
		int i1, i2, ind, ct, n;
					// データの設定
		ct   = 1000;
		eps  = 1.0e-10;
		n    = 3;
		A    = new double [n][n];
		A1   = new double [n][n];
		A2   = new double [n][n];
		X1   = new double [n][n];
		X2   = new double [n][n];

		A[0][0] = 1.0;
		A[0][1] = 0.0;
		A[0][2] = -1.0;

		A[1][0] = 0.0;
		A[1][1] = 1.0;
		A[1][2] = -1.0;

		A[2][0] = -1.0;
		A[2][1] = -1.0;
		A[2][2] = 2.0;
					// 計算
		ind = Jacobi(n, ct, eps, A, A1, A2, X1, X2);
					// 出力
		if (ind > 0)
			System.out.println("収束しませんでした！");
		else {
			System.out.println("固有値");
			for (i1 = 0; i1 < n; i1++)
				System.out.print(" " + A1[i1][i1]);
			System.out.println("\n固有ベクトル（各列が固有ベクトル）");
			for (i1 = 0; i1 < n; i1++) {
				for (i2 = 0; i2 < n; i2++)
					System.out.print(" " + X1[i1][i2]);
				System.out.println();
			}
		}
	}

	/*************************************************************/
	/* 実対称行列の固有値・固有ベクトル（ヤコビ法）              */
	/*      n : 次数                                             */
	/*      ct : 最大繰り返し回数                                */
	/*      eps : 収束判定条件                                   */
	/*      A : 対象とする行列                                   */
	/*      A1, A2 : 作業域（nxnの行列），A1の対角要素が固有値   */
	/*      X1, X2 : 作業域（nxnの行列），X1の各列が固有ベクトル */
	/*      return : =0 : 正常                                   */
	/*               =1 : 収束せず                               */
	/*      coded by Y.Suganuma                                  */
	/*************************************************************/
	static int Jacobi(int n, int ct, double eps, double A[][], double A1[][], double A2[][],
                  double X1[][], double X2[][])
	{
		double max, s, t, v, sn, cs;
		int i1, i2, k = 0, ind = 1, p = 0, q = 0;
					// 初期設定
		for (i1 = 0; i1 < n; i1++) {
			for (i2 = 0; i2 < n; i2++) {
				A1[i1][i2] = A[i1][i2];
				X1[i1][i2] = 0.0;
			}
			X1[i1][i1] = 1.0;
		}
					// 計算
		while (ind > 0 && k < ct) {
						// 最大要素の探索
			max = 0.0;
			for (i1 = 0; i1 < n; i1++) {
				for (i2 = 0; i2 < n; i2++) {
					if (i2 != i1) {
						if (Math.abs(A1[i1][i2]) > max) {
							max = Math.abs(A1[i1][i2]);
							p   = i1;
							q   = i2;
						}
					}
				}
			}
						// 収束判定
							// 収束した
			if (max < eps)
				ind = 0;
							// 収束しない
			else {
								// 準備
				s  = -A1[p][q];
				t  = 0.5 * (A1[p][p] - A1[q][q]);
				v  = Math.abs(t) / Math.sqrt(s * s + t * t);
				sn = Math.sqrt(0.5 * (1.0 - v));
				if (s*t < 0.0)
					sn = -sn;
				cs = Math.sqrt(1.0 - sn * sn);
								// Akの計算
				for (i1 = 0; i1 < n; i1++) {
					if (i1 == p) {
						for (i2 = 0; i2 < n; i2++) {
							if (i2 == p)
								A2[p][p] = A1[p][p] * cs * cs + A1[q][q] * sn * sn -
                                           2.0 * A1[p][q] * sn * cs;
							else if (i2 == q)
								A2[p][q] = 0.0;
							else
								A2[p][i2] = A1[p][i2] * cs - A1[q][i2] * sn;
						}
					}
					else if (i1 == q) {
						for (i2 = 0; i2 < n; i2++) {
							if (i2 == q)
								A2[q][q] = A1[p][p] * sn * sn + A1[q][q] * cs * cs +
                                           2.0 * A1[p][q] * sn * cs;
							else if (i2 == p)
								A2[q][p] = 0.0;
							else
								A2[q][i2] = A1[q][i2] * cs + A1[p][i2] * sn;
						}
					}
					else {
						for (i2 = 0; i2 < n; i2++) {
							if (i2 == p)
								A2[i1][p] = A1[i1][p] * cs - A1[i1][q] * sn;
							else if (i2 == q)
								A2[i1][q] = A1[i1][q] * cs + A1[i1][p] * sn;
							else
								A2[i1][i2] = A1[i1][i2];
						}
					}
				}
								// Xkの計算
				for (i1 = 0; i1 < n; i1++) {
					for (i2 = 0; i2 < n; i2++) {
						if (i2 == p)
							X2[i1][p] = X1[i1][p] * cs - X1[i1][q] * sn;
						else if (i2 == q)
							X2[i1][q] = X1[i1][q] * cs + X1[i1][p] * sn;
						else
							X2[i1][i2] = X1[i1][i2];
					}
				}
								// 次のステップへ
				k++;
				for (i1 = 0; i1 < n; i1++) {
					for (i2 = 0; i2 < n; i2++) {
						A1[i1][i2] = A2[i1][i2];
						X1[i1][i2] = X2[i1][i2];
					}
				}
			}
		}

		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 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 A1 = new Array();
			for (let i1 = 0; i1 < n; i1++)
				A1[i1] = new Array();
			let A2 = new Array();
			for (let i1 = 0; i1 < n; i1++)
				A2[i1] = new Array();
			let X1 = new Array();
			for (let i1 = 0; i1 < n; i1++)
				X1[i1] = new Array();
			let X2 = new Array();
			for (let i1 = 0; i1 < n; i1++)
				X2[i1] = new Array();

			ind = jacobi(n, ct, eps, A, A1, A2, X1, X2);
					// 出力
			if (ind > 0)
				document.getElementById("ans").value = "解を求めることができません！";
			else {
				let str = "固有値：\n";
				for (let i1 = 0; i1 < n; i1++) {
					if (i1 == n-1)
						str = str + A1[i1][i1] + "\n";
					else
						str = str + A1[i1][i1] + " ";
				}
				str = str + "固有ベクトル：\n";
				for (let i1 = 0; i1 < n; i1++) {
					str = str + "   ";
					for (let i2 = 0; i2 < n; i2++) {
						if (i2 == n-1)
							str = str + X1[i1][i2] + "\n";
						else
							str = str + X1[i1][i2] + " ";
					}
				}
				document.getElementById("ans").value = str;
			}
		}

		/*************************************************************/
		/* 実対称行列の固有値・固有ベクトル（ヤコビ法）              */
		/*      n : 次数                                             */
		/*      ct : 最大繰り返し回数                                */
		/*      eps : 収束判定条件                                   */
		/*      A : 対象とする行列                                   */
		/*      A1, A2 : 作業域（nxnの行列），A1の対角要素が固有値   */
		/*      X1, X2 : 作業域（nxnの行列），X1の各列が固有ベクトル */
		/*      return : =0 : 正常                                   */
		/*               =1 : 収束せず                               */
		/*      coded by Y.Suganuma                                  */
		/*************************************************************/
		function jacobi(n, ct, eps, A, A1, A2, X1, X2)
		{
			let max;
			let s;
			let t;
			let v;
			let sn;
			let cs;
			let i1;
			let i2;
			let k = 0;
			let ind = 1;
			let p = 0;
			let q = 0;
						// 初期設定
			for (i1 = 0; i1 < n; i1++) {
				for (i2 = 0; i2 < n; i2++) {
					A1[i1][i2] = A[i1][i2];
					X1[i1][i2] = 0.0;
				}
				X1[i1][i1] = 1.0;
			}
						// 計算
			while (ind > 0 && k < ct) {
							// 最大要素の探索
				max = 0.0;
				for (i1 = 0; i1 < n; i1++) {
					for (i2 = 0; i2 < n; i2++) {
						if (i2 != i1) {
							if (Math.abs(A1[i1][i2]) > max) {
								max = Math.abs(A1[i1][i2]);
								p   = i1;
								q   = i2;
							}
						}
					}
				}
							// 収束判定
								// 収束した
				if (max < eps)
					ind = 0;
								// 収束しない
				else {
									// 準備
					s  = -A1[p][q];
					t  = 0.5 * (A1[p][p] - A1[q][q]);
					v  = Math.abs(t) / Math.sqrt(s * s + t * t);
					sn = Math.sqrt(0.5 * (1.0 - v));
					if (s*t < 0.0)
						sn = -sn;
					cs = Math.sqrt(1.0 - sn * sn);
									// Akの計算
					for (i1 = 0; i1 < n; i1++) {
						if (i1 == p) {
							for (i2 = 0; i2 < n; i2++) {
								if (i2 == p)
									A2[p][p] = A1[p][p] * cs * cs + A1[q][q] * sn * sn -
                    	                       2.0 * A1[p][q] * sn * cs;
								else if (i2 == q)
									A2[p][q] = 0.0;
								else
									A2[p][i2] = A1[p][i2] * cs - A1[q][i2] * sn;
							}
						}
						else if (i1 == q) {
							for (i2 = 0; i2 < n; i2++) {
								if (i2 == q)
									A2[q][q] = A1[p][p] * sn * sn + A1[q][q] * cs * cs +
                                	           2.0 * A1[p][q] * sn * cs;
								else if (i2 == p)
									A2[q][p] = 0.0;
								else
									A2[q][i2] = A1[q][i2] * cs + A1[p][i2] * sn;
							}
						}
						else {
							for (i2 = 0; i2 < n; i2++) {
								if (i2 == p)
									A2[i1][p] = A1[i1][p] * cs - A1[i1][q] * sn;
								else if (i2 == q)
									A2[i1][q] = A1[i1][q] * cs + A1[i1][p] * sn;
								else
									A2[i1][i2] = A1[i1][i2];
							}
						}
					}
									// Xkの計算
					for (i1 = 0; i1 < n; i1++) {
						for (i2 = 0; i2 < n; i2++) {
							if (i2 == p)
								X2[i1][p] = X1[i1][p] * cs - X1[i1][q] * sn;
							else if (i2 == q)
								X2[i1][q] = X1[i1][q] * cs + X1[i1][p] * sn;
							else
								X2[i1][i2] = X1[i1][i2];
						}
					}
									// 次のステップへ
					k++;
					for (i1 = 0; i1 < n; i1++) {
						for (i2 = 0; i2 < n; i2++) {
							A1[i1][i2] = A2[i1][i2];
							X1[i1][i2] = X2[i1][i2];
						}
					}
				}
			}

			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="jacobi.gif"></P>
	</DL>

	<DIV STYLE="text-align:center">
		次数：<INPUT ID="order" STYLE="font-size: 100%" TYPE="text" SIZE="2" VALUE="3"> 
		最大繰り返し回数：<INPUT ID="trial" STYLE="font-size: 100%;" TYPE="text" SIZE="4" VALUE="1000"> 
		<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%">1 0 -1
0 1 -1
-1 -1 2</TEXTAREA><BR><BR>
		<TEXTAREA ID="ans" COLS="60" ROWS="15" STYLE="font-size: 100%"></TEXTAREA>
	</DIV>

</BODY>

</HTML>

  4. PHP 

    <?php

/************************************************/
/* 実対称行列の固有値・固有ベクトル（ヤコビ法） */
/*      coded by Y.Suganuma                     */
/************************************************/

					// データの設定
	$ct   = 1000;
	$eps  = 1.0e-10;
	$n    = 3;
	$A    = array($n);
	$A1   = array($n);
	$A2   = array($n);
	$X1   = array($n);
	$X2   = array($n);
	for ($i1 = 0; $i1 < $n; $i1++) {
		$A[$i1]  = array($n);
		$A1[$i1] = array($n);
		$A2[$i1] = array($n);
		$X1[$i1] = array($n);
		$X2[$i1] = array($n);
	}

	$A[0][0] = 1.0;
	$A[0][1] = 0.0;
	$A[0][2] = -1.0;

	$A[1][0] = 0.0;
	$A[1][1] = 1.0;
	$A[1][2] = -1.0;

	$A[2][0] = -1.0;
	$A[2][1] = -1.0;
	$A[2][2] = 2.0;
					// 計算
	$ind = Jacobi($n, $ct, $eps, $A, $A1, $A2, $X1, $X2);
					// 出力
	if ($ind > 0)
		printf("収束しませんでした！\n");
	else {
		printf("固有値\n");
		for ($i1 = 0; $i1 < $n; $i1++)
			printf(" %f", $A1[$i1][$i1]);
		printf("\n固有ベクトル（各列が固有ベクトル）\n");
		for ($i1 = 0; $i1 < $n; $i1++) {
			for ($i2 = 0; $i2 < $n; $i2++)
				printf(" %f", $X1[$i1][$i2]);
			printf("\n");
		}
	}

/*************************************************************/
/* 実対称行列の固有値・固有ベクトル（ヤコビ法）              */
/*      n : 次数                                             */
/*      ct : 最大繰り返し回数                                */
/*      eps : 収束判定条件                                   */
/*      A : 対象とする行列                                   */
/*      A1, A2 : 作業域（nxnの行列），A1の対角要素が固有値   */
/*      X1, X2 : 作業域（nxnの行列），X1の各列が固有ベクトル */
/*      return : =0 : 正常                                   */
/*               =1 : 収束せず                               */
/*      coded by Y.Suganuma                                  */
/*************************************************************/
function Jacobi($n, $ct, $eps, $A, &$A1, $A2, &$X1, $X2)
{
					// 初期設定
	$k   = 0;
	$ind = 1;
	$p   = 0;
	$q   = 0;

	for ($i1 = 0; $i1 < $n; $i1++) {
		for ($i2 = 0; $i2 < $n; $i2++) {
			$A1[$i1][$i2] = $A[$i1][$i2];
			$X1[$i1][$i2] = 0.0;
		}
		$X1[$i1][$i1] = 1.0; 
	}
					// 計算
	while ($ind > 0 && $k < $ct) {
						// 最大要素の探索
		$max = 0.0;
		for ($i1 = 0; $i1 < $n; $i1++) {
			for ($i2 = 0; $i2 < $n; $i2++) {
				if ($i2 != $i1) {
					if (abs($A1[$i1][$i2]) > $max) {
						$max = abs($A1[$i1][$i2]);
						$p   = $i1;
						$q   = $i2;
					}
				}
			}
		}
						// 収束判定
							// 収束した
		if ($max < $eps)
			$ind = 0;
							// 収束しない
		else {
								// 準備
			$s  = -$A1[$p][$q];
			$t  = 0.5 * ($A1[$p][$p] - $A1[$q][$q]);
			$v  = abs($t) / sqrt($s * $s + $t * $t);
			$sn = sqrt(0.5 * (1.0 - $v));
			if ($s*$t < 0.0)
				$sn = -$sn;
			$cs = sqrt(1.0 - $sn * $sn);
								// Akの計算
			for ($i1 = 0; $i1 < $n; $i1++) {
				if ($i1 == $p) {
					for ($i2 = 0; $i2 < $n; $i2++) {
						if ($i2 == $p)
							$A2[$p][$p] = $A1[$p][$p] * $cs * $cs + $A1[$q][$q] * $sn * $sn -
                                          2.0 * $A1[$p][$q] * $sn * $cs;
						else if ($i2 == $q)
							$A2[$p][$q] = 0.0;
						else
							$A2[$p][$i2] = $A1[$p][$i2] * $cs - $A1[$q][$i2] * $sn;
					}
				}
				else if ($i1 == $q) {
					for ($i2 = 0; $i2 < $n; $i2++) {
						if ($i2 == $q)
							$A2[$q][$q] = $A1[$p][$p] * $sn * $sn + $A1[$q][$q] * $cs * $cs +
                                          2.0 * $A1[$p][$q] * $sn * $cs;
						else if ($i2 == $p)
							$A2[$q][$p] = 0.0;
						else
							$A2[$q][$i2] = $A1[$q][$i2] * $cs + $A1[$p][$i2] * $sn;
					}
				}
				else {
					for ($i2 = 0; $i2 < $n; $i2++) {
						if ($i2 == $p)
							$A2[$i1][$p] = $A1[$i1][$p] * $cs - $A1[$i1][$q] * $sn;
						else if ($i2 == $q)
							$A2[$i1][$q] = $A1[$i1][$q] * $cs + $A1[$i1][$p] * $sn;
						else
							$A2[$i1][$i2] = $A1[$i1][$i2];
					}
				}
			}
								// Xkの計算
			for ($i1 = 0; $i1 < $n; $i1++) {
				for ($i2 = 0; $i2 < $n; $i2++) {
					if ($i2 == $p)
						$X2[$i1][$p] = $X1[$i1][$p] * $cs - $X1[$i1][$q] * $sn;
					else if ($i2 == $q)
						$X2[$i1][$q] = $X1[$i1][$q] * $cs + $X1[$i1][$p] * $sn;
					else
						$X2[$i1][$i2] = $X1[$i1][$i2];
				}
			}
								// 次のステップへ
			$k++;
			for ($i1 = 0; $i1 < $n; $i1++) {
				for ($i2 = 0; $i2 < $n; $i2++) {
					$A1[$i1][$i2] = $A2[$i1][$i2];
					$X1[$i1][$i2] = $X2[$i1][$i2];
				}
			}
		}
	}

	return $ind;
}

?>

  5. Ruby 

    #***********************************************/
# 実対称行列の固有値・固有ベクトル（ヤコビ法） */
#      coded by Y.Suganuma                     */
#***********************************************/

#************************************************************/
# 実対称行列の固有値・固有ベクトル（ヤコビ法）              */
#      n : 次数                                             */
#      ct : 最大繰り返し回数                                */
#      eps : 収束判定条件                                   */
#      a : 対象とする行列                                   */
#      a1, a2 : 作業域（nxnの行列），a1の対角要素が固有値   */
#      x1, x2 : 作業域（nxnの行列），x1の各列が固有ベクトル */
#      return : =0 : 正常                                   */
#               =1 : 収束せず                               */
#      coded by Y.Suganuma                                  */
#************************************************************/
def Jacobi(n, ct, eps, a, a1, a2, x1, x2)
					# 初期設定
	k   = 0
	ind = 1
	p   = 0
	q   = 0
	for i1 in 0 ... n
		for i2 in 0 ... n
			a1[i1][i2] = a[i1][i2]
			x1[i1][i2] = 0.0
		end
		x1[i1][i1] = 1.0
	end
					# 計算
	while ind > 0 && k < ct
						# 最大要素の探索
		max = 0.0
		for i1 in 0 ... n
			for i2 in 0 ... n
				if i2 != i1
					if a1[i1][i2].abs() > max
						max = a1[i1][i2].abs()
						p   = i1
						q   = i2
					end
				end
			end
		end
						# 収束判定
							# 収束した
		if max < eps
			ind = 0
							# 収束しない
		else
								# 準備
			s  = -a1[p][q]
			t  = 0.5 * (a1[p][p] - a1[q][q])
			v  = t.abs() / Math.sqrt(s * s + t * t)
			sn = Math.sqrt(0.5 * (1.0 - v))
			if s*t < 0.0
				sn = -sn
			end
			cs = Math.sqrt(1.0 - sn * sn)
								# akの計算
			for i1 in 0 ... n
				if i1 == p
					for i2 in 0 ... n
						if i2 == p
							a2[p][p] = a1[p][p] * cs * cs + a1[q][q] * sn * sn -
                                       2.0 * a1[p][q] * sn * cs
						elsif i2 == q
							a2[p][q] = 0.0
						else
							a2[p][i2] = a1[p][i2] * cs - a1[q][i2] * sn
						end
					end
				elsif i1 == q
					for i2 in 0 ... n
						if (i2 == q)
							a2[q][q] = a1[p][p] * sn * sn + a1[q][q] * cs * cs +
                                       2.0 * a1[p][q] * sn * cs
						elsif i2 == p
							a2[q][p] = 0.0
						else
							a2[q][i2] = a1[q][i2] * cs + a1[p][i2] * sn
						end
					end
				else
					for i2 in 0 ... n
						if i2 == p
							a2[i1][p] = a1[i1][p] * cs - a1[i1][q] * sn
						elsif i2 == q
							a2[i1][q] = a1[i1][q] * cs + a1[i1][p] * sn
						else
							a2[i1][i2] = a1[i1][i2]
						end
					end
				end
			end
								# xkの計算
			for i1 in 0 ... n
				for i2 in 0 ... n
					if i2 == p
						x2[i1][p] = x1[i1][p] * cs - x1[i1][q] * sn
					elsif i2 == q
						x2[i1][q] = x1[i1][q] * cs + x1[i1][p] * sn
					else
						x2[i1][i2] = x1[i1][i2]
					end
				end
			end
								# 次のステップへ
			k += 1
			for i1 in 0 ... n
				for i2 in 0 ... n
					a1[i1][i2] = a2[i1][i2]
					x1[i1][i2] = x2[i1][i2]
				end
			end
		end
	end

	return ind
end

					# データの設定
ct   = 1000
eps  = 1.0e-10
n    = 3
a    = Array.new(n)
a1   = Array.new(n)
a2   = Array.new(n)
x1   = Array.new(n)
x2   = Array.new(n)
for i1 in 0 ... n
	a[i1]  = Array.new(n)
	a1[i1] = Array.new(n)
	a2[i1] = Array.new(n)
	x1[i1] = Array.new(n)
	x2[i1] = Array.new(n)
end

a[0][0] = 1.0
a[0][1] = 0.0
a[0][2] = -1.0

a[1][0] = 0.0
a[1][1] = 1.0
a[1][2] = -1.0

a[2][0] = -1.0
a[2][1] = -1.0
a[2][2] = 2.0
				# 計算
ind = Jacobi(n, ct, eps, a, a1, a2, x1, x2)
				# 出力
if ind > 0
	printf("収束しませんでした！\n")
else
	printf("固有値\n")
	for i1 in 0 ... n
		printf(" %f", a1[i1][i1])
	end
	printf("\n固有ベクトル（各列が固有ベクトル）\n")
	for i1 in 0 ... n
		for i2 in 0 ... n
			printf(" %f", x1[i1][i2])
		end
		printf("\n")
	end
end

  6. Python 

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

############################################
# 実対称行列の固有値・固有ベクトル（ヤコビ法）
#      n : 次数
#      ct : 最大繰り返し回数
#      eps : 収束判定条件
#      A : 対象とする行列
#      A1, A2 : 作業域（nxnの行列），A1の対角要素が固有値
#      X1, X2 : 作業域（nxnの行列），X1の各列が固有ベクトル
#      return : =0 : 正常
#               =1 : 収束せず
#      coded by Y.Suganuma
############################################

def Jacobi(n, ct, eps, A, A1, A2, X1, X2) :
			# 初期設定
	k   = 0
	ind = 1
	p   = 0
	q   = 0
	for i1 in range(0, n) :
		for i2 in range(0, n) :
			A1[i1][i2] = A[i1][i2]
			X1[i1][i2] = 0.0
		X1[i1][i1] = 1.0
			# 計算
	while ind > 0 and k < ct :
				# 最大要素の探索
		max = 0.0
		for i1 in range(0, n) :
			for i2 in range(0, n) :
				if i2 != i1 :
					if abs(A1[i1][i2]) > max :
						max = fabs(A1[i1][i2])
						p   = i1
						q   = i2
				# 収束判定
					# 収束した
		if max < eps :
			ind = 0
					# 収束しない
		else :
						# 準備
			s  = -A1[p][q]
			t  = 0.5 * (A1[p][p] - A1[q][q])
			v  = fabs(t) / sqrt(s * s + t * t)
			sn = sqrt(0.5 * (1.0 - v))
			if s*t < 0.0 :
				sn = -sn
			cs = sqrt(1.0 - sn * sn)
						# Akの計算
			for i1 in range(0, n) :
				if i1 == p :
					for i2 in range(0, n) :
						if i2 == p :
							A2[p][p] = A1[p][p] * cs * cs + A1[q][q] * sn * sn - 2.0 * A1[p][q] * sn * cs
						elif i2 == q :
							A2[p][q] = 0.0
						else :
							A2[p][i2] = A1[p][i2] * cs - A1[q][i2] * sn
				elif i1 == q :
					for i2 in range(0, n) :
						if i2 == q :
							A2[q][q] = A1[p][p] * sn * sn + A1[q][q] * cs * cs + 2.0 * A1[p][q] * sn * cs
						elif i2 == p :
							A2[q][p] = 0.0
						else :
							A2[q][i2] = A1[q][i2] * cs + A1[p][i2] * sn
				else :
					for i2 in range(0, n) :
						if i2 == p :
							A2[i1][p] = A1[i1][p] * cs - A1[i1][q] * sn
						elif i2 == q :
							A2[i1][q] = A1[i1][q] * cs + A1[i1][p] * sn
						else :
							A2[i1][i2] = A1[i1][i2]
						# Xkの計算
			for i1 in range(0, n) :
				for i2 in range(0, n) :
					if i2 == p :
						X2[i1][p] = X1[i1][p] * cs - X1[i1][q] * sn
					elif i2 == q :
						X2[i1][q] = X1[i1][q] * cs + X1[i1][p] * sn
					else :
						X2[i1][i2] = X1[i1][i2]
						# 次のステップへ
			k += 1
			for i1 in range(0, n) :
				for i2 in range(0, n) :
					A1[i1][i2] = A2[i1][i2]
					X1[i1][i2] = X2[i1][i2]

	return ind

############################################
# 実対称行列の固有値・固有ベクトル（ヤコビ法）
#      coded by Y.Suganuma
############################################
			# データの設定
ct  = 1000
eps = 1.0e-10
n   = 3
A   = np.array([[1, 0, -1],[0, 1, -1],[-1, -1, 2]], np.float)
A1  = np.empty((n, n), np.float)
A2  = np.empty((n, n), np.float)
X1  = np.empty((n, n), np.float)
X2  = np.empty((n, n), np.float)
			# 計算
ind = Jacobi(n, ct, eps, A, A1, A2, X1, X2)
			# 出力
if ind > 0 :
	print("収束しませんでした！")
else :
	print("固有値")
	print(" ", A1[0][0], A1[1][1], A1[2][2])
	print("固有ベクトル（各列が固有ベクトル）")
	for i1 in range(0, n) :
		for i2 in range(0, n) :
			print(" " + str(X1[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       = 1000;
		int n        = 3;
		double eps   = 1.0e-10;
		double[][] A = new double [][] {
			new double[] {1.0, 0.0, -1.0},
			new double[] {0.0, 1.0, -1.0},
			new double[] {-1.0, -1.0, 2.0}
		};
		double[][] A1 = new double [n][];
		double[][] A2 = new double [n][];
		double[][] X1 = new double [n][];
		double[][] X2 = new double [n][];
		for (int i1 = 0; i1 < n; i1++) {
			A1[i1] = new double [n];
			A2[i1] = new double [n];
			X1[i1] = new double [n];
			X2[i1] = new double [n];
		}
					// 計算
		int ind = Jacobi(n, ct, eps, A, A1, A2, X1, X2);
					// 出力
		if (ind > 0)
			Console.WriteLine("収束しませんでした！");
		else {
			Console.WriteLine("固有値");
			for (int i1 = 0; i1 < n; i1++)
				Console.Write(" " + A1[i1][i1]);
			Console.WriteLine();
			Console.WriteLine("固有ベクトル（各列が固有ベクトル）");
			for (int i1 = 0; i1 < n; i1++) {
				for (int i2 = 0; i2 < n; i2++)
					Console.Write(" " + X1[i1][i2]);
				Console.WriteLine();
			}
		}
	}

	/*************************************************************/
	/* 実対称行列の固有値・固有ベクトル（ヤコビ法）              */
	/*      n : 次数                                             */
	/*      ct : 最大繰り返し回数                                */
	/*      eps : 収束判定条件                                   */
	/*      A : 対象とする行列                                   */
	/*      A1, A2 : 作業域（nxnの行列），A1の対角要素が固有値   */
	/*      X1, X2 : 作業域（nxnの行列），X1の各列が固有ベクトル */
	/*      return : =0 : 正常                                   */
	/*               =1 : 収束せず                               */
	/*      coded by Y.Suganuma                                  */
	/*************************************************************/
	int Jacobi(int n, int ct, double eps, double[][] A, double[][] A1, double[][] A2,
	           double[][] X1, double[][] X2)
	{
					// 初期設定
		for (int i1 = 0; i1 < n; i1++) {
			for (int i2 = 0; i2 < n; i2++) {
				A1[i1][i2] = A[i1][i2];
				X1[i1][i2] = 0.0;
			}
			X1[i1][i1] = 1.0;
		}
					// 計算
		int k = 0, ind = 1;
		while (ind > 0 && k < ct) {
						// 最大要素の探索
			double max = 0.0;
			int p = 0, q = 0;
			for (int i1 = 0; i1 < n; i1++) {
				for (int i2 = 0; i2 < n; i2++) {
					if (i2 != i1) {
						if (Math.Abs(A1[i1][i2]) > max) {
							max = Math.Abs(A1[i1][i2]);
							p   = i1;
							q   = i2;
						}
					}
				}
			}
						// 収束判定
							// 収束した
			if (max < eps)
				ind = 0;
							// 収束しない
			else {
								// 準備
				double s  = -A1[p][q];
				double t  = 0.5 * (A1[p][p] - A1[q][q]);
				double v  = Math.Abs(t) / Math.Sqrt(s * s + t * t);
				double sn = Math.Sqrt(0.5 * (1.0 - v));
				if (s*t < 0.0)
					sn = -sn;
				double cs = Math.Sqrt(1.0 - sn * sn);
								// Akの計算
				for (int i1 = 0; i1 < n; i1++) {
					if (i1 == p) {
						for (int i2 = 0; i2 < n; i2++) {
							if (i2 == p)
								A2[p][p] = A1[p][p] * cs * cs + A1[q][q] * sn * sn -
                                           2.0 * A1[p][q] * sn * cs;
							else if (i2 == q)
								A2[p][q] = 0.0;
							else
								A2[p][i2] = A1[p][i2] * cs - A1[q][i2] * sn;
						}
					}
					else if (i1 == q) {
						for (int i2 = 0; i2 < n; i2++) {
							if (i2 == q)
								A2[q][q] = A1[p][p] * sn * sn + A1[q][q] * cs * cs +
                                           2.0 * A1[p][q] * sn * cs;
							else if (i2 == p)
								A2[q][p] = 0.0;
							else
								A2[q][i2] = A1[q][i2] * cs + A1[p][i2] * sn;
						}
					}
					else {
						for (int i2 = 0; i2 < n; i2++) {
							if (i2 == p)
								A2[i1][p] = A1[i1][p] * cs - A1[i1][q] * sn;
							else if (i2 == q)
								A2[i1][q] = A1[i1][q] * cs + A1[i1][p] * sn;
							else
								A2[i1][i2] = A1[i1][i2];
						}
					}
				}
								// Xkの計算
				for (int i1 = 0; i1 < n; i1++) {
					for (int i2 = 0; i2 < n; i2++) {
						if (i2 == p)
							X2[i1][p] = X1[i1][p] * cs - X1[i1][q] * sn;
						else if (i2 == q)
							X2[i1][q] = X1[i1][q] * cs + X1[i1][p] * sn;
						else
							X2[i1][i2] = X1[i1][i2];
					}
				}
								// 次のステップへ
				k++;
				for (int i1 = 0; i1 < n; i1++) {
					for (int i2 = 0; i2 < n; i2++) {
						A1[i1][i2] = A2[i1][i2];
						X1[i1][i2] = X2[i1][i2];
					}
				}
			}
		}

		return ind;
	}
}

  8. VB 

    ''''''''''''''''''''''''''''''''''''''''''''''''
' 実対称行列の固有値・固有ベクトル（ヤコビ法） '
'      coded by Y.Suganuma                     '
''''''''''''''''''''''''''''''''''''''''''''''''
Module Test

	Sub Main()
					' データの設定
		Dim ct As Integer  = 1000
		Dim n As Integer   = 3
		Dim eps As Double  = 1.0e-10
		Dim A(,) As Double = {{1.0, 0.0, -1.0}, {0.0, 1.0, -1.0}, {-1.0, -1.0, 2.0}}
		Dim A1(n,n) As Double
		Dim A2(n,n) As Double
		Dim X1(n,n) As Double
		Dim X2(n,n) As Double
					' 計算
		Dim ind As Integer = Jacobi(n, ct, eps, A, A1, A2, X1, X2)
					' 出力
		If ind > 0
			Console.WriteLine("収束しませんでした！")
		Else
			Console.WriteLine("固有値")
			For i1 As Integer = 0 To n-1
				Console.Write(" " & A1(i1,i1))
			Next
			Console.WriteLine()
			Console.WriteLine("固有ベクトル（各列が固有ベクトル）")
			For i1 As Integer = 0 To n-1
				For i2 As Integer = 0 To n-1
					Console.Write(" " & X1(i1,i2))
				Next
				Console.WriteLine()
			Next
		End If

	End Sub

	'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
	' 実対称行列の固有値・固有ベクトル（ヤコビ法）              '
	'      n : 次数                                             '
	'      ct : 最大繰り返し回数                                '
	'      eps : 収束判定条件                                   '
	'      A : 対象とする行列                                   '
	'      A1, A2 : 作業域（nxnの行列），A1の対角要素が固有値   '
	'      X1, X2 : 作業域（nxnの行列），X1の各列が固有ベクトル '
	'      return : =0 : 正常                                   '
	'               =1 : 収束せず                               '
	'      coded by Y.Suganuma                                  '
	'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
	Function Jacobi(n As Integer, ct As Integer, eps As Double, A(,) As Double,
	                A1(,) As Double, A2(,) As Double, X1(,) As Double, X2(,) As Double)
					' 初期設定
		For i1 As Integer = 0 To n-1
			For i2 As Integer = 0 To n-1
				A1(i1,i2) = A(i1,i2)
				X1(i1,i2) = 0.0
			Next
			X1(i1,i1) = 1.0
		Next
					' 計算
		Dim k As Integer   = 0
		Dim ind As Integer = 1
		Do While ind > 0 and k < ct
						' 最大要素の探索
			Dim max As Double = 0.0
			Dim p As Integer  = 0
			Dim q As Integer  = 0
			For i1 As Integer = 0 To n-1
				For i2 As Integer = 0 To n-1
					If i2 <> i1
						If Math.Abs(A1(i1,i2)) > max
							max = Math.Abs(A1(i1,i2))
							p   = i1
							q   = i2
						End If
					End If
				Next
			Next
						' 収束判定
							' 収束した
			If max < eps
				ind = 0
							' 収束しない
			Else
								' 準備
				Dim s As Double  = -A1(p,q)
				Dim t As Double  = 0.5 * (A1(p,p) - A1(q,q))
				Dim v As Double  = Math.Abs(t) / Math.Sqrt(s * s + t * t)
				Dim sn As Double = Math.Sqrt(0.5 * (1.0 - v))
				If s*t < 0.0
					sn = -sn
				End If
				Dim cs As Double = Math.Sqrt(1.0 - sn * sn)
								' Akの計算
				For i1 As Integer = 0 To n-1
					If i1 = p
						For i2 As Integer = 0 To n-1
							If i2 = p
								A2(p,p) = A1(p,p) * cs * cs + A1(q,q) * sn * sn -
                                           2.0 * A1(p,q) * sn * cs
							ElseIf i2 = q
								A2(p,q) = 0.0
							Else
								A2(p,i2) = A1(p,i2) * cs - A1(q,i2) * sn
							End If
						Next
					ElseIf i1 = q
						For i2 As Integer = 0 To n-1
							If i2 = q
								A2(q,q) = A1(p,p) * sn * sn + A1(q,q) * cs * cs +
                                           2.0 * A1(p,q) * sn * cs
							ElseIf i2 = p
								A2(q,p) = 0.0
							Else
								A2(q,i2) = A1(q,i2) * cs + A1(p,i2) * sn
							End If
						Next
					Else
						For i2 As Integer = 0 To n-1
							If i2 = p
								A2(i1,p) = A1(i1,p) * cs - A1(i1,q) * sn
							ElseIf i2 = q
								A2(i1,q) = A1(i1,q) * cs + A1(i1,p) * sn
							Else
								A2(i1,i2) = A1(i1,i2)
							End If
						Next
					End If
				Next
								' Xkの計算
				For i1 As Integer = 0 To n-1
					For i2 As Integer = 0 To n-1
						If i2 = p
							X2(i1,p) = X1(i1,p) * cs - X1(i1,q) * sn
						ElseIf i2 = q
							X2(i1,q) = X1(i1,q) * cs + X1(i1,p) * sn
						Else
							X2(i1,i2) = X1(i1,i2)
						End If
					Next
				Next
								' 次のステップへ
				k += 1
				For i1 As Integer = 0 To n-1
					For i2 As Integer = 0 To n-1
						A1(i1,i2) = A2(i1,i2)
						X1(i1,i2) = X2(i1,i2)
					Next
				Next
			End If
		Loop

		Return ind

	End Function

End Module

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