数值分析之幂法及反幂法C语言程序实例
- 格式:doc
- 大小:388.00 KB
- 文档页数:10
矩阵的特征值与特征向量的计算摘要物理,力学,工程技术中的很多问题在数学上都归结于求矩阵特征值的问题,例如振动问题(桥梁的振动,机械的振动,电磁振动等)、物理学中某些临界值的确定问题以及理论物理中的一些问题。
矩阵特征值的计算在矩阵计算中是一个很重要的部分,本文使用幂法和反幂法分别求矩阵的按模最大,按模最小特征向量及对应的特征值。
幂法是一种计算矩阵主特征值的一种迭代法,它最大的优点是方法简单,对于稀疏矩阵比较合适,但有时收敛速度很慢。
其基本思想是任取一个非零的初始向量。
由所求矩阵构造一向量序列。
再通过所构造的向量序列求出特征值和特征向量。
反幂法用来计算矩阵按模最小特征向量及其特征值,及计算对应于一个给定近似特征值的特征向量。
本文中主要使用反幂法计算一个矩阵的按模最小特征向量及其对应的特征值。
计算矩阵按模最小特征向量的基本思想是将其转化为求逆矩阵的按模最大特征向量。
然后通过这个按模最大的特征向量反推出原矩阵的按模最小特征向量。
关键词:矩阵;特征值;特征向量;冥法;反冥法THE CALCULATIONS OF EIGENVALUE AND EIGENVECTOR OF MATRIXABSTRACTPhysics, mechanics, engineering technology in a lot of problems in mathematics are attributed to matrix eigenvalue problem, such as vibration (vibration of the bridge, mechanical vibration, electromagnetic vibration, etc.) in physics, some critical values determine problems and theoretical physics in some of the problems. Matrix eigenvalue calculation is a very important part in matrix computation. In this paper, we use the power method and inverse power method to calculate the maximum of the matrix, according to the minimum characteristic vector and the corresponding characteristic value.Power method is an iterative method to calculate the eigenvalues of a matrix. It has the advantage that the method is simple and suitable for sparse matrices, but sometimes the convergence rate is very slow. The basic idea is to take a non - zero initial vector. Construct a vector sequence from the matrix of the matrix. Then the eigenvalues and eigenvectors are obtained by using the constructed vector sequence.The inverse power method is used to calculate the minimum feature vectors and their eigenvalues of the matrix, and to calculate the eigenvalues of the matrix. In this paper, we use the inverse power method to calculate the minimum eigenvalue of a matrix and its corresponding eigenvalues. The basic idea of calculating the minimum characteristic vector of a matrix is to transform it to the maximum characteristic vector of the modulus of the inverse matrix. Then, according to the model, the minimum feature vector of the original matrix is introduced.Key words: Matrix;Eigenvalue;Eigenvector;Iteration methods;目录1 引言 (1)2 相关定理。
1.实验目的:1熟练掌握C 语言程序设计,编程求解问题。
2.运用幂法求解住特征值和特征向量。
2.实验内容:例题:用幂法求 A=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡0.225.05.025.00.10.15.00.10.1 的特征值和特征向量。
完整代码以及截图如下:#include "stdio.h"#include "math.h"#define M 3void main(){float fan(),max(),e1,e2,r1,r2;void au(),ex(),print_x(),std();static float a[M][M]={{1.0,1.0,0.5},{1.0,1.0,0.25},{0.5,0.25,2.0}}; static float u0[M],u1[M],maxn0,maxn1;int i;printf("*********************************\n");printf("****** 幂法*********\n");printf("******求特征值与特征向量*********\n");printf("*********************************\n\n");printf("input precision e1,e2:");scanf("%f,%f",&e1,&e2);printf("\ninput u(%d):",M);for (i=0;i<M;++i){scanf("%f",&u0[i]);}std(u0);maxn0=max(u0);i=0;printf("\n- - - - - - - - - - - - - - - - - -\n");printf(" ............NMD\n");do{au(a,u0,u1);maxn1=max(u1);std(u1);r1=fan(u0,u1);r2=(float)fabs(maxn0-maxn1);maxn0=maxn1;if (r1>e1 || r2>e2){printf("%4d",i++);print_x(u0);printf("\n");ex(u0,u1);}elsebreak;} while (1);}void au(a,u0,u1)float a[][M],u0[],u1[];{int i,j;for (i=0;i<M;i++){u1[i]=0;for (j=0;j<M;j++){u1[i]+=a[i][j]*u0[j];}}}void std(u)float u[];{int i;float t,max();t=max(u);for (i=0;i<M;i++){u[i]=u[i]/t;}}float fan(u0,u1)float u0[],u1[];{float max();int i;float uu[M];for (i=0;i<M;i++){uu[i]=u0[i]-u1[i];}return max(uu);}float max(u)float u[];{int i;float m;m=u[0];for (i=0;i<M;i++){if (u[i]>m){m=u[i];}}return m;}void ex(u0,u1)float u0[],u1[];{int i;for (i=0;i<M;i++){u0[i]=u1[i];}}void print_x(u)float u[];{int i;for (i=0;i<M;i++){printf("%12.6f",u[i]);}}3.运行结果:。
数值计算解矩阵的按模最大最小特征值及对应的特征向量一.幂法1. 幂法简介:当矩阵A 满足一定条件时,在工程中可用幂法计算其主特征值(按模最大)及其特征向量。
矩阵A 需要满足的条件为:(1) 的特征值为A i n λλλλ,0||...||||21≥≥≥>(2) 存在n 个线性无关的特征向量,设为n x x x ,...,,21 1.1计算过程:i ni i i u x x αα,1)0()0(∑==,有对任意向量不全为0,则有111111*********1111011)()(...u u a u a u λu λαu αA x A Ax x k n n k n k k ni i k i i n i i i k )(k (k))(k αλλλλλα++++=+=+++≈⎥⎦⎤⎢⎣⎡+++======∑∑ 可见,当||12λλ越小时,收敛越快;且当k 充分大时,有1)1111)11111λαλαλ=⇒⎪⎩⎪⎨⎧==+++(k )(k k (k k )(k x x u x u x ,对应的特征向量即是)(k x 1+。
2 算法实现.,, 3,,1 , ).5()5(,,,,||).4();max(,).3()(max(;0,1).2(,).1()()()(停机否则输出失败信息转置若转否则输出若计算最大迭代次数,误差限,初始向量输入矩阵βλβεβλβλε←+←<<-←←=←←k k N k y x Ay x x abs x y k Nx A k k k 3 matlab 程序代码function [t,y]=lpowerA,x0,eps,N) % t 为所求特征值,y 是对应特征向量k=1;z=0; % z 相当于λy=x0./max(abs(x0)); % 规范化初始向量x=A*y; % 迭代格式b=max(x); % b 相当于 βif abs(z-b)<eps % 判断第一次迭代后是否满足要求t=max(x);return;endwhile abs(z-b)>eps && k<Nk=k+1;z=b;y=x./max(abs(x));x=A*y;b=max(x);end[m,index]=max(abs(x)); % 这两步保证取出来的按模最大特征值t=x(index); % 是原值,而非其绝对值。
数值计算解矩阵的按模最大最小特征值及对应的特征向量一.幂法1. 幂法简介:当矩阵A 满足一定条件时,在工程中可用幂法计算其主特征值(按模最大)及其特征向量。
矩阵A 需要满足的条件为: (1) 的特征值为A i n λλλλ,0||...||||21≥≥≥>(2) 存在n 个线性无关的特征向量,设为n x x x ,...,,21 1.1计算过程:i ni i i u xx αα,1)0()0(∑==,有对任意向量不全为0,则有1111112211211111111011)()(...u u a u a u λu λαu αA x A Ax x k n n k n k k ni ik i i ni i i k )(k (k))(k αλλλλλα++++=+=+++≈⎥⎦⎤⎢⎣⎡+++======∑∑ 可见,当||12λλ越小时,收敛越快;且当k 充分大时,有1)1111)11111λαλαλ=⇒⎪⎩⎪⎨⎧==+++(k )(k k(k k )(k x x u x u x ,对应的特征向量即是)(k x 1+。
2 算法实现.,, 3,,1 , ).5()5(,,,,||).4();max(,).3()(max(;0,1).2(,).1()()()(停机否则输出失败信息转置若转否则输出若计算最大迭代次数,误差限,初始向量输入矩阵βλβεβλβλε←+←<<-←←=←←k k N k y x Ay x x abs x y k N x A k k k3 matlab 程序代码function [t,y]=lpowerA,x0,eps,N) % t 为所求特征值,y是对应特征向量k=1;z=0; % z 相当于λy=x0./max(abs(x0)); % 规范化初始向量x=A*y; % 迭代格式b=max(x); % b 相当于βif abs(z-b)<eps % 判断第一次迭代后是否满足要求t=max(x);return;endwhile abs(z-b)>eps && k<Nk=k+1;z=b;y=x./max(abs(x));x=A*y;b=max(x);end[m,index]=max(abs(x)); % 这两步保证取出来的按模最大特征值t=x(index); % 是原值,而非其绝对值。
数值分析之幂法及反幂法C 语言程序实例1、算法设计方案:①求1λ、501λ和s λ的值:s λ:s λ表示矩阵的按模最小特征值,为求得s λ直接对待求矩阵A 应用反幂法即可。
1λ、501λ:已知矩阵A 的特征值满足关系 1n λλ<<,要求1λ、及501λ时,可按如下方法求解: a . 对矩阵A 用幂法,求得按模最大的特征值1m λ。
b . 按平移量1m λ对矩阵A 进行原点平移得矩阵1m BA I λ=+,对矩阵B 用反幂法求得B 的按模最小特征值2m λ。
c . 321m m m λλλ=-则:113min(,)m m λλλ=,13max(,)n m m λλλ=即为所求。
②求和A 的与数5011140k k λλμλ-=+最接近的特征值ik λ(k=0,1,…39):求矩阵A 的特征值中与k μ最接近的特征值的大小,采用原点平移的方法:先求矩阵 B=A-k μI 对应的按模最小特征值k β,则k β+k μ即为矩阵A 与k μ最接近的特征值。
重复以上过程39次即可求得ik λ(k=0,1,…39)的值。
③求A 的(谱范数)条件数2cond()A 和行列式det A :在(1)中用反幂法求矩阵A 的按模最小特征值时,要用到Doolittle 分解方法,在Doolittle 分解完成后得到的两个矩阵分别为L 和U ,则A 的行列式可由U 阵求出,即:det(A)=det(U)。
求得det(A)不为0,因此A 为非奇异的实对称矩阵,则: max 2()scond A λλ=,max λ和s λ分别为模最大特征值与模最小特征值。
2、程序源代码:#include<stdio.h>#include<stdio.h>#include<math.h>#define N 501 //列#define M 5 //行#define R 2 //下带宽#define S 2 //上带宽#define K 39#define e 1.0e-12 //误差限float A[M][N]; //初始矩阵float u[N]; //初始向量float y[N],yy[N];float maximum,value1,value2,value_1,value_N,value_s,value_abs_max;const float b=0.16f,c=-0.064f;int max_sign,max_position;void Init_matrix_A() //初始化矩阵A{int i;for(i=2;i<N;i++){A[0][i]= c;}for(i=1;i<N;i++){A[1][i]= b;}for(i=0;i<N;i++){A[2][i]= (1.64-0.024*(i+1))*sin(0.2*(i+1))-0.64*exp(0.1/(i+1));}for(i=0;i<N-1;i++){A[3][i]= b;}for(i=0;i<N-2;i++){A[4][i]= c;}}void Init_u() //初始化迭代向量{int i;for(i=0;i<N;i++)u[i]=1.0;}void Get_max() //获得绝对值最大的数值的模{int i;max_position=0;maximum=fabs(u[0]);for(i=1;i<N;i++){if(maximum<fabs(u[i])){max_position=i;maximum=fabs(u[i]);}}if(u[max_position]<0)max_sign=-1;else max_sign=1;}void Get_y() //单位化迭代向量{int i;for(i=0;i<N;i++)y[i]=u[i]/maximum;}void Get_u() //获得新迭代向量{int i;u[0]=A[2][0]*y[0]+A[1][1]*y[1]+A[0][2]*y[2];u[1]=A[3][0]*y[0]+A[2][1]*y[1]+A[1][2]*y[2]+A[0][3]*y[3];u[N-2]=A[4][N-4]*y[N-4]+A[3][N-3]*y[N-3]+A[2][N-2]*y[N-2]+A[1][N-1]*y[N-1];u[N-1]=A[4][N-3]*y[N-3]+A[3][N-2]*y[N-2]+A[2][N-1]*y[N-1];for(i=2;i<N-2;i++)u[i]=A[4][i-2]*y[i-2]+A[3][i-1]*y[i-1]+A[2][i]*y[i]+A[1][i+1]*y[i+1]+A[0][i+2]*y[i+2]; }void Get_value() //获得迭代后特征值{value2=value1;value1=max_sign*u[max_position];}void Check_value() //幂法第二迭代格迭代{Init_u();Get_max();Get_y();Get_u();Get_value();while(1){Get_max();Get_y();Get_u();Get_value();if(fabs((value2-value1)/value1)<e)break;}}void The_value() //获取绝对值最大的特征值λ_501 {Check_value();value_abs_max=value1;}void The_Other_value() //获取特征值λ_1{int i;float value_temp=value1;for(i=0;i<N;i++){A[2][i]-=value_temp;}Check_value();value1+=value_temp;if(value1<value_temp){value_1=value1;value_N=value_temp;}else{value_N=value1;value_1=value_temp;}}int min(int a,int b) //两值中取最小{if(a<b)return a;elsereturn b;}int max(int a,int b) //两值中取最大{if(a<b)return b;elsereturn a;}void Resolve_LU(){int k,i,j,t;float temp;for(k=1;k<=N;k++){for(j=k;j<=min(k+S,N);j++){temp=0;for(t=max(max(1,k-R),j-S);t<=k-1;t++)temp+=A[k-t+S][t-1]*A[t-j+S][j-1];A[k-j+S][j-1]=A[k-j+S][j-1]-temp;}for(i=k+1;i<=min(k+R,N);i++){temp=0;for(t=max(max(1,i-R),k-S);t<=k-1;t++)temp+=A[i-t+S][t-1]*A[t-k+S][k-1];A[i-k+S][k-1]=(A[i-k+S][k-1]-temp)/A[S][k-1];}}}void Back_substitution()//方程组回代过程{int i,t;float temp=0;for(i=2;i<N+1;i++){for(t=max(1,i-R);t<i;t++)y[i-1]-=A[i-t+S][t-1]*y[t-1];}u[N-1]=y[N-1]/A[S][N-1];for(i=N-1;i>0;i--){temp=0;for(t=i+1;t<=min(i+S,N);t++)temp+=A[i-t+S][t-1]*u[t-1];u[i-1]=(y[i-1]-temp)/A[S][i-1];}}double Det_matrix() //求矩阵行列式值{int i;double det=1;Init_matrix_A();Resolve_LU();for(i=0;i<N;i++)det=det*A[2][i];return det;}float Get_norm() //获得迭代向量模{int i;float normal=0;for(i=0;i<N;i++)normal+=u[i]*u[i];normal=sqrt(normal);return normal;}void Get_yy(float normal) //迭代向量单位化{int i;for(i=0;i<N;i++){y[i]=u[i]/normal;yy[i]=y[i];}}void Get_value_s() //获得绝对值最小的特征值{int i;value2=value1;value1=0;for(i=0;i<N;i++)value1+=yy[i]*u[i];value1=1/value1;}void Value_min() //反幂法求绝对值最小的特征值{float norm=0;int count=0;value1=0,value2=0;Init_u();norm=Get_norm();Get_yy(norm);Back_substitution();Get_value_s();while(count<10000){count++;norm=Get_norm();Get_yy(norm);Back_substitution();Get_value_s();if(fabs((value2-value1)/value1)<e)break;}value_s=value1;}float Get_cond_A() //求矩阵条件数{float cond1;cond1=fabs(value_abs_max/value_s);return cond1;}void Value_translation_min() //偏移条件下反幂法求特征值{int i,k;float tr;for(k=1;k<K+1;k++){tr=value_1+k*(value_N-value_1)/40;Init_matrix_A();for(i=0;i<N;i++)A[2][i]-=tr;Resolve_LU();Value_min();value_s+=tr;printf("k=%d =>>>λi%d=%.13e\n",k,k,value_s);}}void main(){float cond;double value_det;printf("Contactme:****************\n");Init_matrix_A(); //初始化矩阵AThe_value(); //获取绝对值最大的特征值λ_501 The_Other_value(); //获取特征值λ_1printf("λ1=%.13e\n",value_1);printf("λ501=%.13e\n",value_N);value_det=Det_matrix(); //求矩阵行列式值Value_min(); //反幂法求绝对值最小的特征值printf("λs=%.13e\n",value_s);cond=Get_cond_A(); //求矩阵条件数Value_translation_min();//偏移条件下反幂法求特征值printf("cond_A=%.13e\n",cond);printf("value_det=%.13e\n",value_det);}3、程序运行结果:4、迭代初始向量的选取对计算结果的影响:本次计算实习求矩阵A的具有某些特征的特征值,主要用到的方法是幂法和反幂法,这两种方法从原理上看都是迭代法,因此迭代初始向量的选择对计算结果会产生一定影响,主要表现在收敛速度上。