本文所採用的是冪演算法計算矩陣最大特徵值（只計算最大特徵值），具體的細節可以參考：http://blog.csdn.net/whucv/article/details/7636135

這兒只給出虛擬碼和原始碼，方便各位同學使用。

假設待計算的矩陣為A[N][N]，一個N維的方陣.

申明一個輔助向量（一維陣列）v[N].並初始化

//對應函式InitialV()

for i=1:n

v[i]=1

end

然後進行不停的迭代，當迭代的次數趨於無窮時，max的值越接近於真實的最大特徵值。

在實際中，可以考察前後兩次迭代的值的差距，如果該差距在可以接受的範圍內（比如0.00001），則可以結束迭代

迭代過程：

while(true)

1.v=matrix * v; // 對應函式ComputeNewV()

2.max=v(1:n); // 對應函式GetMaxV();

3.v=v/max; // 對應函式NormalizeV();

end

解釋：

1.矩陣相乘。用matrix和v向量進行相乘，結果儲存在v中。相乘的過程和矩陣乘法一樣（需要假象將v倒置）

2.獲取上一步後產生的最大v值，記為max（即一維陣列v中的最大值），用該值去除原向量v，保證v陣列的元素值不超過1.（注，在除法前所獲得的max值才是我們的特徵估計值）

迴圈迭代上述過程，直到前後兩次的估計值（max）在我們的誤差範圍內，結束迭代。

// Matrix.cpp : 定義控制檯應用程式的入口點。

//

#include "stdafx.h"

#include<math.h>

using namespace std;

#define N 3 //行列 NxN方陣

#define e 0.0000001 //誤差限

double matrix[N][N];

double vector[N];

void ComputeNewV();

double GetMaxV();

void NormalizeV();

void InitialV();

int main()

{

///初始化待求矩陣

matrix[0][0]=1.0;

matrix[0][1]=1.0;

matrix[0][2]=0.5;

matrix[1][0]=1.0;

matrix[1][1]=1.0;

matrix[1][2]=0.25;

matrix[2][0]=0.5;

matrix[2][1]=0.25;

matrix[2][2]=2.0;

InitialV();

///分別記錄上一次和這一次的迭代值，即我們需要的最大特徵值

///當誤差連續K次小於一個範圍時，便結束迭代

///關於這一點我也不太確定，沒仔細閱讀教材

///一個可行的辦法是根據可以接受的時間直接迭代10000次甚至更高

double lastValue;

double currentValue=0.0;

int k=0;

for(int i=0;i<10000;i++)

{

ComputeNewV();

lastValue=currentValue;

currentValue=GetMaxV();

if(abs(currentValue-lastValue)<e)

k++;

else

k=0;

if(k>10)

break;

NormalizeV();

}

cout<<"wanted:"<<currentValue<<endl;

return 0;

}

///初始化向量，可以賦予非0任何值

void InitialV()

{

for(int i=0;i<N;i++)

vector[i]=1.0;

}

///迭代計算vector=matrix*vector

void ComputeNewV()

{

double temp[N];//申明臨時變數儲存結果

for(int i=0; i<N; i++)

{

temp[i]=0.0;

for(int j=0;j<N;j++)

temp[i]+=matrix[i][j]*vector[j];

}

///複製結果到vector向量

for(int i=0;i<N;i++)

vector[i]=temp[i];

}

///獲取向量vector的最大值,也是這次迭代之後的一個最大特徵估計值

double GetMaxV()

{

double result=vector[0];

for(int i=1;i<N;i++)

if(vector[i]>result)

result=vector[i];

return result;

}

///用當前vector的最大值對向量做除法，歸一化

void NormalizeV()

{

double max=GetMaxV();

for(int i=0;i<N;i++)

vector[i]/=max;

}