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

最大固有値と固有ベクトル(べき乗法)

    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;
    	}
    }
    			
    アプレット版
    /****************************************/
    /* 最大固有値と固有ベクトル(べき乗法) */
    /*      coded by Y.Suganuma             */
    /****************************************/
    import java.awt.*;
    import java.awt.event.*;
    import java.applet.*;
    import java.util.StringTokenizer;
    
    public class Power extends Applet implements ActionListener {
    	TextField order, mod, max, ini;
    	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("修正:"));
    		mod = new TextField("15", 3);
    		pn1.add(mod);
    
    		pn1.add(new Label("最大繰返し:"));
    		max = new TextField("200", 5);
    		pn1.add(max);
    
    		pn1.add(new Label("初期値:"));
    		ini = new TextField("1 1 1", 20);
    		pn1.add(ini);
    
    		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("11 6 -2\n-2 18 1\n-12 24 13", 10, 30);
    		pn2.add(ta1);
    					// 下のパネル(答え)
    		Panel pn3 = new Panel();
    		add(pn3, BorderLayout.SOUTH);
    
    		ta2 = new TextArea(10, 70);
    		pn3.add(ta2);
    	}
    
    	/******************************/
    	/* ボタンが押されたときの処理 */
    	/******************************/
    	public void actionPerformed(ActionEvent e)
    	{
    		if (e.getSource() == bt) {
    			double A[][], B[][], C[][], r[], v[][], v0[], v1[], v2[], eps;
    			int ind, ct, m, n;
    					// 問題の設定
    			eps = 1.0e-10;
    			n   = Integer.parseInt(order.getText());
    			m   = Integer.parseInt(mod.getText());
    			ct  = Integer.parseInt(max.getText());
    			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];
    
    			String s = ini.getText();
    			StringTokenizer str = new StringTokenizer(s, " \n");
    			int k = 0;
    			while (str.hasMoreTokens()) {
    				v0[k] = Double.parseDouble(str.nextToken());
    				k++;
    			}
    
    			s   = ta1.getText();
    			str = new StringTokenizer(s, " \n");
    			k   = 0;
    			while (str.hasMoreTokens()) {
    				A[k/n][k%n] = Double.parseDouble(str.nextToken());
    				k++;
    			}
    
    			for (int i1 = 0; i1 < n; i1++) {
    				for (int i2 = 0; i2 < n; i2++)
    					B[i1][i2] = 0.0;
    				B[i1][i1] = 1.0;
    			}
    					// 計算
    			ind = power(0, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
    					// 結果
    			if (ind == 0) {
    				ta2.setForeground(Color.red);
    				ta2.setText(" 解を求めることができません");
    			}
    			else {
    				ta2.setForeground(Color.black);
    				ta2.setText("");
    				for (int i1 = 0; i1 < ind; i1++) {
    					ta2.append("固有値: " + String.format("%.5f",r[i1]) + "\n");
    					ta2.append("  固有ベクトル:");
    					for (int i2 = 0; i2 < n; i2++)
    						ta2.append(" " + String.format("%.5f",v[i1][i2]));
    					ta2.append("\n");
    				}
    			}
    		}
    	}
    
    	/****************************************/
    	/* 最大固有値と固有ベクトル(べき乗法) */
    	/*      i : 何番目の固有値かを示す      */
    	/*      n : 次数                        */
    	/*      m : 丸め誤差の修正回数          */
    	/*      ct : 最大繰り返し回数           */
    	/*      eps : 収束判定条件              */
    	/*      A : 対象とする行列              */
    	/*      B : nxnの行列,最初は,単位行列 */
    	/*      C : 作業域,nxnの行列           */
    	/*      r : 固有値                      */
    	/*      v : 各行が固有ベクトル(nxn)   */
    	/*      v0 : 固有ベクトルの初期値       */
    	/*      v1,v2 : 作業域(n次元ベクトル) */
    	/*      return : 求まった固有値の数     */
    	/*      coded by Y.Suganuma             */
    	/****************************************/
    	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
    			

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