[数值算法]求矩阵的最大特征值的幂法

[数值算法]求矩阵的最大特征值的幂法.

对于工程计算而言，矩阵的特征值和特征向量都是相当重要和常见的数据，这里给出的幂法是一种常见的求解方法，用的是迭代的思想。

符号说明:

1A为待求的矩阵，

2Uk，Vk为迭代用的列向量。

3最后的最大特征值maxLamda由最后一次的max(Uk)-----求Uk中的绝对值最大的元素的绝对值.所决定。

而maxLamda所对应的特征向量由最后一次迭代的Vk所决定.

 

 

主要的想法就是先选一个不为0的初始向量U0!=0，然后按下面的式子迭代。

U0=V0!=0

Do

{

       Uk=AVk-1

       Vk=Uk/max( Uk)

}while(abs(max(Uk)-max(Uk-1))>=e)//e为精度.

好了，就这样,更多的细节请去参考相关的数值算法书籍.

 

 

在贴出程序之前，先对一部分我新增加的实用函数进行说明：

如:

void twoDArrMemApply(Type*** inArr,int rowNum,int colNum)

{

       int i=0;/*iterator vaule*/

       (*inArr)=(Type**)malloc(sizeof(Type*)*rowNum);

       for(i=0;i<rowNum;i++)

              (*inArr)[i]=(Type*)malloc(sizeof(Type)*colNum);

       assertF(*inArr!=NULL,"in twoDArrMemApply,inArr at last is null/n");

}

 

 

void twoDArrMemFree(Type*** inArr,int rowNum)

{

       int i=0;/*iterator value*/

       assertF((*inArr)!=NULL,"in 2d arr mem free,in arr is null/n");

      

       for(i=0;i<rowNum;i++)

              free((*inArr)[i]);

      

       free((*inArr));

}

这两个函数的作用相信大家一看就明白，是实现二维指针的申请内存和释放内存的，这样，以后再主程序里的工作量就会小多了。

 

 

还有，我在写一些程序段的时候对待外部传进的指针采用如下处理手段(纯C条件下)

除非这个函数有特殊的作用，如申请内存，或要读入外部文本内容到二维指针等。 其余的情况，一律不对外部指针进行任何申请或释放内存的处理。

对于要保护数据的外部传入指针，则在函数内部再做一个局部指针，在函数结尾释放.

对局部指针的操作，也仅限于赋值，而绝对不要用外部传入针指去指向它（即赋一个临时区的地址给外部的指针变量），这当然是错误的。

 

 

好了，下面是程序段:

/*for max lamda resolve*/

Type powerMethodForLamda(Type** matrixA,int size,char* outputFileName)

{

       Type maxLamda;

       Type* listV;

       Type* listU;

       FILE* outputFile;/*the outputFile for the data output*/

       Type preMax;/*a tween data*/

       float e=(float)0.0001;/*the precise controller*/

       Type tmpData;/*temp data for program*/

       int i=0;/*iterator times*/

       int iteratorNum=0;/*iterator number*/

      

       /*assertion*/

       assertF(matrixA!=NULL,"in powerMethodFor lamda,matrixA is null/n");

       assertF(outputFileName!=NULL,"in readList,listFileName is null/n");

      

       /*open file*/

       assertF((outputFile=fopen(outputFileName,"wb"))!=NULL,"output file open error/n");

      

       /*mem apply*/

       listArrMemApply(&listV,size);

       listArrMemApply(&listU,size);

      

       /*initialization*/

       for(i=0;i<size;i++)

       {

              listV[i]=0;

              listU[i]=0;

       }

       listV[size-1]=1;

       listU[size-1]=1;

      

       /*core program*/

       fprintf(outputFile,"iteratorTime   maxUk/r/n");

       do

       {

assertF(listNotZero(listU,size),"in the core of  powerMethodForLamda list U is NULL/n");

assertF(listNotZero(listV,size),"in the core of  powerMethodForLamda list V is NULL/n");

              preMax=maxAbsValInList(listU,size);

              matirxBy2DWith1DCol(matrixA,listV,listU,size,size);

              tmpData=1/maxAbsValInList(listU,size);

              numByList(tmpData,listU,listV,size);

             

              fprintf(outputFile,"%-16d% -16f/r/n",iteratorNum,maxAbsValInList(listU,size));

       }

       while(fabs(preMax-maxAbsValInList(listU,size))>=e);

       fprintf(outputFile,"charactstic vector is:/r/n");     

       outputListArrFloat(listV,0,size,outputFile);

       /*End of the Core Program*/

       maxLamda=maxAbsValInList(listU,size);

       fprintf(outputFile,"the max lamda is:/r/n  %f./r/n",maxLamda);

       /*mem free*/

       free(listV);

       free(listU);

      

       /*close the file*/

       fclose(outputFile);

      

       return maxLamda;

}

 

 

/*相应的辅助函数*/

/* matirxBy2DWith1DCol 一个n*n的矩阵和一个n*1的列向量作乘法*/

void matirxBy2DWith1DCol(Type** matrixA,Type* matrixListIn,Type* matrixListAns,int rowNum,int mNum)

