北航数值分析大作业第一题幂法与反幂法

1、数值分析计算实习题目第一题：1. 算法设计方案(1) 1， 501 和 s的值。1)首先通过幂法求出按模最大的特征值t1，然后根据 t1 进行原点平移求出另一特征值t2，比较两值大小， 数值小的为所求最小特征值 1，数值大的为是所求最大特征值501。2)使用反幂法求 s，其中需要解线性方程组。因为A 为带状线性方程组，此处采用 LU分解法解带状方程组。( 2)与 k = 1+k 501 1 最接近的特征值 ik 。k 1 40 ik通过带有原点平移的反幂法求出与数 k 最接近的特征值 ik 。(3) cond(A) 2和 detA 。1) cond(A)2 = 1 ，其中 1和 n 分别是按模

2、最大和最小特征值。2)利用步骤( 1)中分解矩阵 A 得出的 LU矩阵， L为单位下三角阵 ,U 为上三角阵，其 中 U 矩阵的主对角线元素之积即为 detA 。由于 A 的元素零元素较多，为节省储存量，将 A 的元素存为 6501的数组中，程序中 采用 get_an_element() 函数来从小数组中取出 A 中的元素。2. 全部源程序#include #include void init_a();/ 初始化 Adouble get_an_element(int,int);/ 取 A 中的元素函数double powermethod(double);/ 原点平移的幂法double inve

3、rsepowermethod(double);/ 原点平移的反幂法int presolve(double);/ 三角 LU 分解int solve(double ,double );/ 解方程组int max(int,int);int min(int,int);double (*u)502=new double502502;/ 上三角 U 数组double (*l)502=new double502502;/ 单位下三角 L 数组double a6502;/ 矩阵 Aint main()int i,k;double lambdat1,lambdat2,lambda1,lambda501,lam

4、bdas,mu40,det;double lambda40; init_a();/ 初始化 A lambdat1=powermethod(0); lambdat2=powermethod(lambdat1); lambda1=lambdat1lambdat2?lambdat1:lambdat2; presolve(0);lambdas=inversepowermethod(0);det=1; for(i=1;i=501;i+) det=det*uii;for (k=1;k=39;k+) muk=lambda1+k*(lambda501-lambda1)/40; presolve(muk); l

5、ambdak=inversepowermethod(muk);printf( 所有特征值如下 n);printf( =%1.11e =%1.11en,lambda1,lambda501); printf( s=%1.n11,leambdas);printf(cond(A)=%1.11en,fabs(lambdat1/lambdas); printf(detA=%1.11e n,det);for (k=1;k=39;k+)printf( i%d=%1.11e ,k,lambdak);if(k % 3=0) prinntf();delete u;delete l;/ 释放堆内存return 0;v

6、oid init_a()/ 初始化 Aint i;for (i=3;i=501;i+) a1i=a5502-i=-0.064;for (i=2;i=501;i+) a2i=a4502-i=0.16;for (i=1;i=501;i+) a3i=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); double get_an_element(int i,int j)/ 从 A 中节省存储量的提取元素方法 if (fabs(i-j)=2) return ai-j+3j;else return 0;double powermethod(double offset)/

7、幂法int i,x1;double u502,y502;double beta=0,prebeta=-1000,yita=0;for (i=1;i=501;i+)ui=1,yi=0;/ 设置初始向量 ufor (int k=1;k=10000;k+)yita=0;for (i=1;i=501;i+) yita=sqrt(yita*yita+ui*ui);for (i=1;i=501;i+) yi=ui/yita;for (x1=1;x1=501;x1+)ux1=0;for (int x2=1;x2=501;x2+)ux1=ux1+(x1=x2)?(get_an_element(x1,x2)-o

8、ffset):get_an_element(x1,x2)*yx2; prebeta=beta;beta=0;for (i=1;i=501;i+) beta=beta+ yi*ui;if (fabs(prebeta-beta)/beta)=1e-12) printf(offset=%f lambda=%f k=%dn,offset,(beta+offset),fabs(prebeta-beta)/beta),k);break;/ 输出中间过程 ,包括偏移量 ,误差 ,迭代次数return (beta+offset);double inversepowermethod(double offset)

9、/ 反幂法int i;double u502,y502;double beta=0,prebeta=0,yita=0;for (i=1;i=501;i+)ui=1,yi=0; / 设置初始向量 ufor (int k=1;k=10000;k+)yita=0;for (i=1;i=501;i+) yita=sqrt(yita*yita+ui*ui);for (i=1;i=501;i+) yi=ui/yita;solve(u,y); prebeta=beta;beta=0;for (i=1;i=501;i+) beta=beta+ yi*ui;beta=1/beta;if (fabs(prebet

10、a-beta)/beta)=1e-12) printf(offset=%f lambda=%f k=%dn,offset,(beta+offset),fabs(prebeta-beta)/beta),k);break;/ 输出中间过程 ,包括偏移量 ,误差 ,迭代次数return (beta+offset);err=%eerr=%eint presolve(double offset)/ 三角 LU 分解int i,k,j,t;double sum;for (k=1;k=501;k+)for (j=1;j=501;j+) ukj=lkj=0; if (k=j) lkj=1; / 初始化 LU

11、矩阵for (k=1;k=501;k+)for (j=k;j=min(k+2,501);j+)sum=0;for (t=max(1,max(k-2,j-2) ; t=(k-1) ; t+) sum=sum+lkt*utj;ukj=(k=j)?(get_an_element(k,j)-offset):get_an_element(k,j)-sum; if (k=501) continue;for (i=k+1;i=min(k+2,501);i+)sum=0;for (t=max(1,max(i-2,k-2);t=(k-1);t+) sum=sum+lit*utk;lik=(i=k)?(get_a

12、n_element(i,k)-offset):get_an_element(i,k)-sum)/ukk; return 0;int solve(double x,double b)/ 解方程组int i,t;double y502;double sum;y1=b1;for (i=2;i=501;i+)sum=0;for (t=max(1,i-2);t=1;i-)sum=0;for (t=i+1;ty?x:y);int min(int x,int y)return (xy?x:y);3. 计算结果 结果如下图所示部分中间结果：给出了偏移量 (offset) ，误差 (err) ，迭代次数 (k)

13、4. 讨论迭代初始向量的选取对计算结果的影响,并说明原因使用 ui=1( i=1,2,.,501)作为初始向量进行迭代， 可得出以上结果。 经过 Mathematica 计算验证结果正确。现修改初始向量 u1=1 ，ui=0，(i=2,3,.,501)。得出结果此结果与正确结果相差较多。令初始向量 um=1 ，un=0 ，(m=1,2,.,250 n=251,252,501) ，得出结果：此结果也与正确结果相差较多。但与上次结果相比，更加靠近准确值一些。再增加初始向量 u 中等于 1 的元素个数，可以发现其结果更加靠近准确值。经验证， 只有当不为 0 元素的个数达到比较高的一个值时， 才能得到精确的收敛结果， 与元素的绝对 值大小无关。分析算法，设 1为按模最大特征值，观察式子u0 1x1