求矩阵特征值算法及程序简介

求矩阵特征值算法及程序简介1.幂法1、幂法规范化算法（1）输入矩阵A、初始向量)0(，误差eps；（2）1k；（3）计算)1()(kkAV；（4）)max(,)max()1(1)(kkkkVmVm；（5）kkkmV/)()(；（6）如果epsmmkk1,则显示特征值1和对应的特征向量)1(x),终止；（7）1kk,转（3）注：如上算法中的符号)max(V表示取向量V中绝对值最大的分量。本算法使用了数据规范化处理技术以防止计算过程中出现益出错误。2、规范化幂法程序Clear[a,u,x];a=Input[系数矩阵A=];u=Input[初始迭代向量u(0)=];n=Length[u];eps=Input[误差精度eps=];nmax=Input[迭代允许最大次数nmax=];fmax[x_]:=Module[{m=0,m1,m2},Do[m1=Abs[x[[k]]];If[m1m,m2=x[[k]];m=m1],{k,1,Length[x]}];m2]v=a.u;m0=fmax[u];m1=fmax[v];t=Abs[m1-m0]//N;k=0;While[teps&&knmax,u=v/m1;v=a.u;k=k+1;m0=m1;m1=fmax[v];t=Abs[m1-m0]//N;Print[k=,k,特征值=,N[m1,10],误差=,N[t,10]];Print[特征向量=,N[u,10]]];If[knmax,Print[迭代超限]]说明：本程序用于求矩阵A按模最大的特征值及其相应特征向量。程序执行后，先通过键盘输入矩阵A、迭代初值向量)0(、精度控制eps和迭代允许最大次数maxn，程序即可给出每次迭代的次数和对应的迭代特征值、特征向量及误差序列，它们都按10位有效数输出。其中最后输出的结果即为所求的特征值和特征向量序列。如果迭代超出maxn次还没有求出满足精度的根则输出迭代超限提示，此时可以根据输出序列判别收敛情况。程序中变量说明a:存放矩阵A；u:初始向量)0(和迭代过程中的向量)(k及所求特征向量；v:存放迭代过程中的向量)(kV；m1:存放所求特征值和迭代过程中的近似特征值；nmax:存放迭代允许的最大次数；eps:存放误差精度；fmax[x]:给出向量x中绝对值最大的分量；k:记录迭代次数；t1:临时变量；注：迭代最大次数可以修改为其他数字。3、例题与实验例1.用幂法求矩阵90688465441356133A的按模最大的特征值及其相应特征向量，要求误差410eps。解：执行幂法程序后在输入的四个窗口中按提示分别输入:{{133,6,135},{44,5,46},{-88,-6,-90}}、{1,1,1}、0.0001、20每次输入后用鼠标点击窗口的“OK”按扭，得如下输出结果：此结果说明迭代6次，求得误差为0.0000101442的按模最大的特征值为44.99999952,及其对应的一个特征向量:{1.000000000,0.3333333371,-0.6666666704}2.反幂法1、反幂法规范化算法（1）输入矩阵A、初始向量)0(，误差eps；（2）1k；（3）计算)1()(kkAV求出解)(kV；（4）)max(,)max()1(1)(kkkkVmVm；（5）kkkmV/)()(；（6）如果epsmmkk1,则显示特征值1和对应的特征向量)1(x),终止；（7）1kk,转（3）注：如上算法中解方程)1()(kkAV可以使用Dololittle分解法。本算法使用了数据规范化处理技术以防止计算过程中出现益出错误。2、规范化反幂法程序Clear[a,u,x];a=Input[系数矩阵A=];u=Input[初始迭代向量u(0)=];n=Length[u];eps=Input[误差精度eps=];nmax=Input[迭代允许最大次数nmax=];fmax[x_]:=Module[{m=0,m1,m2},Do[m1=Abs[x[[k]]];If[m1m,m2=x[[k]];m=m1],{k,1,Length[x]}];m2];v=a.u;a1=Inverse[a];m0=fmax[u];m1=fmax[v];t=Abs[m1-m0]//N;k=0;While[teps&&knmax,u=v/m1;v=a1.u;k=k+1;m0=m1;m1=fmax[v];t=Abs[m1-m0]//N;t1=Abs[1/m1-1/m0]//N;Print[k=,k,特征值=,N[1/m1,10],误差=,N[t1,10]];Print[特征向量=,N[u,10]]];If[knmax,Print[迭代超限]]说明：本程序用于求矩阵A按模最小的特征值及其相应特征向量。程序执行后，先通过键盘输入矩阵A、迭代初值向量)0(、精度控制eps和迭代允许最大次数maxn，程序即可给出每次迭代的次数和对应的迭代特征值、特征向量及误差序列，它们都按10位有效数输出。其中最后输出的结果即为所求的特征值和特征向量序。如果迭代超出maxn次还没有求出满足精度的根则输出迭代超限提示，此时可以根据输出序列判别收敛情况。程序中变量说明a：存放矩阵Au：初始向量)0(和迭代过程中的向量)(k及所求特征向量v:存放迭代过程中的向量)(kVa1：存放逆矩阵1Am1:存放所求特征值和迭代过程中的近似特征值nmax:存放迭代允许的最大次数eps:存放误差精度fmax[x]:给出向量x中绝对值最大的分量k:记录迭代次数t1:临时变量注：迭代最大次数可以修改为其他数字。3、例题与实验例3.用反幂法求矩阵131111322A的按模最小的特征值及其相应特征向量，要求误差510eps。解：执行幂法程序后在输入的四个窗口中按提示分别输入:{{2.,-2.,3.},{1,1.,1},{1.,3,-1}}、{1,0,1}、0.00001、100每次输入后用鼠标点击窗口的“OK”按扭，得如下输出结果：注意到本题按模最小的特征值为1，因此求解效果较满意。3.Jacobi方法1、Jacobi旋转法算法1、输入矩阵A，误差；2、,2,1kFor；（2.1）选取)1()1(maxkijjikpqaa记录qp,（2.2）由4,22tan)1()1()1(kqqkppkpqaaa确定旋转角，获得旋转矩阵,,qpJ；（2.3）qpjaaaaakpjkjpkqjkpjkpj,,sincos)()()1()1()(；（2.4）qpjaaaaakqjkjpkqjkpjkqj,,cossin)()()1()1()(；（2.5）qpjiaakijkij,,)1()(（2.6）cossin2sincos)1(2)1(2)1()(kpqkqqkppkppaaaa（2.7）cossin2cossin)1(2)1(2)1()(kpqkqqkppkqqaaaa（2.8）计算jikijkaAE2)()(（2.9）如果)(kAE，输出对角矩阵),,,(21ndiagD和特征向量矩阵J，停止注：如上算法中)()(kijkaA，)()()0(0ijijaaA。2、Jacobi算法程序Clear[a,bb];a=Input[矩阵A=];n=Input[矩阵阶数n=];eps=Input[误差精度eps=];nmax=Input[迭代允许最大次数nmax=];k=0;bb=IdentityMatrix[n];ea=Sum[a[[i,j]]^2,{i,1,n},{j,1,n}]-Sum[a[[i,i]]^2,{i,1,n}]//N;While[eaeps&&knmax,m=0;Print[迭代次数k=,k];Do[If[Abs[a[[i,j]]]m,m=Abs[a[[i,j]]];p=i;q=j],{i,1,n},{j,i+1,n}];mu=a[[p,p]]-a[[q,q]];If[mu0,thi=Pi/4,thi=ArcTan[2*a[[p,q]]/mu]/2];s=Sin[thi]//N;c=Sqrt[1-s^2];a1=bb[[p]];bb[[p]]=c*bb[[p]]+s*bb[[q]];bb[[q]]=-s*a1+c*bb[[q]];pp=a[[p,p]]*c*c+a[[q,q]]*s*s+2a[[p,q]]*s*c;qq=a[[p,p]]*s*s+a[[q,q]]*c*c-2a[[p,q]]*s*c;Do[a1=a[[p,j]];a[[p,j]]=c*a[[p,j]]+s*a[[q,j]];a[[j,p]]=a[[p,j]];a[[q,j]]=c*a[[q,j]]-s*a1;a[[j,q]]=a[[q,j]],{j,1,n}];a[[p,p]]=pp;a[[q,q]]=qq;a[[p,q]]=0;a[[q,p]]=0;ea=Sum[a[[i,j]]^2,{i,1,n},{j,1,n}]-Sum[a[[i,i]]^2,{i,1,n}]//N;k=k+1;Print[误差=,ea];Print[相似矩阵A=];Print[MatrixForm[a]];Print[特征向量J];Print[MatrixForm[Transpose[bb]]]];说明本程序用于求对称矩阵A的所有特征值及其相应特征向量。程序执行后，先通过键盘输入矩阵A、矩阵阶数n、精度控制eps和迭代允许最大次数maxn，程序即可给出每次迭代的次数和对应的迭代特征值、特征向量及误差序列。其中最后输出的结果即为所求的特征值和特征向量。如果迭代超出maxn次还没有求出满足精度的根则输出迭代超限提示。此外，输出的特征值矩阵可以不是真正的对角矩阵，但它们的主对角元素就是满足要求的所有特征值。程序中变量说明a：存放矩阵A及其相似变换过程中的kAbb:存放特征向量矩阵J的转置v:存放迭代过程中的向量)(kVm1:存放所求特征值和迭代过程中的近似特征值nmax:存放迭代允许的最大次数eps:存放误差精度k:记录迭代次数t1,mu,s,c,ea,p,q,m:临时变量a1,qq,pp：临时向量3、例题与实验例4.用Jacobi方法求矩阵210121012A的所有特征值及其相应特征向量，要求误差610eps。解：执行Jacobi方法程序后在输入的四个窗口中按提示分别输入:{{2,-1,0},{-1,2,-1},{0,-1,2}}、3、0.000001、30每次输入后用鼠标点击窗口的“OK”按扭，得如下输出结果：此结果说明迭代6次，求得误差为1.1668510-9的矩阵A的所有特征值：2,41421.30.585786,321及对应的特征向量{0.500012，0.707107，0.499988}，{-0.5，0.707107，-0.5}，{-0.707098，0.0000120684，0.707115}