{

       /*variable declare*/

       int i,k;/*iterator number*/

       Type sum;

      

       /*assertion*/

       assertF(matrixA!=NULL,"in twoMatrixBy matrixA is null/n");

       assertF(matrixListIn!=NULL,"in twoMatrixBy matrixB is null/n");

       assertF(matrixListAns!=NULL,"in twoMatrixBy matrixAns is null/n");

      

       /*core program*/

       for(i=0;i<rowNum;i++)

       {

                     sum=0;

                     for(k=0;k<mNum;k++)

                            sum+=matrixA[i][k]*matrixListIn[k];

                     matrixListAns[i]=sum;

       }

}

/*求一个一维向量中绝对值的最大值*/

Type maxAbsValInList(Type* inList,int len)

{

       int i;/*iterator num*/

       Type maxData;

       assertF(inList!=NULL,"in maxValInList,inList is NULL/n");

       maxData=(Type)fabs(inList[0]);

       for(i=1;i<len;i++)

              if(fabs(inList[i])>maxData)maxData=(Type)fabs(inList[i]);

       return maxData;

}

 

 

/*test program*/

/*maxLamda resolve  test program*/

#include "Global.h"

#include "Ulti.h"

#include "Matrix.h"

#include "MyAssert.h"

#include <time.h>

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

 

 

 

 

char *inFileName="inputData.txt";

/*

       input data specification

       row,col;

       //Arr

              a11,a12,...;

              , , , ,     

              , , , ,

              an1,an2,...;

       }

*/

 

 

char *outFileName="outputData.txt";

#define DEBUG 1

 

 

void main(int argc,char* argv[])

{

       FILE *inputFile;/*input file*/

       FILE *outputFile;/*output file*/

 

 

       double startTime,endTime,tweenTime;/*time callopsed info*/

      

       int rowNum,colNum;

       Type** wArr;

       Type maxLamda;

 

 

       int n;/*arr deminision for squre matrix*/

 

 

       /*default input file open*/

       if(argc>1)strcpy(inFileName,argv[1]);

       assertF((inputFile=fopen(inFileName,"rb"))!=NULL,"input file error");

       printf("input file open success/n");

      

       /*default outpout file open*/

       if(argc>2)strcpy(outFileName,argv[2]);

       assertF((outputFile=fopen(outFileName,"wb"))!=NULL,"output file error");

       printf("output file open success/n");

      

       /*This function,automatically fullfill the task of

       apply the mem for the 2d pointers. Perfect,right? :)*/

       read2DArrFloat(&wArr,&rowNum,&colNum,"inputData2.txt");

 

 

#if  DEBUG

       printf("/n*******start of test program******/n");

       printf("now is runnig,please wait.../n");

       startTime=(double)clock()/(double)CLOCKS_PER_SEC;

       /******************Core program code*************/

              /*argu prepare*/

              assertF(colNum==rowNum,"in test colNum!=rowNum");

              n=rowNum;/*get the size of square matrix*/

              maxLamda=powerMethodForLamda(wArr,n,"output2.txt");

              printf("the max lamda is:%f./r/n",maxLamda);

              fprintf(outputFile,"the max lamda is:%f./r/n",maxLamda);

       /******************End of Core program**********/

       endTime=(double)clock()/(double)CLOCKS_PER_SEC;

       tweenTime=endTime-startTime;/*Get the time collapsed*/

       /*Time collapsed output*/

       printf("the collapsed time in this algorithm implement is:%f/n",tweenTime);

       fprintf(outputFile,"the collapsed time in this algorithm implement is:%f/r/n",tweenTime); 

       printf("/n*******end of test program******/n");

#endif

       twoDArrMemFree(&wArr,rowNum);

       printf("program end successfully,/n you have to preess any key to clean the buffer area to output,otherwise,you wiil not get the total answer./n");

       getchar();/*Screen Delay Control*/

       return;

}

测试结果:

输入：

3,3;

2,3,2;

10,3,4;

3,6,1;

 

 

输出:

iteratorTime   maxUk

0               4.000000       

0               9.000000       

0               11.444445      

0               10.922329      

0               11.014222      

0               10.997417       

0               11.000469      

0               10.999914      

0               11.000015      

0               10.999997      

charactstic vector is:

0.500000,1.000000,0.750000;

the max lamda is:

  10.999997.

 

 

PS:

最近在数模班里，也的确看到了别的专业一些对程序和算法在理解上比较强能力的人，虽然他们可以告诉你解决了某个问题，但你看一眼他写的程序，充斥在你眼里的，都是大量的“蛮力”算法，像10个层甚至更高层的for循环比比皆是，基本上就是针对这个问题的特型算法。当他们把这个问题这样“解决”的时候，他们会认为事情已经OK了，而我却还不得不承认自己没想出一个比较好的通用算法。呵呵~~~随他去呢，反正我水平本来就不高，但作为一个计算机专业的学生而言，保证算法的通用性和函数的模块性，是首要的任务，如果这个程序里要用什么10个层的循环来做，我还不如不写。还有一些自认为数学感觉好的人，认为经典算法的源程序都在那儿了，自己只要能将它在程序中实现，就了事了.算法在他们的视线里，更多的是只算一种公式，甚至可以不必去了解它的细节.这样的观点和想法，在我这个计算机专业的人看来，同样是不可取的。

路漫漫其修远兮,吾将上下而求索.