静岡理工科大学 菅沼ホーム 目次 索引

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

    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;
    	}
    }
    			
    アプレット版
    /************************************************/
    /* 実対称行列の固有値・固有ベクトル(ヤコビ法) */
    /*      coded by Y.Suganuma                     */
    /************************************************/
    import java.awt.*;
    import java.awt.event.*;
    import java.applet.*;
    import java.util.StringTokenizer;
    
    public class Jacobi extends Applet implements ActionListener {
    	TextField order, max;
    	TextArea ta1, ta2;
    	Button bt;
    	Font f = new Font("TimesRoman", Font.BOLD, 20);
    
    	/************/
    	/* 初期設定 */
    	/************/
    	public void init()
    	{
    					// レイアウト,背景色,フォント
    		setLayout(new BorderLayout(5, 5));
    		setBackground(new Color(225, 255, 225));
    		setFont(f);
    					// 上のパネル
    		Panel pn1 = new Panel();
    		add(pn1, BorderLayout.NORTH);
    
    		pn1.add(new Label("次数:"));
    		order = new TextField("3", 3);
    		pn1.add(order);
    
    		pn1.add(new Label(" 最大繰り返し回数:"));
    		max = new TextField("1000", 5);
    		pn1.add(max);
    
    		pn1.add(new Label("  "));
    		bt = new Button("OK");
    		bt.setBackground(Color.pink);
    		bt.addActionListener(this);
    		pn1.add(bt);
    					// 中央のパネル(問題)
    		Panel pn2 = new Panel();
    		add(pn2, BorderLayout.CENTER);
    
    		ta1 = new TextArea("1 0 -1\n0 1 -1\n-1 -1 2", 10, 30);
    		pn2.add(ta1);
    					// 下のパネル(答え)
    		Panel pn3 = new Panel();
    		add(pn3, BorderLayout.SOUTH);
    
    		ta2 = new TextArea(10, 40);
    		pn3.add(ta2);
    	}
    
    	/******************************/
    	/* ボタンが押されたときの処理 */
    	/******************************/
    	public void actionPerformed(ActionEvent e)
    	{
    		if (e.getSource() == bt) {
    			double A[][], A1[][], A2[][], X1[][], X2[][], eps;
    			int ind, ct, n;
    					// 問題の設定
    			eps = 1.0e-10;
    			n   = Integer.parseInt(order.getText());
    			ct  = Integer.parseInt(max.getText());
    			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];
    
    			String s = ta1.getText();
    			StringTokenizer str = new StringTokenizer(s, " \n");
    			int k = 0;
    			while (str.hasMoreTokens()) {
    				A[k/n][k%n] = Double.parseDouble(str.nextToken());
    				k++;
    			}
    					// 計算
    			ind = jacobi(n, ct, eps, A, A1, A2, X1, X2);
    					// 結果
    			if (ind > 0) {
    				ta2.setForeground(Color.red);
    				ta2.setText(" 解を求めることができません");
    			}
    			else {
    				ta2.setForeground(Color.black);
    				ta2.setText("固有値:");
    				for (int i1 = 0; i1 < n; i1++)
    					ta2.append(" " + String.format("%.5f",A1[i1][i1]));
    				ta2.append("\n固有ベクトル:\n");
    				for (int i1 = 0; i1 < n; i1++) {
    					ta2.append("   ");
    					for (int i2 = 0; i2 < n; i2++)
    						ta2.append(" " + String.format("%.5f",X1[i1][i2]));
    					ta2.append("\n");
    				}
    			}
    		}
    	}
    
    	/*************************************************************/
    	/* 実対称行列の固有値・固有ベクトル(ヤコビ法)              */
    	/*      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[][])
    	{
    		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
    			

静岡理工科大学 菅沼ホーム 目次 索引